# Optimal Probabilistic Constellation Shaping for Covert Communications

Shuai Ma, Yunqi Zhang, Haihong Sheng, Hang Li, Jia Shi, Long Yang, Youlong Wu,
Naofal Al-Dhahir, and Shiyin Li Shuai Ma,Yunqi Zhang, Haihong Sheng and Shiyin
Li are with the School of Information and Control Engineering, China
University of Mining and Technology, Xuzhou 221116, China (e-mail:
mashuai001@cumt.edu.cn).

###### Abstract

In this paper, we investigate the optimal probabilistic constellation shaping
design for covert communication systems from a practical view. Different from
conventional covert communications with equiprobable constellations
modulation, we propose non-equiprobable constellations modulation schemes to
further enhance the covert rate. Specifically, we derive covert rate
expressions for practical discrete constellation inputs for the first time.
Then, we study the covert rate maximization problem by jointly optimizing the
constellation distribution and power allocation. In particular, an approximate
gradient descent method is proposed for obtaining the optimal probabilistic
constellation shaping. To strike a balance between the computational
complexity and the transmission performance, we further develop a framework
that maximizes a lower bound on the achievable rate where the optimal
probabilistic constellation shaping problem can be solved efficiently using
the Frank-Wolfe method. Extensive numerical results show that the optimized
probabilistic constellation shaping strategies provide significant gains in
the achievable covert rate over the state-of-the-art schemes.

###### Index Terms:

Covert communications, probabilistic constellation shaping, achievable rate.

##  I Introduction

Radio frequency (RF) based wireless communication is inevitably susceptible to
eavesdropping due to the broadcast nature of the electromagnetic waves. With
the ever-growing Internet of Things (IoT) applications, the security issue has
become more and more crucial in the future sixth-generation (6G) wireless
networks [1]. For example, in large enterprise buildings, hospitals,
factories, communication may be very sensitive to a hostile adversary. The
conventional cryptography approaches [2] usually focus on protecting the
transmission content or increasing message decoding complexity, while the
physical-layer security [3, 4] approaches exploit the intrinsic wireless
fading channels properties to minimize the information leakage to the
eavesdroppers. In fact, a higher level of security is to hide the existence of
the communication, which not only can be applied in all aforementioned
application scenarios [5, 6, 7, 8], but also meets more critical demands from
military or security agencies. To address this high level security, covert
communications [9, 10], which shields the existence of message transmissions
against the detection of a warden, are emerging as a cutting-edge wireless
communication security technique, and have recently attracted significant
research attention. Note that, covert communication aims to hide the
communication behavior from the eavesdropper, while physical layer security
tries to reduce the interception information of the eavesdropper.

The basic idea of covert communications is as follows. The legitimate
transmitter (Alice) transmits messages to the paired receiver (Bob), while
guaranteeing a low detection probability for a Warden (Willie). Although such
idea has been realized by spread-spectrum techniques [11] for several decades,
the information-theoretic limits of covert communications, which is also
referred to as low probability of detection (LPD) communications in some
literature, were only recently derived [12, 10, 13, 14]. In particular, the
authors in [12] firstly demonstrated that in additive white Gaussian noise
(AWGN) channels, Alice can reliably send at most
$\mathcal{O}\left(\sqrt{n}\right)$ bits to Bob in $n$ channel usages under the
covert requirement. Such result is also called the square root law (SRL).
Subsequent works have extended this result to various channel models such as
binary symmetric channels [15], broadcast communications [16], multiple access
channels [17], and interference channels[18].

Although the SRL indicates that the asymptotic achievable rate under the
covert requirement approaches zero, many researchers have shown that the SRL
limit can be beaten by exploiting additional techniques in the considered
covert communication scenario. These methods include: taking advantage of the
ignorance of the transmission time at Willie [19]; applying an intelligent
reflecting surface [20, 21, 22]; exploring the molecular absorption or
scattering feature of the Terahertz spectrum [23, 24]; cooperative jamming
[25] or uninformed jamming [26, 27, 28, 29]; jointly optimizing the beam
training and data transmission for millimeter-wave communication [30];
exploiting the uncertainty noise power or channel state information (CSI) at
Willie [31, 32, 33, 34, 35]; robust beamforming design [36, 37], over a finite
number of channel uses[38]; applying a full-duplex transceiver [39, 40, 41,
42]; intermediate relay [43, 44] or exploring unmanned aerial vehicle (UAV) as
mobile relay [45, 46].

To be more specific, it is shown in [19] that Alice can covertly transmit
${\cal O}\left({\min\left\\{{n,\sqrt{n\log T\left(n\right)}}\right\\}}\right)$
bits to Bob. When Willie lacks the knowledge of his noise power, Alice can
reliably transmit $\mathcal{O}\left(n\right)$ bits [31, 32, 33]. With the aid
of an uninformed jammer [26], Alice can also achieve positive transmission
rate. With a finite number of channel uses, a uniformly distributed power
allocation scheme was proposed [38] to enhance the covert transmission. In
[39], the authors showed that, the effective throughput under delay
constraints can be improved by adding artificial noise (AN) at the full-duplex
receiver. For a one-way relay network, the authors in [43] studied the
performance limits of convert communications of an energy harvesting relay. In
[44], a multiple-relay network was considered, and both the maximum throughput
and the minimum end-to-end delay routing algorithms were developed with
multiple Willies. In [36], two probabilistic metrics, called the covert outage
probability and the connectivity probability, were analyzed for multi-antenna
covert communications with randomly located wardens and interferers. Using
Kullback-Leibler divergence, i.e., $D\left({{p_{1}}||{p_{0}}}\right)$ or
$D\left({{p_{0}}||{p_{1}}}\right)$ to measure the covertness, a Gaussian input
distribution was shown to be optimal for the covert metric
$D\left({{p_{1}}||{p_{0}}}\right)$, and not optimal for the covert metric
$D\left({{p_{0}}||{p_{1}}}\right)$[9], where ${p_{0}}$ and ${p_{1}}$ represent
Willie’s received signal distributions when covert communications occur and
not occur, respectively.

The aforementioned research advances in covert communications mainly make the
assumption of a Gaussian input distribution at the transmitter side, which can
hardly be realized in practical communication systems. In fact, the
information symbols in practical communication systems are realized in the
form of discrete constellation points, i.e., finite alphabet inputs, such as
pulse amplitude modulation (PAM) and multiple-quadrature amplitude modulation
(M-QAM). In [47], both the lower bound of covert transmission probability and
throughput maximization have been analyzed with discrete constellation inputs.
The discrete constellation points are assumed to be equally likely, which is
not optimal for practical covert communications, especially for high-order
modulation schemes. So far, the optimal discrete constellation inputs of
covert communication are still not well discussed in the literature.

Motivated by the above background, we develop information-theoretic limits of
covert communications with probabilistic constellation shaping. First, we
derive the achievable rate expression of the system with the discrete
constellation input signals, rather than the Gaussian inputs adopted in most
of the existing works. Then, we investigate the performance with the optimal
input distribution. Our results provide a practical design framework for
covert communication systems. The main contributions of this paper are
summarized as follows:

  * •

Generally, the inputs of practical communications systems follow a finite-set
discrete distribution rather than a Gaussian distribution. To evaluate
performance, we derive the achievable rate expressions for an arbitrary
discrete distributed input. Comparing to the existing rate expressions with
equiprobable discrete constellation points, the derived expressions are more
general and practical. Since the derived rate expression is not in closed-
form, we further derive both lower and upper bounds. All these results can be
used as performance metrics for the considered covert communication system.

  * •

Furthermore, we design optimal discrete constellation inputs to maximize the
exact covert rate under the covertness constraints, transmit power
limitations, and the signal distribution requirements, which is a challenging
problem since neither the exact covert rate nor the covertness constraint has
an analytical expression. To efficiently solve it, we conservatively transform
the covertness constraint into its upper bound with closed-form expression.
Then, we adopt the numerical integration method to approximate the covert rate
objective function and its gradient. Afterwards, the optimal probability
distributions of the discrete constellation are calculated by the approximate
gradient descent method, where the step sizes are calculated by the
backtracking line search.

  * •

To reduce the computation complexity of the design problem, we further adopt
the derived lower bound as the covert rate performance metric. To overcome the
non-convexity challenge, this problem is iteratively solved by the proposed
Frank-Wolfe method.

The rest of this paper is organized as follows. The system model and the
derivation of Bob’s achievable rate are presented in Section II. The optimal
probabilistic constellation shaping design for covert communications is
provided in Section III. The probabilistic constellation shaping design and
its approximations are presented in Section IV. In Section V, we evaluate the
proposed probabilistic constellation shaping design using numerical results.
Finally, we conclude the paper in Section VI.

_Notations_ : The vectors and matrices are represented by boldfaced lowercase
and uppercase letters, respectively. The notations
${\left(\cdot\right)^{\rm{*}}}$, ${{\mathbb{E}}}\left\\{\cdot\right\\}$,
$\left\|\cdot\right\|$, ${\rm{Tr}}\left(\cdot\right)$, ${\mathop{\rm
Re}\nolimits}\left(\cdot\right)$ and ${\mathop{\rm
Im}\nolimits}\left(\cdot\right)$ represent the conjugate, the expectation,
Frobenius norm, trace, the real part and imaginary part of its argument,
respectively. And $\odot$ is Hadamard Product,i.e.${A_{m\times
n}}\left[{{a_{ij}}}\right]\odot{B_{m\times
n}}\left[{{b_{ij}}}\right]={C_{m\times n}}\left[{{a_{ij}}{b_{ij}}}\right]$.
The operator ${\bf{A}}\underline{\succ}{\bf{0}}$ means ${\bf{A}}$ is positive
semidefinite. The notation $\mathcal{CN}\left({\mu,{\sigma^{2}}}\right)$
denotes a complex-valued circularly symmetric Gaussian distribution with mean
$\mu$ and variance ${\sigma^{2}}$.

##  II System model

Figure 1:  The system model of covert communication.

Consider a typical covert communication scenario as illustrated in Fig. 1,
where Alice and Bob are a legitimate communication pair, and Willie is the
eavesdropper. Each one of them is equipped with a single antenna. Let
${{{g}}_{\rm{b}}}\sim{\cal CN}\left({{{0}},\sigma_{1}^{2}}\right)$ and
${{{g}}_{\rm{w}}}\sim{\cal CN}\left({{{0}},\sigma_{2}^{2}}\right)$ denote the
Rayleigh flat fading channel from Alice to Bob and Willie, respectively [48],
where $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$ are the variances of
${{{g}}_{\rm{b}}}$ and ${{{g}}_{\rm{w}}}$. Let $x\left[i\right]$ denote
Alice’s transmitted symbol at the $i$-th channel use, where $i=1,...,N$, and
$N$ is the total number of channel uses.

###  II-A Achievable Rate of Bob

As it is the case in practice, $x\left[i\right]\in\Omega$ follows a discrete
constellation distribution instead of Gaussian distribution. Here, $\Omega$
denotes a discrete constellation set with $K$ discrete points
${\left\\{{{x_{k}}}\right\\}_{1\leq k\leq K}}$, i.e.,

$\displaystyle\Omega=\left\\{{X\left|{\begin{array}[]{*{20}{l}}\begin{array}[]{l}\Pr\left({X={x_{k}}}\right)={p_{k}}\geq{\rm{0,}}\sum\limits_{k=1}^{K}{{p_{k}}}=1,\\\
\sum\limits_{k=1}^{K}{{p_{k}}}{\left|{{x_{k}}}\right|^{2}}\leq{P_{\rm{A}}},{x_{k}}\in{\mathbb{C}},k=1,...,K\end{array}\end{array}}\right.}\right\\},$
(4)

where ${{x_{k}}}$ denotes the $k$th discrete point, ${{p_{k}}}$ denotes the
corresponding probability, and ${P_{\text{A}}}$ denotes the average power.

For the $i$-th channel use, the received signal at Bob
${y_{\rm{b}}}\left[i\right]$ is given as

$\displaystyle{y_{\rm{b}}}\left[i\right]={g_{\rm{b}}^{\rm{*}}}{x}\left[i\right]+{z_{\rm{b}}}\left[i\right],$
(5)

where
${z_{\rm{b}}}\left[i\right]\sim\mathcal{CN}\left({0,\sigma_{\rm{b}}^{2}}\right)$
denotes the received noise at Bob. Since $x\left[i\right]\in\Omega$, then the
likehood functions of ${y_{\rm{b}}}\left[i\right]$ is given as

$\displaystyle{p}\left({{y_{\rm{b}}}}\right)=\frac{1}{{{\pi}{{\sigma_{\rm{b}}^{2}}}}}\sum\limits_{k=1}^{K}{{p_{k}}\exp\left({-\frac{{{{\left|{{y_{\rm{b}}}-{g_{\rm{b}}^{\rm{*}}}{x_{k}}}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)}.$
(6)

Therefore, given the discrete constellation, the achievable rate of Bob
${R_{\rm{b}}}$ is given by

$\displaystyle{R_{\rm{b}}}=I\left({{y_{\rm{b}}}\left[i\right];{x}\left[i\right]}\right)$
(7a)
$\displaystyle=h\left({{y_{\rm{b}}}\left[i\right]}\right)-h\left({{z_{\rm{b}}}\left[i\right]}\right)$
(7b)
$\displaystyle=-\int_{-\infty}^{\infty}{{p}\left({{y_{\rm{b}}}}\right)}{\log_{2}}{p}\left({{y_{\rm{b}}}}\right)d{y_{\rm{b}}}-{\log_{2}}\pi
e\sigma_{b}^{2}$ (7c)
$\displaystyle=-\sum\limits_{k=1}^{K}{{p_{k}}}{\mathbb{E}_{{z_{\rm{b}}}}}\left\\{{{{\log}_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)+{z_{\rm{b}}}}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)}\right\\}$
$\displaystyle\quad-\frac{1}{{\ln 2}},$ (7d)

where $h\left(X\right)=-\int{f\left(x\right)}\log f\left(x\right)dx$ denotes
differential entropy, and ${f\left(x\right)}$ represents the probability
density function (PDF).

Based on the achievable rate expression in (7), we will further investigate
the optimal probability of discrete constellations for covert communications.

###  II-B Hypothesis Testing

According to the received signals, Willie attempts to decide whether Alice
covertly transmits information to Bob or not by performing an optimal
statistical hypothesis test (such as the Neyman-Pearson test). Specifically,
Willies needs to distinguish between two hypotheses: 1) the null hypothesis
${{\cal{H}}_{0}}$ indicating no transmission; 2) the hypothesis
${{\cal{H}}_{1}}$ indicating the transmission. Let $y_{\rm{w}}\left[i\right]$
denote the received signal at Willie in the $i$-th channel use. Under the two
hypotheses, the signal received at Willie is given as

$\displaystyle\quad{{\cal{H}}_{0}}:y_{\rm{w}}\left[i\right]={z_{\rm{w}}}\left[i\right],$
(8a)
$\displaystyle\quad{{\cal{H}}_{1}}:y_{\rm{w}}\left[i\right]=g_{\rm{w}}^{\rm{*}}x\left[i\right]+{z_{\rm{w}}}\left[i\right],$
(8b)

where
${z_{\rm{w}}}\left[i\right]\sim\mathcal{CN}\left({0,\sigma_{\rm{w}}^{2}}\right)$
denotes the received noise at Willie. Let ${{\cal{D}}_{1}}$ and
${{\cal{D}}_{0}}$, respectively, denote the binary decisions of Willie. Thus,
the total detection error probability of Willie is defined as [9, 49, 50]

$\displaystyle\xi=\Pr\left({{{\cal{D}}_{1}}|{{\cal{H}}_{0}}}\right)+\Pr\left({{{\cal{D}}_{0}}|{{\cal{H}}_{1}}}\right).$
(9)

Note that, $\Pr\left({{{\cal{D}}_{1}}|{{\cal{H}}_{0}}}\right)$ denotes the
false alarm probability that Willie believes ${{\cal{H}}_{1}}$ when Alice does
not transmit, and $\Pr\left({{{\cal{D}}_{0}}|{{\cal{H}}_{1}}}\right)$ denotes
the missed detection probability that Willie decides ${{\cal{H}}_{0}}$ when
Alice transmits. Moreover, let
${p_{y,0}}=f\left({{y_{\rm{w}}}\left|{{{\cal{H}}_{0}}}\right.}\right)$ and
${p_{y,1}}=f\left({{y_{\rm{w}}}\left|{{{\cal{H}}_{1}}}\right.}\right)$ denote
the likehood functions of $y_{\rm{w}}$ under ${{\cal{H}}_{0}}$ and
${{\cal{H}}_{1}}$, respectively. According to (8), we have

$\displaystyle{p_{y,0}}=\frac{1}{{{\pi}{{\sigma_{\rm{w}}^{2}}}}}\exp\left({-\frac{\left|{{y_{\rm{w}}}}\right|^{2}}{{\sigma_{\rm{w}}^{2}}}}\right),$
(10a)
$\displaystyle{p_{y,1}}=\frac{1}{{{\pi}{{\sigma_{\rm{w}}^{2}}}}}\sum\limits_{k=1}^{K}{{p_{k}}\exp\left({-\frac{{{{\left|{y_{\rm{w}}-{g_{\rm{w}}^{\rm{*}}}{x_{k}}}\right|}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}.$
(10b)

Let
${V_{T}}\left({{{p_{y,0}}\left\|{{p_{y,1}}}\right.}}\right)=\frac{1}{2}{\left\|{{p_{y,0}}-{p_{y,1}}}\right\|_{1}}$
denote the total variation distance between ${p_{y,0}}$ and ${p_{y,1}}$.
According to Theorem 13.1.1 in [50], the optimal detection error probability
of Willie is given as

$\displaystyle{\xi^{\rm{opt}}}=1-{V_{T}}\left({{{p_{y,0}}\left\|{{p_{y,1}}}\right.}}\right)=1-\frac{1}{2}{\left\|{{p_{y,0}}-{p_{y,1}}}\right\|_{1}}.$
(11)

However, in general,
${V_{T}}\left({{{p_{y,0}}\left\|{{p_{y,1}}}\right.}}\right)$ is difficult to
analyze. To address this issue, we apply Pinsker’s inequality [49] to obtain
an upper bound

$\displaystyle{V_{T}}\left({{{p_{y,0}}\left\|{{p_{y,1}}}\right.}}\right)\leq\sqrt{\frac{1}{2}D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)},$
(12)
$\displaystyle{V_{T}}\left({{{p_{y,0}}\left\|{{p_{y,1}}}\right.}}\right)\leq\sqrt{\frac{1}{2}D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)},$
(13)

where
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)=\int_{y}{{p_{y,0}}{{\log}_{2}}\frac{{{p_{y,0}}}}{{{p_{y,1}}}}}d{y_{\rm{w}}}$
denotes the Kullback-Leibler (KL) divergence from ${p_{y,0}}$ to ${p_{y,1}}$,
and
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)=\int_{y}{{p_{y,1}}{\log_{2}}\frac{{p_{y,1}}}{{p_{y,0}}}}d{y_{\rm{w}}}$
denotes the KL divergence from ${p_{y,1}}$ to ${p_{y,0}}$.

Based on the likehood functions of $y_{\rm{w}}$ in (10),
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ are, respectively, given
as

$\displaystyle
D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)=-\frac{1}{{\ln
2}}-{\mathbb{E}_{{z_{\rm{w}}}}}\left\\{{{{\log}_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}}\right.$
$\displaystyle\quad\times\left.{\exp\left({-\frac{{{{\left|{{z_{\rm{w}}}-g_{\rm{w}}^{\rm{*}}{x_{k}}}\right|}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right\\},$
(14a) $\displaystyle
D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)=\sum\limits_{k=1}^{K}{{p_{k}}}{\mathbb{E}_{{z_{\rm{w}}}}}\left\\{{{{\log}_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}}\right.$
$\displaystyle\quad\left.{\times\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}\left({{x_{k}}-{x_{j}}}\right)+{z_{\rm{w}}}}\right|}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right\\}+\frac{{\left|{g_{\rm{w}}^{2}{P_{\rm{A}}}}\right|}}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}+\frac{1}{{\ln 2}}.$ (14b)

Covert communication is achieved for a given $\varepsilon$ if the detection
error probability $\xi$ is no less than $1-\varepsilon$, i.e.,

$\displaystyle\xi\geq 1-\varepsilon,\varepsilon\in\left[{0,1}\right],$ (15)

where $\varepsilon$ is a small number determining the required covertness
level.

Therefore, to achieve covert communication with the given $\varepsilon$, i.e.,
$\xi\geq 1-\varepsilon$, the KL divergences of the likelihood functions should
satisfy one of the following constraints:

$\displaystyle D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}},$ (16a) $\displaystyle
D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq 2{\varepsilon^{2}}.$
(16b)

Both of the above constraints can meet the requirements of covert
communication, but these two constraints are not exactly the same, as we
discuss next.

##  III Optimal Signaling Design under Covert Constraints

In this section, we investigate the design of optimal probability of discrete
constellation points for covert transmission with covertness constraints
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq 2{\varepsilon^{2}}$ or
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq 2{\varepsilon^{2}}$
[12, 13, 10, 9].

###  III-A Case of ${D}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}}$

For the probabilistic constellation shaping scheme, we aim to maximize the
achievable rate of Bob ${R_{\rm{b}}}$ by optimizing the distribution of
discrete constellation inputs, while satisfying both the covert transmission
constraint and the discrete distribution constraint. Mathematically, the
covert discrete constellation input optimization problem can be formulated as
follows

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}~{}$
$\displaystyle{R_{\rm{b}}}$ (17a) $\displaystyle{\rm{s.t.}}~{}$ $\displaystyle
D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq 2{\varepsilon^{2}},$
(17b) $\displaystyle\Pr\left({X={x_{k}}}\right)={p_{k}}\geq{\rm{0}},$ (17c)
$\displaystyle\sum\limits_{k=1}^{K}{{p_{k}}{{\left|{{x_{k}}}\right|}^{2}}}\leq{P_{\text{A}}},$
(17d) $\displaystyle\sum\limits_{k=1}^{K}{{p_{k}}}=1,k=1,...,K.$ (17e)

Since the KL-divergence $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$
in (14b) is not an analytical expression, constraint (17b) is intractable. To
circumvent this, we derive an explicit upper bound for the KL-divergence
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$, which is given by

$\displaystyle{D_{{\rm{U}}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)=-{\log_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}{x_{k}}}\right|}^{\rm{2}}}}}{{\sigma_{\rm{w}}^{2}}}}\right).$
(18)

The details of derivations for (18) can be found in Appendix A. Based on (18),
problem (17) can be rewritten as

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}~{}$
$\displaystyle{R_{\rm{b}}}$ (19a) $\displaystyle{\rm{s.t.}}~{}$
$\displaystyle-{\log_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}{x_{k}}}\right|}^{\rm{2}}}}}{{\sigma_{\rm{w}}^{2}}}}\right)\leq
2{\varepsilon^{2}},$ (19b)
$\displaystyle\eqref{orig1b},\eqref{orig1c},\eqref{orig1d}.$

In order to put problem (19) in a more concise form, we define the following
variables

$\displaystyle{\bf{x}}\buildrel\Delta\over{=}{[{x_{1}},...,{x_{K}}]^{T}},$
(20a)
$\displaystyle{\bf{p}}\buildrel\Delta\over{=}{[{p_{1}},...,{p_{K}}]^{T}},$
(20b)
$\displaystyle{{{\bf{q}}}}\buildrel\Delta\over{=}{\left[{{{\log}_{2}}{{\bf{p}}^{T}}{{{\bf{\hat{q}}}}_{1}},...,{{\log}_{2}}{{\bf{p}}^{T}}{{{\bf{\hat{q}}}}_{K}}}\right]^{T}},$
(20c)
$\displaystyle{{\hat{\bf{q}}}_{l}}\buildrel\Delta\over{=}\left[\begin{array}[]{l}\exp\left({-\frac{{{{\left|{g_{\rm{b}}^{\rm{*}}\left({{x_{l}}-{x_{1}}}\right)+{z_{\rm{b}}}}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}...\\\
\exp\left({-\frac{{{{\left|{g_{\rm{b}}^{\rm{*}}\left({{x_{l}}-{x_{K}}}\right)+{z_{\rm{b}}}}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)\end{array}\right],\forall
l\in K,$ (20g)
$\displaystyle{\bf{t}}\buildrel\Delta\over{=}{\left[{\exp\left({-\frac{{{{{\left|{g_{\rm{w}}^{\rm{*}}{x_{1}}}\right|}}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right),...,\exp\left({-\frac{{{{{\left|{g_{\rm{w}}^{\rm{*}}{x_{K}}}\right|}}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right]^{T}},$
(20h)
$\displaystyle\phi\left({\bf{p}}\right)\buildrel\Delta\over{=}{\mathbb{E}_{{z_{\rm{b}}}}}\left\\{{{{\bf{p}}^{T}}{\bf{q}}}\right\\}.$
(20i)

Furthermore, the rate of Bob ${R_{\rm{b}}}$ and the upper bound of the KL-
divergence ${D_{\rm{U}}}\left({{p_{0}}||{p_{1}}}\right)$ can be, respectively,
rewritten as follows

$\displaystyle{R_{\rm{b}}}=-\phi\left({\bf{p}}\right)-\frac{1}{{\ln 2}},$
(21a)
$\displaystyle{D_{\rm{U}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)=-{\log_{2}}{{\bf{p}}^{T}}{\bf{t}}.$
(21b)

Since $\frac{1}{{\ln 2}}$ is constant and maximizing
$-\phi\left({\bf{p}}\right)$ is equivalent to minimizing
$\phi\left({\bf{p}}\right)$, problem (19) can be reformulated as

$\displaystyle\mathop{\min}\limits_{{\bf{p}}}$
$\displaystyle~{}\phi\left({\bf{p}}\right)$ (22a)
$\displaystyle{\rm{s}}.{\rm{t}}.$
$\displaystyle-{\log_{2}}{{\bf{p}}^{T}}{\bf{t}}\leq 2{\varepsilon^{2}},$ (22b)
$\displaystyle{{\bf{p}}^{T}}{{\bf{1}}_{K}}=1,$ (22c)
$\displaystyle{\bf{p}}^{T}\left({{\bf{x}}\odot{\bf{x}}}\right)\leq{P_{\rm{A}}},$
(22d) $\displaystyle{\bf{p}}\geq{\bf{0}},$ (22e)

where ${{{\bf{1}}_{K}}}$ is a $K\times 1$ vector with all elements equal to
$1$.

Problem (22) is now a convex problem, and we adopt the gradient projection
method to solve it. Specifically, let $\nabla\phi\left({\bf{p}}\right)$ denote
the gradient of the objective function (22a), which is given by

$\displaystyle\nabla\phi\left({\bf{p}}\right)$
$\displaystyle={\mathbb{E}_{{z_{\rm{b}}}}}\left\\{{{\bf{q}}+{\bf{Qp}}}\right\\}=\int_{-\infty}^{\infty}{{f_{{z_{\rm{b}}}}}\left({{z_{\rm{b}}}}\right)\left({{\bf{q}}+{\bf{Qp}}}\right)d}{z_{\rm{b}}}.$
(23a)

Here,
${f_{{z_{\rm{b}}}}}\left({{z_{\rm{b}}}}\right){\rm{=}}\frac{{\rm{1}}}{{\pi\sigma_{\rm{b}}^{2}}}\exp\left({-\frac{{{{\left|{{z_{\rm{b}}}}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)$
denotes the probability density function of ${z_{\rm{b}}}$,
${\bf{Q}}\buildrel\Delta\over{=}\left[{{Q_{i,j}}}\right]$, where
${Q_{i,j}}\buildrel\Delta\over{=}\frac{{{\bf{\hat{q}}}_{j}^{T}{{\bf{e}}_{i}}}}{{{\bf{\hat{q}}}_{j}^{T}{\bf{p}}\ln
2}}$, and ${{\bf{e}}_{i}}$ is the unit vector where the $i$th element is $1$
and the other elements are $0$.

However, neither the objective function $\phi\left({\bf{p}}\right)$ or the
gradient $\nabla\phi\left({\bf{p}}\right)$ has an analytic expression. To
tackle this challenge, we adopt the numerical integration method to
approximate $\phi\left({\bf{p}}\right)$ and $\nabla\phi\left({\bf{p}}\right)$,
i.e.,

$\displaystyle\quad\tilde{\phi}\left({\bf{p}}\right)$
$\displaystyle=\int_{-{\tau_{1}}}^{{\tau_{1}}}{{f_{{z_{\rm{b}}}}}\left({{z_{\rm{b}}}}\right)\left({{{\bf{p}}^{T}}{\bf{q}}}\right)d}{z_{\rm{b}}},$
(24a) $\displaystyle\quad\nabla\tilde{\phi}\left({\bf{p}}\right)$
$\displaystyle=\int_{-{\tau_{2}}}^{{\tau_{2}}}{{f_{{z_{\rm{b}}}}}\left({{z_{\rm{b}}}}\right)\left({{\bf{q}}+{\bf{Qp}}}\right)d}{z_{\rm{b}}},$
(24b)

where $\left[{-{\tau_{1}},{\tau_{1}}}\right]$ and
$\left[{-{\tau_{2}},{\tau_{2}}}\right]$ denote the integration intervals, by
defining ${\tau_{1}}>0$, ${\tau_{2}}>0$. Furthermore,
$\tilde{\phi}\left({\bf{p}}\right)$ and
$\nabla\tilde{\phi}\left({\bf{p}}\right)$ denote the approximations of the
objective function and its gradient, respectively.

Furthermore, let ${{\bf{p}}_{0}}$ denote a feasible starting point, and
${{\bf{p}}_{n}}$ denote the $n$th iteration feasible point. With approximate
gradient $\nabla\tilde{\phi}\left({\bf{p}}\right)$, the gradient descent
iteration step is given as

$\displaystyle{{{\bf{\hat{p}}}}_{n+1}}={{\bf{p}}_{n}}-{\alpha_{n}}\nabla\tilde{\phi}\left({{\bf{p}}_{n}}\right),$
(25)

where ${\alpha_{n}}\in\left({0,1}\right]$ is the stepsize of the $n$th
iteration. To choose a proper step size ${\alpha_{n}}$ with a sufficient
decrease, we adopt the backtracking line search algorithm, given in Algorithm
$1$.

Then, we project ${{{{\bf{\hat{p}}}}_{n+1}}}$ in the feasible region of
problem (22). Specifically, when ${{{{\bf{\hat{p}}}}_{n+1}}}$ satisfies
constraints (22b)-(22e), ${{{{\bf{p}}}_{n+1}}}={{{{\bf{\hat{p}}}}_{n+1}}}$ is
in the feasible region. Otherwise, we need to find the closest point
${{{{\bf{p}}}_{n+1}}}$ in the feasible region as the projection of
${{{{\bf{\hat{p}}}}_{n+1}}}$.

Mathematically, the projection operation of ${{{{\bf{\hat{p}}}}_{n+1}}}$ can
be formulated as follows

$\displaystyle\mathop{\min}\limits_{{{\bf{p}}_{n+1}}}$
$\displaystyle~{}{\left\|{{{\bf{p}}_{n+1}}-{{{\bf{\hat{p}}}}_{n+1}}}\right\|^{2}}$
(26a) $\displaystyle{\rm{s}}.{\rm{t}}.$
$\displaystyle~{}-{\log_{2}}{{{\bf{p}}_{n+1}^{T}}}{\bf{t}}\leq
2{\varepsilon^{2}},$ (26b)
$\displaystyle~{}{{{\bf{p}}_{n+1}^{T}}}{{\bf{1}}_{K}}=1,$ (26c)
$\displaystyle~{}{{\bf{p}}_{n+1}^{T}}\left({{\bf{x}}\odot{\bf{x}}}\right)\leq{P_{\rm{A}}},$
(26d) $\displaystyle~{}{{\bf{p}}_{n+1}}\geq{\bf{0}}.$ (26e)

Algorithm 1 : Backtracking Line Search for Stepsize ${\alpha_{n}}$.

1:Input: ${{{\bf{p}}}_{n}}$, $\tilde{\phi}{\left({{{\bf{p}}_{n}}}\right)}$,
$\nabla\tilde{\phi}\left({{\bf{p}}_{n}}\right)$ and $\bar{\alpha}>0$,
$\rho,c\in\left({0,1}\right)$;

2:Update
${{{\bf{\hat{p}}}}_{n+1}}={{\bf{p}}_{n}}-{\alpha_{n}}\nabla\tilde{\phi}\left({{\bf{p}}_{n}}\right)$;

3:if ${{{\bf{\hat{p}}}}_{n+1}}$ satisfies constraints (22b)-(22e)

4: ${{{{\bf{p}}}_{n+1}}}={{{{\bf{\hat{p}}}}_{n+1}}}$;

5:else

6: Solve problem (26) over ${{{\bf{\hat{p}}}}_{n+1}}$ to obtain
${{{{\bf{p}}}_{n+1}}}$;

7:end;

8:While
$\tilde{\phi}\left({{{{\bf{p}}}_{n+1}}}\right)\leq\tilde{\phi}\left({{{\bf{p}}_{n}}}\right)+c\bar{\alpha}\nabla\tilde{\phi}{\left({{{\bf{p}}_{n}}}\right)^{T}}\left({{{{{\bf{p}}}_{n+1}}}-{{\bf{p}}_{n}}}\right)$;

9: $\bar{\alpha}\leftarrow\rho\bar{\alpha}$;

10: Update
${{{\bf{\hat{p}}}}_{n+1}}={{\bf{p}}_{n}}-{\alpha_{n}}\nabla\tilde{\phi}\left({{\bf{p}}_{n}}\right)$;

11: if ${{{\bf{\hat{p}}}}_{n+1}}$ satisfies constraints (22b)-(22e)

12: ${{{{\bf{p}}}_{n+1}}}={{{{\bf{\hat{p}}}}_{n+1}}}$;

13: else

14: solve problem (26) over ${{{\bf{\hat{p}}}}_{n+1}}$ to obtain
${{{{\bf{p}}}_{n+1}}}$;

15: end;

16:end;

17:return ${\alpha_{n}}=\bar{\alpha}$.

Therefore, we propose the approximate gradient descent projection method to
efficiently solve problem (22), which is summarized in Algorithm 2.

Algorithm 2 Inexact Gradient Descent Projection Method.

1:Input:  choose $K\geq 2$ and choose a random starting point ${{\bf{p}}_{0}}$
which satisfies constraints (22b)-(22e), set ${c_{2}}$ as the stopping
parameter and $n=0$;

2:Repeat

3: Let $n\leftarrow n+1$;

4: Update
$\tilde{\phi}\left({{{\bf{p}}_{n-1}}}\right)=\int_{-{\tau_{1}}}^{{\tau_{1}}}{{f_{{z_{\rm{b}}}}}\left({{z_{\rm{b}}}}\right)\left({{{{\bf{p}}_{n-1}^{T}}}{\bf{q}}}\right)d}{z_{\rm{b}}}$;

5: Update
${\nabla\tilde{\phi}\left({{{\bf{p}}_{n-1}}}\right)}=\int_{-{\tau_{2}}}^{{\tau_{2}}}{{f_{{z_{\rm{b}}}}}\left({{z_{\rm{b}}}}\right)\left({{\bf{q}}+{\bf{Q}}{{\bf{p}}_{n-1}}}\right)d}{z_{\rm{b}}}$;

6: Compute stepsize ${\alpha_{n-1}}$ by Algorithm $1$;

7: Update
${{{\bf{\hat{p}}}}_{n}}={{\bf{p}}_{n-1}}-{\alpha_{n-1}}\nabla\tilde{\phi}\left({{\bf{p}}_{n-1}}\right)$;

8: if ${{{\bf{\hat{p}}}}_{n}}$ satisfies constraints (22b)-(22e)

9: ${{{{\bf{p}}}_{n}}}={{{{\bf{\hat{p}}}}_{n}}}$;

10: else

11: solve problem (26) over ${{{\bf{\hat{p}}}}_{n}}$ to obtain
${{{{\bf{p}}}_{n}}}$;

12: end;

13:Until $\left\|{{{\bf{p}}_{n}}-{{\bf{p}}_{n-1}}}\right\|\leq{c_{2}}$;

14:Output ${{\bf{P}}^{{\rm{opt}}}}={{\bf{p}}_{n}}$.

###  III-B Case of ${D}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}}$

In this subsection, we further consider the other covert constraint
${D}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq 2{\varepsilon^{2}}$,
and the corresponding covert rate optimization problem can be formulated as

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}~{}$
$\displaystyle{R_{\rm{b}}}$ (27a) $\displaystyle{\rm{s.t.}}~{}$
$\displaystyle{D}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}},$ (27b)
$\displaystyle\Pr\left({X={x_{k}}}\right)={p_{k}}\geq{\rm{0}},$ (27c)
$\displaystyle\sum\limits_{k=1}^{K}{{p_{k}}{{\left|{{x_{k}}}\right|}^{2}}}\leq{P_{\rm{A}}},$
(27d) $\displaystyle\sum\limits_{k=1}^{K}{{p_{k}}}=1,k=1,...,K.$ (27e)

To handle intractable constraint (27b), we first derive an upper bound on
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ denoted by
${D_{{\rm{U}}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$, which is
given by

$\displaystyle{D_{{\rm{U}}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)=$
$\displaystyle\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)\right|}^{\rm{2}}}}}{{2\sigma_{\rm{w}}^{2}}}}\right)$
$\displaystyle+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}-1.$ (28)

The details of derivations for (III-B) are given in Appendix B. Then, problem
(27) can be reformulated as

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}~{}$
$\displaystyle{R_{\rm{b}}}$ (29a) $\displaystyle{\rm{s.t.}}~{}$
$\displaystyle{D_{{\rm{U}}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}},$ (29b)
$\displaystyle\eqref{orig2b},\eqref{orig2c},\eqref{orig2d}.$

By defining the following equations to simplify
${D_{{\rm{U}}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$

$\displaystyle{{\bf{s}}_{{\rm{w}},k}}\buildrel\Delta\over{=}\left[\begin{array}[]{l}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}\left({{x_{k}}-{x_{1}}}\right)}\right|}^{2}}}}{{2\sigma_{\rm{w}}^{2}}}}\right),...,\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}\left({{x_{k}}-{x_{K}}}\right)}\right|}^{2}}}}{{2\sigma_{\rm{w}}^{2}}}}\right)\end{array}\right]^{T},$
(30b)
$\displaystyle{{\bf{v}}_{\rm{w}}}\left({\bf{p}}\right)\buildrel\Delta\over{=}{\left[{{{\log}_{2}}{{\bf{p}}^{T}}{{\bf{s}}_{{\rm{w}},1}},...,{{\log}_{2}}{{\bf{p}}^{T}}{{\bf{s}}_{{\rm{w}},K}}}\right]^{T}},$
(30c)

we obtain

$\displaystyle{D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)={{\bf{p}}^{T}}{{\bf{v}}_{\rm{w}}}\left({\bf{p}}\right)+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{x}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}-1.$ (31)

Unfortunately, ${D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$
is non-convex in ${\bf{p}}$, and the covert constraint is also non-convex. To
handle this issue, we apply the first order Taylor expansion to
${D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$. Specifically,
the derivative of
${D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ at
${{{{\bf{\bar{p}}}}_{n}}}$ is given by

$\displaystyle\nabla{D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\left|{{}_{{\bf{p}}={{{{\bf{\bar{p}}}}_{n}}}}}\right.={{\bf{v}}_{\rm{w}}}\left({{{\bf{\bar{p}}}}_{n}}\right)+\nabla{{\bf{v}}_{\rm{w}}}{{{\bf{\bar{p}}}}_{n}},$
(32)

where
$\nabla{\bf{v}}_{\rm{w}}={\left[{\frac{{{{\bf{s}}_{{\rm{w}},1}}}}{{{{{\bf{\bar{p}}}}_{n}^{T}}{{\bf{s}}_{{\rm{w}},1}}}},...,\frac{{{{\bf{s}}_{{\rm{w}},K}}}}{{{{{\bf{\bar{p}}}}_{n}^{T}}{{\bf{s}}_{{\rm{w}},K}}}}}\right]_{K\times
K}}$. Then, the first order Taylor expansion of
${{\bf{p}}^{T}}{{\bf{v}}_{\rm{w}}}\left({\bf{p}}\right)$ is given as follows

$\displaystyle
L\left({\bf{p}}\right)\approx{\bf{\bar{p}}}_{n}^{T}{\bf{v}}_{\rm{w}}\left({{{\bf{\bar{p}}}}_{n}}\right)+{\left({{\bf{v}}_{\rm{w}}\left({{{\bf{\bar{p}}}}_{n}}\right)+\nabla{\bf{v}}_{\rm{w}}{{{\bf{\bar{p}}}}_{n}}}\right)^{T}}\left({{\bf{p}}-{{{\bf{\bar{p}}}}_{n}}}\right).$
(33)

Then, constraint (29b) can be recast to a convex form as

$\displaystyle L\left({\bf{p}}\right)+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}-1\leq 2{\varepsilon^{2}}.$ (34)

Thus, problem (27) can be reformulated as

$\displaystyle\mathop{\min}\limits_{{{\bf{p}}}}$
$\displaystyle~{}\phi\left({\bf{p}}\right)$ (35a)
$\displaystyle\quad{\rm{s}}{\rm{.t}}{\rm{.}}$
$\displaystyle~{}L\left({\bf{p}}\right)+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}-1\leq 2{\varepsilon^{2}},$ (35b)
$\displaystyle~{}{\bf{p}}^{T}{{\bf{1}}_{K}}=1,$ (35c)
$\displaystyle~{}{\bf{p}}^{T}\left({{\bf{x}}\odot{\bf{x}}}\right)\leq{P_{\rm{A}}},$
(35d) $\displaystyle~{}{{\bf{p}}}\geq 0,$ (35e)

which is convex.

Similarly, problem (35) can be efficiently solved by the inexact gradient
descent projection method. The details are omitted due to space limitation.

##  IV Signaling Design with Approximate Covert Rate Expression

Due to the expectation operation, the achievable rate in (7) does not have a
closed-form expression, and can only be computed numerically using the
approximate gradient descent method at the expense of high computational
complexity. To strike a balance between complexity and performance, we further
derive analytical upper bound and lower bound on the achievable rate in (7).

Lemma 1: An upper bound ${R_{\rm{b}}^{\text{U}}}$ on the rate ${R_{\rm{b}}}$
is given by

$\displaystyle{R_{\rm{b}}^{\text{U}}}=-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right),$
(36)

while a lower bound ${R_{\rm{b}}^{\text{L}}}$ on the rate ${R_{\rm{b}}}$ is
given as

$\displaystyle{R_{\rm{b}}^{\text{L}}}=$
$\displaystyle-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)\right|}^{\rm{2}}}}}{{2\sigma_{\rm{b}}^{2}}}}\right)$
$\displaystyle-\frac{1}{{\ln 2}}+1.$ (37)

Please find the derivation in Appendices C and D.

In this section, we adopt the upper bound and lower bound on the achievable
rate ${R_{\rm{b}}}$ in our following analysis.

###  IV-A Maximizing ${R_{\rm{b}}^{\text{U}}}$

In this subsection, we consider the upper bound on the achievable rate for Bob
${R_{\rm{b}}^{\text{U}}}$ as the objective function to find the optimal
probability of discrete constellation points set. Specifically, we study
beamforming design with the objective of maximizing ${R_{\rm{b}}^{\text{U}}}$,
subject to the covert transmission constraint, and the discrete constellation
set.

####  IV-A1 ${D_{\rm{U}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}}$

Finding the optimal probability of discrete constellation set can be
equivalently written as the following optimization problem

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}$
$\displaystyle~{}{R_{\rm{b}}^{\text{U}}}$ (38a)
$\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle~{}{D_{\rm{U}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}},$ (38b)
$\displaystyle~{}\Pr\left({X={x_{k}}}\right)={p_{k}}\geq{\rm{0}},$ (38c)
$\displaystyle~{}\sum\limits_{k=1}^{K}{{p_{k}}{{\left|{{x_{k}}}\right|}^{2}}}\leq{P_{\rm{A}}},$
(38d) $\displaystyle~{}\sum\limits_{k=1}^{K}{{p_{k}}}=1,k=1,...,K.$ (38e)

In order to solve problem (38), we first define the following variables

$\displaystyle{{\bf{r}}_{k}}\buildrel\Delta\over{=}\left[\begin{array}[]{l}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{1}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right),...,\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{K}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)\end{array}\right]^{T},$
(39b)
$\displaystyle{\bf{u}}\left({\bf{p}}\right)\buildrel\Delta\over{=}{\left[{{{\log}_{2}}{{\bf{p}}^{T}}{{\bf{r}}_{1}},...,{{\log}_{2}}{{\bf{p}}^{T}}{{\bf{r}}_{K}}}\right]^{T}}.$
(39c)

In this case, we can obtain

$\displaystyle{R_{\rm{b}}^{\text{U}}}=-{{\bf{p}}^{T}}{\bf{u}}\left({\bf{p}}\right).$
(40)

Therefore, problem (38) can be reformulated as follows

$\displaystyle\mathop{\min}\limits_{{\bf{p}}}$
$\displaystyle~{}{f_{\rm{U}}}\left({\bf{p}}\right)$ (41a)
$\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle~{}-{\log_{2}}{{\bf{p}}^{T}}{\bf{t}}\leq 2{\varepsilon^{2}},$
(41b) $\displaystyle~{}{{\bf{p}}^{T}}{{\bf{1}}_{K}}=1,$ (41c)
$\displaystyle~{}{\bf{p}}^{T}\left({{\bf{x}}\odot{\bf{x}}}\right)\leq{P_{\rm{A}}},$
(41d) $\displaystyle~{}{\bf{p}}\geq{\bf{0}},$ (41e)

where
${f_{\rm{U}}}\left({\bf{p}}\right)={{{\bf{p}}^{T}}{\bf{u}}}\left({\bf{p}}\right)$.
Since the Frank-Wolf method is an algorithm for solving linearly-constrained
problems, it makes a linear approximation of the objective function, obtains
the feasible descending direction by solving the linear programming, and
conducts a one-dimensional search in the feasible region along this direction.
Therefore, we will apply the Frank-Wolf method to solve the optimization
problem.

We use Taylor’s expansion to make a linear approximation of the objective
function ${f_{\rm{U}}}\left({\bf{p}}\right)$. The first order Taylor expansion
at ${{{\bf{p}}_{n}}}$ is as follows

$\displaystyle{f_{\rm{U}}}\left({\bf{p}}\right)\approx{\bf{p}}_{n}^{T}{\bf{u}}\left({{\bf{p}}_{n}}\right)+\nabla{f_{\rm{U}}}{\left({{{\bf{p}}_{i}}}\right)^{T}}\left({{\bf{p}}-{{\bf{p}}_{n}}}\right),$
(42a)
$\displaystyle\nabla{f_{\rm{U}}}\left({{{\bf{p}}_{n}}}\right)={\bf{u}}\left({{\bf{p}}_{n}}\right)+\nabla{\bf{u}}{{\bf{p}}_{n}},$
(42b)

where
$\nabla{\bf{u}}={\left[{\frac{{{{\bf{r}}_{1}}}}{{{{\bf{p}}_{n}^{T}}{{\bf{r}}_{1}}}},...,\frac{{{{\bf{r}}_{K}}}}{{{{\bf{p}}_{n}^{T}}{{\bf{r}}_{K}}}}}\right]_{K\times
K}}$, and ${{\bf{p}}_{n}}$ denotes the current iteration point. Then, we
reformulate the optimization problem of (41) as follows

$\displaystyle\mathop{\min}\limits_{{\bf{p}}}$
$\displaystyle~{}\nabla{f_{\rm{U}}}{\left({{{\bf{p}}_{n}}}\right)^{T}}{\bf{p}}$
(43) $\displaystyle{\rm{s}}{\rm{.t}}.$ $\displaystyle~{}\eqref{approxi_pro
01a},\eqref{approxi_pro 01b},\eqref{approxi_pro 01c},\eqref{approxi_pro 01d}.$

By applying the Frank-Wolf method, the detailed procedures for solving (40)
are summarized in Algorithm $3$. Note that ${\lambda_{n}}$ is the stepsize of
the $n$th iteration and ${{\bf{d}}_{n}}$ denotes the feasible descending
direction of the $n$th iteration.

Algorithm 3 : Solving (43) by Frank-Wolf method.

1:Initialization: Choose a feasible starting point ${{\bf{p}}_{0}}$, set
$\delta>0$ as the stopping parameter, let $n=0$;

2:While
$\left\|{\nabla{f}{{\left({{{\bf{p}}_{n}}}\right)}^{T}}{{\bf{d}}_{n}}}\right\|\leq\delta$;

3: Solve the linear programming problems (43) and obtain optimal solution
${{\bf{\bar{p}}}_{n}}$ ;

4: Construct the feasible descending direction
${{\bf{d}}_{n}}={\bf{\bar{p}}}_{n}-{\bf{p}}_{n}$;

5: Obtain optimal solution
${\lambda_{n}}=\mathop{\arg\min}\limits_{0\leq\lambda\leq
1}{f}\left({{{\bf{p}}_{n}}+\lambda{{\bf{d}}_{n}}}\right)$;

6: Let ${{\bf{p}}_{n+1}}={{\bf{p}}_{n}}+{\lambda_{n}}{{\bf{d}}_{n}}$,
$n\leftarrow n+1;$

7:end

8:Output ${{\bf{p}}_{n}}$

####  IV-A2 ${D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}}$

Furthermore, we consider the optimal probabilistic constellation shaping for
covert communications with covert constraint
${D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}}$, such that

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}$
$\displaystyle~{}{R_{\rm{b}}^{\text{U}}}$ (44a)
$\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle~{}{D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}},$ (44b)
$\displaystyle~{}\eqref{up1b},\eqref{up1c},\eqref{up1d},$

which is non-convex.

Similar to problem (29), by replacing constraint (44b) by constraint (34), we
can obtain the optimization problem

$\displaystyle\mathop{\min}\limits_{{\bf{p}}}$
$\displaystyle~{}\nabla{f_{U}}{\left({{{\bf{p}}_{n}}}\right)^{T}}{\bf{p}}$
(45a) $\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle~{}L\left({\bf{p}}\right)+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}-1\leq 2{\varepsilon^{2}},$ (45b)
$\displaystyle~{}\eqref{approxi_pro 01a},\eqref{approxi_pro
01b},\eqref{approxi_pro 01c},\eqref{approxi_pro 01d}.$

Then, we apply the Frank-Wolf method to solve problem (38), and the details
are omitted since the corresponding algorithm is similar to Algorithm 3.

###  IV-B Maximizing ${R_{\rm{b}}^{\text{L}}}$

In this subsection, we further study the lower bound beamforming design for
covert communication by maximizing the lower bound ${R_{\rm{b}}^{\text{L}}}$,
while satisfying the covert transmission requirement and the discrete
constellation set with $K$.

####  IV-B1 ${D_{\rm{U}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}}$

Under the covert constraint
${D_{\rm{U}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}}$, the optimal probabilistic constellation shaping for
covert communications is formulated as follows

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}$
$\displaystyle~{}{R_{\rm{b}}^{\text{L}}}$ (46a)
$\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle~{}{D_{\rm{U}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}},$ (46b)
$\displaystyle~{}\Pr\left({X={x_{k}}}\right)={p_{k}}\geq{\rm{0}},$ (46c)
$\displaystyle~{}\sum\limits_{k=1}^{K}{{p_{k}}{{\left|{{x_{k}}}\right|}^{2}}}\leq{P_{\rm{A}}},$
(46d) $\displaystyle~{}\sum\limits_{k=1}^{K}{{p_{k}}}=1,k=1,...,K.$ (46e)

To solve the problem, we define the following equations

$\displaystyle{{\bf{s}}_{{\rm{b}},k}}\buildrel\Delta\over{=}\left[\begin{array}[]{l}\exp\left({-\frac{{{{\left|{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{1}}}\right)\right|}^{\rm{2}}}}}{{2\sigma_{\rm{b}}^{2}}}}\right),...,\exp\left({-\frac{{{{\left|{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{K}}}\right)\right|}^{\rm{2}}}}}{{2\sigma_{\rm{b}}^{2}}}}\right)\end{array}\right]^{T},$
(47b)
$\displaystyle{{\bf{v}}_{\rm{b}}}\left({\bf{p}}\right)\buildrel\Delta\over{=}{\left[{{{\log}_{2}}{{\bf{p}}^{T}}{{\bf{s}}_{{\rm{b}},1}},...,{{\log}_{2}}{{\bf{p}}^{T}}{{\bf{s}}_{{\rm{b}},K}}}\right]^{T}},$
(47c)

and then the lower bound ${R_{\rm{b}}^{\text{L}}}$ can be transformed to

$\displaystyle{R_{\rm{b}}^{\text{L}}}=-{{\bf{p}}^{T}}{{\bf{v}}_{\rm{b}}}\left({\bf{p}}\right)-\frac{1}{{\ln
2}}+1.$ (48)

Thus, the covert optimization problem is recast as follows

$\displaystyle\mathop{\max}\limits_{\bf{p}}$
$\displaystyle-{{\bf{p}}^{T}}{{\bf{v}}_{\rm{b}}}\left({\bf{p}}\right)$ (49a)
$\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle~{}-{\log_{2}}{{\bf{p}}^{T}}{\bf{t}}\leq 2{\varepsilon^{2}},$
(49b) $\displaystyle~{}{{\bf{p}}^{T}}{{\bf{1}}_{K}}=1,$ (49c)
$\displaystyle~{}{\bf{p}}^{T}\left({{\bf{x}}\odot{\bf{x}}}\right)\leq{P_{\rm{A}}},$
(49d) $\displaystyle~{}{\bf{p}}\geq{\bf{0}},$ (49e)

Similar to problem (38), we can apply the Frank-Wolf method to solve problem
(46), and the details are omitted.

####  IV-B2 ${D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}}$

With the covert constraint
${D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}}$, the optimal probabilistic constellation shaping for
covert communications is given as

$\displaystyle\mathop{\max}\limits_{\left\\{{{p_{k}}}\right\\}}$
$\displaystyle~{}{R_{\rm{b}}^{\text{L}}}$ (50a)
$\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle{D_{\rm{U}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq
2{\varepsilon^{2}},$ (50b)
$\displaystyle\eqref{lowb},\eqref{lowc},\eqref{lowd}.$

By replacing the constraint (50b) with (34), we obtain the optimization
problem

$\displaystyle\mathop{\max}\limits_{\bf{p}}$
$\displaystyle~{}-{{\bf{p}}^{T}}{{\bf{v}}_{\rm{b}}}\left({\bf{p}}\right)$
(51a) $\displaystyle{\rm{s}}{\rm{.t}}.$
$\displaystyle~{}L\left({\bf{p}}\right)+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}-1\leq 2{\varepsilon^{2}},$ (51b)
$\displaystyle~{}\eqref{low02b},\eqref{low02c},\eqref{low02d}.$

Similar to problem (44), we can apply a similar method to solve problem (50).

##  V Numerical Results

In this section, we present and discuss numerical results to assess the
performance of the proposed probabilistic constellation shaping designs. In
our simulations, the discrete constellation input is QAM modulation, the total
transmit power of Alice is ${P_{\rm{A}}}=10\rm{W}$, the noise variance of
Willie is $\sigma_{\rm{w}}^{2}=1\rm{W}$, and the noise variance of Bob is
$\sigma_{\rm{b}}^{\rm{2}}{\rm{=}}\frac{{{P_{\rm{A}}}}}{{{{10}^{{{{\rm{SNR}}}\mathord{\left/{\vphantom{{{\rm{SNR}}}{10}}}\right.\kern-1.2pt}{10}}}}}}$.
Moreover, we assume that all channels experience Rayleigh flat fading, and
$\sigma_{1}=\sigma_{2}=1$ [48].

We first compare the proposed optimal probabilistic constellation shaping
design with the equiprobable design, starting from the empirical CDF of KL
divergence and the rate comparison.

(a)

(b)

Figure 2:  The empirical CDF of a)
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and (b)
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ with the covertness
threshold $2{\varepsilon^{2}}=0.02$ for the proposed optimal probabilistic
constellation shaping design and the equiprobable design.

Fig. 2 shows the empirical CDF of the achieved
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$, respectively, for both
the proposed optimal probabilistic constellation shaping design and the
equiprobable design with ${\rm{SNR}}=10\rm{dB}$ and $K=8$, where the
covertness threshold is $2{\varepsilon^{2}}=0.02$, i.e.,
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq 0.02$ and
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq 0.02$. As observed
from Fig. 2, the proposed optimal probabilistic constellation shaping design
satisfies the covertness constraint. On the other hand, the equiprobable
design cannot satisfy the covertness constraints.

(a)

(b)

Figure 3:  The covert rate of a)
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and (b)
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ for the proposed optimal
probabilistic constellation shaping design and the equiprobable design.

Fig. 3 (a) and (b) show the achievable covert rate of Bob ${R_{\rm{b}}}$ for
the proposed optimal probabilistic constellation shaping design and the
equiprobable design with $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$
and $D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$, respectively. We
observe that the optimal probabilistic constellation shaping design is
superior when the signal-to-noise ratio is low. Therefore, in practical
applications, the proposed optimal probabilistic constellation shaping design
is advantageous in medium and low SNR, and the equiprobable design is suitable
for high SNR.

Next, we evaluate the performance of the proposed optimal probabilistic
constellation shaping design.

(a)

(b)

Figure 4:  The optimal probability of discrete constellation points with (a)
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and (b)
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ for the proposed optimal
probabilistic constellation shaping design.

Fig. 4 shows the optimal probability distribution of input
$\left\\{{p_{{i,j}}}\right\\}$ with $\rm{SNR=12dB}$ of the proposed optimal
probabilistic constellation shaping design, for
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ in Fig. 4 (a), and for
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ in Fig. 4 (b). As it can
be seen from Fig. 4, for the proposed optimal probabilistic constellation
shaping design, the optimal probability distribution is not equiprobable, and
the symmetrical points have equal probabilities. Specifically, when the number
of discrete constellation points is sixteen, the probability of each
constellation point for $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$
and $D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ is given in Table I.
From the table, we can clearly see that when the coordinates of the
constellation points are symmetrical, their probabilities are the same, and
vice versa.

TABLE I: The optimal probability distribution of the optimal probabilistic constellation shaping design for $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and $D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ with $\rm{SNR=12dB}$ and $K=16$ ${\rm{Im}}\left\\{{x_{{i,j}}}\right\\}$ ${\rm{Re}}\left\\{{x_{{i,j}}}\right\\}$ $\left\\{{p_{{i,j}}}\right\\}$ | $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ | $D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$
---|---|---
$-3$ |  $-1$ |  $1$ |  $3$ |  $-3$ |  $-1$ |  $1$ |  $3$
-3 |  $0.0484$ |  $0.0524$ |  $0.0524$ |  $0.0484$ |  $0.0454$ |  $0.0502$ |  $0.0502$ |  $0.0454$
-1 |  $0.0524$ |  $0.0968$ |  $0.0968$ |  $0.0524$ |  $0.0502$ |  $0.1041$ |  $0.1041$ |  $0.0502$
1 |  $0.0524$ |  $0.0968$ |  $0.0968$ |  $0.0524$ |  $0.0502$ |  $0.1041$ |  $0.1041$ |  $0.0502$
3 |  $0.0484$ |  $0.0552$ |  $0.0552$ |  $0.0484$ |  $0.0454$ |  $0.0502$ |  $0.0502$ |  $0.0454$
Figure 5:  The achievable covert rate of Bob versus $\rm{SNR}$ with different
number of points $K$ for $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$.

Fig. 5 considers the proposed optimal probabilistic constellation shaping
design and depicts, the achievable covert rate of Bob ${R_{\rm{b}}}$ versus
$\rm{SNR}$ with different number of constellation points $K=2,4,8,16$ for
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$. It can be seen from Fig.
5 that when ${\rm{SNR}}$ increases, the covert rate of Bob ${R_{\rm{b}}}$
increases. In addition, we observe that larger number of points $K$ results in
higher covert rate of Bob ${R_{\rm{b}}}$, especially for high SNRs. Thus, as
the modulation order increases, the rate of covert communication increases.

Fig. 6 shows the covert rate versus $\varepsilon$ for the proposed optimal
probabilistic constellation shaping design for
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$, where $K=8$,
${\rm{SNR}}=6{\rm{dB}}$. It can be observed from the figure that as
$\varepsilon$ increases, the covert constraint becomes loose, resulting in a
covert rate increase. The rate of the covert constraint
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq 2{\varepsilon^{2}}$ is
higher than that of the case
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq 2{\varepsilon^{2}}$.
This is because $D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ is less
than $D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ for the same
probability distribution, and thus
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq 2{\varepsilon^{2}}$ is
more stringent than $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq
2{\varepsilon^{2}}$.

Figure 6:  The covert rate versus $\varepsilon$ for the proposed optimal
probabilistic constellation shaping design with the cases of
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$.

Finally, we compare the performance and complexity of three objective
functions.

(a)

(b)

Figure 7:  The achievable rate of Bob versus $\rm{SNR}$ for (a)
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and (b)
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ with the covertness
threshold $2{\varepsilon^{2}}=0.1$.

Fig. 7 depicts the achievable rate of Bob ${R_{\rm{b}}}$ with the proposed
optimal probabilistic constellation shaping design, as well as the objective
functions ${R_{\rm{b}}^{\text{L}}}$ and ${R_{\rm{b}}^{\text{U}}}$ versus the
${\rm{SNR}}$ for the case of
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ and
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$, respectively. It can be
observed that the mutual information of Bob ${R_{\rm{b}}}$ increases as the
${\rm{SNR}}$ increases, while ${R_{\rm{b}}}$ of the proposed optimal
probabilistic constellation shaping design is between the objective functions
${R_{\rm{b}}^{\text{L}}}$ and ${R_{\rm{b}}^{\text{U}}}$, and the objective
function ${R_{\rm{b}}^{\text{U}}}$ is higher than the objective function
${R_{\rm{b}}^{\text{L}}}$.

From Fig. 7 we observe that the proposed optimal probabilistic constellation
shaping design is between the objective functions ${R_{\rm{b}}^{\text{L}}}$
and ${R_{\rm{b}}^{\text{U}}}$. Here, we compare the computational complexity
of the three designs by computational time in Table II, and all simulations of
the three methods are performed using MATLAB 2016b with 2.30GHz, 2.29GHz dual
CPUs and a 128GB RAM, where $K=8$. Specifically, Table III shows that the
computational time of objective functions ${R_{\rm{b}}}$,
${R_{\rm{b}}^{\text{L}}}$ and ${R_{\rm{b}}^{\text{U}}}$ for the covert
constraint condition $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ is
55.04, 2.845 and 3.012 seconds, respectively. The computational times of the
latter two cases is approximately 95 percent shorter than that of the
probabilistic constellation shaping design. Under the covert constraint
condition $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ the
computational time of objective functions ${R_{\rm{b}}}$,
${R_{\rm{b}}^{\text{L}}}$ and ${R_{\rm{b}}^{\text{U}}}$ is 107.77, 3.219 and
3.438 seconds, respectively. The computational times of the latter two cases
is approximately improved by 97 percent compared to that of the probabilistic
constellation shaping design. Moreover, the computational time of the design
for $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ is less than that of
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$.

TABLE II:  Computational time comparison among the objective functions ${R_{\rm{b}}}$, ${R_{\rm{b}}^{\text{L}}}$ and ${R_{\rm{b}}^{\text{U}}}$ Constraint Objective Time/second | ${R_{\rm{b}}}$ | ${R_{\rm{b}}^{\text{L}}}$ | ${R_{\rm{b}}^{\text{U}}}$
---|---|---|---
$D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ | $55.04$ | $2.845$ | $3.012$
$D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ | $107.77$ | $3.219$ | $3.438$

##  VI Conclusions

In this paper, we propose an optimal probabilistic constellation shaping
design for covert communications, where Alice covertly sends a message to Bob
while avoiding being discovered by Willie. We derive the achievable rate
expressions of the covert communications system, and we study the covert rate
maximization problem via optimizing the constellation distribution. In
addition, to strike a balance between the computational complexity and the
transmission performance, we further develop a framework that maximizes the
upper and lower bounds of the achievable rate. Numerical results quantify the
gains of the proposed beamformers design over state-of-the-art schemes in
terms of the achievable covert rate.

##  Appendix A Derivation of the formulation (18)

The upper bound
${D_{{\rm{U}}}}\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ on the KL
divergence $D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)$ is derived as
follows

$\displaystyle
D\left({{p_{y,0}}\left\|{{p_{y,1}}}\right.}\right)\leq-\frac{1}{{\ln 2}}$
$\displaystyle-{\log_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}\exp\left({-{\mathbb{E}_{{z_{\rm{w}}}}}\left\\{{\frac{{{{\left|{{z_{\rm{w}}}-{g_{\rm{w}}^{\rm{*}}}{x_{k}}}\right|}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right\\}}\right)$
(52a) $\displaystyle=-\frac{1}{{\ln
2}}-{\log_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}\exp\left({-\frac{{{\mathbb{E}_{{z_{\rm{w}}}}}\left\\{{z_{\rm{w}}^{2}}\right\\}+{{\left|{g_{\rm{w}}^{\rm{*}}{x_{k}}}\right|}^{\rm{2}}}}}{{\sigma_{\rm{w}}^{2}}}}\right)$
(52b) $\displaystyle=-\frac{1}{{\ln
2}}-{\log_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}\exp\left({-1-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}{x_{k}}}\right|}^{\rm{2}}}}}{{\sigma_{\rm{w}}^{2}}}}\right)$
(52c) $\displaystyle=-\frac{1}{{\ln 2}}-\left({-\frac{1}{{\ln
2}}+{{\log}_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}{x_{k}}}\right|}^{\rm{2}}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right)$
(52d)
$\displaystyle=-{\log_{2}}\sum\limits_{k=1}^{K}{{p_{k}}}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}{x_{k}}}\right|}^{\rm{2}}}}}{{\sigma_{\rm{w}}^{2}}}}\right),$
(52e)

where inequality (52a) holds due to Jensen’s Inequality.

##  Appendix B Derivation of the formulation (III-B)

The upper bound
${D_{{\rm{U}}}}\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ of the KL
divergence $D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)$ is given as

$\displaystyle
D\left({{p_{y,1}}\left\|{{p_{y,0}}}\right.}\right)\leq\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}$
$\displaystyle\times{\mathbb{E}_{{z_{\rm{w}}}}}\left\\{{\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}\left({{x_{k}}-{x_{j}}}\right)+{z_{\rm{w}}}}\right|}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right\\}+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}$ (53a)
$\displaystyle=\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}{\mathbb{E}_{{z_{{\rm{w}},R}}}}\left\\{{\exp\left({-\frac{{{{\left({{c_{R}}+{z_{{{\rm{w}},R}}}}\right)}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right\\}$
$\displaystyle\times{\mathbb{E}_{{z_{{\rm{w}},I}}}}\left\\{{\exp\left({-\frac{{{{\left({{c_{I}}+{z_{{\rm{w}},I}}}\right)}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right\\}+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}$ (53b)
$\displaystyle=\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\left[{\int\limits_{{\rm{-}}\infty}^{\infty}{\frac{{d{z_{{\rm{w}},R}}}}{{\sqrt{\pi}{\sigma_{\rm{w}}}}}\exp\left({-\frac{{{{\left({{c_{R}}+{z_{{\rm{w}},R}}}\right)}^{2}}}}{{\sigma_{\rm{w}}^{2}}}}\right.}}\right.$
$\displaystyle\left.{\left.{+\frac{{z_{{\rm{w}},R}^{2}}}{{\sigma_{\rm{w}}^{2}}}}\right)}\right]\times\left[{\int\limits_{{\rm{-}}\infty}^{\infty}{\frac{{d{z_{{\rm{w}},I}}}}{{\sqrt{\pi}{\sigma_{\rm{w}}}}}\exp\left({-\frac{{{{\left({{c_{R}}+{z_{{\rm{w}},I}}}\right)}^{2}}+z_{{\rm{w}},I}^{2}}}{{\sigma_{\rm{w}}^{2}}}}\right)}}\right]$
$\displaystyle+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}$ (53c)
$\displaystyle=\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\left[{\frac{{\rm{1}}}{{\sqrt{\pi}{\sigma_{\rm{w}}}}}\frac{1}{2}\sqrt{\frac{{\sigma_{\rm{w}}^{2}\pi}}{2}}\exp\left({-\frac{{c_{R}^{2}}}{{2\sigma_{\rm{w}}^{2}}}}\right)2}\right]$
$\displaystyle\times\left[{\frac{{\rm{1}}}{{\sqrt{\pi}{\sigma_{\rm{w}}}}}\frac{1}{2}\sqrt{\frac{{\sigma_{\rm{w}}^{2}\pi}}{2}}\exp\left({-\frac{{c_{I}^{2}}}{{2\sigma_{\rm{w}}^{2}}}}\right)2}\right]+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}$ (53d)
$\displaystyle=\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\frac{1}{2}\exp\left({-\frac{{c_{R}^{2}+c_{I}^{2}}}{{2\sigma_{\rm{w}}^{2}}}}\right)+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}$ (53e)
$\displaystyle=\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\frac{1}{2}+\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)\right|}^{\rm{2}}}}}{{2\sigma_{\rm{w}}^{2}}}}\right)$
$\displaystyle+\frac{1}{{\ln
2}}+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}$ (53f)
$\displaystyle=\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{g_{\rm{w}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)\right|}^{\rm{2}}}}}{{2\sigma_{\rm{w}}^{2}}}}\right)+\frac{1}{{\ln
2}}$
$\displaystyle+\frac{\left|{{{{g_{\rm{w}}^{2}}}}{P_{\rm{A}}}}\right|}{{\left({\ln
2}\right)\sigma_{\rm{w}}^{2}}}-1,$ (53g)

where inequality (53a) is true due to Jensen’s Inequality. Equality (53b)
holds because of the definitions
${z_{{\rm{w}},R}}\buildrel\Delta\over{=}{\rm{Re}}\left\\{{{z_{\rm{w}}}}\right\\}$
and
${z_{{\rm{w}},I}}\buildrel\Delta\over{=}{\rm{Im}}\left\\{{{z_{\rm{w}}}}\right\\}$,
${z_{{\rm{w}},R}},{z_{{\rm{w}},I}}\sim{\cal{N}}\left({0,\frac{1}{2}\sigma_{\rm{w}}^{2}}\right)$,
and ${c_{R}}\buildrel\Delta\over{=}{\mathop{\rm
Re}\nolimits}\left({{g_{\rm{w}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right)$,
${c_{I}}\buildrel\Delta\over{=}{\mathop{\rm
Im}\nolimits}\left({{g_{\rm{w}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right)$.

##  Appendix C Derivation of the formulation (36)

The upper bound ${{R_{\rm{b}}^{{\rm{U}}}}}$ of the covert rate ${R_{\rm{b}}}$
is derived as follows

$\displaystyle{R_{\rm{b}}}\leq-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}$
$\displaystyle\times\exp\left({-{\mathbb{E}_{{z_{\rm{b}}}}}\left\\{{\frac{{{{\left|{g_{\rm{b}}^{\rm{*}}\left({{x_{k}}-{x_{j}}}\right)+{z_{\rm{b}}}}\right|}^{2}}}}{{\sigma_{b}^{2}}}}\right\\}}\right)-\frac{1}{{\ln
2}}$ (54a)
$\displaystyle=-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}-\frac{{{\mathbb{E}_{{z_{\rm{b}}}}}\left\\{{\left|z_{\rm{b}}\right|^{2}}\right\\}}}{{\sigma_{\rm{b}}^{2}}}}\right)$
$\displaystyle-\frac{1}{{\ln 2}}$ (54b) $\displaystyle=-\frac{1}{{\ln
2}}-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}-1}\right)$
(54c)
$\displaystyle=-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)\exp\left({-1}\right)$
$\displaystyle-\frac{1}{{\ln 2}}$ (54d)
$\displaystyle=-\sum\limits_{k=1}^{K}{{p_{k}}}\left({{{\log}_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)-\frac{1}{{\ln
2}}}\right)$ $\displaystyle-\frac{1}{{\ln 2}}$ (54e)
$\displaystyle=-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right|}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right),$
(54f)

where inequality (54a) is true due to Jensen’s Inequality.

##  Appendix D Derivation of the formulation (IV)

The lower bound ${{R_{\rm{b}}^{{\rm{L}}}}}$ on the covert rate ${R_{\rm{b}}}$
is given as

$\displaystyle{R_{\rm{b}}}\geq-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}$
$\displaystyle\times{\mathbb{E}_{{z_{\rm{b}}}}}\left\\{{\exp\left({-\frac{{{{\left({g_{\rm{b}}^{\rm{*}}\left({{x_{k}}-{x_{j}}}\right)+{z_{\rm{b}}}}\right)}^{2}}}}{{\sigma_{\rm{b}}^{2}}}}\right)}\right\\}-\frac{1}{{\ln
2}}$ (55a)
$\displaystyle=-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\left[{\int\limits_{{\rm{-}}\infty}^{\infty}{\frac{{{d{z_{{\rm{b}},R}}}}}{{\sqrt{\pi}{\sigma_{\rm{b}}}}}\exp\left({-\frac{{{{\left({{a_{R}}+{z_{{\rm{b}},R}}}\right)}^{2}}+z_{{\rm{b}},R}^{2}}}{{\sigma_{\rm{b}}^{2}}}}\right)}}\right]$
$\displaystyle\times\left[{\int\limits_{{\rm{-}}\infty}^{\infty}{\frac{{{d{z_{{\rm{b}},I}}}}}{{\sqrt{\pi}{\sigma_{\rm{b}}}}}\exp\left({-\frac{{{{\left({{a_{R}}+{z_{{\rm{b}},I}}}\right)}^{2}}+z_{{\rm{b}},I}^{2}}}{{\sigma_{\rm{b}}^{2}}}}\right)}}\right]-\frac{1}{{\ln
2}}$ (55b)
$\displaystyle=-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\left[{\frac{{\rm{1}}}{{\sqrt{\pi}{\sigma_{\rm{b}}}}}\frac{1}{2}\sqrt{\frac{{\sigma_{\rm{b}}^{2}\pi}}{2}}\exp\left({-\frac{{a_{R}^{2}}}{{2\sigma_{\rm{b}}^{2}}}}\right)2}\right]$
$\displaystyle\times\left[{\frac{{\rm{1}}}{{\sqrt{\pi}{\sigma_{\rm{b}}}}}\frac{1}{2}\sqrt{\frac{{\sigma_{\rm{b}}^{2}\pi}}{2}}\exp\left({-\frac{{a_{I}^{2}}}{{2\sigma_{\rm{b}}^{2}}}}\right)2}\right]-\frac{1}{{\ln
2}}$ (55c) $\displaystyle=-\frac{1}{{\ln
2}}-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\frac{1}{2}\exp\left({-\frac{{a_{R}^{2}+a_{I}^{2}}}{{2\sigma_{\rm{b}}^{2}}}}\right)$
(55d) $\displaystyle=-\frac{1}{{\ln
2}}+1-\sum\limits_{k=1}^{K}{{p_{k}}}{\log_{2}}\sum\limits_{j=1}^{K}{{p_{j}}}\exp\left({-\frac{{{{\left|{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)\right|}^{\rm{2}}}}}{{2\sigma_{\rm{b}}^{2}}}}\right),$
(55e)

where inequality (55a) is true due to Jensen’s inequality, equality (55b)
holds because of the definitions
${z_{{\rm{b}},R}}\buildrel\Delta\over{=}{\rm{Re}}\left\\{{{z_{\rm{b}}}}\right\\}$,
${z_{{\rm{b}},I}}\buildrel\Delta\over{=}{\rm{Im}}\left\\{{{z_{\rm{b}}}}\right\\}$,
${z_{{\rm{b}},R}},{z_{{\rm{b}},I}}\sim{\cal{N}}\left({0,\frac{1}{2}\sigma_{\rm{b}}^{2}}\right)$,
and ${a_{R}}\buildrel\Delta\over{=}{\mathop{\rm
Re}\nolimits}\left({{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right)$,
${a_{I}}\buildrel\Delta\over{=}{\mathop{\rm
Im}\nolimits}\left({{g_{\rm{b}}^{\rm{*}}}\left({{x_{k}}-{x_{j}}}\right)}\right)$.

## References

  * [1] V. L. Nguyen, P. C. Lin, B. C. Cheng, R. H. Hwang, and Y. D. Lin,  “Security and privacy for 6G: A survey on prospective technologies and challenges,”  IEEE Commun. Surveys Tuts., vol. 23, no. 4, pp. 2384–2428, 4th Quart. 2021.
  * [2] X. Chen, D. W. K. Ng, W. H. Gerstacker, and H. H. Chen,  “A survey on multiple-antenna techniques for physical layer security,”  IEEE Commun. Surveys Tuts., vol. 19, no. 2, pp. 1027–1053, 2nd Quart. 2017.
  * [3] M. Bloch and J. Barros,  Physical-Layer Security: From Information Theory to Security Engineering,  U.K.: Cambridge Univ., 2011.
  * [4] H. Wang, T. Zheng, J. Yuan, D. Towsley, and M. H. Lee,  “Physical layer security in heterogeneous cellular networks,”  IEEE Trans. Commun., vol. 64, no. 3, pp. 1204–1219, Mar. 2016.
  * [5] H. Du, D. Niyato, J. Kang, D. I. Kim, and C. Miao,  “Optimal targeted advertising strategy for secure wireless edge metaverse,”  arXiv preprint arXiv:2202.00973, Nov. 2021.
  * [6] H. Du, J. Kang, T. D. Niyato, J. Zhang, and D. I. Kim,  “Reconfigurable intelligent surface-aided joint radar and covert communications: Fundamentals, optimization, and challenges,”  IEEE Vehicular Technology Magazine, pp. 2–12, May 2022.
  * [7] H. Du, D. Niyato, Y.-a. Xie, Y. Cheng, J. Kang, and D. I. Kim,  “Performance analysis and optimization for jammer-aided multi-antenna uav covert communication,”  arXiv preprint arXiv:2202.00973, Feb. 2022.
  * [8] Y.-A. Xie, J. Kang, D. Niyato, N. T. T. Van, N. C. Luong, Z. Liu, and H. Yu,  “Securing federated learning: A covert communication-based approach,”  IEEE Netw., pp. 1–8, Jul. 2022.
  * [9] S. Yan, Y. Cong, S. V. Hanly, and X. Zhou,  “Gaussian signalling for covert communications,”  IEEE Trans. Wireless Commun., vol. 18, no. 7, pp. 3542–3553, Jul. 2019.
  * [10] M. R. Bloch,  “Covert communication over noisy channels: A resolvability perspective,”  IEEE Trans. Inf. Theory, vol. 62, no. 5, pp. 2334–2354, May. 2016\.
  * [11] M. K. Simon, J. K. Omura, R. A. Scholtz, and B. K. Levitt,  Spread Spectrum Communications Handbook,  McGaw-Hill, Apr. 1994.
  * [12] B. A. Bash, D. Goeckel, and D. Towsley,  “Limits of reliable communication with low probability of detection on AWGN channels,”  IEEE J. Sel. Areas Commun., vol. 31, no. 9, pp. 1921–1930, Sept. 2013.
  * [13] L. Wang, W. Wornell, and L. Zheng,  “Fundamental limits of communication with low probability of detection,”  IEEE Trans. Inf. Theory, vol. 62, no. 6, pp. 3493–3503, Jun. 2016\.
  * [14] H. Wu, X. Liao, Y. Dang, Y. Shen, and X. Jiang,  “Limits of covert communication on two-hop AWGN channels,”  in Int. Conf. Netw. Netw. Appl., pp. 42–47, Oct. 2017.
  * [15] P. H. Che, M. Bakshi, and S. Jaggi,  “Reliable deniable communication: Hiding messages in noise,”  in 2013 IEEE International Symposium on Information Theory, pp. 2945–2949, 2013.
  * [16] K. S. K. Arumugam and M. R. Bloch,  “Embedding covert information in broadcast communications,”  IEEE Trans. Inf. Forensics Security, vol. 14, no. 10, pp. 2787–2801, Oct. 2019.
  * [17] K. S. K. Arumugam and M. R. Bloch,  “Covert communication over a $k$ -user multiple-access channel,”  IEEE Trans. Inf. Theory, vol. 65, no. 11, pp. 7020–7044, Nov. 2019\.
  * [18] K. H. Cho and S. H. Lee,  “Treating interference as noise is optimal for covert communication over interference channels,”  IEEE Trans. Inf. Forensics Security, vol. 16, pp. 322–332, Jul. 2021.
  * [19] B. A. Bash, D. Goeckel, and D. Towsley,  “Covert communication gains from adversary’s ignorance of transmission time,”  IEEE Trans. Wireless Commun., vol. 15, no. 12, pp. 8394–8405, Dec. 2016.
  * [20] C. Wang, Z. Li, J. Shi, and D. W. K. Ng,  “Intelligent reflecting surface-assisted multi-antenna covert communications: Joint active and passive beamforming optimization,”  IEEE Trans. Commun., vol. 69, no. 6, pp. 3984–4000, Jun. 2021.
  * [21] X. Chen, T. X. Zheng, L. Dong, M. Lin, and J. Yuan,  “Enhancing MIMO covert communications via intelligent reflecting surface,”  IEEE Wireless Commun. Lett., 2021.
  * [22] J. Si, Z. Li, Y. Zhao, J. Cheng, L. Guan, J. Shi, and N. AL-Dhahir,  “Covert transmission assisted by intelligent reflecting surface,”  IEEE Trans. Commun., vol. 69, no. 8, pp. 5394–5408, Aug. 2021.
  * [23] W. Gao, Y. Chen, C. Han, and Z. Chen,  “Distance-adaptive absorption peak modulation (DA-APM) for terahertz covert communications,”  IEEE Trans. Wireless Commun., vol. 20, no. 3, pp. 2064–2077, Mar. 2021.
  * [24] Z. Liu, J. Liu, Y. Zeng, and J. Ma,  “Covert wireless communication in iot network: From AWGN channel to THz band,”  IEEE IoT-J, vol. 7, no. 4, pp. 3378–3388, Apr. 2020.
  * [25] T. X. Zheng, Z. Yang, C. Wang, Z. Li, J. Yuan, and X. Guan,  “Wireless covert communications aided by distributed cooperative jamming over slow fading channels,”  IEEE Trans. Wireless Commun., vol. 20, no. 11, pp. 7026–7039, Nov. 2021.
  * [26] T. V. Sobers, B. A. Bash, S. Guha, D. Towsley, and D. Goeckel,  “Covert communication in the presence of an uninformed jammer,”  IEEE Trans. Wireless Commun., vol. 16, no. 9, pp. 6193–6206, Sept. 2017.
  * [27] K. Li, P. A. Kelly, and D. Goeckel,  “Optimal power adaptation in covert communication with an uninformed jammer,”  IEEE Trans. Wireless Commun., vol. 19, no. 5, pp. 3463–3473, May. 2020.
  * [28] O. Shmuel, A. Cohen, and O. Gurewitz,  “Multi-antenna jamming in covert communication,”  IEEE Trans. Commun., vol. 69, no. 7, pp. 4644–4658, Jul. 2021.
  * [29] K. Li, T. Sobers, D. Towsley, and D. Goeckel,  “Covert communication in continuous-time systems in the presence of a jammer,”  IEEE Trans. Wireless Commun., 2021.
  * [30] J. Zhang, M. Li, S. Yan, C. Liu, X. Chen, M. Zhao, and P. Whiting,  “Joint beam training and data transmission design for covert millimeter-wave communication,”  IEEE Trans. Inf. Forensics Security, vol. 16, pp. 2232–2245, Jan. 2021.
  * [31] D. Goeckel, B. Bash, S. Guha, and D. Towsley,  “Covert communications when the warden does not know the background noise power,”  IEEE Commun. Lett., vol. 20, no. 2, pp. 236–239, Feb. 2016.
  * [32] S. Lee, R. J. Baxley, M. A. Weitnauer, and B. Walkenhorst,  “Achieving undetectable communication,”  IEEE J. Sel. Topics Signal Process., vol. 9, no. 7, pp. 1195–1205, 2015.
  * [33] B. He, S. Yan, X. Zhou, and V. K. N. Lau,  “On covert communication with noise uncertainty,”  IEEE Commun. Lett., vol. 21, no. 4, pp. 941–944, Apr. 2017.
  * [34] K. Shahzad and X. Zhou,  “Covert wireless communications under quasi-static fading with channel uncertainty,”  IEEE Trans. Inf. Forensics Security, vol. 16, pp. 1104–1116, Oct. 2021.
  * [35] Z. Cheng, J. Si, Z. Li, L. Guan, Y. Zhao, D. Wang, J. Cheng, and N. Al-Dhahir,  “Covert surveillance via proactive eavesdropping under channel uncertainty,”  IEEE Trans. Commun., vol. 69, no. 6, pp. 4024–4037, Jun. 2021.
  * [36] T. Zheng, H. Wang, D. W. K. Ng, and J. Yuan,  “Multi-antenna covert communications in random wireless networks,”  IEEE Trans. Wireless Commun., vol. 18, no. 3, pp. 1974–1987, 2019\.
  * [37] S. Ma, Y. Zhang, H. Li, S. Lu, N. Al-Dhahir, S. Zhang, and S. Li,  “Robust beamforming design for covert communications,”  IEEE Trans. Inf. Forensics Security, vol. 16, pp. 3026–3038, Apr. 2021.
  * [38] S. Yan, B. He, X. Zhou, Y. Cong, and A. L. Swindlehurst,  “Delay-intolerant covert communications with either fixed or random transmit power,”  IEEE Trans. Inf. Forensics Security, vol. 14, no. 1, pp. 129–140, Jan. 2019.
  * [39] F. Shu, T. Xu, J. Hu, and S. Yan,  “Delay-constrained covert communications with a full-duplex receiver,”  IEEE Wireless Commun. Lett., vol. 8, no. 3, pp. 813–816, Jun. 2019\.
  * [40] C. Wang, Z. Li, and D. W. K. Ng,  “Covert rate optimization of millimeter wave full-duplex communications,”  IEEE Trans. Wireless Commun., 2021.
  * [41] R. Sun, B. Yang, S. Ma, Y. Shen, and X. Jiang,  “Covert rate maximization in wireless full-duplex relaying systems with power control,”  IEEE Trans. Commun., vol. 69, no. 9, pp. 6198–6212, Sept. 2021\.
  * [42] M. Zheng, A. Hamilton, and C. Ling,  “Covert communications with a full-duplex receiver in non-coherent rayleigh fading,”  IEEE Trans. Commun., vol. 69, no. 3, pp. 1882–1895, Mar. 2021.
  * [43] J. Hu, S.Yan, F. Shu, and J. Wang,  “Covert transmission with a self-sustained relay,”  IEEE Trans. Wireless Commun., vol. 18, no. 8, pp. 4089–4102, Aug. 2019.
  * [44] A. Sheikholeslami, M. Ghaderi, D. Towsley, B. A. Bash, S. Guha, and D. Goeckel,  “Multi-hop routing in covert wireless networks,”  IEEE Trans. Wireless Commun., vol. 17, no. 6, pp. 3656–3669, Jun. 2018.
  * [45] H. M. Wang, Y. Zhang, X. Zhang, and Z. Li,  “Secrecy and covert communications against UAV surveillance via multi-hop networks,”  IEEE Trans. Commun., vol. 68, no. 1, pp. 389–401, Jan. 2020.
  * [46] S. Yan, S. V. Hanly, and I. B. Collings,  “Optimal transmit power and flying location for uav covert wireless communications,”  IEEE J. Sel. Areas Commun., vol. 39, no. 11, pp. 3321–3333, Nov. 2021.
  * [47] O. A. Topal and G. K. Kurt,  “A countermeasure for traffic analysis attacks: Covert communications with digital modulation,”  IEEE Wireless Commun. Lett., vol. 10, no. 2, pp. 441–445, Feb. 2021\.
  * [48] K. Shahzad, X. Zhou, and S. Yan,  “Covert communication in fading channels under channel uncertainty,”  in Proc. IEEE 85th Veh. Technol. Conf. (VTC Spring), pp. 1–5, Jun. 2017.
  * [49] T.M.Cover and J.A.Thomas,  Elements of Information Theory,  Wiley, Jul. 2006.
  * [50] E. Lehmann and J. Romano,  Testing Statistical Hypotheses,  Springer, Jan. 2005.