# Could SBND-PRISM probe Lepton Flavor Violation?

Gustavo F. S. Alves<EMAIL_ADDRESS>Renata Zukanovich Funchal
<EMAIL_ADDRESS>Departamento de Física Matemática, Instituto de Física
Universidade de São Paulo, 05315-970 São Paulo, Brasil  Pedro A. N. Machado
<EMAIL_ADDRESS>Particle Theory Department, Fermi National Accelerator
Laboratory, Batavia, IL 60510, USA

###### Abstract

We investigate the possibility of using the Short-Baseline Near Detector
(SBND) at Fermilab to constrain lepton flavor violating decays of pions and
kaons. We study how to leverage SBND-PRISM, the use of the neutrino beam
angular spread to mitigate systematic uncertainties, to enhance this analysis.
We show that SBND-PRISM can put stringent limits on the flavor violating
branching ratios $\rm{BR}(\pi^{+}\to\mu^{+}\nu_{e})=8.9\times 10^{-4}$,
$\rm{BR}(K^{+}\to\mu^{+}\nu_{e})=3.2\times 10^{-3}$, improving previous
constraints by factors 9 and 1.25, respectively. We also estimate the SBND-
PRISM sensitivity to lepton number violating decays,
$\rm{BR}(\pi^{+}\to\mu^{+}\overline{\nu}_{e})=2.1\times 10^{-3}$ and
$\rm{BR}(K^{+}\to\mu^{+}\overline{\nu}_{e})=7.4\times 10^{-3}$, though not
reaching previous BEBC limits. Last, we identify several ways how the SBND
collaboration could improve this analysis.111The work and conclusions
presented in this publication are not to be considered as results from the
SBND collaboration.

Neutrino Physics

††preprint: FERMILAB-PUB-24-0119-T

##  I Introduction

Neutrino flavor oscillations represent to this day the sole evidence of lepton
flavor violation (LFV). They also may be connected, if neutrinos are massive
Majorana fermions, to lepton number violation (LNV). The lack of any
conclusive signals from the Large Hadron Collider experiments triggers the
interest on low energy precision measurements. Besides the classic searches
for LFV processes involving charged leptons such as $\mu\to e\gamma$, $\mu\to
3e$ or $\mu\to e$ – conversion in nuclei, and for the quintessential LNV
neutrinoless double-$\beta$ decay, meson decays can also play a role in
scrutinizing these symmetries.

In fact, light pseudoscalar meson leptonic decays, the main sources of the
conventional neutrino fluxes in accelerator neutrino oscillation experiments,
have already been used to this purpose in the past. The best bounds come from
the Big European Bubble Chamber (BEBC) [1, 2, 3], and are, correspondingly,
$\rm{BR}(\pi^{+}\to\mu^{+}\nu_{e})<8.0\times 10^{-3}$,
$\rm{BR}(K^{+}\to\mu^{+}\nu_{e})<4\times 10^{-3}$ for the lepton number
conserving (LNC) and $\rm{BR}(\pi^{+}\to\mu^{+}\overline{\nu}_{e})<1.5\times
10^{-3}$, $\rm{BR}(K^{+}\to\mu^{+}\overline{\nu}_{e})<3.3\times 10^{-3}$ for
the LNV modes. Moreover, the combination of very intense beams and highly
sensitive detectors, now available at some particle accelerator facilities,
allows to use them once again as important probes of these Beyond the Standard
Model (BSM) manifestations of new physics.

From the theoretical point of view LNV and LFV contributions to rare meson
decays are expected to be connected to UV complete models with new degrees of
freedom appearing above the electroweak scale. A useful way to tackle the
problem is by means of the effective field theory approach. One constructs an
effective Lagrangian consisting of the SM part and a series of all possible
non-renormalizable operators of dimension $d\geq 5$ invariant under
$SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ involving the SM fields

${\cal L}_{\rm eff}={\cal L}_{\rm SM}+\sum_{i}\sum_{d\geq
5}\frac{c_{i}^{(d)}}{(\Lambda^{(d)})^{d-4}}{\cal O}_{i}^{d}\,,$ (1)

where $c_{i}^{(d)}$ is the Wilson coefficient of operator
$\mathcal{O}_{i}^{d}$, which will take it to be 1 for simplicity, and
$\Lambda^{(d)}$ is the BSM scale. The lowest dimension of these operators
($d=5$) is the well-known LNV Weinberg operator [4] that lead to Majorana
neutrino masses. There are four (two vector and two scalar) $d=6$ quark-lepton
charged current operators that may induce LFV, either conserving or violating
lepton number, in light meson leptonic decays [5]. Allowing for a single of
these operators to saturate the experimental bound one can probe a BSM scale
$\Lambda^{(6)}$ up to $(0.6-4.7)$ TeV using
$\pi^{+}\to\mu^{+}\nu_{e}(\bar{\nu}_{e})$ or up to $(1.3-7.7)$ TeV using
$K^{+}\to\mu^{+}\nu_{e}(\bar{\nu}_{e})$. In Ref. [6] a classification of
$\Delta L=2$ operators with $d<12$ have been provided. Using this
classification the authors of Ref. [7] have investigated constraints on LNV
operators using rare kaon decays. In particular, they have shown that the
current experimental limits on $\pi^{+}(K^{+})\to\mu^{+}\bar{\nu}_{e}$ can
discard $\Lambda^{(7)}\lesssim 1.9\,(2.4)$ TeV.

There have been several analyses on LFV and LNV connected to hadronic decays
using the EFT formalism; for example Refs. [8, 9] studied the impact on LFV of
adding right-handed neutrinos to the standard model EFT, and Refs. [10, 11]
estimated the beta and pion decays. Some of the operators discussed above can
also be realized in several UV complete models, such as supersymmetric
R-parity violating models [12] or leptoquark models [13, 14, 15, 7].
Additionally, Refs. [16, 17] have discussed how UV complete scenarios of LFV
can affect muon decays and alleviate the LSND anomaly. For a comprehensive
study on how to relate neutrino mass mechanisms and LFV interactions, see Ref.
[18]. Note that the Lepton Flavor Universality observables defined for
pseudoscalar mesons decays as

$R^{M}_{e/\mu}\equiv\frac{\Gamma(M^{+}\to e^{+}+\rm
invisible)}{\Gamma(M^{+}\to\mu^{+}+\rm invisible)}\,,\quad\quad M=\pi,K\;,$
(2)

allow to obtain from data and SM calculations the double ratios
$R^{\pi}_{e/\mu}[\rm exp]/R^{\pi}_{e/\mu}[\rm SM]=0.9994\pm 0.0024$ and
$R^{K}_{e/\mu}[\rm exp]/R^{K}_{e/\mu}[\rm SM]=1.0044\pm 0.0041$ [19]. The
experimental ratios are, however, not very sensitive to the type of LNV or LFV
processes we have discussed above as the final state neutral
particle/antiparticle nature and flavor are not determined.

The Short Baseline Near Detector (SBND) is an ideal candidate to search for
these type of processes. It will collect unprecedented statistics of neutrino
events and set the beam content by probing the unoscillated neutrino flux.
This suppresses natural contamination of $\nu_{e}$ from neutrino oscillations,
mitigating background for the flavor violating meson decays searches 222Note
that the baseline of SBND (110 m) is much shorter than the shortest neutrino
oscillation length $L_{\rm osc}^{\rm atm}=(4\pi E/\Delta m^{2}_{31})\sim
2\times 10^{3}\;\rm km$, so there is no contribution from neutrino oscillation
to the flux. . Before proceed with more details, we can make a simple estimate
of the sensitivity of SBND in constraining light meson decays to the wrong
neutrino. If we consider $M\to\mu\nu_{e}$, we can compare the BSM $\nu_{e}$
flux with the uncertainty in the intrinsic $\nu_{e}$ contamination, namely,

$\text{BR}(M^{+}\to\mu^{+}\nu_{e})\phi_{\nu_{\mu}}^{M}=\sigma_{e}\phi_{\nu_{e}},$
(3)

where $\phi_{\nu_{\alpha}}^{M}$ is the beam flux of $\nu_{\alpha}$ coming from
the decays of meons $M$ and $\sigma_{e}$ is the uncertainty in the electron
neutrino flux for SBND. For a significant part of the neutrino spectrum at
SBND $\phi_{\nu_{\mu}}\sim\mathcal{O}(100)\ \phi_{\nu_{e}}$. We take a
$\sigma_{\phi}\sim 10\%$ uncertainty on the $\nu_{e}$ contamination [21]. This
means that a $\text{BR}(M\to\mu\nu_{e})\sim\mathcal{O}(10^{-3})$ could be
achieved, demonstrating that SBND can potentially provide competitive
constraints on LFV pion decays.

Building from this simple argument, we also note that the SBND collaboration
has identified a potential use angular spread of the Booster Neutrino Beam
–the SBND-PRISM technique [22], similar to previous studies [23, 24]. SBND can
measure neutrino events in off-axis positions by identifying the vertex of the
event. Since the distance to the target is about 110 meters, compared to SBND
largest dimension of 5 meters, the nontrivial angular spread translates into
specific $\nu_{\mu}$ and $\nu_{e}$ energy spectra that are off-axis angle
dependent. Therefore, as the ratio of $\nu_{\mu}$ to $\nu_{e}$ events change
with the off-axis angle, we can leverage that as an additional handle to
control background uncertainties. In this work we aim to show how SBND can
improve the searches for these exotic $\pi$ and $K$ decays.

In section II we discuss the details of the SBND experiment and the PRISM
technique that are pertinent for our study. We discuss how we simulate the
signal and our statistical procedure in section III, followed by our results
in section IV. To conclude, section V summarizes our key findings. We also
discuss potential ways to improve the experimental sensitivity beyond what we
have considered.

##  II The Short Baseline Neutrino Program

The Booster Neutrino Beam (BNB) line at Fermilab provides the neutrino flux
for the Short Baseline Neutrino (SBN) program [25]. This is achieved by
colliding a beam of 8 GeV protons with a beryllium target, producing a
secondary beam of hadrons. In sequence, magnetic horns select the hadron
charge allowing it to function in two different modes: the neutrino mode,
where the positively charged hadrons are focused and the antineutrino mode
where negatively charged hadrons are focused instead. The selected flux of
particles propagate in the 50 m long decay pipe giving rise to the neutrino
beam through their subsequent decays. The SBN program is expected to run
mostly in the neutrino mode where the flux is mainly comprised of $\nu_{\mu}$
(96%) originated predominantly from $\pi^{+}\to\mu^{+}\nu_{\mu}$ for energies
up to approximately $2\ \rm GeV$, beyond which $K^{+}\to\mu^{+}\nu_{\mu}$
decays become the dominant source of $\nu_{\mu}$ production. The BNB flux
contains a small fraction of $\overline{\nu}_{\mu}$ ($5\%$), stemming from the
contamination of wrong-sign hadrons and muon decays. Last, the $\nu_{e}$ and
$\overline{\nu}_{e}$ beam contamination of about 0.5% originates from both
kaon and muon decays.

The SBND detector is the first experiment in the BNB line, located 110 m away
from the BNB target. It is a 112 ton liquid argon time projection chamber
(LArTPC) with an active volume of 5.0 m (L) $\times$ 4.0 m (W) $\times$ 4.0 m
(H). With an initially planned exposure of $6.6\times 10^{20}$ protons on
target (POT), SBND will collect about 7 million $\nu_{\mu}$ charged current
(CC) interactions and 50 thousand $\nu_{e}$ CC interactions, in a three years
period. The LArTPC technology can provide a 3D reconstruction of neutrino-
argon interactions events by reconstructing the tracks of charged particles as
they travel throught the detector ionizing the medium.

The term SBND-PRISM denotes the experiment’s capability to conduct
measurements at off-axis positions. The detector can be logically segmented
into eight slices, each corresponding to an increasing off-axis angle within
the range of $[0^{\circ},1.6^{\circ}]$. The kinematics of two- and three-body
meson decays dictate that the energy spectra of $\nu_{\mu}$ coming from
$\pi\to\mu\nu$ tend to become narrower and shifts to lower energies as we move
off-axis. In contrast, the $\nu_{e}$ from kaons and muons three-body decays
are significantly less affected by the off-axis angle. This distinction
implies a layer-dependent profile in the distribution of $\nu_{\mu}$ and
$\nu_{e}$ events [22], which can be leveraged to mitigate systematic
uncertainties, in particular those related to neutrino-nucleus cross sections.
In Fig. 1 we present the fraction of events for $\nu_{\mu}$ from $\pi$ decays,
$\nu_{\mu}$ from $K$ decays, and intrinsic $\nu_{e}$ expected at each layer.
In the next section we will describe in detail how we extracted these
fractions.

Figure 1: Fraction of CC events for $\nu_{\mu}$ from $\pi$ decay (blue),
$\nu_{\mu}$ from $K$ decay (gray) and intrinsic $\nu_{e}$ (salmon) expected at
each layer labeled L1 to L8 of SBND-PRISM. The insets pictured above each
fraction represent the front view of the SBND detector, segmented in eight
layers, in which the corresponding layer is highlighted in color.

##  III Analysis

In order to study the SBND capability to improve the constraints on LFV meson
decays, we rely on the distribution of simulated events and fluxes publicly
provided by the SBND collaboration [22], for $6.6\times 10^{20}$ POT. This
information allows us to know the total count of $\nu_{\mu}$ and $\nu_{e}$ CC
events as well as their energy spectrum at SBND for each off-axis layer.
However, it lacks the parent particle information, crucial for our purposes.
In the Booster Neutrino Beam, muon neutrinos below $E_{\rm flip}$ of a few GeV
are dominantly produced by pion decays, while above $E_{\rm flip}$ neutrino
production via kaon decay dominates [26]. This flip between parent particle is
fairly sharp. Therefore, to estimate the corresponding $E_{\rm flip}^{i}$ for
each layer $i$, we do the following.

We simulate an 8 GeV proton beam with PYTHIA [27], generating the meson flux.
To mock the magnetic horn effect, we align all positive mesons in the beam
direction, and decay them to neutrinos. First, we select the neutrinos that
cross the front face of SBND, and calculate at which energy the neutrino flux
from pion decay cross the flux from kaon decays. This yields $E_{\rm
flip}=2.35$ GeV in comparison with the official BNB simulation $E_{\rm
flip}^{\rm official}\simeq 2.3$ GeV [26], which serves to validate this
approach as a decent zero order approximation. Then we repeat this procedure
applying angular cuts to select the neutrinos that cross each layer $i$,
obtaining $E_{\rm flip}^{i}$. For simplicity, we use $E_{\rm flip}^{i}$ as a
sharp cut to define the parent particle, that is, we approximate that all muon
neutrinos below $E_{\rm flip}^{i}$ come from pion decays, while above that
energy come from kaon decays. Since we will not explicitly use spectral
information, we do not expect the sharp cut to impact our analysis in any
significant way.

With that information at hand, we can estimate our signal coming from the
flavor violating pions and kaon decays. Since the two-body pion and kaon decay
is dominated by the branching ratio into $\mu\nu_{\mu}$, we can obtain our
signal flux $\phi_{\nu_{e},{\rm sig}}^{i}$ by multiplying the flux of muon
neutrinos coming from the decays of meson $M$ by the nonstandard branching,
that is,

$\phi_{\nu_{e},{\rm sig}}^{i}=\phi_{\nu_{\mu}}^{M,i}\,{\rm
BR}(M\to\mu\nu_{e}),$ (4)

where $i$ denotes the SBND-PRISM layer. The key of the impact of SBND-PRISM
lies in the different SBND-PRISM layer dependence of the signal and background
fluxes, as the signal originates in two-body decays but the intrinsic
$\nu_{e}$ background comes from three-body decays. To make this point more
clear, we present in fig. 2 the number of events for signal (gray) and
background (dark green) in each layer for two representative values of the
nonstandard LFV decays: ${\rm BR}(\pi^{+}\to\mu^{+}\nu_{e})=2\times 10^{-3}$
(left panel) and ${\rm BR}(K^{+}\to\mu^{+}\nu_{e})=2\times 10^{-2}$ (right
panel). For reference, we show a 10% fully correlated systematic error in the
prediction, as it will be discussed later. We also display in the lower panels
the signal-to-background ratio $S/B$ for these benchmarks. We see clearly that
while for nonstandard pion decays $S/B$ decreases with opening angle, while
for kaon decays the opposite happens. By leveraging this geometrical
dependence of $S/B$, SBND-PRISM will be less affected by the large cross
section uncertainties on the background.

Figure 2: Number of signal (gray) and background(dark green) $\nu_{e}$ events,
in each SBND-PRISM layer for two representative branching ratios of LFV pion
(left panel) and kaon (right panel) decays. For reference, we show a 10% fully
correlated uncertainty error bar. In the lower panels, we show the signal-to-
background ratio $S/B$.

Now, we outline our statistical approach. As we lack an understanding of
systematic uncertainties that can lead to spectral distortions on neutrino
energy distributions, and their correlations, we will do simple cut and count
analyses. We focus in two strategies: utilizing the total event rate at SBND
and applying the SBND-PRISM concept. The SBND sensitivity is evaluated through
a $\chi^{2}$ test,

$\chi_{\rm SBND}^{2}=\frac{N_{\nu_{e},\rm
sig}^{2}}{N_{\nu_{e}}+\sigma_{\nu_{e}}^{2}N_{\nu_{e}}^{2}}$ (5)

where $N_{\nu_{e},\rm sig}$ is the number of signal events, while
$N_{\nu_{e}}$ is the number of intrinsic $\nu_{e}$ background events. We take
the systematic uncertainty on the background events to be
$\sigma_{\nu_{e}}=10\%$ [26]. The number of signal events coming from
nonstandard decays of a meson $M=\pi,\,K$ has been calculated as

$N_{\nu_{e},\rm sig}=N_{\nu_{\mu}}^{M}\,{\rm BR}(M^{+}\to\mu^{+}\nu_{e}),$ (6)

where $N_{\nu_{\mu}}^{M}$ is the standard prediction for the number of
$\nu_{\mu}$ charged current events coming from parent meson $M$, as described
in the previous section. Note that this implicitly neglects differences in the
CC cross sections for $\nu_{e}$ and $\nu_{\mu}$. We will also estimate SBND’s
sensitivity to LNV decays, namely $\pi^{+}\to\mu^{+}\overline{\nu}_{e}$ and
$K^{+}\to\mu^{+}\overline{\nu}_{e}$. To that end we substitute the branching
ratio appropriately in the definition of the signal, and we rescale the number
of signal events by the ratio of $\overline{\nu}_{e}$ to $\nu_{e}$ cross
sections. Backgrounds are the same since LArTPCs cannot distinguish electrons
from positrons.

For SBND-PRISM, we build a covariance matrix to correlate event rate in
different layers as

$C_{ij}=\left[(\sigma_{\nu_{e}}^{\rm c})^{2}+\delta_{i,j}(\sigma^{\rm
u}_{\nu_{e}})^{2}\right]N^{i}_{\nu_{e}}N^{j}_{\nu_{e}}+\delta_{i,j}N^{i}_{\nu_{e}},$
(7)

where $i,j$ denote the SBND-PRISM layer, $\sigma_{\nu_{e}}^{\rm c}$ and
$\sigma^{\rm u}_{\nu_{e}}$ are the correlated and uncorrelated systematic
uncertainties across layers, and $N_{\nu_{e}}^{i}$ is the number of intrinsic
$\nu_{e}$ background events in layer $i$. We take the correlated uncertainty
to be $\sigma_{\nu_{e}}^{\rm c}=10\%$, consistently with the previous
analysis, while we study two possibilities for the uncorrelated systematics,
$\sigma^{\rm u}_{\nu_{e}}=5\%$ or $2\%$, corresponding to a pessimistic and an
optimistic case. We then build a $\chi^{2}$ as

$\chi^{2}_{\rm PRISM}=\sum_{i,j=1}^{8}N_{\nu_{e},\rm
sig}^{i}(C^{-1})_{ij}N_{\nu_{e},\rm sig}^{j},$ (8)

where $i,j$ denote the SBND-PRISM layer and the number of signal events per
layer is defined akin to what was described previously. For both analyses we
will calculate the 95% C.L. on the nonstandard branching ratio using
$\chi^{2}=3.841$.

##  IV Results and discussion

We proceed to present our results for the SBND and SBND-PRISM sensitivity on
nonstandard pion and kaon decays. Starting with the lepton flavor violating
decays, $\pi^{+}\to\mu^{+}\nu_{e}$ and $K^{+}\to\mu^{+}\nu_{e}$, we present in
the left panel of fig. 3 and in table 1 the current BEBC constraints [1, 2,
3], as well as the projected SBND sensitivities for an exposure of $10^{21}$
POT [28] for three scenarios: SBND only with 10% uncertainty; SBND-PRISM with
10% (correlated) and 5% (uncorrelated) uncertainties; and SBND-PRISM with 10%
(correlated) and 2% (uncorrelated) uncertainties. For reference, we also show
the statistics only limit. To facilitate the interpretation of this result in
terms of new physics, we also show in table 1 a range of constraints on the
scale of dimensions-6 operators that can lead to such nonstandard transitions.

Figure 3: Sensitivity of SBND and SBND-PRISM to $\pi^{+}\to\mu^{+}\nu_{e}$ and
$K^{+}\to\mu^{+}\nu_{e}$ lepton flavor violating decays. For comparison, the
BEBC [1, 2, 3] results and the statistics only SBND-PRISM limits are shown.

We can see that, even without PRISM, SBND can improve the BEBC constraint on
$\pi^{+}\to\mu^{+}\nu_{e}$ by a factor $\sim 5$; while when employing the
SBND-PRISM technique an improvement of $5-9$ can be achieved. We can
understand these improvements by noting a few differences between BEBC and
SBND. First, the number of background CC $\nu_{e}$ events in BEBC is a factor
two larger than in SBND; namely, $(\nu_{e}:\nu_{\mu})^{\rm BEBC}=(30/1800)\sim
0.016$ [1], compared to $(\nu_{e}:\nu_{\mu})^{\rm SBND}=(3.3\times
10^{4}:4\times 10^{6})\sim 0.008$ [26] for SBND. Besides that, BEBC had a much
smaller statistical sample of events: $6060\pm 440$ candidates $\nu_{\mu}$ CC
events versus $\mathcal{O}(10^{7})$ in SBND. SBND-PRISM presents further
improvement due to the change in the fractions of $\nu_{e}:\nu_{\mu}$ across
different layers as discussed in section III, see fig. 2.

For kaons, the situation is slightly different, as can be seen from the right
panel of fig. 3 and table 1. Since BEBC was on the CERN SPS beamline, which
accelerates protons to 450 GeV, the kaon production was significantly higher
than in SBND which is in the 8 GeV Booster Neutrino Beam. With a lower signal
rate, the systematic uncertainties in the $\nu_{e}$ background become dominant
in SBND, which limit the sensitivity to nonstandard kaon decays. Nevertheless,
SBND-PRISM can lead to improved constraints over BEBC, for the aggressive
scenario with 2% bin-to-bin uncorrelated systematics, by leveraging the larger
asymmetry of muon neutrinos from kaons across each layer in comparison to the
$\nu_{e}$ background. This reduces the impact of systematic uncertainties,
leading to a constraint that is a factor of $1.25$ better than BEBC’s.

For the lepton number violating decays, $\pi^{+}\to\mu^{+}\overline{\nu}_{e}$
and $K^{+}\to\mu^{+}\overline{\nu}_{e}$, we summarize our results in table 1.
For both pion and kaon LNV decays, our findings show that SBND does not
improve the BEBC constraints, even in the optimistic SBND-PRISM scenario.
Although SBND will have higher statistics and smaller $\nu_{e}$ contamination
than BEBC, the latter has an important advantage for this type of search.
BEBC, being a bubble chamber, can disentangle the $\nu_{e}$ and
$\overline{\nu}_{e}$ backgrounds. The $\overline{\nu}_{e}$ contamination there
is found to be below the $0.1\%$, leading to an essentially background free
analysis.

Let us now discuss possible improvements to our proposed analysis. First, we
have not considered any spectral information due to the fact that we do not
have access to realistic systematic uncertainties affecting the spectrum,
particularly in the SBND-PRISM case where the correlations among spectra in
different layers are relevant. Nevertheless, neutrinos coming from two-body
$\pi$ and $K$ decays have very distinctive spectral features, the former being
lower energy than the latter. This features can be used to further improve
background rejection, and most likely would lead to further improvement of all
constraints presented here. A handle that could be used to enhance the
sensitivity to LNV decays is the event topology. It has been observed that
neutrino CC events tend to have more visible protons than antineutrino ones
[29]. Zero-proton events are expected to have a better signal-to-background
ratio than one-proton events, due to the smaller $\overline{\nu}_{e}$
contamination, which could ameliorate the LNV constraints in table 1. Last, a
more direct way of distinguishing $\nu_{e}$ from $\overline{\nu}_{e}$ would be
to magnetize liquid argon time projection chambers [30]. This could also make
feasible doing the LNV analysis when running in antineutrino mode.

Lepton Flavor Violating Decays
---
Experiment | $\rm{BR}(\pi^{+}\to\mu^{+}\nu_{e})$ | $\rm{BR}(K^{+}\to\mu^{+}\nu_{e})$ |  $\Lambda^{(6)}$(TeV)
BEBC | $8\times 10^{-3}$ | $4\times 10^{-3}$ | 0.59 - 4.9
SBND ($10\%$) | $1.5\times 10^{-3}$ | $1.7\times 10^{-2}$ | 0.89 - 4.6
SBND-PRISM ($10\%$, $5\%$) | $1.2\times 10^{-3}$ | $5.8\times 10^{-3}$ | 0.94 - 5.8
SBND-PRISM ($10\%$, $2\%$) | $8.9\times 10^{-4}$ | $3.2\times 10^{-3}$ | 1 - 6.8
Statistics only | $8.5\times 10^{-5}$ | $8.6\times 10^{-4}$ | 1.8 - 9.4
Lepton Number Violating Decays
Experiment | $\rm{BR}(\pi^{+}\to\mu^{+}\overline{\nu}_{e})$ | $\rm{BR}(K^{+}\to\mu^{+}\overline{\nu}_{e})$ |  $\Lambda^{(6)}$ (TeV)
BEBC | $1.5\times 10^{-3}$ | $3.3\times 10^{-3}$ | 0.89 - 6.7
SBND ($10\%$) | $4.0\times 10^{-3}$ | $3.9\times 10^{-2}$ | 0.7 - 3.6
SBND-PRISM ($10\%$,$5\%$) | $3.1\times 10^{-3}$ | $1.3\times 10^{-2}$ | 0.74 - 4.8
SBND-PRISM ($10\%$,$2\%$) | $2.1\times 10^{-3}$ | $7.4\times 10^{-3}$ | 0.82 - 5.5
Statistics only | $2.1\times 10^{-4}$ | $1.9\times 10^{-3}$ | 1.5 - 7.7
Table 1: Summary of the results for the pion and kaon LFV (upper part) and LNV
(lower part) searches considering SBND and SBND-PRISM. The values in
parenthesis indicate the uncertainties (total flux for SBND, and correlated
and uncorrelated for SBND-PRISM, respectively). For comparison we also include
the previous bounds from BEBC [1, 2, 3] and the SBND reach considering only
statistical errors. We also include the range constraints that can be set on
the scale of dimension-6 operators that can lead to those transitions.

##  V Conclusions

In this work we have estimated how the Short-Baseline Near Detector can
improve constraints on lepton flavor violating decays of pions and kaons. In
doing so, we have also shown how to leverage the Booster Neutrino Beam angular
spread to mitigate the impact of systematic uncertainties –a technique known
as SBND-PRISM [22]. We find that SBND-PRISM can constrain ${\rm
BR}(\pi^{+}\to\mu^{+}\nu_{e})<(0.89-1.2)\times 10^{-3}$ and ${\rm
BR}(K^{+}\to\mu^{+}\nu_{e})<(3.2-5.8)\times 10^{-3}$, for aggressive and
conservative assumptions on systematic uncertainties. This represents
improvement factors on these nonstandard decays between 6 to 9 for pion
nonstandard decays and up to 1.25 for kaons with respect to previous bounds
from the Big European Bubble Chamber. We have also estimated SBND’s
sensitivity to lepton number violating decays
$\pi^{+}(K^{+})\to\mu^{+}\overline{\nu}_{e}$. However, we have not found any
improvement over previous BEBC constraints due larger backgrounds in SBND. We
have identified several ways this analysis may be further improved by the SBND
collaboration, including leveraging spectral information and event topology.
The interesting possibility of magnetizing LArTPCs would also bolster the
experiment sensitivity to LNV decays.

## Acknowledgement

We thank Gustavo Burdman, Marco del Tutto, Ornella Palamara and Jure Zupan for
invaluable discussions. G.F.S.A. is fully financially supported by Fundação de
Amparo à Pesquisa do Estado de São Paulo (FAPESP) under Contracts No.
2022/10894-8 and No. 2020/08096-0. G.F.S.A. would like to thank the
hospitality of the Fermilab Theory Group. R. Z. F. is partially supported by
FAPESP under Contract No. 2019/04837-9, and by Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq). P.M. is supported by Fermi
Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S.
Department of Energy, Office of Science, Office of High Energy Physics.

## References

  * Cooper-Sarkar  _et al._ [1982] A. M. Cooper-Sarkar, J. G. Guy, A. G. Michette, M. Tyndel, and W. Venus, Phys. Lett. B 112, 97 (1982).
  * Lyons  _et al._ [1981] L. Lyons, C. Albajar, and G. Myatt, Z. Phys. C 10, 215 (1981).
  * Workman and Others [2022] R. L. Workman and Others (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Weinberg [1979] S. Weinberg, Phys. Rev. Lett. 43, 1566 (1979).
  * Li  _et al._ [2020] T. Li, X.-D. Ma, and M. A. Schmidt, JHEP 10, 115 (2020), arXiv:2007.15408 [hep-ph] .
  * Babu and Leung [2001] K. S. Babu and C. N. Leung, Nucl. Phys. B 619, 667 (2001), arXiv:hep-ph/0106054 .
  * Deppisch  _et al._ [2020] F. F. Deppisch, K. Fridell, and J. Harz, JHEP 12, 186 (2020), arXiv:2009.04494 [hep-ph] .
  * Bhattacharya and Wudka [2016] S. Bhattacharya and J. Wudka, Phys. Rev. D 94, 055022 (2016), [Erratum: Phys.Rev.D 95, 039904 (2017)], arXiv:1505.05264 [hep-ph] .
  * Bischer and Rodejohann [2019] I. Bischer and W. Rodejohann, Nucl. Phys. B 947, 114746 (2019), arXiv:1905.08699 [hep-ph] .
  * Bhattacharya  _et al._ [2012] T. Bhattacharya, V. Cirigliano, S. D. Cohen, A. Filipuzzi, M. Gonzalez-Alonso, M. L. Graesser, R. Gupta, and H.-W. Lin, Phys. Rev. D 85, 054512 (2012), arXiv:1110.6448 [hep-ph] .
  * Cirigliano  _et al._ [2013] V. Cirigliano, S. Gardner, and B. Holstein, Prog. Part. Nucl. Phys. 71, 93 (2013), arXiv:1303.6953 [hep-ph] .
  * Hall and Suzuki [1984] L. J. Hall and M. Suzuki, Nucl. Phys. B 231, 419 (1984).
  * Davidson  _et al._ [1994] S. Davidson, D. C. Bailey, and B. A. Campbell, Z. Phys. C 61, 613 (1994), arXiv:hep-ph/9309310 .
  * Babu and Mohapatra [1995] K. S. Babu and R. N. Mohapatra, Phys. Rev. Lett. 74, 2418 (1995), arXiv:hep-ph/9410326 .
  * Babu and Barr [1997] K. S. Babu and S. M. Barr, Phys. Rev. D 56, 2614 (1997), arXiv:hep-ph/9512389 .
  * Babu and Pakvasa [2002] K. S. Babu and S. Pakvasa, (2002), arXiv:hep-ph/0204236 .
  * Babu  _et al._ [2016] K. S. Babu, D. W. McKay, I. Mocioiu, and S. Pakvasa, Phys. Rev. D 93, 113019 (2016), arXiv:1605.03625 [hep-ph] .
  * Jana [2022] S. Jana, Moscow Univ. Phys. Bull. 77, 371 (2022).
  * Cirigliano and Rosell [2007] V. Cirigliano and I. Rosell, Phys. Rev. Lett. 99, 231801 (2007), arXiv:0707.3439 [hep-ph] .
  * Note [1] Note that the baseline of SBND (110 m) is much shorter than the shortest neutrino oscillation length $L_{\rm osc}^{\rm atm}=(4\pi E/\Delta m^{2}_{31})\sim 2\times 10^{3}\tmspace+{.2777em}\rm km$, so there is no contribution from neutrino oscillation to the flux.
  * Mic [2018] (2018), 10.2172/1573216.
  * Tutto [2021] M. D. Tutto, in  _APS April Meeting 2021, Volume 66, Number 5_ (American Physical Society, Online, 2021) presentation at APS April Meeting 2021.
  * Bhadra  _et al._ [2014] S. Bhadra _et al._ , “Letter of intent to construct a nuprism detector in the j-parc neutrino beamline,” (2014), arXiv:1412.3086 [physics.ins-det] .
  * Abud  _et al._ [2021] A. A. Abud _et al._ , “Deep underground neutrino experiment (dune) near detector conceptual design report,” (2021), arXiv:2103.13910 [physics.ins-det] .
  * Machado  _et al._ [2019] P. A. Machado, O. Palamara, and D. W. Schmitz, Ann. Rev. Nucl. Part. Sci. 69, 363 (2019), arXiv:1903.04608 [hep-ex] .
  * Aguilar-Arevalo  _et al._ [2009] A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. D 79, 072002 (2009), arXiv:0806.1449 [hep-ex] .
  * Bierlich  _et al._ [2022] C. Bierlich, S. Chakraborty, N. Desai, L. Gellersen, I. Helenius, P. Ilten, L. Lönnblad, S. Mrenna, S. Prestel, C. T. Preuss, T. Sjöstrand, P. Skands, M. Utheim, and R. Verheyen, “A comprehensive guide to the physics and usage of pythia 8.3,” (2022), arXiv:2203.11601 [hep-ph] .
  * Palamara [2023] O. Palamara (INT, Seattle, 2023) presentation at the Theoretical Physics Uncertainties to Empower Neutrino Experiments.
  * Palamara [2016] O. Palamara (ArgoNeuT), JPS Conf. Proc. 12, 010017 (2016).
  * Acciarri  _et al._ [2021] R. Acciarri _et al._ (online, 2021) Snowmass 2021 Letter of Interest: Instrumentation Frontier.


solutionSolutionsolutionfile solutionfile[SIM_MISO_SIM_AssistedSystems]

# Performance of Double-Stacked Intelligent Metasurface-Assisted Multiuser
Massive MIMO Communications in the Wave Domain

Anastasios Papazafeiropoulos, Pandelis Kourtessis, Symeon Chatzinotas A.
Papazafeiropoulos is with the Communications and Intelligent Systems Research
Group, University of Hertfordshire, Hatfield AL10 9AB, U. K., and with SnT at
the University of Luxembourg, Luxembourg. P. Kourtessis is with the
Communications and Intelligent Systems Research Group, University of
Hertfordshire, Hatfield AL10 9AB, U. K. S. Chatzinotas is with the SnT at the
University of Luxembourg, Luxembourg. A. Papazafeiropoulos was supported by
the University of Hertfordshire’s 5-year Vice Chancellor’s Research
Fellowship. S. Chatzinotas was supported by the National Research Fund,
Luxembourg, under the project RISOTTI. E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Although reconfigurable intelligent surface (RIS) is a promising technology
for shaping the propagation environment, it consists of a single-layer
structure within inherent limitations regarding the number of beam steering
patterns. Based on the recently revolutionary technology, denoted as stacked
intelligent metasurface (SIM), we propose its implementation not only on the
base station (BS) side in a massive multiple-input multiple-output (mMIMO)
setup but also in the intermediate space between the base station and the
users to adjust the environment further as needed. For the sake of
convenience, we call the former BS SIM (BSIM), and the latter channel SIM
(CSIM). Hence, we achieve wave-based combining at the BS and wave-based
configuration at the intermediate space. Specifically, we propose a channel
estimation method with reduced overhead, being crucial for SIM-assisted
communications. Next, we derive the uplink sum spectral efficiency (SE) in
closed form in terms of statistical channel state information (CSI). Notably,
we optimize the phase shifts of both BSIM and CSIM simultaneously by using the
projected gradient ascent method (PGAM). Compared to previous works on SIMs,
we study the uplink transmission, a mMIMO setup, channel estimation in a
single phase, a second SIM at the intermediate space, and simultaneous
optimization of the two SIMs. Simulation results show the impact of various
parameters on the sum SE, and demonstrate the superiority of our optimization
approach compared to the alternating optimization (AO) method.

###### Index Terms:

Reconfigurable intelligent surface (RIS), stacked intelligent metasurfaces
(SIM), gradient projection, 6G networks.

##  I Introduction

In the pursuit of meeting the stringent demands of future wireless networks
systems [1], massive multiple-input multiple-output (mMIMO) systems and
millimeter-wave (mmWave) communications, suggested the previous years, have
failed to meet the low energy consumption and the ubiquitous wireless
connectivity requirements [2]. For example, mMIMO systems exhibit low
performance in poor scattering conditions despite their implementation with a
large number of active elements which may also consume excessive energy [3,
2]. In the case of mmWave communications, these require costly and power-
hungry transceivers that can be easily blocked by common obstacles such as
walls [4]. It is thus evident that it is crucial for next-generation networks
to introduce a radical paradigm, which will give priority to energy
sustainability and pervasive connectivity through the emergence of
programmable wireless environments.

In this direction, a prominent technology that can dynamically control the
wireless environment in an almost passive manner while enhancing significantly
the wireless channel quality is the reconfigurable intelligent surface (RIS)
[5, 6, 7, 8, 9]. Generally, a RIS consists of an artificial surface
implemented by a large number of cost-efficient nearly passive element, which
can induce independently phase shifts on the impinging electromagnetic (EM)
waves through a smart controller. Thus, the RIS is capable of shaping
dynamically the propagation environment and reducing the energy consumption
[10].

Despite that RIS has been suggested for various communication scenarios
because of its numerous advantages [5, 6, 11, 12, 13, 14, 8, 15, 16, 17],
existing research works have mostly considered single-layer metasurface
structures, which limit the degrees of freedom for managing the beams.
Moreover, it has been shown that RIS cannot suppress the inter-user
interference because of their single-layer configuration and hardware
limitations [18]. The identification of these gaps made the authors in [18]
propose a stacked intelligent metasurface (SIM) with remarkable signal
processing capabilities compared to conventional RIS with a single layer. The
proposed SIM is not a pure mathematical abstraction but it is based on
tangible hardware prototypes [19, 20] and the technological advancements in
wave-based computing. In particular, [19] proposed a deep neural network
($\mathrm{D}^{2}\mathrm{NN}$), which has a similar structure to a SIM to
perform parallel calculations at the speed of light by employing three-
dimensional (3D) printed optical lenses and taking advantage of the wave
properties of photons. However, this fabrication is fixed and cannot be
retrained to perform other tasks. To this end, in [20], a programmable
($\mathrm{D}^{2}\mathrm{NN}$), which is closer to a SIM was designed with each
reprogrammable artificial neuron acting as the SIM meta-atoms. Hence, they
demonstrated the execution of various tasks, e.g., image classification
through the flexible manipulation of the propagation of the EM waves through
its numerous layers.

These experimental results led to propose a SIM-based transceiver for point-
to-point MIMO communications [18], where two SIMs are implemented at the
transmitter and the receiver, respectively to perform precoding and combining
without any digital hardware, which reduced the need for a large number of
antennas while the EM signals propagate through them. Therein, the phase
shifts of each SIM were optimised in an alternating optimisation (AO) manner
based on instantaneous channel state information (CSI) conditions. In [16], we
proposed a simultaneous optimization of the phase shifts of both the
transmitter and the receiver. In [21], the authors deployed a SIM-enabled
(base station) to perform dowlink beamforming in the EM wave domain to
multiple users. However, these works relied on perfect CSI to design the SIMs.
In [22], a digital wave-based channel estimator was proposed, where multiple
subphases were first precoded in the wave domain, and then, in the digital
domain.

###  I-A Contributions

The reflections above led us to the crux of this work, which is to introduce
two SIMs in a mMIMO system with multiple users, one at the BS and one at the
intermediate space under the assumptions of correlated fading and imperfect
CSI, and study the uplink sum spectral efficiency (SE).

Our main contributions are summarised as follows:

  * •

Under the general setup of a mMIMO BS, we propose a new architecture including
two SIMs, and we obtain the uplink achievable SE with imperfect CSI and
maximum-ratio combining (MRC) in closed form that hinges only on large-scale
statistics based on the two time-scale protocol. Specifically, one SIM is
integrated into the BS to enable decoding in the electromagnetic wave domain,
and the other SIM is aimed to assist the communication between the users and
the BS in the same domain.

  * •

We perform channel estimation based on the linear minimum mean-squared error
(LMMSE) method. We also remark that channel estimation of BSIM and CSIM-
assisted mMIMO systems has not taken place before.

  * •

Different from [21], not only do we consider a SIM-based BS, but we also
consider a SIM in the intermediate space to aid the communication between the
users and the BS. Also, we rely on statistical CSI instead of instantaneous-
based analysis in [21]. Moreover, we consider imperfect CSI, which is of
practical interest. Compared to [22], which considered several subfaces to
estimate the channel in SIM-assisted multi-user systems, we manage to update
the estimated channel in a single phase. Thus, the approach in [22], requiring
too much overhead, reduces the SE, while our approach does not appear this
significant disadvantage. Note that [22] assumes a simple SIM at the BS, while
we consider an additional SIM at the intermediate space.

  * •

The proposed optimization is based on statistical CSI, which means it can be
executed at every several coherence intervals. As a result, we achieve reduced
overhead and computational complexity. Previous works on SIM-based
communications relied on instantaneous CSI, which changes at each coherence
interval, and is not suitable for practical SIM-based systems that include
many surfaces, and thus bring high burden.

  * •

We conceive the problem optimizing the phase shifts of both SIMs. Despite the
nonconvexity of this problem, we propose a projected gradient ancient method
(PGAM), which optimises the phase shifts of both SIMs simultaneously. This is
a noteworthy contribution since [18] optimized the phase shifts in an AO way.
We perform a comparison between the two approaches, and show the superiority
of our approach.

  * •

Simulations show the superiority of the proposed architecture, examine which
SIM has the greatest impact, and shed light on the effect of various system
parameters on the sum SE.

Paper Outline: The remainder of this paper is organized as follows. Section II
describes the system and channel models of a BSIM and CSIM-assisted mMIMO
system. Section III elaborates on the CE and the uplink data transmission with
the obtained uplink sum SE. In Section IV, we present the problem formulation
and the simultaneous optimization of both BSIM and CSIM. The numerical results
are discussed in Section LABEL:Numerical. Section VI concludes the paper.

Notation: Vectors and matrices are described by boldface lower and upper case
symbols, respectively. The notations $(\cdot)^{\scriptscriptstyle\mathsf{T}}$,
$(\cdot)^{\scriptscriptstyle\mathsf{H}}$, and $\tr\\!\left({\cdot}\right)$
denote the transpose, Hermitian transpose, and trace operators, respectively.
The notation $\mathbb{E}\left[\cdot\right]$ expresses the expectation
operator, and $\text{diag}\left(\right)\left({\mathbf{A}}\right)$ describes a
vector with elements equal to the diagonal elements of ${\mathbf{A}}$. The
notation $\text{diag}\left(\right)\left({\mathbf{x}}\right)$ denotes a
diagonal matrix whose elements are ${\mathbf{x}}$, while
${\mathbf{b}}\sim{\cal C}{\cal N}{({\mathbf{0}},\mathbf{\Sigma})}$ denotes a
circularly symmetric complex Gaussian vector with zero mean and a covariance
matrix $\mathbf{\Sigma}$. In the case of complex-valued $\mathbf{x}$ and
${\mathbf{y}}$, we denote
$\langle\mathbf{x},\mathbf{y}\rangle=2\Re{\mathbf{x}^{{\scriptscriptstyle\mathsf{H}}}\mathbf{y}}$.

##  II System and Channel Models

Figure 1:  A BSIM and CSIM-assisted mMIMO system with multiple users.

###  II-A System Model

We consider the uplink of a multiuser MISO wireless system, where the BS is
equipped with a large number of $M_{\mathrm{BS}}$ antennas that serve $K$
single-antennas users under the same time-frequency resources. In this
context, we employ the BSIM, integrated into the BS to enable receive decoding
in the EM wave domain, and the CSIM to enhance the communication between the
users and the BS in the wave domain, as shown in Fig. 1. For the fabrication
of BSIM, we utilize an array of $L$ metasurfaces, where each one of them
consists of $M$ meta-atoms. The corresponding sets are
$\mathcal{L}=\\{1,\ldots,L\\}$ and $\mathcal{M}=\\{1,\ldots,M\\}$. Similarly,
the CSIM includes an array of $S$ metasurfaces, each including $N$ meta-
atoms.111Without any loss of generality, we assume that the metasurfaces at
each SIM include an identical number of meta-atoms. The corresponding sets are
$\mathcal{S}=\\{1,\ldots,S\\}$ and $\mathcal{N}=\\{1,\ldots,N\\}$. Both SIMs
are connected to independent smart controllers, which can adjust the
corresponding phase shifts of the EM waves that are transmitted through each
of their meta-atoms. Note that the forward propagation process in each SIM
resembles a fully connected ANN, where training the interconnection
architecture of the BSIM enables uplink decoding, while the respective
training of the CSIM enables shaping the propagation environment. We denote
${\bm{\Phi}}^{l}=\text{diag}\left(\right)({\bm{\phi}}^{l})\in\mathbb{C}^{M\times
M}$, where
${\bm{\phi}}^{l}=[\phi^{l}_{1},\dots,\phi^{l}_{M}]^{{\scriptscriptstyle\mathsf{T}}}\in\mathbb{C}^{N\times
1}$ is the transmission coefficient matrix of the BSIM. Note that
$\phi_{m}^{l}=e^{j\theta_{m}^{l}}$, where
$\theta_{m}^{l}\in[0,2\pi),m\in\mathcal{M},l\in\mathcal{L}$ is the phase shift
by the $m$-th meta-atom on the $k$-th transmit metasurface layer. Similarly,
for the CSIM, we denote
${\bm{\Lambda}}^{s}=\text{diag}\left(\right)({\bm{\lambda}}^{s})\in\mathbb{C}^{M\times
N}$, where
${\bm{\lambda}}^{s}=[\lambda^{s}_{1},\dots,\lambda^{s}_{M}]^{{\scriptscriptstyle\mathsf{T}}}\in\mathbb{C}^{N\times
1}$. Here, $\lambda_{n}^{s}=e^{j\xi_{n}^{s}}$ with
$\xi_{n}^{s}\in[0,2\pi),n\in\mathcal{N},s\in\mathcal{S}$ being the phase shift
by the $n$-th meta-atom on the $s$-th CSIM metasurface layer.222We assume that
the phase shifts are continuously-adjustable. Also, we assume that the modulus
for the coefficient matrices of both SIMs equals to $1$ to evaluate the
maximization of the achievable rate, e.g., see [5]. The study of more
practical designs such as the consideration of coupled phase and magnitude
[23] could be the topic of future work.

###  II-B Channel Model

The transmission coefficient from the $\tilde{m}$th meta-atom on the $(l-1)$st
BSIM metasurface layer to the $m$th meta-atom on the $l$th BSIM metasurface
layer, provided by the Rayleigh-Sommerfeld diffraction theory [19], is written
as

$\displaystyle w_{m,\tilde{m}}^{l}=\frac{A_{t}\cos
x_{m,\tilde{m}}^{l}}{r_{m,\tilde{m}}^{l}}\left(\frac{1}{2\pi
r^{l}_{m,\tilde{m}}}-j\frac{1}{\lambda}\right)e^{j2\pi
r_{m,\tilde{m}}^{l}/\lambda},l\in\mathcal{L},$ (1)

where $A_{t}$ expresses the area of each meta-atom, $x_{m,\tilde{m}}^{l}$ is
the angle between the propagation direction and the normal direction of the
$(l-1)$th BSIM metasurface layer, and $r_{m,\tilde{m}}^{l}$, is the
corresponding transmission distance. As a result, the effect of the BSIM can
be modeled as

$\displaystyle{\mathbf{P}}={\bm{\Phi}}^{L}{\mathbf{W}}^{L}\cdots{\bm{\Phi}}^{2}{\mathbf{W}}^{2}{\bm{\Phi}}^{1}\in\mathbb{C}^{M\times
M},$ (2)

where ${\mathbf{W}}^{l}\in\mathbb{C}^{M\times M},l\in\mathcal{L}/\\{1\\}$ is
the transmission coefficient matrix between the $(l-1)st$ BSIM metasurface
layer and the $l$th transmit metasurface layer. Regarding
${\mathbf{W}}^{1}\in\mathbb{C}^{M\times M_{\mathrm{BS}}}$, it is the
transmission coefficient matrix between the $M_{\mathrm{BS}}$ transmit
antennas and the first metasurface layer of the BSIM.

In the CSIM, the transmission coefficient from the $n$th meta-atom on the
$s$th metasurface layer to the $\tilde{n}$th meta-atom on the $(s-1)$st
metasurface layer can be written as

$\displaystyle
u_{\tilde{n},n}^{s}=\frac{A_{r}\cos\zeta_{\tilde{n},n}^{s}}{t_{\tilde{n},n}^{s}}\left(\frac{1}{2\pi
t^{k}_{\tilde{n},n}}-j\frac{1}{\lambda}\right)e^{j2\pi
t_{\tilde{n},n}^{s}/\lambda},s\in\mathcal{S},$ (3)

where $A_{r}$ denotes the area of each meta-atom in the CSIM,
$\zeta_{\tilde{n},n}^{s}$ expresses the angle between the propagation
direction and the normal direction of the $(s-1)$th CSIM metasurface layer,
and $t_{\tilde{n},n}^{s}$ is the corresponding transmission distance. In this
case, the effect of the CSIM is described by

$\displaystyle{\mathbf{Z}}={\bm{\Lambda}}^{S}{\mathbf{U}}^{S}\cdots{\bm{\Lambda}}^{2}{\mathbf{U}}^{2}{\bm{\Lambda}}^{1}\in\mathbb{C}^{N\times
N},$ (4)

where ${\mathbf{U}}^{s}\in\mathbb{C}^{N\times N},s\in\mathcal{S}/\\{1\\}$ is
the transmission coefficient matrix between the $s$th CSIM metasurface layer
to the $(s-1)$st CSIM metasurface layer, and
${\mathbf{U}}^{1}\in\mathbb{C}^{N\times 1}$ is the transmission coefficient
matrix from the output metasurface layer of the CSIM SIM to user $k$.

Based on narrowband quasi-static block fading channel modeling, we assume that
the duration of each block is $\tau_{\mathrm{c}}$ channel uses. Also, by
adopting the time-division-duplex (TDD) protocol as the recommended protocol
for next-generation systems such as mMIMO systems, $\tau$ channel uses are
allocated for the uplink training phase and $\tau_{\mathrm{c}}-\tau$ channel
uses for the uplink data transmission phase.

Now, we denote
${\mathbf{G}}=[{\mathbf{g}}_{1}\ldots,{\mathbf{g}}_{N}]\in\mathbb{C}^{M\times
N}$ the channel matrix between the BSIM and the CSIM, where
${\mathbf{g}}_{i}\in\mathbb{C}^{M\times 1}$ for $i\in\mathcal{N}$. Moreover,
${\mathbf{q}}_{k}\in\mathbb{C}^{N\times 1}$ expresses the channel between the
CSIM and user $k$. Also, we assume the presence of a direct link between the
BSIM and user $k$. We denote the corresponding link as ${\mathbf{d}}_{k}$. In
this work, we assume the presence of realistic correlated Rayleigh fading,
which occurs in practice [7]. The analysis with correlated Rician fading that
includes a line-of-sight (LoS) component could be the topic of future work.
Specifically, we have

$\displaystyle{\mathbf{G}}$
$\displaystyle=\sqrt{\tilde{\beta}_{g}}{\mathbf{R}}_{\mathrm{BSIM}}^{1/2}{\mathbf{D}}{\mathbf{R}}_{\mathrm{CSIM}}^{1/2},$
(5) $\displaystyle{\mathbf{q}}_{k}$
$\displaystyle=\sqrt{\tilde{\beta}_{k}}{\mathbf{R}}_{\mathrm{CSIM}}^{1/2}{\mathbf{c}}_{k},$
(6) $\displaystyle{\mathbf{d}}_{k}$
$\displaystyle=\sqrt{\bar{\beta}_{k}}{\mathbf{R}}_{\mathrm{BSIM}}^{1/2}\bar{{\mathbf{c}}}_{k},$
(7)

where ${\mathbf{R}}_{\mathrm{BSIM}}\in\mathbb{C}^{M\times M}$ and
${\mathbf{R}}_{\mathrm{CSIM}}\in\mathbb{C}^{N\times N}$ describe the
correlation matrices at the BSIM and the CSIM respectively. Note that these
matrices are deterministic Hermitian-symmetric positive semi-definite and are
assumed to be known by the network. The spatial correlation matrices at the
BSIM and CSIM are given by [7]

$\displaystyle[{\mathbf{R}}_{\mathrm{BSIM}}]_{m,\tilde{m}}$
$\displaystyle=\mathrm{sinc}(2r_{m,\tilde{m}}/\lambda),\tilde{m}\in\mathcal{M},m\in\mathcal{M},$
(8) $\displaystyle[{\mathbf{R}}_{\mathrm{CSIM}}]_{\tilde{n},n}$
$\displaystyle=\mathrm{sinc}(2t_{\tilde{n},n}/\lambda),n\in\mathcal{N},\tilde{n}\in\mathcal{N},$
(9)

where $r_{m,\tilde{m}}$ and $t_{\tilde{n},n}$ are the corresponding meta-atom
spacings. The path-losses $\tilde{\beta}_{g}$, $\bar{\beta}_{k}$, and
$\tilde{\beta}_{k}$ correspond to BSIM-CSIM, CSIM-user $k$, and BSIM-user $k$
links, while
${\mathbf{c}}_{k}\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{N}\right)$,
$\bar{{\mathbf{c}}}_{k}\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{M}\right)$,
and
$\mathrm{vec}({\mathbf{D}})\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{MN}\right)$
represent the small-scale fading parts. Contrary to most previous works, which
assume that one of the links of the cascaded channels is deterministic, i.e.,
it represents a LoS component [24, 8], we consider a more general analysis,
since we consider small-scale fading among all links.

By assuming that the phase shift matrices ${\bm{\Phi}}^{l},l\in\mathcal{L}$
and ${\bm{\Lambda}}^{s},s\in\mathcal{S}$ are given, the aggregated channel
vector from the last BSIM layer to user $k$, includes the cascaded and direct
channels, i.e.,
${\mathbf{h}}_{k}={\mathbf{G}}{\mathbf{Z}}{\mathbf{U}}^{1}{\mathbf{q}}_{k}+{\mathbf{d}}_{k}\in\mathbb{C}^{M\times
1}$, which has a covariance matrix given by
${\mathbf{R}}_{k}=\mathbb{E}\\{{\mathbf{h}}_{k}{\mathbf{h}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}$.
Specifically, we have

$\displaystyle{\mathbf{R}}_{k}$
$\displaystyle=\hat{\beta}_{k}\mathbb{E}\\{{\mathbf{G}}{\mathbf{Z}}{\mathbf{U}}^{1}\mathbb{E}\\{{\mathbf{q}}_{k}{\mathbf{q}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}{\mathbf{U}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{Z}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}\\}+\bar{\beta}_{k}{\mathbf{R}}_{\mathrm{BSIM}}$
(10)
$\displaystyle=\hat{\beta}_{k}\mathbb{E}\\{{\mathbf{G}}{\mathbf{Z}}{\mathbf{U}}^{1}{\mathbf{R}}_{\mathrm{CSIM}}{\mathbf{U}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{Z}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{G}}^{{\scriptscriptstyle\mathsf{H}}}\\}+\bar{\beta}_{k}{\mathbf{R}}_{\mathrm{BSIM}}$
(11)
$\displaystyle=\hat{\beta}_{k}\tr({\mathbf{R}}_{\mathrm{CSIM}}{\mathbf{Z}}\bar{{\mathbf{R}}}_{\mathrm{CSIM}}{\mathbf{Z}}^{{\scriptscriptstyle\mathsf{H}}}){\mathbf{R}}_{\mathrm{BSIM}}+\bar{\beta}_{k}{\mathbf{R}}_{\mathrm{BSIM}},$
(12)

where, in (10), we have applied the independence between ${\mathbf{G}}$ and
${\mathbf{q}}_{k}$, denoted
$\hat{\beta}_{k}=\tilde{\beta}_{g}\tilde{\beta}_{k}$, and used
$\mathbb{E}\\{{\mathbf{d}}_{k}{\mathbf{d}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}=\bar{\beta}_{k}{\mathbf{R}}_{\mathrm{BSIM}}$.
Next, we have used that
$\mathbb{E}\\{{\mathbf{q}}_{k}{\mathbf{q}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}=\tilde{\beta}_{k}{\mathbf{R}}_{\mathrm{CSIM}}$.
In (12), we have denoted
$\bar{{\mathbf{R}}}_{\mathrm{CSIM}}={\mathbf{U}}^{1}{\mathbf{R}}_{\mathrm{CSIM}}{\mathbf{U}}^{1^{{\scriptscriptstyle\mathsf{H}}}}$,
and we have used the property
$\mathbb{E}\\{{\mathbf{V}}{\mathbf{U}}{\mathbf{V}}^{{\scriptscriptstyle\mathsf{H}}}\\}=\tr({\mathbf{U}}){\bm{\mathrm{I}}}_{M}$,
where ${\mathbf{U}}$ is a deterministic square matrix, and ${\mathbf{V}}$ is
any matrix with independent and identically distributed (i.i.d.) entries of
zero mean and unit variance.

##  III Channel Estimation and Uplink Data Transmission

###  III-A Channel Estimation

Given that perfect CSI is not attainable in practice, we employ the TDD
protocol that considers an uplink training phase. In this phase, pilot symbols
are transmitted to estimate the channel. However, a RIS has the peculiar
characteristic that it cannot process any signals, i.e., it cannot obtain
receive pilots and it cannot transmit them because it is implemented without
any RF chains.

Herein, we propose a channel estimation method, where we estimate the
aggregated channel as in [24, 8, 25] instead of the individual channels as in
[13, 26, 27]. The latter approach comes with extra power and hardware cost.
Moreover, the large dimension BSIM-CSIM channel would require a pilot overhead
that would be prohibitively high if we estimated the individual channels.
Also, we manage to express the estimated channel in closed form. Hence, the
proposed method follows.

During the training phase, which lasts $\tau$ channel uses, it is assumed that
all users send orthogonal pilot sequences. In particular,
${\mathbf{x}}_{k}=[x_{k,1},\ldots,x_{k,\tau}]^{{\scriptscriptstyle\mathsf{H}}}\in\mathbb{C}^{\tau\times
1}$ with
${\mathbf{x}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{x}}_{l}=0~{}\forall
k\neq l$ and
${\mathbf{x}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{x}}_{k}=\tau\rho$
joules denotes the pilot sequence of user $k$, where all users use the same
normalized signal-to-noise ratio (SNR) $\rho$ for transmitting each pilot
symbol during the training phase.

Overall, the BS receives

$\displaystyle{\mathbf{Y}}^{\tr}=\sum_{i=1}^{K}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{h}}_{i}{\mathbf{x}}_{i}^{{\scriptscriptstyle\mathsf{H}}}+{\mathbf{Z}}^{\tr},$
(13)

with ${\mathbf{Z}}^{\tr}\in\mathbb{C}^{M\times\tau}$ being the received AWGN
matrix having independent columns with each one distributed as
$\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{M}\right)$.

By multiplying (13) with ${\mathbf{x}}_{k}$, we obtain

$\displaystyle{\mathbf{r}}_{k}={\mathbf{c}}_{k}+\frac{{\mathbf{z}}_{k}}{\tau\rho},$
(14)

where
${\mathbf{c}}_{k}={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{h}}_{k}$
and ${\mathbf{z}}_{k}={\mathbf{Z}}^{\tr}{\mathbf{x}}_{k}$.

Based on linear MMSE estimation, the overall perfect channel can be written as

$\displaystyle{\mathbf{c}}_{k}=\hat{{\mathbf{c}}}_{k}+\tilde{{\mathbf{c}}}_{k},$
(15)

where $\hat{{\mathbf{c}}}_{k}$ and $\tilde{{\mathbf{c}}}_{k}$ are channel
estimate and the estimation channel error vectors, respectively. The following
lemma provides the linear MMSE estimated channel. Note that
$\hat{{\mathbf{c}}}_{k}$ and $\tilde{{\mathbf{c}}}_{k}$ are uncorrelated and
each of them has zero mean, but they are not independent [28].

###### Lemma 1

The linear MMSE estimated channel ${\mathbf{c}}_{k}$ of user $k$ at the BS is
obtained as

$\displaystyle\hat{{\mathbf{c}}}_{k}=\hat{{\mathbf{R}}}_{k}{\mathbf{Q}}_{k}{\mathbf{r}}_{k},$
(16)

where
$\hat{{\mathbf{R}}}_{k}={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{R}}_{k}{\mathbf{W}}^{1}{\mathbf{P}}$,
${\mathbf{Q}}_{k}\\!=\\!\left(\\!\hat{{\mathbf{R}}}_{k}\\!+\\!\frac{\sigma^{2}}{\tau\rho}{\bm{\mathrm{I}}}_{M}\\!\right)^{\\!-1}$,
and ${\mathbf{r}}_{k}$ is the noisy channel given by (14).

###### Proof:

Please see Appendix A. ∎

The covariances of $\hat{{\mathbf{g}}}_{k}$ and $\tilde{{\mathbf{g}}}_{k}$ are
easily obtained as

$\displaystyle{\bm{\Psi}}_{k}\\!$
$\displaystyle=\\!\hat{{\mathbf{R}}}_{k}{\mathbf{Q}}_{k}\hat{{\mathbf{R}}}_{k},$
(17) $\displaystyle\tilde{{\bm{\Psi}}}_{k}$
$\displaystyle=\hat{{\mathbf{R}}}_{k}-{\bm{\Psi}}_{k}.$ (18)

As can be seen, in this method, the length of the pilot sequence depends only
on the number of users $K$, and is independent of the number of surfaces and
elements at the BSIM and CSIM, which are generally large. Thus, we achieve
substantial overhead reduction compared to estimating the individual channel,
where the complexity would increase with increasing SIMs sizes. Also, we
achieved to obtain the estimated overall channel in closed form instead of
works such as [22] that cannot result in analytical expressions.

Another benefit of the proposed two-timescale method is that the channel is
required to be estimated when the large-scale statistics (statistical CSI)
varies, which is at every several coherence intervals. On the contrary,
instantaneous CSI-based approaches would require channel estimation at every
coherence interval that would result in high power consumption and
computational complexity.

###  III-B Uplink Data Transmission

In this subsection, we provide the derivation of the uplink achievable sum SE
of a SIM-assisted mMIMO system.

During the uplink data transmission phase, the BS receives

$\displaystyle{\mathbf{y}}=\sqrt{\rho}\sum_{i=1}^{K}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{h}}_{i}s_{i}+{\mathbf{n}},$
(19)

where
${\mathbf{n}}\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{M}\right)$
is the AWGN at the BS, $s_{i}\sim\mathcal{CN}(0,1)$ is the transmit signal
vector, $\rho$ is the uplink SNR, and ${\mathbf{h}}_{i}$ is the aggregated
channel from the last BSIM layer to user $i$.

After application of a receive combining vector ${\mathbf{v}}_{k}$ as
$\hat{s}_{k}={\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{y}}$,
the BS obtains $s_{k}$ from user $k$ as

$\displaystyle\hat{s}_{k}$
$\displaystyle=\sqrt{\rho}\mathbb{E}\left\\{{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{{\mathbf{h}}}_{k}\right\\}s_{k}$
$\displaystyle+\sqrt{\rho}\left({\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{{\mathbf{h}}}_{k}-\mathbb{E}\left\\{{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{{\mathbf{h}}}_{k}\right\\}\right)s_{k}$
$\displaystyle+\sum_{i\neq
k}^{K}\sqrt{\rho}{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{{\mathbf{h}}}_{i}s_{i}+{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{n}},$
(20)

where we have also accounted for channel hardening, which is applicable in
mMIMO systems with $M>8$ [28]. In (20), the first term is the desired signal,
the second term is the beamforming gain uncertainty, the third term is
multiuser interference, and the last term is the post-processed AWGN noise.

Next, we derive the use-and-then-forget bound on the uplink average SE, which
is a lower bound [28]. For this reason, we assume the worst-case uncorrelated
additive noise scenario for the inter-user interference and the AWGN noise. In
this case, the achievable SE can be written as

$\displaystyle\mathrm{SE}_{k}=\frac{\tau_{\mathrm{c}}-\tau}{\tau_{\mathrm{c}}}\log_{2}\left(1+\gamma_{k}\right)\\!,$
(21)

where the fraction represents the percentage of samples transmitted during the
uplink data transmission, and $\gamma_{k}=\frac{S_{k}}{I_{k}}$ denotes the
uplink signal-to-interference-plus-noise ratio (SINR). Here, $S_{k}$ and
$I_{k}$ describe the desired signal power and the interference plus noise
power. In terms of ${\mathbf{g}}_{k}$, we have

$\displaystyle S_{k}$
$\displaystyle=\rho|\mathbb{E}\left\\{{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}\right\\}|^{2},$
(22) $\displaystyle I_{k}$
$\displaystyle=\mathbb{E}\left\\{|{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}-\mathbb{E}\left\\{{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}\right\\}|^{2}\right\\}+\sum_{i\neq
k}^{K}\mathbb{E}\left\\{|{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{i}|^{2}\right\\}$
$\displaystyle+\\!\mathbb{E}\left\\{|{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{n}}|^{2}\right\\}\\!.$
(23)

In this work, we assume MRC decoding to provide a closed-form SE, i.e.,
${\mathbf{v}}_{k}=\hat{{\mathbf{c}}}_{k}$.333The application of the MMSE
decoder, which is another common linear receiver used for mMIMO systems due to
its better performance, is the topic of ongoing work. The reason is that the
corresponding analysis requires the use of deterministic equivalent tools [29,
30, 31], which result in complicated expressions that would result in an
intractable optimization regarding the phase shifts matrices including in BSIM
and CSIM.

###### Theorem 1

The uplink SINR of user $k$ for given SIMs in a BSIM and CSIM assisted mMIMO
system is tightly approximated as

$\gamma_{k}\ {\color[rgb]{0,0,0}\approx}\
\frac{S_{k}}{\color[rgb]{0,0,0}\tilde{I}_{k}},$ (24)

where

$\displaystyle{S}_{k}$ $\displaystyle=\tr^{2}\left({\bm{\Psi}}_{k}\right)\\!,$
(25) $\displaystyle\tilde{I}_{k}$
$\displaystyle=\sum_{i=1}^{K}\tr\\!\left(\hat{{\mathbf{R}}}_{i}{\bm{\Psi}}_{k}\right)-\tr\left({\bm{\Psi}}_{k}^{2}\right)+\frac{1}{\rho}\tr({\bm{\Psi}}_{k}).$
(26)

###### Proof:

Please see Appendix B. ∎

Theorem 1 provides a tight approximation as the numerical results in Section
LABEL:Numerical reveal. The tightness is also expected because it comes as a
result of the exploitation of channel hardening appearing in mMIMO systems
[28].

##  IV Problem Formulation and Optimization

It is crucial to optimize the sum SE in terms of the phase shift matrices of
the SIMs included in the proposed architecture.

###  IV-A Problem Formulation

By assuming infinite-resolution phase shifters, the formulation of the
optimization problem follows as

(27a) (27b) (27c) (27d) (27e)

In problem $(\mathcal{P})$, we have denote
$f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$ the objective function that includes
the approximate of the SINR provided by Theorem (1).

###  IV-B Simultaneous Optimization of Both SIMs

The non-convexity of problem $(\mathcal{P})$ and the coupling between the
optimization variables of the two SIMs, which also have constant modulus
constraints lead us to seek a locally optimal solution. Note that many
previous optimization methods on RIS-aided systems have applied alternating
optimization (AO) to optimize the relevant variables in each case [5, 32].
Although AO-based methods can be implemented easily, their convergence may
require many iterations, which increase with increasing the size of RIS [33].
Since SIM-assisted systems are expected to include a large number of elements,
the AO methods is not suggested.

Given that the projection operators can be obtained in closed form, we can
apply the projected gradient ascent method (PGAM) [34, Ch. 2] for the
simultaneous optimization of ${\bm{\phi}}_{l}$ and ${\bm{\lambda}}_{s}$
instead of the AO method. Numerical results in Section 27bk will demonstrate
the superiority of the proposed approach.

###  IV-C Proposed Algorithm-PGAM

According to the proposed PGAM, starting from
$({\bm{\phi}}_{l}^{0},{\bm{\lambda}}_{l}^{0})$, we shift towards the gradient
of the objective function, i.e., $\nabla
f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$, where the rate of change of
$f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$ becomes maximum. The newly computed
points ${\bm{\phi}}_{l}$ and ${\bm{\lambda}}_{s}$ are projected onto
$\Phi_{l}$ and $\Lambda_{s}$, respectively to keep the updated points inside
the feasible set.

The step size of this shift is determined by a parameter $\mu>0$. For the sake
of description, we define the following sets.

$\displaystyle\Phi_{l}$
$\displaystyle=\\{{\bm{\phi}}_{l}\in\mathbb{C}^{M\times
1}:|\phi^{l}_{i}|=1,i=1,\ldots,M\\},$ (27ab) $\displaystyle\Lambda_{s}$
$\displaystyle=\\{{\bm{\lambda}}_{s}\in\mathbb{C}^{N\times
1}:|\lambda^{s}_{i}|=1,i=1,\ldots,N\\}.$ (27ac)

The proposed PGAM solving (27) is outlined in Algorithm 1. As can be seen, the
algorithm includes the following two iterations.

$\displaystyle{\bm{\phi}}_{l}^{i+1}$
$\displaystyle=P_{\Phi_{l}}({\bm{\phi}}_{l}^{i}+\mu_{i}\nabla_{{\bm{\phi}}_{l}}f({\bm{\phi}}_{l}^{i},{\bm{\lambda}}_{s}^{i}))$
(27ad) $\displaystyle{\bm{\lambda}}_{s}^{i+1}$
$\displaystyle=P_{\Lambda_{s}}({\bm{\lambda}}_{s}^{i}+\mu_{i}\nabla_{{\bm{\lambda}}_{s}}f({\bm{\phi}}_{l}^{i},{\bm{\lambda}}_{s}^{i})),$
(27ae)

where $P_{\Phi_{l}}(\cdot)$ and $P_{\Lambda_{s}}(\cdot)$ are the projections
onto $\Phi_{l}$ and $\Lambda_{s}$, respectively. Also, to find the step size,
we apply the Armijo-Goldstein backtracking line search. For this reason, we
define the quadratic approximation of $f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$
as

$\displaystyle
Q_{\mu}({\bm{\phi}}_{l},{\bm{\lambda}}_{s};{\mathbf{x}},{\mathbf{y}})=f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})+\langle\nabla_{{\bm{\phi}}_{l}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s}),{\mathbf{x}}-{\bm{\phi}}_{l}\rangle$
$\displaystyle-\frac{1}{\mu}\|{\mathbf{x}}-{\bm{\phi}}_{l}\|^{2}_{2}+\langle\nabla_{{\bm{\lambda}}_{s}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s}),{\mathbf{y}}-{\bm{\lambda}}_{s}\rangle$
$\displaystyle-\frac{1}{\mu}\|{\mathbf{y}}-{\bm{\lambda}}_{s}\|^{2}_{2}.$
(27af)

The step size in (27ad) and (27ae) can be found as
$\mu_{i}=L_{i}\kappa^{m_{i}}$ with $L_{i}>0$, and $\kappa\in(0,1)$, where the
parameter $m_{n}$ is the smallest nonnegative integer satisfying

$\displaystyle f({\bm{\phi}}_{l}^{i+1},{\bm{\lambda}}_{s}^{i+1})\geq
Q_{L_{i}\kappa^{m_{i}}}({\bm{\phi}}_{l}^{i},{\bm{\lambda}}_{s}^{i};{\bm{\phi}}_{l}^{i+1},{\color[rgb]{0,0,0}{\bm{\lambda}}_{s}^{i+1}}),$
(27ag)

which can be performed by an iterative procedure.

Algorithm 1 Simultaneous PGAM for both SIMS

1: Input: ${\bm{\phi}}_{l}^{0},{\bm{\lambda}}_{s}^{0},\mu_{1}>0$,
$\kappa\in(0,1)$

2: $i\leftarrow 1$

3: repeat

4: repeat

5:
${\bm{\phi}}_{l}^{i+1}=P_{\Phi_{l}}({\bm{\phi}}_{l}^{i}+\mu_{ni}\nabla_{{\bm{\phi}}_{l}}f({\bm{\phi}}_{l}^{i},{\bm{\lambda}}_{s}^{i}))$

6:
${\bm{\lambda}}_{s}^{i+1}=P_{\Lambda_{s}}({\bm{\lambda}}_{s}^{i}+\mu_{i}\nabla_{{\bm{\lambda}}_{s}}f({\bm{\phi}}_{l}^{i},{\bm{\lambda}}_{s}^{i}))$

7: if $f({\bm{\phi}}_{l}^{i+1},{\bm{\lambda}}_{s}^{i+1})\leq
Q_{\mu_{i}}({\bm{\phi}}_{l}^{i},{\bm{\lambda}}_{s}^{i};{\bm{\phi}}_{l}^{i+1},{\color[rgb]{0,0,0}{\bm{\lambda}}_{s}^{i+1}})$
then

8: $\mu_{i}=\mu_{i}\kappa$

9: end if

10: until
$f({\bm{\phi}}_{l}^{i+1},{\bm{\lambda}}_{s}^{i+1})>Q_{\mu_{i}}({\bm{\phi}}_{l}^{i},{\bm{\lambda}}_{s}^{i};{\bm{\phi}}_{l}^{i+1},{\bm{\lambda}}_{s}^{i+1})$

11: $\mu_{i+1}\leftarrow\mu_{i}$

12: $i\leftarrow i+1$

13: until convergence

14: Output: ${\bm{\phi}}_{l}^{i+1},{\bm{\lambda}}_{s}^{i+1}$

###  IV-D Complex-valued Gradients

Below, we present the gradients
$\nabla_{{\bm{\phi}}_{l}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$ and
$\nabla_{{\bm{\lambda}}_{s}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$.

###### Proposition 1

The gradients of $f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$ regarding
${\bm{\phi}}_{l}^{*}$ and ${\bm{\lambda}}_{s}^{*}$ are provided closed-forms
by

$\displaystyle\nabla_{{\bm{\phi}}_{l}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$
$\displaystyle=\frac{\tau_{c}-\tau}{\tau_{c}\log_{2}(e)}\sum_{k=1}^{K}\frac{\tilde{I}_{k}\nabla_{{\bm{\phi}}_{l}}S_{k}-S_{k}\nabla_{{\bm{\phi}}_{l}}\tilde{I}_{k}}{(1+\gamma_{k})\tilde{I}_{k}},$
(27ah)
$\displaystyle\nabla_{{\bm{\lambda}}_{s}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$
$\displaystyle=\frac{\tau_{c}-\tau}{\tau_{c}\log_{2}(e)}\sum_{k=1}^{K}\frac{\tilde{I}_{k}\nabla_{{\bm{\lambda}}_{s}}S_{k}-S_{k}\nabla_{{\bm{\lambda}}_{s}}\tilde{I}_{k}}{(1+\gamma_{k})\tilde{I}_{k}},$
(27ai)

where

Figure 2:

###  V-C Sum SE Performance Evaluation

First, Figs. V-C.(a) and V-C.(b) illustrate the sum SE versus the numbers of
meta-atoms $M$ and $N$ on each metasurface layer of the BSIM and CSIM,
respectively. In both figures, we observe an increase in the sum SE with the
sizes of the corresponding surfaces. However, the impact of the CSIM is
greater since a double increase in $N$ results in $170\%$ improvement compared
to $200\%$ improvement when $M$ becomes almost double (from $40$ to $100$).
Also, in Fig.V-C.(b), we observe the impact of the CSIM. In particular, we
have shown the SE when the CSIM is absent, and only the direct signal exists,
the performance is quite low.

(a) Achievable sum SE of the proposed double-SIM mMIMO architecture versus the
number of meta-atoms of the BSIM $M$ while varying the number of meta-atoms of
the CSIM $N$. Figure 3:

Figs. V-C.(a) and V-C.(b) depict the sum SE versus the numbers of layers $L$
and $S$ of the BSIM and CSIM, respectively. In Fig. V-C.(a), we observe that
the sum SE improves until $L=4$ because of the ability of the BSIM to mitigate
the inter-user interference in the EM wave domain. In Fig. V-C.(b), it is
shown that the sum SE increases with the number of metasurfaces $S$ due to
increasing array and multiplexing gains. In both figures, a significant
improvement is observed compared to the single-layer BSIM and CSIM,
respectively. However, this improvement in the BSIM is greater.

(a) Achievable sum SE of the proposed double-SIM mMIMO architecture versus the
number of metasurfaces of the BSIM $L$ while varying the number of
metasurfaces of the CSIM $S$. Figure 4:

In Fig. 7, we show the sum SE versus the number of users $K$. We observe that
the sum SE grows with $K$ at the beginning but decreases after a certain value
of $K$. The reason behind this observation is that the increase in the
beginning comes from the increase of the spatial multiplexing gain with $K$.
Later, the increase of inter-user interference due to increasing $K$ cannot be
mitigated by the BSIM and its wave-based decoding. Based on our setup, the
BSIM can mitigate the multiuser interference of maximum $K=5$ users. Moreover,
we observe that an increase in the number of meta-atoms $M$ of the BSIM can
serve a larger number of users.

Figure 5: Achievable sum SE of the proposed double-SIM mMIMO architecture
versus the number of users $K$ while varying the number of BCSIM meta-atoms
$M$.

Fig. 6 depicts the sum SE versus the transmit SNR for various cases regarding
the values of the phase shifts of the BSIM. The optimal phase shifts design
achieves the best performance as expected. The designs with random phase
shifts and equal phase shifts follow with lower performance. In parallel,
Monte Carlo simulations verify the analytical results.

Figure 6: Achievable sum SE of the proposed double-SIM mMIMO architecture
versus the SNR for different cases regarding the phase shifts.

In Fig. 7, we present a comparison between the BCIM and CSIM-assisted mMIMO
system with the case where the BSIM is replaced by a conventional BS that
employs a conventional zero-forcing (ZF) decoder for mitigating the multiuser
interference. For the sake of a fair comparison, we assume that the path loss
corresponding to user $k$ is
$d_{k}=\sqrt{H_{\mathrm{BS}^{2}}+(d_{0}(k-1))^{2}}$. The conventional BS
requires $M_{\mathrm{BS}}=16$ antennas to exhibit better performance than the
BSIM with $M_{\mathrm{BS}}=64$ and $M=49$ elements. In other words, we can
remove completely the digital decoding at the BS while equipping each antenna
with low-resolution ADCs/DACs that can reduce overall the energy consumption
and hardware cost.

Figure 7: Achievable sum SE of the proposed double-SIM mMIMO architecture
versus the number of BS antennas $M_{\mathrm{BS}}$ while varying the number of
receiver metasurface layers $K$.

##  VI Conclusion

In this paper, a novel SIM-assisted mMIMO architecture was proposed. In
particular, a SIM enables wave-based decoding in the uplink, and another SIM
is used at the intermediate space between the users and the BS to increase
further the ability of the proposed configuration to reshape the surrounding
EM environment. We first suggested a channel estimation scheme, which obtains
the estimated channel in closed form in a single phase. Then, we derived the
uplink sum SE, and optimized both SIMs simultaneously. Simulations showed the
superiority of the proposed architecture compared to counterparts,
demonstrated the impact of various parameters on the sum SE, and exhibited
outperformance against the AO method, where the SIMs are optimized separately.
For example, an increase in the numbers of layers and meta-atoms increases the
performance. Especially, the BSIM exhibits interference cancellation. Also,
the proposed architecture is superior to the conventional digital MISO
counterpart, which means that SIM-based designs could lead to next-generation
energy-efficient networks.

##  Appendix A Proof of Lemma 1

The LMMSE estimator of
${\mathbf{c}}_{k}={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{h}}_{k}$
is derived according to according to [35, Ch. 12] as

$\displaystyle\hat{{\mathbf{c}}}_{k}=\mathbb{E}\\{{\mathbf{r}}_{k}{\mathbf{g}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}\left(\mathbb{E}\\{{\mathbf{r}}_{k}{\mathbf{r}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}\right)^{-1}.$
(27at)

The first term of (27at) is written as

$\displaystyle\mathbb{E}\\{{\mathbf{r}}_{k}{\mathbf{c}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}$
$\displaystyle=\mathbb{E}\big{[}\big{(}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{c}}_{k}+\frac{{\mathbf{z}}_{k}}{\tau\rho}\big{)}{\mathbf{c}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{W}}^{1}{\mathbf{P}}\big{]}$
(27au)
$\displaystyle={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}\mathbb{E}\\{{\mathbf{c}}_{k}{\mathbf{c}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}{\mathbf{W}}^{1}{\mathbf{P}}={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{R}}_{k}{\mathbf{W}}^{1}{\mathbf{P}}.$
(27av)

The second term is obtained as

$\displaystyle\mathbb{E}\\{{\mathbf{r}}_{k}{\mathbf{r}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}$
$\displaystyle={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{R}}_{k}{\mathbf{W}}^{1}{\mathbf{P}}+\frac{\sigma^{2}}{\tau\rho}{\bm{\mathrm{I}}}_{M}.$
(27aw)

By inserting (27av) and (27aw) into (27at), we obtain

$\displaystyle\hat{{\mathbf{c}}}_{k}=\hat{{\mathbf{R}}}_{k}{\mathbf{Q}}_{k}{\mathbf{r}}_{k},$
(27ax)

where
$\hat{{\mathbf{R}}}_{k}={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{R}}_{k}{\mathbf{W}}^{1}{\mathbf{P}}$
and
${\mathbf{Q}}_{k}\\!=\\!\left(\\!\hat{{\mathbf{R}}}_{k}\\!+\\!\frac{\sigma^{2}}{\tau\rho}{\bm{\mathrm{I}}}_{M}\\!\right)^{\\!-1}$.

##  Appendix B Proof of Theorem 1

Regarding the desired signal power part, we have

$\displaystyle\mathbb{E}\left\\{{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}\right\\}$
$\displaystyle=\tr(\mathbb{E}\\{{\mathbf{c}}_{k}\hat{{\mathbf{c}}}_{k}\\})$
$\displaystyle=\tr(\mathbb{E}\\{{\mathbf{c}}_{k}\hat{{\mathbf{R}}}_{k}{\mathbf{Q}}_{k}{\mathbf{r}}_{k}\\})$
(27ay) $\displaystyle=\tr({\bm{\Psi}}_{k}),$ (27az)

where, in (27ay), we have inserted (20), while in the last equation, we have
computed the expectation between ${\mathbf{c}}_{k}$ and ${\mathbf{r}}_{k}$,
and we have used (17).

In the case of $I_{k}$, the first term can be written as

$\displaystyle\\!\\!\mathbb{E}\left\\{|\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}-\mathbb{E}\left\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}\right\\}|^{2}\right\\}\\!=\\!\mathbb{E}\big{\\{}\big{|}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}\big{|}^{2}\big{\\}}\\!-\\!\big{|}\mathbb{E}\left\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}\right\\}\big{|}^{2}$
(27ba)
$\displaystyle=\mathbb{E}\big{\\{}\big{|}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}+\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\big{|}^{2}\big{\\}}-\big{|}\mathbb{E}\big{\\{}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\big{\\}}\big{|}^{2}$
(27bb)
$\displaystyle=\mathbb{E}\big{\\{}\big{|}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\big{|}^{2}\big{\\}}+2\real\\{\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\\}\\}$
(27bc)

where in (27bb), we have substituted (15). The second term in (27bc) goes to
zero according to the channel hardening property in mMIMO, which gives with
high accuracy
$\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\approx\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\\}$
[28]. Specifically, the second term becomes

$\displaystyle\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\\}\approx\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\\}\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\\}=0.$
(27bd)

Hence, (27bc) becomes

$\displaystyle\\!\\!\mathbb{E}\big{\\{}\big{|}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}-\mathbb{E}\big{\\{}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{k}\big{\\}}\big{|}^{2}\big{\\}}\\!\approx\\!\mathbb{E}\big{\\{}|\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}|^{2}\big{\\}}$
(27be)
$\displaystyle=\tr\\!\left(\hat{{\mathbf{R}}}_{k}{\bm{\Psi}}_{k}\right)-\tr\left({\bm{\Psi}}_{k}^{2}\right).$
(27bf)

For the derivation of (27bf), we apply the approximations
$\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\approx\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\\}$
and
$\tilde{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\approx\mathbb{E}\\{\tilde{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\\}$
due to channel hardening. Thus, we result in

$\displaystyle\mathbb{E}\\{|\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}|^{2}\\}$
$\displaystyle=\tr(\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{k}\tilde{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\})$
(27bg)
$\displaystyle\approx\tr(\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}\mathbb{E}\\{\tilde{{\mathbf{c}}}_{k}\tilde{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\})$
(27bh)
$\displaystyle=\tr((\hat{{\mathbf{R}}}_{k}-{\bm{\Psi}}_{k}){\bm{\Psi}}_{k}).$
(27bi)

The multiuser interference term in (23) becomes

$\displaystyle\mathbb{E}\left\\{|\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{c}}_{i}|^{2}\right\\}=\mathbb{E}\big{\\{}\big{|}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}(\hat{{\mathbf{c}}}_{i}+\tilde{{\mathbf{c}}}_{i})\big{|}^{2}\big{\\}}$
(27bj)
$\displaystyle=\mathbb{E}\big{\\{}\big{|}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{i}\big{|}^{2}\big{\\}}+\mathbb{E}\big{\\{}\big{|}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\tilde{{\mathbf{c}}}_{i}\big{|}^{2}\big{\\}}+2\real\\{\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{i}\tilde{{\mathbf{c}}}_{i}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\\}\\}$
(27bk)
$\displaystyle=\tr\\!\left({\bm{\Psi}}_{k}\hat{{\mathbf{R}}}_{i}\right),$
(27bl)

where the third term in (27bl) is zero because it holds that

$\displaystyle
2\real\\{\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{i}\tilde{{\mathbf{c}}}_{i}^{{\scriptscriptstyle\mathsf{H}}}\hat{{\mathbf{c}}}_{k}\\}\\}$
$\displaystyle\\!=\\!2\real\left\\{\mathbb{E}\\{\tr\\!\left(\\!\left(\hat{{\mathbf{c}}}_{k}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\right)\\!\\!\left(\tilde{{\mathbf{c}}}_{i}\hat{{\mathbf{c}}}_{i}^{{\scriptscriptstyle\mathsf{H}}}\right)\\!\right)\\!\right\\}$
(27bm)
$\displaystyle\\!=\\!2\real\left\\{\tr\\!\left(\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\}\mathbb{E}\\{\tilde{{\mathbf{c}}}_{i}\hat{{\mathbf{c}}}_{i}^{{\scriptscriptstyle\mathsf{H}}}\\}\right)\\!\right\\}$
(27bn) $\displaystyle=0,$ (27bo)

where, first, we used that $\hat{{\mathbf{c}}}_{k}$ and
$\hat{{\mathbf{c}}}_{i}$ are independent, and next that
$\tilde{{\mathbf{c}}}_{i}$ and $\hat{{\mathbf{c}}}_{i}$ are uncorrelated.

The last term in (23) is written as

$\displaystyle\mathbb{E}\left\\{|\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{n}}|^{2}\right\\}$
$\displaystyle=\tr(\mathbb{E}\\{\hat{{\mathbf{c}}}_{k}\hat{{\mathbf{c}}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\})$
(27bp) $\displaystyle=\tr({\bm{\Psi}}_{k}),$ (27bq)

where, in (27bp), we have used that $\hat{{\mathbf{c}}}_{k}$ and
${\mathbf{n}}$ are uncorrelated.

The combination of (27bf), (27bl), and (27bp) provides the approximate of
${I}_{k}$, which is $\tilde{I}_{k}$, and concludes the proof.

##  Appendix C Proof of Proposition 1

First, we focus on the derivation of
$\nabla_{{\bm{\phi}}_{l}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$. From (27a),
we can easily obtain

$\displaystyle\nabla_{{\bm{\phi}}_{l}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})=\frac{\tau_{c}-\tau}{\tau_{c}\log_{2}(e)}\sum_{k=1}^{K}\frac{\tilde{I}_{k}\nabla_{{\bm{\phi}}_{l}}S_{k}-S_{k}\nabla_{{\bm{\phi}}_{l}}\tilde{I}_{k}}{(1+\gamma_{k})\tilde{I}_{k}}.$
(27br)

For the computation of $\nabla_{{\bm{\phi}}_{l}}S_{k}$, we obtain its
differential as

$\displaystyle d(S_{k})$
$\displaystyle=d\bigl{(}\tr(\bm{\Psi}_{k})^{2}\bigr{)}$
$\displaystyle=2\tr(\bm{\Psi}_{k})d\tr(\bm{\Psi}_{k})$
$\displaystyle=2\tr(\bm{\Psi}_{k})\tr(d\bm{\Psi}_{k}).$ (27bs)

Application of [36, Eq. (3.35)] gives

$\displaystyle
d(\bm{\Psi}_{k})=d(\hat{\mathbf{R}}_{k}\mathbf{Q}_{k}\hat{\mathbf{R}_{k}})$
$\displaystyle=d(\hat{\mathbf{R}_{k}})\mathbf{Q}_{k}\hat{\mathbf{R}}_{k}+\hat{\mathbf{R}}_{k}d(\mathbf{Q}_{k})\hat{\mathbf{R}}_{k}+\hat{\mathbf{R}}_{k}\mathbf{Q}_{k}d(\hat{\mathbf{R}}_{k}).$
(27bt)

To this end, we have to derive $d(\hat{\mathbf{R}_{k}})$ and
$d(\mathbf{Q}_{k})$. The former can be written as

$\displaystyle d(\hat{\mathbf{R}}_{k})$
$\displaystyle={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}d({\mathbf{P}}^{H}){\mathbf{R}}_{k}{\mathbf{W}}^{1}{\mathbf{P}}+{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{R}}_{k}{\mathbf{W}}^{1}d({\mathbf{P}}).$
(27bu)

The latter is derived based on [36, eqn. (3.40)] as

$\displaystyle
d(\mathbf{Q}_{k})=d\bigl{(}\hat{\mathbf{R}}_{k}+\frac{\sigma}{\tau
P}\mathbf{I}_{M}\bigr{)}^{-1}$
$\displaystyle=-\bigl{(}\hat{\mathbf{R}}_{k}+\frac{\sigma^{2}}{\tau
P}\mathbf{I}_{M}\bigr{)}^{-1}d\bigl{(}\hat{\mathbf{R}}_{k}+\frac{\sigma^{2}}{\tau
P}\mathbf{I}_{M}\bigr{)}\bigl{(}\hat{\mathbf{R}}_{k}+\frac{\sigma}{\tau
P}\mathbf{I}_{M}\bigr{)}^{-1}$
$\displaystyle=-\mathbf{Q}_{k}d(\hat{\mathbf{R}}_{k})\mathbf{Q}_{k}.$ (27bv)

Substitution of (27bu) and (27bv) into (27bt) allows to obtain $d(S_{k})$ from
(27bs) as

where
${\mathbf{B}}_{k}=\mathbf{Q}_{k}\hat{\mathbf{R}}_{k}-\mathbf{Q}_{k}\hat{\mathbf{R}}^{2}_{k}\mathbf{Q}_{k}+\hat{\mathbf{R}}_{k}\mathbf{Q}_{k}$,
${\mathbf{A}}_{l}={\bm{\Phi}}^{L}{\mathbf{W}}^{L}\cdots{\bm{\Phi}}^{l+1}{\mathbf{W}}^{l+1}$,
and
${\mathbf{C}}_{l}={\mathbf{W}}^{l}{\bm{\Phi}}^{l-1}{\mathbf{W}}^{l-1}\cdots{\bm{\Phi}}^{1}$.
In (C), we have also used (2) by writing its differential as
$d({\mathbf{P}})={\mathbf{A}}_{l}{\bm{\Phi}}_{l}{\mathbf{C}}_{l}$

It can be easily obtained from (C) that

Regarding $\nabla_{{\bm{\phi}}_{l}}\tilde{I}_{k}$, we obtain the differential
from (26) as

$\displaystyle
d(\tilde{I}_{k})=\tr(\bm{\Psi}d(\hat{\mathbf{R}}_{k}))+\sum\nolimits_{i=1}^{K}\tr(\hat{\mathbf{R}}_{k}d(\bm{\Psi}_{i}))+\tr\bigl(\bar{\bm{\Psi}}_{k}d(\bm{\Psi}_{k})\bigr{missing}),$
(27bw)

where $\bm{\Psi}=\sum\nolimits_{i=1}^{K}\bm{\Psi}_{i}$,
$\bar{\bm{\Psi}}_{k}=\frac{1}{\rho}{\bm{\mathrm{I}}}-2\bm{\Psi}_{k}$. Having
obtained the differentials in (27bw) previously and after several algebraic
manipulations, we obtain

$\displaystyle\nabla_{{\bm{\phi}}_{l}}\tilde{I}_{k}$
$\displaystyle=\frac{\partial}{\partial{\bm{\phi}}_{l}^{*}}\tilde{I}_{k}$
(27bx)
$\displaystyle=\text{diag}\left(\right)\big{(}{\mathbf{A}}_{l}^{{\scriptscriptstyle\mathsf{H}}}\big{(}{\mathbf{R}}_{k}{\mathbf{W}}^{1}{\mathbf{P}}\big{(}{{\bm{\Psi}}}+\bar{{\bm{\Psi}}}_{k}{\mathbf{B}}_{k}$
$\displaystyle+\sum\nolimits_{i=1}^{K}{\mathbf{R}}_{i}{\mathbf{W}}^{1}{\mathbf{P}}{\mathbf{B}}_{i}\big{)}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{C}}_{l}^{{\scriptscriptstyle\mathsf{H}}}\big{)}.$
(27by)

In the case of
$\nabla_{{\bm{\lambda}}_{s}}f({\bm{\phi}}_{l},{\bm{\lambda}}_{s})$, we obtain
similar expressions, but now $d(\hat{\mathbf{R}}_{k})$ is derived as

$\displaystyle
d(\hat{\mathbf{R}}_{k})={\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}d({\mathbf{R}}_{k}){\mathbf{W}}^{1}{\mathbf{P}}$
$\displaystyle=\hat{\beta}_{k}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}$R_BSIMW^1P
$\displaystyle=\hat{\beta}_{k}{\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}$R_BSIMW^1P,

where, in (27bk), we have substituted (12). We have also denoted
${\mathbf{D}}_{s}={\mathbf{U}}^{s-1}{\bm{\Lambda}}^{s-1}\cdots{\mathbf{U}}^{2}{\bm{\Lambda}}^{1}$
and
${\mathbf{F}}_{s}={\bm{\Lambda}}^{S}{\mathbf{U}}^{S}\cdots{\bm{\Lambda}}^{s+1}{\mathbf{U}}^{s+1}$.
In (C), we have used (4) by writing its differential as
$d({\mathbf{Z}})={\mathbf{D}}_{s}{\bm{\Lambda}}_{s}{\mathbf{F}}_{s}$.

Substitution of (C) and (27bv) into (27bt) allows to obtain $d(S_{k})$ from
(27bs) as

Thus, we obtain

$\displaystyle\nabla_{{\bm{\lambda}}_{s}}S_{k}$
$\displaystyle=2\hat{\beta}_{k}\tr(\bm{\Psi}_{k})\tr({\mathbf{W}}^{1^{{\scriptscriptstyle\mathsf{H}}}}{\mathbf{P}}^{H}{\mathbf{R}}_{\mathrm{BSIM}}{\mathbf{W}}^{1}{\mathbf{P}}{\mathbf{B}}_{k})$
$\displaystyle\times\text{diag}\left(\right)({\mathbf{F}}_{s}\bar{{\mathbf{R}}}_{\mathrm{CSIM}}{\mathbf{Z}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{R}}_{\mathrm{CSIM}}{\mathbf{D}}_{s})^{{\scriptscriptstyle\mathsf{H}}}$

Similar to (27by), we obtain from (27bw) and (C) and after some algebraic
manipulations

$\displaystyle\nabla_{{\bm{\lambda}}_{s}}\tilde{I}_{k}$ $\displaystyle=2$
$\displaystyle\times\text{diag}\left(\right)({\mathbf{F}}_{s}\bar{{\mathbf{R}}}_{\mathrm{CSIM}}{\mathbf{Z}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{R}}_{\mathrm{CSIM}}{\mathbf{D}}_{s})^{{\scriptscriptstyle\mathsf{H}}}$

which concludes the proof.

## References

  * [1] F. Boccardi _et al._ , “Five disruptive technology directions for 5G,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 74–80, 2014.
  * [2] J. Zhang _et al._ , “Prospective multiple antenna technologies for beyond 5G,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1637–1660, 2020\.
  * [3] F. Sohrabi and W. Yu, “Hybrid digital and analog beamforming design for large-scale antenna arrays,” _IEEE J. Sel. Top. Signal Proc._ , vol. 10, no. 3, pp. 501–513, 2016.
  * [4] T. S. Rappaport _et al._ , “Wideband millimeter-wave propagation measurements and channel models for future wireless communication system design,” _IEEE Trans. Commun._ , vol. 63, no. 9, pp. 3029–3056, 2015.
  * [5] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [6] E. Basar _et al._ , “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, 2019\.
  * [7] E. Björnson and L. Sanguinetti, “Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 4, pp. 830–834, 2021.
  * [8] A. Papazafeiropoulos _et al._ , “Intelligent reflecting surface-assisted MU-MISO systems with imperfect hardware: Channel estimation and beamforming design,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 3, pp. 2077–2092, 2021.
  * [9] A. Papazafeiropoulos, “Ergodic capacity of IRS-assisted MIMO systems with correlation and practical phase-shift modeling,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 2, pp. 421–425, 2022.
  * [10] C. Huang _et al._ , “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Transa. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, 2019.
  * [11] M. Di Renzo _et al._ , “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [12] Q. U. A. Nadeem _et al._ , “Asymptotic max-min SINR analysis of reconfigurable intelligent surface assisted MISO systems,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 7748–7764, 2020.
  * [13] Y. Yang _et al._ , “Intelligent reflecting surface meets OFDM: Protocol design and rate maximization,” _IEEE Trans. Commun._ , vol. 68, no. 7, pp. 4522–4535, 2020.
  * [14] C. Pan _et al._ , “Multicell MIMO communications relying on intelligent reflecting surfaces,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 8, pp. 5218–5233, 2020.
  * [15] A. Papazafeiropoulos _et al._ , “Coverage probability of distributed IRS systems under spatially correlated channels,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 8, pp. 1722–1726, 2021.
  * [16] ——, “Achievable rate of a STAR-RIS assisted massive MIMO system under spatially-correlated channels,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2023.
  * [17] A. Papazafeiropoulos, P. Kourtessis, and S. Chatzinotas, “Max-Min SINR analysis of STAR-RIS assisted massive MIMO systems with hardware impairments,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2023.
  * [18] J. An _et al._ , “Stacked intelligent metasurfaces for efficient holographic MIMO communications in 6G,” _IEEE J. Sel. Areas Commun._ , 2023.
  * [19] X. Lin _et al._ , “All-optical machine learning using diffractive deep neural networks,” _Science_ , vol. 361, no. 6406, pp. 1004–1008, 2018.
  * [20] C. Liu _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.
  * [21] J. An _et al._ , “Stacked intelligent metasurfaces for multiuser downlink beamforming in the wave domain,” _arXiv preprint arXiv:2309.02687_ , 2023\.
  * [22] Q.-U.-A. Nadeem, J. An, and A. Chaaban, “Hybrid digital-wave domain channel estimator for stacked intelligent metasurface enabled multi-user MISO systems,” _arXiv preprint arXiv:2309.16204_ , 2023.
  * [23] S. Abeywickrama _et al._ , “Intelligent reflecting surface: Practical phase shift model and beamforming optimization,” _IEEE Trans. Commun._ , vol. 68, no. 9, pp. 5849–5863, 2020.
  * [24] Q. Nadeem _et al._ , “Intelligent reflecting surface-assisted multi-user MISO Communication: Channel estimation and beamforming design,” _IEEE Open J. Commun. Soc._ , vol. 1, pp. 661–680, 2020.
  * [25] N. V. Deshpande _et al._ , “Spatially-correlated irs-aided multiuser FD mMIMO systems: Analysis and optimization,” _IEEE Trans. Commun._ , vol. 70, no. 6, pp. 3879–3896, 2022.
  * [26] C. Wu _et al._ , “Channel estimation for STAR-RIS-aided wireless communication,” _IEEE Commun. Letters_ , vol. 26, no. 3, pp. 652–656, 2021\.
  * [27] B. Zheng _et al._ , “A survey on channel estimation and practical passive beamforming design for intelligent reflecting surface aided wireless communications,” _IEEE Commun. Sur. & Tut._, vol. 24, no. 2, pp. 1035–1071, 2022.
  * [28] E. Björnson _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.
  * [29] J. Hoydis, S. ten Brink, and M. Debbah, “Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?” _IEEE J. Select. Areas Commun._ , vol. 31, no. 2, pp. 160–171, 2013.
  * [30] A. K. Papazafeiropoulos and T. Ratnarajah, “Deterministic equivalent performance analysis of time-varying massive MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 10, pp. 5795–5809, 2015.
  * [31] A. K. Papazafeiropoulos, “Impact of general channel aging conditions on the downlink performance of massive MIMO,” _IEEE Trans. Veh. Tech._ , vol. 66, no. 2, pp. 1428–1442, Feb 2017.
  * [32] S. Zhang and R. Zhang, “Capacity characterization for intelligent reflecting surface aided MIMO communication,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1823–1838.
  * [33] N. S. Perović _et al._ , “Achievable rate optimization for MIMO systems with reconfigurable intelligent surfaces,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 6, pp. 3865–3882, 2021.
  * [34] D. Bertsekas, _Nonlinear Programming_ , 2nd ed., M. A. Scientific, Ed., 1999\.
  * [35] S. M. Kay, _Fundamentals of statistical signal processing: Estimation theory_. Upper Saddle River: Prentice Hall PTR, 1993.
  * [36] A. Hjørungnes, _Complex-Valued Matrix Derivatives: With Applications in Signal Processing and Communications_. Cambridge University Press, 2011.

solutionfile