Stacked Intelligent Metasurfaces for Holographic MIMO Aided Cell-Free Networks

Qingchao Li, Graduate Student Member, IEEE, Mohammed El-Hajjar, Senior Member, IEEE,
Chao Xu, Senior Member, IEEE, Jiancheng An, Member, IEEE, Chau Yuen, Fellow, IEEE,
and Lajos Hanzo, Life Fellow, IEEE The work of Chau Yuen was supported by the Ministry of Education, Singapore, under its Ministry of Education (MOE) Tier 2 (Award number MOE-T2EP50220-0019). The work of Lajos Hanzo was supported by the Engineering and Physical Sciences Research Council projects EP/W016605/1, EP/X01228X/1, EP/Y026721/1 and EP/W032635/1 as well as of the European Research Council’s Advanced Fellow Grant QuantCom (Grant No. 789028). (Corresponding author: Lajos Hanzo.) Qingchao Li, Mohammed El-Hajjar, Chao Xu and Lajos Hanzo are with the School of Electronics and Computer Science, University of Southampton, Southampton SO17 1BJ, U.K. (e-mail: [email protected]; [email protected]; [email protected]; [email protected]). Jiancheng An and Chau Yuen are with the School of Electrical and Electronics Engineering, Nanyang Technological University, Singapore 639798 (e-mail: [email protected]; [email protected]).

Abstract

Large-scale multiple-input and multiple-output (MIMO) systems are capable of achieving high date rate. However, given the high hardware cost and excessive power consumption of massive MIMO systems, as a remedy, intelligent metasurfaces have been designed for efficient holographic MIMO (HMIMO) systems. In this paper, we propose a HMIMO architecture based on stacked intelligent metasurfaces (SIM) for the uplink of cell-free systems, where the SIM is employed at the access points (APs) for improving the spectral- and energy-efficiency. Specifically, we conceive distributed beamforming for SIM-assisted cell-free networks, where both the SIM coefficients and the local receiver combiner vectors of each AP are optimized based on the local channel state information (CSI) for the local detection of each user equipment (UE) information. Afterward, the central processing unit (CPU) fuses the local detections gleaned from all APs to detect the aggregate multi-user signal. Specifically, to design the SIM coefficients and the combining vectors of the APs, a low-complexity layer-by-layer iterative optimization algorithm is proposed for maximizing the equivalent gain of the channel spanning from the UEs to the APs. At the CPU, the weight vector used for combining the local detections from all APs is designed based on the minimum mean square error (MMSE) criterion, where the hardware impairments (HWIs) are also taken into consideration based on their statistics. The simulation results show that the SIM-based HMIMO outperforms the conventional single-layer HMIMO in terms of the achievable rate. We demonstrate that both the HWI of the radio frequency (RF) chains at the APs and the UEs limit the achievable rate in the high signal-to-noise-ratio (SNR) region.

Index Terms:

Holographic multiple-input and multiple-output (HMIMO), stacked intelligent metasurface (SIM), cell-free network, hardware impairment (HWI).

I Introduction

In the fifth generation (5G) wireless systems, large-scale multiple-input and multiple-output (MIMO) systems have been harnessed for providing significantly increased throughput by employing a large number of antennas at the base station (BS) [1], [2], [3], [4], [5]. However, they require numerous active radio frequency (RF) chains, which results in excessive hardware cost and energy consumption. Hence, some authors have focused their attention on the design of low-cost and energy-efficient solutions [6], [7], [8], [9], [10]. As an attractive design alternative, holographic MIMO (HMIMO) solutions rely on an intelligent software reconfigurable paradigm in support of improved hardware efficiency and energy efficiency. They achieve this ambitious objective by utilizing a spatially near-continuous aperture and holographic radios having reduced power consumption and fabrication cost [11], [12], [13], [14], [15], [16]. The recent progress in channel modeling and efficient channel estimation conceived for HMIMO systems is reported in [17].

Currently, the hardware architectures of HMIMO are mainly based on reconfigurable refractive surfaces (RRS) [18], reconfigurable holographic surfaces (RHS) [19], [20], [21], [22], [23], [24], [25], [26], and dynamic metasurface antennas (DMA) [27], [28], [29], [30].

I-1 Reconfigurable refractive surfaces

It is infeasible to realize HMIMO schemes relying on a large number of conventional RF chains and active antennas due to the excessive power consumption [16]. As a remedy, Zeng et al. [18] employed a RRS illuminated by a single RF chain at the BS for creating an energy-efficient HMIMO. A substantial beamforming gain can be achieved by adjusting the coefficient of each RRS element. Both their theoretical analysis and simulation results show that the RRS-aided HMIMO has higher energy efficiency than the conventional MIMO systems using phased arrays.

I-2 Reconfigurable holographic surfaces

In [19], [20], [21], Deng et al. proposed a RHS based architecture. Their digital beamformer relying on the state-of-the-art (SoA) zero-forcing (ZF) transmit precoding method and a holographic beamformer are jointly optimized at the BS and the RHS, respectively. The simulation results showed that the RHS-based hybrid beamformer achieves a higher sum-rate than the SoA massive MIMO based hybrid beamformer relying on phase shift arrays. To reduce the configuration complexity, Hu et al. [22] proposed a holographic beamforming scheme based on amplitude-controlled RHS elements having limited resolution. They showed that the holographic beamformer associated with 2-bit quantized amplitude resolution achieves a similar sum-rate as that associated with continuous amplitude values. Furthermore, in [23] the RHS beamformer is employed for ultra-dense low-Earth-orbit (LEO) satellite communications to compensate for the severe path-loss of satellite communications. The simulation results showed that the RHS provides a more cost-effective solution for pursuing higher data rate than the phased array architecture. The impact of the number of radiation elements on the sum-rate of RHS-based LEO satellite communications is further investigated in [24]. Explicitly, the authors theoretically analyzed the minimum number of RHS elements required for the sum-rate of the RHS-aided system to exceed that of the phased array system. Furthermore, Wei et al. [25] proposed a low-complexity ZF beamformer based on utilizing Neumann series expansion to replace the matrix inversion operation in multi-user RHS-based MIMO communications. To reduce the channel state information (CSI) estimation overhead, Wu et al. [26] minimized the transmit power of the RHS-based holographic MIMO based on a two-time scale beamformer. Specifically, this holographic beamformer was designed based on the statistical CSI, and then the instantaneous CSI of the equivalent channel links was estimated and utilized for designing the digital transmit precoding matrix.

I-3 Dynamic metasurface antennas

The dynamic metasurface antenna consists of multiple microstrips and each microstrip is composed of a multitude of sub-wavelength and frequency-selective resonant metamaterial radiating elements, which can be employed to realize low-cost, power-efficient and size compact antenna arrays [27], [28], [29], [30]. In the DMA, the information is beamformed by linearly combining the radiation observation from all metamaterial elements in each microstrip. The mathematical model of DMA-based massive MIMO systems were firstly proposed by Shlezinger et al. in [27], where the fundamental limits of DMA-aided uplink communications was also investigated. To approach these limits, the alternating optimization algorithm was employed for designing practical DMAs for arbitrary multipath channels and frequency selectivity profiles. Furthermore, the achievable sum-rate of DMA-based downlink massive MIMO systems was characterized in [28], and an efficient alternating algorithm was proposed for dynamically configuring the DMA weights to maximize the achievable sum-rate. It was shown that the fundamental limits of DMA-based massive MIMO systems are comparable to those of conventional MIMO systems based on ideal antenna arrays. In [29], You et al. employed DMA to realize large-scale antenna arrays having reduced physical size, hardware cost, and power consumption. Specifically, the energy efficiency of the DMA-based massive MIMO system was maximized by the Dinkelbach transform, alternating optimization, and deterministic equivalent methods. The simulation results showed that higher energy efficiency can be achieved by the DMA-based massive MIMO architecture than by the conventional fully-digital and hybrid massive MIMO systems. Furthermore, Li et al. [30] proposed the power-efficient DMA, which operate at high frequencies and realize extremely large-scale MIMO (XL-MIMO) schemes. The DMA-based XL-MIMO architecture is composed of the conventional digital beamformer and the DMA-based holographic beamformer. Specifically, in the holographic beamformer, the DMAs can be configured based on three different modes, including continuous-amplitude configurations, binary-amplitude configurations and Lorentzian-constrained phase configurations [30]. An efficient successive convex approximation based alternating direction method of multipliers (ADMM) aided algorithm was proposed for optimizing the digital beamformer and the DMA-based holographic beamformer in an alternating manner. The associated simulation results showed that the DMA-based array has lower hardware overhead and power consumption than the conventional hybrid massive MIMO beamformer.

The above holographic MIMO architectures are based on a single-layer metasurface. To further improve both the spatial-domain gain and the beamformer’s degree-of-freedom, An et al. [31] proposed a HMIMO system based on stacked intelligent metasurfaces (SIM), which is composed of stacked multi-layer reconfigurable surfaces to carry out advanced signal processing directly in the native electromagnetic (EM) wave regime without digital beamformer. The gradient descent algorithm is employed for optimizing the phase shifts of the elements in all layers of the metasurfaces to maximize the sum-rate. The simulation results show that the SIM architecture outperforms its single-layer metasurface counterparts. The wave-based beamforming relying on the SIM is benefit of simplifying hardware architecture and improving computational efficiency. Furthermore, an SIM has the application of performing interference cancellation for multiple access and enabling integrated sensing and communications [32].

However, the above contributions assume idealized perfect RF hardware both at the BSs and the UEs, which is impractical. Furthermore, the above HMIMO architectures are tailored for cellular networks, where the cell-edge users (UEs) suffer from low data-rate both due to the increased BS-UE distance and owing to the inter-cell interference. As a design alternative, the cell-free concept mitigates the signal path loss and the inter-cell interference by deploying a set of distributed access points (APs) for cooperatively serving UEs without cell boundaries [33]. Specifically, the ZF method and the minimum mean square error (MMSE) method can be employed for designing the beamformers in the centralized algorithm. To reduce the required overhead of CSI-sharing between APs, a distributed algorithm can be employed based on the maximum ratio transmission (MRT) or the maximum ratio combining (MRC) criteria, albeit while at the cost of a performance degradation. Furthermore, in the distributed optimization algorithm of the cell-free system, the cooperation between APs is promising in terms of harnessing parallel computing resources and achieving almost the same data rate as the centralized algorithm [34]. To deal with these challenges, in this paper we propose an SIM assisted HMIMO architecture for cell-free networks, while directly considering the signal distortion resulting from the realistic hardware impairments (HWIs) of the RF chains both at the APs and the UEs. Against this background, our contributions are detailed as follows, while Table I explicitly contrasts them to the literature at a glance.

TABLE I: Contrasting the novelty of our paper to the existing MIMO techniques in the literature [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34].

	Metasurface technique	Multi-user system	Cell-free network	Transceiver hardware impairments
Our paper	SIM	✔	✔	✔
\hdashline[18]	RRS
\hdashline[19, 20, 21, 22, 23, 24]	RHS	✔
\hdashline[25]	RHS	✔
\hdashline[26]	RHS	✔
\hdashline[27]	DMA	✔
\hdashline[28]	DMA	✔
\hdashline[29]	DMA	✔
\hdashline[30]	DMA	✔
\hdashline[31, 32]	SIM	✔
\hdashline[33]	Conventional antennas	✔	✔
\hdashline[34]	Conventional antennas	✔	✔

•

We conceive an SIM-based HMIMO architecture for the uplink of a cell-free network, where the SIM is employed at the APs to attain high spectral- and energy-efficient information transfer. The distributed operation is employed for the SIM-based cell-free network. Specifically, at each AP the hybrid beamformer coefficients and the receiver combining (RC) vectors are jointly optimized to acquire a local estimate of the information arriving from the UEs. Afterwards, the central processing unit (CPU) uses the data detected by all APs to carry out the final detection of each UE’s data.
•

Since designing the hybrid beamformer coefficient of the SIM coefficients and the RC vectors at each AP is a non-convex problem, we propose a low-complexity layer-by-layer iterative optimization algorithm. Specifically, when the RC vectors are given, the coefficients of the intelligent metasurface are optimized on a layer-by-layer basis. By contrast, when the coefficients of all layers of the intelligent metasurface are given, the RC vectors are designed based on the MRC criterion. The SIM coefficients and the RC vectors of each AP are alternately optimized until reaching convergence.
•

For recovering the information gleaned from the UEs, the weight vector of the CPU used for combining the local detections arriving from all APs is designed based on the MMSE criterion harnessed for maximizing the signal-to-noise-ratio (SNR) of the received signal. Furthermore, since having RF hardware impairments at APs and UEs is inevitable, we take them into account in the weight vector design by exploiting their statistics.
•

Our numerical results show that the average achievable rate of our SIM-based HMIMO architecture in the cell-free network outperforms the conventional single-layer intelligent surface aided HMIMO. Furthermore, owing to the HWIs at APs and UEs, the achievable rate saturates at high SNRs without reaching its theoretical maximum.

The rest of this paper is organized as follows. In Section II, we present the system model, while the beamformer design is described in Section III. Our simulation results are presented in Section IV, while we conclude in Section V.

Notations: Vectors and matrices are denoted by boldface lower and upper case letters, respectively, $(\cdot)^{\text{T}}$ , $(\cdot)^{{\dagger}}$ and $(\cdot)^{\text{H}}$ represent the operation of transpose, conjugate and Hermitian transpose, respectively, $\odot$ represents the Hadamard product operation, $|a|$ and $\angle a$ denote the amplitude and angle of the complex scalar $a$ , respectively, $\|\mathbf{a}\|$ denotes the norm of the vector $\mathbf{a}$ , $\mathbb{C}^{m\times n}$ denotes the space of $m\times n$ complex-valued matrices, $\mathbf{0}_{N}$ is the $N\times 1$ zero vector, $\mathbf{I}_{N}$ represents the $N\times N$ identity matrix, $\mathbf{Diag}\{\mathbf{a}\}$ denotes a diagonal matrix having elements of $\mathbf{a}$ in order, $[\mathbf{a}]_{n}$ is the $n$ th element in the vector $\mathbf{a}$ , $\mathcal{CN}(\boldsymbol{\mu},\mathbf{\Sigma})$ is a circularly symmetric complex Gaussian random vector with the mean $\boldsymbol{\mu}$ and covariance matrix $\mathbf{\Sigma}$ .

Refer to caption — Figure 1: System model of the SIM-aided holographic MIMO in cell-free networks.

II System Model

In this section, we describe our proposed SIM-aided HMIMO architecture, and then present the channel model of the cell-free network considered.

The system model of the SIM-aided HMIMO operating in narrowband cell-free networks is shown in Fig. 1¹¹1In this paper, we investigate the SIM-based hybrid beamforming design of narrowband cell-free networks. In wideband networks, both spatial-wideband effects and frequency-selective effects should be considered [35]. The SIM-based hybrid beamforming design of wideband networks considering the spatial-wideband effect and the frequency-selective effect is set aside for our future work.. In contrast to the conventional cellular network, where a central AP is deployed in each cell to support the surrounding UEs, the cell-free network deploys multiple distributed APs. We focus our attention on the uplink scenario, where we consider $L$ distributed multi-antenna APs and $K$ single-antenna UEs. AP- $l$ ( $l=1,2,\cdots,L$ ) has $M$ antennas and a stacked intelligent metasurface containing $T_{l}$ layers, with each layer composed of $N$ reconfigurable passive elements. Specifically, the information is sent from the $K$ UEs to $L$ APs. In each AP, the received RF signals go through the passive beamformer of the multi-layer SIM, and then they are converted to baseband signals via the RF chains. To recover the desired information, we rely on distributed operation in order to reduce the detection complexity, where each AP locally detects the data of the UEs supported by it, and then the CPU fuses the data detected at all APs via the fronthaul links to carry out the final detection of the desired information.

II-A SIM-Aided Holographic MIMO Architecture

Before describing the SIM-aided holographic MIMO architecture, we briefly review the conventional fully-digital massive MIMO architecture and the hybrid digital-analog hybrid massive MIMO architecture. In the conventional fully-digital massive MIMO architecture, the number of the RF chains is the same as that of the transmit antennas. To reduce the power consumption of the RF chains, a hybrid digital-analog beamforming architecture is formed, where only a few RF chains is employed, and a phase shift array (PSA) is harnessed between the RF chains and AP antennas.

The SIM-aided HMIMO architecture of cell-free networks is shown in Fig. 1. Firstly, the signals received at each AP go through an SIM-based beamformer. The output signals of the SIM are received by the antennas, and then converted to baseband signals via multiple RF chains. The SIM is constituted by a closed vacuum container having several layers of stacked reconfigurable metasurfaces, each of which is composed of a large number of passive densely-spaced elements [31], [32]. By appropriately configuring the phase shifts of the elements in each layer of the metasurface with the aid of a software controllers, such as a field programmable gate array (FPGA) [31], [32], a substantial beamforming gain may be achieved. In practice, the SIM is enclosed in a supporting structure surrounded by wave-absorbing materials, to prevent interferences from undesired diffraction, scattering, and environmental noise [36], [37]. The AP antennas are $M_{x}\times M_{y}$ uniform rectangular planar arrays (URPA), while the intelligent surfaces in each layer are $N_{x}\times N_{y}$ URPAs, satisfying $M=M_{x}M_{y}$ and $N=N_{x}N_{y}$ . We denote the size of each reconfigurable element as $\delta_{x}\times\delta_{y}$ , and the distances between the adjacent antennas are $d_{x}$ and $d_{y}$ along the $x$ and $y$ axis, respectively.

Explicitly, we contrast our SIM-aided HMIMO to various MIMO technologies in Table II.

TABLE II: Comparison of our SIM-aided HMIMO to various MIMO technologies.

	Beamforming structure	Number of RF chains	Number of PSA (or metasurface) elements	Hardware cost	Energy efficiency
Fully digital MIMO	Fully digital beamformer	Large	No	High	Low
\hdashlinePSA-aided hybrid MIMO	Digital beamformer + Analog beamformer relying on PSA	Moderate	Large	Moderate	Moderate
\hdashlineSIM-aided HMIMO in [31], [32]	Wave-based beamformer relying on SIM	No	Very large	Low	High
\hdashlineOur SIM-aided HMIMO	Digital beamformer + Wave-based beamformer relying on SIM	Moderate	Large	Low	High

II-B Channel Model

In this section, we describe the channel between the UEs and the AP antennas. As shown in Fig. 1, we have to consider the channel between the SIM and the antennas at the AP as well as the channel between the UEs and the SIM at the AP.

II-B1 Channel links between the SIM and the antennas at APs

We assume that the antenna array and the SIM layers are parallel to the $xoy$ plane of the Cartesian coordinate. We denote the Cartesian coordinates of the $M$ antennas at the $l$ th AP as $\mathbf{p}_{1}^{(l,0)}=(x_{1}^{(l,0)},y_{1}^{(l,0)},z^{(l,0)})$ , $\mathbf{p}_{2}^{(l,0)}=(x_{2}^{(l,0)},y_{2}^{(l,0)},z^{(l,0)}),\cdots,\mathbf{p}_{M}^{(l,0)}=(x_{M}^{(l,0)},y_{M}^{(l,0)},z^{(l,0)})$ , and that of the $N$ elements in the SIM layer- $t$ as $\mathbf{p}_{1}^{(l,t)}=(x_{1}^{(l,t)},y_{1}^{(l,t)},z^{(l,t)})$ , $\mathbf{p}_{2}^{(l,t)}=(x_{2}^{(l,t)},y_{2}^{(l,t)},z^{(l,t)}),\cdots,\mathbf{p}_{N}^{(l,t)}=(x_{N}^{(l,t)},y_{N}^{(l,t)},z^{(l,t)})$ with $t=1,2,\cdots,T_{l}$ . We employ the near-field model to describe the channel response between the antennas and the SIM layers.

As shown in Fig. 1, we denote the channel response between AP- $l$ of the SIM layer- $1$ and the antennas as $\mathbf{A}^{(l,1)}\in\mathbb{C}^{M\times N}$ , with the $(m,n)$ th entry $a_{m,n}^{(l,1)}$ representing the response of the $n$ th element in the SIM layer- $1$ to the $m$ th antenna, while $a_{m,n}^{(l,1)}$ is given by

\displaystyle a_{m,n}^{\left(l,1\right)}=\sqrt{\Gamma_{m,n}^{\left(l,1\right)}}\mathrm{e}^{-\jmath\frac{2\pi}{\lambda}\left\|\mathbf{p}_{n}^{\left(l,1\right)}-\mathbf{p}_{m}^{\left(l,0\right)}\right\|}.

(1)

In (1), $\lambda$ is the carrier wavelength, $\|\mathbf{p}_{n}^{(l,1)}-\mathbf{p}_{m}^{(l,0)}\|$ is the distance between the $n$ th element in the SIM layer- $1$ and the $m$ th antenna, while $\Gamma_{m,n}^{(l,1)}$ is the power radiated from the $n$ th element in the SIM layer- $1$ to the $m$ th antenna, which can be represented as [38]

\displaystyle\Gamma_{m,n}^{\left(l,1\right)}=\iint_{\mathcal{D}_{n}}\frac{\kappa\left(\frac{z^{\left(l,1\right)}-z^{\left(l,0\right)}}{\left\|\mathbf{p}_{n}^{\left(l,1\right)}-\mathbf{p}_{m}^{\left(l,0\right)}\right\|}\right)^{\frac{\kappa}{2}}}{4\pi\left\|\mathbf{p}_{n}^{\left(l,1\right)}-\mathbf{p}_{m}^{\left(l,0\right)}\right\|^{2}}\mathrm{d}x\mathrm{d}y,

(2)

where $\kappa$ is the directional radiation gain of the reconfigurable elements and the integration region is formulated as:

	$\displaystyle\mathcal{D}_{n}=$	$\displaystyle\left\{\left(x,y\right)\in\mathbb{R}^{2}:x_{n}-\frac{\delta_{x}}{2}\leq x\leq x_{n}+\frac{\delta_{x}}{2},\right.$
		$\displaystyle\left.y_{n}-\frac{\delta_{y}}{2}\leq y\leq y_{n}+\frac{\delta_{y}}{2}\right\}.$		(3)

Next, we focus our attention on the channel model between the SIM layers. For AP- $l$ we denote the response of the channel spanning from the SIM layer- $t$ to layer- $(t-1)$ as $\mathbf{A}^{(l,t)}\in\mathbb{C}^{N\times N}$ ( $t=2,3,\cdots,T_{l}$ ) with the $(n_{2},n_{1})$ th entry representing the response from the $n_{1}$ th element in the SIM layer- $t$ to the $n_{2}$ th element in the SIM layer- $(t-1)$ , given by

\displaystyle a_{n_{2},n_{1}}^{\left(l,t\right)}=\sqrt{\Gamma_{n_{2},n_{1}}^{\left(l,t\right)}}\mathrm{e}^{-\jmath\frac{2\pi}{\lambda}\left\|\mathbf{p}_{n_{1}}^{\left(l,t\right)}-\mathbf{p}_{n_{2}}^{\left(l,t-1\right)}\right\|}.

(4)

In (4), $\|\mathbf{p}_{n_{1}}^{(l,t)}-\mathbf{p}_{n_{2}}^{(l,t-1)}\|$ is the distance between the $n_{1}$ th element in the SIM layer- $t$ and the $n_{2}$ th element in the SIM layer- $(t-1)$ , while $\Gamma_{n_{2},n_{1}}^{(l,t)}$ is the power radiated from the $n_{1}$ th element in the SIM layer- $t$ to the $n_{2}$ th element in the SIM layer- $(t-1)$ , which can be represented as [18]

\displaystyle\Gamma_{n_{2},n_{1}}^{\left(l,t\right)}=\iint_{\mathcal{D}_{n}}\frac{\kappa\left(\frac{z^{\left(l,t\right)}-z^{\left(l,t-1\right)}}{\left\|\mathbf{p}_{n_{1}}^{\left(l,t\right)}-\mathbf{p}_{n_{2}}^{\left(l,t-1\right)}\right\|}\right)^{\frac{\kappa}{2}}}{4\pi\left\|\mathbf{p}_{n_{1}}^{\left(l,t\right)}-\mathbf{p}_{n_{2}}^{\left(l,t-1\right)}\right\|^{2}}\mathrm{d}x\mathrm{d}y.

(5)

Furthermore, the phase shifts of the reconfigurable elements in each layer can be adjusted to achieve a beamforming gain. In AP- $l$ , we denote the phase shift of the $n$ th element in layer- $t$ as $\theta_{n}^{(l,t)}$ . Upon considering the signal power attenuation resulting from the EM waves travelling through a SIM, we denote the radiation coefficient of the $n$ th element in layer- $t$ as $\upsilon_{n}^{(l,t)}$ satisfying $\upsilon_{n}^{(l,t)}\in[0,1]$ . Thus, by defining $\mathbf{\Upsilon}^{(l,t)}=\mathbf{Diag}\{\upsilon_{1}^{(l,t)},\upsilon_{2}^{(l,t)},\cdots,\upsilon_{N}^{(l,t)}\}$ and $\mathbf{\Theta}^{(l,t)}=\mathbf{Diag}\{\mathrm{e}^{\jmath\theta_{1}^{(l,t)}},\mathrm{e}^{\jmath\theta_{2}^{(l,t)}},\cdots,\mathrm{e}^{\jmath\theta_{N}^{(l,t)}}\}$ , the response matrix of the SIM’s layer- $t$ can be represented as

$\displaystyle\mathbf{\Xi}^{\left(l,t\right)}=$	$\displaystyle\sqrt{\mathbf{\Upsilon}^{\left(l,t\right)}}\mathbf{\Theta}^{\left(l,t\right)}$
$\displaystyle=$	$\displaystyle\mathbf{Diag}\left\{\sqrt{\upsilon_{1}^{\left(l,t\right)}}\mathrm{e}^{\jmath\theta_{1}^{\left(l,t\right)}},\sqrt{\upsilon_{2}^{\left(l,t\right)}}\mathrm{e}^{\jmath\theta_{2}^{\left(l,t\right)}},\cdots,\right.$
	$\displaystyle\left.\sqrt{\upsilon_{N}^{\left(l,t\right)}}\mathrm{e}^{\jmath\theta_{N}^{\left(l,t\right)}}\right\}.$	(6)

Therefore, the equivalent channel response of the SIM-based beamformer at the $l$ th AP, denoted as $\mathbf{G}^{(l)}\in\mathbb{C}^{M\times N}$ , is²²2In this paper, we assume that each SIM layer can be configured independently without mutual coupling. In reality, having mutual coupling among the metamaterial elements is inevitable due to their dense packaging. The holographic beamforming design considering the effect of mutual coupling is set aside for our future work.

$\displaystyle\mathbf{G}^{(l)}=$	$\displaystyle\mathbf{A}^{\left(l,1\right)}\mathbf{\Xi}^{\left(l,1\right)}\mathbf{A}^{\left(l,2\right)}\mathbf{\Xi}^{\left(l,2\right)}\cdots\mathbf{A}^{\left(l,T_{l}\right)}\mathbf{\Xi}^{\left(l,T_{l}\right)}$
$\displaystyle=$	$\displaystyle\mathbf{A}^{\left(l,1\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,1\right)}}\mathbf{\Theta}^{\left(l,1\right)}\mathbf{A}^{\left(l,2\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,2\right)}}\mathbf{\Theta}^{\left(l,2\right)}\cdots\mathbf{A}^{\left(l,T_{l}\right)}$
	$\displaystyle\sqrt{\mathbf{\Upsilon}^{\left(l,T_{l}\right)}}\mathbf{\Theta}^{\left(l,T_{l}\right)}.$	(7)

II-B2 Channel links between UEs and APs

We denote the large-scale fading and the small-scale fading between UE- $k$ and the SIM at AP- $l$ as $\varrho_{k}^{(l)}$ and $\mathbf{h}_{k}^{(l)}\in\mathbb{C}^{N\times 1}$ , respectively. We assume the knowledge of $\mathbf{h}_{1}^{(l)},\mathbf{h}_{2}^{(l)},\cdots,\mathbf{h}_{K}^{(l)}$ can be attained at the AP- $l$ . In practical metasurface-based holographic MIMO systems, this has to be acquired by CSI acquisition methods, such as the subspace-based channel estimator of [39] or the sparse channel estimator of [40]. We adopt the millimeter-wave (mmWave) channel model to characterize the propagation environment between the APs and each user³³3Note that here we consider a mmWave channel model as an example for describing the signal propagation between the APs and the UEs. But indeed, our proposed SIM-based beamforming architecture is also applicable to other channel models.. Specifically, the mmWave uplink channel of UE- $k$ is assumed to be the superposition of all propagation paths that are scattered in $\zeta_{c}$ clusters and each cluster contributes $\zeta_{p}$ paths, expressed as [41]

\displaystyle\mathbf{h}_{k}^{(l)}=\sqrt{\frac{1}{\zeta_{c}\zeta_{p}}}\sum_{c=1}^{\zeta_{c}}\sum_{p=1}^{\zeta_{p}}\alpha_{c,p}^{\left(l,k\right)}\mathbf{f}(\psi_{c,p}^{\left(l,k\right)},\varphi_{c,p}^{\left(l,k\right)}),

(8)

where $\alpha_{c,p}^{(l,k)}$ is the complex gain of the $p$ th path in the $c$ th cluster following $\alpha_{c,p}^{(l,k)}\sim\mathcal{CN}(0,1)$ , and $\mathbf{f}(\psi_{c,p}^{(l,k)},\varphi_{c,p}^{(l,k)})$ is

		$\displaystyle\mathbf{f}\left(\psi_{c,p}^{\left(l,k\right)},\varphi_{c,p}^{\left(l,k\right)}\right)=\left[1,\cdots,\right.$
		$\displaystyle\left.\mathrm{e}^{\jmath\frac{2\pi}{\lambda}\left(\delta_{x}n_{x}\sin\psi_{c,p}^{\left(l,k\right)}\cos\varphi_{c,p}^{\left(l,k\right)}+\delta_{y}n_{y}\sin\psi_{c,p}^{\left(l,k\right)}\sin\varphi_{c,p}^{\left(l,k\right)}\right)},\cdots\right.$
		$\displaystyle\left.\mathrm{e}^{\jmath\frac{2\pi}{\lambda}\left(\delta_{x}\left(N_{x}-1\right)\sin\psi_{c,p}^{\left(l,k\right)}\cos\varphi_{c,p}^{\left(l,k\right)}+\delta_{y}\left(N_{y}-1\right)\sin\psi_{c,p}^{\left(l,k\right)}\sin\varphi_{c,p}^{\left(l,k\right)}\right)}\right]^{\mathrm{H}},$		(9)

where $\psi_{c,p}^{(l,k)}$ and $\varphi_{c,p}^{(l,k)}$ are the elevation and azimuth angles of departure (AoD) from UE- $k$ to the SIM in AP- $l$ in the $p$ th path of the $c$ th cluster, respectively. Within the $c$ th cluster, the random variables $\psi_{c,p}^{(l,k)}$ and $\varphi_{c,p}^{(l,k)}$ have uniformly distributed mean values of $\mu_{\psi_{c}^{k}}$ and $\mu_{\varphi_{c}^{k}}$ , respectively, and have angular spreads, i.e. standard deviations, of $\sigma_{\psi_{c}^{k}}$ and $\sigma_{\varphi_{c}^{k}}$ , respectively.

The signal received by the antennas at AP- $l$ is given by

	$\displaystyle\mathbf{y}^{(l)}=$	$\displaystyle\sum_{k=1}^{K}\left(\sqrt{\rho_{k}\varepsilon_{u_{k}}\varepsilon_{\mathbf{v}^{(l)}}}\mathbf{q}_{k}^{(l)}s_{k}+\sqrt{\rho_{k}\left(1-\varepsilon_{u_{k}}\right)\varepsilon_{\mathbf{v}^{(l)}}}\mathbf{q}_{k}^{(l)}u_{k}\right.$
		$\displaystyle\left.+\sqrt{\rho_{k}(1-\varepsilon_{\mathbf{v}^{(l)}})}\mathbf{q}_{k}^{(l)}\odot\mathbf{v}_{k}^{(l)}\right)+\mathbf{w}^{(l)},$		(10)

where $\mathbf{q}_{k}^{(l)}=\sqrt{\varrho_{k}^{(l)}}\mathbf{G}^{(l)}\mathbf{h}_{k}^{(l)}$ is the equivalent channel spanning from the UE- $k$ to the antennas at AP- $l$ , $s_{k}$ is the desired information of UE- $k$ , $\rho_{k}$ denotes the transmit power of UE- $k$ , and $\mathbf{w}^{(l)}\sim\mathcal{CN}(\mathbf{0}_{M},\sigma_{w^{(l)}}^{2}\mathbf{I}_{M})$ is the additive noise at AP- $l$ . Furthermore, $u_{k}\sim\mathcal{CN}(0,1)$ represents the contamination of the information symbol $s_{k}$ due to HWIs at UE- $k$ , resulting from the power amplifier non-linearities, amplitude/phase imbalance in the In-phase/Quadrature mixers, phase noise in the local oscillator, sampling jitter and finite-resolution quantization in the analog-to-digital converters. Furthermore, $\mathbf{v}_{k}^{(l)}\sim\mathcal{CN}(\mathbf{0}_{M},\mathbf{I}_{M})$ is the distortion of the information symbol $s_{k}$ due to HWIs of the RF chains at AP- $l$ . Finally, $\varepsilon_{u_{k}}$ and $\varepsilon_{\mathbf{v}^{(l)}}$ represents the hardware quality factors of UE- $k$ and AP- $l$ satisfying $0\leq\varepsilon_{u_{k}}\leq 1$ and $0\leq\varepsilon_{\mathbf{v}^{(l)}}\leq 1$ , respectively [42, 43]. Explicitly, a hardware quality factor of 1 indicates that the hardware is ideal, while 0 means that the hardware is completely inadequate.

III Beamforming Design

As shown in the SIM-based cell-free network of Fig. 1, we rely on the distributed processing philosophy for reducing the required overhead of CSI sharing between APs. Specifically, for the information transmitted from the UE- $k$ , each AP carries out a local detection based on its received signal, where the SIM coefficient matrices and the RC vectors of each AP are optimized based on the corresponding local CSI. Then, the CPU gathers the locally detected signals from all APs to generate a final estimate of the UE- $k$ information. In the following, we first present the beamformer design in terms of the SIM coefficient matrices and combining vectors at each AP, followed by the optimization of the weight vectors at the CPU for generating the final detection of the UE information.

III-A SIM Coefficient Matrices and RC Vectors Design at APs

To recover the information transmitted from UE- $k$ , we denote the corresponding RC vector at AP- $l$ as $\mathbf{b}_{k}^{(l)}$ satisfying $\|\mathbf{b}_{k}^{(l)}\|^{2}=1$ . Thus, the locally recovered information of $s_{k}$ at AP- $l$ , denoted as $\hat{s}_{k}^{(l)}$ , is given by

$\displaystyle\hat{s}_{k}^{(l)}=$	$\displaystyle\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{y}^{(l)}$
$\displaystyle=$	$\displaystyle\sum_{k=1}^{K}\left(\sqrt{\rho_{k}\varepsilon_{u_{k}}\varepsilon_{\mathbf{v}^{(l)}}}\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{q}_{k}^{(l)}s_{k}\right.$
	$\displaystyle\left.+\sqrt{\rho_{k}\left(1-\varepsilon_{u_{k}}\right)\varepsilon_{\mathbf{v}^{(l)}}}\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{q}_{k}^{(l)}u_{k}\right.$
	$\displaystyle\left.+\sqrt{\rho_{k}\left(1-\varepsilon_{\mathbf{v}^{(l)}}\right)}\mathbf{b}_{k}^{(l)\mathrm{H}}\left(\mathbf{q}_{k}^{(l)}\odot\mathbf{v}_{k}^{(l)}\right)\right)+\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{w}^{(l)},$	(11)

and the corresponding signal-to-interference-plus-noise ratio (SINR), denoted as $\gamma_{k}^{(l)}$ , can be formulated as

\displaystyle\gamma_{k}^{(l)}=\frac{\rho_{k}\varepsilon_{u_{k}}\varepsilon_{\mathbf{v}^{(l)}}\left|\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{q}_{k}^{(l)}\right|^{2}}{\varsigma_{k,k}+\mathop{\sum\limits_{k^{\prime}=1}^{K}}\limits_{k^{\prime}\neq k}\left(\rho_{k^{\prime}}\varepsilon_{u_{k^{\prime}}}\varepsilon_{\mathbf{v}^{(l)}}\left|\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{q}_{k^{\prime}}^{(l)}\right|^{2}+\varsigma_{k,k^{\prime}}\right)+\sigma_{w^{(l)}}^{2}}.

(12)

In (12) $\varsigma_{k,{k^{\prime}}}$ represents the interference resulting from the distortion imposed on UE- $k^{\prime}$ signal for the recovery of $s_{k}$ due to the HWIs, given by

	$\displaystyle\varsigma_{k,{k^{\prime}}}=$	$\displaystyle\rho_{k^{\prime}}\left(1-\varepsilon_{u_{k^{\prime}}}\right)\varepsilon_{\mathbf{v}^{(l)}}\left\|\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{q}_{k^{\prime}}^{(l)}\right\|^{2}+\rho_{k^{\prime}}\left(1-\varepsilon_{\mathbf{v}^{(l)}}\right)\cdot$
		$\displaystyle\mathbf{b}_{k}^{(l)\mathrm{H}}\left(\left(\mathbf{q}_{k^{\prime}}^{(l)}\mathbf{q}_{k^{\prime}}^{(l)\mathrm{H}}\right)\odot\mathbf{I}_{M}\right)\mathbf{b}_{k}^{(l)}.$		(13)

For a specific UE- $k$ , we aim for jointly optimizing the active beamformers $\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)}$ and the passive SIM-based beamformers $\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ to maximize the SINR of $\gamma_{k}^{(l)}$ with $l=1,2,\cdots,L$ . The corresponding optimization problem can be formulated as

(P1)	$\displaystyle\max_{\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)},\mathbf{\Theta}^{\left(l,1\right)},\mathbf{\Theta}^{\left(l,2\right)},\cdots,\mathbf{\Theta}^{\left(l,T_{l}\right)}}\gamma_{k}^{(l)},\ l=1,2,\cdots,L$	(14)
s.t.	$\displaystyle\quad\mathbf{\Theta}^{\left(l,t\right)}\mathbf{\Theta}^{\left(l,t\right)\mathrm{H}}=\mathbf{I}_{N},\quad t=1,2,\cdots,T_{l},$	(15)
	$\displaystyle\quad\left\\|\mathbf{b}_{k^{\prime}}^{(l)}\right\\|^{2}=1,\quad k^{\prime}=1,2,\cdots,K.$	(16)

Since the $L$ APs jointly support all $K$ UEs, the beamforming of all APs focused on a specific UE- $k$ to maximize the SINR $\gamma_{k}^{(l)}$ ( $l=1,2,\cdots,L$ ) is carried out at the cost of disregarding the other $K-1$ UEs. The main advantage of the cell-free architecture is that of reducing the path-loss from each AP to its nearest UE. Thus, at each AP we aim for maximizing the channel gain between the AP and its nearest user within the SIM-based holographic beamformer, while the digital beamformer is optimized for all $K$ UEs to get their local information estimates. Specifically, we denote the AP set associated with the SIM-based beamforming focused on UE- $k$ as $\mathcal{L}_{k}=\{l:D_{k}^{(l)}\leq D_{k}^{(l^{\prime})},l^{\prime}=1,2,\cdots,L\}$ with $D_{k}^{(l)}$ being the distance from UE- $k$ to AP- $l$ . Therefore, problem (P1) can be formulated as

(P2)	$\displaystyle\max_{\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)},\mathbf{\Theta}^{\left(l,1\right)},\mathbf{\Theta}^{\left(l,2\right)},\cdots,\mathbf{\Theta}^{\left(l,T_{l}\right)}}\gamma_{k}^{(l)},\ l\in\mathcal{L}_{k}$	(17)
s.t.	$\displaystyle\quad\mathbf{\Theta}^{\left(l,t\right)}\mathbf{\Theta}^{\left(l,t\right)\mathrm{H}}=\mathbf{I}_{N},\quad t=1,2,\cdots,T_{l},$	(18)
	$\displaystyle\quad\left\\|\mathbf{b}_{k^{\prime}}^{(l)}\right\\|^{2}=1,\quad k^{\prime}=1,2,\cdots,K.$	(19)

Since (P2) is a non-convex problem, we can decouple it into a pair of sub-problems and optimize them iteratively as follows.

III-A1 Design of RC vectors at APs

Once the SIM-based beamformer $\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ is given, we can estimate the equivalent channels impinging from all $K$ UEs to AP- $l$ , i.e $\mathbf{q}_{1}^{(l)},\mathbf{q}_{2}^{(l)},\cdots,\mathbf{q}_{K}^{(l)}$ . Then, the optimal active beamformer $\mathbf{b}_{k^{\prime}}^{(l)}$ designed for the recovery of the information $s_{k^{\prime}}$ at the AP- $l$ can be optimized by using the MRC criterion as

\displaystyle\mathbf{b}_{k^{\prime}}^{(l)}=\frac{\mathbf{q}_{k^{\prime}}^{(l)}}{\left\|\mathbf{q}_{k^{\prime}}^{(l)}\right\|},\quad k^{\prime}=1,2,\cdots,K.

(20)

III-A2 Design of SIM coefficient matrices at APs

When the active beamformer $\mathbf{b}_{k}^{(l)}$ is given, the problem (P2) aims for optimizing the channel gain $\|\mathbf{G}^{(l)}\mathbf{h}_{k}^{(l)}\|^{2}$ , and it is given by

	(P3)	$\displaystyle\max_{\mathbf{\Theta}_{1}^{(l)},\mathbf{\Theta}_{2}^{(l)},\cdots,\mathbf{\Theta}_{T_{l}}^{(l)}}\left\\|\mathbf{G}^{(l)}\mathbf{h}_{k}^{(l)}\right\\|^{2},\ l\in\mathcal{L}_{k}$		(21)
	s.t.	$\displaystyle\quad\mathbf{\Theta}_{t}^{(l)}\mathbf{\Theta}_{t}^{(l)\mathrm{H}}=\mathbf{I}_{N},\quad t=1,2,\cdots,T_{l}.$		(22)

Since the sub-problem (P3) is still a non-convex one, we propose a layer-by-layer iterative optimization algorithm. Firstly, we optimize $\mathbf{\Theta}^{(l,1)}$ by fixing all the other $T_{l}-1$ layers of the SIM. Therefore, the channel gain can be represented as $\overline{\mathbf{b}}_{k}^{(l,1)\mathrm{H}}\mathbf{\Theta}^{(l,1)}\overline{\mathbf{h}}_{k}^{(l,1)}$ in conjunction with

\displaystyle\overline{\mathbf{b}}_{k}^{\left(l,1\right)\mathrm{H}}=\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{A}^{\left(l,1\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,1\right)}}

(23)

and

\displaystyle\overline{\mathbf{h}}_{k}^{\left(l,1\right)}=\mathbf{A}^{\left(l,2\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,2\right)}}\mathbf{\Theta}^{\left(l,2\right)}\cdots\mathbf{A}^{\left(l,T_{l}\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,T_{l}\right)}}\mathbf{\Theta}^{\left(l,T_{l}\right)}\mathbf{h}_{k}^{(l)},

(24)

and the passive beamformer $\mathbf{\Theta}^{(l,1)}$ can be optimized as

\displaystyle\mathbf{\Theta}^{\left(l,1\right)}=\mathbf{Diag}\left\{\mathrm{e}^{\jmath\left(\angle\overline{\mathbf{b}}_{k}^{\left(l,1\right)}-\angle\overline{\mathbf{h}}_{k}^{\left(l,1\right)}\right)}\right\}.

(25)

Afterwards, we optimize $\mathbf{\Theta}^{(l,2)}$ by fixing all the other $T_{l}-1$ layers of the SIM with the channel gain represented as $\overline{\mathbf{b}}_{k}^{(l,2)\mathrm{H}}\mathbf{\Theta}^{(l,2)}\overline{\mathbf{h}}_{k}^{(l,2)}$ along with

\displaystyle\overline{\mathbf{b}}_{k}^{\left(l,2\right)\mathrm{H}}=\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{A}^{\left(l,1\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,1\right)}}\mathbf{\Theta}^{\left(l,1\right)}\mathbf{A}^{\left(l,2\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,2\right)}}

(26)

and

\displaystyle\overline{\mathbf{h}}_{k}^{\left(l,2\right)}=\mathbf{A}^{\left(l,3\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,3\right)}}\mathbf{\Theta}^{\left(l,3\right)}\cdots\mathbf{A}^{\left(l,T_{l}\right)}\sqrt{\mathbf{\Upsilon}^{\left(l,T_{l}\right)}}\mathbf{\Theta}^{\left(l,T_{l}\right)}\mathbf{h}_{k}^{(l)},

(27)

and the passive beamformer $\mathbf{\Theta}^{(l,2)}$ can be optimized as

\displaystyle\mathbf{\Theta}^{\left(l,2\right)}=\mathbf{Diag}\left\{\mathrm{e}^{\jmath\left(\angle\overline{\mathbf{b}}_{k}^{\left(l,2\right)}-\angle\overline{\mathbf{h}}_{k}^{(l,2)}\right)}\right\}.

(28)

Then, $\mathbf{\Theta}^{(l,3)},\mathbf{\Theta}^{(l,4)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ can be optimized in turn.

The details of the layer-by-layer iterative optimization of the hybrid beamformer is presented in Algorithm 1.

Algorithm 1 Layer-by-layer iterative optimization algorithm for hybrid beamforming at AP-

l

(

l\in\mathcal{L}_{k}

0: The power attenuation matrices of SIM elements

\mathbf{\Upsilon}^{(l,1)},\mathbf{\Upsilon}^{(l,2)},\cdots,\mathbf{\Upsilon}^{(l,T_{l})}

, and the channel links

\mathbf{A}^{(l,1)},\mathbf{A}^{(l,2)},\cdots,\mathbf{A}^{(l,T_{l})}

and

\mathbf{h}_{k}^{(l)}

1: Set the random initial passive SIM-based beamforming matrices

\mathbf{\Theta}^{(l,1)}

\mathbf{\Theta}^{(l,2)}

\cdots

\mathbf{\Theta}^{(l,T_{l})}

, satisfying

\mathbf{\Theta}^{(l,1)}\mathbf{\Theta}^{(l,1)\mathrm{H}}=\mathbf{\Theta}^{(l,2)}\mathbf{\Theta}^{(l,2)\mathrm{H}}=\mathbf{\Theta}^{(l,T_{l})}\mathbf{\Theta}^{(l,T_{l})\mathrm{H}}=\mathbf{I}_{N}

2: repeat

3: The equivalent channel from the UE-

k

to AP-

l

antennas is

\mathbf{q}_{k}^{(l)}=\sqrt{\varrho_{k}^{(l)}}\mathbf{A}^{(l,1)}\sqrt{\mathbf{\Upsilon}^{(l,1)}}\mathbf{\Theta}^{(l,1)}\mathbf{A}^{(l,2)}\sqrt{\mathbf{\Upsilon}^{(l,2)}}\newline \mathbf{\Theta}^{(l,2)}\cdots\mathbf{A}^{(l,T_{l})}\sqrt{\mathbf{\Upsilon}^{(l,T_{l})}}\mathbf{\Theta}^{(l,T_{l})}\mathbf{h}_{k}^{(l)}

4: Design the active beamforming vector by the MRC method as

\mathbf{b}_{k^{\prime}}^{(l)}=\frac{\mathbf{q}_{k^{\prime}}^{(l)}}{\|\mathbf{q}_{k^{\prime}}^{(l)}\|}

for

k^{\prime}=1,2,\cdots,K

5: for

t=1

T_{l}

\overline{\mathbf{b}}_{k}^{(l,t)}=\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{A}^{(l,1)}\sqrt{\mathbf{\Upsilon}^{(l,1)}}\mathbf{\Theta}^{(l,1)}\cdots\mathbf{A}^{(l,t)}\sqrt{\mathbf{\Upsilon}^{(l,t)}}

\overline{\mathbf{h}}_{k}^{(l,t)}=\mathbf{A}^{(l,t+1)}\sqrt{\mathbf{\Upsilon}^{(l,t+1)}}\mathbf{\Theta}^{(l,t+1)}\cdots\mathbf{A}^{(l,T_{l})}\newline \sqrt{\mathbf{\Upsilon}^{(l,T_{l})}}\mathbf{\Theta}^{(l,T_{l})}\mathbf{h}_{k}^{(l)}

8: The optimal SIM-based passive beamforming matrix for layer-

t

\mathbf{\Theta}^{(l,t)}=\mathbf{Diag}\{{\mathrm{e}^{\jmath(\angle\overline{\mathbf{b}}_{k}^{(l,t)}-\angle\overline{\mathbf{h}}_{k}^{(l,t)})}}\}

9: end for

10: until reaching the iteration times.

10: The optimized active beamforming vectors

\mathbf{b}_{k}^{(l)}

and the optimized passive SIM-based beamforming matrices

\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}

III-B Weight Vector Design at the CPU

The CPU computes its estimate as a linear combination of the local estimates as [44]

\displaystyle\hat{s}_{k}=\sum_{l=1}^{L}\eta_{k}^{(l){\dagger}}\hat{s}_{k}^{(l)},

(29)

where $\boldsymbol{\eta}_{k}=[\eta_{k}^{(l)},\eta_{k}^{(2)},\cdots,\eta_{k}^{(L)}]^{\mathrm{T}}\in\mathbb{C}^{L\times 1}$ is the weight vector that the CPU assigns to the local signal estimate of the signal arriving from UE- $k$ satisfying $\|\boldsymbol{\eta}_{k}\|^{2}=1$ . Therefore, according to (III-A), $\hat{s}_{k}$ can be written as in (III-B).

$\displaystyle\hat{s}_{k}=$	$\displaystyle\sum_{l=1}^{L}\eta_{k}^{(l){\dagger}}\mathbf{b}_{k}^{(l)\mathrm{H}}\mathbf{y}^{(l)}$
$\displaystyle=$	$\displaystyle\underbrace{\sum_{l=1}^{L}\eta_{k}^{(l){\dagger}}\sqrt{\rho_{k}\varepsilon_{u_{k}}\varepsilon_{\mathbf{v}^{(l)}}}\left\\|\mathbf{q}_{k}^{(l)}\right\\|s_{k}}_{\text{Desired signal for $s_{k}$}}+\underbrace{\sum_{l=1}^{L}\eta_{k}^{(l){\dagger}}\sqrt{\rho_{k}\left(1-\varepsilon_{u_{k}}\right)\varepsilon_{\mathbf{v}^{(l)}}}\left\\|\mathbf{q}_{k}^{(l)}\right\\|u_{k}}_{\text{Signal distortion resulting from the HWI at UE-$k$}}+\underbrace{\sum_{l=1}^{L}\eta_{k}^{(l){\dagger}}\sqrt{\rho_{k}\left(1-\varepsilon_{\mathbf{v}^{(l)}}\right)}\frac{\mathbf{q}_{k}^{(l)\mathrm{H}}\left(\mathbf{q}_{k}^{(l)}\odot\mathbf{v}_{k}^{(l)}\right)}{\left\\|\mathbf{q}_{k}^{(l)}\right\\|}}_{\text{Signal distortion resulting from the HWI at APs}}$
	$\displaystyle+\underbrace{\mathop{\sum_{k^{\prime}=1}^{K}}_{k^{\prime}\neq k}\sum_{l=1}^{L}\eta_{k}^{(l){\dagger}}\Bigg{(}\sqrt{\rho_{k^{\prime}}\varepsilon_{u_{k^{\prime}}}\varepsilon_{\mathbf{v}^{(l)}}}\frac{\mathbf{q}_{k}^{(l)\mathrm{H}}\mathbf{q}_{k^{\prime}}^{(l)}s_{k^{\prime}}}{\left\\|\mathbf{q}_{k}^{(l)}\right\\|}+\sqrt{\rho_{k^{\prime}}(1-\varepsilon_{u_{k^{\prime}}})\varepsilon_{\mathbf{v}^{(l)}}}\frac{\mathbf{q}_{k}^{(l)\mathrm{H}}\mathbf{q}_{k^{\prime}}^{(l)}u_{k^{\prime}}}{\left\\|\mathbf{q}_{k}^{(l)}\right\\|}+\sqrt{\rho_{k^{\prime}}\left(1-\varepsilon_{\mathbf{v}^{(l)}}\right)}\frac{\mathbf{q}_{k}^{(l)\mathrm{H}}\left(\mathbf{q}_{k^{\prime}}^{(l)}\odot\mathbf{v}_{k^{\prime}}^{(l)}\right)}{\left\\|\mathbf{q}_{k}^{(l)}\right\\|}\Bigg{)}}_{\text{Inter-user interference}}$
	$\displaystyle+\underbrace{\sum_{l=1}^{L}\eta_{k}^{(l){\dagger}}\frac{\mathbf{q}_{k}^{(l)\mathrm{H}}\mathbf{w}^{(l)}}{\left\\|\mathbf{q}_{k}^{(l)}\right\\|}}_{\text{Additive noise}}$	(30)

The SINR of $\hat{s}_{k}$ , denoted as $\gamma_{k}$ , can be derived as

\displaystyle\gamma_{k}=\frac{\rho_{k}\left|\boldsymbol{\eta}_{k}^{\mathrm{H}}\mathbf{z}_{k,k}\right|^{2}}{\boldsymbol{\eta}_{k}^{\mathrm{H}}\mathbf{R}_{k}\boldsymbol{\eta}_{k}},

(31)

where $\mathbf{R}_{k}$ is given by

	$\displaystyle\mathbf{R}_{k}=$	$\displaystyle\rho_{k}\left(\mathbf{z}_{u_{k,k}}\mathbf{z}_{u_{k,k}}^{\mathrm{H}}+\left(\mathbf{z}_{\mathbf{v}_{k,k}}\mathbf{z}_{\mathbf{v}_{k,k}}^{\mathrm{H}}\right)\odot\mathbf{I}_{L}\right)+\mathop{\sum_{k^{\prime}=1}^{K}}_{k^{\prime}\neq k}\rho_{k^{\prime}}\left(\mathbf{z}_{k,k^{\prime}}\right.$
		$\displaystyle\left.\mathbf{z}_{k,k^{\prime}}^{\mathrm{H}}+\mathbf{z}_{u_{k,k^{\prime}}}\mathbf{z}_{u_{k,k^{\prime}}}^{\mathrm{H}}+\left(\mathbf{z}_{\mathbf{v}_{k,k^{\prime}}}\mathbf{z}_{\mathbf{v}_{k,k^{\prime}}}^{\mathrm{H}}\right)\odot\mathbf{I}_{L}\right)+\mathbf{W}$		(32)

in conjunction with

\displaystyle\mathbf{z}_{k,i}=\left[\begin{array}[]{c}\frac{\sqrt{\varepsilon_{u_{k}}\varepsilon_{\mathbf{v}^{(1)}}}}{\left\|\mathbf{q}_{k}^{(1)}\right\|}\mathbf{q}_{k}^{(1)\mathrm{H}}\mathbf{q}_{i}^{(1)}\\ \vdots\\ \frac{\sqrt{\varepsilon_{u_{k}}\varepsilon_{\mathbf{v}^{(L)}}}}{\left\|\mathbf{q}_{k}^{(L)}\right\|}\mathbf{q}_{k}^{(L)\mathrm{H}}\mathbf{q}_{i}^{(L)}\end{array}\right],

(36)

\displaystyle\mathbf{z}_{u_{k,i}}=\left[\begin{array}[]{c}\frac{\sqrt{\left(1-\varepsilon_{u_{k}}\right)\varepsilon_{\mathbf{v}^{(1)}}}}{\left\|\mathbf{q}_{k}^{(1)}\right\|}\mathbf{q}_{k}^{(1)\mathrm{H}}\mathbf{q}_{i}^{(1)}\\ \vdots\\ \frac{\sqrt{\left(1-\varepsilon_{u_{k}}\right)\varepsilon_{\mathbf{v}^{(L)}}}}{\left\|\mathbf{q}_{k}^{(L)}\right\|}\mathbf{q}_{k}^{(L)\mathrm{H}}\mathbf{q}_{i}^{(L)}\end{array}\right],

(40)

\displaystyle\mathbf{z}_{\mathbf{v}_{k,i}}=\left[\begin{array}[]{c}\frac{\sqrt{1-\varepsilon_{\mathbf{v}^{(1)}}}}{\left\|\mathbf{q}_{k}^{(L)}\right\|}\left\|\mathbf{q}_{k}^{(1){\dagger}}\odot\mathbf{q}_{i}^{(1)}\right\|\\ \vdots\\ \frac{\sqrt{1-\varepsilon_{\mathbf{v}^{(L)}}}}{\left\|\mathbf{q}_{k}^{(L)}\right\|}\left\|\mathbf{q}_{k}^{(L){\dagger}}\odot\mathbf{q}_{i}^{(L)}\right\|\end{array}\right],

(44)

and

\displaystyle\mathbf{W}=\mathbf{Diag}\{\sigma_{w^{(1)}}^{2},\sigma_{w^{(2)}}^{2},\cdots,\sigma_{w^{(L)}}^{2}\}.

(45)

Based on the generalized Rayleigh quotient, the maximum of $\gamma_{k}$ in (31) can be attained as follows:

$\displaystyle\gamma_{k}=$	$\displaystyle\rho_{k}\mathbf{z}_{k,k}^{\mathrm{H}}\mathbf{R}_{k}^{-1}\mathbf{z}_{k,k}$
$\displaystyle=$	$\displaystyle\rho_{k}\mathbf{z}_{k,k}^{\mathrm{H}}\Big{(}\rho_{k}\left(\mathbf{z}_{u_{k,k}}\mathbf{z}_{u_{k,k}}^{\mathrm{H}}+\left(\mathbf{z}_{\mathbf{v}_{k,k}}\mathbf{z}_{\mathbf{v}_{k,k}}^{\mathrm{H}}\right)\odot\mathbf{I}_{L}\right)$
	$\displaystyle+\mathop{\sum_{k^{\prime}=1}^{K}}_{k^{\prime}\neq k}\rho_{k^{\prime}}\left(\mathbf{z}_{k,k^{\prime}}\mathbf{z}_{k,k^{\prime}}^{\mathrm{H}}+\mathbf{z}_{u_{k,k^{\prime}}}\mathbf{z}_{u_{k,k^{\prime}}}^{\mathrm{H}}\right.$
	$\displaystyle\left.+\left(\mathbf{z}_{\mathbf{v}_{k,k^{\prime}}}\mathbf{z}_{\mathbf{v}_{k,k^{\prime}}}^{\mathrm{H}}\right)\odot\mathbf{I}_{L}\right)+\mathbf{W}\Big{)}^{-1}\mathbf{z}_{k,k},$	(46)

by setting

\displaystyle\boldsymbol{\eta}_{k}=\mathbf{R}_{k}^{-1}\mathbf{z}_{k,k}.

(47)

III-C Computational Complexity

III-C1 Computational complexity of the hybrid beamformer at each AP

The computational complexity of the proposed layer-by-layer iterative hybrid beamformer at AP- $l$ ( $l\in\mathcal{L}_{k}$ ) in Algorithm 1 depends both on the number of iterations in the alternating maximization, which is denoted as $\tau$ , and on the computational complexity required in each iteration. As for each iteration, the sub-problem of optimizing the combining vectors $\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)}$ is solved in line 3 and 4 of Algorithm 1, while the sub-problem of optimizing the coefficients of all layers in the SIM $\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ is evaluated in line 5 to 9 of Algorithm 1. The complexity of additions is neglected, since its operation is easily implemented in hardware. Hence, we quantify the computational complexity by counting the number of floating-point multiplication operations that are required. Specifically, the calculation of $\mathbf{q}_{k}^{(l)}$ in line 3 of Algorithm 1 requires

\displaystyle c_{1}^{\prime}=\left(T_{l}-1\right)\left(N^{2}+2N\right)+MN+2N+M

(48)

floating-point multiplication operations by exploiting the property that $\mathbf{\Upsilon}^{(l,1)},\mathbf{\Upsilon}^{(l,2)},\cdots,\mathbf{\Upsilon}^{(l,T_{l})}$ and $\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ are diagonal matrices. The calculation of $\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)}$ in line 4 of Algorithm 1 requires

\displaystyle c_{2}^{\prime}=2KM

(49)

floating-point multiplications. Furthermore, the calculation of $\overline{\mathbf{b}}_{k}^{(l,t)}$ and $\overline{\mathbf{h}}_{k}^{(l,t)}$ in line 6 and 7 of Algorithm 1 requires

\displaystyle c_{3,t1}^{\prime}=\left(t-1\right)(N^{2}+2N)+\left(M+1\right)N

(50)

and

\displaystyle c_{3,t2}^{\prime}=\left(T_{l}-t\right)(N^{2}+2N)

(51)

floating-point multiplications, respectively. The loop between line 8 of Algorithm 1 has no floating-point multiplications. Thus, the loop between line 5 to 9 of Algorithm 1 entails

	$\displaystyle c_{3}^{\prime}=$	$\displaystyle\sum\nolimits_{t=1}^{T_{L}}\left(c_{3,t1}^{\prime}+c_{3,t2}^{\prime}\right)$
	$\displaystyle=$	$\displaystyle\left(T_{l}^{2}-T_{l}\right)(N^{2}+2N)+T_{l}(M+1)N$		(52)

floating-point multiplications. Therefore, according to (48), (49) and (III-C1), the total number of floating-point multiplications in each iteration is

$\displaystyle c^{\prime}=$	$\displaystyle c_{1}^{\prime}+c_{2}^{\prime}+c_{3}^{\prime}$
$\displaystyle=$	$\displaystyle\left(T_{l}^{2}-1\right)N^{2}+\left(2T_{l}^{2}+T_{l}M+T_{L}+M\right)N$
	$\displaystyle+\left(2K+1\right)M.$	(53)

Hence for $N>M$ and $N>K$ , the overall computational complexity of our proposed layer-by-layer iterative optimization of the hybrid beamformer at all APs is

\displaystyle\mathcal{O}\left(\tau\left(\sum_{l=1}^{L}T_{l}^{2}\right)N^{2}\right).

(54)

This shows that the proposed layer-by-layer iterative optimization algorithm conceived for the SIM-based hybrid beamformer is of polynomial time-complexity with respect to the number of APs, the number of SIM layers at each AP, and the number of reconfigurable elements in each SIM layer.

III-C2 Computational complexity of the CPU processing

To recover the information $s_{k}$ , the computational complexity at the CPU includes the calculation of the matrix $\mathbf{R}_{k}$ defined in (III-B), of the linear combination vector $\boldsymbol{\eta}_{k}$ defined in (47) and of the information recovery $\hat{s}_{k}$ defined in (29). Specifically, the calculation of $\mathbf{R}_{k}$ requires

\displaystyle c_{1}^{\prime\prime}=\left(L^{2}+2L+3M\right)K

(55)

floating-point multiplication operations according to (III-B), (36), (40) and (44). The calculation of the linear combining vector $\boldsymbol{\eta}_{k}$ requires

\displaystyle c_{2}^{\prime\prime}=\frac{1}{3}\left(L^{3}-L\right)+L^{2}

(56)

floating-point multiplication operations by employing the Cholesky decomposition of the Hermitian positive-definite matrix $\mathbf{R}_{k}^{-1}$ . Furthermore, the calculation of the information recovery $\hat{s}_{k}$ requires

\displaystyle c_{3}^{\prime\prime}=L

(57)

floating-point multiplication operations according to (29). Therefore, according to (55), (56) and (57), the total number of floating-point multiplications required for recovering the information $s_{k}$ is

	$\displaystyle c^{\prime\prime}=$	$\displaystyle c_{1}^{\prime\prime}+c_{2}^{\prime\prime}+c_{3}^{\prime\prime}$
	$\displaystyle=$	$\displaystyle\frac{1}{3}L^{3}+\left(K+1\right)L^{2}+\left(2K+\frac{2}{3}\right)L+3MK.$		(58)

Therefore, the overall computational complexity of recovering the information $s_{1},s_{2},\cdots,s_{K}$ at the CPU is

\displaystyle\mathcal{O}\left(KL^{3}\right)+\mathcal{O}\left(K^{2}L^{2}\right)+\mathcal{O}\left(K^{2}M\right).

(59)

This shows that the CPU processing conceived for information recovery in the cell-free systems is of polynomial time-complexity with respect to the number of UEs, the number of APs and the number of RF chains at each AP.

III-D Convergence Analysis

First, in line 3 and 4 of Algorithm 1 convinced for the digital beamforming design, the received SINR $\gamma_{k}$ becomes non-decreasing after the digital beamformer $\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)}$ has been optimized, given the holographic beamformer $\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ , i.e.,

\displaystyle\gamma_{k}(\ddot{\mathbf{b}}^{(l,i+1)},\ddot{\mathbf{\Theta}}^{(l,i)})\geq\gamma_{k}(\ddot{\mathbf{b}}^{(l,i)},\ddot{\mathbf{\Theta}}^{(l,i)}),

(60)

where $\gamma_{k}(\ddot{\mathbf{b}}^{(l,i_{1})},\ddot{\mathbf{\Theta}}^{(l,i_{2})})$ represents the received SINR based on the digital beamformer $\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)}$ in the $i_{1}$ th iteration and on the holographic beamformer $\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ in the $i_{2}$ th iteration. Secondly, in line 5 to 9 of Algorithm 1, the received SINR $\gamma_{k}$ becomes non-decreasing after the holographic beamformer $\mathbf{\Theta}^{(l,1)},\mathbf{\Theta}^{(l,2)},\cdots,\mathbf{\Theta}^{(l,T_{l})}$ has been optimized, given the digital beamformer $\mathbf{b}_{1}^{(l)},\mathbf{b}_{2}^{(l)},\cdots,\mathbf{b}_{K}^{(l)}$ , i.e.,

\displaystyle\gamma_{k}(\ddot{\mathbf{b}}^{(l,i)},\ddot{\mathbf{\Theta}}^{(l,i+1)})\geq\gamma_{k}(\ddot{\mathbf{b}}^{(l,i)},\ddot{\mathbf{\Theta}}^{(l,i)}).

(61)

Therefore, (60) and (61) imply that in each iteration of the proposed layer-by-layer iterative optimization algorithm at each AP, the objective function value of the received SINR $\gamma_{k}$ is non-decreasing. Additionally, the objective function value sequence obtained throughout the iteration process is monotonic and it is also bounded, hence the overall algorithm is guaranteed to converge.

TABLE III: Simulation Parameters.

Parameters	Values
Carrier frequency	$f_{c}=30$ GHz
\hdashlineNumber of APs	$L=16$
\hdashlineNumber of UEs	$K=8$
\hdashlineNumber of AP antennas	$M=4\times 4$
\hdashlineNumber of SIM elements in each layer	$N=16\times 16$
\hdashlineNumber of SIM layers	$T_{l}=4$ ( $l=1,2,\cdots L$ )
\hdashlineInter-layer distance	$z^{(l,0)}=z^{(l,1)}=\cdots=z^{(l,T_{l}-1)}=\lambda$ ( $l=1,2,\cdots L$ )
\hdashlineDirectional radiation gain	$\kappa=10\text{dB}$
\hdashlineDistance between AP antennas	$d_{x}=d_{y}=\frac{\lambda}{2}$
\hdashlineSize of the SIM elements	$\delta_{x}=\delta_{y}=\frac{\lambda}{4}$
\hdashlineNoise power	$\sigma_{w^{(l)}}^{2}=-80\text{dBm}$ ( $l=1,2,\cdots L$ )
\hdashlineHardware quality factor	$\varepsilon_{u_{k}}=\varepsilon_{\mathbf{v}^{(l)}}=\varepsilon=1$ ( $k=1,2,\cdots,K$ , and $l=1,2,\cdots,L$ )
\hdashlineTransmit power at UEs	$\rho=\rho_{1}=\rho_{2}=\cdots=\rho_{K}$
\hdashlinePath loss exponent	$\beta=3.5$
\hdashlinePath loss at the reference distance of 1 meter	$\mathrm{C}_{0}=10^{-3}$
\hdashlineNumber of clusters	$\zeta_{c}=4$ ( $c=1,2,\cdots,\zeta_{c}$ )
\hdashlineNumber of paths in each cluster	$\zeta_{p}=8$ ( $p=1,2,\cdots,\zeta_{p}$ )
\hdashlineMean values of the elevation angles	$\mu_{\psi_{c}^{k}}$ are randomly distributed in $[0^{\circ},180^{\circ}]$ ( $c=1,2,\cdots,\zeta_{c}$ , and $k=1,2,\cdots,K$ )
\hdashlineMean values of the azimuth angles	$\mu_{\varphi_{c}^{k}}$ are randomly distributed in $[0^{\circ},360^{\circ}]$ ( $c=1,2,\cdots,\zeta_{c}$ , and $k=1,2,\cdots,K$ )
\hdashlineElevation angle spreads	$\sigma_{\psi_{c}^{k}}=7.5^{\circ}$ ( $c=1,2,\cdots,\zeta_{c}$ , and $k=1,2,\cdots,K$ )
\hdashlineAzimuth angle spreads	$\sigma_{\varphi_{c}^{k}}=7.5^{\circ}$ ( $c=1,2,\cdots,\zeta_{c}$ , and $k=1,2,\cdots,K$ )
\hdashlineSIM power radiation coefficients	$\upsilon_{n}^{(l,t)}=1$ ( $n=1,2,\cdots,N$ , $t=1,2,\cdots,T_{l}$ , and $l=1,2,\cdots,L$ )

IV Numerical and Simulation Results

In this section, the average achievable rate of the SIM-based holographic MIMO of cell-free networks is quantified. A total of $L$ APs are uniformly distributed in the area $\mathcal{S}$ with the Cartesian coordinate of $\{(x,y):-100\text{m}\leq x\leq 100\text{m},-100\text{m}\leq y\leq 100\text{m}\}$ . Furthermore, $K$ UEs are uniform-randomly distributed in the area $\mathcal{S}$ . We employ a distance-dependent path-loss model for the UE-AP channel, given by $\varrho_{k}^{(l)}=\min\{\mathrm{C}_{0},\mathrm{C}_{0}(d_{k}^{(l)})^{-\beta}\}$ , where $d_{k}^{(l)}$ is the length of the link spanning from AP- $l$ to UE- $k$ , $\beta$ is the path loss exponent, and $\mathrm{C}_{0}$ is the path loss at the reference distance of 1 meter [45]. Referring to [32], [46], the simulation parameters are those given in Table III, unless otherwise specified.

Fig. 2 compares the average achievable rate $\overline{R}=\frac{1}{K}\sum_{k=1}^{K}R_{k}=\frac{1}{K}\sum_{k=1}^{K}\log_{2}(1+\gamma_{k})$ versus the transmit power $\rho$ for different number of APs in cell-free networks, with the number of AP $L=16$ and $L=64$ , respectively. Observe from Fig. 2 that the average achievable rate can be improved by increasing the number of APs when the average distance between the APs and UEs can be reduced. When the hardware is non-ideal, i.e. $\varepsilon<1$ , the average achievable rate tends to a constant value upon increasing of the transmit power $\rho$ . Furthermore, as seen in Fig. 2, an obvious performance degradation exists when the HWIs are ignored during information recovery, especially with the increase of the number of APs. Fig. 3 compares the average achievable rate $\overline{R}$ versus the transmit power $\rho$ for different number of SIM layers in each AP, where the number $N$ of elements in each SIM layer is set to $N=16\times 16$ , $N=8\times 8$ and $N=4\times 4$ when the number of SIM layers $T_{l}$ are 1, 2 and 4 respectively for ensuring that the same total number of SIM elements are employed. Observe in Fig. 3 that the average achievable rate is improved upon increasing the number of SIM layers. It inspires our future research on the SIM design to determine how many SIM layers is optimal in terms of maximizing the performance gain given the total number of SIM elements.

To provide further insights, Fig. 4 characterizes the average achievable rate of both our SIM-based hybrid beamforming architecture and of the conventional cell-free full-digital beamforming architecture, for the hardware quality factors of $\varepsilon=1$ , $1-10^{-4}$ and $1-10^{-2}$ , respectively. In terms of the perfect hardware quality factor, i.e., $\varepsilon=1$ , Fig. 4 (a) shows that the SIM-based hybrid beamformer has higher average achievable rate than the conventional full-digital beamformer, even when the SIM has the signal radiation attenuation associated with the power radiation coefficients of $\upsilon_{n}^{(l,t)}=0.5$ . This is a benefit of the beamforming gain attained by the configuration of the SIM elements. Furthermore, the average achievable rate can be improved upon increasing the number of SIM layers. When the hardware quality is imperfect, i.e., $\varepsilon<1$ , it is illustrated in Fig. 4 (b) and Fig. 4 (c) that the SIM-based hybrid beamformer outperforms the full-digital beamformer in the low-SNR region, while the achievable rate in the high-SNR region is limited by the hardware quality. Moreover, in the high-SNR region, increasing the number of SIM layers even degrades the achievable rate when the SIM power radiation coefficients $\upsilon_{n}^{(l,t)}<1$ , since the signals are attenuated in each layer of the SIM and it cannot be compensated by the configuration of SIM elements due to the limitation of the hardware quality. It can be observed that although the SIM can increase the equivalent channel gain, but it cannot compensate for the deleterious effects of HWIs in the high-SNR region.

The performance comparison of the average achievable rate $\overline{R}$ versus the number of antennas $M$ at each AP is presented in Fig. 5, for different hardware quality factors $\varepsilon$ . It shows that the average achievable rate can be improved upon increasing the number of antennas at each AP, but only at the cost of higher energy consumption. However, the signal distortion resulting from the imperfect hardware quality can be compensated by employing more AP antennas.

Fig. 6 portrays the average achievable rate versus the inter-layer distance, where we have $z=z^{(l,0)}=z^{(l,1)}=\cdots=z^{(l,T_{l}-1)}$ for $l=1,2,\cdots,L$ . The figure shows that there exists an optimal inter-layer distance in terms of the highest average achievable rate. Furthermore, it also shows that the optimal inter-layer distance increases upon employing more intelligent surface elements. Observe that when the number of intelligent surface elements is small, the beamforming gain attained by the SIM is predominantly limited by the path loss between layers, and thus a higher average achievable rate can be attained by reducing the inter-layer distance. By contrast, upon increasing the number of intelligent surface elements, the inter-layer path loss effect becomes negligible and increasing the inter-layer distance improves the signal radiation between adjacent layers, increasing the achievable rate.

To investigate the impact of the signal attenuation caused by the signal travelling through each layer of the SIM on the achievable rate performance, Fig. 7 presents the average achievable rate $\overline{R}$ versus the number of SIM layers at each AP, characterize by different power radiation coefficients $\upsilon_{n}^{(l,t)}$ . It shows that although the average achievable rate degrades with the reduction of the radiation coefficients $\upsilon_{n}^{(l,t)}$ , this effect can be compensated by increasing the number of SIM layers.

In Fig. 8 we portray the average achievable rate $\overline{R}$ versus the number of iterations $\tau$ attained by the proposed layer-by-layer iterative optimization algorithm used by the hybrid beamformer. Observe that although the convergence speed is reduced as the number of SIM layers increases, it achieves convergence within 10 iterations.

In the above simulations, we assume that the reconfigurable SIM elements have continuous phase shifts in the range of $[0,2\pi)$ . However, in practical hardware implementations, the phase shift of each reconfigurable element is limited to a finite number of discrete values. For simplicity, we assume that the discrete phase shift set is obtained by uniformly quantizing the interval $[0,2\pi)$ into $2^{b}$ levels, with $b$ respecting the number of phase shift control bits. Thus, we have $\theta_{n}^{(l,t)}\in\{2\pi\cdot\frac{0}{2^{B}},2\pi\cdot\frac{1}{2^{B}},\cdots,2\pi\cdot\frac{2^{b}-1}{2^{b}}\}$ . Fig. 9 characterizes the achievable sum-rate, denoted as $R_{\mathrm{sum}}=\sum_{k=1}^{K}R_{k}=\sum_{k=1}^{K}\log_{2}(1+\gamma_{k})$ versus the number of phase shift control bits $b$ , with different number of UEs $k$ . It shows that 4-bit phase shift quantization can approach the rate of the infinite phase shift resolution.

V Conclusions

We conceived an uplink SIM-based cell-free HMIMO architecture, where the distributed signal processing was employed. Each AP carries out a local detection of the UE information, where the hybrid beamformer and the RC vectors of each distributed AP are optimized based on our low-complexity layer-by-layer iterative optimization algorithm. The CPU recovers the final UE information by fusing the local detections gleaned from all APs, where the RC weight vector used for combining the local detections is designed based on the MMSE criterion, taking into account the HWIs of the RF chains at the APs and of the UEs. The simulation results showed that the achievable rate of the SIM-based cell-free HMIMO network improves upon increasing the number of SIM layers as well as the number of elements in each layer. Furthermore, due to the HWI of the RF chains at the APs and the UEs, the achievable rate saturates in the high-SNR region.

References

[1] P. Wang, J. Fang, L. Dai, and H. Li, “Joint transceiver and large intelligent surface design for massive MIMO mmWave systems,” IEEE Trans. Wireless Commun., vol. 20, no. 2, pp. 1052–1064, 2020.
[2] L. Dai, J. Tan, Z. Chen, and H. V. Poor, “Delay-phase precoding for wideband THz massive MIMO,” IEEE Trans. Wireless Commun., vol. 21, no. 9, pp. 7271–7286, 2022.
[3] J. Du, S. Ye, L. Jin, X. Li, H. Q. Ngo, and O. A. Dobre, “Tensor-based joint channel estimation for multi-way massive MIMO hybrid relay systems,” IEEE Trans. Veh. Technol., vol. 71, no. 9, pp. 9571–9585, 2022.
[4] X. Li, M. Zhang, H. Chen, C. Han, L. Li, D.-T. Do, S. Mumtaz, and A. Nallanathan, “UAV-enabled multi-pair massive MIMO-NOMA relay systems with low-resolution ADCs/DACs,” IEEE Trans. Veh. Technol., vol. 73, no. 2, pp. 2171–2186, 2024.
[5] J. Gao, C. Zhong, G. Y. Li, J. B. Soriaga, and A. Behboodi, “Deep learning-based channel estimation for wideband hybrid mmWave massive MIMO,” IEEE Trans. Commun., vol. 71, no. 6, pp. 3679–3693, 2023.
[6] C.-J. Wang, C.-K. Wen, S. Jin, and S.-H. Tsai, “Finite-alphabet precoding for massive MU-MIMO with low-resolution DACs,” IEEE Trans. Wireless Commun., vol. 17, no. 7, pp. 4706–4720, 2018.
[7] L. Yang, J. Yang, W. Xie, M. O. Hasna, T. Tsiftsis, and M. Di Renzo, “Secrecy performance analysis of RIS-aided wireless communication systems,” IEEE Trans. Veh. Technol., vol. 69, no. 10, pp. 12 296–12 300, 2020.
[8] L. Yang, F. Meng, J. Zhang, M. O. Hasna, and M. Di Renzo, “On the performance of RIS-assisted dual-hop UAV communication systems,” IEEE Trans. Veh. Technol., vol. 69, no. 9, pp. 10 385–10 390, 2020.
[9] L. V. Nguyen, A. L. Swindlehurst, and D. H. Nguyen, “Linear and deep neural network-based receivers for massive MIMO systems with one-bit ADCs,” IEEE Trans. Wireless Commun., vol. 20, no. 11, pp. 7333–7345, 2021.
[10] Q. Li, M. El-Hajjar, I. Hemadeh, A. Shojaeifard, A. A. Mourad, B. Clerckx, and L. Hanzo, “Reconfigurable intelligent surfaces relying on non-diagonal phase shift matrices,” IEEE Trans. Veh. Technol., vol. 71, no. 6, pp. 6367–6383, 2022.
[11] C. Huang, S. Hu, G. C. Alexandropoulos, A. Zappone, C. Yuen, R. Zhang, M. Di Renzo, and M. Debbah, “Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends,” IEEE Wireless Commun., vol. 27, no. 5, pp. 118–125, 2020.
[12] Q. Li, M. El-Hajjar, I. Hemadeh, D. Jagyasi, A. Shojaeifard, E. Basar, and L. Hanzo, “The reconfigurable intelligent surface-aided multi-node IoT downlink: Beamforming design and performance analysis,” IEEE Internet Things J., vol. 10, no. 7, pp. 6400–6414, 2022.
[13] Q. Li, M. El-Hajjar, I. Hemadeh, A. Shojaeifard, A. A. Mourad, and L. Hanzo, “Reconfigurable intelligent surface aided amplitude-and phase-modulated downlink transmission,” IEEE Trans. Veh. Technol., vol. 72, no. 6, pp. 8146–8151, 2023.
[14] I. Yoo and D. R. Smith, “Sub-6-Ghz uplink massive MIMO system using holographic beamforming metasurfaces: A conceptual development,” IEEE Wireless Commun. Lett., vol. 12, no. 4, pp. 644–648, 2023.
[15] R. Deng, Y. Zhang, H. Zhang, B. Di, H. Zhang, and L. Song, “Reconfigurable holographic surface: A new paradigm to implement holographic radio,” IEEE Veh. Technol. Mag., vol. 18, no. 1, pp. 20–28, 2023.
[16] T. Gong, P. Gavriilidis, R. Ji, C. Huang, G. C. Alexandropoulos, L. Wei, Z. Zhang, M. Debbah, H. V. Poor, and C. Yuen, “Holographic MIMO communications: Theoretical foundations, enabling technologies, and future directions,” IEEE Commun. Surv. Tutor., vol. 26, no. 1, pp. 196–257, 2024.
[17] J. An, C. Yuen, C. Huang, M. Debbah, H. V. Poor, and L. Hanzo, “A tutorial on holographic MIMO communication–Part I: Channel modeling and channel estimation,” IEEE Commun. Lett., vol. 27, no. 7, pp. 1664–1668, 2023.
[18] S. Zeng, H. Zhang, B. Di, H. Qin, X. Su, and L. Song, “Reconfigurable refractive surfaces: An energy-efficient way to holographic MIMO,” IEEE Commun. Lett., vol. 26, no. 10, pp. 2490–2494, 2022.
[19] R. Deng, B. Di, H. Zhang, Y. Tan, and L. Song, “Reconfigurable holographic surface: Holographic beamforming for metasurface-aided wireless communications,” IEEE Trans. Veh. Technol., vol. 70, no. 6, pp. 6255–6259, 2021.
[20] R. Deng, B. Di, H. Zhang, and L. Song, “HDMA: Holographic-pattern division multiple access,” IEEE J. Sel. Areas Commun., vol. 40, no. 4, pp. 1317–1332, 2022.
[21] R. Deng, B. Di, H. Zhang, D. Niyato, Z. Han, H. V. Poor, and L. Song, “Reconfigurable holographic surfaces for future wireless communications,” IEEE Wireless Commun., vol. 28, no. 6, pp. 126–131, 2021.
[22] X. Hu, R. Deng, B. Di, H. Zhang, and L. Song, “Holographic beamforming for ultra massive MIMO with limited radiation amplitudes: How many quantized bits do we need?” IEEE Commun. Lett., vol. 26, no. 6, pp. 1403–1407, 2022.
[23] R. Deng, B. Di, H. Zhang, H. V. Poor, and L. Song, “Holographic MIMO for LEO satellite communications aided by reconfigurable holographic surfaces,” IEEE J. Sel. Areas Commun., vol. 40, no. 10, pp. 3071–3085, 2022.
[24] X. Hu, R. Deng, B. Di, H. Zhang, and L. Song, “Holographic beamforming for LEO satellites,” IEEE Commun. Lett., vol. 27, no. 10, pp. 2717–2721, 2023.
[25] L. Wei, C. Huang, G. C. Alexandropoulos, E. Wei, Z. Zhang, M. Debbah, and C. Yuen, “Multi-user holographic MIMO surfaces: Channel modeling and spectral efficiency analysis,” IEEE J. Sel. Top. Signal Process., vol. 16, no. 5, pp. 1112–1124, 2022.
[26] H. Wu, Y. Chen, Y. Ming, and Z. Wang, “Two-timescale beamforming optimization for downlink multi-user holographic MIMO surfaces,” IEEE Trans. Veh. Technol., vol. 73, no. 3, pp. 4476–4481, 2023.
[27] N. Shlezinger, O. Dicker, Y. C. Eldar, I. Yoo, M. F. Imani, and D. R. Smith, “Dynamic metasurface antennas for uplink massive MIMO systems,” IEEE Trans. Commun., vol. 67, no. 10, pp. 6829–6843, 2019.
[28] H. Wang, N. Shlezinger, S. Jin, Y. C. Eldar, I. Yoo, M. F. Imani, and D. R. Smith, “Dynamic metasurface antennas based downlink massive MIMO systems,” in 2019 IEEE 20th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC). IEEE, 2019, pp. 1–5.
[29] L. You, J. Xu, G. C. Alexandropoulos, J. Wang, W. Wang, and X. Gao, “Energy efficiency maximization of massive MIMO communications with dynamic metasurface antennas,” IEEE Trans. Wireless Commun., vol. 22, no. 1, pp. 393–407, 2022.
[30] Y. Li, S. Gong, H. Liu, C. Xing, N. Zhao, and X. Wang, “Near-field beamforming optimization for holographic XL-MIMO multiuser systems,” IEEE Trans. Commun., vol. 72, no. 4, pp. 2309–2323, 2024.
[31] J. An, C. Xu, D. W. K. Ng, G. C. Alexandropoulos, C. Huang, C. Yuen, and L. Hanzo, “Stacked intelligent metasurfaces for efficient holographic MIMO communications in 6G,” IEEE J. Sel. Areas Commun., vol. 41, no. 8, pp. 2380–2396, 2023.
[32] J. An, C. Yuen, C. Xu, H. Li, D. W. K. Ng, M. Di Renzo, M. Debbah, and L. Hanzo, “Stacked intelligent metasurface-aided MIMO transceiver design,” IEEE Wireless Commun., 2024.
[33] S. Chen, J. Zhang, E. Björnson, J. Zhang, and B. Ai, “Structured massive access for scalable cell-free massive MIMO systems,” IEEE J. Sel. Areas Commun., vol. 39, no. 4, pp. 1086–1100, 2020.
[34] W. Xu, J. An, H. Li, L. Gan, and C. Yuen, “Algorithm unrolling-based distributed optimization for RIS-assisted cell-free networks,” IEEE Internet Things J., vol. 11, no. 1, pp. 944–957, 2024.
[35] J. Xu, L. You, G. C. Alexandropoulos, X. Yi, W. Wang, and X. Gao, “Near-field wideband extremely large-scale MIMO transmissions with holographic metasurface-based antenna arrays,” IEEE Trans. Wireless Commun., 2024.
[36] C. Liu, Q. Ma, Z. J. Luo, Q. R. Hong, Q. Xiao, H. C. Zhang, L. Miao, W. M. Yu, Q. Cheng, L. Li et al., “A programmable diffractive deep neural network based on a digital-coding metasurface array,” Nature Electronics, vol. 5, no. 2, pp. 113–122, 2022.
[37] J. An, C. Yuen, Y. L. Guan, M. Di Renzo, M. Debbah, H. V. Poor, and L. Hanzo, “Two-dimensional direction-of-arrival estimation using stacked intelligent metasurfaces,” arXiv preprint arXiv:2402.08224, 2024.
[38] H. Lu and Y. Zeng, “Communicating with extremely large-scale array/surface: Unified modeling and performance analysis,” IEEE Trans. Wireless Commun., vol. 21, no. 6, pp. 4039–4053, 2021.
[39] Ö. T. Demir, E. Björnson, and L. Sanguinetti, “Channel modeling and channel estimation for holographic massive MIMO with planar arrays,” IEEE Wireless Commun. Lett., vol. 11, no. 5, pp. 997–1001, 2022.
[40] M. Cui and L. Dai, “Channel estimation for extremely large-scale MIMO: Far-field or near-field?” IEEE Trans. Commun., vol. 70, no. 4, pp. 2663–2677, 2022.
[41] I. A. Hemadeh, K. Satyanarayana, M. El-Hajjar, and L. Hanzo, “Millimeter-wave communications: Physical channel models, design considerations, antenna constructions, and link-budget,” IEEE Commun. Surv. Tutor., vol. 20, no. 2, pp. 870–913, 2017.
[42] E. Björnson, J. Hoydis, L. Sanguinetti et al., “Massive MIMO networks: Spectral, energy, and hardware efficiency,” Foundations and Trends® in Signal Processing, vol. 11, no. 3-4, pp. 154–655, 2017.
[43] Q. Li, M. El-Hajjar, Y. Sun, I. Hemadeh, A. Shojaeifard, Y. Liu, and L. Hanzo, “Achievable rate analysis of the STAR-RIS aided NOMA uplink in the face of imperfect CSI and hardware impairments,” IEEE Trans. Commun., vol. 71, no. 10, pp. 6100–6114, 2023.
[44] Ö. T. Demir, E. Björnson, L. Sanguinetti et al., “Foundations of user-centric cell-free massive MIMO,” Foundations and Trends® in Signal Processing, vol. 14, no. 3-4, pp. 162–472, 2021.
[45] M. Haenggi, Stochastic geometry for wireless networks. Cambridge University Press, 2012.
[46] R. Deng, B. Di, H. Zhang, Y. Tan, and L. Song, “Reconfigurable holographic surface-enabled multi-user wireless communications: Amplitude-controlled holographic beamforming,” IEEE Trans. Wireless Commun., vol. 21, no. 8, pp. 6003–6017, 2022.