text:
\section{Introduction}
Owing to their excellent features, the reconfigurable intelligent surfaces (RISs) have recently attracted a noticeable attention as a promising technique for future wireless communication networks.
An RIS is an artificial surface, made of electromagnetic
material, that is capable of customizing the propagation of
the radio waves impinging upon it \cite{Renzo1}.
It has been proposed as a new low-cost and less complicated solution to realize wireless
communication with high energy and spectrum efficiencies \cite{Basar1}.

In \cite{Wu}, it was shown that RIS outperforms
conventional massive multiple-input multiple-output
systems and multi-antenna
amplify-and-forward relaying networks with smaller
number of antennas, while reducing the system complexity and
cost. The performances of relay-assisted and RIS-assisted
wireless networks from coverage, probability of signal-to-noise ratio (SNR) gain, and delay outage rate aspects have been compared in \cite{Yang2}. In \cite{Yang4}, the authors utilized RIS to improve
the quality of a source signal that is sent to destination through an unmanned aerial vehicle.

The authors in \cite{Yang6} and \cite{Boulogeorgos} derived accurate approximations for channel distributions of RIS-assisted networks assuming Raleigh fading. Recently, some works on RIS-assisted networks over Nakagami-$m$ fading channels started to appear in literature \cite{Ferreira}, \cite{Anas1}. Most recently, limited number of papers have considered the interference phenomenon in a RIS context \cite{Li}, \cite{Hou}. All these works on interference effect considered the scenario, where an RIS is used as a replacement for relay.

There exists another important scenario where the RIS could be used as a part of the transmitter itself \cite{Basar1}, \cite{Yang6}. In such scenario, the transmitter includes an RF signal generator and an RIS, where the RIS is used as a transmitter along with the RF signal generator \cite{Basar1}. Recently, this idea was validated with the aid of a testbed platform \cite{Tang}. Combining RIS with relays is expected to achieve very high diversity order. In addition, addressing the impact of interference on the behavior of such scenario is of a noticeable importance.

In this letter, we consider the scenario, where the source and relays are internally utilizing RISs to enhance the quality of their transmitted signals in a decode-and-forward (DF) relaying network with interference. We derive accurate closed-form approximation for the system outage probability assuming Rayleigh fading channels for both intended users and interferers. The derived results are valid for arbitrary number of reflecting elements. We also derive closed-form expression for the asymptotic outage probability at high SNR values, where the system diversity order and coding gain are provided. To the best of authors' knowledge, the derived expressions are new and the achieved key findings on the impact of interference on the performance of RIS-assisted sources and relays network are being reported for the first time here.

\section{System and Channel Model}\label{SM}
Consider a dual-hop relay system with RIS-assisted source of $N_{1}$ reflecting elements, $K$ RIS-assisted DF relay nodes each of $N_{2}$ reflecting elements, one destination, and arbitrary number of interferers at both the relays and destination with opportunistic relaying. The entire communication takes place in two phases. In the first phase, the source {\sf S} transmits its signal to $K$ relays. In the second phase, only the best relay among all other relays who succeeded in decoding the source signal in the first phase is selected to forward it to destination {\sf D}. We assume that the signal at the $k^{\mathrm{th}}$ relay is corrupted by interfering signals from
$I_{k}$ co-channel interferers $\{x_{i}\}_{i=1}^{I_{k}}$.

With assuming that only the signals reflected by the RIS once are considered while those reflected two times and more are ignored, the received signal at the $k^{\mathrm{th}}$ relay can be expressed as
\begin{equation}\label{Eq.1}
y_{{\sf r}_{k}}= \sum_{i = 1}^{N_{1}} h_{{\sf s},k,i}e^{j\phi_{{\sf s},k,i}}x_{0}+\sum_{i_{k}=1}^{I_{k}}h_{ i_{k},k}^{I}x_{i_{k},k}^{I}+n_{{\sf
s},k},
\end{equation}
where $h_{{\sf s},k,i}=\alpha_{{\sf s},k,i}e^{-j\theta_{{\sf s},k,i}}$ is the channel coefficient between the $i^{\mathrm{th}}$ reflecting element at $\sf S$ and the $k^{\mathrm{th}}$ relay, where $\alpha_{{\sf s},k,i}$ is the amplitude and $\theta_{{\sf s},k,i}$ is the phase shift, $\phi_{{\sf s},k,i}$ is the adjustable phase induced by the $i^{\mathrm{th}}$ reflecting
element, $x_{0}$ is the transmitted symbol with $\mathbb{E}\{|{x_{0}}|^{2}\}=P_{0}$, $h_{i_{k},k}^{I}$ is the channel coefficient between the $i_{k}^{\mathrm{th}}$ interferer and $k^{\mathrm{th}}$ relay, $x_{i_{k},k}^{I}$ is the transmitted symbol from the $i_{k}^{\mathrm{th}}$ interferer with
$\mathbb{E}\{|{x_{i_{k},k}^{I}}|^{2}\}=P_{i_{k},k}^{I}$, $n_{{\sf s},k}\thicksim\mathcal{CN}(0, N_{0})$ is an additive white Gaussian noise (AWGN) of zero mean and power $N_{0}$, and $\mathbb{E}\{\cdot\}$ denotes the expectation operation. Let us define $h_{k,{\sf d},i}=\beta_{k,{\sf d},i}e^{-j\theta_{k,{\sf d},i}}$, where $\beta_{k,{\sf d},i}$ is the amplitude and $\theta_{k,{\sf d},i}$ is the phase shift, and ${h_{i_{d},\sf d}^{I}}$ as the channel coefficients between the $i^{\mathrm{th}}$ reflecting
element at the $k^{\mathrm{th}}$ relay and {\sf
D}, and the $i_{d}^{\mathrm{th}}$ interferer and {\sf D}, respectively. The amplitudes of the channels between {\sf S} and the $k^{\mathrm{th}}$ relay ${\alpha_{{\sf s},k,i}}, k=1...K$ and between the $k^{\mathrm{th}}$ relay and {\sf D} ${\beta_{k,{\sf d},i}}, k=1...K$ are assumed to be Rayleigh distributed with mean $\frac{\sqrt{\pi}}{2}$ and variance $\frac{(4-\pi)}{4}$. That is, their mean powers $\mathbb{E}\{|\alpha_{{\sf s},k,i}|^{2}\}=\mathbb{E}\{|\beta_{k,{\sf d},i}|^{2}\}=1$. In addition, the interferers channel coefficients are assumed to follow Rayleigh distribution. That is, the channel powers denoted by $|{h_{i_{k},k}^{I}}|^2$ and $|{h_{i_{d},\sf d}^{I}}|^2$ are
exponential distributed random variables (RVs) with parameters $\sigma_{I,i_{k},k}^{2}$ and
$\sigma_{I,i_{d},\sf d}^{2}$, respectively. Using (\ref{Eq.1}),
the signal-to-interference-plus-noise ratio (SINR) at the $k^{\mathrm{th}}$ relay can be written as
\begin{equation}
\gamma_{{\sf s},k}=\frac{P_{0}}{N_{0}}{\left(\sum_{i=1}^{N_{1}}|\alpha_{{\sf
s},k,i}|\right)^{2}}\Big/\left(\sum_{i_{k}=1}^{I_{k}}\frac{P_{i_{k},k}^{I}}{N_{0}}|{h_{
i_{k},k}^{I}}|^2+1\right).
\end{equation}

Let $B_{L}$ denote a decoding set defined by the set of relays who
successfully decoded the source message at the first phase. It is
defined as
\begin{align}\label{Eq.3}
B_{L}&\triangleq\left\{k\in \mathcal{S}_{r}:\gamma_{{\sf s},k}\geq
2^{2R}-1\right\},
\end{align}
where $\mathcal{S}_{r}$ is the set of all relays and $R$ denotes a
fixed spectral efficiency threshold.

In the second phase and after decoding the received signal, only the best relay in $B_{L}$ forwards the
re-encoded signal to the destination. The best
relay is the relay with the maximum
$\gamma_{l,\sf d}$, where $\gamma_{l,\sf d}$ is the SINR at the destination resulting from the $l^{\mathrm{th}}$ relay being the
relay, which forwarded the source information. It can be written as
\begin{equation}
\gamma_{l,\sf d}=\frac{P_{l}}{N_{0}}\left(\sum_{k=1}^{N_{2}}|{\beta_{l,{\sf
d},k}}|\right)^{2}\Big/\left(\sum_{i_{d}=1}^{I_{d}}\frac{P_{i_{d},\sf
d}^{I}}{N_{0}}\big|{h_{ i_{d},\sf d}^{I}}\big|^2+1\right),
\end{equation} where
$P_{l}$, $P_{i_{d},\sf d}^{I}$, and $N_{0}$ are the transmit power
of the $l^{\mathrm{th}}$ active relay, the transmit power of the $i_{d}^{\mathrm{th}}$ interferer, and the AWGN power at the
destination, respectively, and $I_{d}$ is the number of
interferers at the destination node. Since the denominator is
common to the SINRs from all relays belonging to $B_{L}$, the
best relay is the relay with the
maximum $\frac{P_{l}}{N_{0}}\left(\sum_{k=1}^{N_{2}}|{\beta_{l,{\sf
d},k}}|\right)^{2}$.

The end-to-end (e2e) SINR
at {\sf D} can be written as
\begin{equation}\label{dsinr}
\gamma_{\sf
d}=\frac{P_{b}}{N_{0}}\left(\sum_{k=1}^{N_{2}}|{\beta_{b,{\sf
d},k}}|\right)^{2}\bigg/\left(\sum_{i_{d}=1}^{I_{d}}\frac{P_{i_{d},\sf
d}^{I}}{N_{0}}\big|{h_{ i_{d},\sf d}^{I}}\big|^2+1\right),
\end{equation}
where the subscript $b$ is used to denote the best selected relay.
\section{Outage Performance Analysis}\label{PA}

\subsection{Preliminary Study}
The probability of the decoding set defined in (\ref{Eq.3}) can be
written as
\begin{align}\label{Eq.8}
\mathrm{P_{r}}\left[B_{L}\right]=\prod_{l\in
B_{L}}\mathrm{P_{r}}\left[\gamma_{{\sf s},l}\geq
u\right]\prod_{m\notin B_{L}}\mathrm{P_{r}}\left[\gamma_{{\sf
s},m}< u\right],
\end{align}
where $u=\left(2^{2R}-1\right)$ is the outage SNR threshold.
The outage probability of the system can be achieved by averaging
over all the possible decoding sets as follows \cite{JKim}
\begin{align}\label{Eq.9}
P_{\sf out}&\triangleq\mathrm{P_{r}}\left[\frac{1}{2}\
\mathrm{log}_{2}\left(1+\gamma_{\sf d}\right)<R\right]\nonumber\\
&=\sum_{L=0}^{K}\sum_{B_{L}}\mathrm{P_{r}}\left[\gamma_{\sf
d}<u|B_{L}\right]\mathrm{P_{r}}\left[B_{L}\right],
\end{align}
where the internal summation is taken over all of
${{K}\choose{L}}$ possible subsets of size $L$ from the set with
$K$ relays. In order to evaluate (\ref{Eq.9}), we need first to
derive $\mathrm{P_{r}}\left[\gamma_{\sf d}<u|B_{L}\right]$ and
$\mathrm{P_{r}}\left[B_{L}\right]$, which are presented in the
following section.

Throughout the analysis below, it assumed that $\rho=
P_{0}/N_{0}=P_{l}/N_{0}$ and $\rho_{I}=
P_{i_{k},k}^{I}/N_{0}=P_{i_{d},\sf d}^{I}/N_{0}$. The terms
$\rho_{I}|h_{i_{k},k}^{I}|^{2}$ and
$\rho_{I}|h_{i_{d},\sf d}^{I}|^{2}$ are exponential distributed
with parameters $\lambda_{i_{k},k}^{I}=1/\rho_{I}\sigma_{I,i_{k},k}^{2}$ and
$\lambda_{i_{d},\sf d}^{I}=1/\rho_{I}\sigma_{I,i_{d},\sf d}^{2}$. As the channels between {\sf S} and all relays and from relays to {\sf D} are assumed to have the same mean power, which equals 1, their average SNRs will be equal to $\rho=
P_{0}/N_{0}=P_{l}/N_{0}$ and their parameters will be equal to $\lambda_{{\sf s},k}=\lambda_{k,\sf d}=1/\rho, k=1,...,K$.

\subsection{Outage Probability}
The approximate outage probability of the considered system is summarized in the
following key result.

\begin{theorem}\label{theorem:1}
The outage probability of RIS-assisted source and relays DF relay network with independent identically distributed (i.i.d.) second hop cumulative distribution functions (CDFs) ($\lambda_{1,{\sf d}}=\lambda_{2,{\sf d}}=...=\lambda_{{\sf r,d}}$) and i.i.d. interferers' powers at {\sf D} ($\lambda^{I}_{1,{\sf d}}=\lambda^{I}_{2,{\sf d}}=...=\lambda^{I}_{I_{d},{\sf d}}=\lambda^{I}_{{\sf d}}$) and at relays ($\lambda^{I}_{1,k}=\lambda^{I}_{2,k}=...=\lambda^{I}_{I_{k},k}=\lambda^{I}_{k}, k=1,...,K$) can be obtained in a closed-form expression by using \eqref{Eq.9}, after evaluating the
terms $\mathrm{P_{r}}\left[\gamma_{\sf d}<u|B_{L}\right]$ and
$\mathrm{P_{r}}\left[\gamma_{{\sf s},k}<u\right]$ as follows
\end{theorem}
\begin{align}\label{Eq.27az}
&\mathrm{P_{r}}\left[\gamma_{\sf
d}<u|B_{L}\right]= - \frac{\left(\lambda_{\sf
 d}^{I} \right)^{I_{d}}}{\left(I_{d} - 1 \right)!}e^{\lambda_{\sf r,d}^{I}} \left(-1 \right)^{I_{d}}\sum_{g=0}^{I_{d} - 1}\ \begin{pmatrix}
I_{d} - 1 \\ g
\end{pmatrix} (-1)^{g} \nonumber\\&\times\sum_{k=0}^{L} \begin{pmatrix}
L \\ k \end{pmatrix} (-1)^{k}  \sum_{j_{0}=0}^{N_{2}-1}...\sum_{j_{k}=0}^{N_{2}-1} \frac{u^{\sum_{n=0}^{k}j_{n}}}{(C/\lambda_{\sf r,d})^{\sum_{n=0}^{k}j_{n}}\prod_{n=0}^{k}j_{n}!}\nonumber\\
&\times\left( \lambda^{I}_{\sf d}+ u\lambda_{\sf r,d} k \right)^{-g - \sum_{n=0}^{k}j_{n} - 1} \Gamma\left(g+\sum_{n=0}^{k}j_{n}+1, \lambda^{I}_{\sf d}+ \frac{u\lambda_{\sf r,d} k}{C}\right),
\end{align}
\begin{align}\label{Eq.27bz}
&\mathrm{P_{r}}\left[\gamma_{{\sf
s},k}<u\right]\approx \frac{\left(\lambda_{k}^{I} \right)^{I_{k}}}{\left(I_{k} - 1 \right)!}e^{\lambda_{k}^{I}} \left(-1 \right)^{I_{k}}\sum_{g=0}^{I_{k} - 1}\ \begin{pmatrix}
I_{k} - 1 \\ g
\end{pmatrix} (-1)^{g} \nonumber\\
 &\times \bBigg@{3.2}(- (\lambda^{I}_{k})^{-g-1} \Gamma\left(g+1,\lambda^{I}_{k}\right)+ \sum_{i = 0}^{N_{1} - 1} \frac{u^{i}  (\lambda_{{\sf s},k})^{i}}{B^i i!\left(\lambda^{I}_{k} + \frac{u  \lambda_{{\sf s},k}}{B}\right)^{g+i+1}} \nonumber\\
 &\times\Gamma\left(g+i+1,\lambda^{I}_{k} + \frac{u \lambda_{{\sf s},k}}{B }\right)\bBigg@{3.2}),
\end{align}
where $\Gamma(.,.)$ is the upper incomplete gamma function \cite[Eq. (8.352.2)]{Grad.}.

\noindent\hspace{2em}{\itshape Proof: }  See Appendix \ref{appendix:1}.\hspace*{\fill}~$\blacksquare$\par\endtrivlist\unskip

\section{Asymptotic Outage Behavior}\label{AOB}
Here, $P_{\sf out}$ can be written as $P_{\sf
out}\approx\left(G_{c}\rho\right)^{-G_{d}}$, where $G_{c}$ is the coding gain and $G_{d}$ is the diversity order \cite{Alouinib}.

\begin{theorem}\label{Theorem:2}
The asymptotic outage probability for RIS-assisted source and relays DF relay network with interference can be obtained in a closed-form expression by using \eqref{Eq.9}, after evaluating the
terms $\mathrm{P_{r}}\left[\gamma_{\sf d}<u|B_{L}\right]$,
$\mathrm{P_{r}}\left[\gamma_{{\sf s},k}<u\right]$ as follows
\begin{align}\label{Eq.28a}
&\mathrm{P_{r}}\left[\gamma_{\sf
d}<u|B_{L}\right]\approx - \frac{\left(\lambda_{\sf
 d}^{I} \right)^{I_{d}}e^{\lambda_{\sf d}^{I}}\left(-1 \right)^{I_{d}}}{\left(I_{d} - 1 \right)!C^{N_{2}L}(N_{2}!)^{L}} \sum_{g=0}^{I_{d} - 1}\ \begin{pmatrix}
I_{d} - 1 \\ g
\end{pmatrix}\nonumber\\
&\times (-1)^{g} (\lambda_{\sf r,d})^{N_{2}L} u^{N_{2}L} (\lambda_{\sf d}^{I})^{-g-N_{2}L-1} \Gamma\left(g+N_{2}L+1, \lambda^{I}_{\sf d}\right),
\end{align}
\begin{align}\label{Eq.28b}
&\mathrm{P_{r}}\left[\gamma_{{\sf
s},k}<u\right]\approx - \frac{\left(\lambda_{k}^{I} \right)^{I_{k}}e^{\lambda_{k}^{I}} \left(-1 \right)^{I_{k}}}{\left(I_{k} - 1 \right)!B^{N_{1}}N_{1}!}\sum_{g=0}^{I_{k} - 1}\ \begin{pmatrix}
I_{k} - 1 \\ g
\end{pmatrix} (-1)^{g} \nonumber\\&\times (\lambda_{{\sf s},k})^{N_{1}} u^{N_{1}} (\lambda_{k}^{I})^{-g-N_{1}-1} \Gamma\left(g+N_{1}+1, \lambda^{I}_{k}\right).
\end{align}
\end{theorem}
\noindent\hspace{2em}{\itshape Proof: } See Appendix \ref{appendix:3}.\hspace*{\fill}~$\blacksquare$\par\endtrivlist\unskip

Two extreme cases have been noticed here as follows.
\begin{itemize}
\item Case 1: $N_{1}<N_{2}$\\
Here, $P_{\sf out}$ is dominated by the first hop as follows
\begin{align}\label{Eq.28bss}
P_{\sf out}^{\infty}&\approx \bBigg@{3.2}(\bBigg@{3.2}[- \frac{\left(\lambda_{\sf r}^{I} \right)^{I_{r}}e^{\lambda_{\sf r}^{I}} \left(-1 \right)^{I_{r}}}{\left(I_{r} - 1 \right)!B^{N_{1}}N_{1}!}\sum_{g=0}^{I_{r} - 1}\ \begin{pmatrix}
I_{r} - 1 \\ g
\end{pmatrix} (-1)^{g}\nonumber\\
&\times u^{N_{1}} \frac{\Gamma\left(g+N_{1}+1, \lambda^{I}_{\sf r}\right)}{(\lambda_{\sf r}^{I})^{g+N_{1}+1}}\bBigg@{3.2}]^{-1/N_{1}}\rho\bBigg@{3.2})^{-N_{1}K}.
\end{align}
Clearly, the diversity order of the system for this case is $G_{d}=N_{1}K$ and the coding gain $G_{c}$ is a function of several parameters, including $N_{1}$, $I_{r}$, $\lambda_{\sf r}^{I}$, and $u$.

\item Case 2: $N_{2}<N_{1}$\\
Here, $P_{\sf out}$ is dominated by the second hop as follows
\begin{align}\label{Eq.28ass}
P_{\sf out}^{\infty}&\approx \bBigg@{3.2}(\bBigg@{3.2}[ - \frac{\left(\lambda_{\sf
 d}^{I} \right)^{I_{d}}e^{\lambda_{\sf d}^{I}}\left(-1 \right)^{I_{d}}u^{N_{2}L}}{\left(I_{d} - 1 \right)!C^{N_{2}L}(N_{2}!)^{L}} \sum_{g=0}^{I_{d} - 1}\ \frac{\begin{pmatrix}
I_{d} - 1 \\ g
\end{pmatrix}}{(-1)^{-g}} \nonumber\\
&\times \frac{\Gamma\left(g+N_{2}L+1, \lambda^{I}_{\sf d}\right)}{(\lambda_{\sf d}^{I})^{g+N_{2}L+1}}\bBigg@{3.2}]^{-1/N_{2}L}\rho\bBigg@{3.2})^{-N_{2}L}.
\end{align}
The second hop is dominant when the number of active relays is equal to number of all available relays, that is $L=K$. This means that the diversity order of this system for this case is equal to $G_{d}=N_{2}K$, whereas it is clear from \eqref{Eq.28ass} that the coding gain $G_{c}$ is a function of several parameters, including $N_{2}$, $I_{d}$, $\lambda_{\sf d}^{I}$, and $u$.
\end{itemize}

Therefore, the diversity order of the system when $N_{1}$ is not equal to $N_{2}$ is $G_{d}=\mathrm{min}(N_{1},N_{2})K$. On the other hand, when $N_{1}=N_{2}=N$, the interference value at the relays and destination specifies, which hop dominates the system performance. The hop that suffers from higher interference dominates the system performance by determining its coding gain. When the interference at both hops is equal, we noticed that the second hop dominates the system performance through determining its coding gain as in Case 2 above. For the last two cases, the diversity order is $G_{d}=NK$.
\section{Numerical Results}\label{NR}
We assume here that the parameters $\sigma_{k,\sf d}^{2}$, $\sigma_{I,i_{k},k}^{2}$, and
$\sigma_{I,i_{d},\sf d}^{2}$ are all equal to 1. In addition, we call the outage SNR threshold as $\gamma_{\sf out}=u$. A good matching between the analytical and asymptotic results and simulations is clear in Fig. \ref{Pout_SNR_N1_K}. In addition, the diversity order of the system $G_{d}$ is increasing when both $N_{1}$ and $K$ increase. As the system performance is dominated here by the first hop, $G_{d}=N_{1}K$ and this coincides with Case 1 in the asymptotic analysis section.
\begin{figure}[htb!]
\centering
\includegraphics[scale=0.34]{Pout_SNR_N1_K.eps}
\caption{$P_{\sf out}$ vs SNR for various values of $N_{1}$ and $K$.}\label{Pout_SNR_N1_K}
\end{figure}

Same insights can be drawn from Fig. \ref{Pout_SNR_N2_K}, but here $G_{d}=N_{2}K$ as the performance is dominated by the second hop. Results match Case 2 in the asymptotic analysis section.
\restylefloat{figure}
\begin{figure}[htb!]
\centering
\includegraphics[scale=0.34]{Pout_SNR_N2_K.eps}
\caption{$P_{\sf out}$ vs SNR for various values of $N_{2}$ and $K$.}\label{Pout_SNR_N2_K}
\end{figure}

In Fig. \ref{Pout_SNR_K_N2}, as $N_{2}$ is smaller than $N_{1}$, $G_{d}=3$. This figure informs us that the system with one relay node and multiple reflecting elements ($K=1,N_{2}=3$), gives better performance than the system with multiple relays each with one reflecting element ($K=3,N_{2}=1$).
\begin{figure}[htb!]
\centering
\includegraphics[scale=0.34]{Pout_SNR_K_N2.eps}
\caption{$P_{\sf out}$ vs SNR for various values of $K$, $N_{2}$, $\bar{\gamma}_{\sf r}^{I}$, and $\bar{\gamma}_{\sf d}^{I}$.}\label{Pout_SNR_K_N2}
\end{figure}

Fig. \ref{Pout_SNR_gIr_gId} shows two main cases, (solid line: $N_{1}=4, N_{2}=1$) where the behavior is dominated by the second hop and (dash line: $N_{1}=1, N_{2}=4$) where the behavior is dominated by the first hop. For all curves here, $G_{d}=2$. For each main case, decreasing the interference power enhances $G_{c}$ and the system performance. In addition, this figure informs us that the interference at the destination is more impactful/severe on the system performance than the interference at the relays. This is clear as the case where ($\bar{\gamma}_{\sf r}^{I}=10\ \mathrm{dB},\bar{\gamma}_{\sf d}^{I}=20\ \mathrm{dB}$) gives worse results than the case where ($\bar{\gamma}_{\sf r}^{I}=20\ \mathrm{dB},\bar{\gamma}_{\sf d}^{I}=10\ \mathrm{dB}$).
\begin{figure}[htb!]
\centering
\includegraphics[scale=0.34]{Pout_SNR_gIr_gId.eps}
\caption{$P_{\sf out}$ vs SNR for various values of $\bar{\gamma}_{\sf r}^{I}$ and $\bar{\gamma}_{\sf d}^{I}$.}\label{Pout_SNR_gIr_gId}
\end{figure}

Fig. \ref{Pout_SNR_N1_N2} shows two cases, first hop dominating the system performance ($N_{1}<N_{2}$) and second hop dominating the system performance ($N_{2}<N_{1}$). For each case, increasing the number of reflecting elements $N_{1}$ or $N_{2}$ increases $G_{d}$ and enhances the performance. For the first case, $G_{d}=N_{1}K$ and for the second case, $G_{d}=N_{2}K$. Also, having more reflecting elements at the relays than the source gives better performance. This is because when the first hop is dominating the performance, the interference power at the relays determines $G_{c}$, which is better than that when the second hop is dominating and the interference at the destination is more severe on the behavior.
\restylefloat{figure}
\begin{figure}[htb!]
\centering
\includegraphics[scale=0.34]{Pout_SNR_N1_N2.eps}
\caption{$P_{\sf out}$ vs SNR for various values of $N_{1}$ and $N_{2}$.}\label{Pout_SNR_N1_N2}
\end{figure}

Fig. \ref{Pout_SNR_N1_N2_gIr_gId_Ir_Id} studies the case when $N_{1}=N_{2}$, where the interference power is determining which hop dominates the system performance. When $N_{1}=N_{2}=15$, $G_{d}=30$, whereas when $N_{1}=N_{2}=20$, $G_{d}=40$. Clearly, when the interference power is reduced at the relays, the enhancement on $G_{c}$ is unnoticeable, whereas when the interference power is reduced at the destination, $G_{c}$ increases and the system gain becomes noticeable. This behavior is expected as the interference at the destination is more impactful on the performance than the relays interference. We can also see that when the first hop is dominating the performance ($\bar{\gamma}^{I}_{\sf r}=15\ \mathrm{dB}, \bar{\gamma}^{I}_{\sf d}=5\ \mathrm{dB}$), decreasing the number of interferers at the relays from 3 to 1 has a noticeable impact on $G_{c}$. Same thing applies when the second hop is dominating the performance ($\bar{\gamma}^{I}_{\sf r}=5\ \mathrm{dB}, \bar{\gamma}^{I}_{\sf d}=15\ \mathrm{dB}$), where decreasing the number of interferers at the destination from 3 to 1 also noticeably enhances $G_{c}$.
\restylefloat{figure}
\begin{figure}[htb!]
\centering
\includegraphics[scale=0.34]{Pout_SNR_N1_N2_gIr_gId_Ir_Id.eps}
\caption{$P_{\sf out}$ vs SNR for various values of $N_{1}$ and $N_{2}$ and various interference parameters.}\label{Pout_SNR_N1_N2_gIr_gId_Ir_Id}
\end{figure}

\section{Conclusion}\label{C}
This letter considered the scenario of DF relay network with RIS-assisted source and relays with interference. It derived approximate closed-form and asymptotic expressions for the outage probability assuming Rayleigh fading channels. Results showed that the system can achieve a diversity order of $G_{d}=\mathrm{min}(N_{1},N_{2})K$. Findings illustrated that for the same $G_{d}$, utilizing one relay with multiple reflecting elements gives better performance than utilizing multiple relays with a single reflecting element. The interference at the destination is more severe on the performance than the interference at the relays. Therefore, under the same interference and for a given $K$, the case $N_{1}<N_{2}$ outperforms the case $N_{2}<N_{1}$ in terms of $G_{c}$.

\appendices
\section{Proof of Theorem \ref{theorem:1}}\label{appendix:1}
In this Appendix, we evaluate the first term in (\ref{Eq.9})
$\mathrm{P_{r}}\left[\gamma_{\sf d}<u|B_{L}\right]$. First, the e2e SINR
can be written as $\gamma_{\sf d} = Y/Z$. The CDF of
$\gamma_{\sf d}$ given a decoding set $B_{L}$ is given by
$\mathrm{P_{r}}\left[\gamma_{\sf
d}<u|B_{L}\right]=\int_{1}^{\infty}f_{Z}(z)\int_{0}^{uz}f_{Y}(y)dydz = \int_{1}^{\infty}f_{Z}(z)F_{Y}(uz)dz$. With RIS-assisted relays of $N_{2}$ reflecting elements, the CDF of the second hop SNR is given by \cite{Yang6}
\begin{align}\label{sese}
F_{\gamma_{\sf r,d}}(\gamma) = 1- e^{-\frac{\lambda_{\sf r,d} \gamma}{C}} \sum_{k = 0}^{N_{2} - 1} \frac{(\lambda_{\sf r,d} \gamma)^k}{C^k k!},
\end{align}
where $C = 1 + (N_{2} - 1)\Gamma^{2}\left(\frac{3}{2} \right)$. Upon utilizing the Binomial rule, the CDF of the best selected relay can be rewritten as
\begin{align*}
&F_{Y}(uz)= (F_{\gamma_{\sf r,d}}(uz))^L=\left(1- e^{-\frac{\lambda_{\sf r,d} uz}{C}} \sum_{k = 0}^{N_{2} - 1} \frac{(\lambda_{\sf r,d} uz)^k}{C^k k!}\right)^L\\
&= \sum_{k=0}^{L} \begin{pmatrix}
L \\ k \end{pmatrix} (-1)^{k} \exp\left(-\lambda_{\sf r,d}u z k\right) \sum_{j_{1}=0}^{N_{2}-1}...\sum_{j_{k}=0}^{N_{2}-1}u^{\sum_{n=1}^{k}j_{n}} \nonumber\\
&\times \frac{1}{(C/\lambda_{\sf r,d})^{\sum_{n=1}^{k}j_{n}}\prod_{n=1}^{k}j_{n}!} .
\end{align*}
With Rayleigh distributed interference, the interference-to-noise ratio (INR) CDFs at {\sf D} will be following exponential distribution as $F_{\gamma_{\sf d}^{I}}(\gamma)=1-e^{-\lambda_{\sf d}^{I}\gamma}$. Hence, the RV $Z$, which equal to $Z=X+1$, where $X$ is the interference at the destination has the following PDF \cite{Anas}
\begin{align}\label{Eq.13w}
f_{X}(x)= - \frac{\left(\lambda_{\sf
 d}^{I} \right)^{I_{d}}}{\left(I_{d} - 1 \right)!}x^{I_{d}-1}e^{-\lambda_{\sf d}^{I}x}.
\end{align}
Using transformation of RVs and then the binomial rule, the PDF of $Z$ can be obtained as
\begin{align}\label{Eq.13}
f_{Z}(z)= - \frac{\left(\lambda_{\sf
 d}^{I} \right)^{I_{d}}}{\left(I_{d} - 1 \right)!}\frac{\left(-1 \right)^{I_{d}}}{e^{-\lambda_{\sf d}^{I}}}\sum_{g=0}^{I_{d} - 1}\ \begin{pmatrix}
I_{d} - 1 \\ g
\end{pmatrix} \frac{(-1)^{g} z^{g}}{e^{\lambda_{\sf d}^{I}z}}.
\end{align}
Now, using the integral $\int_{1}^{\infty}f_{Z}(z)F_{Y}(uz)dz$ and with the help of \cite[Eq. (8.351.2)]{Grad.}, the first term in \eqref{Eq.9} $\mathrm{P_{r}}\left[\gamma_{\sf d}<u|B_{L}\right]$ can be obtained as in \eqref{Eq.27az}.

To obtain the second term in (\ref{Eq.9}) $\mathrm{P_{r}}\left[\gamma_{{\sf s},k}<u\right]$, the first hop SINR can be written as $Y'/Z'$, where $Y'$ is the first hop SNR with RIS-aided transmitter and $Z'=X'+1$, where $X'$ is the interference at the relays. Again, the CDFs of the interferers at the relays follow exponential distribution as $F_{\gamma_{k}^{I}}(\gamma)=1-e^{-\lambda_{k}^{I}\gamma}, k=1...K$. With RIS-assisted source of $N_{1}$ reflecting elements, the CDF of first hop SNR is given by \cite{Yang6}
\begin{align}
F_{Y'}(uz) = 1- e^{-\frac{\lambda_{{\sf s},k} u z}{B}} \sum_{i = 0}^{N_{1} - 1} \frac{(\lambda_{{\sf s},k} u z)^i}{B^i i!},
\end{align}
where $B = 1 + (N_{1} - 1)\Gamma^{2}\left(\frac{3}{2} \right)$.\\
The PDF of $Z'$ is similar to that of $Z$, but with replacing $\lambda^{I}_{\sf d}$ with $\lambda^{I}_{k}$ and $I_{d}$ with $I_{k}$. Again, upon using the integral $\int_{1}^{\infty}f_{Z'}(z)F_{Y'}(uz)dz$ and with the help of \cite[Eq. (8.351.2)]{Grad.}, the second term in \eqref{Eq.9} $\mathrm{P_{r}}\left[\gamma_{{\sf s},k}<u\right]$ can be obtained as in \eqref{Eq.27bz}.

\section{Proof of Theorem \ref{Theorem:2}}\label{appendix:3}
To find the asymptotic outage probability, we first need to
obtain $\mathrm{P_{r}}\left[\gamma_{\sf d}<u|B_{L}\right]$. As
$\rho\rightarrow\infty$ and with constant values of $\rho_{I}$,
$I_{k}$, and $I_{d}$, the CDF of the second hop channels with opportunistic relaying can be approximated as
\begin{align}\label{C.1}
F_{\gamma_{\sf r,d}}(\gamma)\approx
\frac{\gamma^{N_{2}}}{(C/ \lambda_{\sf r,d})^{N_{2}} N_{2}!}.
\end{align}
Now, the CDF of best selected relay can be written as
\begin{align}\label{C.1w}
F_{Y}(u z)\approx
\frac{(u z)^{LN_{2}}}{(C/ \lambda_{\sf r,d})^{LN_{2}} (N_{2}!)^{L}}.
\end{align}
With no change in the PDF of $Z$ and upon following the same procedure as in Appendix \ref{appendix:1},
the term $\mathrm{P_{r}}\left[\gamma_{\sf d}<u|B_{L}\right]$ in
(\ref{Eq.9}) can be evaluated at high SNR values as in \eqref{Eq.28a}.

Now, the second term in \eqref{Eq.9}
$\mathrm{P_{r}}\left[B_{L}\right]$ can be obtained after
evaluating the CDF of $\gamma_{{\sf s},k}$, which can be
approximated at high SNR values as
\begin{align}\label{C.1}
F_{\gamma_{{\sf s},k}}(u z)\approx
\frac{(u z)^{N_{1}}}{(B/ \lambda_{{\sf s},k})^{N_{1}} N_{1}!}.
\end{align}
Again, with no change in the PDF of $Z'$ and upon following the same procedure as in Appendix \ref{appendix:1},
the term $\mathrm{P_{r}}\left[B_{L}\right]$ in
(\ref{Eq.9}) can be evaluated at high SNR values as in \eqref{Eq.28b}.

\balance


text:
\section{Introduction}
\label{s1}

Acceleration of particles in plasmas, or plasma wakefield acceleration, offers the possibility to reach multi-GeV and, potentially, TeV range electron and positron energies in facilities that are orders of magnitude smaller than modern high-energy accelerators \cite{Nat.445-741,Nat.524-442,PPCF58-034017, PRL122-084801,RMP81-1229,RAST9-19,RAST9-63,NatPhys5-363,PoP18-103101,RAST9-85}.
Studies of novel accelerators significantly benefit from numerical simulations that go in parallel with them \cite{RAST9-165}. These simulations complement fragmentary experimental data \cite{RMP90-035002} and form a complete picture of physical processes that occur at tiny spatial and temporal scales inside the plasma and, therefore, cannot be measured.

Recently, the Advanced WAKefield Experiment (AWAKE) at CERN \cite{NIMA-829-3,NIMA-829-76,PPCF60-014046} demonstrated seeded self-modulation of a long proton beam in the plasma \cite{PRL122-054801,PRL122-054802} and electron acceleration by the wakefield of this beam \cite{Nat.561-363}. This milestone achievement opens the way for using proton beams from modern synchrotrons as drivers for plasma wakefield acceleration \cite{RAST9-85,PPCF56-084013,EPJC76-463}. Because of the high proton energy, the micro-bunching and acceleration occurs in a single plasma cell, avoiding difficulties related to staged acceleration \cite{PRST-AB3-071301,PRST-AB15-111303,PRL116-124801,IPAC16-2561}.

In this paper, we show that simulations of beam self-modulation agree with related AWAKE measurements. The axisymmetric self-modulation process is thus well understood in the sense that simulations include the most important effects, they can help in understanding details of the process and can serve as a valuable starting point for future predictions. We also discuss how to conduct these simulations properly, and show what the results would look like otherwise. The AWAKE data used for the studies reported in this paper were taken during the 2017 and early 2018 running periods.

\section{Methods}

AWAKE uses the proton beam from the SPS \cite{IPAC17-1747}. The proton beam co-propagates with a short laser pulse through a 10-meter-long cylindrical cell filled with rubidium vapor (Fig.\,\ref{fig1}). The laser pulse quickly ionizes the vapor, so that the newly appearing plasma interacts only with the rear part of the beam, as if the beam would be cut at the pulse position. The leading edge of this ``sharply cut'' beam seeds a small amplitude wakefield in the plasma that grows in space and time and eventually converts the beam into a train of short micro-bunches \cite{PRL104-255003,PoP22-103110}. The self-modulation occurs due to defocusing of protons from the space between the micro-bunches. The radial momentum gained by the defocused protons is an integral quantitative measure of the self-modulation process: the stronger the longitudinal wakefield micro-bunch formation, the larger the transverse momentum of the defocused protons \cite{PRL122-054801,NIMA-829-314,NAPAC16-707}.

\begin{figure*}[t]
\includegraphics{fig1.pdf}
\caption{Setup of AWAKE experiments on beam self-modulation.}\label{fig1}
\end{figure*}

\begin{figure}[tb]
\includegraphics[width=\columnwidth]{Fig2.pdf}
\caption{Images of beam core (a) and halo (b) at the second (downstream) imaging station. The red solid line in (b) shows the detected halo boundary used for determination of the maximum deflection radius. Note that the intense light emitted by the core of the beam is blocked by a mask.}\label{fig2-images}
\end{figure}

\begin{figure}[tb]
\includegraphics{fig3-profile.pdf}
\caption{Angle-averaged radial distributions of the beam flux in several individual shots (grey lines), additionally averaged over 30 shots (thin dark gray line), and simulated: with optimum parameters from Table~\ref{t1} (solid red line), with a reduced angular spread and the same beam radius (dotted line) and with a reduced beam radius and the same angular spread (dashed line). The vertical black lines mark the maximum proton radius calculated with the contour method \cite{NIMA-909-123}.}\label{fig3-profile}
\end{figure}

The defocused protons form a halo around the micro-bunched part of the beam. Its time-integrated flux is measured at two imaging stations located 1.66\,m and 9.75\,m downstream the plasma cell exit \cite{IPAC17-1682,NIMA-854-100} (Fig.\,\ref{fig1}). Each imaging station separately measures the flux distribution in the halo and in the central part (core) of the beam  with two cameras of different sensitivity (Fig.\,\ref{fig2-images}). Separate registration of the core and halo, followed by stitching two parts of the image, provides a wide dynamic range of measurements (approximately $10^4$) and absolute flux calibration through equating the integral flux and the known beam charge (Fig.\,\ref{fig3-profile}). The halo images processed with a contour method \cite{NIMA-909-123} yield the maximum radius of defocused protons, which is a convenient scalar to be compared with simulations.

The simulations are performed with the quasistatic two-dimensional (axisymmetric) particle-in-cell code LCODE \cite{PRST-AB6-061301,NIMA-829-350}. The simulation window encloses the whole beam and moves with the speed of light $c$ in the positive $z$ direction. The radial window size is always $30c/\omega_p$, where $\omega_p = \sqrt{4 \pi n_0 e^2/m}$ is the plasma frequency, $n_0$ is the plasma density, $e$ is the elementary charge, and $m$ is the electron mass. The grid cell has equal sizes $0.01c/\omega_p$ in both longitudinal and radial directions. There are 10~radius-weighted plasma macro-particles of each species (electrons and ions) in each cell and about $10^6$ or more (proportionally to $\omega_p$) equal beam macro-particles in total. The time step for the beam is $100\omega_p^{-1}$, and the plasma state is updated every $200c/\omega_p$. The plasma has sharp boundaries in both longitudinal and radial directions, as simulations \cite{JPD51-025203,NIMA-740-197,PRA99-063423} have suggested. After the plasma cell, the protons propagate to the screens along straight lines.

\begin{table}[b]
 \caption{Beam and plasma parameters taken as input for simulations. Parameters in the first group are the same in all runs. Parameters in the second group (below the line) may vary in different regimes, and the listed values are those for illustrative cases (Figs.\,\ref{fig3-profile}, \ref{fig5-wall}, \ref{fig6-expl}).}\label{t1}
 \begin{indented}
 \item[]\begin{tabular}{@{}ll}\br
  Parameter, notation & Value \\ \mr
  Plasma length, $L$ & 10.3\,m \\
  Plasma radius, $R_p$ & 1.5\,mm \\
  Plasma ion-to-electron mass ratio, $M_i$ & 157\,000\\
  Plasma initial temperature, $T_e$ & 0\,eV \\
  Beam energy, $W_b$ & 400\,GeV \\
  Beam energy spread, $\delta W_b$ & 0.035\,\% \\
  \mr
  Plasma density, $n_0$ \quad & $1.965\times 10^{14}\,\text{cm}^{-3}$ \\
  Beam population, $N_b$ & $2.86\times 10^{11}$ \\
  Beam length, $\sigma_z$ & 6.57\,cm \\
  Beam radius, $\sigma_r$ & $200\,\mu$m \\
  Normalized beam emittance, $\epsilon_b$ & 3.6\,mm\,mrad \\
  Seed pulse position, $\xi_s$ & 3.75\,cm \\
  \br
 \end{tabular}
 \end{indented}
\end{table}

\begin{figure}[tb]
\includegraphics{fig4-streak.pdf}
\caption{Streak camera image of the proton beam (a) and Gaussian approximation of the beam current (b).}\label{fig4-streak}
\end{figure}

Most of the input parameters for simulations (Table~\ref{t1}) are known with sufficient precision and their uncertainties have little effect on quantitative characteristics of beam self-modulation. Less well known are the beam length, radius and emittance and these are determined from the data used in this analysis. The proton beam was modelled with Gaussian distributions in both the longitudinal and transverse directions. The standard deviation of the proton beam in the longitudinal direction, $\sigma_z$, was determined prior to the start of data taking using a streak camera \cite{RSI88-025110}  (Fig.\,\ref{fig4-streak}) and correlated to parameters measured with beam quality monitors in the SPS \cite{sz_at_SPS}. Measurements performed in the SPS during data taking were then used to calculate the bunch length on a shot-by-shot basis.

The initial beam radius $\sigma_r$ and emittance $\epsilon_b$ were not measured with the required precision during the discussed parameter scans, so their input values are fitted to the experimental data in four steps. First, we determine the radial beam flux distributions by separate averaging of left and right halves of the images (Fig.\,\ref{fig3-profile}). This allows evaluation of the symmetry of the distribution and verifies that the axisymmetric model is applicable. The systematic asymmetry that we see in the halo is relatively low and can be caused by transverse misalignments between the proton beam and the plasma column. Second, we average the measured distributions over 30 shots. Third, in simulations, we vary the beam angular spread $\delta \alpha \approx \epsilon_b / (\gamma \sigma_r)$, where $\gamma$ is the relativistic factor of the beam. It determines the width of the beam head that propagates in the neutral gas and therefore the shape of the core part of beam flux distribution (compare red solid and magenta dotted lines in Fig.\,\ref{fig3-profile}). Fourth, we adjust the beam radius to match the height and width of the flux profile ``shoulders'' (compare red solid and blue dashed lines in Fig.\,\ref{fig3-profile}). The height and width of the shoulders depend on the strength of the fields produced in the plasma and are therefore sensitive to the initial beam density, which scales as $\sigma_r^{-2}$. Even with the best parameter fit (red line), the simulated beam profile is sharply peaked at the axis and therefore locally differs from the measured profile. This happens because both radial and angular distributions of beam micro-bunches in the plasma are strongly peaked in the perfectly axisymmetric case \cite{PoP24-023119}, while in the experiment this narrow spike is blurred. The obtained values of $\sigma_r$ and $\epsilon_b$ depend weakly on beam population. We approximate them as
\begin{equation}\label{em}
    \epsilon_b [\mu\mathrm{m}] = 2.3 + 4.5 \times 10^{-12} N_b, \ \ \sigma_r [\mu\mathrm{m}] = 105 \sqrt{\epsilon_b [\mu\mathrm{m}]},
\end{equation}
so the beam parameters vary from $\epsilon_b = 2.75$\,mm\,mrad and $\sigma_r = 175\,\mu$m at $N_b = 10^{11}$ to the values listed in Table~\ref{t1} at higher $N_b$.

The simulated maximum radius of deflected protons $r_\text{max}$ depends on the number of beam macro-particles used. The larger the number, the better we resolve the ``wings'' of the Gaussian angular distribution, the larger the maximum initial transverse momentum beam macro-particles can have, and the larger the maximum deflection we observe. To avoid this ambiguity, we consider $r_\text{max}$ as the radius at which the simulated beam flux density equals the noise level in experiments.

\begin{figure}[tb]
\includegraphics{fig5-wall.pdf}
\caption{(a) Density of plasma electrons (color map) and trajectories of individual electrons (lines) after 4 meters of beam propagation in the plasma. The vertical black line shows the cross-section detailed in Fig.\,\ref{fig6-expl}. (b) Wakefield potential $\Phi$ on the axis after the same propagation distance calculated with wide ($30c/\omega_p \approx 11.4\,\text{mm}$) and narrow ($1.7\,\text{mm} \approx 4.5c/\omega_p$) simulation windows.}\label{fig5-wall}
\end{figure}

\begin{figure}[tb]
\includegraphics{fig6-expl.pdf}
\caption{Radial dependencies of plasma-period-averaged ion density $\langle n_i \rangle$, electron density $\langle n_e \rangle$, and radial electric field $\langle E_r \rangle$ 8\,cm behind the seed laser pulse after 4 meters of beam propagation in the plasma. The electron density outside the plasma column is 100 times magnified for visibility. Arrows schematically show the field effect on plasma electrons and beam protons.}\label{fig6-expl}
\end{figure}

\begin{figure}[tb]
\includegraphics{fig7-charge.pdf}
\caption{Beam radius on the screens for low (a) and high (b) density plasmas and varying beam population. Top and bottom groups of lines in each panel are for two imaging stations. Black dots with error bars (``expt'') are experimentally measured values (individual shots), lines are simulation results. The error bars were determined the same way as in \cite{NIMA-909-123} and show the standard deviation of the determined radii along different directions of the image. Solid blue lines are for optimum parameter sets. The green dashed line in (a) shows the effect of a narrow simulation window ($R_w = 1.7$\,mm). The additional lines in (b) are for narrower beam (with $\sigma_r$ reduced by 20\% at the same $\delta \alpha$, green dashed line) or lower-emittance beam ($\delta \alpha$ reduced by 20\% at the same $\sigma_r$, orange dotted line). Dotted red lines are scalings $\propto N_b^{1/3}$ in (a) and $\propto N_b^{1/2}$ in (b).}\label{fig7-charge}
\end{figure}

The wide simulation window ($30c/\omega_p$) is necessary for correct calculation of high-energy electron trajectories. These electrons appear as a result of wavebreaking, which often occurs in the simulated regimes after beam micro-bunching [Fig.\,\ref{fig5-wall}(a)]. The electrons are ejected from the plasma column radially and leave an unbalanced positive charge behind \cite{PRL112-194801, PoP25-063108}. The charge separation electric field keeps them from traveling far, so these electrons form a negatively charged halo around the positively charged plasma column (Fig.\,\ref{fig6-expl}). The halo electrons move predominantly in the positive $z$ direction \cite{PRL112-194801} and efficiently interact with the wakefield when crossing the near-axis area: they can take energy from some regions or release it in others. As a consequence, the wakefield damps when a substantial number of halo electrons returns to the near-axis area. This is seen when comparing trajectories of halo electrons [Fig.\,\ref{fig5-wall}(a)] and the wakefield potential on the axis [Fig.\,\ref{fig5-wall}(b)]. The wakefield potential $\Phi$ characterizes the force exerted on axially moving ultra-relativistic particles:
\begin{equation}\label{e1}
    E_z = - \frac{\partial \Phi}{\partial z}, \qquad E_r - B_\varphi = - \frac{\partial \Phi}{\partial r}.
\end{equation}
A reduced amplitude of potential oscillations on the axis evidences reduction of both accelerating and defocusing properties of the wakefield.

\section{Comparison of simulations and experiments}

We compare simulations and experiments in three parameter scans. The first two are beam population scans at low ($1.965\times 10^{14}\,\text{cm}^{-3}$) and high ($7.7\times 10^{14}\,\text{cm}^{-3}$) plasma densities (Fig.\,\ref{fig7-charge}). The best agreement variants (blue lines) correspond to the approximation (\ref{em}) and other parameters from Table~\ref{t1}, with the exception that the high density scan was with $\sigma_z = 7$\,cm and $\xi_s = 1.9$\,cm.

The green dashed line in Fig.\,\ref{fig7-charge}(a) shows the effect of a narrow simulation window. Disagreement with experiments is expected, but the sign of the effect may seem counter-intuitive. If the simulation wall is located close to the plasma boundary, the positive electric field around the plasma (Fig.\,\ref{fig6-expl}) does not additionally accelerate beam protons in the radial direction. However, the narrow window also ``switches off'' the effect of halo electrons, the wakefield becomes stronger, and we observe larger proton deflection in simulations.

The additional lines in Fig.\,\ref{fig7-charge}(b) illustrate the sensitivity to several beam parameters.  The simulations agree well with the measured data over most of the range of proton beam population for variations of 20~\% in the beam radius and the beam angular divergence. The largest discrepancy comes from varying the beam radius at small bunch populations.  In this case, a 20\% reduction in the beam radius results in visible disagreement with measurements (green dashed line) for the second imaging station. The sensitivity to the angular spread is lower: a 20\% reduction changes $r_{\rm max}$ almost imperceptibly (orange dotted line). Variation of other beam parameters also have small effects on proton defocusing, except the seed position discussed in Ref.\,\cite{Marlene}.

Theory \cite{PoP20-083119} suggests that the maximum wakefield at AWAKE may be limited by nonlinear elongation of the plasma wave period. The wave goes out of resonance with the bunch train, and electric field growth saturates at about $0.4 mc\omega_p/e$. If this mechanism is indeed the main limiting factor, then the maximum wakefield amplitude should scale as the contribution of an individual bunch to the power of 1/3, or as $N_b^{1/3}$ (see Eq.\,(12) of \cite{PoP20-083119}). The dependence of maximum proton deflection on beam population at the low plasma density follows this scaling remarkably well [dotted line in Fig.\,\ref{fig7-charge}(a)]. At the high density, however, the maximum radius scales as $N_b^{1/2}$ [dotted line in Fig.\,\ref{fig7-charge}(b)], which indicates that the microbunch train is not dense or long enough to excite the nonlinear wakefields at this plasma density. In this case, the radial momentum gained by the deflected protons $p_\text{max}$ is defined by the depth of the wakefield potential well $e\Phi_\text{max}$ that scales as $N_b$. As a result, $p_\text{max} \propto \sqrt{2 e \Phi_\text{max} W_b}/c \propto \sqrt{N_b}$.

\begin{figure}[htb]
\includegraphics{fig8-grad.pdf}
\caption{Beam radius on the screens for different plasma density gradients (characterized by the relative density increase at the downstream plasma end): black dots are the experimentally measured values; error bars show the standard deviation of all individual measurements (about 30 shots) at a given gradient setting, lines are simulation results for different beam radii.}\label{fig8-grad}
\end{figure}

For the density gradient scan, the simulations also follow the measured trend (Fig.\,\ref{fig8-grad}). Here the plasma density changes linearly along the cell due to a temperature difference between the two Rb reservoirs attached to the cell at the ends \cite{JPD51-025203}. The density at the upstream end is always $1.8\times 10^{14}\,\text{cm}^{-3}$. Beam parameters are those from Table~\ref{t1}, except $\sigma_z = 7.4$\,cm and $\xi_s = 2.4$\,cm.
The simulations show significant sensitivity to the beam radius for positive gradients larger than $+5$\%.  Overall, the best agreement is for $\sigma_r \approx 220\,\mu$m.

\begin{figure}[htb]
\includegraphics{fig9-gmaxf.pdf}
\caption{Maximum amplitude of the wakefield potential $\Phi_{\rm max}$ versus beam propagation distance $z$ for $\sigma_r = 220\,\mu$m and different magnitudes of relative density increase along the plasma cell.}\label{fig9-gmaxf}
\end{figure}
Small positive gradients contribute to beam self-modulation by increasing the total charge of microbunches \cite{NIMA-829-63}. The effect is similar to that of the density step \cite{PoP18-024501,PoP22-103110}: larger fractions of initially formed micro-bunches remain focused by the wave because of plasma wavelength shortening. With small positive gradients (density increase about 2\% over 10 meters), the longitudinal wakefield is stronger \cite{NIMA-829-63}, and externally injected electrons gain a higher energy \cite{Nat.561-363}. Higher positive gradients, however, are destructive and reduce the maximum wakefield (Fig.\,\ref{fig9-gmaxf}), which results in weaker deflection of defocused protons (Fig.\,\ref{fig8-grad}). High negative gradients also reduce the wakefield, but not the maximum deflection radius. This happens because the strongest defocused protons do not always experience the strongest wakefield generated by the beam. Negative density gradients change the wavelength so that the micro-bunches mostly fall into the defocusing wave phase, and more efficient coupling with the transverse field compensates the lower wakefield amplitude.

\section{Summary}

To conclude, simulations with two-dimensional (axisymmetric) particle-in-cell code LCODE agree with the measurements of the maximum radius of deflected protons in three different experimental scans: proton beam population scans ($1 \div 3 \times 10^{11}$ particles) for high and low plasma densities, and gradient scan ($-20 \div 20\%$). The difference between calculated and measured values is about or below 5\%. This means that particle-in-cell codes can be used to interpret the physics underlying experimental data on SSM. However, the agreement with experiments can be achieved only for wide enough simulation windows ($\sim 30c/\omega_p$ for AWAKE parameters) that are necessary to correctly account for ejected plasma electrons that charge the plasma column positively and then return to the axis destroying the wakefield. This additionally increases the computational power required for simulations of long microbunch trains in a plasma.

\ack

This work was supported in parts by the Foundation for the Development of Theoretical Physics and Mathematics ``BASIS''; a Leverhulme Trust Research Project Grant RPG-2017-143 and by STFC (AWAKE-UK, Cockcroft Institute core, John Adams Institute core, and UCL consolidated grants), United Kingdom; a Deutsche Forschungsgemeinschaft project grant PU 213-6/1 ``Three-dimensional quasi-static simulations of beam self-modulation for plasma wakefield acceleration''; the National Research Foundation of Korea (Nos.\ NRF-2016R1A5A1013277 and NRF-2019R1F1A1062377); the Portuguese FCT---Foundation for Science and Technology, through grants CERN/FIS-TEC/0032/2017, PTDC-FIS-PLA-2940-2014, UID/FIS/50010/2013 and SFRH/IF/01635/2015; NSERC and CNRC for TRIUMF's contribution; the U.S.\ National Science Foundation under grant PHY-1903316; the Wolfgang Gentner Programme of the German Federal Ministry of Education and Research (grant no.\ 05E15CHA); and the Research Council of Norway. A.A.\,Gorn and K.V.\,Lotov acknowledge the support of the Russian Foundation for Basic Research (grant 19-32-90125). M. Wing acknowledges the support of the Alexander von Humboldt Stiftung and DESY, Hamburg. Support of the Wigner Datacenter Cloud facility through the
``Awakelaser'' project and the support of P\'eter L\'evai is acknowledged. The work of V. Hafych has been supported by the European Union's Framework Programme for Research and Innovation Horizon 2020 (2014--2020) under the Marie Sklodowska-Curie Grant Agreement No.\ 765710.  The AWAKE collaboration acknowledge the SPS team for their excellent proton delivery. LCODE simulations were performed on HPC-cluster ``Akademik V.M. Matrosov'' \cite{matrosov}.

\section*{References}

redpajama_set_name: RedPajamaArXiv