{'fraction_non_alphanumeric': 0.05054287525756855, 'fraction_numerical': 0.013175622127119988, 'mean_word_length': 4.231446932006634, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Jamming Aided Covert Communication with Multiple Receivers

We consider that a transmitter covertly communicates with multiple receivers under the help of a friendly jammer. The messages intended for different receivers are transmitted in mutually orthogonal frequency bands. An adversary observes all these frequency bands aiming at detecting whether or not communication occurs, while the friendly jammer broadcasts jamming signals to degrade the detection performance of the adversary. We consider a block Rayleigh fading channel model and evaluate the performance of covert communication in two situations: 1) the wireless channels vary slowly such that the transmission ends within one channel coherent time block, and 2) the wireless channels vary fast such that the wireless channels have changed several times before the whole transmission is finished. In the former case, subject to a covertness constraint, we maximize the sum of the effective rates by optimizing the transmit power allocation and the transmission rate for each receiver. In the latter case, we take the channel training process into consideration, and subject to a covertness constraint, we maximize the sum of the ergodic rates by optimizing the power allocation and the pilot length. Though both of the two optimization problems are non-convex, we presented methods to find their global optimal solutions. Besides, we also present methods to find sub-optimal solutions with lower computational complexities. Numerical results are presented to evaluate the performance under the two situations.


I. INTRODUCTION
Providing endogenous network security represents one of the new paradigm shifts of next generation wireless networks [1], [2]. Recently, covert communication has gained considerable attention for its ability to hide the occurrence of the communication itself [3], [4]. Consider a wireless system wherein a transmitter sends its message to an intended receiver. An adversary listens to the wireless channels aiming at detecting whether or not transmission occurs, posing a threat on system security. The technique of covert communication enables the transmitter to reliably communicate with the receiver while ensuring a high probability that the detector of the adversary produces an incorrect result, which greatly enhances wireless security.
The information-theoretic performance limits of covert communication have been studied in [5]- [8]. It was revealed in [5] that covert communication over additive white Gaussian noise (AWGN) channels is subject to the so called square root law (SRL). Specifically, a transmitter is able to reliably and covertly transmit at most O( √ n) bits to a receiver in n channel uses as n → ∞. [6]- [8] further extended the SRL to binary symmetric channels and general discrete memoryless channels. The SRL indicates that the covert communication rate is asymptotically zero, i.e., lim n→∞ √ n n = 0. Some recently efforts have been devoted to improving the covert communication performance. It has been revealed that by exploiting the adversary's uncertainty on the statistical information of its channel outputs, the covert communication performance can be greatly improved, and in some cases, non-vanishing covert communication rates exist, for example, when the adversary has uncertainty on its noise power [9], [10] or when the adversary is uncertain about the transmission time [11]- [13].

A. Wireless covert communication
Covert communication over practical wireless channels has also been extensively investigated. Due to the openness of wireless media and the randomness of wireless environment, wireless transmissions inevitably suffer from co-channel interference. Though co-channel interference is harmful to normal wireless communication process, it also leads to a poor detection performance at the adversary, and thus its impacts have been studied in many existing works [14]- [27].
Covert communication under co-channel interference has been studied in [14]- [19]. In [14], [15], the locations of the interferers were modeled as a Possion point process. Subject to a covert outage probability upper bound, the covert throughput was maximized by optimizing the transmit power and the transmission rate. Covert communication in device-to-device (D2D) underlaying cellular networks was studied in [16], where the wireless signals of the cellular users were treated as interference and used to hide the communications between D2D pairs. [17], [18] studied covert communication in oneway relay networks, and the authors maximized the effective covert rate achieved by the relay. In [19], the dynamicity of the interference environment was considered, and the authors studied the covert throughput scaling law with respect to the codeword length and the variation rate of the background interference environment.
The performance of covert communication can be improved by letting a friendly jammer deliberately broadcast jamming signal to degrade the detection performance of the adversary [20]- [24]. Remarkably, [20] theoretically showed the existence of constant covert communication rates provided that the adversary does not know the jamming power or the instantaneous realization of the jamming channel. In [21], a truncated channel inversion power adaption scheme was proved to be optimal in term of minimizing the outage probability subject to a covertness constraint. In [22], the covert throughput in a random network was studied, wherein the jammer who is closest to the adversary broadcasts jamming signal. In [23] and [24], a friendly jammer was assumed to have multiple antennas and use beamforming to maximize its ability to degrade the detection performance of the adversary.
Instead of relying on external friendly jammers, the works in [25]- [27] assumed that the receiver operates in full-duplex mode, i.e., simultaneously receiving the signal from the transmitter and broadcasting jamming signal to the deteriorate detection performance of the adversary. Due to the imperfect self interference cancellation, the trade-offs between the covertness and the reliability should be carefully designed. In particular, [25] and [26] designed the transmit power of the transmitter and the jamming power of the receiver to maximize the detection error probability at the adversary subject to a lower bound on the effective throughput, and to minimize the outage probability subject to a lower bound on the detection error probability at the adversary, respectively. In [27], an onoff transmission strategy were proposed to optimize the covert communication performance.

B. Motivations, challenges, and contributions
Existing works have presented important insights on the achievable performance of covert communication over wireless channels, however, most of them have only discussed the transmissions from a single transmitter to a single receiver. Motivated by this observation, in this paper, we study a covert communication scenario wherein a transmitter simultaneously communicates with multiple receivers over wireless fading channels. Under the help of a friendly jammer, the transmitter communicates with different receivers in different and mutually orthogonal frequency bands (namely the system works in a frequency-division multiplexing manner). The adversary is able to observe all the frequency bands to make a decision on whether communication occurs or not, whereas we let the friendly jammer broadcast jamming signals in all these frequency bands to degrade the detection performance of the adversary.
It is worth noting that compared to the single-receiver case investigated in literature, analyzing the covertness of the communication becomes more challenging in the considered multiple-receiver case. Specifically, in existing works such as [14]- [18], [21]- [25], [27], the optimal detector of the adversary was shown to be an energy-based detector [20], i.e., comparing the received energy to a predesigned threshold. By exploiting the simple mathematical form of the energy-based detector, the optimal detection performance of the adversary can be accurately characterized by analyzing the false alarm and missed detection probabilities of the optimal energy-based detector. However, in our case, due to the fact that the fading channel coefficients in different frequency bands are different, the energy-based detector becomes strictly suboptimal, which prevents us from analyzing the covertness of the communication by using the methods adopted in existing works.
In addition, the transmit power allocation problem is also a key issue in the considered multiple-receiver scenario. Since the adversary is able to observe all the frequency bands, increasing the transmit power in each single frequency band will lead to an improved detection performance at the side of the adversary. Then, a natural question is that subject to a certain constraint on the covertness, how to allocate the transmit power to the multiple receivers in different frequency bands in order to achieve the optimal communication performance. Based on this observation, the issue of power allocation at the side of the transmitter constitutes the main problem that will be studied in this paper.
The contributions of this paper are summarized as follows, 1) We consider that a transmitter covertly communicates with multiple receivers in mutually orthogonal frequency bands under the help of a jammer. Depending on the temporal dynamic properties of the wireless channels, we consider two different situations: 1) the wireless channels vary slowly over time and a single transmission of the transmitter terminates before the wireless channels have changed, and 2) the wireless channels vary fast and during a single transmission period, the wireless channels change several times. For convenience, the wireless channels are said to be quasi-static and fast-varying in these two situations, respectively. For both situations, analytically tractable upper bounds on the total variation distance between f 1 and f 0 are derived to characterize the covertness of the communication, where f 1 and f 0 denote the probability density functions (PDFs) of the adversary's channel outputs given that covert communication occurs and does not occur, respectively. 2) Under the condition that the channels are quasi-static, we maximize the sum of the effective communication rates of the multiple receivers subject to a constraint on covertness by jointly optimizing the power allocation and the transmission rate for each receiver. The established optimization problem is non-convex, however, by exploiting its monotonic properties, we can obtain the global optimum by using the polyblock outer approximation (POA) method. Besides, we also present a computationally efficient method based on successive convex approximation (SCA) to search for a sub-optimal solution. 3) Under the condition that the channels are fast-varying, due to the limited channel coherent time, we take the channel training process into consideration. Accordingly, we maximize the sum of the ergodic rates subject to a constraint on covertness by jointly optimizing the pilot length and the power allocation. An exhaustive search (ES) based method is proposed to calculate the global optimal solution. In addition, we also present an alternating optimization based method to obtain a sub-optimal solution, which has a lower computational complexity than the ES based method but causes little performance loss.
The rest of this paper is organized as follows: in Section II, we introduce the system model; in Section III and Section IV, we study the covert communication performance under the conditions that the channels are quasi-static and fast-varying, respectively; Numerical results are presented in Section V; and finally, Section VI concludes the paper.
Notations: C, R, and R + denote the set of complex, real, and non-negative numbers, respectively. (·) T and (·) H denote transpose and conjugate transpose, respectively. E(·) and P(·) denote mathematical expectation and probability, respectively. | · | and · denote the absolute value and the norm, respectively. CN (·, ·) and E(·) denote the complex Gaussian and exponential distributions, respectively. Diagonal matrix is denoted by diag(·). For any two real vectors x, y ∈ R n , x ≥ y if x i ≥ y i for ∀i = 1, 2, · · · , n where x i and y i are the i-th elements of x and y, respectively. I m denotes the m-by-m identity matrix.

II. SYSTEM MODEL
Consider that a transmitter communicates with K receivers in K mutually orthogonal frequency bands (K ≥ 1). An adversary aims to detect the occurrence of the communication by observing the K frequency bands. The goal of the transmitter is to reliably communicates with the K receivers while ensuring that the adversary is unable to effectively detect the existence of the communication. Following the works in [20]- [24], we consider that a friendly jammer broadcasts jamming signals in order to deteriorate the detection performance of the adversary, which potentially helps the transmitter to accomplish its transmission covertly. A comprehensive system model is depicted in Fig. 1.

A. Signal model and basic assumption
We assume that the transmitter has M (M ≥ 1) antennas and that the K receivers, the adversary, and the jammer each have a single antenna. For convenience, assume that the transmitter uses the k-th orthogonal frequency band to communicate with the k-th receiver where 1 ≤ k ≤ K. We consider a block fading channel model. Specifically, the channel coefficients remain unchanged in each channel coherent time block and are identically and independently distributed (i.i.d.) in different time blocks. Besides, the channel coefficients in different frequency bands are also i.i.d. For simplicity, we assume that all the channels share the same channel coherent time. In the i-th (i ≥ 1) time block, denote by h k,i ∼ CN (0, I M ) and φ k,i ∼ CN (0, 1) the smallscale fading channels from the k-th receiver to the transmitter and the jammer, respectively. Denote by g k,i ∼ CN (0, I M ) and ϕ k,i ∼ CN (0, 1) the small-scale fading channels from the transmitter and the jammer to the adversary in the k-th frequency band, respectively.
We consider that the communication system works in time division duplex mode. In each channel coherent time block, which consists of N (N > 1) symbol periods, each receiver first sends a pilot sequence of length N t to the transmitter, and the transmitter estimates the channels by exploiting the channel reciprocity. After that, the transmitter uses the remaining N d = N − N t symbol periods to transmit its data to the receivers. In the i-th time block, the pilot sequence received by the transmitter in the k-th frequency band is where P R k is the transmit power of the k-th receiver, S R k ,T is the distance-based path loss between the k-th receiver and the transmitter, x (p) k ∈ C Nt×1 is the pilot sequence sent by the kth receiver, and Z k,i ∈ C M×Nt is the AWGN with each of its elements distributed as CN (0, σ 2 k,T ). For notation simplicity, we assume that P R1 = · · · = P RK = P R and σ 2 1,T = · · · = σ 2 K,T = σ 2 T . Based on (1), the minimum mean square error estimation of h k,i , denoted byĥ k,i , and the estimation error, denoted byh k,i , are respectively given bŷ where µ k σ 2 T PRSR k ,T . After obtainingĥ k,i , the transmitter adopts the maximum ratio transmission beamforming scheme to transmits its data to the k-th receiver, and the beamforming vector is given by b k,i =ĥ k,i ||ĥ k,i || . Accordingly, the k-th receiver receives where P k is the transmit power allocated for the k-th receiver, x is the data sequence, Q k is the jamming power that the jammer emits in the k-th frequency band, S R k ,J is the path loss between the k-th receiver and the jammer, v k,i ∼ CN (0, I N d ) is a sequence of the jamming signals transmitted by the jammer, and z k,i ∼ CN (0, σ 2 R k I N d ) is the AWGN. It is worth noting that by (1) and (4), we assume that the jammer keeps silent during the channel training process and jams the K frequency bands during the data transmission process.
Remark 1: In this paper, we aim to keep the adversary unaware of the transmission from the transmitter while the receivers are allowed to broadcast wireless signals. Similar assumptions have also been adopted in existing works such as [25]- [27] wherein full-duplex receivers are exploited to broadcast jamming signals to improve covert communication performance. In our situation, we let the receivers broadcast pilot sequences so that the transmitter can estimate the wireless channels and design beamforming vectors.

B. Adversary model and metric for covertness
From the perspective of the adversary, it does not know whether communication occurs or not, and the signals received by the adversary in the i-th time block can be modeled as where k ∈ {1, 2, · · · , K}, H 1 and H 0 stand for the hypotheses that communication occurs and does not occur, respectively, w k (i) ∈ C N d ×1 is the signal received in the k-th frequency band, S A,T and S A,J are the path losses from the transmitter and the jammer to the adversary, respectively, andz k,i ∼ CN 0, σ 2 A k I M is the AWGN. We assume that the adversary does not know the instantaneous channel realizations, i.e., ϕ k,i and g k,i for 1 ≤ k ≤ K, but knows their statistic distributions. Besides, we assume that the jamming powers, i.e., Q k for 1 ≤ k ≤ K, are fixed constants.
The goal of the adversary is to determine whether or not communication occurs based on its received signals in the K frequency bands. We consider the following two scenarios depending on the temporal dynamic properties of the wireless channels: 1) The channels are quasi-static: The channel coherent time block is sufficiently long, and the transmission from the transmitter terminates before the end of a channel coherent time block. In this case, we assume that the adversary makes a decision to determine whether communication occurs at the end of each time block.
2) The channels are fast-varying: The wireless channels vary fast as compared to the length of a single transmission. Specifically, we consider that a single transmission is comprised of L > 1 consecutive short channel coherent time blocks. In this case, we assume that the adversary collects the signals received in the L consecutive time blocks together to make a decision.
For a given detector of the adversary, denoted by d, let P FA (d) and P MD (d) be the false alarm and the missed detection probabilities, respectively. In this paper, the communication between the transmitter and the multiple receivers is said to be (1 − ǫ)-covert if it satisfies that where ǫ ∈ (0, 1) is a pre-fixed number. Note that the minimization operation in (6) is with respect to the detector adopted by the adversary. This means that the covertness is guaranteed no matter how the adversary performs the detection [5]. For convenience, we refer to (6) as the covertness constraint in the subsequent part of this paper. Remark 2: From the viewpoint of the adversary, it solves a binary hypothesis testing problem which tests H 1 against H 0 defined in (5). By [34, Theorem 13.1.1], the optimal detector that minimizes the sum of the false alarm and the missed detection probabilities is the likelihood ratio test, which makes a decision by comparing the likelihood ratio with a pre-designed threshold. We assume that the adversary knows the exact values ofp k = P k S A,T ,q k = Q k S A,J , and σ 2 A k , 1 ≤ k ≤ K, so that it can perform the optimal detection. It is also worth noting that the optimal detector is not equivalent to the energy-based detector due to the different fading channel coefficients in different frequency bands and time blocks, see e.g., [19], [20].

C. Covert communication problem formulation 1) Covert communication over quasi-static channels:
In this situation, we assume that the pilot sequence is sufficiently long so that the transmitter can accurately estimate the channel coefficients. As a result, according to (4), the signal-to-noise ratio (SNR) of the k-th user is Here, we have omitted the time block index i because in this case, the transmission terminates within a single time block. Suppose that the transmission rate for the k-th receiver is r k > 0. Due the randomness of the interfering channel φ k , an outage event occurs if ln(1 + SNR k ) ≤ r k , and the outage probability is given by We define the effective rate of the k-th receiver as r The goal of the transmitter is to maximize the sum of the effective rates of the K receivers while ensuring that its transmission is at least (1−ǫ)-covert for some pre-fixed ǫ ∈ (0, 1), i.e.

2) Covert communication over fast-varying channels:
In this case, the wireless channel changes several times before the transmission is finished. Without loss of generality, we assume that the transmission starts in the 1-st time block and ends at the L-th time block. Due to the limited channel coherent time, the pilot length, N t , cannot be neglected as compared to the block length N . As a result, the channel estimation error becomes non-negligible. To evaluate the communication performance, we reformulated (4) as follow, where k ∈ {1, 2, · · · , K} and i ∈ {1, 2, · · · , L}. Following [28, Theorem 1], an achievable ergodic rate of the k-th receiver isr Note that in (11), |E{h H k,i b k,i }| 2 and Var(h H k,i b k,i ) are independent of the time block index i due to the fact that the channel coefficients are i.i.d. over different time blocks. Based on (11), we consider to maximize the sum of the ergodic rates by optimizing the power allocation and the pilot length, i.e, In summary, in this section, we have formulated two optimization problems to enhance the covert communication performance. Subject to the communication being at least (1 − ǫ)-covert, we maximize the sum of the effective rates under the condition that the channels are quasi-static, and we maximize the sum of the ergodic rates under the condition that the channels are fast-varying. In the subsequent two sections, we present our methods to solve these two problems, respectively.
Remark 3: In this paper, we do not impose any transmit power constraint at the transmitter. In fact, to avoid being detected by the adversary, the transmit power is usually very small. This means that in many situations, the power consumption does not constitute a performance bottleneck of covert communication. Based on this viewpoint, in this paper, we assume that the transmitter has enough power budget to accomplish its transmission.

CHANNELS
In this section, we study the covert communication performance under the condition the channels are quasi-static. In the following, we first analyze the covertness constraint and transform it into a mathematically tractable form. Then, we present our method to solve the effective sum-rate maximization problem.

A. Covertness constraint analysis
By assumption, the adversary independently runs its detection process in each single time block. In the i-th time block, the detection problem at the adversary can be formulated as the following binary hypothesis testing problem, , where are the PDFs of W (i) under the condition that transmission occurs and does not occur in the i-th time block, respectively. For notation simplicity, we omit the time block index i in the subsequent part of this section.
can be written as where f is the PDF of a circular symmetrical complex Gaussian random vector with covariance matrix xI N d . The PDFs of U k and V k are presented in the following lemma.
Lemma 1: The PDFs of U k and V k are given by Besides, due to the fact that g k ∼ CN (0, I M ) and that b k is independent of g k with ||b k || 2 = 1, it is straight that g H k b k ∼ CN (0, 1) and thus |g H k b k | 2 ∼ E(1). Note that in covert communication, the signal power is usually very small compared to the jamming-plus-noise power in order to guarantee the covertness of the transmission. Therefore, in the following, we assume thatp k <q k holds.
By [34, Theorem 13.1.1], the optimal detector of the adversary that minimize the sum pf the false alarm and missed detection probabilities satisfies min d P FA (d) ). Therefore, the covertness constraint (6) can be re-formulated as V(f are presented in (15a) and (15b), respectively. As N d → ∞, the total variation distance between f Proof: Please refer to Appendix A.
For the special case of K = 1, V (f U , f V ) can be written in a closed form, where which involves a K-fold integral. Therefore, in the following, we derive an analytically tractable upper bound on V (f U , f V ) for the ease of numerically optimizing the covert communication performance.
is due to the triangle inequality, i.e., |u + v| ≤ |u| + |v| for any real-valued u and v, step (b) is obtained by calculating the integral with respect to x 2 , x 3 , · · · , x K for the first term and calculating the integral with respect to x 1 for the second term, and using the definition of V(·, ·), step (c) follows from (17), step (d) is obtained by applying the steps Based on Proposition 2, in order to make sure that the communication is at least (1 − ǫ)-covert, we consider the following constraint, where η(x) x 1 1−x . The tightness of using K k=1 η(χ k ) as an upper bound on V (f U , f V ) will be numerically evaluated in Section V-A.

B. Covert communication performance optimization
In previous subsection, we have presented our method to reformulate the covertness constraint. In this subsection, we maximize the effective sum-rate subject to the covertness constraint (20).
For the case with K = 1, based on (17), we ensure the covertness of the communication by letting χ 1−x = ǫ in (0, 1). Substituting P * 1 into (9), we obtain the following problem, where γ 1 e r1 − 1 is viewed as the optimization variable, The optimal solution to (21), denoted by γ * 1 , is presented below.

C. An SCA based method
In this subsection, we present an alternative method to search for a suboptimal solution to (22). Define t = [t 1 , t 2 , · · · , t K ] T with t k χ k γ k for 1 ≤ k ≤ K. We view t and γ as the optimization variables and recast (22) as the following problem, The main difficulty on solving (24) lies in the fact that (24b) is not a convex constraint and that the objective of (24) is non-concave. Note that η(x) satisfies that η(0) = 0, lim x→0 η ′ (x) = 1. By checking the high order derivatives of η(x), if can be proved that η ′′ (x) ≤ 0 for ∀x ∈ (0, 1). Therefore, η(x) is a monotonically increasing concave function, and satisfies Based on (25), we can replace the inequality constraint (24b) with t T γ ≤ ǫ to obtain a lower bound on optimal covert communication performance. We further simplify (24), by introducing two vector-valued slack variables α = [α 1 , α 2 , · · · , α K ] T ∈ R K + and β = [β 1 , β 2 , · · · , β K ] T ∈ R K + , and obtain the following problem, where R(α, β) (26) is still a non-convex problem due to the non-convex objective function R s (α, β) and the non-convex constraint (26b). However, (26), in its current form, can be efficiently handled by the SCA method. In brief, the SCA method handles a non-convex optimization problem by transforming it into a series of parameterized convex problem. By iteratively solving the obtained convex problems with the parameters at each iteration being the optimal solution obtained in previous iteration, the SCA method generates a series of solutions which converges to a Karush-Kuhn-Tucker (KKT) point of the original non-convex problem. A brief introduction of the SCA method is presented in Appendix C. In our case, in the (j + 1)-th iteration of the SCA method, we need to solve the following convex problem, min t,γ,α,β∈R+ where ρ (j) = α (j) + β (j) , Λ , Γ (j) = diag(γ (j) ), T (j) = diag(t (j) ), and {α (j) , β (j) , t (j) , γ (j) } consists of the optimal solution obtained in the j-th iteration. The derivation of (27) is explained in Appendix C. To start the SCA iteration, {α (0) , β (0) , t (0) , γ (0) } can be initialized to be any feasible solution to (26). Note that as (27) is a convex problem, it can be efficiently solved by using software such as CVX [38].
Remark 4: In this section, we have proposed the POA and the SCA methods to maximize the effective sum-rate subject to a covertness constraint. Though the POA method provisdes us the globally optimal solution, it generally exhibits a high computational complexity. Specifically, after a vertex is searched, K new vertexes will be added into the vertex set (see Appendix B for more details), meaning that the worstcase computational complexity of the POA method scales with K exponentially. As for the SCA method, it solves a sequence of convex optimization problems, i.e., (27), by using standard convex optimization algorithms, for example, the interior-point method combined with Newton's method [37], and the computational complexity grows with K in a polynomial manner. Based on these observations, we conclude that the POA method better suits to the case where the number of receivers, K, is small, whereas if K is large, the SCA is computationally more efficient.

IV. COVERT COMMUNICATION OVER FAST-VARYING CHANNELS
In this section, we consider the situation where the wireless channels are fast-varying. In the following, we first analyze the covertness constraint, based on which we then discuss the ergodic sum-rate maximization problem.

A. Covertness constraint analysis
Without loss of generality, assume that the transmission starts in the 1-st and ends at the L-th time blocks (if it occurs). Accordingly, the detection problem of the adversary can be formulated as the following hypothesis testing problem, where W (i), 1 ≤ i ≤ L, are defined in (13),f ). Note that (28) differs from (13) in the sense that in (28), the signal matrices obtained in different time blocks are collected together to make a decision, while in (13), the decision is made based solely on the signal matrix obtained in a single time block.
In principle, we can use the method in Proposition 1 and 2 to obtain an upper bound on V(f ). However, given that the channel coherent time is short, such an upper bound, which requires N d → ∞, becomes inaccurate. In the following, we use Pinsker's inequality, see e.g., [36, Lemma 11.6.1] to derive a mathematically tractable bound on V(f where step (a) is because W (i) for 1 ≤ i ≤ L are i.i.d. random matrices. Based on (29), we use the following constraint to ensure the covertness of the transmission, The mathematical expression of D f 1,k ) can be simplified as (31) in the top of next page. where in step (a), we express the complex-valued integral variables w 1 , w 2 , · · · , w n in polar coordinates, i.e., w i = x i e jθi with x i ∈ R + and θ i ∈ (−π, π], step (b) is obtained by first calculating the integral w.r.t. θ and then changing the integral variables by letting y i = x 2 i for 1 ≤ i ≤ n; and finally, in step (c), we let z = n j=1 y j and z i = i j=1 y j for 1 ≤ i ≤ n − 1, and calculate the integral w.r.t. (z 1 , z 2 , · · · , z n−1 ). We further have that Similarly, we obtain We complete the proof by inserting (32) and (33) into (31) and letting z ← z Based on Proposition 4, for a given power allocation vector p = [P 1 , P 2 , · · · , P K ] T , we can check its feasibility subject to constraint (30). However, (30) is too complicated to facilitate to solve optimization problem (12). In view of the fact that p k is generally small to ensure the covertness, the following proposition provides a way to approximate D(f We further have that Based on Proposition 5, we can approximate the covertness constraint (30) by, where χ k for 1 ≤ k ≤ K are defined in (22). Remark 5: For q k > 0, it can be verified through numerical calculation that 1,k ) dp 3 k p k =0 < 0, meaning that the lefthand-side (LHS) of (34) is an upper bound on the LHS of (30) when p k is sufficiently small. Note that in (34), with N and N t fixed, the transmit power P k is upper bounded by O( 1 √ L ), which is analogous to the SRL. Note that here, L stands for the number of time blocks but not the number of channel uses. Similar result has also been presented in [19].
In the following, we present our method to solve (37). We observe that if τ is fixed,R(χ, τ ) is a concave function of χ, and thus (37) becomes a convex optimization problem, the global optimum of which can be efficiently obtained. In light of this, we can solve (37) by first optimizing χ with τ fixed and then calculating the optimal τ by an exhaustively search over the finite discrete set Z N . For a given value of τ , the optimal χ is presented in the following proposition.
Proposition 6: Let τ be a fixed constant in (0, 1). For 1 ≤ k ≤ K, the optimal χ k , denoted by χ * k , is the unique solution in (0, ∞) to the following equation, whereG k τ G k ,Ẽ k τ E k +μ k ,F k τ F k,1 + F k,2 , and λ * > 0 is the optimal Lagrange multiplier satisfying that K k=1 ζ k (τ ) 2 (χ * k ) 2 = 2 L ǫ 2 . Proof: Let λ be the Lagrange multiplier associated with constraint (37b). The KKT conditions of (37) state that ∂L(χ,λ) is the Lagrangian function, which results in (38). SinceR(χ, τ ) and (38) is active at the optimum, meaning that the optimal Lagrange multiplier λ * satisfies that Proposition 6,(38) can be solved in a closedform due to that fact that (38) is a cubic equation. As for λ * , it can be efficiently searched by using the bisection method. Based on Proposition 6, we can obtain the optimal χ for each fixed value of τ . Then, the optimal τ can be obtained by an exhaustive search over the set Z N .

C. An alternating optimization based method
In this subsection, we present a method to obtain a suboptimal solution of (37) with a lower computational complexity than the ES based method proposed in previous subsection. The basic idea is that we replace constraint (37b) with the following one In fact, by setting τ = 0 in (37b), we equivalently let the adversary observe N > N d samples in each time block. Therefore, this improves the detection performance of the adversary and leads to a lower bound on the achievable communication performance. Since (39) is independent of τ , we replace (37b) with (39) and solve the resulting problem by using the AO method. Specifically, we alternatingly optimize χ with τ fixed and optimize τ with χ fixed until the value ofR(χ, τ ) converges. For the sub-problem of optimizing χ, the optimal solution is presented in Proposition 6. Given that χ is fixed, the subproblem of optimizing τ is given by, To solve (40), we relax τ to be a real number in (0, 1), and the optimal real-valued τ is presented in the following proposition.
Remark 6: In this section, the ES and the AO methods are proposed to maximize the ergodic sum-rate subject to a covertness constraint. The ES method treats τ = N d N as a discrete variable and searches the optimal τ in Z N exhaustively. The AO method relaxes τ as a continuous variable in (0, 1) and optimizes τ and χ alternately. Theoretically, the AO method only produces a sub-optimal solution. Note that for each fixed value of τ , both the ES and the AO methods need to compute the optimal χ, denoted by χ * (τ ), by using the bisection method to solve for the optimal dual variable according to Proposition 6. Conceptually, the ES method needs to compute χ * (τ ) for N − 1 times as |Z N | = N − 1. For the AO method, there is generally no guarantee on the maximal number of iterations before convergence. However, our simulation results in Section V-B show that only a few iterations are sufficient for the AO method to converge, and thus the AO method is computationally more efficient than the ES method.

V. NUMERICAL RESULT
We numerically evaluate the covert communication performance in this section. We let the transmitter, the adversary, the jammer, and the k-th receiver (1 ≤ k ≤ K) be in the same two-dimensional plane with their positions given by s T , s A , s J , and s R k , respectively. The distance-based path loss between node a and node b is S a,b = ||s a − s b || −4 where a, b ∈ {T, A, J, R 1 , · · · , R K }. We set s T = (0, 0), s A = (−d A , 0), and s J = (−d J , 0). The positions of each receivers, s R k for 1 ≤ k ≤ K, are independently and uniformly distributed in a circular region with the center and the radius being (d R , 0) and r c , respectively. Unless specified, we set d A = 150 m, d J = 250 m, d R = 150 m, and r c = 30 m, σ 2 A1 = · · · = σ 2 AK = −80 dBm, σ 2 R1 = · · · = σ 2 RK = −80 dBm, and σ 2 T = −80 dBm, P R1 = · · · = P RK = P R = 5 dBm, and Q 1 = · · · = Q K = Q = 25 dBm.

A. Covert communication over quasi-static channels
In this subsection, we evaluate the covert communication performance under the condition that the channels are quasistatic. Unless specified, we set M = 20 and ǫ = 0.005.
In Fig. 2, we compare V(f U , f V ) with its upper bound derived in (18). Here, we set K = 2 and χ 1 = χ 2 = χ for illustrative simplicity. As there is no closed-form expression for V(f U , f V ), we calculate it via numerical integration. Fig. 2 reveals that the upper bound derived in (18) tightly approximates V(f U , f V ), especially when χ → 0. In Fig.  2, we also illustrate two widely used upper bounds on total variation distance, namely, the Hellinger distance based upper bound, see e.g., [34,Theorem 13.1.2], and the KL-divergence based upper bound obtained by using Pinsker's inequality [36, Lemma 11.6.1]. However, Fig. 2 shows that in the considered scenario, these bounds are not as tight as the proposed upper bound in (18). Fig. 3 illustrates the convergence behavior of the SCA method. We run the SCA method under four groups of randomly generated system parameters. Specifically, for each curve in Fig. 3, the locations of the receivers are randomly generated as introduced at the beginning of this section, and the channel coefficients are randomly generated by using the complex Gaussian distribution as assumed in Section II-A. From Fig. 3, we can see that the SCA method converges within a few iterations. In fact, through a large number of numerical experiments, we observe that under the considered system settings, the SCA method converges within 10 iterations with the relative error smaller than 0.01 in most cases.
In Fig. 4, we plot the optimized effective sum-rate, R (e) , versus the jamming power, Q, wherein the results obtained by using the SCA method are compared to those obtained by using the POA method. Fig. 4 indicates that the performance obtained by the SCA method is nearly optimal. Besides, Fig.  4 also shows that the covert communication performance gets improved as the jamming power increases. This is because  whether or not the covertness constraint (20) is satisfied only depends on the values of χ 1 , χ 2 , · · · , χ K . By definition, χ k is proportional to P k Q k , meaning that we can proportionally increase P k with Q k without sacrificing the covertness. And thus if P 1 , P 2 , · · · , P K are chosen such that P1 Q1 , P2 Q2 , · · · , PK QK are fixed constants and satisfy the covertness constraint, then it can be checked that the effective rate r (e) k increases with Q k for 1 ≤ k ≤ K. Therefore, the optimal effective sum-rate also increases with the jamming power.
In Fig. 5, we plot the maximal effective sum-rate as a function of M . Fig. 5 inspires us that increasing M is a promising method to enhance the covert communication performance. In fact, in the considered scenario, with the knowledge about the channel state information (CSI), the transmitter can design beamforming vectors to enhance the SNRs at the receivers. Note that due to the mismatch between the beamforming vector of the transmitter and the CSI of the adversary, increasing M does not increase the signal power received by the adversary, meaning that the covert communication performance can be improved by using more antennas without sacrificing the covertness.

B. Covert communication over fast-varying channels
In this subsection, we evaluate the ergodic sum-rate under the condition that the channels are fast-varying. Unless spec-  We first study the convergence behavior of the AO method introduced in Section IV-C. For each curve in Fig. 6, we randomly generate the locations of the receivers and then apply the AO method to optimize the ergodic sum-rate. We set the initial value of τ as 0.5. The value ofR obtained in each iteration of the AO method are plotted. As we can see from Fig. 6, the AO method converges very fast. In fact, through a large number of numerical experiments, we observe that the improvement onR becomes negligible after two iterations of the AO method in most cases.
In Fig. 7, we showR as a function of P R . Both the AO and the ES methods are implemented to solve (37). Note that the solutions obtained by the ES method are global optimal. From Fig. 7, it can be seen that the performance achieved by the AO method is nearly optimal. Combining Fig. 6 and 7, we conclude that the AO method is computationally efficient while causes little performance loss. Fig. 7 also reveals that the communication performance gets significantly improved as the P R increases, which is due to that fact that with a larger value of P R , the transmitter can estimate the CSIs more accurately.
The influences of the jamming power on the ergodic sumrate are plotted in Fig. 8. From Fig. 8, we observe that increasing the jamming power does not always improve the communi- cation performance, and there seems to exist an optimal value of the jamming power. In fact, through numerical calculation, we observe that ζ(q k , N −N t ) in constraint (34) increases with q k (recall that q k ∝ Q k ). This means that subject to (34), the increase of Q k inevitably leads to the decrease of χ k , which finally reducesR. We notice that the effects of increasing the jamming power are quiet different in Fig. 4 and 8. Specifically, in Fig. 4, the covert communication performance improves as the jamming power increases, whereas in Fig. 8, the covert communication performance becomes degraded when the jammer power is sufficiently large. This inspires us that in practice, the jamming power should be carefully designed according to the temporal dynamic properties of the wireless channel of the wireless channels. In Fig. 9, we illustrate the ergodic sum-rate,R, as a function of the lower bound on the error probability of the adversary, i.e., 1 − ǫ, wherein we fix N × L = 1.5 × 10 3 . As 1 − ǫ approaches 1,R monotonically decreases, which is the direct result of the reduced transmit power in order to satisfy the covertness constraint. It is also worth noting that R decreases with L. In fact, as L increases, the adversary obtains more signal samples that are under different random channel realizations. This helps the adversary to average out the unknown fading channel coefficients, and thus leads to a more reliable detection result.

VI. CONCLUSION
Covert communication between a transmitter and multiple receivers under the help of a friendly jammer has been studied. When the channels are quasi-static, we maximized the sum of the effective rates. When the channels are fast-varying, we maximized the sum of the ergodic rates. The optimal solutions to both of the two optimization problems can be numerically obtained. And we have also presented methods to search for sub-optimal solutions with relatively lower computational complexities. Numerical results have been presented to show the covert communication performance. We revealed that increasing the number of antenna at the transmitter is promising in improving the covert communication performance. Besides, it has also been revealed that the existence of a friendly jamming greatly helps improving the covert communication performance.

A. The proof of Proposition 1
Let n ≥ 1 be an integer, by the definition of the total variation distance, we obtain that for some sufficiently large positive constant G, step (a) is obtained by using step (a), (b), and (c) in the derivation of (31), and in step (b), we make a change of variable y k ← z k /n for 1 ≤ k ≤ K.
almost everywhere. Besides, for n ≥ 2, there exists some constraint c > 0 (independent of n) such that for ∀x > 0, f Lemma 2 indicates that for ∀n ≥ 1, V k (y k ) is bounded above by some constant. Therefore, by using the bounded convergence theorem, see e.g., [33,Theorem 1.4], we have that 2 , Proposition 1 is proved by using Lemma 3 and letting G → ∞.
Proof of Lemma 2: For simplicity, we only prove for the case of f (n) V k (x), and the proof can be directly extended to the case of f (n) where step (a) is because n! ≥ √ 2πnn n /e −n (Stirling's formula) and ̺(z) 1 − z + ln z, and step (b) follows from the fact that ̺(z) ≤ ̺(1) = 0 for any z > 0, and that ̺(z) is increasing in z ∈ (0, 1).
where x ′ is some real number that lies between (1 − ξ)x and (1 + ξ)x, step (a) is due to the fact that . Besides, by using the steps in (46) and (47), it can be shown that for ∀x ≥ 0, f q k and thus is bounded above by some constant that is independent of n.
Proof of Lemma 3: By changing the order of integration, X k,1 satisfies X k,1 = (1), which completes the proof.

B. A brief introduction on the POA method
In this part, we briefly introduce the POA method. We first present some basic definitions.
Definition 1 (Monotonic functions): A function h : R n → R is monotonically increasing if h(x) ≥ h(y) when x ≥ y.
Definition 2 (Hyper-rectangles (a.k.a. Boxes)): Let a, b ∈ R n with x ≤ y, then the set of x ∈ R n such that a ≤ x ≤ b is a box in R n , which is denoted by [a, b].

Definition 3 (Normal sets):
A set X is normal if ∀x ∈ X , then the box [0, x] ⊆ X .
Definition 4 (Polyblocks): A set P ∈ R n + is a polyblock if P = ∪ v∈V [0, v], where V is refer to as the vertex set of P with |V| < +∞.
Definition 5 (Monotonic optimization): Let h : R n → R be a monotonically increasing function, and G ⊆ [0, v] for some v ∈ R n + is a compact normal set with nonempty interior. Then, max x∈G h(x) is a monotonic optimization problem.
Based on the definitions above, let's consider a simple monotonic optimization problem in the form of max x∈G h(x), where h : R n + → R is continuous and monotonically increasing, and G ⊆ R n + is a compact normal set. The basic idea behind the POA method to solve such a monotonic optimization problem is to iteratively generate a sequence of polyblocks to approximate the feasible set G, i.e., P 1 ⊇ P 2 ⊇ · · · ⊇ P j ⊇ · · · ⊇ G, where P j is the polyblock used to approximate G in the j-th iteration. Denote by V j the vertex set of P j . The j-th iteration of the POA method consists of the following steps: 1) Find the best vertex: x Fig. 10: Illustration of the POA method in a two-dimensional space, where the shadow area is the feasible set G. Here, in the first iteration, G is approximated by P1 whose vertex set is V1 = {v1,1}, see Fig. 10(a). After the first iteration, we obtain an refined polyblock outer approximation, i.e., P2, whose vertex set if V2 = {v2,1, v2,2}, see Fig. 10(b). Following similar steps, we obtain P3, whose vertex set is V3 = {v3,1, v3,2, v3,3}, see Fig. 10(c). The iteration continues until some pre-defined termination criterion is satisfied.
2) Project x and e i is the i-th standard basis of R n . The iteration continues until some certain termination criterion is satisfied, which we introduce later. We illustrate the POA iteration introduced above in Fig. 10 for better understanding. As we can see that the POA method iteratively generates a sequence of feasible points. Denote the set of the feasible points obtained after the j-th iteration as j }, and among X j , denote the best feasible point as x * j = max x∈Xj h(x). Then, at the (j + 1)-th iteration, it is reasonable to terminate the iteration process if h(x (v) j+1 ) − h(x * j ) ≤ δ where δ > 0 a given tolerance. In fact, this means that h(x * j ) ≤ max x∈G h(x) ≤ max x∈Pj+1 h(x) = max x∈Vj+1 h(x) = h(x (v) j+1 ) ≤ h(x * j ) + δ, and thus x * j is a δ-optimal solution.
Note that in general, the number of points in the vertex set V j increases exponentially with j, leading to a high computational complexity. There are some methods to reduce the points in V j , for example, [31,Proposition 2.7] and [31,Section 3.3], which is helpful to reduce the computational burden. We also note that in this part, we have only considered a simplified version of the canonical monotonic optimization formulation, which is, however, enough to handle the optimization problem in (23). For more details about the framework of monotonic optimization and the POA method, please refer to [30], [31] and references therein.

C. A brief introduction on the SCA method and the derivation of (27)
Consider a general non-convex optimization problem min x∈X g 0 (x) + h 0 (x), where X {x : g t (x) + h t (x) ≤ 0, 1 ≤ t ≤ T } is the feasible set with T being the number of the constraints, g t (x) for 0 ≤ t ≤ T are convex functions, and h t (x) for 0 ≤ t ≤ T are non-convex functions. The SCA method handles such a general non-convex optimization problem by replacing the non-convex function h t (x) with a convex approximation near some feasible pointx, denoted bŷ h t (x;x), and iteratively solving the resultant convex problem. Specifically, in the j-th iteration, the SCA method solves min x∈Xj g 0 (x) + h 0 (x; x j−1 ) where X j {x : g t (x) + h t (x; x j−1 ) ≤ 0, 1 ≤ t ≤ T }, where x j−1 is set to be the optimal solution obtained in the (j − 1)-th iteration. Consequently, the SCA method generates a sequence of solutions {x 1 , x 2 , · · · }. It has been shown in literature that for ∂x x=x , and 4)ĥ t (x;x) is continuous in (x;x), then under some mild assumptions, the limiting point generated by the SCA method is a KKT solution to the original non-convex problem. For more details about the SCA method, please see, e.g. [32], [33] and references therein. In our problem in (26), the objective and the constraint in (26b) are non-convex. To use the SCA method, in the (j+1)-th iteration, we approximate the non-convex part of the objective, i.e., −||α + β|| 2 , by using its first order Taylor expansion. For the non-convex constraint (26b), we approximate it using are presented below (27). It can be verified that the convex approximations adopted here satisfy the conditions mentioned above, and thus enable the SCA method.