Optimal Downlink Transmission for Cell Free SWIPT Massive MIMO Systems with Active Eavesdropping

TL;DR

This paper investigates secure SWIPT in cell-free massive MIMO systems with active eavesdroppers, deriving energy and secrecy rate bounds, and proposing a power control method that enhances security and energy harvesting performance.

Contribution

It introduces a novel power control approach using relaxed SDP for secure downlink transmission in cell-free MIMO with active eavesdroppers, considering nonlinear power constraints.

Findings

Cell-free MIMO outperforms colocated MIMO at low AHE constraints.

The proposed SDP-based power control achieves rank-one globally optimal solutions.

Cell-free MIMO is more resilient to active eavesdropping power increases.

Abstract

This paper considers secure simultaneous wireless information and power transfer (SWIPT) in cell-free massive multiple-input multiple-output (MIMO) systems. The system consists of a large number of randomly (Poisson-distributed) located access points (APs) serving multiple information users (IUs) and an information-untrusted dual-antenna active energy harvester (EH). The active EH uses one antenna to legitimately harvest energy and the other antenna to eavesdrop information. The APs are networked by a centralized infinite backhaul which allows the APs to synchronize and cooperate via a central processing unit (CPU). Closed-form expressions for the average harvested energy (AHE) and a tight lower bound on the ergodic secrecy rate (ESR) are derived. The obtained lower bound on the ESR takes into account the IUs' knowledge attained by downlink effective precoded-channel training. Since the transmit power constraint is per AP, the ESR is nonlinear in terms of the transmit power elements of the APs and that imposes new challenges in formulating a convex power control problem for the downlink transmission. To deal with these nonlinearities, a new method of balancing the transmit power among the APs via relaxed semidefinite programming (SDP) which is proved to be rank-one globally optimal is derived. A fair comparison between the proposed cell-free and the colocated massive MIMO systems shows that the cell-free MIMO outperforms the colocated MIMO over the interval in which the AHE constraint is low and vice versa. Also, the cell-free MIMO is found to be more immune to the increase in the active eavesdropping power than the colocated MIMO.

Tables

Table 1. TABLE I: List of Notations

Notation	Description
$𝒂$, $𝑨$	vectors and matrices are denoted by boldface lowercase and boldface uppercase letters, respectively
${𝑰}_N$	denotes the $N\times N$ identity matrix
$\text{diag}\left(𝒔\right)$	a matrix whose diagonal entries are the entries of vector $𝒔$ and zeros elsewhere
$\text{diag}\left(𝑺\right)$	a column vector whose entries are the diagonal entries of matrix $𝑺$
$𝑺⪰0$	indicates that $𝑺$ is a positive semidefinite matrix
${(\cdot )}^T\text{and}$ ${(\cdot )}^H$	the transpose and the conjugate transpose, respectively
$\text{tr}(\cdot )\text{and}$ ${\text{log}}_2(\cdot )$	the trace of a matrix and logarithm to base 2, respectively
$|\cdot |\text{and}$ $\parallel \cdot \parallel$	the absolute value of scalars and the Euclidean norm, respectively
$ℛ,{ℛ}_{+}^n$, ${𝒮}_{+}^n$ and ${𝒞}^{m\times n}$	sets of real numbers, nonnegative real numbers, symmetric positive semidefinite $n\times n$ real matrices and complex $m\times n$ matrices, respectively
$𝒞𝒩(\mathrm{𝟎},𝚺)$	circularly symmetric complex Gaussian distribution of a random vector with zero mean and covariance matrix $𝚺$
$\text{cov}(x,y)$ $\text{and}\text{var}\left(x\right)$	the covariance between the random variables (RVs) $x$ and $y$, and the variance of $x$, respectively
$\left\{{𝒂}_n\right\}\text{and}$ ${\left\{a_{m,n}\right\}}_m$	a set of all vectors indexed by $n$ and a set of all scalars indexed by $m$, respectively
${\left[𝒂\right]}_n\text{and}$ ${\left[𝑨\right]}_{n,m}$	the $n$th entry of $𝒂$ and the $(n,m)$th entry of $𝑨$, respectively
$𝑩=\text{null}\left(𝑨\right)$	means $𝑨𝑩=\mathrm{𝟎}$ and $𝑩{𝑩}^H=𝑰$
${\left[x\right]}^{+}$	is equivalent to $\max (x,0)$

Table 2. TABLE II: Selected Values of System Parameters

Parameter	Value
$A_a$, ${\lambda }_a$	$1\times 1$ ${\text{Km}}^2$, $1.5\times {10}^{-4}{m}^{-1}$
$A_u$, $M$	$300\times 300$ ${\text{m}}^2$, 3
$\alpha$, $\sigma$	2.5, $8\text{dB}$ [37]
$P$, $P_E$, $P_t$	$1W$, $1W$, $500mW$
$\tau$, ${\tau }_d$, $\zeta$	10, 10, 0.5 [38]

Table 3. TABLE III: Relative values of Δ 1 subscript Δ 1 \Delta_{1} and Δ ¯ 2 subscript ¯ Δ 2 \overline{\Delta}_{2}

Realization	1st	2nd	3rd	4th
${\Delta }_1$	$9.8\times {10}^{-3}$	$7.3\times {10}^{-3}$	$1.6\times {10}^{-2}$	$4.6\times {10}^{-3}$
${\overline{\Delta }}_2$	$6.9\times {10}^{-5}$	$1.1\times {10}^{-4}$	$3.8\times {10}^{-5}$	$1.2\times {10}^{-4}$
$\frac{{\Delta }_1}{{\overline{\Delta }}_2}$	$1.4\times {10}^2$	$0.65\times {10}^2$	$4.1\times {10}^2$	$0.3\times {10}^2$

Equations

$$\mathbf{Y}=$$

$$\sum _{i=1}^M\sqrt{P_I}{\mathbf{h}}_i{\mathbf{\psi }}_i^T+\sqrt{\varphi P_E}{\mathbf{g}}_E{\mathbf{\psi }}_k^T+\sqrt{(1-\varphi )P_E}\mathbf{g}{\mathbf{\psi }}_E^T+\mathbf{N},$$

$${\hat{\mathbf{h}}}_i={\mathbf{C}}_i{\mathbf{y}}_i,{\mathbf{C}}_i=\sqrt{P_I}{\mathbf{\Gamma }}_i{\left(\tau P_I{\mathbf{\Gamma }}_i+{\delta }_{ik}\tau \varphi P_E\mathbf{\Gamma }+{\sigma }_n^2{\mathbf{I}}_N\right)}^{-1},$$

$${\mathbf{y}}_i=\mathbf{Y}{\mathbf{\psi }}_i^{\ast }=\tau \sqrt{P_I}{\mathbf{h}}_i+{\delta }_{ik}\tau \sqrt{\varphi P_E}{\mathbf{g}}_E+\mathbf{N}{\mathbf{\psi }}_i^{\ast },$$

$$\hat{\mathbf{g}}=\mathbf{C}\mathbf{y},\mathbf{C}=\sqrt{(1-\varphi )P_E}\mathbf{\Gamma }{\left(\tau (1-\varphi )P_E\mathbf{\Gamma }+{\sigma }_n^2{\mathbf{I}}_N\right)}^{-1},$$

$$\mathbf{y}=\mathbf{Y}{\mathbf{\psi }}_E^{\ast }=\tau \sqrt{(1-\varphi )P_E}\mathbf{g}+\mathbf{N}{\mathbf{\psi }}_E^{\ast },$$

$$\frac{{\mathbf{y}}_i^H{\mathbf{y}}_i-{\tau }^2{P}_I\text{tr}\left({\mathbf{\Gamma }}_i\right)-N\tau {\sigma }_n^2}{{\tau }^2\text{tr}\left(\mathbf{\Gamma }\right)}\stackrel{N\to \infty }{\to }{\delta }_{ik}\varphi P_E,$$

$${\mathbf{x}}_k=\sum _{i=1}^M{\mathbf{w}}_i{q}_i+{\bar{\mathbf{w}}}_{k}z+\mathbf{w},$$

$${\mathbf{w}}_i=\text{diag}\left({\mathbf{p}}_i\right){\hat{\mathbf{h}}}_i^{\ast },{\mathbf{p}}_i={\left[\sqrt{p_{i,1}},\dots ,\sqrt{p_{i,N}}\right]}^T,$$

$${\bar{\mathbf{w}}}_k=\text{diag}\left(\bar{\mathbf{p}}\right){\hat{\mathbf{h}}}_k^{\ast },\bar{\mathbf{p}}={\left[\sqrt{{\bar{p}}_1},\dots ,\sqrt{{\bar{p}}_N}\right]}^T,$$

$$\mathbf{w}=\text{diag}\left(\mathbf{p}\right){\hat{\mathbf{g}}}^{\ast },\mathbf{p}={\left[\sqrt{p_1},\dots ,\sqrt{p_N}\right]}^T.$$

$$y_{k,i}={\mathbf{h}}_i^T{\mathbf{x}}_k+n_i,$$

$$y_k={\mathbf{g}}^T{\mathbf{x}}_k+\hat{n},$$

$$y_{E_k}={\mathbf{g}}_E^T{\mathbf{x}}_k+\bar{n}.$$

$${\mathbf{y}}_{I_i}={\mathbf{h}}_i^T{\mathbf{X}}_d+{\mathbf{n}}_i=\sum _{j=1}^M{a}_{i,j}{\mathbf{\psi }}_{d_j}^T+{\mathbf{n}}_i,$$

$$\frac{{\mathbf{p}}_i^T{\mathbf{\Gamma }}_i{\mathbf{R}}_i{\mathbf{p}}_i}{{\mathbf{p}}_i^T{\mathbf{\Gamma }}_i{\mathbf{R}}_i{\mathbf{p}}_i+{\tau }_d{\sigma }_n^2}{y}_{I_i},$$

$${\hat{a}}_{i,i}=\frac{y_{I_i}}{{\tau }_d}=a_{i,i}+{\tilde{a}}_{i,i},$$

$$y_{k,i}$$

$$=\mathbb{E}\left[a_{i,i}|{\hat{a}}_{i,i}\right]q_i+\left(a_{i,i}-\mathbb{E}\left[a_{i,i}|{\hat{a}}_{i,i}\right]\right)q_i+Z_{k,i},$$

$$Z_{k,i}=\sum _{j\ne i}{a}_{i,j}{q}_j+{\mathbf{h}}_i^T\left({\bar{\mathbf{w}}}_{k}z+\mathbf{w}\right)+n_i.$$

$$\mathbb{E}\left[a_{i,i}|{\hat{a}}_{i,i}\right]$$

$$=\mathbb{E}\left[{\hat{a}}_{i,i}|{\hat{a}}_{i,i}\right]+\mathbb{E}\left[{\tilde{a}}_{i,i}|{\hat{a}}_{i,i}\right]={\hat{a}}_{i,i}+\mathbb{E}\left[{\tilde{a}}_{i,i}\right]={\hat{a}}_{i,i},$$

$$a_{i,i}-\mathbb{E}\left[a_{i,i}|{\hat{a}}_{i,i}\right]={\tilde{a}}_{i,i},$$

$$R_k=\mathbb{E}\left\{{\text{log}}_2\left(1+{\text{SINR}}_k\right)\right\},$$

$${\text{SINR}}_k$$

$$=\frac{{\left|{\hat{a}}_{k,k}\right|}^2}{\text{var}\left({\tilde{a}}_{k,k}\right)+\text{var}\left(Z_{k,k}\right)},$$

$${\underline{\text{SINR}}}_k=\frac{{\tau }^2{P}_I{c}_k^2}{\sum _{j\ne k}{c}_{k,j}+{\tau }^2{P}_I{\bar{c}}_k^2+{\bar{c}}_k^{(1)}+{\tilde{c}}_k+{\sigma }_n^2\frac{{\tau }_d+1}{{\tau }_d}},$$

$$c_k={\mathbf{p}}_k^T\text{diag}\left({\mathbf{\Gamma }}_k{\mathbf{C}}_k\right),c_{k,j}={\mathbf{p}}_j^T{\mathbf{\Gamma }}_k{\mathbf{R}}_j{\mathbf{p}}_j,{\tilde{c}}_k={\mathbf{p}}^T{\mathbf{\Gamma }}_k\mathbf{R}\mathbf{p},$$

$${\bar{c}}_k={\bar{\mathbf{p}}}^T\text{diag}\left({\mathbf{\Gamma }}_k{\mathbf{C}}_k\right),{\bar{c}}_k^{(1)}={\bar{\mathbf{p}}}^T{\mathbf{\Gamma }}_k{\mathbf{R}}_k^{(1)}\bar{\mathbf{p}},\text{and}$$

$${\mathbf{R}}_k^{(1)}={\bar{\mathbf{R}}}_k-{\tau }^2{P}_I{\mathbf{C}}_k^2{\mathbf{\Gamma }}_k.$$

$${\underline{R}}_k={\text{log}}_2\left(1+{\underline{\text{SINR}}}_k\right)\stackrel{N\to \infty }{<}{R}_k$$

$$y_{E_k}=b_k{q}_k+\sum _{j\ne k}{b}_j{q}_j+{\hat{b}}_{k}z+b+\bar{n},$$

$$b_j={\mathbf{g}}_E^T{\mathbf{w}}_j,{\hat{b}}_k={\mathbf{g}}_E^T{\bar{\mathbf{w}}}_k,b={\mathbf{g}}_E^T\mathbf{w}.$$

$${\overline{R}}_{E_k}={\text{log}}_2\left(1+\mathbb{E}\left[{\text{SINR}}_{E_k}\right]\right)\ge$$

$$R_{E_k}=\mathbb{E}\left[{\text{log}}_2\left(1+{\text{SINR}}_{E_k}\right)\right],$$

$$\mathbb{E}\left[{\text{SINR}}_{E_k}=\frac{{\left|b_k\right|}^2}{|{\hat{b}}_k{|}^2+|b{|}^2+{\sigma }_n^2}\right]\stackrel{N\to \infty }{\to }$$

$$\frac{\mathbb{E}\left[{\left|b_k\right|}^2\right]}{\mathbb{E}\left[|{\hat{b}}_k{|}^2+|b{|}^2+{\sigma }_n^2\right]}=\frac{{\tau }^2\varphi P_E{d}_k^2+d_k^{(1)}}{{\tau }^2\varphi P_E{\bar{d}}_k^2+{\bar{d}}_k^{(1)}+d+{\sigma }_n^2},$$

$$d_k={\mathbf{p}}_k^T\text{diag}\left(\mathbf{\Gamma }{\mathbf{C}}_k\right),d_k^{(1)}={\mathbf{p}}_k^T\mathbf{\Gamma }{\mathbf{R}}_k^{(2)}{\mathbf{p}}_k,d={\mathbf{p}}^T\mathbf{\Gamma }\mathbf{R}\mathbf{p},$$

$${\bar{d}}_k={\bar{\mathbf{p}}}^T\text{diag}\left(\mathbf{\Gamma }{\mathbf{C}}_k\right),{\bar{d}}_k^{(1)}={\bar{\mathbf{p}}}^T\mathbf{\Gamma }{\mathbf{R}}_k^{(2)}\bar{\mathbf{p}},\text{and}$$

Full text

I Introduction

In contrast to multi-cell massive multiple-input multiple-output (MIMO) systems in which the users in each cell (of a confined area) are served by an array of colocated antennas, cell-free massive MIMO is an architecture in which the users over a large area are served by a large number of distributed antennas (access points (APs)) [1]. Given the provision of backhaul phase-coherent cooperation between the APs [2, 3, 4], the distributed deployment of the APs offers many advantages such as: eliminating the correlation between the transmitting antennas, the ability to overcome deep shadow fading, and more importantly, the large freedom in balancing the simultaneous transmissions of information, jamming and energy signals.

In massive MIMO systems, the asymptotic orthogonality between independent users’ channels makes downlink transmission very robust against passive eavesdropping attacks [5]. Therefore, the active eavesdropping attack in massive MIMO systems (which introduces correlation between the estimated channels of both the attacker and the attacked user) is relevant. Active information-eavesdropping relies on attacking the uplink channel estimation phase by sending an identical training sequence as the legitimate information user (IU), such that the estimated IU’s channel is correlated with the channel of the attacking eavesdropper (EV). Therefore, the active EV benefits from the downlink transmission which is beamformed based on the estimated IU’s channel [5, 6].

The broadcast nature of the wireless channel imposes challenges in securing wireless communication systems, particularly, in the presence of adversarial EVs [7]. One example of such systems is simultaneous wireless information and power transfer (SWIPT) systems that comprise information-untrusted EHs. The secrecy issue in SWIPT massive MIMO systems, particularly under active attack, has previously lacked in-depth study in the literature. The main body of research concerning the secrecy problems in SWIPT systems has considered the colocated massive MIMO architecture [8, 9, 10, 11, 12, 13]. The large dimensionality of transmit antennas in massive MIMO systems allows the use of random matrix theory to simplify the system design and performance analysis. Moreover, the asymptotic orthogonality between independent users’ channels encourages the use of artificial noise (AN) jamming against any potential information eavesdropping. In [8], an asymptotic expression for the ergodic secrecy rate (ESR) of one IU and one passive information-untrusted energy harvester (EH) (both have multiple antennas) is derived in terms of the covariance matrix of the downlink signal vector. This asymptotic ESR is maximized by optimizing the covariance matrix subject to some average harvested energy (AHE) constraints. The AN jamming can be deployed in the downlink transmission phase to provide direct power transfer and to degrade the information signal quality at the EHs [9]. In [12], the use of AN is extended for both the downlink training and payload data transmission phases to further degrade the eavesdropping capabilities of the information EV. The authors in [13] considered joint enhancement of the secrecy and power transfer in the presence of an active dual-antenna information-untrusted EH. Asymptotic expressions for a lower bound on the ESR and the AHE are derived. Then, these results are used to optimize the power allocation for the downlink SWIPT transmission. Throughout the literature, much of the research regarding optimizing the performance of cell-free MIMO systems deals with the spectral efficiency [2, and the references therein], the energy efficiency [14, 15, 16, 17], and the secrecy rate of wire-taped systems [18].

This paper investigates the design and the performance evaluation of SWIPT in cell-free massive MIMO, particularly, the secrecy of the information transmission under an active attack from a dual-antenna information-untrusted EH. From the service provider (cooperative APs) point of view, the dual-antenna active EH’s request for service equivalently appears as a separate legitimate EH using a training power $\phi P_{E}$ (where $0<\phi<1$ and $P_{E}$ is the total available training power) via the energy harvesting antenna, and illegitimate active EV attacking a certain IU with training power $(1-\phi)P_{E}$. However, the cooperative APs can rely on their large dimensionality to monitor the levels of training powers, therefore, they can blame the legitimate EH for the active attack. Upon the detection of the active attack, the cooperative APs have no option but to deal with this attack, and only two possible actions might be taken: 1) Dropping the IU under attack from service, i.e., stop sending information to the IU being attacked. With an exception for IUs receiving information with a high degree of importance, such an action seems impractical. Therefore, there is no secrecy design for the downlink transmission; 2) Dealing with the case by optimising the secrecy of the downlink transmission. Taking this action is useful and practical, particularly with the advantage of the large dimensionality of the APs.

Contributions: We are motivated by the lack of literature on the security of cell-free MIMO systems to provide a new globally optimal solution to the problem of joint power and data transfer in a cell-free massive MIMO system. The proposed system established by a large number of randomly (Poisson-distributed) located APs which cooperate via a central processing unit (CPU). The communication links between the APs and the IUs are vulnerable to be wire-tapped by an information-untrusted dual-antenna active EH. Since the transmit power constraint is per AP, the secrecy rate is nonlinear in terms of the transmit power elements of the APs and that imposes new challenges in formulating a convex power control problem for the downlink transmission. The main contributions of our work are: 1) To jointly improve the ESR and the AHE (of the legitimate EH), we propose optimized downlink transmissions of three different signals: information, AN and energy signals beamformed towards the IUs, legitimate and illegitimate antennas of the EH, respectively; 2) We derive closed-form expressions for the AHE and a tight lower bound on the ESR. The derived expressions are deterministic at the CPU and take into account the IUs’ knowledge attained by downlink effective precoded-channel training; 3) Knowing that the ESR is nonlinear in terms of the transmit power elements of the APs, a new globally optimal iterative method for cooperatively balancing the transmit powers at the APs via relaxed semidefinite programming (SDP) is derived; 4) We provide a proof for the rank-one global optimality of our SDP solution (Theorem 3) and the convergence of our iterative SDP problem (Subsection IV-C2); 5) Finally, a fair performance comparison between the proposed cell-free and colocated massive MIMO systems is performed. The comparison shows informative results of the secrecy performance with respect to the active eavesdropping training power and the range of the AHE constraint values.

Related Work: To the best of the authors’ knowledge, the secrecy performance in cell-free massive MIMO systems has only been studied in [18] where the focus was on maximizing the secrecy rate of a given IU when being attacked by an active EV under constraints on the individual rates of all IUs. We can compare the work in this paper to the work in [18] from two perspectives: 1) From system and signal design perspectives, our work considers the worst-case SWIPT problem by optimizing three different downlink signals: information, AN and energy signals beamformed towards the IUs, legitimate and illegitimate antennas of the dual-antenna EH, respectively; while work in [18] considers the secrecy problem of a certain IU by optimizing the downlink information signals (no jamming or power transfer are considered); 2) From a problem-solving perspective, the employed lower bound on the secrecy rate in [18] imposes constraints on the domain of the linear programming (LP) optimization variables (the allocated power of the downlink information vectors) [18, (23)], i.e., the values of allocated power vectors are feasible on a sub-region of $\mathcal{R}_{+}^{N}$, $N$ is the total number of APs. Since the update in the proposed iterative algorithm does not include the power vector of the considered IU, the obtained solution is locally optimal, or at least, the globally optimal solution is not guaranteed. In contrast, in our work, both the objective function and constraints of the SDP formulation are differentiable and there are no constraints on the domain of the optimization variables which implies the satisfaction of Slater’s condition. Therefore, by proving the optimal rank requirements (please see Theorem 3 and its proof) and the convergence of the employed iterative problem (please see Subsection IV-C2), we claim the global optimality of our solution. In our early work in [13], an active dual-antenna information-untrusted EH (equivalent to the proposed EH in this paper) has been considered for a colocated SWIPT massive MIMO system. However, considering such a secrecy problem for cell-free massive MIMO will result in a non-linear objective function in terms of the allocated power elements at the APs. Inevitably, this problem can not be solved by the LP method used for a colocated massive MIMO in [13], and this leads to a completely different SDP optimization challenge.

Notation: For referencing convenience, the notations used in this paper are listed in Table I at the top of the next page.

II System Model

As illustrated in Fig. 1, we consider the downlink of a cell-free massive MIMO system consisting of a large number of APs which are randomly located on a two dimensional Euclidean area $\text{A}_{a}$ based on an homogeneous Poisson point process (PPP) $\Phi_{a}$ with an intensity $\lambda_{a}$; $M$ single antenna IUs interested in information decoding, $\{\text{IU}_{i}\}$, $i=1,2,...,M$; and an active information-untrusted EH, equipped with two antennas, where one antenna is used to legitimately harvest energy, while the other antenna is used to illegitimately and actively eavesdrop and decode an information signal intended for a certain IU, $\text{IU}_{k}$, $k\in\{1,2,...,M\}$. Unless otherwise stated, the IUs and the EH are randomly located on a two dimensional Euclidean area $\text{A}_{u}<\text{A}_{a}$111Since each user (IU or EH) is dominantly served by a subset of the APs. Therefore, the assumption $\text{A}_{u}<\text{A}_{a}$ introduces an overlap between the dominant AP groups serving different users. From the secure SWIPT design point of view, this case is more severe than the case when the users are widely apart, i.e., $\text{A}_{u}=\text{A}_{a}$.. The origins of both $\text{A}_{u}$ and $\text{A}_{a}$ coincide. The APs are networked by a centralized infinite backhaul which allows them to synchronize and cooperate via a CPU.

Let $\{\text{AP}_{1},\;\dots,\text{AP}_{N}\}$ be the set of the adopted realization of APs. ${\boldsymbol{h}_{i}=[h_{i,1},\;\dots,h_{i,N}]^{T}=\boldsymbol{\Gamma}_{i}^{\frac{1}{2}}\boldsymbol{\bar{h}}_{i}}$ denotes the uplink channel vector between $\text{IU}_{i}$ and the set of APs, where ${\boldsymbol{\bar{h}}_{i}\sim\mathcal{CN}(\boldsymbol{0},\boldsymbol{I}_{N})}$ is the small-scale fading vector and ${\boldsymbol{\Gamma}_{i}=\text{diag}([\gamma_{i,1},\;\dots,\gamma_{i,N}])}$, $\gamma_{i,j}$ is the large-scale fading coefficient of the channel between $\text{IU}_{i}$ and $\text{AP}_{j}$. ${\boldsymbol{g}=[g_{1},\;\dots,g_{N}]^{T}=\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{\bar{g}}}$ and ${\boldsymbol{g}_{E}=[g_{E_{1}},\;\dots,g_{E_{N}}]^{T}=\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{\bar{g}}_{E}}$ denote the uplink channel vectors between the legitimate and the illegitimate (eavesdropping) antennas of the EH and the set of APs, respectively, where ${\boldsymbol{\bar{g}}}=[\bar{g}_{1},\;\dots,\bar{g}_{N}]^{T},~{}{\boldsymbol{\bar{g}}_{E}}=[\bar{g}_{E_{1}},\;\dots,\bar{g}_{E_{N}}]^{T}\sim\mathcal{CN}(\boldsymbol{0},\boldsymbol{I}_{N})$ are independent, uncorrelated small-scale fading vectors. ${\boldsymbol{\Gamma}=\text{diag}([\gamma_{1},\;\dots,\gamma_{N}])}$ where $\gamma_{j}$ is the large-scale fading coefficient of the channel between the EH and $\text{AP}_{j}$. The large-scale fading coefficients $\{\gamma_{i,j},~{}\gamma_{j}\}$ change very slowly compared to the small-scale fading coefficients, therefore, we assume that $\{\gamma_{i,j},~{}\gamma_{j}\}$ are perfectly known at the APs [19].

II-A Uplink Channel Estimation

The user small-fading channels manifest block fading, i.e., they remain constant over one time block, but change independently from one block to another. Each time block is divided into three time slots of lengths: $\tau$ transmission samples for uplink training, $\tau_{d}$ transmission samples for downlink training and $\tau_{s}$ samples for downlink data transmission. Without loss of generality, we assume a unit time slot for the downlink data transmission $\tau_{s}T_{s}=1s$, where $T_{s}$ is the duration of the transmitted data symbol [8, 20]. During the uplink training phase, a training sequence is sent from each IU with an average power $P_{I}$. Pessimistically, we assume that the EH has the potential to acquire the training sequence of a certain IU (made possible by overhearing the leaking electromagnetic signalling between the APs and the IUs [21]). Therefore, the EH sends a copy of the training sequence of the attacked IU, $\text{IU}_{k},~{}k\in\{1,2,~{}\dots,M\}$, via its eavesdropping antenna using part of its total average power $\phi P_{E}$, $0<\phi<1$, such that the cooperative APs estimate the uplink composite channel coefficients of both $\text{IU}_{k}$ and the eavesdropping antenna of the EH. Consequently, the estimated channel of $\text{IU}_{k}$ will be corrupted and correlated with the illegitimate channel of the EH [5, 22]. The remaining training power $(1-\phi)P_{E}$ is used for transmitting the legitimate uplink training sequence via the energy harvesting antenna. The uplink training sequences of the IUs and legitimate EH are assumed to be orthogonal. The signal at the APs received across $\tau$ training transmissions is

[TABLE]

where $\boldsymbol{N}\in\mathcal{C}^{N\times\tau}$ is the additive noise matrix with entries following the distribution $\mathcal{CN}(0,\sigma_{n}^{2})$. $k$ is the index of the attacked IU, $\text{IU}_{k}$. $\boldsymbol{\psi}_{i},~{}\boldsymbol{\psi}_{k},~{}\boldsymbol{\psi}_{E}\in\mathcal{C}^{\tau\times 1}$ are the uplink training sequences of $\text{IU}_{i}$, the IU under attack, $\text{IU}_{k}$, and the legitimate antenna of the EH, respectively. $\boldsymbol{\psi}_{i}^{H}\boldsymbol{\psi}_{j\neq i},\;\boldsymbol{\psi}_{i}^{H}\boldsymbol{\psi}_{E}=0$; and $\boldsymbol{\psi}_{i}^{H}\boldsymbol{\psi}_{i},\;\boldsymbol{\psi}_{E}^{H}\boldsymbol{\psi}_{E}=\tau$. We assume centralized channel estimation via the CPU. Given that $\text{IU}_{k}$ is the attacked IU, the minimum mean square error (MMSE) estimate of $\boldsymbol{h}_{i}$, $\boldsymbol{\hat{h}}_{i}=[\hat{h}_{i,1},\;\dots,\hat{h}_{i,N}]^{T}$, and of $\boldsymbol{g}$, $\boldsymbol{\hat{g}}=[\hat{g}_{1},\;\dots,\hat{g}_{N}]^{T}$, are given as

[TABLE]

where $\delta_{ik}=1$ if $i=k$ (i.e., $\text{IU}_{i}$ is the attacked IU) and $\delta_{ik}=0$ if $i\neq k$. The covariance matrices $\mathbb{E}[\boldsymbol{\hat{h}}_{i}\boldsymbol{\hat{h}}_{i}^{H}]$ and $\mathbb{E}[\boldsymbol{\hat{g}}\boldsymbol{\hat{g}}^{H}]$ are equal to $\boldsymbol{R}_{i}=\tau\sqrt{P_{I}}\boldsymbol{\Gamma}_{i}\boldsymbol{C}_{i}$ and $\boldsymbol{R}=\tau\sqrt{(1-\phi)P_{E}}\boldsymbol{\Gamma}\boldsymbol{C}$, respectively. To emphasize whether $\text{IU}_{i}$ is being attacked or not, we use $\boldsymbol{R}_{i}$ to describe the covariance matrix of $\text{IU}_{i}$ if not being attacked and $\boldsymbol{\bar{R}}_{i}$ to describe the covariance matrix of $\text{IU}_{i}$ if being attacked. Both $\boldsymbol{R}_{i}$ and $\boldsymbol{\bar{R}}_{i}$ are calculated by the same aforementioned formula, but with $k\neq i$ for $\boldsymbol{R}_{i}$ and with $k=i$ for $\boldsymbol{\bar{R}}_{i}$. The results in (2a) and (2c) follow from standard channel estimation theory [23, 24]. Active eavesdropping attack detection and the identification of the attacked IU, $\text{IU}_{k}$, are possible and have been studied in [25, 26, 27]. Alternatively, the cooperative APs can exploit their large dimensionality to detect the active eavesdropping attack by monitoring the values of training powers which have been proven to be accurate as $N\to\infty$. The CPU can calculate the eavesdropping (illegitimate) and the legitimate training powers of the EH, $\phi P_{E}$ and $(1-\phi)P_{E}$, respectively, by using the following lemma222Since the cooperative APs are able to monitor the changes in the training powers of the IUs and the EH using Lemma 1, we assume that the cooperative APs blame the information-untrusted EH for the active eavesdropping attack.

Lemma 1

For a large density of APs as $\lambda_{a}\to\infty$, which leads to a large number of APs as $N\to\infty$, any illegitimate active training power can be identified and calculated as

[TABLE]

where $\text{IU}_{i}$ is under attack if $\delta_{ik}=1$, i.e., $k=i$, and $\text{IU}_{i}$ is not being attacked if $\delta_{ik}=0$, i.e., $k\neq i$. All the scalars, vector and matrices in the left-hand side of (3) are deterministic at the CPU.

Proof:

See Appendix A. ∎

II-B Downlink Transmission

The APs cooperate via the CPU to control the power allocation of the downlink data, AN, and energy signal transmissions. From the service provider (cooperative APs) point of view, the EH’s request for service equivalently appears to the cooperative APs as a separate legitimate EH which uses a training power $\phi P_{E}$ and illegitimate active eavesdropper attacking a certain IU, $\text{IU}_{k}$, with a training power $(1-\phi)P_{E}$. However, the CPU relies on the large dimensionality of the APs to monitor the levels of training powers, and based on Lemma 1, it can blame the legitimate EH for the active attack. Upon the detection of the active attack, the CPU has no option but to deal with this attack, and only two possible actions might be taken:

• •

Dropping the IU under attack from service, i.e., stop sending information to the IU being attacked. With an exception for IUs receiving information with a high degree of importance, such an action seems impractical. Therefore, there is no secrecy design for the downlink transmission.

• •

Dealing with the case by optimizing the secrecy of downlink transmission (by employing controlled transmissions of information, jamming and energy signals). Taking this action is useful and practical, particularly with the advantage of the large number of randomly located APs. Compared to the case of collocated APs (conventional MIMO), the average path-loss from an AP to the active EH and the attacked IU varies from one AP to another. This property of randomly distributed APs would increase the efficiency of power control in tackling the active eavesdropping.

Given that the $\text{IU}_{k}$ is the attacked IU, the APs employ the matched filter (MF) precoder to transmit the downlink signal vector

[TABLE]

where the $j$th entry of $\boldsymbol{x}_{k}$, $[\boldsymbol{x}_{k}]_{j}$, is the signal transmitted by $\text{AP}_{j}$, $\boldsymbol{w}_{i}q_{i}$ is the information signal vector directed towards $\text{IU}_{i}$, $\boldsymbol{\bar{w}}_{k}z$ is the AN signal vector directed towards the eavesdropping antenna of the EH, and $\boldsymbol{w}$ is the energy signal vector directed towards the legitimate antenna of the EH. $\{q_{i}\}$ and $z$ are the information signal symbols intended for $\{\text{IU}_{i}\}$ and the AN symbol, respectively, and they are mutually independent and follow the distribution $\mathcal{CN}(0,1)$. The MF beamforming vectors in (4) are defined as333Please note, due to active attack, the $\boldsymbol{\hat{h}}_{k}$ used to design $\boldsymbol{\bar{w}}_{k}$ in (5b) is the estimate of the composite channel of both $\boldsymbol{h}_{k}$ and $\boldsymbol{g}_{E}$. By optimizing the per AP AN power factors $\{\bar{p}_{1},\;\dots,\;\bar{p}_{N}\}$, the AN power can be maximized at the EH and minimized at the $\text{IU}_{k}$.

[TABLE]

For example, $|\left[\boldsymbol{w}_{i}\right]_{j}|^{2}=p_{i,j}|\hat{h}_{i,j}|^{2}$, $|\left[\boldsymbol{\bar{w}}_{k}\right]_{j}|^{2}=\bar{p}_{j}|\hat{h}_{k,j}|^{2}$ and $|\left[\boldsymbol{w}\right]_{j}|^{2}=p_{j}|\hat{g}_{j}|^{2}$ are the allocated powers at $\text{AP}_{j}$ for $\text{IU}_{i}$’s data, AN and energy signals, respectively. Power allocation is controlled via the factors $\{p_{i,j}\}$, $\{\bar{p}_{j}\}$ and $\{p_{j}\}$. Referring to (2a) and (5b), it can be noticed that the received AN signal power at the eavesdropping antenna of the EH, $|\boldsymbol{g}_{E}^{T}\boldsymbol{\bar{w}}_{k}|^{2}$, is directly proportional to the eavesdropping training power, $\phi P_{E}$, i.e., the larger the eavesdropping training power, the larger the jamming received power by the EH. Therefore, although the AN is aligned to the $\text{IU}_{k}$’s estimated channel coefficients, the cooperative APs can improve the information secrecy by exploit the nature of the cell-free system – in which $\text{IU}_{k}$ and the EH experience different path-losses to a single AP – by optimizing the per AP per user power control.

Given that $\text{IU}_{k}$ is the attacked IU. The received signals at $\text{IU}_{i}$, $y_{k,i}$; the legitimate antenna of the EH, $y_{k}$; and at the eavesdropping antenna of the EH, $y_{E_{k}}$, are

[TABLE]

where $n_{i}$, $\hat{n}$ and $\bar{n}$ are zero mean $\sigma_{n}^{2}$ variance complex Gaussian noises at $\text{IU}_{i}$, the legitimate and eavesdropping antennas of the EH, respectively.

II-C Downlink Effective Precoded-Channel Estimation

With a large number of APs, the channel estimation at all IUs requires training sequences of a length $\geq N$ which is practically infeasible. Alternatively, we propose the estimation of the effective precoded-channels, $\{a_{i,i}=\boldsymbol{h}_{i}^{T}\boldsymbol{w}_{i}\}$ at the IUs444The EH has the potential to estimate the precoded channel for the attacked IU, $b_{k}=\boldsymbol{g}_{E}^{T}\boldsymbol{w}_{k}$, however, as will be seen in Subsection III-B, the worst case in which the EH can perfectly estimate $b_{k}$ is assumed.. The downlink estimation of the effective precoded-channels at the IUs requires $M$ orthogonal training sequences that can be of a finite length, $\geq M$. Therefore, such a downlink estimation is practically possible. Notice that $\text{IU}_{i}$ needs to estimate its effective precoded-channel $a_{i,i}$ which includes the values of power control factors $\{p_{i,j}\}$, $\{\bar{p}_{j}\}$ and $\{p_{j}\}$. Therefore the values of $\{p_{i,j}\}$, $\{\bar{p}_{j}\}$ and $\{p_{j}\}$ to be used for downlink data transmission are employed for downlink training. The cooperative APs transmit the downlink training signal matrix $\boldsymbol{X}_{d}=\sum_{i=1}^{M}\boldsymbol{w}_{i}~{}\boldsymbol{\psi}_{d_{i}}^{T}$, where $\{\boldsymbol{\psi}_{d_{i}}\}$ – $\boldsymbol{\psi}_{d_{i}}^{H}\boldsymbol{\psi}_{d_{i}}=\tau_{d}$ and $\boldsymbol{\psi}_{d_{i}}^{H}\boldsymbol{\psi}_{d_{j\neq i}}=0$ – are the downlink training sequences of the IUs555 The same training sequences could be used in the uplink and downlink.. The received training signal vector at $\text{IU}_{i}$, $\boldsymbol{y}_{I_{i}}\in\mathcal{C}^{1\times\tau_{d}}$ is

[TABLE]

where $a_{i,j}=\boldsymbol{h}_{i}^{T}\boldsymbol{w}_{j}$ and $\boldsymbol{n}_{i}\sim\mathcal{CN}(0,\sigma_{n}^{2}\boldsymbol{I}_{\tau_{d}})$ is the noise vector at $\text{IU}_{i}$. First, let us examine the MMSE estimate of $a_{i,i}$ at $\text{IU}_{i}$ which can be calculated as [23, 24]

[TABLE]

where $y_{I_{i}}=\boldsymbol{y}_{I_{i}}\boldsymbol{\psi}_{d_{i}}^{*}=\tau_{d}a_{i,i}+\boldsymbol{n}_{i}\boldsymbol{\psi}_{d_{i}}^{*}$. However, since the allocated power control factors in $\boldsymbol{p}_{i}$ are not available at $\text{IU}_{i}$, the calculation of (8) is not possible, and instead, we assume that $\text{IU}_{i}$ performs a simple least square error (LSE) estimate of $a_{i,i}$, $\hat{a}_{i,i}$ which is given as

[TABLE]

where $\tilde{a}_{i,i}=\frac{\boldsymbol{n}_{i}\boldsymbol{\psi}_{d_{i}}^{*}}{\tau_{d}}$ is the estimation error which is statistically independent from the effective precoded channel $a_{i,i}$.

III Secrecy Analysis

III-A Lower Bound on the IU Rate

The received signal at $\text{IU}_{i}$ given in (6a) can be recast as follows

[TABLE]

where

[TABLE]

$\mathbb{E}\left[a_{i,i}|\hat{a}_{i,i}\right]q_{i}$ is the desired information signal received through a deterministic precoded channel $\mathbb{E}\left[a_{i,i}|\hat{a}_{i,i}\right]$, while $\left(a_{i,i}-\mathbb{E}\left[a_{i,i}|\hat{a}_{i,i}\right]\right)q_{i}$ is the desired information signal received through a non-deterministic precoded channel $a_{i,i}-\mathbb{E}\left[a_{i,i}|\hat{a}_{i,i}\right]$. $\mathbb{E}\left[a_{i,i}|\hat{a}_{i,i}\right]q_{i}$ and $\left(a_{i,i}-\mathbb{E}\left[a_{i,i}|\hat{a}_{i,i}\right]\right)q_{i}$ are statistically dependent. $Z_{k,i}$ is the equivalent noise666$Z_{k,i}$ is considered as an equivalent noise since $Z_{k,i}$ and $\mathbb{E}\left[a_{i,i}|\hat{a}_{i,i}\right]q_{i}$ are independent and that follows since $\{q_{j}\}$, $z$, $n_{i}$ and $\boldsymbol{w}$ are statistically independent. which accounts for inter user interference, energy signal interference and the thermal noise. Referring to (9), we can see that $a_{i,i}$ is explicitly decoupled and therefore $a_{i,i}$ and $\tilde{a}_{i,i}$ are uncorrelated and statistically independent. Since $\hat{a}_{i,i}$ is deterministic at $\text{IU}_{i}$, then

[TABLE]

where $\mathbb{E}\left[\hat{a}_{i,i}|\hat{a}_{i,i}\right]=\hat{a}_{i,i}$ follows as an expectation over a deterministic value; $\mathbb{E}\left[\tilde{a}_{i,i}|\hat{a}_{i,i}\right]=\mathbb{E}\left[\tilde{a}_{i,i}\right]$ follows from the statistical independence between $\tilde{a}_{i,i}$ and $\hat{a}_{i,i}$; and $\mathbb{E}\left[\tilde{a}_{i,i}\right]=0$ follows since $\mathbb{E}\left[\boldsymbol{n}_{i}\boldsymbol{\psi}_{d_{i}}^{*}\right]=0$. Using the results in [28, Theorem 1] and in [29, (22)], the downlink information rate at the attacked user $\text{IU}_{k}$, $R_{k}$ (given in (13)) is achievable and forms a lower bound on the ergodic information rate

[TABLE]

where

[TABLE]

Theorem 1

For $N\rightarrow\infty$, the value of $\text{SINR}_{k}$ is tightly lower bounded by a deterministic value $\underline{\text{SINR}}_{k}\overset{N\rightarrow\infty}{<}\text{SINR}_{k}$ which is given by

[TABLE]

where

[TABLE]

Since $\underline{\text{SINR}}_{k}$ is deterministic (independent of the small-fading randomness, $\mathbb{E}[\underline{\text{SINR}}_{k}]=\underline{\text{SINR}}_{k}$), and based on (13) and (15), $\underline{R}_{k}=\text{log}_{2}(1+\underline{\text{SINR}}_{k})$ is a tight lower bound on the ergodic rate of the attacked user $\text{IU}_{k}$, and known at the CPU.

[TABLE]

Proof:

See Appendix A. ∎

III-B Upper Bound on the EH Ergodic Rate

The received signal at the eavesdropping antenna of the EH in (6c) can be recast as follows

[TABLE]

In the following, we assume the worst-case scenario in which the EH has full knowledge of its own channel vectors, $\boldsymbol{g}_{E}$ and $\boldsymbol{g}$; and the beamforming vectors $\{\boldsymbol{w}_{i}\}$. With this worst-case assumption, an upper bound on the ergodic information rate at the EH is given in the following theorem.

Theorem 2

With a worst-case scenario assumption that the EH has full knowledge of its own channel and the beamforming vectors of the IUs, the EH is capable of cancelling the inter-user interference [30, Chapter 8]. Since the information, $\{q_{i}\}$, the AN signal, $z$, and the energy signal, $\boldsymbol{w}$, are statistically independent, we have the following upper bound, $\overline{R}_{E_{k}}$, on the ergodic rate of the EH intending to eavesdrop $\text{IU}_{k}$, $R_{E_{k}}$, given by

[TABLE]

for which

[TABLE]

where

[TABLE]

Proof:

See Appendix A. ∎

Such a worst-case scenario is commonly employed by much of the current research to guarantee maximum information security [31, 20]. Ensuring the confidentiality of the information for the worst-case scenario design ensures confidentiality for more optimistic scenarios.

III-C Lower Bound on the Ergodic Secrecy Rate of $\text{IU}_{k}$

Using the lower bound and the upper bound on the information rates at the attacked user $\text{IU}_{k}$ and the EH given in (16) and (18), we assess the secrecy of information at $\text{IU}_{k}$ in terms of ESR which has the following lower bound

[TABLE]

III-D Average Harvested Energy at the EH

The EH relies on the dual functionality of its antennas to harvest energy and eavesdrop information simultaneously. The whole signal received via the legitimate antenna is devoted for energy harvesting, while the signal received via the illegitimate antenna is used for information decoding. However, since the CPU blames the EH for the active attack, the received signals via both antennas are accounted for the CPU for energy harvesting. The AHE by the EH intending to eavesdrop $\text{IU}_{k}$ is777Detailed derivation of the results in (LABEL:eeq22) are in Appendix A.

[TABLE]

where

[TABLE]

IV Power Control of Downlink Transmission

IV-A Problem Formulation

In our system, a single AP, $\text{AP}_{j}$, transmits a set of $M+2$ different types of signals, $\{\{[\boldsymbol{w}_{i}]_{j}q_{i}\}_{i},~{}[\boldsymbol{\bar{w}}_{k}]_{j}z,~{}[\boldsymbol{w}]_{j}\}$. With the random geometric distribution of the APs with respect to the IUs and the EH, the power control in the cell-free MIMO system has an advantage over the conventional MIMO that different users have different subsets of dominant serving APs. In the long-term, the CPU can achieve a fair and secured SWIPT transmission towards the IUs and the EH by balancing the average levels of transmit powers at the APs within the power limits of each AP. The power control aims to maximize the worst-case ESR, $\min_{k}R_{S_{k}}$, with a constraint on the minimum AHE requirement of the legitimate EH. Therefore, our constrained problem is

[TABLE]

where $P_{t}$ is the available power budget at each AP. The constraint (22b) guarantees the average power consumption at each AP is within the limit, $P_{t}$. Problem (22) is non-convex since the objective function is a logarithm of multiplicative fractional functions. Without loss of generality, we assume that (22) is always feasible and focus on solving it. We use the exponential variable substitution method used in [32] and [33] to transform the logarithmic objective function of (22) into an equivalent linear function. By using the properties of logarithmic and exponential functions, the objective function of (22) can be expressed as $\text{log}_{e}2~{}\text{ln}(e^{u_{k}-s_{k}}e^{v_{k}-t_{k}})$ where

[TABLE]

Since the logarithmic functions are monotonically increasing in their arguments, then (22) can be recast as

[TABLE]

Our new objective in (24) is monotonically increasing with $\min_{k}R_{S_{k}}$. The constraints (24a)–(24e) bound the slack variables $u_{k},~{}s_{k},~{}t_{k},~{}v_{k}$ of the objective function within their limits defined in (23a)–(23d). The exponential variables $e^{s_{k}}$ and $e^{t_{k}}$ are linearized as $e^{\bar{s}_{k}}(s_{k}-\bar{s}_{k}+1)$ and $e^{\bar{t}_{k}}(t_{k}-\bar{t}_{k}+1)$. $\bar{s}_{k},~{}\bar{t}_{k}$ are the initial values around which $e^{s_{k}}$ and $e^{t_{k}}$ are linearized.

The formulation in (24) is still non-convex since the right-hand sides of the constraints (24a)–(24e) contain expressions which are nonlinear in the optimization variables (the power control factors $\{\{\boldsymbol{p}_{i}\},~{}\boldsymbol{\bar{p}},~{}\boldsymbol{p}\}$), such as $c_{k}^{2}=\left(\boldsymbol{p}_{k}^{T}\text{diag}\left(\boldsymbol{\Gamma}_{k}\boldsymbol{C}_{k}\right)\right)^{2}$. These nonlinearities arise from the per AP per user power control (specific for cell-free massive MIMO systems) where each AP has its own transmit power constraint. In comparison, these nonlinearities do not exist in the power control for the conventional (collocated) massive MIMO systems in which the constraint is on the total transmit power from all collocated antennas [13]. To deal with these nonlinearities, we introduce a new method of cooperative balancing of the transmit powers at the APs via relaxed SDP formulation which has been proved to be optimal as will be described in the next subsection.

IV-B SDP Formulation for Optimal Power Control

In this subsection, we reformulate the non-convex problem (24) into a relaxed SDP convex problem. To achieve this, the nonlinear expressions in the power control factors $\{\{\boldsymbol{p}_{i}\},~{}\boldsymbol{\bar{p}},~{}\boldsymbol{p}\}$ are represented as linear expressions in terms of new rank-one positive semidefinite matrix variables $\{\{\boldsymbol{P}_{i}=\boldsymbol{p}_{i}\;\boldsymbol{p}_{i}^{H}\},~{}\boldsymbol{\bar{P}}=\boldsymbol{\bar{p}}\;\boldsymbol{\bar{p}}^{H},~{}\boldsymbol{P}=\boldsymbol{p}\;\boldsymbol{p}^{H}\}$. For instance, given that $k$ is the index of the IU under attack, the expression of $c_{k}^{2}$ can be recast in an SDP form as

[TABLE]

where ${\boldsymbol{A}_{k}=\text{diag}(\boldsymbol{\Gamma}_{k}\boldsymbol{C}_{k})~{}\text{diag}(\boldsymbol{\Gamma}_{k}\boldsymbol{C}_{k})^{T}}$. In a comparable way, the rest of the expressions $\{c_{k,j},~{}\bar{c}_{k}^{2},~{}\bar{c}_{k}^{(1)},~{}\tilde{c}_{k}\}$, $\{d_{k}^{2},~{}d_{k}^{(1)},~{}d,~{}\bar{d}_{k}^{2},~{}\bar{d}_{k}^{(1)}\}$ and $\{d_{k,j},~{}\tilde{d}_{k},~{}\tilde{d}^{2},~{}\tilde{d}^{(1)}\}$ in (24) can be transformed into linear expressions in terms of $\{\{\boldsymbol{P}_{i}\},~{}\boldsymbol{\bar{P}},~{}\boldsymbol{P}\}$. With these transformations, we can recast the non-convex problem in (24) into a convex relaxed888The formulation in (26) does not impose any constraints on the rank of $\{\{\boldsymbol{P}_{i}\},~{}\boldsymbol{\bar{P}},~{}\boldsymbol{P}\}$, i.e, $\{\text{rank}(\boldsymbol{P}_{i})\},~{}\text{rank}(\boldsymbol{\bar{P}}),~{}\text{rank}(\boldsymbol{P})\leq N$. SDP formulation as in (26) at the top of the next page, where $\mathbb{S}=\{\{\boldsymbol{P}_{i}\},~{}\boldsymbol{\bar{P}},~{}\boldsymbol{P},~{}\{u_{k},~{}s_{k},~{}t_{k},~{}v_{k}\}\}$ is the set of optimization variables and

[TABLE]

The constraints (26e) and (26f) are an SDP recast of (22a) and (22b), respectively. The constraint (26f) is equivalent to (22b), where $\boldsymbol{D}_{l}\in\mathcal{R}^{N\times N}$ has zero entries except $[\boldsymbol{D}_{l}]_{l,l}=1$. This equivalent representation in (26f) is required to facilitate the proof of Theorem 3 presented in Appendix A.

The formulation in (26) is convex and can be solved iteratively based on the initial value update method given in Algorithm 1. It can be shown that the complex-valued SDP problem (27) (which is equivalent to (26)) contains: $M+2$ semidefinite complex-valued $N\times N$ matrix variables, $5M+1$ real scalar variables, $6M+NM+2$ constraints on matrix variable of size $N\times N$, and $M$ constraints on scalar variables. The complexity (in terms of number of complex operations) of obtaining a per iteration solution of (26) within accuracy $\epsilon$ is asymptotically upper bounded by $\mathcal{O}(M^{4}N^{\frac{9}{2}}\log(\frac{1}{\epsilon}))$ [34]. This result assumes unstructured input data matrices. However, the optimization solver (such as SeDuMi employed by CVX software [35]) can exploit the structure of input data matrices – for example, the structure of single non-zero element matrices $\{\boldsymbol{D}_{l}\}$ – to reduce the computational complexity [34].

IV-C Global Optimality of the SDP Formulation

To investigate the optimality of the solution obtained by (26), let us rewrite (26) in the equivalent form in (27) by replacing the objective $\min_{k}\left(u_{k}-s_{k}+v_{k}-t_{k}\right)$ by a new slack viable $\pi$ and $K$ linear constraints as

[TABLE]

By examining (27) with the first-order and the second-order conditions of convexity, we have

[TABLE]

This means that (26) is convex with an affine objective function. Since the constraints of (26) are differentiable and there are no constraints on the domain of the optimization variables $\{\boldsymbol{P}_{i}\},~{}\boldsymbol{\bar{P}},~{}\boldsymbol{P}\in\mathcal{S}^{+}$, $\{u_{k},~{}s_{k},~{}t_{k},~{}v_{k}\},~{}\pi\in\mathcal{R}$, then Slater’s condition holds and the solution obtained by solving (26) is globally optimal subject to: 1) satisfying the rank requirement of $\{\boldsymbol{P}_{i}\},~{}\boldsymbol{\bar{P}}~{}~{}\text{and}~{}~{}\boldsymbol{P}$; 2) and the convergence of the constraints (26b) and (26c) (which results in the convergence of the iterative problem (26)).

IV-C1 Rank-one Optimality

Generally, the optimality of the solutions obtained via SDP programming might require a rank higher than one. The rank requirement for the optimality of the solutions obtained by SDP problems has been investigated in [36, Lemma 3.1] which is quoted as:

Lemma 2

Suppose that the separable SDP (P1) and its dual (D1) are solvable. Then, problem (P1) has always an optimal solution $\{\boldsymbol{X}_{1}^{\star},~{}\dots,~{}\boldsymbol{X}_{L}^{\star}\}$ such that

[TABLE]

where $\{\boldsymbol{X}_{1},~{}\dots,~{}\boldsymbol{X}_{L}\}$ are the semi-definite matrix variables of (P1), $\{\boldsymbol{X}_{1}^{\star},~{}\dots,~{}\boldsymbol{X}_{L}^{\star}\}$ are their optimal values and $M$ (for (29) only) is the number of constraints. Nevertheless, for our problem (26), the obtained solution $\{\boldsymbol{P}_{i}^{\star}\},~{}\boldsymbol{\bar{P}}^{\star}$ and $\boldsymbol{P^{\star}}$ needs to satisfy the rank-one structure $\{\text{rank}(\boldsymbol{P}_{i}^{\star})\},~{}\text{rank}(\boldsymbol{\bar{P}}^{\star}),~{}\text{rank}(\boldsymbol{P}^{\star})=1$ which complies with the optimality condition (29) in Lemma 2. The compliance of $\{\boldsymbol{P}_{i}^{\star}\},~{}\boldsymbol{\bar{P}}^{\star}$ and $\boldsymbol{P^{\star}}$ with rank-one requirement is given in the following theorem

Theorem 3

Given that $\mathbb{S}^{\star}=\{\{\boldsymbol{P}_{i}^{\star}\},~{}\boldsymbol{\bar{P}}^{\star},~{}\boldsymbol{P^{\star}},~{}\{u_{k}^{\star},~{}s_{k}^{\star},~{}t_{k}^{\star},$ $v_{k}^{\star}\}\}$ is the solution obtained by solving (26), then, the optimized power control factor matrices $\{\boldsymbol{P}_{i}^{\star}\},~{}\boldsymbol{\bar{P}}^{\star},~{}\boldsymbol{P^{\star}}$ always satisfy the rank-one constraint, i.e., $\{\text{rank}(\boldsymbol{P}_{i}^{\star})\},~{}\text{rank}(\boldsymbol{\bar{P}}^{\star}),~{}\text{rank}(\boldsymbol{P}^{\star})=1$.

Proof:

See Appendix B. ∎

IV-C2 Convergence of the Iterative Problem

Here, we prove that the iterative optimization (26) converges to a globally optimal value, and the objective value (which is monotonically increasing with $\min_{k}R_{S_{k}}$) is increasing with the iterations. To facilitate our proof, let us introduce the following results:

Lemma 3

For arbitrary real values of $x$ and $\bar{x}\neq x$, the first order approximation $e^{\bar{x}}\left(x-\bar{x}+1\right)$ is always an underestimate of $e^{x}$, i.e.

[TABLE]

Proof:

See Appendix C. ∎

Lemma 4

For arbitrary real values of $x$ and $\bar{x}\neq x$, the successive first order approximations of $e^{x}$; $f^{[n]}=e^{\bar{x}^{[n]}}\left(x-\bar{x}^{[n]}+1\right)$ and $f^{[n+1]}=e^{\bar{x}^{[n+1]}}\left(x-\bar{x}^{[n+1]}+1\right)$, $e^{\bar{x}^{[n+1]}}=f{[n]}$; always satisfy $f^{[n+1]}>f^{[n]}$ for $\bar{x}^{[n]}\neq x$.

Proof:

See Appendix C. ∎

Without loss of generality, we assume that (26) is feasible in its first iteration. Since our problem is convex and Slater’s condition holds (see (28) and the paragraph that follows), constraints (26b) and (26c) can strictly hold. With the first order linearisation in (26b) and (26c), and according to Lemma 3, the constraints (26b) and (26c) are tighter than their original formulations in (24b) and (24c), i.e., the feasibility region of (24) is smaller than and a subregion of the feasibility region of (26). Therefore, any non-converged solution is suboptimal.

According to Lemma 4, and since the constraints (26b) and (26c) are initialized in the $n\text{th}$ iteration by the optimal values obtained at the $(n-1)\text{th}$ preceding iteration such as $e^{\bar{s}_{k}^{[n]}}=e^{\bar{s}_{k}^{[n-1]}}(s_{k}^{\star^{[n-1]}}-\bar{s}_{k}^{[n-1]}+1)$ and $e^{\bar{t}_{k}^{[n]}}=e^{\bar{t}_{k}^{[n-1]}}(t_{k}^{\star^{[n-1]}}-\bar{t}_{k}^{[n-1]}+1)$, $\forall~{}k$, the feasibility of the $(n-1)\text{th}$ iteration will ensure the feasibility of the succeeding $n\text{th}$ iteration. Furthermore, the feasibility region at the $n\text{th}$ iteration is larger than the feasibility region at the $(n-1)\text{th}$ iteration and contains it. This ensures that the optimized objective value is monotonically increasing with the successive iteration. Given that the constrained values in (26b) and (26c) are finitely bounded (since both constraints are linear and the power budget at every AP is finite, $\leq P_{t}$), therefore, it can be concluded that the increasing optimized objective value will certainly converge, let us say at the $n$th iteration, i.e.

[TABLE]

By solving the updating method, $e^{\bar{s}_{k}^{[n]}}=e^{\bar{s}_{k}^{[n-1]}}(s_{k}^{\star^{[n-1]}}-\bar{s}_{k}^{[n-1]}+1)$ and (31), we have $s_{k}^{\star^{[n]}}=\bar{s}_{k}^{{[n]}}$, and $e^{\bar{s}_{k}^{{[n]}}}(s_{k}^{\star^{[n]}}-\bar{s}_{k}^{{[n]}}+1)=e^{s_{k}^{\star^{[n]}}}$. This indicates that the constraint (26b) converges to its original nonlinearized form. Likewise, the constraint (26c) converges to its original nonlinearized form.

V Evaluations

In this section, we evaluate the asymptotic performance of our SWIPT cell-free massive MIMO system. The APs are randomly located on a two dimensional Euclidean area $A_{a}$ based on an homogeneous PPP $\Phi_{a}$ with an intensity $\lambda_{a}$. The IUs and the EH are randomly located on a two dimensional Euclidean area $A_{u}$ with the origins of $A_{a}$ and $A_{u}$ coincident (please refer to footnote 1 regarding this assumption). The large-scale fading coefficients $\{\gamma_{i,j},~{}\gamma_{j}\}$ are modeled with the standard distance-based model as $\gamma_{i,j}\triangleq d_{i,j}^{-\alpha}~{}10^{\frac{\nu_{i,j}}{10}}$ and $\gamma_{j}\triangleq d_{j}^{-\alpha}~{}10^{\frac{\nu}{10}}$, where $d_{i,j}$ and $d_{j}$ are the distances from $\text{IU}_{i}$ and the EH to $\text{AP}_{j}$, respectively. $\alpha$ is the pathloss exponent and $\nu,~{}\nu_{i,j}\sim\mathcal{CN}(0,\sigma^{2})$ are the shadow fading coefficients with standard division $\sigma$. All users experience independent shadow fading, i.e., $\nu_{i,j}$ and $\nu_{i,j}$s are independent random variables (RVs). $P$, $P_{E}$ and $P_{t}$ are the training power budget at every IU, the training power budget at the EH, and the transmit power budget at every AP, respectively. $\tau$ and $\zeta$ are the length of the training sequences and the energy harvesting efficiency at the EH, respectively. Unless otherwise stated, and for referencing convenience, the selected values of system parameters are listed in Table II.

Fig. 2 shows the AP-user deployment geometry of a realization in which the number of APs is $N=145$ ($\mathbb{E}(N)=\lambda_{a}A_{a}=150$), $M=3$ IUs and one EH zoomed into the central $500\times 500$ $\text{m}^{2}$.

The SWIPT secrecy performance is presented by the E-R plot which relates the achievable worst-case ESR, $\min_{k}R_{S_{k}}$, to the constraint on the minimum AHE by the EH, $\bar{E}$. The larger the area under the E-R curve, the better the SWIPT performance. Our design analyses are made based on the asymptotic assumption $N\to\infty$, then, the system’s performance can be examined for a realistic scenario of a large but finite number of APs.

In colocated MIMO systems, the user exhibits a constant average path-loss to all of the base station’s (BS’s) colocated antennas, and that average path-loss varies from one user to another based on the user’s location. In contrast, in cell-free MIMO systems, the average path-loss of a given user varies from one AP to another. Intuitively, this property of randomly distributed APs is anticipated to increase the efficiency of power control in tackling the active eavesdropping. For fair comparisons between the performance of cell-free and colocated MIMO systems, a comparable model of single-cell colocated massive MIMO system is derived such that: 1) the total number of colocated antennas at the BS is equal to the total number of APs, $N$; 2) the average value of a user’s pathloss to the BS in colocated MIMO (all users experience equal pathlosses) is equal to the average value of the users’ pathlosses in cell-free MIMO; 3) the total transmit power is equal for both system, and the power limits at the colocated MIMO is per antenna; 4) the antennas of the BS are uncorrelated. Defining $\bar{\gamma}_{i}$ and $\bar{\gamma}$ as the pathlosses of $\text{IU}_{i}$ and the EH in the colocated MIMO system, respectively, we have $\bar{\gamma}_{i}=\sum_{j}\frac{\gamma_{i,j}}{N}$ and $\bar{\gamma}=\sum_{j}\frac{\gamma_{j}}{N}$. The downlink beamfoming and power control of the colocated MIMO can be performed by the same methods used for cell-free MIMO.

Fig. 3 shows the E-R regions of the colocated MIMO system for two different values of active eavesdropping power corresponding to training power splitting factors $\phi=0.3$ and $\phi=0.4$. It can be noticed that there is a tradoff between the ESR and the constraint on the AHE. As the AHE constraint increases, more downlink transmission resources are optimized to satisfy the increase in AHE constraint at the expense of the ESR which tends to decrease. Also, there is a clear gap between the ESR performances at different levels of active eavesdropping powers. The larger the eavesdropping power the lower the ESR.

Fig. 4 shows the E-R regions of the cell-free MIMO system for the same values of active eavesdropping powers used for colocated MIMO system. By comparing Fig. 3 and Fig. 4, it can be noticed that the cell-free MIMO outperforms the colocated MIMO within the interval in which the harvested energy constraint is low and vice versa. The cell-free MIMO is also found to be more immune to the increase in the active eavesdropping power than the colocated MIMO. In colocated MIMO, all antennas contribute to the AHE by an equal average value which is not the case for the cell-free MIMO. Therefore, the colocated MIMO is more efficient at power transfer than the cell-free MIMO. The difference between channel gains of the IU and the EH in the cell-free system offers the optimizer more freedom to balance the tradeoff between the information, AN and the energy signal powers than in the case of the colocated MIMO system. That justifies the advantage of cell-free MIMO over the colocated MIMO in the feasible region (the low AHE constraint region).

Fig. 5 shows how the density of APs affects the secrecy performance. The achievable worst-case ESR is measured versus a set of practically large values of AP densities $\lambda_{a}=10^{-4}\times[0.2,~{}0.4,~{}...,~{}2.6]~{}m^{-1}$. The values of the worst-case ESR in Fig. 5 are obtained by Monte Carlo simulation averaged over 50 independent realizations of AP deployments, with $\bar{E}=0$ and $\phi=0.3$. As expected, as the AP density (which is directly proportional to $\mathbb{E}[N]$) increases, the worst-case ESR increases.

The secrecy performance is affected by the relative location of the attacked IU with respect to the EH. Fig. 6 shows the ESR performance for the case where the system comprises one IU and one EH lying on the $x$-axis symmetrically around the origin of the APs’ deployment given in Fig 2. The results represent the achieved ESR for different separation distances between the IU and the EH, $\Delta=[0,~{}100,~{}\dots,500]$. As the separation $\Delta$ increases, the ESR performance improves. This can be justified since as the separation increases, the APs subsets that dominantly serve the IU and the EH become more distinctive. Beyond a certain value of $\Delta>200$, the achieved ESR starts to saturate since the dominant subsets of the APs that serve the IU and the EH remain unchanged, but the position of each user within its set varies. The value $\Delta=0$ means that the IU and the EH are colocated, i.e., $\boldsymbol{\Gamma}_{1}=\boldsymbol{\Gamma}$.

Fig. 7 shows the convergence speed of Algorithm 1 at $\phi=0.4$, $\bar{E}=2~{}mW$ and the initial values are selected arbitrarily as $\bar{s}_{i}^{[1]}=\bar{t}_{i}^{[1]}=[-8,~{}-6]~{}\forall~{}i$. As discussed in Subsection IV-C2, the optimized objective value is increasing with iterations until convergence.

VI Conclusions

In this paper, relaxed SDP programming has been proposed to optimize a nonlinear power control of the downlink transmission in a SWIPT cell-free massive MIMO system in the presence of an information-untrusted dual-antenna active EH. The downlink SWIPT transmissions include: information, AN and energy signals beamformed towards the IUs, legitimate and illegitimate antennas of the EH, respectively. Analytic expressions for the AHE and a tight lower bound on the ESR were derived with taking into account the IUs’ knowledge attained by downlink effective precoded-channel training. It has been proved that the proposed SDP iterative problem can always achieve a converged rank-one globally optimal solution. A fair comparison between the proposed cell-free and the colocated massive MIMO systems showed that the cell-free MIMO outperformed the colocated MIMO over the interval in which the AHE constraint is low and vice versa. Also, cell-free MIMO was more immune to the increase in the active eavesdropping power than the colocated MIMO.

Appendix A Proofs of Lemma 1, Theorems 1 and 2

A-A Proof of Lemma 1

Since the spectral radius of the diagonal matrices $\boldsymbol{\Gamma}_{i}$, $\boldsymbol{\Gamma}$ and $\sqrt{\boldsymbol{\Gamma}_{i}\boldsymbol{\Gamma}}$ are bounded [13, Lemma 2], then by expanding $\boldsymbol{y}_{i}^{H}\boldsymbol{y}_{i}$ followed by applying Corollary 1 in [39] we get

[TABLE]

which satisfies the asymptotic convergence in (3). This concludes the proof.

A-B Proof of Theorem 1

Before commencing our proof, let us introduce the following result.

Lemma 5

For a non-negative bounded RV $X_{1}\leq U$, $U$ is a positive real value, and a symmetrical zero mean RV $X_{2}$, the non-negative RV $Y=X_{1}+X_{2}$ is upper bounded as $Y\leq 2U$

Proof:

We have $P((X_{1}+X_{2})<0)=0$, then $P(U+X_{2}<0)=0$, i.e., $P(X_{2}<-U)=0$. By symmetry of distribution, $P(X_{2}>U)=0$ which implies $X_{2}\leq U$. Therefore, $Y\leq 2U$. This concludes the proof. ∎

Let $\text{IU}_{k}$ be the attacked IU. Based on (2), (5a), (7) and (9) we have

[TABLE]

where $\boldsymbol{\tilde{h}}_{k}=\boldsymbol{\hat{h}}_{k}-\tau\sqrt{P_{I}}\boldsymbol{C}_{k}\boldsymbol{h}_{k}=\tau\sqrt{\phi P_{E}}\;\boldsymbol{g}_{E}+\boldsymbol{N}\boldsymbol{\psi}_{k}^{*}$, and the entries of $\boldsymbol{\tilde{h}}_{k}$, $\boldsymbol{h}_{k}$ and $\boldsymbol{n}_{k}$ are statistically independent. Using Corollary 1 in [39] and Lemma 2 in [13], the first term in (LABEL:ApA21) asymptotically converges to the deterministic value

[TABLE]

where $\mathbb{I}=\{\{k,j\}_{j}\times\{k,m\}_{m}|\{k,j\}\neq\{k,m\}\}$ and $\bar{c}_{k,j}=[\boldsymbol{C}_{k}]_{j,j}$. With the assumption that the noise variance $\sigma_{n}^{2}\ll\tau\phi P_{E}$, we can approximate the sum of the second and the third terms in (LABEL:ApA21) as $\Delta_{2}=\tau^{2}\phi P_{E}|\boldsymbol{\bar{h}}_{k}^{H}\boldsymbol{\Gamma}_{k}^{\frac{1}{2}}\text{diag}\left(\boldsymbol{p}_{k}\right)\boldsymbol{C}_{k}\boldsymbol{g}_{E}|^{2}$, which is equivalent to

[TABLE]

where $\kappa_{j}=|\hat{h}_{k,j}\hat{g}_{E_{j}}|^{2}$ (equivalent to the product of two independent exponential RVs of parameter 1) is a non-negative RV with the mean value $\mathbb{E}[\kappa_{j}]=1$. $\kappa_{j,m}=\hat{h}_{k,j}\hat{g}_{E_{j}}\hat{h}_{k,m}\hat{g}_{E_{m}}$, $j\neq m$, is a zero mean RV with a symmetric distribution [40, 41]. Since $\Delta_{2}$ is always positive, i.e., $P(\Delta_{2}<0)=0$, then, by applying Theorem A in [42] (which defines an upper bound on the sum of non-negative RVs) and Lemma 5, $\Delta_{2}$ is upper bounded by a deterministic value as

[TABLE]

Given that the additive terms that constitute $\Delta_{1}$ in (LABEL:ApA22) and the upper bound of $\Delta_{2}$, $\overline{\Delta}_{2}$, in (36) are of a finite order of magnitude, then, asymptotically, we have $\Delta_{1}\overset{N\to\infty}{\to}\mathcal{O}\left(N^{2}\right)$ and $\overline{\Delta}_{2}\overset{N\to\infty}{\to}\mathcal{O}\left(N\right)$. Therefore, as $N\to\infty$, $\Delta_{1}$ and $\Delta_{2}$ differ by $\mathcal{O}\left(N\right)$ order of magnitude which implies that the bound $\left|\hat{a}_{k,k}\right|^{2}\geq\Delta_{1}=\tau^{2}P_{I}\text{tr}\left(\boldsymbol{P}_{k}\boldsymbol{A}_{k}\right)$ is tight. Based on this result, (14), (15), and since $\text{SINR}_{k}$ and $\underline{\text{SINR}}_{k}$ share the same denominator, then $\text{SINR}_{k}>\underline{\text{SINR}}_{k}$ in (15) is of the same degree of tightness. To validate the tightness of $\text{SINR}_{k}\overset{N\rightarrow\infty}{>}\underline{\text{SINR}}_{k}$ numerically, Table III presents the values of $\Delta_{1}$, $\overline{\Delta}_{2}$ and $\frac{\Delta_{1}}{\overline{\Delta}_{2}}$ for different realizations of $\{\boldsymbol{\Gamma}_{i}\}$ and $\boldsymbol{\Gamma}$ at a large average value of $N=100$, and at $\bar{E}=5~{}mW$. The optimized values of $\{\boldsymbol{P}_{i}\}$, $\boldsymbol{\bar{P}}$ and $\boldsymbol{P}$ used to obtain the values of $\Delta_{1}$ are used to calculate corresponding value of $\overline{\Delta}_{2}$. The obtained results validate our analysis.

The values of $\text{var}(Z_{k,k})=\sum_{j\neq k}c_{k,j}+\tau^{2}P_{I}\bar{c}_{k}^{2}+\bar{c}_{k}^{(1)}+\tilde{c}_{k}$ and $\text{var}(\tilde{a}_{k,k})=\sigma_{n}^{2}\frac{\tau_{d}+1}{\tau_{d}}$ (as in (14)-(15)) can be calculated as follows. We have

[TABLE]

The second equality follows from (12a). The third and fifth equalities follow from substituting (9) and (7), respectively. Since $\mathbb{E}(\boldsymbol{n}_{k})=0$, then $\mathbb{E}(\tilde{a}_{k,k})=0$ and therefore

[TABLE]

$\text{var}\left(\tilde{a}_{k,k}\right)+\sigma_{n}=\frac{\sigma_{n}}{\tau_{d}^{2}}+\sigma_{n}^{2}=\sigma_{n}^{2}\frac{\tau_{d}+1}{\tau_{d}}$ is equivalent to the last term in the denominator of (15). Regarding $Z_{k,k}$, the additive terms that constitute $Z_{k,k}$ in (11) are zero mean statistically independent RVs since the entries of $\{a_{k,j}\}_{j\neq k}$, $\boldsymbol{\bar{w}}_{k}$, $\boldsymbol{w}$ and $n_{i}$ are zero mean statistically independent RVs. Therefore,

[TABLE]

The expectations in (LABEL:Aprte1) are calculated as follows.

[TABLE]

The third equality in (LABEL:Aprffss1) is obtained by substituting the values of $\boldsymbol{h}_{k}$ and $\boldsymbol{w}_{j}$ from (5a). In the fourth equality, the expectation is moved to $\boldsymbol{\hat{h}}_{j}^{*}\boldsymbol{\hat{h}}_{j}^{T}$ based on the statistical independence between $\boldsymbol{\bar{h}}_{k}$ and $\boldsymbol{\hat{h}}_{j}$. The fifth equality follows since the entries of $\boldsymbol{\bar{h}}_{k}$ are zero mean unit variance independent RVs. The sixth and the seventh equalities follow since the matrices $\boldsymbol{\Gamma}_{k}$, $\text{diag}\left(\boldsymbol{p}_{j}\right)$ and $\boldsymbol{R}_{j}$ are diagonal. The forms $\boldsymbol{p}_{j}^{T}\boldsymbol{\Gamma}_{k}\boldsymbol{R}_{j}\boldsymbol{p}_{j}$ and $\text{tr}(\boldsymbol{P}_{j}\boldsymbol{\Gamma}_{k}\boldsymbol{R}_{j})$ are identical to those used in (15) and (26a)-(26b), respectively. For $\mathbb{E}[|\boldsymbol{h}_{k}^{T}\boldsymbol{\bar{w}}_{k}|^{2}]$, we substitute the value of $\boldsymbol{\bar{w}}_{k}$ from (5b), (2a) and (2b) as

[TABLE]

where $\boldsymbol{\tilde{h}}_{k}=\delta_{ik}\;\tau\sqrt{\phi P_{E}}\;\boldsymbol{g}_{E}+\boldsymbol{N}\boldsymbol{\psi}_{i}^{*}$. The second equality follows from the statistical independent between $\boldsymbol{h}_{k}$ and $\boldsymbol{\tilde{h}}_{k}$. The first term after the third equality follows from applying Corollary 1 in [39]. In the second term in the third equality, the expectation is moved to $\boldsymbol{\tilde{h}}_{k}^{*}\boldsymbol{\tilde{h}}_{k}^{T}$ based on the statistical independence between $\boldsymbol{\bar{h}}_{k}$ and $\boldsymbol{\tilde{h}}_{k}$. The first term in the fourth equality follows since the matrices $\boldsymbol{\Gamma}_{k}$, $\boldsymbol{C}_{k}$ and $\text{diag}\left(\boldsymbol{\bar{p}}\right)$ are diagonal. The second term in the fourth equality follows since the entries of $\boldsymbol{\bar{h}}_{k}$ are zero mean unit variance independent RVs and the matrices $\boldsymbol{\Gamma}_{k}$ and $\boldsymbol{R}_{k}^{(1)}$ are diagonal. The form which is after the fourth equality is identical to that used in (15), while the SDP form which is after the sixth equality is identical to that used in (26a)-(26b). For $\mathbb{E}[|\boldsymbol{h}_{k}^{T}\boldsymbol{w}|^{2}]$, we substitute the value of $\boldsymbol{w}$ from (5c) as

[TABLE]

In the second equality, the expectation is moved to $\boldsymbol{\hat{g}}^{*}\boldsymbol{\hat{g}}^{T}$ based on the statistical independence between $\boldsymbol{\bar{h}}_{k}$ and $\boldsymbol{\hat{g}}$. The third equality follows since the entries of $\boldsymbol{\bar{h}}_{k}$ are zero mean unit variance independent RVs and the matrices $\boldsymbol{\Gamma}_{k}$ and $\boldsymbol{R}$ are diagonal. The form which is after the third equality is identical to that used in (15), while the SDP form which is after the fifth equality is identical to that used in (26a)-(26b).

A-C Proof of Theorem 2

Based on the assumption that the EH has a full knowledge of the IUs’ beamforming vectors and its own channel, the EH is capable of cancelling the inter-user interference [30]. Furthermore, the information, AN and energy signals; $\{q_{i}\}$, $z$ and $\boldsymbol{w}$; are statistically independent. Therefore, based on the concavity of the logarithmic function, applying Jensen’s inequality (which has been proven to be tight and suitable for characterizing the performance of massive MIMO systems [43]) will result in the following upper bound on the ergodic rate at the EH

[TABLE]

where $\text{SINR}_{E_{k}}$ is the SINR at the EH when attacking $\text{IU}_{k}$. As defined in (LABEL:eeq20), ${\text{SINR}_{E_{k}}}$ ${=\frac{X_{k}}{Y_{k}}}$, ${X_{k}=\mathbb{E}[|\boldsymbol{g}_{E}^{T}\boldsymbol{w}_{k}q_{k}|^{2}]=|b_{k}|^{2}}$ and ${Y_{k}}$ ${=\mathbb{E}[|\boldsymbol{g}_{E}^{T}(\boldsymbol{\bar{w}}_{k}z+\boldsymbol{w})+\bar{n}|^{2}]}$ ${=|\hat{b}_{k}|^{2}+|b|^{2}+\sigma_{n}^{2}}$ . Using the multivariate Taylor expansion, ${\mathbb{E}[\text{SINR}_{E_{k}}]}$ can be expanded as [44]

[TABLE]

where $R=f(\text{var}(Y_{k}),\text{cov}(X_{k},Y_{k}))$ is the reminder of the series expansion. We have

[TABLE]

where

[TABLE]

[TABLE]

where $\boldsymbol{\tilde{h}}_{k}^{(2)}=\boldsymbol{\hat{h}}_{k}^{*}-\tau\sqrt{\phi P_{E}}\boldsymbol{C}_{k}\boldsymbol{g}_{E}^{*}$ and $\boldsymbol{R}_{k}^{(2)}=\mathbb{E}[\boldsymbol{\tilde{h}}_{k}^{(2)}\boldsymbol{\tilde{h}}_{k}^{(2)^{H}}]=\boldsymbol{\bar{R}}_{k}-\tau^{2}\phi P_{E}\boldsymbol{C}_{k}^{2}\boldsymbol{\Gamma}$. The third equality in (LABEL:ApAc3) is obtained by substituting the value of $\boldsymbol{w}_{k}$ from (5a), (2a) and (2b). The fourth equality in (LABEL:ApAc3) follows from the statistical independence between $\boldsymbol{g}_{E}$ and $\boldsymbol{\tilde{h}}_{k}^{(2)}$. The first equality in (46) follows from applying Corollary 1 in [39] and the diagonality of the matrices $\boldsymbol{\Gamma}$, $\boldsymbol{C}_{k}$ and $\text{diag}\left(\boldsymbol{p}_{k}\right)$. In the first equality in (47), the expectation is moved to $\boldsymbol{\tilde{h}}_{k}^{(2)}\boldsymbol{\tilde{h}}_{k}^{(2)^{H}}$ based on the statistical independence between $\boldsymbol{\bar{g}}_{E}$ and $\boldsymbol{\tilde{h}}_{k}^{(2)}$. The second equality in (47) follows from applying Corollary 1 in [39]. The third and the fourth equalities in in (47) follow from the diagonality of the matrices $\boldsymbol{\Gamma}$, $\boldsymbol{R}_{k}^{(2)}$ and $\text{diag}\left(\boldsymbol{p}_{k}\right)$. By substituting (46) and (47) in (LABEL:ApAc3), we get $\mathbb{E}[X_{k}]=\mathbb{E}[|b_{k}|^{2}]=\tau^{2}\phi P_{E}d_{k}^{2}+d_{k}^{(1)}$, $d_{k}^{2}=\text{tr}(\boldsymbol{P}_{k}\boldsymbol{B}_{k})$ and $d_{k}^{(1)}=\text{tr}(\boldsymbol{P}_{k}\boldsymbol{\bar{B}}_{k})$.

For $\mathbb{E}[Y_{k}]$, based on the statistical independence between $\boldsymbol{\bar{w}}_{k}$ and $\boldsymbol{w}$, we have

[TABLE]

[TABLE]

The first equality in (LABEL:Apxxc05) follows from substituting the value of $\boldsymbol{\bar{w}}_{k}$ from (5b), (2a) and (2b); and the statistical independence between $\boldsymbol{g}_{E}$ and $\boldsymbol{\tilde{h}}_{k}^{(2)}$. The first term in the second equality follows from applying Corollary 1 in [39]. In the second term after the second equality, the expectation is moved to $\boldsymbol{\tilde{h}}_{k}^{(2)}\boldsymbol{\tilde{h}}_{k}^{(2)^{H}}$ based on the statistical independence between $\boldsymbol{\bar{g}}_{E}$ and $\boldsymbol{\tilde{h}}_{k}^{(2)}$. The first term in the third equality follows from the diagonality of the matrices $\boldsymbol{\Gamma}$, $\boldsymbol{C}_{k}$ and $\text{diag}\left(\boldsymbol{\bar{p}}_{k}\right)$, while the second term follows since the entries of $\boldsymbol{\bar{g}}_{E}$ are zero mean unit variance independent RVs and the matrices $\boldsymbol{\Gamma}$ and $\boldsymbol{R}_{k}^{(2)}$ are diagonal. The fourth equality follows since the matrices $\boldsymbol{\Gamma}$, $\boldsymbol{C}_{k}$ and $\boldsymbol{R}_{k}^{(2)}$ are diagonal.

[TABLE]

The first equality in (LABEL:Apjxc05), the expectation is moved to $\boldsymbol{\hat{g}}^{*}\boldsymbol{\hat{g}}^{T}$ based on the statistical independence between $\boldsymbol{\bar{g}}_{E}$ and $\boldsymbol{\hat{g}}$. The second equality follows since the entries of $\boldsymbol{\bar{g}}_{E}$ are zero mean unit variance independent RVs. The third equality follows since the matrices $\boldsymbol{\Gamma}$ and $\boldsymbol{R}$ are diagonal. By substituting the results from (LABEL:Apxxc05) and (LABEL:Apjxc05) in (48) we get

[TABLE]

where $d=\text{tr}(\boldsymbol{P}\boldsymbol{B}),~{}\bar{d}_{k}^{2}=\text{tr}(\boldsymbol{\bar{P}}\boldsymbol{B}_{k})$ and $\bar{d}_{k}^{(1)}=\text{tr}(\boldsymbol{\bar{P}}\boldsymbol{\bar{B}}_{k})$. For $\text{var}(Y_{k})$ we have

[TABLE]

Now let us calculate the first and the second terms in (LABEL:ApAddc5), as follows

[TABLE]

[TABLE]

where $\circ$ denotes the Hadamard power operation and $\mathbb{I}=\{\{j\}\times\{m\}|~{}j\neq m\}$. The second equality in (LABEL:ApAddc5) follows from the statistical independence between $\boldsymbol{w}$, $\boldsymbol{\bar{w}}_{k}$ and $\bar{n}$. The second equality in (LABEL:ApAc6a) is obtaioned by substituting the value of $\mathbb{E}[|\boldsymbol{g}_{E}^{T}\boldsymbol{w}|^{2}]$ from (LABEL:Apjxc05). The second equality in (54) is obtained since the norm $\left|\boldsymbol{g}_{E}^{T}\boldsymbol{\bar{w}}_{k}\right|^{2}$ converges to $\tau^{2}\phi P_{E}\text{tr}(\boldsymbol{\bar{P}}\boldsymbol{B}_{k})+|\boldsymbol{g}_{E}^{T}\text{diag}(\boldsymbol{\bar{p}})\boldsymbol{\tilde{h}}_{k}^{(2)}|^{2}$ in the signal domain. By expanding the norms $|\boldsymbol{g}_{E}^{T}\boldsymbol{w}|^{4}$ and $|\boldsymbol{g}_{E}^{T}\text{diag}(\boldsymbol{\bar{p}})\boldsymbol{\tilde{h}}_{k}^{(2)}|^{4}$ – which are composed of squared exponential RVs – followed by applying the statistical expectation999The expectation is obtained by applying the fact: For an exponential RV $X\sim\mathcal{E}(\lambda)$, the $n$th moment of $X$ is $\mathbb{E}[X^{n}]=\frac{\Gamma(n+1)}{\lambda^{n}}=\frac{\Gamma(n+1)}{\lambda^{n}}=\frac{n!}{\lambda^{n}}$., we obtain the final results in (LABEL:ApAc6a) and (54), respectively.

Given that the entries of the matrix $\boldsymbol{B}_{k}$ have non-zero positive values, the matrices $\boldsymbol{\bar{B}}_{k}$ and $\boldsymbol{B}$ are diagonal, and the additive terms in $(\text{tr}(\boldsymbol{\bar{P}}\boldsymbol{B}_{k})+\text{tr}(\boldsymbol{\bar{P}}\boldsymbol{\bar{B}}_{k})+\text{tr}(\boldsymbol{P}\boldsymbol{B}))^{2}$ are of a finite order of magnitude; then, asymptotically, we have $(\mathbb{E}[Y_{k}])^{2}\overset{N\to\infty}{\to}\mathcal{O}\left(N^{4}\right)$ and $\text{var}(Y_{k})\overset{N\to\infty}{\to}\mathcal{O}\left(2N^{2}\right)$. This implies that $\frac{\text{var}(Y_{k})}{(E[Y_{k}])^{2}}\overset{N\to\infty}{\to}0$ and $R=f(\text{var}(Y_{k}),\text{cov}(X_{k},Y_{k}))\overset{N\to\infty}{\to}0$. Based on (LABEL:ApAc2) and this result, and since $\text{cov}(X_{k},Y_{k})=0$ (follows from the statistical independence between $X_{k}$ and $Y_{k}$), we have

[TABLE]

By substituting (LABEL:ApAc7) in (43), we get (18)–(LABEL:eeq20). This concludes the proof.

A-D Deriving the asymptotic value of AHE in (LABEL:eeq22)

The details of deriving $\mathbb{E}[|b_{k}|^{2}]$ and $\mathbb{E}[|\hat{b}_{k}|^{2}+|b|^{2}]$ are provide (LABEL:ApAc3)-(47) and (48)-(LABEL:ApAcxs05), respectively. The details of deriving the values $\mathbb{E}[|b_{j\neq k}|^{2}]$, $\mathbb{E}[|\tilde{b}_{j}|^{2}]$, $\mathbb{E}[|\tilde{\hat{b}}_{j}|^{2}]$ and $\mathbb{E}[|\tilde{b}|^{2}]$ that constitute $\bar{E}_{k}$ in (LABEL:eeq22) are as follows

[TABLE]

[TABLE]

[TABLE]

In the second equalities in (LABEL:Apjfgc05), (LABEL:Apjxchj5) and (LABEL:Ap89chj5), the expectation is moved to $\boldsymbol{\hat{h}}_{j\neq k}^{*}\boldsymbol{\hat{h}}_{j\neq k}^{T}$, $\boldsymbol{\hat{h}}_{j}^{*}\boldsymbol{\hat{h}}_{j}^{T}$ and $\boldsymbol{\hat{h}}_{k}^{*}\boldsymbol{\hat{h}}_{k}^{T}$ based on the statistical independence between $\boldsymbol{\bar{g}}_{E}$ and $\boldsymbol{\hat{h}}_{j\neq k}$, $\boldsymbol{\bar{g}}$ and $\boldsymbol{\hat{h}}_{j}$, and between $\boldsymbol{\bar{g}}_{E}$ and $\boldsymbol{\hat{h}}_{k}$, respectively. The third equalities in (LABEL:Apjfgc05), (LABEL:Apjxchj5) and (LABEL:Ap89chj5) follows since the entries of $\boldsymbol{\bar{g}}_{E}$ and $\boldsymbol{\bar{g}}$ are zero mean unit variance independent RVs. The fourth equalities follows since the matrices $\boldsymbol{\Gamma}$, $\boldsymbol{R}_{j}$ and $\boldsymbol{\bar{R}}_{k}$ are diagonal.

[TABLE]

The second equality follows from the statistical independent between $\boldsymbol{g}$ and $\boldsymbol{N}$. The first term after the fourth equality follows from applying Corollary 1 in [39]. In the second term in the fourth equality, the expectation is moved to $\boldsymbol{N}^{*}\psi_{E}\psi_{E}^{H}\boldsymbol{N}^{T}$ based on the statistical independence between $\boldsymbol{\bar{g}}$ and $\boldsymbol{N}$. The first term in the fifth equality follows since the matrices $\boldsymbol{\Gamma}$, $\boldsymbol{C}$ and $\text{diag}\left(\boldsymbol{p}\right)$ are diagonal. The second term in the fifth equality follows since $\mathbb{E}[\boldsymbol{N}^{*}\psi_{E}\psi_{E}^{H}\boldsymbol{N}^{T}]=\tau\sigma_{n}^{2}\boldsymbol{I}_{N}$ and the entries of $\boldsymbol{\bar{g}}$ are zero mean unit variance independent RVs. The seventh equality follows since the matrices $\boldsymbol{\Gamma}$ and $\boldsymbol{C}$ are diagonal. The form which is after the fifth equality is identical to that used in (LABEL:eeq22), while the SDP form which is after the seventh equality is identical to that used in (26e).

Appendix B Proof of Theorem 3

To prove that the optimal solution $\{\boldsymbol{P}_{i}^{\star}\},~{}\boldsymbol{\bar{P}}^{\star},~{}\boldsymbol{P^{\star}}$ obtained by solving (26) is always of unity rank, we exploit the boundedness property of the dual Lagrangian function to show that the optimal primal matrices $\{\boldsymbol{P}_{i}^{\star}\},~{}\boldsymbol{\bar{P}}^{\star},~{}\boldsymbol{P^{\star}}$ can satisfy the KKT conditions of optimality at one case in which $\{\text{rank}(\boldsymbol{P}_{i}^{\star})\},~{}\text{rank}(\boldsymbol{\bar{P}}^{\star}),~{}\text{rank}(\boldsymbol{P}^{\star})=1$, and that has been validated by computer simulation.

The Lagrangian of the equivalent problem (27) is

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\mathbb{L}=\{\{\lambda_{1_{k}}\},~{}\dots,~{}\{\lambda_{7_{k}}\},~{}\{\boldsymbol{F}_{k}\},~{}\boldsymbol{\bar{F}},~{}\boldsymbol{F}\}$ are the Lagrange multipliers of the constraints (27a), (26a)–(26f) and the constraints on $\{\boldsymbol{P}_{k}\}$, $\boldsymbol{\bar{P}}$ and $\boldsymbol{P}$ in (26g), respectively, with $\{\lambda_{j_{k}}\}\geq 0$, $\{\boldsymbol{F}_{k}\},~{}\boldsymbol{\bar{F}},~{}\boldsymbol{F}\succeq\boldsymbol{0}$. Now, for the Lagrangian function to exist, the infimum of $\mathcal{L}$ over the primal variable $\mathbb{S}$, $\inf_{\mathbb{S}}\mathcal{L}$, should be bounded from below, therefore, we have

[TABLE]

Given that Slater’s condition holds (see (28) and the paragraph that follows) and based on the non-negativeness of the dual variables (Lagrange multipliers), the satisfaction of the KKT’s complementary slackness condition results in $\text{tr}(\boldsymbol{P}_{k}^{\star}\boldsymbol{F}_{k}^{\star})=0~{}\forall~{}k$. The KKT’s stationarity condition should satisfy $\sum_{k}\frac{\partial\mathcal{L}}{\partial\boldsymbol{P}_{k}^{\star}}=\boldsymbol{0}$, therefore

[TABLE]

Given the fact that if the summation of multiple positive semi-definite matrices is equal to zero, all the matrices are equal to zero. And based on the (B), (65) and (66), we have

[TABLE]

[TABLE]

Let $\text{null}\left(\boldsymbol{F}_{k}^{\star}\right)=\boldsymbol{\Omega}_{k}=\left[\boldsymbol{\omega}_{k,1},...,\boldsymbol{\omega}_{k,N-\text{rank}(\boldsymbol{F}_{k}^{\star})}\right]$ and $\text{null}\left(\boldsymbol{G}_{k}^{\star}\right)=\boldsymbol{\Psi}_{k}=\left[\boldsymbol{\psi}_{k,1},...,\boldsymbol{\psi}_{k,N-\text{rank}(\boldsymbol{G}_{k}^{\star})}\right]$. By making use of the inequality of matrix sum [45, subsection 3.3.4] and (67) we have101010Please note that $\text{rank}(\boldsymbol{A}_{k})=1$, this is understandable from the structure of $\boldsymbol{A}_{k}$. Please refer to the first paragraph in IV-B. $\text{rank}(\boldsymbol{F}_{k}^{\star})\geq\text{rank}(\boldsymbol{G}_{k}^{\star})-1$. Based on this result, and since $\text{rank}(\boldsymbol{\Omega}_{k})=N-\text{rank}(\boldsymbol{F}_{k}^{\star})$ and $\text{rank}(\boldsymbol{\Psi}_{k})=N-\text{rank}(\boldsymbol{G}_{k}^{\star})$, then the following result is true

[TABLE]

Now, let us examine the null space of $\boldsymbol{G}_{k}^{\star}$, $\boldsymbol{\psi}_{k,j}\in\boldsymbol{\Psi}_{k}$, by computing the inner product between $\boldsymbol{\psi}_{k,j}$ and $\boldsymbol{F}_{k}^{\star}$ in (67) as follows

[TABLE]

where the inequality in (70) follows from111111$\boldsymbol{A}_{k}\succeq 0$ follows since it is structured from a vector whose all entries are positive (see (25) and the paragraph that follows). $\boldsymbol{A}_{k}\succeq 0$. However, since $\boldsymbol{F}_{k}^{\star}\succeq 0$, (70) can only hold with equality, i.e., $\boldsymbol{\psi}_{k,j}^{H}\left(\lambda_{2_{k}}^{\star}\tau^{2}P_{I}\boldsymbol{A}_{k}\right)\boldsymbol{\psi}_{k,j}=0$. This result implies that the null space of $\boldsymbol{G}_{k}^{\star}$ always forms null space of $\boldsymbol{F}_{k}^{\star}$, i.e., $\boldsymbol{\Psi}_{k}$ is a sub-matrix of $\boldsymbol{\Omega}_{k}$, therefore, and according to (69), $\boldsymbol{\omega}_{k,j}\in\boldsymbol{\Omega}_{k}$ belongs to one of the following two spaces: 1) the column space of $\boldsymbol{\Psi}_{k}$, $\boldsymbol{\omega}_{k,j}\in\{\boldsymbol{\psi}_{k,j}\}$; 2) 1-dimensional vector space, $\boldsymbol{\omega}_{k,j}=\boldsymbol{a}\in\mathcal{C}^{N\times 1}$ where $\boldsymbol{a}\notin\{\boldsymbol{\psi}_{k,j}\}$.

Since the optimal value of $\boldsymbol{P}_{k}^{\star}$ needs to satisfy the complementary slackness condition, $\text{tr}(\boldsymbol{P}_{k}^{\star}\boldsymbol{F}_{k}^{\star})=0~{}\forall~{}k$, the structure of $\boldsymbol{P}_{k}^{\star}$ is

[TABLE]

where $\{m_{k,j}\}$ are non-negative scaling factors. The $\boldsymbol{P}_{k}^{\star}$’s component $m_{k,j}\boldsymbol{\psi}_{k,j}\boldsymbol{\psi}_{k,j}^{H}$ introduces zero information signal power at $\text{IU}_{k}$ since $\boldsymbol{\psi}_{k,j}^{H}\boldsymbol{A}_{k}\boldsymbol{\psi}_{k,j}=0$, and therefore contributes by a negative ESR. Thus, $m_{k,j}\boldsymbol{\psi}_{k,j}\boldsymbol{\psi}_{k,j}^{H}$ is a non-optimal component of $\boldsymbol{P}_{k}^{\star}$. By this, we can conclude that $\boldsymbol{P}_{k}^{\star}$ is constructed by the single component $\boldsymbol{P}_{k}^{\star}=m_{k,1}\boldsymbol{a}\boldsymbol{a}^{H}$, $\boldsymbol{a}\notin\{\boldsymbol{\psi}_{k,j}\}$, therefore, $\text{rank}(\boldsymbol{P}_{k}^{\star})=1$ is always true. This concludes the proof.

Appendix C Proofs of Lemmas 3 and 4

C-A Proof of Lemma 3

Since the exponential function is convex (has a downward curvature), the tangent line at any point is below the function trajectory. Using triangulation (as depicted in Fig. 8), it can be easily understood that the value of the Taylor first order approximation of $e^{x}$, $e^{\bar{x}}(x-\bar{x}+1)$, always lies at the tangent line ($L_{1}$, black solid line) which is always below the function trajectory (circle-marked line). Fig. 8 shows the case $\bar{x}<x$. Following a comparable reasoning, the previous result can be proved for the other case $\bar{x}>x$. This concludes the proof.

C-B Proof of Lemma 4

Building upon the proof of Lemma 3, the value of successive Taylor approximation $e^{\bar{\bar{x}}}(x-\bar{\bar{x}}+1)$, $e^{\bar{\bar{x}}}=e^{\bar{x}}(x-\bar{x}+1)$, lies at the tangent line ($L_{2}$, star-marked line in Fig. 8) touched at $(\bar{\bar{x}},e^{\bar{\bar{x}}})$. Since the derivative of the exponential function is non-decreasing, $L_{2}$ always lies above $L_{1}$ for $x\geq\bar{\bar{x}}$. Therefore, $e^{\bar{\bar{x}}}(x-\bar{\bar{x}}+1)>e^{\bar{x}}(x-\bar{x}+1)$ is always true. This concludes the proof.