# Cell-Free Massive MIMO SWIPT with Beyond Diagonal Reconfigurable Intelligent
Surfaces

Thien Duc Hua, Mohammadali Mohammadi, Hien Quoc Ngo, and Michail Matthaiou
Centre for Wireless Innovation (CWI), Queen’s University Belfast, U.K.
Email:{dhua01, m.mohammadi, hien.ngo<EMAIL_ADDRESS>
###### Abstract

This paper investigates the integration of beyond-diagonal reconfigurable
intelligent surfaces (BD-RISs) into cell-free massive multiple-input multiple-
output (CF-mMIMO) systems, focusing on applications involving simultaneous
wireless information and power transfer (SWIPT). The system supports
concurrently two user groups: information users (IUs) and energy users (EUs).
A BD-RIS is employed to enhance the wireless power transfer (WPT) directed
towards the EUs. To comprehensively evaluate the system’s performance, we
present an analytical framework for the spectral efficiency (SE) of IUs and
the average harvested energy (HE) of EUs in the presence of spatial
correlation among the BD-RIS elements and for a non-linear energy harvesting
circuit. Our findings offer important insights into the transformative
potential of BD-RIS, setting the stage for the development of more efficient
and effective SWIPT networks. Finally, incorporating a heuristic scattering
matrix design at the BD-RIS results in a substantial improvement compared to
the scenario with random scattering matrix design. ††footnotetext: This work
was supported by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement No.
101001331). The work of H. Q. Ngo was supported by the U.K. Research and
Innovation Future Leaders Fellowships under Grant MR/X010635/1.

##  I Introduction

The telecommunications sector is undergoing rapid expansion, fueled by the
increasing demand for high data traffic. This surge is notably driven by the
deployment of 5G and 6G networks and the widespread adoption of the Internet
of Everything (IoE). As a consequence, there is an anticipated substantial
rise in the energy consumption of wireless networks. Recognizing this
challenge, SWIPT has emerged as a promising technology to address the growing
energy demands in this dynamic landscape [1]. Traditional cellular systems
exhibit inherent limitations, including elevated infrastructure costs, inter-
cell interference, cell-edge effects, and significant path loss attributed to
extended distances. These constraints impede the seamless application of
SWIPT. To tackle these challenges, CF-mMIMO emerges as a promising technology.
In this paradigm, several APs are distributed across the coverage area,
effectively serving multiple users. This approach reduces the distance between
users and their nearby APs, fostering macro-diversity and mitigating path loss
[2]. As a result, CF-mMIMO provides a compelling solution to improve the SE of
wireless networks and enhance the performance of WPT in SWIPT scenarios [3, 4,
5].

In recent advancements, RISs have been seamlessly integrated into wireless
networks, enhancing the overall system performance. Notably, a RIS offers the
capability to shape radio waves at the electromagnetic level, eliminating the
need for digital signal processing methods and power amplifiers [6]. Recent
research efforts have delved into the exploration of RIS-assisted CF-mMIMO
networks [7, 8]. Specifically, Trinh et al. [7] conducted an analysis of the
uplink and downlink SE in RIS-assisted CF-mMIMO systems. Additionally, Dai et
al. [8] investigated and optimized the achievable rate performance of uplink
RIS-aided CF-mMIMO systems.

Although the application of RISs has undeniably enhanced the WPT efficiency in
various scenarios [9], their potential within the context of CF-mMIMO SWIPT
systems has remained vastly unexplored. In their recent work, Shi et al. [10]
investigated potential application scenarios and system architectures for RIS-
aided CF-mMIMO systems in the context of WPT. More precisely, they discussed
the utilization of RISs in these systems to direct energy beams towards
obstructed energy zones. Recently, the BD-RIS concept, which unifies different
RIS modes/architectures, has garnered significant attention from the research
community [11, 12]. This motivates us to investigate, for the first time ever,
the performance of BD-RIS-assisted CF-mMIMO SWIPT systems.

We focus on a CF-mMIMO SWIPT system that serves two distinct groups: IUs and
EUs. To optimize the resource utilization, we categorize the available APS
into information APs (I-APs) and energy APs (E-APs). This ensures that the APs
simultaneously serve IUs and EUs throughout the entire time slot, resulting in
enhanced SE and EH. We employ a BD-RIS to aid in directing energy beams toward
designated energy zones, thereby efficiently supporting the EUs. This novel
architecture allows the EUs to harvest energy from all APs, but it also
creates more interference at the IUs due to simultaneous WPT. To tackle this,
we utilize the protective partial zero-forcing (PPZF) precoder, which applies
partial zero-forcing (ZF) precoding at I-APs and protective maximum ratio
transmission (PMRT) at E-APS. PMRT guarantees full protection of the IUs from
the energy-transmitted interference. We consider a more realistic channel
model by taking into account the spatial correlation among the BD-RIS
elements. Our major contributions are:

  * •

We derive the second-order and fourth-order moments of the estimate of the AP-
to-UE channel, which consists of direct and indirect links. The indirect link
takes into account the spatial correlation among the channels of the BD-RIS
elements.

  * •

We provide closed-form expressions for the ergodic SE of the IUs and the HE of
the EUs with non-linear energy harvester. These expressions are derived based
on the channel moments and the precoding scheme. The obtained results
establish a fundamental basis for designing and optimizing the performance of
the considered CF-mMIMO SWIPT system.

  * •

We conduct extensive simulations to compare the performance of a heuristic
scattering matrix design for BD-RIS against a random design, and no service
from the BD-RIS.

Notation: We use bold upper/lower case letters to denote matrices/vectors. The
superscripts $(\cdot)^{T}$ and $(\cdot)^{H}$ stand for the transpose and the
conjugate-transpose, respectively; $\mathbf{I}_{N}$ denotes the $N\times N$
identity matrix. A circular symmetric complex Gaussian variable with variance
$\sigma^{2}$ is denoted by $\mathcal{CN}(0,\sigma^{2})$. Finally,
$\mathbb{E}\\{\cdot\\}$ denotes the statistical expectation.

##  II System model

Figure 1 illustrates a CF-mMIMO SWIPT system aided by a BD-RIS. In this
configuration, a central processing unit (CPU) connects to $M$ APs. All APs
cooperate to simultaneously serve $K$ IUs and $J$ EUs in the same frequency
bands. We define the sets $\mathcal{K}\triangleq\\{1,\ldots,K\\}$,
$\mathcal{J}\triangleq\\{1,\ldots,J\\}$ and
$\operatorname{\mathcal{M}}\triangleq\\{1,\ldots,M\\}$ as the collections of
indices of the IUs, EUs, and APs, respectively. The APs are equipped with $L$
antennas each, while the IUs and EUs are equipped with a single antenna, so
that $ML\gg(K+J)$. To ensure SWIPT under the same frequency spectrum, we
consider the selection of operation modes for APs. Therefore, the APs are
classified into I-APs and E-APs. The I-APs support wireless information
transfer (WIT) to the IUs, while the E-APs support WPT towards the EUs. To
further improve the HE of the EUs, a BD-RIS is located near the EU zone to
support the transmission. In the BD-RIS, the $n$-th element connects to all
$n^{\prime}$-th elements for $n^{\prime}\neq n,\forall
n,n^{\prime}\in\mathcal{N}$ through a network of impedances, where
$\mathcal{N}=\\{1,\ldots,N\\}$ denotes the set of BD-RIS elements. Due to
element-wise correlation, the scattering matrix
$\operatorname{\boldsymbol{\Theta}}$ is a full matrix, and satisfies
$\operatorname{\boldsymbol{\Theta}}=\operatorname{\boldsymbol{\Theta}}^{T}$,
and
$\operatorname{\boldsymbol{\Theta}}^{H}\operatorname{\boldsymbol{\Theta}}=\textbf{I}_{N}$.
In addition, the condition of perfect matching and no mutual coupling can be
practically achieved by individually matching each antenna to the reference
impedance and maintaining an antenna spacing greater than half-wavelength [11,
12].

###  II-A Channel Model

The considered architecture operates in time division duplex, which implies
that channel reciprocity holds as in the canonical form of CF-mMIMO [2]. We
assume a block-fading channel model, i.e., the channel remains time-invariant
and frequency-flat during each coherence interval containing $\tau_{c}$
symbols, and independently varies between coherence intervals. Particularly,
we denote $\tau$ as the number of symbols per coherence interval spent on
transmission of uplink training, and $(\tau_{c}-\tau)$ is the duration of
downlink transmission.

The channel vector between the $m$-th I-AP and the $k$-th IU is denoted by
${\textbf{g}}_{mk}^{\mathtt{I}}\in\mathbb{C}^{L\times
1},\forall{m}\in\operatorname{\mathcal{M}}$, and
$\forall{k}\in\operatorname{\mathcal{K}}$. Since the BD-RIS is deployed to
improve the EH system, it is located near the EU area which is far from the IU
area. As a result, the indirect links reflected from the BD-RIS to the IUs are
relative weak compared to the direct links, and are neglected in
${\textbf{g}}_{mk}^{\mathtt{I}}$. The aggregated channel between the $m$-th
E-AP and the $j$-th EU, ${\textbf{g}}_{mj}^{\mathtt{E}}\in\mathbb{C}^{L\times
1},\forall{m}\in\operatorname{\mathcal{M}}$, and $\forall{j}\in\mathcal{J}$,
is comprising the direct link, ${\textbf{h}}_{mj}^{\mathtt{E}}$, and cascaded
indirect links reflected from the BD-RIS,
${\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}{\bf z}_{j}$.
Mathematically speaking,

$\leavevmode\nobreak\
{\textbf{g}}_{mj}^{\mathtt{E}}={\textbf{h}}_{mj}^{\mathtt{E}}+{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}{\bf
z}_{j},$ (1)

where ${\textbf{h}}_{mj}^{\mathtt{E}}\in\operatorname{\mathbb{C}}^{L\times
1}$, ${\textbf{H}}_{m}\in\operatorname{\mathbb{C}}^{N\times L}$, and ${\bf
z}_{l}\in\operatorname{\mathbb{C}}^{N\times 1}$ are the channels from the
$m$-th AP to the j-th EU, from the $m$-th AP to the RIS, and from the RIS to
the j-th EU, respectively. In particular, we denote the $l$-th column vector
of ${\textbf{H}}_{m}$ as an $N$-dimensional vector ${\textbf{h}}_{ml}$, such
that ${\textbf{H}}_{m}=[{\textbf{h}}_{m1},\ldots,{\textbf{h}}_{mL}]$ with
${\textbf{h}}_{ml}\in\operatorname{\mathbb{C}}^{N\times 1}$. We assume that
${\textbf{g}}_{mk}^{\mathtt{I}}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\beta_{mk}^{\mathtt{I}}{\textbf{I}}_{L})$
and
${\textbf{h}}_{mj}^{\mathtt{E}}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\beta_{mj}^{\mathtt{E}}{\textbf{I}}_{L})$.
Moreover, as in [7], we consider the spatial correlated Rayleigh fading model
for the BD-RIS-related channels as follows:

$\displaystyle{\textbf{h}}_{ml}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\boldsymbol{\Omega}_{m}),\hskip
10.00002pt\text{and}\hskip 10.00002pt{\bf
z}_{j}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\boldsymbol{\Psi}_{j}),$

where $\boldsymbol{\Omega}_{m}=\alpha_{m}d_{H}d_{V}{\bf R}$ and
$\boldsymbol{\Psi}_{j}=\alpha_{j}d_{H}d_{V}{\bf R}$, while $\beta^{I}_{mk}$,
$\beta_{mj}^{\mathtt{E}}$, $\alpha_{m}$, and $\alpha_{j}$ are the large-scale
fading coefficients; $d_{H}$ and $d_{v}$ are the horizontal and vertical
lengths of a RIS element. Moreover, ${\bf
R}\in\operatorname{\mathbb{C}}^{N\times N}$ is the normalized spatial
correlation among the RIS elements. Each element of the normalized spatial
correlation matrix ${\bf R}$, can be expressed as [7]

$\leavevmode\nobreak\ [{\bf R}]_{nn^{\prime}}=\text{sinc}\bigg{(}\frac{2||{\bf
u}_{n}-{\bf
u}_{n^{\prime}}||}{\lambda}\bigg{)},\forall{n,n^{\prime}}\in\mathcal{N},$ (2)

where

$\leavevmode\nobreak\ {\bf
u}_{n}\\!=\\![0,\text{mod}(n\\!-\\!1)N_{H},\\!\lfloor(n-1)/N_{H}\rfloor
N_{V}]^{T},\forall{n}\in\mathcal{N},$ (3)

is the location vector of the $n_{I}$-th element w.r.t. the origin. In (3),
$N_{H}$ and $N_{V}$ are the respective numbers of RIS elements per row and per
column, by which $N=N_{H}N_{V}$.

Given the complex analytical nature of the aggregated channel, we are
providing the following Lemma to facilitate subsequent derivations.

######  Lemma 1.

The second-order and fourth-order moments of the channel norm
${\textbf{g}}_{mj}^{\mathtt{E}}$ can be written as

$\displaystyle\mathbb{E}\\{\|{\textbf{g}}_{mj}^{\mathtt{E}}\|^{2}\\}$
$\displaystyle=L\delta_{mj},$ (4a)
$\displaystyle\mathbb{E}\\{\|{\textbf{g}}_{mj}^{\mathtt{E}}\|^{4}\\}$
$\displaystyle=L(L+1)\big{(}\delta_{mj}^{2}+\mathrm{tr}((\bar{\operatorname{\boldsymbol{\Theta}}}_{mj})^{2})\big{)},\leavevmode\nobreak\
$ (4b)

where
$\delta_{mj}\triangleq\beta_{mj}^{\mathtt{E}}+\mathrm{tr}(\bar{\operatorname{\boldsymbol{\Theta}}}_{mj})$,
and
$\bar{\operatorname{\boldsymbol{\Theta}}}_{mj}\triangleq\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}$.

###### Proof.

See Appendix A. ∎

Figure 1: The proposed CF-mMIMO RIS-assisted SWIPT system.

###  II-B Uplink Training for Channel Estimation

In the training phase, the channels are estimated at the APs. To this end, all
the IUs and EUs transmit orthogonal pilot sequences
$\operatorname{\boldsymbol{\varphi}}_{k}^{I}$ and
$\operatorname{\boldsymbol{\varphi}}_{j}^{E}$ of length $\tau$ symbols, which
requires $\tau\geq|\mathcal{K}|+|\mathcal{J}|$. Let
$\operatorname{\boldsymbol{\Phi}}_{K}=[\operatorname{\boldsymbol{\varphi}}_{1}^{I},\ldots,\operatorname{\boldsymbol{\varphi}}_{K}^{I}]$
and
$\operatorname{\boldsymbol{\Phi}}_{J}=[\operatorname{\boldsymbol{\varphi}}_{1}^{E},\ldots,\operatorname{\boldsymbol{\varphi}}_{J}^{E}]$
denote the pilot matrices transmitted by the IUs and EUs. By applying similar
analytical methodology as in [13], the minimum mean square error (MMSE)
estimates of ${\textbf{g}}_{mk}^{\mathtt{I}}$ and
${\textbf{g}}_{mj}^{\mathtt{E}}$ are
$\hat{{\textbf{g}}}_{mk}^{\mathtt{I}}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\gamma_{mk}^{\mathtt{I}}{\textbf{I}}_{L})$
and
$\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\gamma_{mj}^{\mathtt{E}}{\textbf{I}}_{L})$,
where

$\displaystyle\gamma_{mk}^{\mathtt{I}}\triangleq\mathbb{E}\Big{\\{}\big{\|}[\hat{{\textbf{g}}}_{mk}^{\mathtt{I}}]_{l}\big{\|}^{2}\Big{\\}}=\frac{\tau\rho_{u}(\beta_{mk}^{\mathtt{I}})^{2}}{\tau\rho_{u}\beta_{mk}^{\mathtt{I}}+1},$
(5)
$\displaystyle\gamma_{mj}^{\mathtt{E}}\triangleq\mathbb{E}\Big{\\{}\big{\|}[\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}]_{l}\big{\|}^{2}\Big{\\}}=\frac{\tau\rho_{u}(\delta_{mj})^{2}}{\tau\rho_{u}\delta_{mj}+1},$
(6)

where $\rho_{u}$ denotes the uplink training signal-to-noise ratio (SNR). Due
to the independence property of MMSE estimation, the estimation errors are
distributed as
$\tilde{{\textbf{g}}}_{mk}^{\mathtt{I}}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},(\beta_{mk}^{\mathtt{I}}-\gamma_{mk}^{\mathtt{I}}){\textbf{I}}_{L})$
and
$\tilde{{\textbf{g}}}_{mj}^{\mathtt{E}}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},(\delta_{mj}-\gamma_{mj}^{\mathtt{E}}){\textbf{I}}_{L})$.

###  II-C Downlink SWIPT

During downlink transmission phase, the $K$ IUs and $J$ EUs are ensured
simultaneous downlink service. Following the channel acquisition from the
training phase, $M$ APs initially select between the WIT or WPT operation
mode. Subsequently, according to the mode of operation, precoding schemes are
appropriately employed for WIT and WPT operations. Specifically, the decision
on which operation mode is assigned to each AP is made based on the network
design requirements. We denote the binary variable $a_{m}$ as the operation
mode indicator, which can be formulated as

$\displaystyle\leavevmode\nobreak\ a_{m}$
$\displaystyle\triangleq\begin{cases}1&\text{if the $m$-th AP operates as
I-AP}\\\ 0&\text{if the $m$-th AP operates as E-AP}.\end{cases}$ (7)

We denote $x_{\mathtt{I},k}$ and $x_{\mathtt{E},j}$, where
$\mathbb{E}\\{|x_{\mathtt{I},k}|^{2}\\}=\mathbb{E}\\{|x_{\mathtt{E},j}|^{2}\\}=1$,
as the data and energy symbol transmitted to the $k$-th IU and $j$-th EU,
respectively. Then, the transmitted signal from the $m$-th AP is

$\displaystyle\leavevmode\nobreak\ {\bf x}_{m}$
$\displaystyle=\sqrt{a_{m}\rho_{d}}\sum\nolimits_{k\in\mathcal{K}}\sqrt{\eta_{mk}^{\mathtt{I}}}{\bf
w}_{\mathrm{I},mk}x_{\mathtt{I},k}$
$\displaystyle+\sqrt{(1-a_{m})\rho_{d}}\sum\nolimits_{j\in\mathcal{J}}\sqrt{\eta_{mj}^{\mathtt{E}}}{\bf
w}_{\mathrm{E},mj}x_{\mathtt{E},j},$ (8)

where $\rho_{d}=\tilde{\rho_{d}}/\sigma_{n}^{2}$ is the maximum downlink SNR,
while $\tilde{\rho_{d}}$ denotes the transmit power and $\sigma_{n}^{2}$ is
the noise power; ${\bf w}_{\mathrm{I},mk}\in\operatorname{\mathbb{C}}^{L\times
1}$ and ${\bf w}_{\mathrm{E},mj}\in\operatorname{\mathbb{C}}^{L\times 1}$ are
the precoding vectors for the $k$-th IU and $j$-th EU, respectively, which
adhere to the constraint $\mathbb{E}\big{\\{}\big{\|}{\bf
w}_{\mathrm{I},mk}\big{\|}^{2}\big{\\}}=\mathbb{E}\big{\\{}\big{\|}{\bf
w}_{\mathrm{E},mj}\big{\|}^{2}\big{\\}}=1$. We denote the transmitted signal
from the $m$-th I-AP and the $m$-th E-AP as ${\bf x}_{\mathtt{I},m}$ and ${\bf
x}_{\mathtt{E},m}$, respectively. Furthermore, $\eta_{mk}^{\mathtt{I}}$ and
$\eta_{mj}^{\mathtt{E}}$ are the power control coefficients chosen to satisfy
the power constraint at each AP, i.e.,

$\displaystyle a_{m}\mathbb{E}\big{\\{}\big{\|}{\bf
x}_{\mathtt{I},m}\big{\|}^{2}\big{\\}}+(1-a_{m})\mathbb{E}\big{\\{}\big{\|}{\bf
x}_{\mathtt{E},m}\big{\|}^{2}\big{\\}}\leq\rho_{d}.$ (9)

With the transmitted signal ${\bf x}_{m}$, the $k$-th IU and the $j$-th EU
then receive

$\displaystyle
y_{\mathtt{I},k}\\!=\\!\sqrt{a_{m}\rho_{d}}\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\sqrt{\eta_{mk}^{\mathtt{I}}}({\textbf{g}}_{mk}^{\mathtt{I}})^{H}{\bf
w}_{\mathrm{I},mk}x_{\mathtt{I},k}$ (10)
$\displaystyle+\sqrt{a_{m}\rho_{d}}\sum\nolimits_{k^{\prime}\neq
k}^{K}\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\sqrt{\eta_{mk^{\prime}}^{\mathtt{I}}}({\textbf{g}}_{mk}^{\mathtt{I}})^{H}{\bf
w}_{\mathrm{I},mk^{\prime}}x_{\mathtt{I},k^{\prime}}$
$\displaystyle+\\!\sqrt{(\\!1\\!-\\!a_{m}\\!)\rho_{d}}\\!\sum\nolimits_{j\in\mathcal{J}}\\!\\!\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\\!\\!\sqrt{\eta_{mj}^{\mathtt{E}}}({\textbf{g}}_{mk}^{\mathtt{I}})^{\\!H}\\!{\bf
w}_{\mathrm{E},mj}x_{\mathtt{E},j}\\!\\!+\\!n_{d,k},$

and,

$\displaystyle
y_{\mathtt{E},j}\\!=\\!\sqrt{(1\\!-\\!a_{m})\rho_{d}}\\!\sum\nolimits_{j^{\prime}\in\mathcal{J}}\\!\\!\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\\!\\!\sqrt{\eta_{mj^{\prime}}^{\mathtt{E}}}({\textbf{g}}_{mj}^{\mathtt{E}})^{\\!H}{\bf
w}_{\mathrm{E},mj^{\prime}}x_{\mathtt{E},j^{\prime}}$
$\displaystyle+\\!\sqrt{a_{m}\rho_{d}}\\!\sum\nolimits_{k\in\mathcal{K}}\\!\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\\!\\!\sqrt{\eta_{mk}^{\mathtt{I}}}({\textbf{g}}_{mj}^{\mathtt{E}})^{\\!H}{\bf
w}_{\mathrm{I},mk}x_{\mathtt{I},k}\\!+\\!n_{d,j},$ (11)

respectively, where $n_{d,k}$ and
$n_{d,j}\sim\operatorname{\mathcal{CN}}(0,1)$ are the corresponding additive
white Gaussian noises.

##  III Precoding Design and Performance Evaluation

###  III-A Precoding Design

Different beamforming designs can be explored at the system design level to
meet the performance requirements of the system. In order to strike a balance
between complexity and performance, we consider the PZF technique. Then, local
PZF precoding is considered at the I-APs and PMRT is applied at the E-APs.
This precoding design capitalizes on the effectiveness of the ZF principle in
suppressing inter-user interference, making it nearly optimal for information
transmission [14]. Additionally, power transfer with MRT proves to be optimal
for WPT systems, especially with a large number of antennas [15]. However, it
is worth noting that in our system, IUs also encounter non-coherent
interference, denoted as $\mathrm{EUI}_{kj}$, from EUs due to the simultaneous
transmission of information and energy signals. To mitigate the non-coherent
interference experienced from energy signals sent to EUs, MRT can be
implemented in the orthogonal complement of the IUs’ channel space, which is
called PMRT. Let
$\hat{{\textbf{G}}}_{m}^{\mathtt{I}}=\big{[}\hat{{\textbf{g}}}_{m1}^{\mathtt{I}},\ldots,\hat{{\textbf{g}}}_{mK}^{\mathtt{I}}\big{]}\in\operatorname{\mathbb{C}}^{L\times
K}$ and
$\hat{{\textbf{G}}}_{m}^{\mathtt{E}}=\big{[}\hat{{\textbf{g}}}_{m1}^{\mathtt{E}},\ldots,\hat{{\textbf{g}}}_{mJ}^{\mathtt{E}}\big{]}\in\operatorname{\mathbb{C}}^{L\times
J}$ be the matrices of the estimated channels between the $m$-th AP and all
IUs, and all EUs, respectively. Accordingly, the PZF and PMRT precoders at the
$m$-th AP towards the $k$-th IU and the $j$-th EU can be designed as

$\displaystyle{\bf w}_{\mathrm{I},mk}^{\mathsf{PZF}}$
$\displaystyle=\alpha_{mk}^{\mathsf{PZF}}{\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\Big{(}\big{(}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\big{)}^{H}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\Big{)}^{-1}{\bf
e}_{k}^{I}},\leavevmode\nobreak\ $ (12a) $\displaystyle{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}$
$\displaystyle=\alpha_{mj}^{\mathsf{PMRT}}{{\bf
B}_{m}\hat{{\textbf{G}}}_{m}^{\mathtt{E}}{\bf e}_{j}^{E}},\leavevmode\nobreak\
$ (12b)

where

$\displaystyle\alpha_{mk}^{\mathsf{PZF}}$
$\displaystyle\triangleq\Big{(}\mathbb{E}\Big{\\{}\big{\|}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\big{(}\big{(}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\big{)}^{H}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\big{)}^{-1}{\bf
e}_{k}^{I}\big{\|}^{2}\Big{\\}}\Big{)}^{-\frac{1}{2}}$
$\displaystyle{\alpha_{mj}^{\mathsf{PMRT}}}$
$\displaystyle\triangleq\Big{(}\mathbb{E}\big{\\{}\big{\|}{\bf
B}_{m}\hat{{\textbf{G}}}_{m}^{\mathtt{E}}{\bf
e}_{j}^{E}\big{\|}^{2}\big{\\}}\Big{)}^{-\frac{1}{2}},$ (13)

with ${\bf e}_{k}^{I}$ and ${\bf e}_{j}^{E}$ being the $k$-th column of
${\textbf{I}}_{K}$ and the $j$-th column of ${\textbf{I}}_{J}$, respectively.
In addition, ${\bf B}_{m}$ denotes the projection matrix onto the orthogonal
complement of $\hat{{\textbf{G}}}_{m}^{\mathtt{I}}$ so that
$\big{(}\hat{{\textbf{g}}}_{mk}^{\mathtt{I}}\big{)}^{H}{\bf
B}_{m}=\boldsymbol{0}$. Thus, ${\bf B}_{m}$ can be computed as

$\displaystyle{\bf
B}_{m}={\textbf{I}}_{L}-\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\Big{(}\big{(}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\big{)}^{H}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\Big{)}^{-1}\big{(}\hat{{\textbf{G}}}_{m}^{\mathtt{I}}\big{)}^{H}.$
(14)

According to [14], under the consideration of independent Rayleigh fading
channels, the analytical terms of the normalization factors in (12a) and (12b)
can be calculated as
$\big{(}\alpha_{mk}^{\mathsf{PZF}}\big{)}^{2}={(L-K)\gamma_{mk}^{\mathtt{I}}}$
and
$\big{(}\alpha_{mj}^{\mathsf{PMRT}}\big{)}^{2}=\big{(}(L-K)\gamma_{mj}^{\mathtt{E}}\big{)}^{-1}$.

###  III-B Spectral Efficiency and Average Harvested Energy

Upon receiving the signal (10), the $k$-th IU detects its desired symbol
$x_{\mathtt{I},k}$. In the absence of a downlink training phase, the IUs avail
of channel statistics, specifically the average effective channel gain, to
detect their desired symbols. For the downlink SE analysis at the IUs, we
apply the hardening bound [13]. Consequently, the received signal at the
$j$-th IU can be expressed as

$\displaystyle\leavevmode\nobreak\ y_{\mathtt{I},k}$
$\displaystyle=\mathrm{DS}_{k}x_{\mathtt{I},k}+\mathrm{BU}_{k}x_{\mathtt{I},k}+\sum\nolimits_{k^{\prime}\in\operatorname{\mathcal{K}}\setminus
k}\mathrm{IUI}_{kk^{\prime}}x_{\mathtt{I},k^{\prime}}$ $\displaystyle\hskip
10.00002pt+\sum\nolimits_{j\in\operatorname{\mathcal{J}}}\mathrm{EUI}_{kj}x_{\mathtt{E},j}+n_{k},\leavevmode\nobreak\
\forall k\in\operatorname{\mathcal{K}},$ (15)

where $\mathrm{DS}_{k}$, $\mathrm{BU}_{k}$, $\mathrm{IUI}_{kk^{\prime}}$, and
$\mathrm{EUI}_{kj}$ represent the desired signal, the beamforming gain
uncertainty, the interference cause by the $k^{\prime}$-th IU, and the
interference caused by the $j$-th EU, respectively, given by

$\displaystyle\leavevmode\nobreak\ \mathrm{DS}_{k}$
$\displaystyle\triangleq\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\sqrt{a_{m}\rho_{d}\eta_{mk}^{\mathtt{I}}}\mathbb{E}\big{\\{}({\textbf{g}}_{mk}^{\mathtt{I}})^{H}{\bf
w}_{\mathrm{I},mk}^{\mathsf{PZF}}\big{\\}},$ (16a)
$\displaystyle\mathrm{BU}_{k}$
$\displaystyle\triangleq\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\sqrt{a_{m}\rho_{d}\eta_{mk}^{\mathtt{I}}}\Bigl{(}({\textbf{g}}_{mk}^{\mathtt{I}})^{H}{\bf
w}_{\mathrm{I},mk}^{\mathsf{PZF}}$ $\displaystyle\hskip
10.00002pt-\mathbb{E}\big{\\{}({\textbf{g}}_{mk}^{\mathtt{I}})^{H}{\bf
w}_{\mathrm{I},mk}^{\mathsf{PZF}}\big{\\}}\Bigl{)}\leavevmode\nobreak\ ,$
(16b) $\displaystyle\mathrm{IUI}_{kk^{\prime}}$
$\displaystyle\triangleq\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\sqrt{a_{m}\rho_{d}\eta_{mk^{\prime}}^{\mathtt{I}}}({\textbf{g}}_{mk}^{\mathtt{I}})^{H}{\bf
w}_{\mathrm{I},mk^{\prime}}^{\mathsf{PZF}},$ (16c)
$\displaystyle\mathrm{EUI}_{kj}$
$\displaystyle\triangleq\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\sqrt{(1-a_{m})\rho_{d}\eta_{mj}^{\mathtt{E}}}({\textbf{g}}_{mk}^{\mathtt{I}})^{H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}.$ (16d)

The corresponding DL SE in [bit/s/Hz] for the $k$-th IU can be obtained as

$\displaystyle\leavevmode\nobreak\ \mathrm{SE}_{k}$
$\displaystyle=\Big{(}1\\!-\\!\frac{\tau}{\tau_{c}}\Big{)}\log_{2}\left(1\\!+\mathrm{SINR}_{k}\right),$
(17)

where $\mathrm{SINR}_{k}$ is the effective signal- to-interference-and-noise
(SINR), given by $\mathrm{SINR}_{k}=$

$\displaystyle\\!\frac{\big{|}\mathrm{DS}_{k}\big{|}^{2}}{\mathbb{E}\big{\\{}\big{|}\mathrm{BU}_{k}\big{|}^{2}\\!\big{\\}}\\!+\\!\\!\sum_{k^{\prime}\in\mathcal{K}\setminus
k}\\!\mathbb{E}\big{\\{}\big{|}\mathrm{IUI}_{kk^{\prime}}\big{|}^{2}\\!\big{\\}}\\!+\\!\\!\sum_{j\in\mathcal{J}}\\!\mathbb{E}\big{\\{}\big{|}\mathrm{EUI}_{kj}\big{|}^{2}\\!\big{\\}}\\!+\\!1}.$
(18)

To characterize the HE, a non-linear energy harvesting model with the
sigmoidal function is used. Therefore, the total HE at EU $j$ is given by [16]

$\displaystyle\leavevmode\nobreak\ \Phi_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}}\big{)}=\frac{\Lambda_{j}\big{(}\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}})\big{)}-\phi\nu}{1-\nu},\leavevmode\nobreak\
\forall j\in\mathcal{J},$ (19)

where ${\bf a}$ is an indicator vector, whose entries are $a_{m}$;
$\operatorname{\boldsymbol{\eta}^{\mathtt{I}}}=[\eta_{m1}^{\mathtt{I}},\ldots,\eta_{mK}^{\mathtt{I}}]$;
$\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}=[\eta_{m1}^{\mathtt{E}},\ldots,\eta_{mJ}^{\mathtt{E}}]$;
$\phi$ is the maximum output DC power; $\nu=\frac{1}{1+\exp(\xi\chi)}$ is a
constant to guarantee a zero input/output response, while $\xi$ and $\chi$ are
constant related parameters that depend on the circuit. Moreover,
$\Lambda\big{(}\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}})\big{)}$
is the traditional logistic function, given by

$\displaystyle\leavevmode\nobreak\ \Lambda_{j}\big{(}\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}})\big{)}$
$\displaystyle\\!=\\!\\!\frac{\phi}{1\\!+\\!\exp\big{(}\\!-\xi\big{(}\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}})\\!-\\!\chi\big{)}\\!\big{)}}.$
(20)

Note that $\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}})$
denotes the received RF energy at EU $j$, $\forall j\in\mathcal{J}$. We
consider the average of the harvested energy as the performance metric of the
WPT operation, given by

$\displaystyle\leavevmode\nobreak\ \mathbb{E}\big{\\{}\Phi_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}}\big{)}\big{\\}}=\frac{\mathbb{E}\big{\\{}\Lambda_{j}\big{(}\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}})\big{)}\big{\\}}-\phi\nu}{1-\nu}.$
(21)

Finding $\mathbb{E}\big{\\{}\Lambda_{j}\big{(}\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}})\big{)}\big{\\}}$
is complicated if not impossible. Nevertheless, since the logistic function in
(20) is a convex function of $\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}})$,
by using Jensen’s inequality, we have

$\displaystyle\mathbb{E}\left\\{\Lambda_{j}\left(\mathrm{E}_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}\right)\right\\}$
$\displaystyle\geq\Lambda_{j}\left(\mathbb{E}\left\\{\mathrm{E}_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}\right\\}\right)$
$\displaystyle=\Lambda_{j}\left(Q_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}\right),$
(22)

where $Q_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}=\mathbb{E}\\{\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}\\}$.
Thus, we can find a tight lower-bound on the average HE at the EU
$j\in\mathcal{J}$. By using (II-C), we have

$\displaystyle\mathbb{E}\big{\\{}\mathrm{E}_{j}({\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}\big{\\}}=(\tau_{c}-\tau)\sigma_{n}^{2}$
$\displaystyle\hskip
5.0pt\times{\rho}_{d}\\!\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\\!\\!{(1\\!-\\!a_{m})\eta_{mj}^{\mathtt{E}}}\mathbb{E}\Big{\\{}\\!\big{|}\big{(}{\textbf{g}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}\big{|}^{2}\\!\Big{\\}}$
$\displaystyle\hskip
5.0pt+{\rho}_{d}\\!\\!\sum_{j^{\prime}\in\mathcal{J}\setminus
j}\\!\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\\!\\!{(1\\!-\\!a_{m})\eta_{mj^{\prime}}^{\mathtt{E}}}\mathbb{E}\Big{\\{}\\!\big{|}\big{(}{\textbf{g}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{E},mj^{\prime}}^{\mathsf{PMRT}}\big{|}^{2}\\!\Big{\\}}$
$\displaystyle\hskip
5.0pt+\\!{\rho}_{d}\sum\nolimits_{k\in\mathcal{K}}\\!\sum\nolimits_{m\in\operatorname{\mathcal{M}}}\\!{a_{m}\eta_{mk}^{\mathtt{I}}}\mathbb{E}\Big{\\{}\big{|}\big{(}{\textbf{g}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{I},mk}^{\mathsf{PZF}}\big{|}^{2}\Big{\\}}\\!+\\!1.$ (23)

We provide closed-form expressions for the SE and average HE under the PPZF
precoding scheme.

######  Proposition 1.

The ergodic SE for the $k$-th IU, achieved by PZF precoding at the I-APs and
PMRT at the E-APs is given by (17), where the effective SINR is given in
closed-form by (24) at the top of next page.

$\displaystyle\mathrm{SINR}_{k}^{\mathsf{PZF}}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}=\\!\frac{(L-K)\Big{(}\sum_{m\in\operatorname{\mathcal{M}}}\sqrt{a_{m}\eta_{mk}^{\mathtt{I}}\gamma_{mk}^{\mathtt{I}}}\Big{)}^{2}}{\sum_{k^{\prime}\in\mathcal{K}}\sum_{m\in\operatorname{\mathcal{M}}}a_{m}\eta_{mk^{\prime}}^{\mathtt{I}}(\beta_{mk}^{\mathtt{I}}-\gamma_{mk}^{\mathtt{I}})+\sum_{j\in\mathcal{J}}\sum_{m\in\operatorname{\mathcal{M}}}{(1-a_{m})\eta_{mj}^{\mathtt{E}}(\beta_{mk}^{\mathtt{I}}-\gamma_{mk}^{\mathtt{I}})}+1/\rho_{d}}.$
(24)

###### Proof.

The proof is omitted due to space limitations. ∎

######  Proposition 2.

The average HE for the $j$-th EU, achieved by PMRT precoding at the E-APs and
PZF at the I-APs is bounded by

$\displaystyle\leavevmode\nobreak\ \mathbb{E}\big{\\{}\Phi_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}}\big{)}\big{\\}}\geq\frac{\Lambda_{j}\left(Q_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}\right)-\phi\nu}{1-\nu},$
(25)

where $Q_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}$
is given by (2) at the top of next page.

$\displaystyle\leavevmode\nobreak\ Q_{j}\big{(}{\bf
a},\operatorname{\boldsymbol{\eta}^{\mathtt{I}}},\operatorname{\boldsymbol{\eta}^{\mathtt{E}}}\big{)}$
$\displaystyle=(\tau_{c}-\tau)\sigma_{n}^{2}\rho_{d}\Big{(}\big{(}L-K+1\big{)}\sum\nolimits_{m\in\operatorname{\mathcal{M}}}(1-a_{m})\eta_{mj}^{\mathtt{E}}\gamma_{mj}^{\mathtt{E}}+\sum\nolimits_{j^{\prime}\in\mathcal{J}\setminus
j}\sum\nolimits_{m\in\operatorname{\mathcal{M}}}(1-a_{m})\eta_{mj^{\prime}}^{\mathtt{E}}\delta_{mj}$
$\displaystyle\hskip
80.00012pt+\sum\nolimits_{k\in\mathcal{K}}\sum\nolimits_{m\in\operatorname{\mathcal{M}}}a_{m}\eta_{mk}^{\mathtt{I}}\delta_{mj}+1/\rho_{d}\Big{)}.$
(26)

###### Proof.

See Appendix B. ∎

By examining the expressions for the SE (24) and HE (25), we note that
utilizing the current degrees-of-freedom for system design, encompassing the
design of the scattering matrix for the BD-RIS, the AP mode selection vector,
and the power control coefficients at the AP, enables us to fulfill various
design objectives in CF-mMIMO SWIPT systems.

Figure 2: Impact of the number of antennas per AP ($L$) on the average HE and
SE ($N=40,ML=480,K=3,J=5$).

##  IV Numerical Results and Discussion

In this section, numerical results are presented to illustrate the performance
of our proposed scenario. We consider a BD-RIS-assisted CF-mMIMO SWIPT
network, where the APs are uniformly distributed in an area of $1\times
1\text{km}^{2}$. The height of the APs, BD-RIS, and users is $15$ m, $30$ m,
and $1.65$ m, respectively. Moreover, we set $\tau_{c}=200$, $\tau=K+J$,
$\sigma_{n}^{2}=-92$ dBm, $\rho_{d}=1$ W and $\rho_{u}=0.2$ W. In addition, we
set the non-linear energy harvesting parameters as $\xi=150$, $\chi=0.014$,
and $\phi=0.024$ W [15]. The large-scale fading coefficients are generated
following the three-slope propagation model from [13]. The shadow fading
follows a log-normal distribution with a standard deviation of $8$ dB. In our
simulations, we compare three different design scenarios:

  * •

Heuristic scattering matrix design at the BD-RIS (Heuristic BD-RIS), where the
symmetric-unitary projection scheme in [17] is adopted. To this end, we first
construct $M$ matrices ${\bf
F}_{m}\triangleq{\textbf{H}}_{m}{\textbf{H}}^{E}_{m}{\textbf{Z}}_{m}^{H}$,
where
${\textbf{H}}^{E}_{m}=[{\textbf{h}}_{m1}^{\mathtt{E}},\ldots,{\textbf{h}}_{mJ}^{\mathtt{E}}]\in\operatorname{\mathbb{C}}^{M\times
J}$ and ${\textbf{Z}}_{m}=[{\bf z}_{1},\ldots,{\bf
z}_{J}]\in\operatorname{\mathbb{C}}^{N\times J}$. Next, the symmetric
projection of ${\bf F}_{m}$ is obtained by averaging ${\bf F}_{m}$ and its
transpose across various small-scale channel realizations. To meet the unitary
constraint, the rank $r_{m}$ and singular value decompositions $[{\bf
U}_{m},{\bf V}_{m}]$ of ${\bf F}_{m}$ are computed. Then, ${\bf U}_{m},{\bf
V}_{m}$ are partitioned into $[{\bf U}_{m}^{r_{m}},{\bf U}_{m}^{N-r_{m}}]$,
and $[{\bf V}_{m}^{r_{m}},{\bf V}_{m}^{N-r_{m}}]$, respectively. This approach
allows for the heuristic design of $\operatorname{\boldsymbol{\Theta}}$ via
the symmetric-unitary projection, which can be computed as
$\operatorname{\boldsymbol{\Theta}}_{m}\triangleq\hat{{\bf U}}_{m}{\bf
V}^{H}$, where $\hat{{\bf U}}_{m}=[{\bf U}_{m}^{r_{m}},\big{(}{\bf
V}_{m}^{N-r_{m}}\big{)}^{*}]$. The trace products
$\mathrm{tr}(\operatorname{\boldsymbol{\Theta}}^{H}{\bf
R}\operatorname{\boldsymbol{\Theta}}{\bf R})$ and the indexes of
$\operatorname{\boldsymbol{\Theta}}_{m}$ of each iteration are stored, such
that

$\displaystyle\operatorname{\boldsymbol{\Theta}}_{m}^{*}=\underset{\operatorname{\boldsymbol{\Theta}}_{m}}{\operatorname*{arg\,max}}{\mathrm{tr}(\operatorname{\boldsymbol{\Theta}}^{H}{\bf
R}\operatorname{\boldsymbol{\Theta}}{\bf R})}.$

This process iteratively refines the design of the scattering matrix at the
BD-RIS.

  * •

BD-RIS with a random scattering matrix (Random BD-RIS), where a discrete
Fourier transformation matrix is used to generate random
$\operatorname{\boldsymbol{\Theta}}$ matrices, which are then normalized by
$N$.

  * •

CF-mMIMO SWIPT without BD-RIS (No BD-RIS), to study the role of BD-RIS in the
SWIPT process.

For both Heuristic BD-RIS and Random BD-RIS, the scattering matrix satisfies
the constraints
$\operatorname{\boldsymbol{\Theta}}=\operatorname{\boldsymbol{\Theta}}^{T}$,
and
$\operatorname{\boldsymbol{\Theta}}^{H}\operatorname{\boldsymbol{\Theta}}=\mathbf{I}_{N}$
[11, 12]. We assume that the APs’ operation mode selection parameters
($a_{m}$, $\forall m\in\operatorname{\mathcal{M}}$) are randomly assigned,
while the downlink power control coefficients are
$\eta_{mk}^{\mathtt{I}}=\big{(}\sum_{k\in\mathcal{K}}\gamma_{mk}^{\mathtt{I}}\big{)}^{-1}$,
$\forall m,k$, and
$\eta_{mj}^{\mathtt{E}}=\big{(}\sum_{j\in\mathcal{J}}\gamma_{mj}^{\mathtt{E}}\big{)}^{-1}$,
$\forall m,j$.

Figure 2 shows the average sum HE and the per-IU SE achieved by three
scenarios versus the number of antennas per AP. We note that, for a fixed
number of service antennas (i.e., $ML=480$), the number of APs decreases as
the number of antennas per AP $L$ increases. From the standpoint of average
sum heuristic HE, it is evident that Heuristic BD-RIS outperforms Random BD-
RIS, achieving a noticeable gain of up to $12\%$. This result showcases the
significance of an optimization scheme for BD-RIS scattering matrix design. On
the other hand, it is evident that without BD-RIS, WPT fails to meet the
energy requirements of the EUs. Finally, we observe that, by increasing $L$
and consequently decreasing $M$, both the per-IU SE and average sum HE are
degraded. The decrease rate of the SE is more dominant compared to the
decrease rate of the average sum HE, which can be interpreted as a consequence
of employing BD-RISs to assist the EUs.

Figure 3 shows the average HE of the EUs and per-IU SE versus the number of
APs for fixed $ML=480$. As the number of APs increases from $M=8$ to $M=24$,
there is a notable enhancement in the average HE, with improvements of $14\%$
and $20\%$ for Random BD-RIS and Heuristic BD-RIS, respectively. Beyond this,
the improvement continues, albeit at a more gradual rate. These results
confirm that both the number of APs and the per-AP antenna count should be
carefully selected to improve both SE and average HE in the considered CF-
mMIMO SWIPT system.

Figure 3: Impact of the number of APs ($M$) on the average HE and SE
($N=40,ML=480,K=3,J=5$).

##  V Conclusion

We evaluated the SE and average HE performance in a novel setup of a CF-mMIMO
SWIPT network assisted by a BD-RIS. In this setup, the APs function as either
I-APs or E-APs, serving two distinct groups of IUs and EUs simultaneously,
with the BD-RIS strategically positioned near the energy zone to facilitate
the WPT towards the EUs. Our simulation results demonstrated the superior
performance of the BD-RIS assisted architectures compared to the conventional
designs without BD-RIS. We further observed that the proposed Heuristic BD-RIS
can achieve up $12\%$ gain over the Random BD-RIS design for a fixed number of
service antennas. For future work, integrating a robust optimization scheme to
jointly optimize the design variables could pave the way for next-generation
communication systems.

##  Appendix A Proof of Lemma 1

The following lemma is useful to prove Lemma 1. The proof is omitted due to
space limitations.

######  Lemma 2.

For an $N\times L$ matrix ${\bf X}=[{\bf x}_{1},...,{\bf x}_{L}]$ with ${\bf
x}_{l}\in\operatorname{\mathbb{C}}^{N\times 1}$ distributed as ${\bf
x}_{l}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\bar{{\bf R}})$ with
${\bf R}\in\operatorname{\mathbb{C}}^{N\times N}$, and two deterministic
matrices ${\bf M}$, ${\bf N}\in\operatorname{\mathbb{C}}^{N\times N}$, it
holds that $\mathbb{E}\bigl{\\{}{\bf X}^{H}{\bf M}{\bf X}{\bf X}^{H}{\bf
N}{\bf X}\bigl{\\}}$ is a diagonal matrix of size $L\times L$ with
$L\mathrm{tr}\big{(}{\bf M}{\bf R}{\bf N}{\bf R}\big{)}+\mathrm{tr}\big{(}{\bf
R}{\bf M}\big{)}\mathrm{tr}({\bf R}{\bf N}\big{)}$ on its main diagonal.

By leveraging the statistical independence of the direct and indirect links,
the second-order moment of ${\textbf{g}}_{mj}^{\mathtt{E}}$ becomes

$\displaystyle\mathbb{E}\\{\|{\textbf{g}}_{mj}^{\mathtt{E}}\|^{2}\\}$
$\displaystyle=\mathbb{E}{\|{\textbf{h}}_{mj}^{\mathtt{E}}\|^{2}}+\mathbb{E}\Big{\\{}\|{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}{\bf
z}_{j}\|^{2}\Big{\\}}$
$\displaystyle=L\beta_{mj}^{\mathtt{E}}+\mathbb{E}\Big{\\{}{\bf
z}^{H}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}{\bf
z}_{j}\Big{\\}}$
$\displaystyle\overset{\mathrm{(a)}}{=}L\beta_{mj}^{\mathtt{E}}+\mathrm{tr}\Big{(}\operatorname{\boldsymbol{\Theta}}^{H}\mathbb{E}\big{\\{}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\big{\\}}\operatorname{\boldsymbol{\Theta}}\mathbb{E}\big{\\{}{\bf
z}_{j}{\bf z}^{H}_{j}\big{\\}}\Big{)}$
$\displaystyle=L\Big{(}\beta_{mj}^{\mathtt{E}}+\mathrm{tr}\Big{(}\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\Big{)}\Big{)},$
(27)

where (a) exploits $\mathbb{E}\\{{\bf
x}^{H}{\textbf{y}}\\}=\mathbb{E}\\{\mathrm{tr}({\textbf{y}}{\bf x}^{H})\\}$,
followed by the property of expectation of independent channels
${\textbf{H}}_{m}$ and ${\bf z}_{j}$.

To derive the fourth-order moment
$\mathbb{E}\\{\|{\textbf{g}}_{mj}^{\mathtt{E}}\|^{4}\\}$, we define
$a\\!=\\!\|{\textbf{h}}_{mj}^{\mathtt{E}}\|^{2}$,
$b\\!=\\!({\textbf{h}}_{mj}^{\mathtt{E}})^{H}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}{\bf
z}_{j}$, $c\\!=\\!{\bf
z}_{j}^{H}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{h}}_{mj}^{\mathtt{E}}$,
and $d\\!=\\!\|{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}{\bf
z}_{j}\|^{2}$. Thus, we have

$\displaystyle\leavevmode\nobreak\
\mathbb{E}\\{\|{\textbf{g}}_{mj}^{\mathtt{E}}\|^{4}\\}\\!=$
$\displaystyle\mathbb{E}\\{|a|^{2}\\}\\!+\\!\mathbb{E}\\{|b|^{2}\\}\\!+\\!\mathbb{E}\\{|c|^{2}\\}\\!+\\!\mathbb{E}\\{|d|^{2}\\}\\!+\\!2\mathbb{E}\\{ad\\}.$
(28)

We notice that

$\leavevmode\nobreak\
\mathbb{E}\\{|a|^{2}\\}=\mathbb{E}\\{\|{\textbf{h}}_{mj}^{\mathtt{E}}\|^{4}\\}=L(L+1)(\beta_{mj}^{\mathtt{E}})^{2}.$
(29)

Moreover, $\mathbb{E}\\{|b|^{2}\\}=\mathbb{E}\\{|c|^{2}\\}$, which can be
calculated as

$\displaystyle\leavevmode\nobreak\ \mathbb{E}\\{|b|^{2}\\}$
$\displaystyle=\beta_{mj}^{\mathtt{E}}\mathrm{tr}(\operatorname{\boldsymbol{\Theta}}^{H}\mathbb{E}\\{{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\\}\operatorname{\boldsymbol{\Theta}}\mathbb{E}\\{{\bf
z}_{j}{\bf z}_{j}^{H}\\})$
$\displaystyle=L\beta_{mj}^{\mathtt{E}}\bigl{(}\mathrm{tr}(\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j})\bigl{)}.$
(30)

In order to derive $\mathbb{E}\big{\\{}|d|^{2}\big{\\}}$, we first express it
as

$\displaystyle\leavevmode\nobreak\ \mathbb{E}\big{\\{}|d|^{2}\big{\\}}$
$\displaystyle=\mathbb{E}\Big{\\{}\mathrm{tr}(\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}})$
$\displaystyle\hskip
10.00002pt+\big{|}\mathrm{tr}(\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}})\big{|}^{2}\Big{\\}},$
(31)

where [18, Lemma 2] is applied for a complex Gaussian random vector ${\bf
z}_{j}\sim\operatorname{\mathcal{CN}}(\boldsymbol{0},\boldsymbol{\Psi}_{j})$
and a given deterministic matrix
$\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}$.

The first expectation term in (A) can be obtained as

$\displaystyle\mathbb{E}\Big{\\{}\mathrm{tr}\big{(}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}\big{)}\Big{\\}}$
$\displaystyle=\mathrm{tr}\Big{(}\mathbb{E}\Big{\\{}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}\Big{\\}}\Big{)}$
$\displaystyle\\!\stackrel{{\scriptstyle(a)}}{{=}}\\!L\bigg{(}L\mathrm{tr}\big{(}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m}\big{)}$
$\displaystyle\hskip
20.00003pt+\mathrm{tr}\big{(}\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}\big{)}\mathrm{tr}(\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}\big{)}\Big{)},$
$\displaystyle\\!=\\!L\Big{(}L\mathrm{tr}\big{(}(\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m})^{\\!2}\big{)}\\!+\\!\big{(}\mathrm{tr}(\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{\\!H})\big{)}^{\\!2}\Big{)},$
(32)

where we have used Lemma 2.

Due to space limitations, we present the conclusive results for the second
terms in (A) as

$\displaystyle\mathbb{E}\Big{\\{}\Big{|}\mathrm{tr}\big{(}{\textbf{H}}_{m}^{H}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}{\textbf{H}}_{m}\big{)}\Big{|}^{2}\Big{\\}}=L\mathrm{tr}\bigl{(}(\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m})^{2}\bigl{)}$
$\displaystyle\hskip
60.00009pt+L^{2}\bigl{(}\mathrm{tr}(\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j}\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m})\bigl{)}^{2}.$
(33)

Finally, using the independence property between the direct and indirect
channels, $\mathbb{E}\\{ad\\}$ can be derived as

$\displaystyle\mathbb{E}\\{2ad\\}=L^{2}\beta_{mj}^{\mathtt{E}}\bigl{(}\mathrm{tr}(\operatorname{\boldsymbol{\Theta}}^{H}\boldsymbol{\Omega}_{m}\operatorname{\boldsymbol{\Theta}}\boldsymbol{\Psi}_{j})\bigl{)}.$
(34)

By substituting (29), (A), (A), (34) into (28), followed by some algebraic
transformations, Lemma 1 is obtained.

##  Appendix B Proof of Proposition 2

First, we compute the first term in (III-B) as

$\displaystyle\mathbb{E}\big{\\{}\big{|}\big{(}{\textbf{g}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}\big{|}^{2}\big{\\}}$
$\displaystyle\\!=\\!\mathbb{E}\big{\\{}\big{|}\big{(}\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}\\!+\\!\tilde{{\textbf{g}}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}\big{|}^{2}\big{\\}}$
$\displaystyle\\!=\\!\mathbb{E}\big{\\{}\big{|}\big{(}\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}\big{|}^{2}\big{\\}}.$ (35)

According to (12b),
$\mathbb{E}\big{\\{}\big{|}\big{(}\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}\big{|}^{2}\big{\\}}$ is computed as

$\displaystyle\mathbb{E}\big{\\{}\\!\big{(}\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}\\!{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}\\!\big{\\}}$
$\displaystyle\\!\\!\stackrel{{\scriptstyle(a)}}{{=}}\\!\\!\frac{\mathbb{E}\big{\\{}\\!{\bf
B}_{m}\\!\big{\\}}\mathbb{E}\big{\\{}\\!\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}({\textbf{g}}_{mj}^{\mathtt{E}})^{H}\\!\big{\\}}}{\sqrt{(L\\!-K)\gamma_{mj}^{\mathtt{E}}}}\\!\\!=\\!\\!\sqrt{(L\\!-\\!K)\gamma_{mj}^{\mathtt{E}}},$

where (a) exploits the independence between ${\textbf{g}}_{mj}^{\mathtt{E}}$
and the idempotent Hermitian matrix ${\bf B}_{m}$, where
$\mathbb{E}\big{\\{}\\!{\bf B}_{m}\\!\big{\\}}=(L-K/L){\textbf{I}}_{L}$ [19,
Lemma 1]. Moreover, by applying the fourth-order moment rule [13, Appendix
A.2.4], we derive

$\displaystyle\mathbb{E}\big{\\{}\\!\big{|}\big{(}\hat{{\textbf{g}}}_{mj}^{\mathtt{E}}\big{)}^{\\!H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}\big{|}^{2}\\!\big{\\}}$
$\displaystyle=\frac{(L-K)(L-K+1)(\gamma_{mj}^{\mathtt{E}})^{2}}{(L-K)\gamma_{mj}^{\mathtt{E}}}$
$\displaystyle=(L-K+1)\gamma_{mj}^{\mathtt{E}}.$ (36)

Following the independence property between
$\tilde{{\textbf{g}}}_{mj}^{\mathtt{E}}$ and ${\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}$, the second expectation term in (III-B),
is obtained as

$\mathbb{E}\big{\\{}(\tilde{{\textbf{g}}}_{mj}^{\mathtt{E}})^{H}{\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}}({\bf
w}_{\mathrm{E},mj}^{\mathsf{PMRT}})^{H}\tilde{{\textbf{g}}}_{mj}^{\mathtt{E}}\big{\\}}=\delta_{mj}.$
(37)

Similarly, the third expectation term can be computed as

$\mathbb{E}\big{\\{}({\textbf{g}}_{mj}^{\mathtt{E}})^{H}{\bf
w}_{\mathrm{I},mk}^{\mathsf{PZF}}({\bf
w}_{\mathrm{I},mk}^{\mathsf{PZF}})^{H}{\textbf{g}}_{mj}^{\mathtt{E}}\big{\\}}=\delta_{mj},$
(38)

due to that ${\textbf{g}}_{mj}^{\mathtt{E}}$ and ${\bf
w}_{\mathrm{I},mk}^{\mathsf{PZF}}$ are independent.

## References

  * [1] B. Clerckx, R. Zhang, R. Schober, D. W. K. Ng, D. I. Kim, and H. V. Poor, “Fundamentals of wireless information and power transfer: From RF energy harvester models to signal and system designs,” _IEEE J. Sel. Areas Commun._ , vol. 37, no. 1, pp. 4–33, Jan. 2019.
  * [2] H. Q. Ngo, L.-N. Tran, T. Q. Duong, M. Matthaiou, and E. G. Larsson, “On the total energy efficiency of cell-free massive MIMO,” _IEEE Trans. Green Commun. Netw._ , vol. 2, no. 1, pp. 25–39, Mar. 2018.
  * [3] O. T. Demir and E. Björnson, “Joint power control and LSFD for wireless-powered cell-free massive MIMO,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 3, pp. 1756–1769, Mar. 2021.
  * [4] D. L. Galappaththige, R. Shrestha, and G. A. Aruma Baduge, “Exploiting cell-free massive MIMO for enabling simultaneous wireless information and power transfer,” _IEEE Trans. Green Commun. Netw._ , vol. 5, no. 3, pp. 1541–1557, Jun. 2021.
  * [5] M. Mohammadi, L.-N. Tran, Z. Mobini, H. Q. Ngo, and M. Matthaiou, “Cell-free massive MIMO and SWIPT: Access point operation mode selection and power control,” in _Proc. IEEE GLOBECOM_ , Dec. 2023.
  * [6] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface-aided wireless communications: A tutorial,” _IEEE Trans. Commun._ , vol. 69, no. 5, pp. 3313–3351, May. 2021.
  * [7] T. V. Chien, H. Q. Ngo, S. Chatzinotas, M. D. Renzo, and B. Ottersten, “Reconfigurable intelligent surface-assisted cell-free massive MIMO systems over spatially-correlated channels,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 7, pp. 5106–5128, Jul. 2022.
  * [8] J. Dai, J. Ge, K. Zhi, C. Pan, Z. Zhang, J. Wang, and X. You, “Two-timescale transmission design for RIS-aided cell-free massive MIMO systems,” _IEEE Trans. Wireless Commun._ , pp. 1–1, Mar. 2023.
  * [9] Y. Zhao, B. Clerckx, and Z. Feng, “IRS-aided SWIPT: Joint waveform, active and passive beamforming design under nonlinear harvester model,” _IEEE Trans. Commun._ , vol. 70, no. 2, pp. 1345–1359, Feb. 2022.
  * [10] E. Shi _et al._ , “Wireless energy transfer in RIS-aided cell-free massive MIMO systems: Opportunities and challenges,” _IEEE Commun. Mag._ , vol. 60, no. 3, pp. 26–32, Mar. 2022.
  * [11] H. Li, S. Shen, M. Nerini, and B. Clerckx, “Reconfigurable intelligent surfaces 2.0: Beyond diagonal phase shift matrices,” _IEEE Commun. Mag._ , pp. 1–7, Nov. 2023.
  * [12] H. Li, S. Shen, and B. Clerckx, “Beyond diagonal reconfigurable intelligent surfaces: A multi-sector mode enabling highly directional full-space wireless coverage,” _IEEE J. Sel. Areas Commun._ , vol. 41, no. 8, pp. 2446–2460, Aug. 2023.
  * [13] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, “Cell-free massive MIMO versus small cells,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 3, pp. 1834–1850, Mar. 2017.
  * [14] G. Interdonato, M. Karlsson, E. Björnson, and E. G. Larsson, “Local partial zero-forcing precoding for cell-free massive MIMO,” _IEEE Tran. Wireless Commun.s_ , vol. 19, no. 7, pp. 4758–4774, Jul. 2020.
  * [15] A. Almradi and K. A. Hamdi, “The performance of wireless powered MIMO relaying with energy beamforming,” _IEEE Trans. Commun._ , vol. 64, no. 11, pp. 4550–4562, Nov. 2016.
  * [16] E. Boshkovska, D. W. K. Ng, N. Zlatanov, and R. Schober, “Practical non-linear energy harvesting model and resource allocation for SWIPT systems,” _IEEE Commun. Lett._ , vol. 19, no. 12, pp. 2082–2085, Dec. 2015.
  * [17] T. Fang and Y. Mao, “A low-complexity beamforming design for beyond-diagonal RIS aided multi-user networks,” _IEEE Commun. Lett._ , vol. 28, no. 1, pp. 203–207, Jan. 2024.
  * [18] E. Björnson, M. Matthaiou, and M. Debbah, “Massive MIMO with non-ideal arbitrary arrays: Hardware scaling laws and circuit-aware design,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 8, pp. 4353–4368, Aug. 2015.
  * [19] M. Elfiatoure, M. Mohammadi, H. Q. Ngo, P. J. Smith, and M. Matthaiou, “Protecting massive MIMO-radar coexistence: Precoding design and power control,” _IEEE Open J. Commun. Society_ , vol. 5, pp. 276–293, Jan. 2024.