*** *** ***

# 3D Spectrum Mapping and Reconstruction under Multi-Radiation Source
Scenarios

Jie Wang 1122 Zhipeng Lin 1122 Qiuming<EMAIL_ADDRESS>1122** Qihui
Wu 1122 Tianxu Lan 1122 Yi Zhao 1122 Yunpeng Bai 1122 Weizhi Zhong 2233 The
Key Laboratory of Dynamic Cognitive System of Electromagnetic Spectrum Space,
Nanjing University of Aeronautics and Astronautics, 211106, Nanjing, China
College of Electronic and Information Engineering, Nanjing University of
Aeronautics and Astronautics, 211106, Nanjing, China  The Key Laboratory of
Dynamic Cognitive System of Electromagnetic Spectrum Space, College of
Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing
211106, China

###### Abstract

Spectrum map construction, which is crucial in cognitive radio (CR) system,
visualizes the invisible space of the electromagnetic spectrum for spectrum-
resource management and allocation. Traditional reconstruction methods are
generally for two-dimensional (2D) spectrum map and driven by abundant
sampling data. In this paper, we propose a data-model-knowledge-driven
reconstruction scheme to construct the three-dimensional (3D) spectrum map
under multi-radiation source scenarios. We firstly design a maximum and
minimum path loss difference (MMPLD) clustering algorithm to detect the number
of radiation sources in a 3D space. Then, we develop a joint location-power
estimation method based on the heuristic population evolutionary optimization
algorithm. Considering the variation of electromagnetic environment, we self-
learn the path loss (PL) model based on the sampling data. Finally, the 3D
spectrum is reconstructed according to the self-learned PL model and the
extracted knowledge of radiation sources. Simulations show that the proposed
3D spectrum map reconstruction scheme not only has splendid adaptability to
the environment, but also achieves high spectrum construction accuracy even
when the sampling rate is very low.

###### keywords:

cognitive radio; 3D spectrum map; map reconstruction; path loss model;
radiation source;

##  1 INTRODUCTION

Figure 1: (a): Realistic 3D spectrum map under an urban scenario; (b):
Reconstruction model of 3D spectrum map.

With the extensive growth of wireless frequency devices, the problem of
electromagnetic spectrum resource deficit is increasing severe. In this case,
the cognitive radio (CR) technique is proposed, which is beneficial for
efficient management and allocation for rational spectrum resources [1, 2, 3,
4, 5, 6, 7, 8]. A key task in CR is using spectrum awareness to accurately
determine the spectrum situation, including the spectrum environments, the
received signal strength, the spectrum busy or the idle states and the
spectrum access protocols, etc [9, 10, 11, 12, 13]. Radio environment map
(REM), also known as spectrum map, comprising time, frequency, signal
strength, location and other multi-dimensional spectrum information can
visualize spectrum state awareness results on a geographical map [14, 15]. REM
is useful for numerous applications in wireless communication, e.g.,
opportunistic spectrum access (OSA), spectrum management, and interference
coordination [9, 16, 17, 18, 19]. However, considering the limited capability
of sensing devices and the costs of measurement, it is necessary to construct
an accurate and robust REM with limited sampling data.

Many spectrum map reconstruction methods have emerged in recent years. In
[20], the authors constructed a REM to process the data collected by
heterogeneous spectrum sensors. A novel machine learning approach based on
generative adversarial networks (GANs) was presented in [21], which considers
the completion of spectrum map from the perspective of image processing. The
authors in [22] converted the estimation task into an image reconstruction
task by image color mapping. Nevertheless, these methods only focus on two-
dimensional (2D) ground spectrum mapping, which cannot meet the needs of
three-dimensional (3D) spectrum monitoring. Recently, electromagnetic spectrum
monitoring is facing a huge challenge that the monitoring space extends from
the land to the air. There have been many high-flexibility air spectrum
sensing devices, such as unmanned aerial vehicle (UAV), which have excellent
spatial freedom and abundant spectrum access opportunities. Compared to fixed
ground sensing devices, they can more easily capture the dynamic changes of
spectrum [23], and can be used to construct the spectrum map of 3D spaces.

Typically, spectrum map reconstruction can be divided into data-driven methods
and model-driven methods. The data-driven methods are usually regarded as a
matrix completion problem by interpolation methods. In [24], the authors
proposed an inverse distance weighted (IDW) interpolation construction method.
It does not take any physical fact into account, degrading the performance in
terms of REM construction of urban environments with severe shadow fading.
Kriging interpolation is another spatial interpolation method [25]. In [26],
the authors proposed an optimal data-driven spectrum mapping method based on
distributed Kriging. Moreover, by extending the low-rank property from
matrices to tensors, a tensor completion method combining the prediction model
was proposed in [27]. However, these methods can only achieve satisfactory
performance when a large amount of data is collected, and they cannot well
reflect the essence of electromagnetic propagations. In addition, sampled
locations have a great influence on the completion performance. If the sampled
location is far away from the radiation source, the completion performance
will be poor.

Model-driven methods based on radiation source information and electromagnetic
propagation model can be addressed in [28, 29]. They can achieve high mapping
accuracy and are applicable to sparse sampling data. The authors in [28]
proposed a LocatIon Estimation based (LIvE) algorithm. The spectrum map is
reconstructed by estimating the location of primary user (PU) based on the RSS
value. However, the power of PU and the PL model, are not easy to obtain. In
[29], authors proposed an indirect model-driven reconstruction method by
taking the location, power and antenna pattern. In [30], a least absolute
shrinkage and selection operator (LASSO) method was proposed by using the PL
model for each radiation source. Nevertheless, these model-driven methods are
not applicable to the sce- narios with multi-radiation sources.

In this paper, we propose a data-model-knowledge-driven 3D spectrum
construction method under multi-radiation source scenarios. We first extract
the knowledge of radiation sources in the region of interest (ROI) based on
the sampling data, and then combine it with the propagation model to build a
spectrum map that truly presents the real-time spectrum situation, which
combines sampling data, propagation model and radiation source knowledge to
achieve spectrum mapping.

The main contributions of this paper are summarized as follows:

  * •

A maximum and minimum path loss difference (MMPLD) clustering algorithm is
developed to detect the number of multi-radiation sources in the 3D space.
Based on the propagation model, the clustering algorithm identifies the number
of radiation sources in the measurement area, and thus, can achieve the
spectrum map construction of 3D multi-radiation source scenarios.

  * •

The shuffled frog leaping algorithm (SFLA) is proposed to estimate the high-
dimensional parameters of the locations and the powers of radiation sources.
SFLA can efficiently execute the process of multi-source fusion RSS data and
achieve the parameters estimation of radiation sources.

  * •

Considering the dynamic nature of the urban environment, we self-learn the
propagation model by modifying the PL model with sparse sampled data. The
model self-learning algorithm can precisely capture the characteristics of
electromagnetic propagation in real time.

The remainder of this paper is organized as follows. Section II introduces and
elaborates the 3D spectrum reconstruction model. In Section III, a MMPLD-SFLA-
based self-learning reconstruction mechanism of 3D spectrum map under multi-
radiation source scenarios is proposed. Then, Section IV presents the
simulation results and analyses. Finally, Section V draws conclusions.

##  2 SYSTEM MODEL

In this paper, we consider to construct the 3D spectrum map of multi-radiation
source urban scenario as shown in Figure 1. In the considered urban scenario,
there are multi-radiation sources that have high priorities or legally
licenses to utilize the frequency bands. Users opportunistically access the
available spectrum imperceptibly. We construct a complete 3D spectrum map for
rational use of spectrum resources, so that we can obtain the spectrum
occupation status in the ROI and achieve efficient and reliable dynamic
spectrum access in the future work.

To construct the 3D spectrum map, we first discretize the whole ROI into small
cubes. Each cube is colored according to the RSS of its central position,
where red cubes represent high RSS values and blue cubes represent low RSS
values. We establish a three-dimensional cartesian coordinate system according
to the starting and ending points of the ROI and set up a spectrum tensor
$\chi\in{\Re^{{N_{1}}\times{N_{2}}\times{N_{3}}}}$, in which ${N_{1}}$,
${N_{2}}$ and ${N_{3}}$ indicate the grid number in $x,y,z$ dimensions,
respectively.

We proceed to mark the variables to be used. Note that it is impossible to
deploy a large number of sensors or spectrum sensing devices for all cubes of
the spectrum map due to high cost. Therefore, we complete the spectrum map by
sampling data in the ROI as shown in Figure 1. Suppose there are $N$ sampled
points in the ROI, and each of them can provide its location and RSS values.
$M$ radiation sources are located at
${\bf{t}}_{j}^{t}=\left[{x_{{}_{j}}^{t},y_{{}_{j}}^{t},z_{{}_{j}}^{t}}\right],j=1,2,\cdots
M$ with the transmission power ${P_{j}}^{t}$. Let
${{\bf{T}}^{r}}\buildrel\Delta\over{=}\left[{{\bf{t}}_{{}_{1}}^{r},{\bf{t}}_{{}_{2}}^{r},\cdots,{\bf{t}}_{{}_{N}}^{r}}\right]$
be the collection of the sampling positions.
${\bf{t}}_{{}_{i}}^{r}=\left[{x_{{}_{i}}^{r},y_{{}_{i}}^{r},z_{{}_{i}}^{r}}\right],i=1,2,\ldots,N$
are the coordinates of the cube in $\chi$.
${{\bf{P}}^{r}}\buildrel\Delta\over{=}\left[{{P_{1}}^{r},{P_{2}}^{r},\cdots,{P_{N}}^{r}}\right]$
is the collection of the sampled RSS with the value ${{\bf{T}}^{r}}$. We
define
$\mathbb{R}=\left\\{{{{\bf{R}}_{1}},{{\bf{R}}_{2}},\ldots,{{\bf{R}}_{N}}}\right\\}$
as a set of samples, where ${{\bf{R}}_{i}}=({\bf{t}}_{i}^{r},P_{i}^{r})$.

##  3 3D SPECTRUM MAPPING AND RECONSTRUCTION

###  3.1 The Procedure of 3D Spectrum Map

Figure 2: The flowchart of 3D spectrum map reconstruction.

Our 3D spectrum map reconstruction consists of four steps: $i)$ Scenario
setting up and ROI modeling; $ii)$ Sampling data collection; $iii)$ Knowledge
extraction of the radiation sources; $iv)$ Self-learning of the propagation
model and reconstruction of the spectrum map. At the first step, we set up and
model the ROI. We collect sampling data from the spectrum sensing devices
within the ROI. We then design the MMPLD-clustering algorithm and SFLA-
optimization algorithm to extract the knowledge from the ROI, including the
total number of radiation sources, and the location and power of each
radiation source. Finally, given the complexity of urban environment, we learn
the PL model based on sampling data to make the algorithm have self-learning
ability. The spectrum map is constructed based on the data, model and
knowledge. The flow chart of the whole algorithm is shown in Figure 2.

###  3.2 3D Radiation Sources Number Detection

The challenge of constructing 3D multi-radiation source spectrum map is to
mine the knowledge of the radiation sources. In this section, we first briefly
review the $K$-means clustering algorithm, and then, design the MMPLD-based
clustering algorithm to detect the number of radiation sources.

####  3.2.1 Review of $K$-means Clustering Algorithm

Clustering data based on the measurements of similarity is critical in data
process and pattern recognition. The basic idea of clustering data is to learn
a set of clustering centers from the data to minimize the sum of squared
errors between the data points and their nearest centers [31].

$K$-means is a classic clustering technique based on geometric centers and
partitions [32, 33]. It is based on Euclidean distance and considers that the
closer the distance between two objects is, the greater the similarity is.
However, the initial clustering centers are difficult to be determined and the
$K$ value needs to be known in advance. Therefore, $K$-means usually needs to
be rerun multiple times with different initializations to find a good
solution. The clustering algorithm based on the maximum and minimum distance
(MMD) criterion, which is developed from $K$-means, selects the farthest point
to be the clustering center. The MMD algorithm can dynamically determine
clustering number and avoid the shortcomings of $K$-means [34]. Note that a
hyperparameter threhold can be set determine whether the clustering process
needs to be continued. Nevertheless, the selection of the hyperparameter
threhold directly affects the partition results. Traditional clustering
algorithms takes the Euclidean distance as the basis, but it is not suitable
for spectrum data processing. Moreover, they cannot accurately determine the
number of clusters.

####  3.2.2 MMPLD-based 3D Radiation Sources Number Detection

Traditional clustering algorithms need to know the number of clusters in
advance. However, in our paper, the number of clusters is the parameter that
we need to estimate. The maximum and minimum distance (MMD) clustering
algorithm, which is developed from K-means, can dynamically determine
clustering number. Considering the traditional the Euclidean distance
clustering criterion is not suitable for spectrum data, we propose a maximum
and minimum path loss difference (MMPLD) criterion by exploiting the
correlation of spectrum data. Besides, a criterion function is designed to
determine the number of clusters at the end of the algorithm to balance the
complexity of parameter settings. Given the electromagnetic propagation law,
the RSS of a particular position is more consistent with the propagation law
of the nearest radiation source than other radiation sources. The main
procedure of the proposed MMPLD algorithm includes three steps: $i)$ The path
loss difference (PLD) is calculated as the basic criterion; $ii)$ The
selection of clustering centers and data classification are executed by
maximum and minimum PLD criterion; $iii)$ The convergence condition is
designed for the algorithm by calculating the criterion function of clustering
number $K$ to determine the optimal clustering number $K$.

The PL model of a free space can be expressed as

${L_{0}}=32.4+20{\log_{10}}\left({{f_{\rm{c}}}}\right)+20{\log_{10}}\left(d\right),$
(1)

where ${L_{0}}$ is the path loss in free space. $d$ and ${f_{c}}$ denote the
distance between sampled position and radiation source and frequency,
respectively.

We define $\mathbb{C}$ as the set of the clustering centers. Firstly, we
choose the ${{\bf{R}}_{n}}\in\mathbb{R}$, which has the highest
$P_{{}^{n}}^{r}$ of ${{\bf{P}}^{r}}$ as the first clustering center, denoted
by ${C_{1}}$, ${C_{1}}\in\mathbb{C}$. According to (1), we can obtain the
theoretical PL value $L_{i}^{{C_{1}}},i=1,2,3\cdots,N$ between ${C_{1}}$ and
the $i-th$ sampling point. Collecting all the values, we have
${{\bf{L}}^{{C_{1}}}}$, as given by

${{\bf{L}}^{{C_{1}}}}=\left\\{{L_{{}^{2}}^{{C_{1}}},L_{{}^{3}}^{{C_{1}}},\ldots,L_{{}^{N}}^{{C_{1}}}}\right\\}.$
(2)

${\bf{P}}_{\Delta}^{{C_{1}}}$ collects all the difference RSS values between
${C_{1}}$ and all sampling points, which can be expressed as

$\begin{array}[]{l}{\bf{P}}_{\Delta}^{{C_{1}}}=\left\\{{\Delta
P_{1}^{{}_{{C_{1}}}},\Delta P_{2}^{{}_{{C_{1}}}},\cdots,\Delta
P_{N}^{{}_{{C_{1}}}}}\right\\}\\\
=\left\\{\begin{array}[]{l}\left|{{P_{{C_{1}}}}^{r}-{P_{1}}^{r}}\right|,\left|{{P_{{C_{1}}}}^{r}-{P_{2}}^{r}}\right|,\\\
\ldots\left|{{P_{{C_{1}}}}^{r}-{P_{N}}^{r}}\right|\end{array}\right\\},\end{array}$
(3)

where ${P_{{C_{1}}}}^{r}$ is the RSS of the first clustering center ${C_{1}}$.
${P_{i}}^{r}\begin{array}[]{*{20}{c}},&{i=1,2,\ldots,N}\end{array}$ is the RSS
value of the $i$-th sampling point, $\Delta P_{i}^{{}_{{C_{1}}}}$ is the
absolute difference values of RSS between the first clustering center
${C_{1}}$ and the $i$-th sampling point, and $\Delta P_{i}^{{}_{{C_{1}}}}$ is
the actual PL value. We define error set ${{\bf{\delta}}^{{C_{1}}}}$ as
clustering criteria, according to the following formula:

$\begin{array}[]{l}{{\bf{\delta}}^{{C_{1}}}}=\left\\{{\delta_{{}_{1}}^{{C_{1}}},\delta_{{}_{2}}^{{C_{1}}},\ldots,\delta_{{}_{N}}^{{C_{1}}}}\right\\}\\\
=\left\\{\begin{array}[]{l}\left|{L_{1}^{{C_{1}}}-\Delta
P_{1}^{{C_{1}}}}\right|,\left|{L_{2}^{{C_{1}}}-\Delta
P_{2}^{{C_{1}}}}\right|,\\\ \cdots,\left|{L_{N}^{{C_{1}}}-\Delta
P_{N}^{{C_{1}}}}\right|\end{array}\right\\},\end{array}$ (4)

where $\delta_{{}_{I}}^{{C_{1}}}$ is the PLD between the first clustering
center ${C_{1}}$ and the $i$-th sampling point, and
$\delta_{{}_{I}}^{{C_{1}}}$ is the absolute error between theoretical PL and
actual PL. We choose the max value of ${{\bf{\delta}}^{{C_{1}}}}$ to be the
second cluster center, denoted as ${C_{2}}$, as given by

${C_{2}}\leftarrow\max(d_{1}^{\min},d_{2}^{\min},\ldots,d_{N}^{\min}),$ (5)

where $d_{i}^{\min}=\delta_{{}_{i}}^{{C_{1}}},i=1,2,\ldots,N$ indicates the
PLD between the point and the center of the class, to which the point belongs.
We initiate $d_{i}^{\min}$ with $\delta_{{}_{i}}^{{C_{1}}}$ and its value is
constantly updated as the clustering process progresses. To be specific, if
$d_{n}^{\min}=\max(d_{1}^{\min},d_{2}^{\min},\ldots,d_{N}^{\min})$, the second
cluster center ${C_{2}}$ is the $n$-th point. Afterwards, we calculate
$L_{i}^{{C_{2}}}$ according to (1). The ${\bf{P}}_{\Delta}^{{C_{2}}}$ of the
second cluster center ${C_{2}}$ are

$\begin{array}[]{l}{\bf{P}}_{\Delta}^{{C_{2}}}=\left\\{{\Delta
P_{1}^{{}_{{C_{2}}}},\Delta P_{2}^{{}_{{C_{2}}}},\cdots,\Delta
P_{N}^{{}_{{C_{2}}}}}\right\\}\\\
=\left\\{\begin{array}[]{l}\left|{{P_{{C_{2}}}}^{r}-{P_{1}}^{r}}\right|,\left|{{P_{{C_{2}}}}^{r}-{P_{2}}^{r}}\right|,\\\
\ldots\left|{{P_{{C_{2}}}}^{r}-{P_{N}}^{r}}\right|\end{array}\right\\},\end{array}$
(6)

where $\Delta P_{i}^{{}_{{C_{2}}}}$ is the actual PL between ${C_{2}}$ and the
$i$-th sampling point. The PLD set ${{\bf{\delta}}^{{C_{2}}}}$ in the second
class can be obtained by

$\begin{array}[]{l}{{\bf{\delta}}^{{C_{2}}}}=\left\\{{\delta_{{}_{1}}^{{C_{2}}},\delta_{{}_{2}}^{{C_{2}}},\ldots,\delta_{{}_{N}}^{{C_{2}}}}\right\\}\\\
=\left\\{\begin{array}[]{l}\left|{L_{1}^{{C_{2}}}-\Delta
P_{1}^{{C_{2}}}}\right|,\left|{L_{2}^{{C_{2}}}-\Delta
P_{2}^{{C_{2}}}}\right|,\\\ \cdots,\left|{L_{N}^{{C_{2}}}-\Delta
P_{N}^{{C_{2}}}}\right|\end{array}\right\\},\end{array}$ (7)

where $\delta_{{}_{I}}^{{C_{2}}}$is the PLD between the second clustering
center ${C_{2}}$ and the $i$-th sampling point. We update $d_{i}^{\min}$ by

$d_{i}^{\min}=\min\left({\delta_{i}^{{C_{1}}},\delta_{i}^{{C_{2}}}}\right),i=1,2,\ldots,N.$
(8)

According to $d_{i}^{\min}$, we categorize the $i$-th sampling point into the
class with the least PLD. Specifically, if
$d_{i}^{\min}=\delta_{i}^{{C_{m}}},m=1,2$, then the $i$-th point is assigned
to class $m$. Furthermore, the $q$-th clustering center ${C_{q}}$ is selected
from the maximum value of $d_{i}^{\min}$, as given by

${C_{q}}\leftarrow\max(d_{1}^{\min},d_{2}^{\min},\ldots,d_{N}^{\min}).$ (9)

If $d_{n}^{\min}=\max(d_{1}^{\min},d_{2}^{\min},\ldots,d_{N}^{\min})$,
${C_{q}}$ is the $n-th$ point. We keep obtaining the PLD
${{\bf{\delta}}^{{C_{q}}}}$ of the sampling points in the $q-th$ class

$\begin{array}[]{l}{{\bf{\delta}}^{{C_{q}}}}=\left\\{{\delta_{{}_{1}}^{{C_{q}}},\delta_{{}_{2}}^{{C_{q}}},\ldots,\delta_{{}_{N}}^{{C_{q}}}}\right\\}\\\
=\left\\{\begin{array}[]{l}\left|{L_{1}^{{C_{q}}}-\Delta
P_{1}^{{C_{q}}}}\right|,\left|{L_{2}^{{C_{q}}}-\Delta
P_{2}^{{C_{q}}}}\right|,\\\ \cdots,\left|{L_{N}^{{C_{q}}}-\Delta
P_{N}^{{C_{q}}}}\right|\end{array}\right\\},\end{array}$ (10)

where $\delta_{{}_{I}}^{{C_{q}}}$ is the PLD between the $q$-th clustering
center ${C_{q}}$ and the $i$-th sampling point; $L_{i}^{{C_{q}}}$ is the
theoretical path loss value between ${C_{q}}$ and the $i$-th sampling point in
the $q$-th class, and $\Delta P_{i}^{{}_{{C_{q}}}}$ is the actual path loss
between ${C_{q}}$ and the $i$-th sampling point.

Meanwhile, we update the $d_{i}^{\min}$ of each point:

$d_{i}^{\min}=\min\left({\delta_{i}^{{C_{1}}},\delta_{i}^{{C_{2}}},\ldots,\delta_{i}^{{C_{q}}}}\right),i=1,2,\ldots,N.$
(11)

We reclassify the sampling points according to the $d_{i}^{\min}$ until the
end condition of the algorithm is satisfied. If
$d_{i}^{\min}=\delta_{i}^{{C_{m}}},m=1,2,\ldots,q$, then the $i$-th point is
assigned to the class $m$. There is a criterion that we select the clustering
centers according to the maximum PLD, and classify points according to the
minimum PLD. The final clustering number $K$ is the number of radiation
sources need to be detected. The selection of $K$ depends on whether the
current clustering status conforms to the maximum similarity within the
cluster and the minimum similarity between the clusters. In order to determine
the optimal value $K$, we define the criterion function $\varpi(*)$ as the
convergence condition of the algorithm.

For the $i$-th point, which is in class $\gamma$, we define the function
$\varsigma(i)$ of the point with in other classes, given by

$\varsigma(i)=\frac{1}{{\left({K-1}\right)}}\sum\limits_{{C_{\rm{q}}},{C_{\rm{q}}}\neq\gamma}^{K-1}{\delta_{i}^{{C_{q}}}},$
(12)

where $\delta_{i}^{{C_{q}}}$ is the PLD calculated by (15), and it does not
belong to the current class $\gamma$. Then for the number $K$ of clusters we
can define the criterion function $\varpi(K)$, given by

$\varpi(K)=\frac{1}{N}\sum\limits_{i=1}^{N}{\frac{{d_{i}^{\min}}}{{\varsigma(i)}}},$
(13)

where $d_{i}^{\min}$ is the PLD between the $i$-th sampling point and the
clustering center ${C_{\gamma}}$ of the class $\gamma$.

We set two hyperparameters as thresholds as the convergence conditions of the
algorithm, given by

$\left\\{\begin{array}[]{l}\left|{(\varpi(K)-\varpi(K-1))}\right|\leq{\sigma_{1}}\begin{array}[]{*{20}{c}},&{\begin{array}[]{*{20}{c}}{{\rm{if}}}&{1\leq
K<3}\end{array},}\end{array}\\\
\left|{\frac{{(\varpi(K)-\varpi(K-1))}}{{(\varpi(K-1)-\varpi(K-2))}}}\right|\leq{\sigma_{2}}\begin{array}[]{*{20}{c}},&{\begin{array}[]{*{20}{c}}{{\rm{if}}}&{K\geq
3}\end{array},}\end{array}\end{array}\right.$ (14)

where ${\sigma_{1}}$ and ${\sigma_{2}}$ are the thresholds of the clustering
algorithm. When $K$ satisfies the above inequality, the clustering algorithm
ends and the $K$ value is the number of radiation sources.

###  3.3 SFLA-based 3D Radiant Location

In this section, we estimate the parameters of the radiation sources including
locations and powers. Suppose the transmission power of the $j$-th radiation
source is ${P_{j}}^{t},j=1,2,\cdots K$ and the location of the $j$-th
radiation source is
${\bf{t}}_{j}^{t}=\left[{x_{{}_{j}}^{t},y_{{}_{j}}^{t},z_{j}^{t}}\right],j=1,2,\cdots
K$. For each radiation source, we define a coefficient ${\eta_{j}}$ related to
its propagation loss. On the basis of the electromagnetic propagation law, the
RSS value ${\hat{P}_{ij}}^{r}$ of the $i$-th of $N$ sampling points from the
$j$-th radiation source is

$\begin{array}[]{l}{{\hat{P}}_{ij}}^{r}={\eta_{j}}{P_{\rm{j}}}^{t}{d_{ij}}^{-\alpha},\\\
i=1,2,\ldots N\begin{array}[]{*{20}{c}},&{j=1,2,\ldots
K}\end{array}\end{array}$ (15)

where ${d_{ij}}$ is the Euclidean distance between the location
${\bf{t}}_{j}^{r}$ of the $i$-th sampling point and the location
${\bf{t}}_{j}^{t}$ of the $j$-th radiation source, and $\alpha$ is the path
loss index.

For each sampled point, we can only measure the combined RSS of all radiation
sources. Then its RSS value from all radiation sources are added up to
generate the measured data ${\hat{P}^{r}}_{i}$, as given by

${\hat{P}_{i}}^{r}=\sum\limits_{j=1}^{K}{{{\hat{P}}_{ij}}^{r}}=\sum\limits_{j=1}^{K}{{\eta_{j}}{P_{\rm{j}}}^{t}{d_{ij}}^{-\alpha}}\begin{array}[]{*{20}{c}},&{i=1,2,\ldots
N}\end{array},$ (16)

where $K$ is the number of radiation sources estimated in Section 3.2.2.

We define ${\mathbf{\theta}}$ as the unknown vector of $5K$ parameters
${\bf{\theta}}=\left[{{{\bf{\theta}}_{1}},{{\bf{\theta}}_{2}},\ldots,{{\bf{\theta}}_{K}}}\right],$
${{\bf{\theta}}_{j}}={\left[{{\eta_{j}},x_{j}^{t},y_{j}^{t},z_{j}^{t},{P_{j}}^{t}}\right]^{T}},j=1,2,\cdots,K$
and ${\bf{\Omega}}$ be the set of all acceptable values of ${\bf{\theta}}$.
Given $K$, the cost function is established to estimate parameter, and the
estimation of ${{\bf{\theta}}^{\bf{*}}}$ is a high dimensional optimization
problem, as given by

$\begin{array}[]{l}{{\bf{\theta}}^{\bf{*}}}=\arg\mathop{\min}\limits_{\theta\in{\bf{\Omega}}}sqrt\left({{{\sum\limits_{i=1}^{N}{\left({{P_{i}}^{r}-{{\hat{P}}_{i}}^{r}}\right)}}^{2}}}\right)\\\
=\arg\mathop{\min}\limits_{\theta\in{\bf{\Omega}}}sqrt\left({{{\sum\limits_{i=1}^{N}{\left({{P_{i}}^{r}-\sum\limits_{j=1}^{K}{{\eta_{j}}{P_{\rm{j}}}^{t}{d_{ij}}^{-\alpha}}}\right)}}^{2}}}\right),\\\
{\rm{s}}{\rm{.t}}.\begin{array}[]{*{20}{c}}{}\hfil&{{P_{j}}^{t}\geq
0,\begin{array}[]{*{20}{c}}{}\hfil&{j=1,2,\cdots,K}\end{array}}\end{array},\end{array}$
(17)

where ${P_{i}}^{r},i=1,2,\cdots,N$ is the RSS of the $i$-th sampled position.

The algorithms can be used to solve this high dimensional parameters
eatimation problem[35]. SFLA is a swarm intelligence optimization algorithm
based on post-heuristics. It is based on the evolution of individual memes and
the global information exchange by memes. SFLA has many advantages, such as
simple concept, few parameters, fast computing speed, strong global
optimization ability, and easy implementation.

We develop SFLA to estimate the parameters ${{\bf{\theta}}^{\bf{*}}}$ in a
$5K$ dimensional target search space. Firstly, $P$ frogs (solutions) are
randomly generated within ${\bf{\Omega}}$ to form the initial population,
denoted as set
${\bf{\Lambda}}=\left\\{{{{\bf{F}}_{1}},{{\bf{F}}_{2}},\ldots,{{\bf{F}}_{P}}}\right\\}$.
The $s$-th frog represents a possible solution
${{\bf{F}}_{s}}=\left[{{\bf{\theta}}_{s1}^{*},{\bf{\theta}}_{s2}^{*},\ldots,{\bf{\theta}}_{sK}^{*}}\right]$
of the problem, where
${\bf{\theta}}_{sj}^{*}={\left[{{\eta_{sj}},x_{sj}^{t},y_{sj}^{t},z_{sj}^{t},{P_{sj}}^{t}}\right]^{T}},j=1,2,\cdots,K$,
is the parameter estimated of the $s$-th frog. Furthermore, the fitness
function of the $s$-th frog can be defined according (24), given by

$\begin{array}[]{l}\varphi(s)=sqrt\left({{{\sum\limits_{i=1}^{N}{\left({{P_{i}}^{r}-\sum\limits_{j=1}^{K}{{\eta_{sj}}{P_{{\rm{sj}}}}^{t}{d_{sij}}^{-\alpha}}}\right)}}^{2}}}\right),\\\
s=1,2,\ldots,P.\end{array}$ (18)

Then we divide the whole population into $M$ subgroups in a descending order
of $\varphi$. The frog ranked first is divided into the first subgroup, the
frog ranked second is allocated into the second subgroup, the frog ranked $M$
is allocated into the $M$ subgroup, the frog ranked $M+1$ is allocated into
the first subgroup, the frog ranked $M+2$ is allocated into the second
subgroup, and so on until all frogs were allocated.

We carry out local depth search for each subpopulation. Specifically, in each
iteration of the subpopulation, we determine the frog with the worst $\varphi$
and the best $\varphi$ as ${{\bf{F}}_{nw}}$ and ${{\bf{F}}_{nb}}$,
respectively. The frog with the best $\varphi$ in global population is denoted
as ${{\bf{F}}_{g}}$. Only the worst individual ${{\bf{F}}_{nw}}$ in the
current subpopulation is updated and the update strategy is

$\begin{array}[]{l}{\varepsilon_{\rm{t}}}=\lambda\cdot\left({{{\bf{F}}_{nb}}-{{\bf{F}}_{nw}}}\right),\\\
\left\|{{\varepsilon_{\min}}}\right\|<\left\|{{\varepsilon_{t}}}\right\|<\left\|{{\varepsilon_{\max}}}\right\|,\end{array}$
(19)

where $\lambda$ is a random number between 0 and 1, and ${\varepsilon_{t}}$
denotes the distance that the frog moves, and $t=1,2,\ldots,Q$ in which $Q$ is
the number of local iteration. ${\varepsilon_{\min}}$ and
${\varepsilon_{\max}}$ are the minimum and the maximum distances allowing
frogs to move. The individual frog ${F_{nw}}$ position is updated by

${{\bf{F}}_{nw}}={{\bf{F}}_{nw}}+{\varepsilon_{t}}.$ (20)

If the updated frog is superior to the original frog in $\varphi$, the frog
${F_{nw}}$ of the original meme group can be replaced. Otherwise,
${{\bf{F}}_{g}}$ is used instead of ${{\bf{F}}_{nb}}$ to perform local
location update, as given by

$\begin{array}[]{l}{\varepsilon_{\rm{t}}}=\lambda\cdot\left({{{\bf{F}}_{g}}-{{\bf{F}}_{nw}}}\right),\\\
\left\|{{\varepsilon_{\min}}}\right\|<\left\|{{\varepsilon_{t}}}\right\|<\left\|{{\varepsilon_{\max}}}\right\|,\end{array}$
(21) ${{\bf{F}}_{nw}}={{\bf{F}}_{nw}}+{\varepsilon_{t}}.$ (22)

If the above steps still do not achieve a satisfactory fit, a new frog is
randomly generated to replace the original ${{\bf{F}}_{nw}}$. When the local
deep search of all subpopulations is completed, all frog individuals are
remixed and sorted and the subpopulation is divided again. Then the local
depth search is carried out again and the process is repeated until
convergence conditions or the maximum number of iterations are met. The final
output global optimal value ${{\bf{F}}_{g}}$ is the parameter
${{\bf{\theta}}^{\bf{*}}}$ we estimate.

###  3.4 Propagation Model Self-Learning and Spectrum Reconstruction

In this section, we propose a model self-learning method to learn the
electromagnetic environment. Considered that fading under the urban scenario
can result in a loss of signal power without reducing the power of noise, it
will cause performance degradation of a communication system. It can be solved
with a semi-empirical propagation PL model to achieve an accurate estimate of
the actual scenarios [36].

The RSS at the $i$-th point from the $j$-th radiation source can be modeled as

${P_{ij}}^{r}={P_{j}}^{t}-{L_{ij}},$ (23)

where ${P_{ij}}^{r}$ and ${P_{j}}^{t}$ are the RSS at the $i$-th grid from the
$j$-th radiation source and the transmitting power of the $j$-th radiation
source, respectively. ${L_{ij}}$ is the PL between the $i$-th point and the
$j$-th radiation source. RSS depends on the PL model and the transmitting
power of all radiation sources in the measurement area. Then, the RSS
${P_{i}}^{r}$ of a certain point can be expressed as

$\begin{array}[]{l}{P_{ij}}^{r}\left[{{\rm{mW}}}\right]={10^{{P_{ij}}^{r}\left[{{\rm{dBm}}}\right]/10}},\\\
{P_{i}}^{r}\left[{{\rm{mW}}}\right]=\sum\limits_{i=1}^{K}{{P_{ij}}^{r}\left[{{\rm{mW}}}\right]},\end{array}$
(24)

where ${P_{i}}^{r}$ is the RSS value at the $i$-th point from all radiation
sources.

If the semi-empirical propagation model is chosen properly, the awareness of
the environment can be accurate by modifying the semi-empirical propagation
model. It should be mentioned that for different scenarios, such as urban,
suburban, indoor, mountain, etc, we can choose different semi-empirical
models. In this paper, we choose an urban scenario in which the
electromagnetic propagation characteristics of 3D space and the influence of
shadow fading are considered. Besides, we reduce the impact of small-scale
fading by averaging multiple RSS samples over time. Then the channel model is
modeled according to the 3D airspace scenario, as given by [37]

$\begin{array}[]{l}L=32.4+20\log 10({f_{\rm{c}}})\\\ +10(A\cdot{(h)^{B}})\log
10(d)+{\chi_{\sigma}},\end{array}$ (25)

where $A$ and $B$ are environment-dependent parameters and they vary
extraordinary in different scenarios. $H$ is the height relative to the
ground. ${\chi_{\sigma}}$ is a zero-mean Gaussian variable with variance of
$\sigma$ which represents the log-normal shadow fading.

According to (3) and the data set we sampled, the parameters $A$, $B$ and
$\sigma$ are estimated by solving a nonlinear fitting problem. For the $i$-th
data, the actual measured path loss is

${L_{i}}^{{\rm{real}}}=\left({10\cdot\log
10\left({\sum\limits_{j=1}^{K}{{P_{j}}^{t}}}\right)}\right)-{P_{i}}^{r},$ (26)

where ${P_{j}}$ and ${P_{i}}$ represents the transmitting power of the $j$-th
radiation source and the RSS value of the $i$-th point.

Theoretical path loss can be obtained by (25)

$\begin{array}[]{l}{L_{i}}^{{\rm{theo}}}=\sum\limits_{j=1}^{K}{{L_{ji}}}\\\
=\sum\limits_{j=1}^{K}{\left({32.4+20\lg{f_{c}}+10(A\cdot{{(h)}^{B}})\lg{d_{ij}}+{\chi_{\sigma}}}\right)},\end{array}$
(27)

where ${d_{ij}}$ is the distance between the $i$-th point and the $j$-th
radiation source. We define ${\bf{\Theta}}=\left[{A,B,\sigma}\right]$ to be
the unknown parameters of PL model. We obtain the estimated parameters
${\bf{\hat{\Theta}}}$ by constructing an optimization function, given by

${\bf{\hat{\Theta}}}=\arg\min{\left\|{{L_{i}}^{{\rm{theo}}}-{L_{i}}^{{\rm{real}}}}\right\|_{2}},$
(28)

where ${\bf{\hat{\Theta}}}=\left[{\hat{A},\hat{B},\hat{\sigma}}\right]$ is the
parameter to be estimated. There are only three parameters to be estimated
here, so we can solve this problem by selecting the appropriate optimization
algorithm.

After solving the ${\bf{\hat{\Theta}}}$, we can utilize ${\bf{\hat{\Theta}}}$
and ${{\bf{\theta}}^{\bf{*}}}$ to construct the 3D spectrum map. According to
(23) and (24), the RSS at any location can be evaluated by the parameters of
radiation sources and the modified propagation model.

##  4 SIMULATION RESULTS AND VALIDATIONS

###  4.1 Experiment Setup

In this chapter, the performance of the proposed 3D spectrum map construction
method is analyzed by simulation data. The ROI is a typical campus scenario.
There are 9 buildings with the heights from 19 m to 35 m, and most of them are
about 20 m. The scope of the entire region is $500{\kern 1.0pt}{\rm{m}}\times
500{\kern 1.0pt}{\rm{m}}\times 100{\kern 1.0pt}{\rm{m}}$, which means there
are about $7.7854\times 10^{-5}$ buildings per square meter. We discretize the
space into a spectrum tensor of ${N_{1}}\times{N_{2}}\times{N_{3}}=100\times
100\times 10$ with interval size of ${\rm{5{\kern 1.0pt}m}}\times{\rm{5{\kern
1.0pt}m}}\times{\rm{10{\kern 1.0pt}m}}$. It should be mentioned that it is
impossible to acquire all spectrum data of the ROI in a real world. However,
it is necessary for the performance evaluation of the spectrum map
construction methods to acquire all spectrum data of the measurement area. In
this case, we leverage the RT technique to obtain the 3D spectrum map
simulation results of the ROI used for the performance evaluation. The RT
technique has been widely used for radio propagation modeling, which has good
performance for the small specific area and high frequency. In RT, the
electromagnetic wave is viewed as a bunch of rays. When these rays interact
with the scatterer surface and wedge, they can be reflected and diffracted.
After all rays are tracked with the forward technique or the reverse
technique, propagation parameters, i.e., electric field intensity or received
signal strength of each ray, can be acquired. In this article, we only
concentrate on the RSS. The simulation parameters of radiation sources are
shown in Table1. The values of other parameters are presented in Table2.

We compare the results of the proposed method based on self-learning
propagation model (SLPM) with that of the proposed method based on free space
propagation model (FSPM). The results of IDW and High Accuracy Low Rank Tensor
Completion (HaLRTC) are also provided. Root mean square error (RMSE), the
Correct Detection Zone Ration (CDZR) and False Alarm Zone Ratios (FAZR) are
selected as performance metrics. We also analysed the impact of radiation
source number $K$. When comparing the performance of the spectrum map
construction methods, we sample the data randomly, so that only a part of
entries is retained in the spectrum tensor. The sampling rate indicates the
percentage of the observed elements in the tensor.

Table 1: The parameters of radiation sources for different number Number | Location (m) |  Transmission power (w)
---|---|---
x | y | z
2 | 310 | -239 | 2 | 1
235 | -105 | 2 | 1
3 | 345 | -365 | 33.77 | 1
205 | -265 | 2 | 1
245 | -95 | 2 | 1
4 | 330 | -370 | 33.77 | 1
400 | -140 | 23.3 | 1
185 | -255 | 2 | 1
245 | -85 | 2 | 1
Table 2: Other simulation parameters Parameter | Value
---|---
Scenario | Urban
| The range of
---
measurement area
| X(m): 0-500
---
Y(m): -450-50
Z(m): 0-100
Frequency | 100 MHz
Antenna type | Omnidirectional

###  4.2 Comparison and Validation

$i)$ Comparison of Received Signal Strength Recovery: In this section, the
RMSE are chosen as the performance evaluation criterion of RSS recovery
average error. The average error is defined as the difference of each grid
point in RSS between the estimated REM and the real REM, given by

${\rm{RMSE}}=\frac{1}{N}\sum\limits_{i=1}^{N}{\sqrt{{{\left|{\frac{{P_{i}^{{}^{r}{\rm{est}}}-P_{i}^{{}^{r}{\rm{real}}}}}{{P_{i}^{{}^{r}{\rm{real}}}}}}\right|}^{2}}}},$
(29)

where $P_{i}^{{}^{r}{\rm{est}}}$ and $P_{i}^{{}^{r}{\rm{real}}}$ are the
estimated value and the actual value respectively.
$N={N_{1}}\times{N_{2}}\times{N_{3}}$ is the total number of grids.

In Figure 3, we compare the results of the proposed method based on SLPM with
that of the proposed method based on FSPM. The results of IDW and HaLRTC are
also provided. Given the sampling rate $r$, the sampling locations are
arbitrarily chosen from the measurement area. It shows that the RMSE of RSS
recovery of all algorithms decrease with the increase of sampling rate. The
proposed method based on SLPM outperforms other three algorithms. The data-
driven method IDW and HaLRTC has poor performance in the complex urban
environment. In addition, the proposed method based on the SLPM can better
guarantee the recovery of the actual electromagnetic environment than the
method directly based on FSPM.

We also compare the influence of the different number $K$ of radiation sources
on the RSS recovery errors of the four algorithms with the same sampling rate.
The results can be seen in Figure 4. We can see that the RMSE of the RSS
recovery error of all algorithms decrease with the decrease of radiation
source number. The recovery performance of the data-driven methods, i.e., IDW
and HaLRTC, is far inferior to the two model-driven methods when multi-
radiation sources exist. This is because the data-driven methods have high
requirements on the number and locations of sampled data. Moreover, the
spectrum recovery performance of the plane with unsampled height is poor but
the sampling rate is high. By comparing the proposed method based on the SLPM
with that based on the FSPM, self-learning of the model can further improve
the recovery performance.

Figures 3 and 4 also show that the construction accuracy of model-driven
spectrum map methods is higher than that of data-driven spectrum map methods
in terms of the sampling rate and the radiation source number. In addition,
although the proposed algorithm based on the FSPM and the SLPM respectively
have similar convergence rates, the recovery accuracy of the latter is high.

Figure 3: The RMSE of RSS recovery error vs. the sampling rate $r$ ($K$=3).
Figure 4:  The RMSE of RSS recovery error vs. the radiation source number $K$
($r$=0.2).

$ii)$ Comparison of CDZR and FAZR: In $i)$, we calculated the RMSE of RSS
recovery results between the real spectrum map and the estimated spectrum map
to quantify the similarity. To determine where spectrum opportunities are lost
and where CR conflicts with ongoing transmissions of radiation sources in the
region, the Correct Detection Zone Ration (CDZR) and False Alarm Zone Ratios
(FAZR) are introduced in [28].

Figure 5: The CDZR/FAZR of RSS recovery error vs. the sampling rate $r$
($K$=3).

The intersection of the zones of true REM and estimated REM defines the zones
that are correctly and incorrectly determined. Regions where both the
estimated spectrum map and the actual spectrum map correctly identify
forbidden is marked as the Correct Detection Zone of type 1
($Z_{1}^{{\rm{C}}{{\rm{D}}_{j}}}$), while identifying permitted is marked as
Correct Detection Zone of type 0 ($Z_{0}^{{\rm{C}}{{\rm{D}}_{j}}}$). False
Alarm Zone $Z_{0}^{{\rm{F}}{{\rm{A}}_{j}}}$ represents the area where
communication is allowed but is incorrectly detected as the prohibitive zone.
Missed Detection Zone $Z_{1}^{{\rm{M}}{{\rm{D}}_{j}}}$ represents the area
where transmission is forbidden but missed detection. The CDZR and FAZR can be
obtained from the following equations

$\begin{array}[]{l}{\rm{FAZ}}{{\rm{R}}_{j}}=\frac{{S\left({Z_{0}^{{\rm{F}}{{\rm{A}}_{j}}}}\right)}}{{S\left({Z_{0}^{{\rm{F}}{{\rm{A}}_{j}}}}\right)+S\left({Z_{0}^{{\rm{C}}{{\rm{D}}_{j}}}}\right)}}\begin{array}[]{*{20}{c}},&{j=1,2,\ldots,K}\end{array},\\\
{\rm{FAZR}}=\sum\limits_{j=1}^{K}{{\rm{FAZ}}{{\rm{R}}_{j}}},\end{array}$ (30)
$\begin{array}[]{l}{\rm{CDZ}}{{\rm{R}}_{j}}=\frac{{S\left({Z_{1}^{{\rm{C}}{{\rm{D}}_{j}}}}\right)}}{{S\left({Z_{1}^{{\rm{C}}{{\rm{D}}_{j}}}}\right)+S\left({Z_{1}^{{\rm{M}}{{\rm{D}}_{j}}}}\right)}}\begin{array}[]{*{20}{c}},&{j=1,2,\ldots,K}\end{array},\\\
{\rm{CDZR}}=\sum\limits_{j=1}^{K}{{\rm{CDZ}}{{\rm{R}}_{j}}},\end{array}$ (31)

where $S\left(Z\right)$ is the size of the zone $Z$.
${\rm{FAZ}}{{\rm{R}}_{j}}$ and ${\rm{CDZ}}{{\rm{R}}_{j}}$ are Correct
Detection Zone Ration and False Alarm Zone Ratios at the $j-th$ radiation
source respectively. Having higher CDZR means disturbing the active user is
less probable. Having less FAZR means white spaces can be utilized better.
Better spectrum map construction performance fit better on the true spectrum
map contour with predetermined threshold value.

We compared the ${\rm{CDZR}}/{\rm{FAZR}}$ curves as shown in Figure 5. The
simulation results show that based on the knowledge of radiation sources and
self-learning of the model, the proposed algorithm can make better use of the
white spaces and reduce the interference to the radiation sources than the
other two algorithms. Figure 5 can also be seen that the data-driven methods
cannot achieve satisfactory performance, since their performance largely
depends on the sampling positions. The data of unsampled position, which is
interpolated by the data-driven methods, is smaller than the maximum value of
the sampled data. Then if the sampling position is far away from the radiation
source, it is difficult to recover its nearby real spectrum data.

Figure 6: The localization errors vs. the sampling rate $r$. Figure 7: The
source signal strength recovery errors vs. the sampling rate $r$. Figure 8:
The RSS recovery errors vs. the sampling rate $r$. Figure 9: The success rate
of radiation source number identification vs. the sampling rate $r$.

(a)

(b)

(c)

(d)

(e)

Figure 10: 3D visualization results of spectrum situation recovery (K=4). (a):
the simulation of spectrum map. (b)-(e): the spectrum situation recovery
results of the proposed method based on the SLPM, FSPM, IDW and Tensor
completion, respectively.

$iii)$ Impact of Radiation Source Number: We study the impact of radiation
source number $K$ on the performance of the proposed algorithm including
localizations, signal strength recovery, the success rate of radiation source
number identification and RSS recovery.

The average localization error is defined as

${\rm{LOC\\_E}}=\frac{1}{K}\sum\limits_{i=1}^{K}{{{\left\|{{\rm{lo}}{{\rm{c}}_{i}}^{{\rm{est}}}-{\rm{lo}}{{\rm{c}}_{i}}^{{\rm{real}}}}\right\|}_{2}}},$
(32)

where ${\rm{lo}}{{\rm{c}}_{i}}^{{\rm{est}}}$ and
${\rm{lo}}{{\rm{c}}_{i}}^{{\rm{real}}}$ denote the estimated and the real
locations of radiation sources in Cartesian coordinate.

The signal strength recovery error is defined as

${\rm{SS\\_E}}=\frac{1}{K}\sum\limits_{i=1}^{K}{{{\left\|{{\rm{si}}{{\rm{g}}_{i}}^{{\rm{est}}}-{\rm{si}}{{\rm{g}}_{i}}^{{\rm{real}}}}\right\|}_{2}}},$
(33)

where ${\rm{si}}{{\rm{g}}_{i}}^{{\rm{est}}}$ and
${\rm{si}}{{\rm{g}}_{i}}^{{\rm{real}}}$ represent estimated and true signal
strengths of radiation sources.

In Figures 6 and 7, the impact of $K$ on the performances of localizations and
signal strength recovery is analyzed. It can be seen that the performances of
three algorithms are all improved as $K$ decreases. The localization error and
signal strength recovery error decrease when the number of radiation sources
is reduced. These factors are the key to determine the performance of the
proposed algorithm. The RMSE of RSS recovery also decreases with the increase
of $K$, as shown in Figure 8. This is because increasing $K$ will aggravate
the computation and lead to more interactions between radiation sources. In
addition, when the sampling rate $r$ grows, the performances of these four
algorithms improve significantly.

Furthermore, we analyze the success rate of radiation source number
identification as shown in Figures 9. It can be seen that when the sampling
rate is over 0.3, our MMPLD clustering algorithm can identify the radiation
source number with the success rate of 0.7.

$iv)$ 3D Spectrum Map Visualization: Figures 10 and 11 present the 3D spectrum
map visualization results of the IDW, the tensor completion, the proposed
method based on the SLPM and FSPM. In order to analyze the construction
situation of spectrum map clearly, we observe two figures from different
angles.

Figure 10 shows the 3D visualization of spectrum situation recovery. By
comparing Figure 10 (b)-(e), we can see that the IDW and tensor completion
perform poorly in 3D spectrum situation recovery, as their performance only
depends on the data. Figure 10 (b) and (c) can be seen that the proposed
method based on the SLPM can recover the original spectrum situation veritably
with the consideration of shadow fading in the urban environment.

Figure 11 (b)-(e) further visually reveal that our proposed method performs
well in the spectrum situation recovery at radiation sources and in the source
signal strength recovery of radiation sources. As shown in Figure 11 (b) and
(c), by exploiting the knowledge of radiation sources, the spectrum situation
can be well recovered around the radiation sources. Besides, considering the
self-learning of the propagation model, the visualization in Figure 11(b) can
capture the shadowing, which is superior than the counterpart in Figure 11
(c). The IDW and tensor completion in Figure 11 (d) and (e) failed to recover
the spectrum data around the radiation sources due to the influence of
sampling positions . It can be seen that the construction performance of the
tensor completion method is very poor when the sampling data are relatively
unevenly distributed in the ROI.

(a)

(b)

(c)

(d)

(e)

Figure 11: Visualization results of spectrum situation recovery at radiation
sources (K=4). (a): the simulation of spectrum map. (b)-(e): the spectrum
situation recovery at radiation sources of the proposed method based on the
SLPM, FSPM, IDW and Tensor completion, respectively.

##  5 CONCLUSION

In this paper, we have investigated 3D spectrum mapping under multi-radiation
source urban scenarios. We have constructed the spectrum map indirectly
through the study of the parameters of radiation sources and the environment
of the 3D space. We have proposed the 3D spectrum mapping and reconstruction
scheme, which is composed of an extraction of radiation source knowledge and
self-learning of the propagation model. We have designed a MMPLD-based
clustering algorithm to detect the number of radiation sources and a SFLA-
based optimizing estimation algorithm to locate the radiation sources.
Propagation model self-learning based on sampling data has also been proposed,
which can describe the real time electromagnetic propagation states. We have
developed the empirical PL model by sampling data to make the proposed
algorithm self learn. We have compared the RMSE of RSS recovery and CDZR/FAZR
curves among four algorithms i.e., the proposed method based on the SLPM,
FSPM, IDW and tensor completion. Furthermore, the impact of radiation source
number has also been studied. Simulation results have shown that the proposed
method has obvious advantages in the RSS recovery under multi-radiation source
urban scenarios. Further research will be conducted to sampling optimization
and dynamic reconstruction of radiation source movement.

## ACKNOWLEDGEMENT

This work was supported in part by the National Key Scientific Instrument and
Equipment Development Project under Grant No. 61827801, in part by the open
research fund of State Key Laboratory of Integrated Services Networks, No.
ISN22-11, in part by Natural Science Foundation of Jiangsu Province, No.
BK20211182, and in part by the open research fund of National Mobile
Communications Research Laboratory, Southeast University, No. 2022D04.

## References

  * [1] S. Haykin, “Cognitive radio: brain-empowered wireless communications,” _IEEE journal on selected areas in communications_ , vol. 23, no. 2, pp. 201–220, 2005.
  * [2] X. Zhang, L. Xu, L. Xu, and D. Xu, “Direction of departure (dod) and direction of arrival (doa) estimation in mimo radar with reduced-dimension music,” _IEEE communications letters_ , vol. 14, no. 12, pp. 1161–1163, 2010.
  * [3] G. Ding, Q. Wu, L. Zhang, Y. Lin, T. A. Tsiftsis, and Y.-D. Yao, “An amateur drone surveillance system based on the cognitive internet of things,” _IEEE Communications Magazine_ , vol. 56, no. 1, pp. 29–35, 2018.
  * [4] B. Wang and K. R. Liu, “Advances in cognitive radio networks: A survey,” _IEEE Journal of selected topics in signal processing_ , vol. 5, no. 1, pp. 5–23, 2010.
  * [5] Y. Zhang, L. Zhao, and Z. He, “A 3-d hybrid dynamic channel model for indoor thz communications,” _China Communications_ , vol. 18, pp. 50–65, 2021.
  * [6] X. Zhou, J. Ma, G. Y. Li, Y. H. Kwon, and A. C. K. Soong, “Probability-based combination for cooperative spectrum sensing,” _IEEE Transactions on Communications_ , vol. 58, no. 2, pp. 463–466, 2010.
  * [7] X. Liu, C. Sun, M. Zhou, C. Wu, B. Peng, and P. Li, “Reinforcement learning-based multislot double-threshold spectrum sensing with bayesian fusion for industrial big spectrum data,” _IEEE Transactions on Industrial Informatics_ , vol. 17, no. 5, pp. 3391–3400, 2021.
  * [8] X. Liu and X. Zhang, “Noma-based resource allocation for cluster-based cognitive industrial internet of things,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 8, pp. 5379–5388, 2020.
  * [9] J. Wang, M. Ghosh, and K. Challapali, “Emerging cognitive radio applications: A survey,” _IEEE Communications Magazine_ , vol. 49, no. 3, pp. 74–81, 2011\.
  * [10] Z. Feng, C. Qiu, Z. Feng, Z. Wei, W. Li, and P. Zhang, “An effective approach to 5g: Wireless network virtualization,” _IEEE Communications Magazine_ , vol. 53, pp. 53–59, 2015.
  * [11] Q. hui Wu, G. Ding, J. Wang, and Y. dong Yao, “Spatial-temporal opportunity detection for spectrum-heterogeneous cognitive radio networks: Two-dimensional sensing,” _IEEE Transactions on Wireless Communications_ , vol. 12, pp. 516–526, 2013.
  * [12] Y. Chen, G. Yu, Z. Zhang, H.-h. Chen, and P. Qiu, “On cognitive radio networks with opportunistic power control strategies in fading channels,” _IEEE Transactions on Wireless Communications_ , vol. 7, no. 7, pp. 2752–2761, 2008\.
  * [13] Z. Quan, S. Cui, A. H. Sayed, and H. V. Poor, “Optimal multiband joint detection for spectrum sensing in cognitive radio networks,” _IEEE Transactions on Signal Processing_ , vol. 57, pp. 1128–1140, 2009.
  * [14] Y. Zhao, L. Morales, J. Gaeddert, K. K. Bae, J.-S. Um, and J. H. Reed, “Applying radio environment maps to cognitive wireless regional area networks,” in _2007 2nd IEEE International Symposium on New Frontiers in Dynamic Spectrum Access Networks_. IEEE, 2007, pp. 115–118.
  * [15] Y. Hu and R. Zhang, “Differentially-private incentive mechanism for crowdsourced radio environment map construction,” in _IEEE INFOCOM 2019 \- IEEE Conference on Computer Communications_ , 2019, pp. 1594–1602.
  * [16] Y. Xu, A. Anpalagan, Q. Wu, L. Shen, Z. Gao, and J. Wang, “Decision-theoretic distributed channel selection for opportunistic spectrum access: Strategies, challenges and solutions,” _IEEE Communications Surveys & Tutorials_, vol. 15, no. 4, pp. 1689–1713, 2013.
  * [17] Z. Lin, T. Lv, W. Ni, J. A. Zhang, and R. P. Liu, “Tensor-based multi-dimensional wideband channel estimation for mmwave hybrid cylindrical arrays,” _IEEE Transactions on Communications_ , vol. 68, pp. 7608–7622, 2020.
  * [18] Q. Zhao and B. M. Sadler, “A survey of dynamic spectrum access,” _IEEE Signal Processing Magazine_ , vol. 24, pp. 79–89, 2007.
  * [19] H. Chen, M. Zhou, L. Xie, and J. Li, “Cooperative spectrum sensing with m-ary quantized data in cognitive radio networks under ssdf attacks,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 8, pp. 5244–5257, 2017\.
  * [20] V. Atanasovski, J. van de Beek, A. Dejonghe, D. Denkovski, L. Gavrilovska, S. Grimoud, P. Mähönen, M. Pavloski, V. Rakovic, J. Riihijarvi _et al._ , “Constructing radio environment maps with heterogeneous spectrum sensors,” in _2011 IEEE International Symposium on Dynamic Spectrum Access Networks (DySPAN)_. IEEE, 2011, pp. 660–661.
  * [21] J. Talvitie, M. Renfors, and E. S. Lohan, “Distance-based interpolation and extrapolation methods for rss-based localization with indoor wireless signals,” _IEEE transactions on vehicular technology_ , vol. 64, no. 4, pp. 1340–1353, 2015.
  * [22] X. Han, L. Xue, F. Shao, and Y. Xu, “A power spectrum maps estimation algorithm based on generative adversarial networks for underlay cognitive radio networks,” _Sensors_ , vol. 20, no. 1, p. 311, 2020.
  * [23] F. Shen, G. Ding, Z. Wang, and Q. Wu, “Uav-based 3d spectrum sensing in spectrum-heterogeneous networks,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 6, pp. 5711–5722, 2019.
  * [24] G. Y. Lu and D. W. Wong, “An adaptive inverse-distance weighting spatial interpolation technique,” vol. 34, no. 9. Elsevier, 2008, pp. 1044–1055.
  * [25] Y. Hu and R. Zhang, “A spatiotemporal approach for secure crowdsourced radio environment map construction,” _IEEE/ACM Transactions on Networking_ , vol. 28, no. 4, pp. 1790–1803, 2020.
  * [26] E. Dall’Anese, S.-J. Kim, and G. B. Giannakis, “Channel gain map tracking via distributed kriging,” _IEEE transactions on vehicular technology_ , vol. 60, no. 3, pp. 1205–1211, 2011.
  * [27] M. Tang, G. Ding, Q. Wu, Z. Xue, and T. A. Tsiftsis, “A joint tensor completion and prediction scheme for multi-dimensional spectrum map construction,” _IEEE Access_ , vol. 4, pp. 8044–8052, 2016.
  * [28] H. B. Yilmaz and T. Tugcu, “Location estimation-based radio environment map construction in fading channels,” _Wireless communications and mobile computing_ , vol. 15, no. 3, pp. 561–570, 2015.
  * [29] M. Pesko, T. Javornik, L. Vidmar, A. Košir, M. Štular, and M. Mohorčič, “The indirect self-tuning method for constructing radio environment map using omnidirectional or directional transmitter antenna,” _EURASIP Journal on wireless communications and networking_ , vol. 2015, no. 1, pp. 1–12, 2015.
  * [30] J. A. Bazerque and G. B. Giannakis, “Distributed spectrum sensing for cognitive radio networks by exploiting sparsity,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 3, pp. 1847–1862, 2009.
  * [31] B. J. Frey and D. Dueck, “Clustering by passing messages between data points,” _science_ , vol. 315, no. 5814, pp. 972–976, 2007.
  * [32] K. P. Sinaga and M.-S. Yang, “Unsupervised k-means clustering algorithm,” _IEEE Access_ , vol. 8, pp. 80 716–80 727, 2020.
  * [33] H. Xiong, J. Wu, and J. Chen, “K-means clustering versus validation measures: A data-distribution perspective,” _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , vol. 39, no. 2, pp. 318–331, 2009.
  * [34] I. Katsavounidis, C.-C. J. Kuo, and Z. Zhang, “A new initialization technique for generalized lloyd iteration,” _IEEE Signal processing letters_ , vol. 1, no. 10, pp. 144–146, 1994.
  * [35] S. Mahdavi, M. E. Shiri, and S. Rahnamayan, “Metaheuristics in large-scale global continues optimization: A survey,” _Information Sciences_ , vol. 295, pp. 407–428, 2015.
  * [36] Q. Zhu, Z. Zhao, K. Mao, X. Chen, W. Liu, and Q. Wu, “A real-time hardware emulator for 3d non-stationary u2v channels,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 68, no. 9, pp. 3951–3964, 2021\.
  * [37] Q. Zhu, M. Yao, F. Bai, X. Chen, W. Zhong, B. Hua, and X. Ye, “A general altitude-dependent path loss model for uav-to-ground millimeter-wave communications,” _Frontiers of Information Technology & Electronic Engineering_, vol. 22, no. 6, pp. 767–776, 2021.

wangjie.epsJie Wang received the B.S. degree in internet of things engineering
from the College of Information Science and Technology, Nanjing Forestry
University of China, Nanjing, China, in 2021. She is currently pursuing the
Ph.D. degree in communications and information systems with the College of
Electronic and Information Engineering, Nanjing University of Aeronautics and
Astronautics. Her current research interests is spectrum mapping.
ZhiPengLin.epsZhipeng Lin received the Ph.D. degrees from the School of
Information and Communication Engineering, Beijing University of Posts and
Telecommunications, Beijing, China, and the School of Electrical and Data
Engineering, University of Technology of Sydney, NSW, Australia, in 2021. He
is currently an Associate Researcher in the College of Electronic and
Information Engineering, Nanjing University of Aeronautics and Astronautics,
Nanjing, China. His current research interests include signal processing,
massive MIMO, spectrum sensing, and UAV communications. zhuqiuming.epsQiuming
Zhu received the B.S. degree in electronic engineering and the M.S. and Ph.D.
degrees in communication and information system from the Nanjing University of
Aeronautics and Astronautics (NUAA), Nanjing, China, in 2002, 2005, and 2012,
respectively. He is currently a Professor in the College of Electronic and
Information Engineering, Nanjing University of Aeronautics and Astronautics,
Nanjing, China. His current research interests include channel sounding,
modeling, and emulation for the fifth/sixth generation (5G/6G) mobile
communication, vehicle-to-vehicle (V2V) communication, and unmanned aerial
vehicles (UAV) communication systems. wuqihui.epsQihui Wu received the B.S.
degree in communications engineering and the M.S. and Ph.D. degrees in
communications and information system from the PLA University of Science and
Technology, Nanjing, China, in 1994, 1997, and 2000, respectively. He is
currently a Professor with the College of Electronic and Information
Engineering, Nanjing University of Aeronautics and Astronautics. His current
research interests include algorithms and optimization for cognitive wireless
networks, soft-defined radio, and wireless communication systems.
lantianxu.epsTianxu Lan received the B.S. degree in information engineering
from Nanjing University of Aeronautics and Astronautics (NUAA) in 2019. He is
currently working towards the master degree in electronic information
engineering, NUAA. His current research is reconstruction of spectrum
situation. zhaoyi.epsYi Zhao received the B.S. degree in information
engineering from Nanjing University of Aeronautics and Astronautics (NUAA) in
2021. He is currently working towards the master degree in electronic
information engineering, NUAA. His current research is temporal and spatial
prediction and reconstruction of spectrum situation. YunPengBai.epsYunpeng Bai
received the B.E. degree in electrical engineering and automation from
Shandong University (SDU) in 2019. He is currently working towards the master
degree in electronic information engineering, NUAA. His current research is
UAV path planning. zhongweizhi.epsWeizhi Zhong received her B.S. and M.S.
degrees in communication and information system from the Jilin University, and
received the Ph.D. degree in communication and information system from the
Harbin Institute of Technology, respectively. She is an associate professor in
the College of Astronautics, Nanjing University of Aeronautics and
Astronautics (NUAA). Her research interests include millimeter wave
communication for the fifth generation (5G), massive MIMO technique and
beamforming and beam tracking technique.


# Basis set calculations of heavy atoms

M. G. Kozlov1,2 Yu. A. Demidov1,2 M. Y. Kaygorodov1 E. V. Triapitsyna1 1
Petersburg Nuclear Physics Institute NRC ``Kurchatov Institute'', Gatchina
188300, Russia
2 St. Petersburg Electrotechnical University ``LETI'', St. Petersburg 197376,
Russia

###### Abstract

Most modern calculations of many-electron atoms use basis sets of atomic
orbitals. An accurate account for the electronic correlations in heavy atoms
is very difficult computational problem and optimization of the basis sets can
reduce computational costs and increase final accuracy. Here we suggest a
simple differential ansatz to form virtual orbitals from the Dirac-Fock
orbitals of the core and valence electrons. We use basis sets with such
orbitals to calculate different properties in Cs including hyperfine structure
constants and QED corrections to the valence energies and to the $E1$
transition amplitudes.

##  I Introduction

Many efficient methods of atomic calculations, such as configuration
interaction (CI), many-body perturbation theory (MBPT), coupled cluster (CC),
or their combinations, employ basis sets of one-electron orbitals. For
example, the methods CI+MBPT and CI+CC Dzuba et al. (1996); Kozlov (2004);
Dzuba (2005); Safronova et al. (2009); Kahl and Berengut (2019) use CI for
valence electrons, where the electronic-correlations effects are strong, and
account for weaker core-valence correlations by means of either the second
order MBPT, or (linearized) CC method. In these methods, the relativistic
effects are treated within the no-pair approximation for the Dirac-Coulomb-
Breit Hamiltonian. On the other hand, QED effects may be approximately
accounted for using the model-QED-potential approach Pyykkö and Zhao (2003);
Flambaum and Ginges (2005); Tupitsyn and Berseneva (2013). For heavy
polyvalent atoms such calculations become very computationally expensive
Cheung et al. (2021). That is why it is very important to develop efficient
basis sets, which provide high accuracy at reasonable length.

The Dirac-Fock (DF) method serves as initial approximation for most of the
atomic calculations. However in practice, it is useful to keep DF orbitals in
the basis set without re-expanding them in some nonphysical basic orbitals. On
the other hand, the highly excited DF orbitals are usually ineffective in
accounting for the correlations between valence electrons, as their radius
grows too rapidly. Because of that the B-splines Johnson et al. (1988), Sturm
orbitals Tupitsyn and Loginov (2003), or other simple orbitals Bogdanovich and
Žukauskas (1983); Kozlov et al. (1996) usually turn out to be more useful.
Therefore, an efficient basis set has to include different types of orbitals.
An effective method to merge two subsets of orbitals in one basis set was
suggested in Ref. Kozlov and Tupitsyn (2019). The first subset consisted of DF
orbitals for core and valence shells and the second subset consisted of
B-splines, whose parameters were optimized to complement the first subset.
Such mixed basis sets were tested for calculations of the valence energies of
Au and Fr and sufficiently fast convergence was observed: saturation was
reached for 25 splines per partial wave Kozlov and Tupitsyn (2019).

It is known, that accurate calculations of the hyperfine structure and parity
non-conservation effects for heavy atoms are much more challenging then
calculations of the energies, or transition amplitudes. These properties
depend on the wave function at small distances. Correlations change the
behaviour of the valence-electrons wave functions in this region and lead to
the large corrections to the matrix elements of these operators. Because of
that the convergence with respect to the number of orbitals in the basis set
is typically much slower. Test calculations of the magnetic hyperfine constant
$A$ for Au did not converge even for the mixed basis set with 45 B-splines per
partial wave Kozlov and Tupitsyn (2019). It was found that small components of
the valence orbitals near the origin were not smooth. Increasing the size of
the basis set decreased the amplitude of these non-physical dips or bumps, but
shifted them closer to the origin, where the hyperfine interaction was
particularly strong. The solution suggested in that paper was to add to the
basis set DF orbitals for the ion AuM+, or, in other words, the orbitals for
the $V^{N-M}$ potential, where N is the number of electrons in the neutral
atom. Such orbitals are similar to those of the neutral atom, but are more
contracted. That allowed to smoothly change the orbitals near the nucleus. The
calculation of the hyperfine constants $A$ for Au already converged for two
such ionic orbitals and 24 B-splines per partial wave. To make this recipe
work one need to choose $M$ such that ionic orbitals are neither too
different, nor too similar to those of the neutral atom. In the former case
such orbitals become less useful, while in the latter one can run into the
linear dependency problem. In Ref. Kozlov and Tupitsyn (2019) the orbitals for
the $V^{N-6}$ potential were used. Here, we suggest another method to form
additional orbitals to supplement DF orbitals and B-splines. The new procedure
is more formal and does not require arbitrary adjustments. We test these basis
sets for calculation of the QED corrections and hyperfine constants in neutral
Cs.

##  II Method

The basis set for atomic calculations consists of subsets for different
partial waves, which are defined by the relativistic quantum number
$\varkappa=(l-j)(2j+1)$, where $l$ and $j$ are the orbital and total angular
momenta. The relativistic two-component radial wave functions have the form

$\displaystyle\psi_{n\varkappa}(r)$
$\displaystyle=\frac{1}{r}\left(\begin{array}[]{c}P_{n\varkappa}(r)\\\
Q_{n\varkappa}(r)\end{array}\right)\,.$ (3)

Most atomic calculations are done in the no-pair approximation. QED
corrections may be included approximately by means of the model operators
Flambaum and Ginges (2005); Shabaev et al. (2013); Tupitsyn and Berseneva
(2013). This means that the Hamiltonian is projected on the positive energy
electronic states and the basis set does not include negative energy
continuum. To this end, the small components $Q$ of the virtual basis orbitals
can not be chosen arbitrary. Most often they are found from the kinetic
balance condition Stanton and Havriliak (1984):

$\displaystyle Q_{n\varkappa}(r)$
$\displaystyle=-\frac{\alpha}{2}\left(\frac{d}{dr}+\frac{\varkappa}{r}\right)P_{n\varkappa}(r)\,,$
(4)

where $\alpha\approx\tfrac{1}{137}$ is the fine structure constant (we use
atomic units throughout the paper). Below we assume that the small components
are always formed with the help of Eq. (4) and focus on the large components
$P$.

When some small perturbation is added to the Hamiltonian, the large component
$P$ of the electronic orbital slightly changes. These changes do not affect
the asymptotic behaviour at large and small distances and the number of nodes
of the function $P$ remains the same. To account for such changes we need to
add basic functions with similar properties. Let us consider a simple
stretching transformation:

$\displaystyle P(r)$ $\displaystyle\to\tilde{P}(r)=P(kr)\,,$ (5)

where $k=1+\varepsilon$, $|\varepsilon|\ll 1$. Expanding $\tilde{P}$ in powers
of $\varepsilon$ we get:

$\displaystyle\tilde{P}(r)$
$\displaystyle=P(kr)=P\left((1+\varepsilon)r\right)\approx P(r)+\varepsilon
r\frac{dP}{dr}\,.$ (6)

Thus, if we want to guarantee that stretched DF orbital
$\tilde{P}_{n\varkappa}$ can be accurately expanded in the basis set, we can
add a virtual orbital:

$\displaystyle P_{n^{\prime}\varkappa}(r)$
$\displaystyle=r\frac{dP_{n\varkappa}}{dr}\,.$ (7)

Note that this function has the same asymptotic behavior at $r\to 0$ and at
$r\to\infty$ as $P_{n\varkappa}$, however it has an extra node.

Figure 1: Comparison of the differential ansatz (7) with the difference between DF orbitals for the Dirac-Coulomb Hamiltonian with and without model QED potential. Table 1:  Comparison of the basis sets B7 and B30 for Cs. Results for the QED corrections to the valence binding energies and reduced E1 transition amplitudes and for the RPA corrections to the hyperfine constants $A$. Column three gives DF values. Next two columns give QED or RPA corrections calculated on the basis set B7 with 7 virtual orbitals per partial wave and basis B30 with 30 virtual orbitals per partial wave. Last column gives the difference between two calculations in percent.  Property | Units | DF value | QED correction
---|---|---|---
|  |  | Basis B7 | Basis B30 | Diff.
$\varepsilon(6s_{1/2})$ | (cm-1) | 27954.1 | $-15.798$ | $-15.782$ | $-0.1\%$
$\varepsilon(6p_{1/2})$ |  | 18790.5 | $0.814$ | $0.805$ | $-1.1\%$
$\varepsilon(6p_{3/2})$ |  | 18388.8 | $0.169$ | $0.165$ | $-2.6\%$
$\varepsilon(5d_{3/2})$ |  | 14138.5 | $2.177$ | $2.151$ | $-1.2\%$
$\varepsilon(5d_{5/2})$ |  | 14162.6 | $2.063$ | $2.032$ | $-1.5\%$
$E1(6s_{1/2}-6p_{1/2})$ | (a.u.) | 5.27769 | $0.00322$ | $0.00322$ | $0.0\%$
$E1(6s_{1/2}-6p_{3/2})$ |  | 7.42643 | $0.00476$ | $0.00475$ | $-0.2\%$
$E1(6p_{1/2}-5d_{3/2})$ |  | 8.97833 | $-0.00159$ | $-0.00175$ | $10.1\%$
$E1(6p_{3/2}-5d_{3/2})$ |  | 4.06246 | $-0.00069$ | $-0.00076$ | $10.1\%$
$E1(6p_{3/2}-5d_{5/2})$ |  | 12.18640 | $-0.00180$ | $-0.00190$ | $5.6\%$
|  | DF value | RPA correction
$A(6s_{1/2})$ | (MHz) | 1423.3 | $297.6$ | $293.3$ | $-1.4\%$
$A(6p_{1/2})$ |  | 160.88 | $40.02$ | $40.56$ | $1.3\%$
$A(6p_{3/2})$ |  | 23.916 | $19.12$ | $18.91$ | $-1.1\%$

##  III Test calculations for Cs

To illustrate usefulness of the ansatz (7) let us consider QED corrections in
neutral Cs. Firstly, we used a sufficiently large basis set to diagonalize the
Dirac-Fock operator together with the model QED potential of Ref. Shabaev et
al. (2013, 2015). After that we calculated the difference between large
components with and without QED corrections and compared them with the
respective scaled orbitals (7). Figure 1 shows that the latter very much
resemble the former. Thus, we can expect that adding orbitals (7) to the basis
set can significantly improve its quality and speeds up the convergence.

We formed the basis set for Cs which included DF orbitals $1s-7s$, $2p-7p$,
and $3d-6d$. For each DF orbital we added respective virtual orbital (7). To
have the same number of virtual orbitals per partial wave, we formed one
additional $p$ orbital per partial wave and three additional $d$ orbitals per
partial wave using the method of Ref. Kozlov et al. (1996). Thus, this basis
set included 7 virtual orbitals per partial wave. This calculation was
compared with calculation on the long basis set, which included these 7 plus
23 additional virtual orbitals per partial wave. Below, we refer to these
basis sets as B7 and B30, respectively. We see that the QED correction for the
valence $6s$ orbital is practically the same for both basis sets. For the
orbital $1s$, both the basis sets produce some non-physical oscillations, but
their amplitude for the large basis set is few times smaller.

Figure 2: Comparison of the QED calculations with the short (B7) and the long
(B30) basis sets. The short basis set includes 7 orbitals per partial wave. In
the $s$ wave all of them are formed using ansatz (7), while in other partial
waves additional orbitals were added using method Kozlov et al. (1996).

In the calculations of neutral atoms we usually explicitly need only wave
functions of the valence electrons, where our short basis set seems to work
nicely. To test it further we calculated QED corrections to the valence
energies and to the E1 transition amplitudes. In both cases we used model QED
potential Shabaev et al. (2013). This model potential accounts for the QED
corrections to the wave functions. For E1 amplitudes there are also QED
corrections to the vertex, which are not included here. They were calculated
in Refs. Sapirstein et al. (2004); Sapirstein and Cheng (2005) and found to be
much smaller. Results obtained with the basis sets B7 and B30 are compared in
Table 1. We see that two calculations give practically identical QED
correction to the binding energy of the $6s$ electron. For the other partial
waves the corrections to the energies are much smaller and the difference
between results is about 1 – 2%. For the $s-p$ transition amplitudes the
agreement is again very good, but for $p-d$ transitions the difference between
two basis sets is about 10%. On one hand, these corrections are smaller. On
the other hand, in $d$ wave we have only four orbitals formed with the ansatz
(7), while other three orbitals of the basis set B7 are apparently less
useful. We mentioned earlier that convergence in calculations of the hyperfine
constants $A$ is often very slow. Because of that we compared random-phase
approximation (RPA) corrections for these constants obtained with basis sets
B7 and B30. For the orbitals $6s$, $6p_{1/2}$, and $6p_{3/2}$ the difference
was around 1%.

Results presented in Table 2 show that the basis sets, which include virtual
orbitals (7) provide fast convergence for calculations of the QED and RPA
corrections, which are particularly sensitive to the quality of the basis set.
Below we compare QED corrections calculated with these basis sets with
calculations of other groups and study potential sources of the differences.
We focus primarily on the QED corrections to the $E1$ transition amplitudes.

###  III.1 QED corrections to the E1 amplitudes in Cs

QED corrections to the $E1$ amplitudes in Cs were calculated several times.
Sapirstein and Cheng (2005) calculated $6s_{1/2}-6p_{1/2}$ amplitude within
one-determinant approximation, while calculations Fairhall et al. (2023); Tran
Tan and Derevianko (2023); Roberts et al. (2023) included correlations, but
neglected the vertex correction. Traditionally QED corrections to the $E1$
amplitude are expressed in terms of the parameter $R$:

$\displaystyle E1$ $\displaystyle=E1_{0}\left(1+\frac{\alpha}{\pi}R\right)\,,$
(8)

where $E1_{0}$ is the amplitude without QED correction. There are two
contributions to the QED correction $R$:

$\displaystyle R$ $\displaystyle=R_{\mathrm{PO}}+R_{\mathrm{V}}\,,$ (9)

where $R_{\mathrm{PO}}$ accounts for the correction to the orbitals (perturbed
orbitals) and $R_{\mathrm{V}}$ is the vertex correction. The former is
supposed to be larger than the latter Flambaum and Ginges (2005) and it can be
approximated by model QED potentials. At present, the vertex correction can be
calculated only within the consistent QED approach.

Table 2: Comparison of the PO QED corrections $R_{\mathrm{PO}}$ (see Eq. (9)) to the $E1$ transition amplitudes in Cs with other calculations. Our calculations were done using model potentials STY Shabaev et al. (2013) and FG Flambaum and Ginges (2005).  Transition | This work | Sapirstein and Cheng (2005) | Fairhall et al. (2023) | Tran Tan and Derevianko (2023)
---|---|---|---|---
Method | STY | FG | Ab initio | FG | FG
$R_{\mathrm{PO}}(6s_{1/2}-6p_{1/2})$ | $0.262$ | $0.267$ | $0.326$ | $0.272$ | $0.26$
$R_{\mathrm{PO}}(6s_{1/2}-6p_{3/2})$ | $0.275$ | $0.280$ |  | $0.286$ | $0.28$
$R_{\mathrm{PO}}(6s_{1/2}-7p_{1/2})$ | $-2.89$ | $-2.97$ |  | $-2.80$ | $-2.89$
$R_{\mathrm{PO}}(6s_{1/2}-7p_{3/2})$ | $-1.84$ | $-1.84$ |  | $-1.76$ | $-1.84$
$R_{\mathrm{PO}}(7s_{1/2}-6p_{1/2})$ | $-0.442$ | $-0.452$ |  | $-0.433$ | $-0.44$
$R_{\mathrm{PO}}(7s_{1/2}-6p_{3/2})$ | $-0.369$ | $-0.376$ |  | $-0.359$ | $-0.37$

In Table 2 we compare our results for $R_{\mathrm{PO}}$ with results of other
groups. Sapirstein and Cheng (2005) made ab initio QED calculation using local
Kohn-Sham potential. They calculated both corrections and found, that the
vertex correction $R_{\mathrm{V}}$ was significantly smaller,
$R_{\mathrm{V}}=-0.065$. Other calculations used non-local DF operator and
model potential Flambaum and Ginges (2005); Ginges and Berengut (2016, 2016).
Here we made calculations with model potentials Shabaev et al. (2013) and
Flambaum and Ginges (2005), which we denote as STY and FG potentials
respectively.

Table 2 shows that the results obtained with two QED potentials agree within 2
– 3%. The difference between different groups is somewhat bigger, about 4 –
6%. There are small differences in the model potentials described in Refs.
Flambaum and Ginges (2005) and Ginges and Berengut (2016, 2016).
Incompleteness of the basis sets can be another reasons of the disagreements.
The only direct QED calculation is about 20% larger, than all three
calculations with model potentials. Note that this calculation was done with
the local Kohn-Sham potential instead of the DF potential used by all other
groups. Because of that, it is important to check, how much QED corrections
depend on the choice of the potential. Note also that parametrization (8)
allows to partly compensate for the differences in the value of the initial
amplitude $E1_{0}$ in different approximations.

###  III.2 Local screening potentials

We made calculations with the model QED operator Shabaev et al. (2013) and
three local screening potentials: the first one is the core-Hartee potential
(CH) induced by the core electrons of Cs+ ground-state configuration, the
other two are Kohn-Sham (KS) Kohn and Sham (1965) and Slater (S) Slater and
Johnson (1972) potentials corresponding to the ground-state configuration of
the neutral Cs. The proper asymptotic behavior of the potentials KS and S is
restored by including the Latter correction Latter (1955). In Table 3 the
results obtained using the local screening potentials are compared with the
Dirac-Fock results. We see that calculations of QED corrections
$\delta\varepsilon$ with local potentials do not correlate with the DF
calculations for all partial waves, but the $s$ wave. The reason for this
discrepancy is obvious. In the DF approximation the energy correction comes
from the direct contribution of the perturbation and the contribution from the
adjustment of the self-consistent field (relaxation of the core) $\delta
V_{\mathrm{SCF}}$. The former is large for the $s$ wave and rapidly decreases
with the orbital quantum number $l$. At the same time, the relaxation of the
core affects all valence orbitals in a similar way. Two contributions become
comparable already for the valence $p$ orbitals. For the local potentials
there is only direct contribution. The KS potential reproduces valence orbital
energies $\varepsilon$ better than two other potentials, while QED corrections
$\delta\varepsilon$ to the $s$ wave energies for the CH and KS potentials are
of the similar accuracy. Slater potential is clearly the least accurate.

Table 3: QED corrections $\delta\varepsilon$ to the binding energies $\varepsilon$ in Cs calculated in DF approximation and with three local potentials, core-Hartree (CH), Kohn-Sham (KS), and Slater (S). All calculations were done using model QED potential Shabaev et al. (2013) (STY). Orbital | Dirac-Fock | core-Hartree | Kohn-Sham | Slater
---|---|---|---|---
| $\varepsilon$ | $\delta\varepsilon$ | $\varepsilon$ | $\delta\varepsilon$ | $\varepsilon$ | $\delta\varepsilon$ | $\varepsilon$ | $\delta\varepsilon$
| (a.u.) | (cm-1) | (a.u.) | (cm-1) | (a.u.) | (cm-1) | (a.u.) | (cm-1)
$6s_{1/2}$ | $0.12737$ | $-15.80$ | $0.11998$ | $-16.13$ | $0.12394$ | $-16.22$ | $0.13573$ | $-23.54$
$6p_{1/2}$ | $0.08562$ | $0.80$ | $0.08130$ | $-0.20$ | $0.08489$ | $-0.17$ | $0.08908$ | $-0.20$
$6p_{3/2}$ | $0.08379$ | $0.16$ | $0.07876$ | $-1.09$ | $0.08278$ | $-0.91$ | $0.08676$ | $-0.96$
$5d_{3/2}$ | $0.06442$ | $2.15$ | $0.05898$ | $0.04$ | $0.06041$ | $0.07$ | $0.07735$ | $0.38$
$5d_{5/2}$ | $0.06453$ | $2.03$ | $0.05886$ | $-0.05$ | $0.06011$ | $-0.09$ | $0.07570$ | $-0.49$
$7s_{1/2}$ | $0.05519$ | $-4.30$ | $0.05327$ | $-4.71$ | $0.05457$ | $-4.85$ | $0.05791$ | $-6.40$
$7p_{1/2}$ | $0.04202$ | $0.29$ | $0.04060$ | $-0.07$ | $0.04184$ | $-0.06$ | $0.04334$ | $-0.07$
$7p_{3/2}$ | $0.04137$ | $0.05$ | $0.03971$ | $-0.39$ | $0.04108$ | $-0.32$ | $0.04251$ | $-0.34$
$6d_{3/2}$ | $0.03609$ | $1.22$ | $0.03318$ | $0.11$ | $0.03422$ | $0.19$ | $0.04152$ | $0.24$
$6d_{5/2}$ | $0.03609$ | $0.88$ | $0.03311$ | $-0.13$ | $0.03402$ | $-0.22$ | $0.04106$ | $-0.32$

Taking into account that core relaxation is very important for the $p$ and $d$
waves, the local approximation is clearly inapplicable for the QED corrections
to the $E1$ transitions between $p$ and $d$ orbitals. The results of the
calculations for the $s\leftrightarrow p$ transitions are presented in Table
4. Again, the KS potential gives better agreement with DF, than two other
potentials. In most cases the difference is less, or about 10%. The only
exception is the $6s_{1/2}-7p_{1/2}$ transition, where the difference is 24%.
For the CH and S potentials the average differences are about 14% and 22%
respectively.

Table 4: Comparison of the PO QED corrections $R_{\mathrm{PO}}$ to the $s\leftrightarrow p$ transition amplitudes in Cs for different potentials. For local potentials we also give the difference with DF results in percent. All calculations were done using model QED potential Shabaev et al. (2013) (STY). Transition | Dirac-Fock | core-Hartree | Kohn-Sham | Slater
---|---|---|---|---
$R_{\mathrm{PO}}(6s_{1/2}-6p_{1/2})$ | $0.266$ | $0.299$ | 12% | $0.281$ | 6% | $0.365$ | 37%
$R_{\mathrm{PO}}(6s_{1/2}-6p_{3/2})$ | $0.274$ | $0.316$ | 15% | $0.296$ | 8% | $0.389$ | 42%
$R_{\mathrm{PO}}(6s_{1/2}-7p_{1/2})$ | $-2.856$ | $-2.646$ | 7% | $-3.541$ | 24% | $-3.000$ | 5%
$R_{\mathrm{PO}}(6s_{1/2}-7p_{3/2})$ | $-1.810$ | $-1.388$ | 23% | $-1.945$ | 7% | $-1.814$ | 0%
$R_{\mathrm{PO}}(7s_{1/2}-6p_{1/2})$ | $-0.446$ | $-0.417$ | 7% | $-0.433$ | 3% | $-0.492$ | 10%
$R_{\mathrm{PO}}(7s_{1/2}-6p_{3/2})$ | $-0.373$ | $-0.311$ | 17% | $-0.341$ | 8% | $-0.403$ | 8%
$R_{\mathrm{PO}}(7s_{1/2}-7p_{1/2})$ | $0.255$ | $0.311$ | 22% | $0.284$ | 11% | $0.372$ | 46%
$R_{\mathrm{PO}}(7s_{1/2}-7p_{3/2})$ | $0.292$ | $0.318$ | 9% | $0.289$ | 1% | $0.380$ | 30%

We conclude that for Cs the KS potential provides higher accuracy for the QED
calculations than CH and S potentials. Note, that the result of Sapirstein and
Cheng (2005) for the $6s_{1/2}-6p_{1/2}$ transition cited in Table 2 was
obtained using KS potential. Our KS value is closer to their value, then the
DF one. Thus, the approximately 20% difference between ab initio calculation
Sapirstein and Cheng (2005) and calculations with model QED potentials
presented in Table 2 can be attributed in part to the local potential error.

##  IV Conclusions

Many numerical methods for atomic calculations require one-electron basis
sets. Here we suggested the differential ansatz to generate virtual orbitals
from the DF orbitals of the core and valence electrons. The number of such
orbitals cannot exceed the number of DF orbitals. More virtual orbitals can be
added using B-splines, or other traditional methods. Our test calculations for
Cs show good convergence of calculations of different properties, including
QED corrections and hyperfine structure constants.

We also tested three local potentials, which are often used for the QED
calculations in many-electron atoms and found out that for the lower partial
waves the Kohn-Sham local potential demonstrated best agreement with the DF
potential. For the higher partial waves the results obtained with all tested
local potentials are very different from the results obtained with the non-
local DF potential. These differences are caused by the absences of the
contribution from the change of the self-consistent field induced by the
perturbation. This effect becomes dominant for the partial waves with $l>1$
making calculation of the respective QED corrections problematic. In order to
improve the accuracy of the ab initio QED calculations in heavy atoms it is
necessary to include core-relaxation effects by going to the next order of the
perturbation theory in the electron-electron interaction. The next step in
quantum-mechanical calculations of the QED corrections requires to construct
effective operators for the vertex contribution.

###### Acknowledgements.

We are grateful to Ilya Tupitsyn and Vladimir Yerokhin for extremely useful
discussions and suggestions. This work was supported by the Russian Science
Foundation grant # 23-22-00079.

## References

  * Dzuba et al. (1996) V. A. Dzuba, V. V. Flambaum, and M. G. Kozlov, Phys. Rev. A 54, 3948 (1996).
  * Kozlov (2004) M. G. Kozlov, Int. J. Quant. Chem. 100, 336 (2004).
  * Dzuba (2005) V. A. Dzuba, Phys. Rev. A 71, 032512 (2005).
  * Safronova et al. (2009) M. S. Safronova, M. G. Kozlov, W. R. Johnson, and D. Jiang, Phys. Rev. A 80, 012516 (2009).
  * Kahl and Berengut (2019) E. Kahl and J. Berengut, Comput. Phys. Commun. 238, 232 (2019).
  * Pyykkö and Zhao (2003) P. Pyykkö and L.-B. Zhao, J. Phys. B: At. Mol. Opt. Phys. 36, 1469 (2003).
  * Flambaum and Ginges (2005) V. V. Flambaum and J. S. M. Ginges, Phys. Rev. A 72, 052115 (2005).
  * Tupitsyn and Berseneva (2013) I. I. Tupitsyn and E. V. Berseneva, Opt. Spectrosc. 114, 682 (2013).
  * Cheung et al. (2021) C. Cheung, M. Safronova, and S. Porsev, Symmetry 13, 621 (2021), ISSN 2073-8994.
  * Johnson et al. (1988) W. R. Johnson, S. A. Blundell, and J. Sapirstein, Phys. Rev. A 37, 307 (1988).
  * Tupitsyn and Loginov (2003) I. I. Tupitsyn and A. V. Loginov, Opt. Spectrosc. 94, 319 (2003).
  * Bogdanovich and Žukauskas (1983) P. Bogdanovich and G. Žukauskas, Sov. Phys. Collection 23, 18 (1983).
  * Kozlov et al. (1996) M. G. Kozlov, S. G. Porsev, and V. V. Flambaum, J. Phys. B 29, 689 (1996).
  * Kozlov and Tupitsyn (2019) M. Kozlov and I. Tupitsyn, Atoms 7, 92 (2019).
  * Shabaev et al. (2013) V. M. Shabaev, I. I. Tupitsyn, and V. A. Yerokhin, Phys. Rev. A 88, 012513 (2013).
  * Stanton and Havriliak (1984) R. E. Stanton and S. Havriliak, J. Chem. Phys. 81, 1910 (1984).
  * Shabaev et al. (2015) V. Shabaev, I. Tupitsyn, and V. Yerokhin, Comp. Phys. Commun. 189, 175 (2015).
  * Sapirstein et al. (2004) J. Sapirstein, K. Pachucki, and K. T. Cheng, Phys. Rev. A 69, 022113 (2004).
  * Sapirstein and Cheng (2005) J. Sapirstein and K. T. Cheng, Phys. Rev. A 71, 022503 (2005).
  * Fairhall et al. (2023) C. J. Fairhall, B. M. Roberts, and J. S. M. Ginges, Phys. Rev. A 107, 022813 (2023).
  * Tran Tan and Derevianko (2023) H. B. Tran Tan and A. Derevianko, Phys. Rev. A 107, 042809 (2023).
  * Roberts et al. (2023) B. M. Roberts, C. J. Fairhall, and J. S. M. Ginges, Phys. Rev. A 107, 052812 (2023).
  * Ginges and Berengut (2016) J. S. M. Ginges and J. C. Berengut, J. Phys. B: At. Mol. Opt. Phys. 49, 095001 (2016).
  * Ginges and Berengut (2016) J. S. M. Ginges and J. C. Berengut, Phys. Rev. A 93, 052509 (2016).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965), ISSN 0031-899X.
  * Slater and Johnson (1972) J. C. Slater and K. H. Johnson, Phys. Rev. B 5, 844 (1972), ISSN 0556-2805.
  * Latter (1955) R. Latter, Phys. Rev. 99, 510 (1955), ISSN 0031-899X.