# Coherent inverse Compton scattering by bunches in fast radio bursts

Bing Zhang  Department of Physics and Astronomy, University of Nevada, Las
Vegas, Las Vegas, NV 89154<EMAIL_ADDRESS>
###### Abstract

The extremely high brightness temperature of fast radio bursts (FRBs) requires
that their emission mechanism must be “coherent”, either through concerted
particle emission by bunches or through an exponential growth of a plasma wave
mode or radiation amplitude via certain maser mechanisms. The bunching
mechanism has been mostly discussed within the context of curvature radiation
or cyclotron/synchrotron radiation. Here we propose a family of model invoking
coherent inverse Compton scattering (ICS) of bunched particles that may
operate within or just outside of the magnetosphere of a flaring magnetar.
Crustal oscillations during the flaring event may excite low-frequency
electromagnetic waves near the magnetar surface. The X-mode of these waves
could penetrate through the magnetosphere. Bunched relativistic particles in
the charge starved region inside the magnetosphere or in the current sheet
outside of the magnetosphere would upscatter these low-frequency waves to
produce GHz emission to power FRBs. The ICS mechanism has a much larger
emission power for individual electrons than curvature radiation. This greatly
reduces the required degree of coherence in bunches, alleviating several
criticisms to the bunching mechanism raised in the context of curvature
radiation. The emission is $\sim 100\%$ linearly polarized (with the
possibility of developing circular polarization) with a constant or varying
polarization angle across each burst. The mechanism can account for a narrow-
band spectrum and a frequency downdrifting pattern, as commonly observed in
repeating FRBs.

fast radio bursts – radio transient sources – magnetars

=1

##  1 Introduction

Fast radio bursts (FRBs) (Lorimer et al., 2007; Petroff et al., 2019; Cordes &
Chatterjee, 2019; Zhang, 2020) have a brightness temperature

$\displaystyle T_{b}$ $\displaystyle\simeq$ $\displaystyle\frac{{\cal
S}_{\nu,p}D_{\rm A}^{2}}{2\pi k_{\rm B}(\nu\Delta t)^{2}}=(1.2\times 10^{36}\
{\rm K})$ (1) $\displaystyle\times$ $\displaystyle D_{\rm A,28}^{2}({\cal
S}_{\nu,p}/{\rm Jy})\nu_{9}^{-2}\Delta t_{-3}^{-2},$

where $k_{\rm B}$ is Boltzmann constant, $S_{\nu,p}$ is specific flux at the
peak time, $\nu$ is observing frequency, $\Delta t$ is variability timescale,
$D_{\rm A}$ is angular distance of the source, and the convention
$Q_{n}=Q/10^{n}$ has been adopted in cgs units throughout the paper. This is
much greater than the maximum brightness temperature for an incoherent
emitting source

$T_{\rm b,max}^{\rm incoh}\simeq\Gamma\gamma m_{e}c^{2}/k_{\rm B}=(5.9\times
10^{13}\ {\rm K})\Gamma_{2}\gamma_{2},$ (2)

where $m_{e}$, $c$, $k$ are the fundamental constants electron mass, speed of
light and Boltzman constant, respectively, $\Gamma$ is the bulk Lorentz factor
of the emitter towards earth (taken as unity if the source is not moving with
a relativistic speed), and $\gamma={\rm max}(\gamma_{m},\gamma_{a})$ is the
characteristic electron Lorentz factor in the emission region, which, for a
synchrotron source, is the greater of the minimum injection Lorentz factor
$\gamma_{m}$ and the Lorentz factor $\gamma_{a}$ corresponding to synchrotron
self-absorption (Kumar & Zhang, 2015). This suggests that FRB emission
mechanism must be coherent.

Many coherent radiation mechanisms have been explored to interpret the
emission of radio pulsars, which also have the brightness temperature
$T_{b}\gg T_{\rm b,max}^{\rm incoh}$ (but is about 10 orders of magnitude
lower than that of FRBs). In general, these mechanisms can be grouped into
three types (e.g. Melrose, 1978): coherent radiation by bunches, intrinsic
growth of plasma wave modes, and maser mechanism. The second and third
mechanisms were also termed as “plasma maser” and “vacuum maser” mechanisms,
respectively. Some of these mechanisms have been reinvented to interpret FRB
coherent emission (e.g. Lu & Kumar, 2018; Zhang, 2020; Xiao et al., 2021;
Lyubarsky, 2021, for surveys of various mechanisms discussed in the
literature).

This paper mainly concerns about the first type of coherent mechanism, namely,
coherent radiation by bunches, also called the “antenna” mechanism. Within
this mechanism, charged particles are clustered in both position and momentum
spaces and emit as a macroscopic charge. The emission power of the bunch is
$P_{b}\simeq N_{e,b}^{2}P_{e}$, where $P_{e}$ is the power of individual
electrons, and $N_{e,b}$ is the total number of net charges in the bunch.
Within the FRB context, a widely discussed possibility is coherent curvature
radiation by bunches (Katz, 2014; Kumar et al., 2017; Yang & Zhang, 2018; Lu &
Kumar, 2018; Wang et al., 2019; Kumar & Bošnjak, 2020; Lu et al., 2020; Yang
et al., 2020; Wang & Lai, 2020; Cooper & Wijers, 2021). Another related
mechanism is the so-called “synchrotron maser” mechanism in 90o magnetized
relativistic shocks, which in fact invoke cyclotron-radiating bunches of
charged particles in momentum space (Lyubarsky, 2014; Beloborodov, 2017, 2020;
Metzger et al., 2019; Plotnikov & Sironi, 2019; Margalit et al., 2020). These
mechanisms can interpret some FRB emission properties but also suffer from
criticisms in theoretical and/or observational aspects (e.g. Zhang, 2020;
Lyubarsky, 2021).

In this paper, we discuss another family of bunching coherent mechanisms,
namely, coherent inverse Compton scattering (ICS) by bunches. Before getting
into the details of the model, it is informative to clarify the meaning of
“bunches” discussed in this paper. In the pulsar literature, “bunches” are
usually treated as giant charges surrounded by a background plasma (e.g. Gil
et al., 2004). The mechanism was criticized for the the formation and
maintenance mechanisms of such bunches (e.g. Melrose, 1978). Also the
surrounding plasma works against the coherence of the bunch and suppresses the
coherent emission (Gil et al., 2004; Lyubarsky, 2021). For the curvature
radiation mechanism widely discussed in the literature, since the radiation
power of individual particles is low, highly coherent bunches are required and
these criticisms may be relevant (but see Melikidze et al., 2000; Qu & Zhang,
2021). However, as will be shown below, in the case of ICS emission, since the
emission power of individual particles is much higher, the required coherence
for bunches is much reduced. Rather than distinct giant charges surrounded by
an ambient medium, the bunches discussed in this paper are merely charge
density fluctuations in an relativistic particle outflow (see also Yang &
Zhang 2018).

Figure 1: The geometric sketch of the proposed model. Charges in the
magnetosphere near the crustal oscillation region send out low-frequency
electromagnetic waves which propagate across the closed field line region to
reach the emission sites, either in the outer magnetosphere in the open field
line region where an $E_{\parallel}$ is developed, or in the current sheet
region outside the magnetosphere. Relativistic particles in the emission
region upscatter the low-frequency waves and produce FRBs. The FRB emission
direction is along the local magnetic field lines in the first scenario and is
out of page in the second scenario.

##  2 Coherent ICS by bunches

###  2.1 General picture

The general picture of the hereby discussed mechanism is the following (Figure
1): Sudden cracking of the curst of a magnetar excites crustal quakes and
plasma oscillations near the surface of the magnetar. Low-frequency
electromagnetic waves with angular frequency $\omega_{0}$ would be generated
by the coherently oscillating charged particles and would propagate across the
magnetosphere. The waves encounter bunches of relativistic particles with a
typical Lorentz factor $\gamma$ and a net charge $N_{e,b}e$, which moves along
magnetic field lines in the outer magnetosphere or in current sheets outside
of the magnetosphere. In the rest-frame of the particle bunch, the incident
wave frequency is boosted to $\sim\gamma(1-\beta\cos\theta_{i})\nu_{0}$. The
bunched particles are induced to oscillate at the same frequency in the
comoving frame, which is transferred to a lab-frame outgoing ICS frequency

$\nu\simeq\gamma^{2}\nu_{0}(1-\beta\cos\theta_{i})=(1\ {\rm GHz})\
\gamma_{2.5}^{2}\nu_{0,4}(1-\beta\cos\theta_{i}),$ (3)

where $\nu_{0}=\omega_{0}/2\pi$ and $\theta_{i}$ is the angle between the
incident photon momentum and the electron momentum. Similar to curvature
radiation, such an ICS mechanism does not depend intrinsically on the
dispersive properties of the plasma and may be treated assuming vacuum wave
properties. Different from the curvature radiation mechanism which appeals to
the curved trajectory of charged particles for acceleration, the ICS mechanism
invokes an oscillating electromagnetic field from the low frequency waves to
accelerate bunched particles.

Such a mechanism is not new in the pulsar literature. Linear acceleration
emission (LAE) invokes an oscillating $E_{\parallel}$ along the direction of
particle motion (e.g. Melrose, 1978; Rowe, 1995). A free-electron laser (FEL)
invokes interaction between an electromagnetic disturbance (wiggler) with a
relativistically moving bunch (e.g. Fung & Kuijpers, 2004) to power pulsar
radio emission or FRB emission Lyutikov (2021a). Lyubarskii (1996) introduced
induced upscattering of longitudinal oscillations to interpret pulsar radio
emission. Qiao & Lin (1998) considered low-frequency electromagnetic waves
generated from regular discharges of an inner vacuum gap in the pulsar polar
cap region and showed that ICS off these waves may reproduce many
observational features of pulsar radio emission (Qiao et al., 2001; Xu et al.,
2000).

###  2.2 Low frequency waves: generation and propagation

Since FRBs are rarely produced from their sources (e.g. the Galactic magnetar
SGR 1935+2154 only produced one detected FRB since its discovery, CHIME/FRB
Collaboration et al. 2020; Bochenek et al. 2020, while many X-ray bursts have
been produced without association with any FRBs, Lin et al. 2020), it is
reasonable to assume that something extraordinary must have happened to
trigger an FRB. We envisage the low-frequency electromagnetic waves generated
near the surface of the neutron star as one condition to produce FRBs. This
may be related to violent oscillations of charged particles near the surface
triggered by, e.g. strong oscillations of the neutron star crust during a
starquake (e.g. Thompson & Duncan, 2001; Beloborodov & Thompson, 2007; Wang et
al., 2018; Dehman et al., 2020; Yang & Zhang, 2021). The bulk of the global
oscillation energy would be carried by Alfvén waves that propagate along the
magnetic field lines in the form a magnetic pulse, which would deposit energy
in the form of particle energy and radiation at a large distance, as envisaged
in most of the FRB models (e.g. Kumar & Bošnjak, 2020; Lyubarsky, 2020; Yuan
et al., 2020). However, the same shearing oscillations of the crust would
induce oscillations of charges in the near-surface magnetosphere in the
direction perpendicular to magnetic field lines, so that a small fraction of
the oscillation energy may be radiated by these charges coherently and nearly
isotropically as X-mode electromagnetic waves with angular frequency
$\omega_{0}$, similar to an antenna in a radio station on Earth. Strictly
speaking, the electromagnetic waves discussed here are electromagnetic modes
propagating in a plasma, which may be equivalent to some forms of fast
magnetosonic waves. However, since the magnetization factor $\sigma$ is
extremely high in the inner magnetosphere, the magnetosonic waves all
propagate at essentially the speed of light and would behave like
electromagnetic waves in vacuum with a proper dispersion relation. Since the
electromagnetic waves do not carry the bulk of FRB energy but rather provide
seed photons for particles to upscatter at a large radius, the brightness
temperature of the low-frequency coherent emission does not need to be very
high. In this model, the low frequency waves should at least last several
milliseconds (the duration of the FRB), but could last longer (in which case
the FRB duration would be defined by the duration of the magnetic pulse
dissipation leading to coherent ICS emission by relativistic particles).

For a cold, non-relativistic plasma, the X-mode electromagnetic waves can
propagate in the frequency regimes $\omega>\omega_{\rm R}$, $\omega_{\rm
L}<\omega<\omega_{\rm uh}$, and $\omega<\omega_{\rm lh}$, where $\omega_{\rm
R}$, $\omega_{\rm L}$ are the cutoff frequencies for the R-mode and L-mode and
$\omega_{\rm uh}$ and $\omega_{\rm lh}$ are the upper and lower hybrid
principle resonant frequencies, respectively (Boyd & Sanderson, 2003). For an
electron-positron pair plasma relevant for a pulsar/magnetar magnetosphere,
one has $\omega_{\rm R}=\omega_{\rm uh}=\omega_{\rm
L}=\sqrt{\omega_{p}^{2}+\omega_{\rm B}^{2}}$ and $\omega_{\rm lh}=\omega_{\rm
B}$, so that the X-mode electromagnetic waves are free to propagate at
$\omega>\sqrt{\omega_{p}^{2}+\omega_{\rm B}^{2}}$ or $\omega<\omega_{\rm B}$
(see Appendix). For a magnetar with surface magnetic field $B_{s}=10^{15}{\rm
G}B_{s,15}$ and rotation period $P$, the Goldreich-Julian charge number
density (Goldreich & Julian, 1969; Ruderman & Sutherland, 1975) is

$n_{\rm GJ}\simeq\frac{\Omega B}{2\pi ce}\simeq(6.9\times 10^{7}\ {\rm
cm^{-3}})B_{s,15}P^{-1}\hat{r}_{2}^{-3},$ (4)

where $\Omega=2\pi/P$ is the angular frequency, and

$\hat{r}\equiv\frac{r}{R}$ (5)

is the radius $r$ normalized to the neutron star radius $R=10^{6}{\rm
cm}R_{6}$. The plasma number density can be estimated as

$n=\xi n_{\rm GJ},$ (6)

where $\xi$ is the multiplicity parameter due to electron-positron pair
production. As a result, the plasma frequency and the Larmor frequency can be
estimated as

$\displaystyle\omega_{p}$ $\displaystyle=$ $\displaystyle\left(\frac{4\pi
ne^{2}}{m_{e}}\right)^{1/2}$ (7) $\displaystyle\simeq$
$\displaystyle(4.7\times 10^{8}\ {\rm rad\ s^{-1}})\
\xi^{1/2}B_{s,15}^{1/2}P^{-1/2}\hat{r}_{2}^{-3/2}$ $\omega_{\rm
B}=\frac{eB}{m_{e}c}\simeq(1.8\times 10^{16}\ {\rm rad\ s^{-1}})\
B_{s,15}\hat{r}_{2}^{-3}.$ (8)

It is also informative to evaluate these frequencies at the light cylinder of
the magnetar, which has $r=R_{\rm lc}=c/\Omega$ and

$\hat{r}_{\rm lc}=(4.8\times 10^{3})PR_{6}^{-1}.$ (9)

This gives

$\displaystyle\omega_{p}^{\rm lc}$ $\displaystyle\simeq$
$\displaystyle(1.4\times 10^{6}\ {\rm rad\ s^{-1}})\
\xi^{1/2}B_{s,15}^{1/2}P^{-2}R_{6}^{3/2}$ (10) $\displaystyle\omega_{\rm
B}^{\rm lc}$ $\displaystyle\simeq$ $\displaystyle(1.6\times 10^{11}\ {\rm rad\
s^{-1}})\ B_{s,15}P^{-3}R_{6}^{3}.$ (11)

One can see that for the low-frequency EM waves we are interested in, the
condition $\omega_{0}\ll\omega_{p}\ll\omega_{\rm B}$ is satisfied throughout
the magnetosphere. Therefore the X-mode low frequency waves are transparent to
the magnetosphere in all directions for essentially all frequencies (see also
Lu et al., 2019).

The above treatment makes the assumption of a cold plasma. Near surface
electromagnetic waves need to penetrate through the closed field line region
to reach the emission sites in our two scenarios (Fig.1). Since the closed
field line region of a magnetar is likely populated with inactive non-
relativistic particles, the cold plasma dispersion relation treatment can
serve the purpose for our discussion. We note that the dispersion relation is
more complicated when relativistic plasmas are considered (e.g. Arons &
Barnard, 1986; Lyubarskii & Petrova, 1998; Melrose, 2017). Nonetheless, the
conclusion that X-mode low-frequency electromagnetic waves are transparent in
a magnetar magnetosphere remains valid for such more complicated situations.

###  2.3 Cross section and ICS power

In the following, we perform a vacuum treatment of the ICS process for
simplicity. A more rigorous treatment should consider various plasma effects
(e.g. Melrose, 1978; Rowe, 1995; Lyubarskii, 1996; Fung & Kuijpers, 2004;
Lyutikov, 2021a). However, a vacuum treatment can catch the essential features
of this family of models, as discussed below.

The Compton scattering (in the rest frame of electron) cross section in a
strong magnetic field is significantly modified from the Thomson cross section
$\sigma_{\rm T}$. The cross section for the two modes of photons read (Herold,
1979; Xia et al., 1985)

$\displaystyle\sigma^{\prime}(1)$ $\displaystyle=$ $\displaystyle\sigma_{\rm
T}\left\\{\sin^{2}\theta^{\prime}_{i}+\frac{1}{2}\cos^{2}\theta^{\prime}_{i}\left[\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}+\omega_{\rm
B})^{2}}+\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}-\omega_{\rm
B})^{2}}\right]\right\\},$ (12) $\displaystyle\sigma^{\prime}(2)$
$\displaystyle=$ $\displaystyle\frac{\sigma_{\rm
T}}{2}\left[\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}+\omega_{\rm
B})^{2}}+\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}-\omega_{\rm
B})^{2}}\right].$ (13)

where 1 denotes the mode that the electric vector ${\bf E}^{\prime}$ of the
incident wave is parallel to the $({\bf k},{\bf B})$ plane, and 2 denotes the
mode that ${\bf E}^{\prime}$ is perpendicular to the plane. Here the primed
symbols denote in the rest frame of the electron (since electrons move along
magnetic field lines, $\omega_{\rm B}$ does not change when the frame
changes). In this frame, the photon incident angle $\theta^{\prime}_{i}$ is
connected to the lab-frame incident angle $\theta_{i}$ through

$\displaystyle\sin\theta^{\prime}_{i}$ $\displaystyle=$
$\displaystyle\frac{\sin\theta_{i}}{\gamma(1-\beta\cos\theta_{i})},$ (14)
$\displaystyle\cos\theta^{\prime}_{i}$ $\displaystyle=$
$\displaystyle\frac{\cos\theta_{i}-\beta}{1-\beta\cos\theta_{i}},$ (15)

where $\gamma$ and $\beta$ are the Lorentz factor and dimensionless speed of
the electron, respectively. In the lab frame, the ICS cross section is related
to the the rest-frame Compton scattering cross section through
$\sigma=(1-\beta\cos\theta_{i})\sigma^{\prime}$ (Pacholczyk, 1970). This
finally gives

$\displaystyle\sigma(1)$ $\displaystyle=$ $\displaystyle\sigma_{\rm
T}\left\\{\frac{\sin^{2}\theta_{i}}{\gamma^{2}(1-\beta\cos\theta_{i})}+\frac{(\cos\theta_{i}-\beta_{i})^{2}}{2(1-\beta\cos\theta_{i})}\right.$
(16) $\displaystyle\times$
$\displaystyle\left.\left[\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}+\omega_{\rm
B})^{2}}+\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}-\omega_{\rm
B})^{2}}\right]\right\\},$ $\displaystyle\sigma(2)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{\rm
T}}{2}\left[\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}+\omega_{\rm
B})^{2}}+\frac{{\omega^{\prime}}^{2}}{(\omega^{\prime}-\omega_{\rm
B})^{2}}\right].$ (17)

At this point, it is informative to compare the order of magnitude of the two
terms in the curly brackets. In the lab frame, $\theta_{i}$ can be an
arbitrary angle, so that all the factors involving $\theta_{i}$ could be of
the order of unity. For our nominal parameters $\gamma\sim 10^{2.5}$ and
$\omega_{0}\sim(2\pi)10^{4}$ and noticing
$\omega^{\prime}\simeq\gamma\omega_{0}(1-\beta\cos\theta_{i})$, the first term
in $\sigma(1)$ is therefore $\propto\gamma^{-2}\sim 10^{-5}$ and the second
term (also $\sigma(2)$) is $\propto(\omega^{\prime}/\omega_{\rm
B})^{2}\sim(10^{6}/10^{11})^{2}\sim 10^{-10}$ even at the light cylinder where
$\omega_{\rm B}$ is the lowest in the magnetosphere. One can therefore ignore
$\sigma(2)$ and the second term of $\sigma(1)$, so that only (see also Qiao &
Lin, 1998)

$\sigma(1)\simeq\sigma_{\rm
T}\frac{\sin^{2}\theta_{i}}{\gamma^{2}(1-\beta\cos\theta_{i})}\simeq(6.65\times
10^{-30}\ {\rm cm^{2}})\gamma_{2.5}^{-2}f(\theta_{i})$ (18)

is relevant, where $f(\theta_{i})=\sin^{2}\theta_{i}/(1-\beta\cos\theta_{i})$
is defined.

The ICS emission power of a single relativistic electron may be estimated as

$\displaystyle P_{e}^{\rm ICS}$ $\displaystyle\simeq$
$\displaystyle\frac{4}{3}\gamma^{2}\sigma(1)cU_{\rm ph}$ (19)
$\displaystyle\simeq$ $\displaystyle(2.1\times 10^{-7}\ {\rm erg\
s^{-1}})f(\theta_{i})(\delta B_{0,6})^{2}\hat{r}_{2}^{-2},$

where

$\displaystyle U_{\rm ph}$ $\displaystyle\simeq$
$\displaystyle\frac{(|{\bf\delta E_{0}\times\delta
B_{0}}|)^{2}}{4\pi}\hat{r}^{-2}$ (20) $\displaystyle\simeq$
$\displaystyle(8.0\times 10^{6}\ {\rm erg\ cm^{-3}})(\delta
B_{0,6})^{2}\hat{r}_{2}^{-2}$

is the photon energy density in the emission region, and ${\bf\delta E_{0}}$
and ${\bf\delta B_{0}}$ are the electric and magnetic vectors of the low
frequency waves near the surface region, which we normalize to a relatively
small value $\delta B_{0}\sim 10^{6}\ {\rm G}$. Comparing with the curvature
radiation (CR) emission power of a single electron

$P_{e}^{\rm CR}=\gamma^{4}\frac{2e^{2}c}{3\rho^{2}}\simeq(4.6\times 10^{-15}\
{\rm erg\ s^{-1}})\gamma_{2.5}^{4}\rho_{8}^{-2},$ (21)

one can see that the ICS power is much larger, i.e. $P_{e}^{\rm ICS}\gg
P_{e}^{\rm CR}$, a known fact in pulsar studies (Qiao & Lin, 1998; Zhang et
al., 1999). This suggests that when coherent ICS operates, the effect of
coherent CR is negligibly small.

###  2.4 Coherent ICS emission by bunches

A bunch of leptons with a net charge value of $N_{e,b}e$ can radiate
coherently in response to the low frequency waves, with an emitted power of
$\sim N_{e,b}^{2}P_{e}^{\rm ICS}$. The observed power is boosted by a factor
of $\gamma^{2}$ because the observer time is shorter by a factor of
$(1-\beta\cos\theta)\sim 1/\gamma^{2}$ (when the angle between the charge
motion direction and the line of sight is $\theta<1/\gamma$) with respect to
the emission time. For a total number of bunch number $N_{b}$, the total
luminosity of the emitter may be written as

$L\simeq N_{b}N_{e,b}^{2}P_{e}^{\rm ICS}\gamma^{2}.$ (22)

One may estimate

$N_{e,b}=\zeta n_{\rm GJ}A_{b}\lambda,$ (23)

where $A_{b}$ is the cross section of the bunch, $\lambda=c/\nu$ is the
wavelength, and $\zeta$ is the factor to denote the net charge density with
respect to $n_{\rm GJ}$. The most conservative estimate is (Kumar et al.,
2017)

$A_{b}^{\rm min}=\pi(\gamma\lambda)^{2},$ (24)

since it describes the causally connected area of a relativistic bunch111The
true $A_{b}$ can be in principle much larger, up to the area defined by the
Fresnel radius i.e. $A_{b}\simeq\pi x\lambda$, where $x$ is the distance
between the emission site and the projected focus of the tangential lines of
the emission region field lines.. Since the ICS mechanism is very efficient,
in the following we use this conservative limit. Plugging Eqs. (4), (19), (23)
, and (24) in Eq. (22), one finally gets

$\displaystyle L$ $\displaystyle\simeq$ $\displaystyle(7.3\times 10^{38}\ {\rm
erg\ s^{-1}})\zeta^{2}N_{b,5}\gamma_{2.5}^{6}\nu_{9}^{-6}$ (25)
$\displaystyle\times$ $\displaystyle B_{s,15}^{2}P^{-2}f(\theta_{i})\delta
B_{0,6}^{2}\hat{r}_{2}^{-8}.$

One may compare Eq.(25) with the true FRB luminosity from the observations,
i.e. $L^{\rm obs}=L_{\rm iso}^{\rm obs}{\rm
max}(\pi/\gamma^{2},\pi\theta_{j}^{2})$, where $\theta_{j}$ is the half
opening angle of the FRB jet. For the narrow jet scenario with solid angle
defined by the $\gamma^{-1}$ cone, one has $L^{\rm obs}\simeq(3\times 10^{38}\
{\rm erg\ s^{-1}})L^{\rm obs}_{\rm iso,43}\gamma_{2.5}^{-2}$. One can see that
the coherent ICS mechanism can easily account for the typical FRB luminosity
by only requiring a moderate number of $N_{b}\sim 10^{5}$ even with the most
conservative estimate of the coherent bunch cross section $A_{b}^{\rm min}$.
From Eq.(22), one can immediately see that the required $N_{b}N_{e,b}^{2}$ for
the coherent ICS model is smaller by a factor $P_{e}^{\rm CR}/P_{e}^{\rm
ICS}\sim 2.2\times 10^{-8}$ to interpret the same FRB luminosity for the same
electron Lorentz factor $\gamma$.

A bunch with $N_{e,b}$ particles cools more efficiently than single particles
since the emission power is $\sim N_{e,b}^{2}$ times of the individual
particle emission power. The bunch quickly loses the total energy, which
itself is not large enough to power FRB emission. In order to power FRB
emission, a sustained $E_{\parallel}$ is needed to continuously pump energy to
the bunch. Such a scenario is relevant to coherent CR by bunches (Kumar et
al., 2017). Below we investigate the case for coherent ICS emission. The
cooling timescale of the coherent ICS bunches may be estimated as

$\displaystyle t_{c,b}$ $\displaystyle=$ $\displaystyle\frac{N_{e,b}\gamma
m_{e}c^{2}}{N_{e,b}^{2}P_{e}^{\rm ICS}}\simeq(2.1\times 10^{-15}\ {\rm s})$
(26) $\displaystyle\times$
$\displaystyle\zeta^{-1}\nu_{9}^{3}\gamma_{2.5}^{-1}B_{s,15}^{-1}P[f(\theta_{i})]^{-1}\delta
B_{0,6}^{-2}\hat{r}_{2}^{5}.$

Since this is much shorter than the typical FRB duration, similar to the
bunched coherent CR mechanism, a parallel electric field is also needed to
continuously pump energy to the bunch. One may estimate the strength of this
$E_{\parallel}$ by requiring $N_{e,b}e(ct_{c,b})\simeq N_{e,b}\gamma
m_{e}c^{2}$, which gives

$\displaystyle E_{\parallel}$ $\displaystyle\simeq$
$\displaystyle\frac{N_{e,b}P_{e}^{\rm ICS}}{ec}\simeq(8.6\times 10^{9}\ {\rm
esu})$ (27) $\displaystyle\times$
$\displaystyle\zeta\nu_{9}^{-3}\gamma_{2.5}^{2}B_{s,15}P^{-1}f(\theta_{i})\delta
B_{0,6}^{2}\hat{r}_{2}^{-5}.$

This is larger than that required for the bunched coherent CR mechanism
because of a larger $P_{e}$. However, in view of the steep dependence on
$\hat{r}$, the required $E_{\parallel}$ comfortably falls into the range
available for magnetars if $\hat{r}$ is somewhat greater than 100. Several
mechanisms may provide this $E_{\parallel}$ at $\hat{r}\gg 1$. First, a slot
gap near the last open field line region may extend from the neutron surface
all the way to light cylinder because the region cannot be filled with
electron-positron pairs produced from a pair production cascade – a
consequence of straight line propagation of $\gamma$-rays and curved magnetic
field lines (e.g. Arons & Scharlemann, 1979; Muslimov & Harding, 2004). Charge
starvation may become more prominent in old, low-twist magnetars (Wadiasingh
et al., 2020). Second, Alfvén waves propagating along magnetic field lines
would enter a charge starved region at a large enough radius because the
charge density may become insufficient to provide the required current (e.g.
Kumar & Bošnjak 2020; Lu et al. 2020, cf. Chen et al. 2020). Third, in the
current sheet outside of the magnetosphere, an $E_{\parallel}$ naturally
exists in the magnetic reconnection layer (e.g. Lyubarsky, 2020;
Kalapotharakos et al., 2018; Philippov & Spitkovsky, 2018).

##  3 Salient features of the model

###  3.1 Two physical scenarios

After discussing some general properties of this model, we discuss two
specific scenarios within the framework of the magnetar repeating FRB models.
The geometric configurations of the two scenarios are presented in Figure 1.

The first scenario is similar to that of Kumar & Bošnjak (2020) and Lu et al.
(2020), but with coherent curvature radiation by bunches replaced by the more
efficient coherent ICS by bunches. Within this scenario, neutron star crust
cracking excites crustal seismic waves, which in turn excite Alfvén waves in
the magnetosphere. Low frequency electromagnetic waves (effectively fast
magnetosonic waves in a high-$\sigma$ medium) are generated and the X-mode
waves propagate with a speed close to the speed of light. The waves are
upscattered at a large enough altitude where $E_{\parallel}$ is developed,
either due to charge starvation in the Alfvén waves (Kumar & Bošnjak, 2020) or
a due to the traditional pulsar gap mechanism (Arons & Scharlemann, 1979;
Cheng et al., 1986; Muslimov & Harding, 2004; Wadiasingh et al., 2020).
Assuming that this emission site is at $\hat{r}\sim 100$, based on the
scalings presented in Section 2.3, the FRB luminosities are easily interpreted
even if the low frequency waves have a relatively low amplitude with $\delta
B_{0}\sim 10^{6}$ G.

In the second scenario, the emission site is in the reconnection current sheet
region just outside of the light cylinder, similar to Lyubarsky (2020).
However, the radiation mechanism is via coherent ICS by bunches rather than
oscillations of colliding magnetic islands within the current sheet as
proposed by Lyubarsky (2020). Here the bunches could be density fluctuations
in the turbulent relativistic particle outflows within the reconnection layer.
At this radius, the amplitude of low frequency waves is degraded significantly
with respect to that in the near-surface region. Since magnetic reconnection
is enhanced by a strong magnetic pulse that is eventually responsible for the
FRB, the magnetic energy density is enhanced by a factor of $b\equiv B_{\rm
pulse}/B_{\rm wind}$ in the current sheet region (Lyubarsky, 2020). Based on
magnetic flux conservation, the cross section of the magnetic tube is also
compressed by the same factor $b$. Since the length scale in the ${\bf B}$
direction is not changed, this leads to a local plasma density increased by a
factor of $\sim b$. According to Eqs.(1) and (4) of Lyubarsky (2020), one
estimates $b\sim 4\times 10^{4}$. As a result, we can estimate the local
charge density using Eq.(6) with $\hat{r}=\hat{r}_{\rm lc}$ (Eq.(9)) and
$\zeta\simeq b=4\times 10^{4}b_{4.6}$. Rewriting the equations presented in
Section 2.3 with $\hat{r}=\hat{r}_{\rm lc}$ (Eq.(9)), one obtains:

$\displaystyle U_{\rm ph}^{\rm lc}$ $\displaystyle\simeq$
$\displaystyle(3.5\times 10^{3}\ {\rm erg\ cm^{-3}})\delta
B_{0,6}^{2}P^{-2}R_{6}^{2}$ (28) $\displaystyle P_{e}^{\rm ICS,lc}$
$\displaystyle\simeq$ $\displaystyle(9.3\times 10^{-11}\ {\rm erg\
s^{-1}})\delta B_{0,6}^{2}P^{-2}f(\theta_{i})R_{6}^{2}$ (29) $\displaystyle
L^{\rm lc}$ $\displaystyle\simeq$ $\displaystyle(4.4\times 10^{38}\ {\rm erg\
s^{-1}})b_{4.6}^{2}N_{b,9}\gamma_{2.5}^{6}\nu_{9}^{-6}$ (30)
$\displaystyle\times B_{s,15}^{2}P^{-10}f(\theta_{i})\delta
B_{0,6}^{2}R_{6}^{8}$ $\displaystyle t_{c,b}^{\rm lc}$ $\displaystyle\simeq$
$\displaystyle(1.3\times 10^{-11}\ {\rm
s})b_{4.6}^{-1}\nu_{9}^{3}\gamma_{2.5}^{-1}B_{s,15}^{-1}$ (31)
$\displaystyle\times P^{6}[f(\theta_{i})]^{-1}\delta B_{0,6}^{-2}R_{6}^{-5}$
$\displaystyle E_{\parallel}^{\rm lc}$ $\displaystyle\simeq$
$\displaystyle(1.4\times 10^{6}\ {\rm
esu})b_{4.6}\nu_{9}^{-3}\gamma_{2.5}^{2}B_{s,15}P^{-6}$ (32) $\displaystyle
f(\theta_{i})\delta B_{0,6}^{2}R_{6}^{5}.$

One can see that the FRB luminosities can be still explained with $\delta
B_{0}\sim 10^{6}$ G but a somewhat larger bunch number $N_{b}\sim 10^{9}$.

Note that in the magnetic reconnection site, the reconnection-driven electric
field is perpendicular to ${\bf B}$. This is different from the first scenario
where $E_{\parallel}$ is along the ${\bf B}$ field direction. However, since
we do not know the geometric configuration of the FRB central engine, the same
FRB may be interpreted by the two scenarios with different viewing geometries:
The inner magnetospheric scenario has the viewing angle sweep the open field
line region, whereas the reconnection scenario has the viewing angle sweep the
equatorial plane (see Fig. 1).

###  3.2 Energy budget

In our model, the energy budget to power FRBs is carried by particles rather
than by low frequency waves. The latter only carries a small amount of energy.
Its role is to provide seed photons for coherent ICS emission to operate.
Since $\delta B\propto r^{-1}$ and $B\propto r^{-3}$, with the nominal
parameters $\delta B_{0}\sim 10^{6}$ G and $B_{s}\sim 10^{15}$ G, one
estimates $\delta B^{\rm lc}/B^{\rm lc}=4.8\times 10^{-6}PR_{6}^{-1}\ll 1$ at
the light cylinder. As a result, the enhancement of the magnetospheric Thomson
cross section (which is relevant when $\delta B>B$, Beloborodov 2021)222In a
weak background $\bf B$ field, the Thomson scattering cross section is
enlarged by a factor of $a^{2}$ for strong waves for $a\gg 1$, where $a$ is
the amplitude parameter (Yang & Zhang, 2020). However, with the presence of a
strong background field, this effect is suppressed until the wave vector
$\delta B$ becomes stronger than the background $B$.  does not occur for the
low-frequency waves so that they can propagate freely within the magnetar
magnetosphere333The constraint on the upscattered waves still applies,
especially if the emission has a very high luminosity (Beloborodov, 2021). The
problem is alleviated if the emission altitude is high (as envisaged in our
scenarios) and when radiation pressure (Ioka, 2020; Wang et al., 2021) or
ponderomotive pre-acceleration of background plasma (Lyutikov, 2021b) by the
intense FRB emission field are considered.. The energy of the emitting
particles comes from the main magnetic pulse driven by the violent event
itself. The energy that powers the FRB is first carried by the Alfvén waves,
which advect particles to a large altitude where an $E_{\parallel}$ is
developed in a charge-starved region in the magnetar magnetosphere (Kumar &
Bošnjak, 2020; Muslimov & Harding, 2004; Wadiasingh et al., 2020) or in
magnetic reconnection sites in the current sheet just outside the
magnetosphere (Lyubarsky, 2020; Kalapotharakos et al., 2018; Philippov &
Spitkovsky, 2018). Particles are accelerated in these $E_{\parallel}$ fields
so that their ultimate energy comes from the magnetic pulse itself.

###  3.3 Polarization properties

Since only the X-mode low frequency waves can propagate in the magnetosphere,
the incident photon is linearly polarized and has an electric field vector
perpendicular to the local magnetic field direction. The outgoing photon in an
ICS process is also linearly polarized. According to Herold (1979), only the
cross section from mode 1 to mode 1 is related to the $\sin^{2}\theta_{i}$
term in Eq.(17). All the other scattering modes (1 to 2, 2 to 1, 2 to 2)
involve the $[{\omega^{\prime}}^{2}/(\omega^{\prime}+\omega_{\rm
B})^{2}+{\omega^{\prime}}^{2}/(\omega^{\prime}-\omega_{\rm B})^{2}]$ term,
which is negligibly small for our problem. The 1 mode has both an X-mode
component and an O-mode component. The latter cannot propagate. As a result,
our model predicts linearly polarized X-mode emission with a $\sim$ 100%
degree of linear polarization in the X mode. Circular polarization may develop
under special conditions (Xu et al., 2000). This is in general consistent with
the FRB observations (e.g. Gajjar et al., 2018; Michilli et al., 2018; Cho et
al., 2020; Day et al., 2020; Luo et al., 2020).

The duration of an FRB is defined by the longer of the emission duration and
the timescale for the line-of-sight to sweep the bundle of the magnetic field
lines where emitting particles flow out (Wang et al., 2019, 2021). For a slow
rotator and an emission site far from the neutron star surface, the
polarization angle would be nearly constant across an individual burst, as
observed in some repeating FRBs (e.g. Michilli et al., 2018). For a rapid
rotator and an emission site closer to the neutron star surface, the
polarization angle may display diverse swing features, as observed in some
other repeating FRBs (Luo et al., 2020). Our model can account for both
patterns.

###  3.4 Narrow spectrum and frequency down-drifting

The characteristic frequency of coherent ICS emission Eq.(3) depends on
$\omega_{0}$, $\gamma$ and $\theta_{i}$. For the coherent mechanism to work,
particles need to continuously tap energy from $E_{\parallel}$ so that they
are likely emitting in the radiation-reaction-limited regime. As a result,
$\gamma$ may remain roughly constant during the emission phase. At the large
emission radius as envisaged in the two scenarios, $\theta_{i}$ varies little
as the line of sight sweeps across different field lines. For a seismic
oscillation originating from the neutron star crust, the oscillation frequency
$\omega_{0}$ depends on the resonance frequency of the crust which may have a
characteristic value. As a result, the characteristic frequency $\nu$ has a
narrow value range, unlike curvature radiation which is an intrinsically broad
band mechanism (Yang & Zhang, 2018)444A relatively narrow spectrum can be
obtained if oppositely charged particles are spatially separated (Yang et al.,
2020).. This is consistent with the narrow-band emission as observed from
repeating FRBs (Spitler et al., 2016; CHIME/FRB Collaboration et al., 2021).

The characteristic Lorentz factor of electrons may slightly decrease with
radius. If so, the outgoing ICS frequency would decrease with radius and give
a radius-to-frequency mapping within the ICS model (Qiao & Lin, 1998).
Alternatively, as the seismic waves damp, $\omega_{0}$ would gradually
decrease with a decreasing amplitude. This would also result in a decreasing
$\nu$ of FRB emission with reduced amplitude. Both factors may offer an
explanation to the observed sub-pulse frequency down-drifting pattern (also
called the “sad trombone” effect) (Hessels et al., 2019; CHIME/FRB
Collaboration et al., 2019a, b)555Other interpretations include radius-to-
frequency mapping for curvature radiation (Wang et al., 2019), decreasing
magnetic field strengths in synchrotron maser shocks (Metzger et al., 2019),
or external free-free absorption (Kundu & Zhang, 2021)..

###  3.5 Required degree of coherence

Since the ICS emission power of individual electrons is much higher than that
of curvature radiation, the required degree of coherence in the bunch is much
lower. In the scalings presented in Section 2.3, the net charge parameter
$\zeta$ is normalized to unity. In other words, a bunch with the Goldreich-
Julian net charge density can produce the high brightness temperature of FRB
emission. In fact, the model even allows the bunch charge density to be below
the Goldreich-Julian density, with the luminosity compensated with a somewhat
larger total number of bunches, $N_{b}$, which is currently normalized to a
small value. The only requirement to make a bunch is that there are
fluctuations of net charge density in space with respect to the background
number density (Yang & Zhang, 2018). This alleviates traditional criticisms on
the bunching mechanisms within the context of curvature radiation regarding
the formation and maintenance of bunches (Melrose, 1978). Since the mechanism
can apply to magnetospheres with a low plasma density (see also Wadiasingh et
al. 2020), another major criticism on the bunching mechanism, i.e. the
suppression of coherent emission flux due to the dense plasma effect (Gil et
al., 2004; Lyubarsky, 2021), is also alleviated.

##  4 Conclusions and discussion

We have proposed a new family of FRB emission model invoking coherent ICS
emission by bunches within or just outside of the magnetosphere of a flaring
magnetar. The key assumption is that the oscillations of the magnetar crust
would excite oscillations of charges in the magnetosphere near the neutron
star surface, which would excite low-frequency electromagnetic waves with an
antenna mechanism. The X-mode of these waves could propagate freely in the
magnetosphere. Bunched relativistic particles in the outer magnetosphere or in
the current sheet region outside of the magnetosphere could upscatter the low
frequency waves coherently to power FRBs observed in the $\sim$ GHz band
(Eq.(3)). For standard parameters, the ICS power of individual particles
(Eq.(19)) is much greater than that of curvature radiation (Eq.(21)). As a
result, the typical FRB luminosity can be easily reproduced with a low charge
density bunch (plasma density of the order of or even lower than the
Goldreich-Julian charge density) and a moderate number of bunches (Eqs.(25)
and (30)). An $E_{\parallel}$ is needed to sustain the coherent radiation
(similar to the curvature radiation mechanism), but such an $E_{\parallel}$ is
expected in charge starved region in the outer magnetar magnetosphere or in
the current sheet outside the magnetosphere. This model can account for
several observational properties of repeating FRBs, including nearly 100%
polarization degree of emission, constant or varying polarization angle across
each burst, narrow emission spectrum, and frequency downdrifting. The low
degree of coherence also alleviates several criticisms to the bunched coherent
curvature radiation mechanism, including the formation and maintenance of
bunches as well as the plasma suppression effect.

The vacuum treatment discussed in this paper shares some essential features
with a broader family of models that invoke relativistic particles scattering
off various plasma waves. For example, Lyubarskii (1996) presented a detailed
treatment of induced scattering of relativistic particles off longitudinal
subluminal plasma waves and showed that superluminal transverse
electromagnetic waves can be efficiently generated to power pulsar radio
emission with the characteristic radius-to-frequency mapping. Lyutikov (2021a)
studied the scattering of relativistic particles off Alfvén wave wigglers and
showed a narrow characteristic frequency and a high efficiency, similar to the
features of the vacuum model presented in this paper. All these suggest that
ICS of relativistic particles against various forms of low-frequency waves
could be in general an attractive mechanism to power magnetospheric FRB
emission from magnetars.

The author acknowledges the referee for critical comments that helped to
improve the presentation of the paper, Yuan-Pei Yang, Pawan Kumar, and Wenbin
Lu for excellent comments, Guojun Qiao for discussion and encouragement, and
Yuanhong Qu for discussion and the help to draw the figure.

##  Appendix A Propagation of low frequency waves in a magnetar magnetosphere

The low frequency electromagnetic waves conjectured in this paper need to
propagate through the closed field line region of the magnetar to reach the
two FRB emission regions (Fig. 1). Since the particles residing in the closed
field line region of a magnetar magnetosphere are believed to have non-
relativistic motion, the dispersion relation of wave propagation may be
treated under the assumption of a cold plasma. For a cold, magnetized medium,
the conductivity is a tensor $\sigma_{ij}$ so that ${\bf
j}=\sigma_{ij}\cdot{\bf E}$, where ${\bf E}$ is the electric field vector of
the wave. Define the dielectric tensor

$\epsilon_{ij}=\delta_{ij}+4\pi i\sigma_{ij}/\omega,$ (A1)

where $\omega$ is the angular frequency of a wave. The fourth Maxwell’s
equation becomes ${\bf n}\times({\bf n}\times{\bf E})=-\epsilon_{ij}\cdot{\bf
E}$, which defines the dispersion relations for wave propagation. In general,
defining $\theta$ as the angle between the wave number vector ${\bf k}$ and
the magnetic field vector ${\bf B}$, one can write the Maxwell response tensor
as (e.g. Boyd & Sanderson, 2003)

$\epsilon_{ij}\equiv\left(\begin{array}[]{ccc}S-n^{2}\cos^{2}\theta&-iD&n^{2}\cos\theta\sin\theta\\\
iD&S-n^{2}&0\\\
n^{2}\cos\theta\sin\theta&0&P-n^{2}\sin^{2}\theta\end{array}\right)$ (A2)

where ${\bf B}$ is defined in the $\hat{\bf z}$ direction. Here

$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{1}{2}(R+L)=1-\frac{\omega_{p}^{2}(\omega^{2}+\Omega_{i}\Omega_{e})}{(\omega^{2}-\Omega_{i}^{2})(\omega^{2}-\Omega_{e}^{2})},$
(A3) $\displaystyle D$ $\displaystyle=$
$\displaystyle\frac{1}{2}(R-L)=\frac{\omega_{p}^{2}\omega(\Omega_{i}+\Omega_{e})}{(\omega^{2}-\Omega_{i}^{2})(\omega^{2}-\Omega_{e}^{2})},$
(A4) $\displaystyle R$ $\displaystyle=$ $\displaystyle
1-\frac{\omega_{p}^{2}}{(\omega+\Omega_{i})(\omega+\Omega_{e})},$ (A5)
$\displaystyle L$ $\displaystyle=$ $\displaystyle
1-\frac{\omega_{p}^{2}}{(\omega-\Omega_{i})(\omega-\Omega_{e})},$ (A6)
$\displaystyle P$ $\displaystyle=$ $\displaystyle
1-\frac{\omega_{p}^{2}}{\omega^{2}},$ (A7)

where $\omega_{p}$ is the plasma frequency,
$\Omega_{e}=q_{e}B/m_{e}c=-eB/m_{e}c=-\omega_{\rm B}$ is the electron gyration
frequency (which is the Larmor frequency $\omega_{\rm B}$ with a negative
sign), and $\Omega_{i}=q_{i}B/m_{i}c$ is the ion gyration frequency (positive
sign). For an ion with atomic number $Z$ and mass number $A$, one has
$\Omega_{i}=ZeB/Am_{p}c\ll\omega_{\rm B}$. For hydogen, one has
$\Omega_{i}=eB/m_{p}c$. For an electron positron plasma, one has
$\Omega_{i}=eB/m_{e}c=-\Omega_{e}=\omega_{\rm B}$. The general dispersion
relation for cold plasma waves is

$An^{4}-Bn^{2}+C=0,$ (A8)

where

$\displaystyle A$ $\displaystyle=$ $\displaystyle
S\sin^{2}\theta+P\cos^{2}\theta,$ (A9) $\displaystyle B$ $\displaystyle=$
$\displaystyle RL\sin^{2}\theta+PS(1+\cos^{2}\theta)$ (A10) $\displaystyle C$
$\displaystyle=$ $\displaystyle PRL.$ (A11)

The propagation of the waves is prohibited at $k=0$ (cut-offs) or at
$k\rightarrow\infty$ (resonances) for certain propagation angles. At principle
resonances ($k\rightarrow\infty$ and $\theta_{\rm res}=0$ or $90^{\rm o}$,
where the resonant angle $\theta_{\rm res}$ is defined by $\tan^{2}\theta_{\rm
res}=-P/S$), propagation in all directions is prohibited.

For the case of electromagnetic waves propagating along a magnetic field line,
i.e. ${\bf k}\parallel{\bf B}$, the dispersion relations become $n^{2}=R$ and
$n^{2}=L$. Making $R=0$ and $L=0$, one can define two cut-off frequencies for
the R-mode and L-mode, respectively, i.e.

$\displaystyle\omega_{\rm R}$ $\displaystyle\equiv$
$\displaystyle\left[\omega_{p}^{2}+\frac{(\Omega_{i}-\Omega_{e})^{2}}{4}\right]^{1/2}-\frac{(\Omega_{i}+\Omega_{e})}{2}=\sqrt{\omega_{p}^{2}+\omega_{\rm
B}^{2}},$ (A12) $\displaystyle\omega_{\rm L}$ $\displaystyle\equiv$
$\displaystyle\left[\omega_{p}^{2}+\frac{(\Omega_{i}-\Omega_{e})^{2}}{4}\right]^{1/2}+\frac{(\Omega_{i}+\Omega_{e})}{2}=\sqrt{\omega_{p}^{2}+\omega_{\rm
B}^{2}}.$ (A13)

Hereafter the second equation in each expression is for a pair plasma, where
$\Omega_{i}=-\Omega_{e}=\omega_{\rm B}$ is adopted. The principle resonances
can be obtained from $R\rightarrow\infty$ and $L\rightarrow\infty$, which
gives

$\displaystyle\omega_{\rm res,R}$ $\displaystyle=$
$\displaystyle-\Omega_{e}=\omega_{\rm B},$ (A14) $\displaystyle\omega_{\rm
res,L}$ $\displaystyle=$ $\displaystyle\Omega_{i}=\omega_{\rm B}.$ (A15)

Electromagnetic waves are transparent above the cut-off frequencies or below
the resonances. So for ${\bf k}\parallel{\bf B}$, the condition for wave
propagation in a cold, magnetized plasma is

$\omega>\sqrt{\omega_{p}^{2}+\omega_{\rm B}^{2}},~{}~{}~{}{\rm
or}~{}~{}~{}\omega<\omega_{\rm B}.$ (A16)

For the case of electromagnetic waves propagating in the direction
perpendicular to the field lines, i.e. ${\bf k}\perp{\bf B}$, one should
consider two modes: the O-mode with ${\bf E}\parallel{\bf B}$ and the X-mode
with ${\bf E}\perp{\bf B}$. The O-mode is similar to the case as if no
magnetic field exists (the electron moving in response of the ${\bf E}$ field
does not feel the Lorentz force from ${\bf B}$). Its cut-off frequency is
defined by $P=0$, i.e. $\omega=\omega_{p}$. The waves cannot propagate at
$\omega<\omega_{p}$.

The X-mode propagation is more complicated. Taking $\theta=90^{\rm o}$, the
dispersion equation (Eq.(A8)) is simplified to

$n^{2}=\frac{RL}{S}.$ (A17)

The cut-off frequencies are again defined by $R=0$ (i.e. $\omega=\omega_{\rm
R}$) and $L=0$ (i.e. $\omega=\omega_{\rm L}$), and the resonance is defined by
$S=0$, with the principle resonances defined with an additional condition
$\theta_{\rm res}=\pi/2$. This demands
$\omega^{4}-\omega^{2}(\omega_{p}^{2}+\Omega_{i}^{2}+\Omega_{e}^{2})-\Omega_{i}\Omega_{e}(\omega^{2}-\Omega_{i}\Omega_{e})=0$,
which defines two principle resonances

$\displaystyle\omega^{2}=\left(\frac{\omega_{p}^{2}+\Omega_{i}^{2}+\Omega_{e}^{2}}{2}\right)\left[1\pm\left(1+\frac{4\Omega_{i}\Omega_{e}(\omega_{p}^{2}-\Omega_{i}\Omega_{e})}{(\omega_{p}^{2}+\Omega_{i}^{2}+\Omega_{e}^{2})^{2}}\right)^{1/2}\right].$
(A18)

For a pair plasma ($\Omega_{i}=-\Omega_{e}=\omega_{\rm B}$), one obtains the
two (upper and lower) hybrid resonances as

$\displaystyle\omega_{\rm uh}^{2}$ $\displaystyle=$
$\displaystyle\omega_{p}^{2}+\omega_{\rm B}^{2},$ (A19)
$\displaystyle\omega_{\rm lh}^{2}$ $\displaystyle=$ $\displaystyle\omega_{\rm
B}^{2}.$ (A20)

The cut-off bands for X-mode in the case of ${\bf k}\perp{\bf B}$ is
$\omega_{\rm uh}<\omega<\omega_{\rm R}$ and $\omega_{\rm
lh}<\omega<\omega_{\rm L}$. As shown above, for a pair plasma, one has
$\omega_{\rm R}=\omega_{\rm uh}=\omega_{\rm
L}=\sqrt{\omega_{p}^{2}+\omega_{\rm B}^{2}}$ and $\omega_{\rm lh}=\omega_{\rm
B}$. One therefore draws the conclusion that the cut-off band for a pair
plasma in ${\bf k}\perp{\bf B}$ X mode is identical to Eq.(A16) derived for
the ${\bf k}\parallel{\bf B}$ case.

Since both $\theta=0$ and $\theta=90^{\rm o}$ cases (the latter for X-mode
only) have the identical wave propagation condition, Eq.(A16) should apply to
the oblique case with an arbitrary $\theta$ value for the X-mode waves. O-mode
waves cannot propagate below $\omega_{p}$.

## References

  * Arons & Barnard (1986) Arons, J., & Barnard, J. J. 1986, ApJ, 302, 120
  * Arons & Scharlemann (1979) Arons, J., & Scharlemann, E. T. 1979, ApJ, 231, 854
  * Beloborodov (2017) Beloborodov, A. M. 2017, ApJ, 843, L26
  * Beloborodov (2020) —. 2020, ApJ, 896, 142
  * Beloborodov (2021) —. 2021, arXiv e-prints, arXiv:2108.07881
  * Beloborodov & Thompson (2007) Beloborodov, A. M., & Thompson, C. 2007, ApJ, 657, 967
  * Bochenek et al. (2020) Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Nature, 587, 59
  * Boyd & Sanderson (2003) Boyd, T. J. M., & Sanderson, J. J. 2003, The Physics of Plasmas
  * Chen et al. (2020) Chen, A. Y., Yuan, Y., Beloborodov, A. M., & Li, X. 2020, arXiv:2010.15619
  * Cheng et al. (1986) Cheng, K. S., Ho, C., & Ruderman, M. 1986, ApJ, 300, 522
  * CHIME/FRB Collaboration et al. (2019a) CHIME/FRB Collaboration, Amiri, M., Bandura, K., et al. 2019a, Nature, 566, 235
  * CHIME/FRB Collaboration et al. (2019b) CHIME/FRB Collaboration, Andersen, B. C., Bandura, K., et al. 2019b, ApJ, 885, L24
  * CHIME/FRB Collaboration et al. (2020) CHIME/FRB Collaboration, Andersen, B. C., Bandura, K. M., et al. 2020, Nature, 587, 54
  * CHIME/FRB Collaboration et al. (2021) CHIME/FRB Collaboration, :, Amiri, M., et al. 2021, arXiv e-prints, arXiv:2106.04352
  * Cho et al. (2020) Cho, H., Macquart, J.-P., Shannon, R. M., et al. 2020, ApJ, 891, L38
  * Cooper & Wijers (2021) Cooper, A. J., & Wijers, R. A. M. J. 2021, arXiv e-prints, arXiv:2108.07818
  * Cordes & Chatterjee (2019) Cordes, J. M., & Chatterjee, S. 2019, ARA&A, 57, 417
  * Day et al. (2020) Day, C. K., Deller, A. T., Shannon, R. M., et al. 2020, arXiv e-prints, arXiv:2005.13162
  * Dehman et al. (2020) Dehman, C., Viganò, D., Rea, N., et al. 2020, ApJ, 902, L32
  * Fung & Kuijpers (2004) Fung, P. K., & Kuijpers, J. 2004, A&A, 422, 817
  * Gajjar et al. (2018) Gajjar, V., Siemion, A. P. V., Price, D. C., et al. 2018, ApJ, 863, 2
  * Gil et al. (2004) Gil, J., Lyubarsky, Y., & Melikidze, G. I. 2004, ApJ, 600, 872
  * Goldreich & Julian (1969) Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869
  * Herold (1979) Herold, H. 1979, Phys. Rev. D, 19, 2868
  * Hessels et al. (2019) Hessels, J. W. T., Spitler, L. G., Seymour, A. D., et al. 2019, ApJ, 876, L23
  * Ioka (2020) Ioka, K. 2020, ApJ, 904, 15
  * Kalapotharakos et al. (2018) Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44
  * Katz (2014) Katz, J. I. 2014, Phys. Rev. D, 89, 103009
  * Kumar & Bošnjak (2020) Kumar, P., & Bošnjak, Ž. 2020, MNRAS, 494, 2385
  * Kumar et al. (2017) Kumar, P., Lu, W., & Bhattacharya, M. 2017, MNRAS, 468, 2726
  * Kumar & Zhang (2015) Kumar, P., & Zhang, B. 2015, Phys. Rep., 561, 1
  * Kundu & Zhang (2021) Kundu, E., & Zhang, B. 2021, MNRAS, 508, L48
  * Lin et al. (2020) Lin, L., Zhang, C. F., Wang, P., et al. 2020, Nature, 587, 63
  * Lorimer et al. (2007) Lorimer, D. R., Bailes, M., McLaughlin, M. A., Narkevic, D. J., & Crawford, F. 2007, Science, 318, 777
  * Lu & Kumar (2018) Lu, W., & Kumar, P. 2018, MNRAS, 477, 2470
  * Lu et al. (2019) Lu, W., Kumar, P., & Narayan, R. 2019, MNRAS, 483, 359
  * Lu et al. (2020) Lu, W., Kumar, P., & Zhang, B. 2020, MNRAS, 498, 1397
  * Luo et al. (2020) Luo, R., Wang, B. J., Men, Y. P., et al. 2020, Nature, 586, 693
  * Lyubarskii (1996) Lyubarskii, Y. E. 1996, A&A, 308, 809
  * Lyubarskii & Petrova (1998) Lyubarskii, Y. E., & Petrova, S. A. 1998, A&A, 333, 181
  * Lyubarsky (2014) Lyubarsky, Y. 2014, MNRAS, 442, L9
  * Lyubarsky (2020) —. 2020, ApJ, 897, 1
  * Lyubarsky (2021) —. 2021, Universe, 7, 56
  * Lyutikov (2021a) Lyutikov, M. 2021a, arXiv e-prints, arXiv:2102.07010
  * Lyutikov (2021b) —. 2021b, arXiv e-prints, arXiv:2110.08435
  * Margalit et al. (2020) Margalit, B., Beniamini, P., Sridhar, N., & Metzger, B. D. 2020, ApJ, 899, L27
  * Melikidze et al. (2000) Melikidze, G. I., Gil, J. A., & Pataraya, A. D. 2000, ApJ, 544, 1081
  * Melrose (1978) Melrose, D. B. 1978, ApJ, 225, 557
  * Melrose (2017) —. 2017, Reviews of Modern Plasma Physics, 1, 5
  * Metzger et al. (2019) Metzger, B. D., Margalit, B., & Sironi, L. 2019, MNRAS, 485, 4091
  * Michilli et al. (2018) Michilli, D., Seymour, A., Hessels, J. W. T., et al. 2018, ArXiv e-prints
  * Muslimov & Harding (2004) Muslimov, A. G., & Harding, A. K. 2004, ApJ, 606, 1143
  * Pacholczyk (1970) Pacholczyk, A. G. 1970, Radio astrophysics. Nonthermal processes in galactic and extragalactic sources
  * Petroff et al. (2019) Petroff, E., Hessels, J. W. T., & Lorimer, D. R. 2019, A&A Rev., 27, 4
  * Philippov & Spitkovsky (2018) Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94
  * Plotnikov & Sironi (2019) Plotnikov, I., & Sironi, L. 2019, MNRAS, 485, 3816
  * Qiao & Lin (1998) Qiao, G. J., & Lin, W. P. 1998, A&A, 333, 172
  * Qiao et al. (2001) Qiao, G. J., Liu, J. F., Zhang, B., & Han, J. L. 2001, A&A, 377, 964
  * Qu & Zhang (2021) Qu, Y., & Zhang, B. 2021, in preparation
  * Rowe (1995) Rowe, E. T. 1995, A&A, 296, 275
  * Ruderman & Sutherland (1975) Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51
  * Spitler et al. (2016) Spitler, L. G., Scholz, P., Hessels, J. W. T., et al. 2016, Nature, 531, 202
  * Thompson & Duncan (2001) Thompson, C., & Duncan, R. C. 2001, ApJ, 561, 980
  * Wadiasingh et al. (2020) Wadiasingh, Z., Beniamini, P., Timokhin, A., et al. 2020, ApJ, 891, 82
  * Wang & Lai (2020) Wang, J.-S., & Lai, D. 2020, ApJ, 892, 135
  * Wang et al. (2018) Wang, W., Luo, R., Yue, H., et al. 2018, ApJ, 852, 140
  * Wang et al. (2019) Wang, W., Zhang, B., Chen, X., & Xu, R. 2019, ApJ, 876, L15
  * Wang et al. (2021) Wang, W.-Y., Yang, Y.-P., Niu, C.-H., et al. 2021, ApJ, submitted
  * Xia et al. (1985) Xia, X. Y., Qiao, G. J., Wu, X. J., & Hou, Y. Q. 1985, A&A, 152, 93
  * Xiao et al. (2021) Xiao, D., Wang, F., & Dai, Z. 2021, Science China Physics, Mechanics, and Astronomy, 64, 249501
  * Xu et al. (2000) Xu, R. X., Liu, J. F., Han, J. L. & Qiao, G. J. 2000, ApJ, 535, 354
  * Yang & Zhang (2018) Yang, Y.-P., & Zhang, B. 2018, ApJ, 868, 31
  * Yang & Zhang (2020) —. 2020, ApJ, 892, L10
  * Yang & Zhang (2021) —. 2021, ApJ, 919, 89
  * Yang et al. (2020) Yang, Y.-P., Zhu, J.-P., Zhang, B., & Wu, X.-F. 2020, ApJ, 901, L13
  * Yuan et al. (2020) Yuan, Y., Beloborodov, A. M., Chen, A. Y., & Levin, Y. 2020, ApJ, 900, L21
  * Zhang (2020) Zhang, B. 2020, Nature, 587, 45
  * Zhang et al. (1999) Zhang, B., Hong, B. H., & Qiao, G. J. 1999, ApJ, 514, L111


# Distributed Expectation Propagation Detection for Cell-Free Massive MIMO

Hengtao He1, Hanqing Wang2, Xianghao Yu1, Jun Zhang3, S.H. Song1, and Khaled
B. Letaief1 This work is supported by the Hong Kong Research Grant Council
under Grant No. 16212120. 1Department of ECE, The Hong Kong University of
Science and Technology,
Kowloon, Hong Kong, E-mail: {eehthe, eeyu, eeshsong<EMAIL_ADDRESS>2Huawei
Technologies Company Ltd., Shanghai 201206, China, E-mail:
<EMAIL_ADDRESS>3Department of EIE, The Hong Kong Polytechnic
University, Hong Kong, E-mail<EMAIL_ADDRESS>
###### Abstract

In cell-free massive MIMO networks, an efficient distributed detection
algorithm is of significant importance. In this paper, we propose a
distributed expectation propagation (EP) detector for cell-free massive MIMO.
The detector is composed of two modules, a nonlinear module at the central
processing unit (CPU) and a linear module at the access point (AP). The turbo
principle in iterative decoding is utilized to compute and pass the extrinsic
information between modules. An analytical framework is then provided to
characterize the asymptotic performance of the proposed EP detector with a
large number of antennas. Simulation results will show that the proposed
method outperforms the distributed detectors in terms of bit-error rate.

##  I Introduction

Massive multiple-input multiple-output (MIMO) systems have been regarded as a
key enabling technology for $5$G because of its high spectral efficiency (SE),
energy efficiency and link reliability [1]. However, the SE gain becomes
marginal for the cell-edge user equipments (UEs). To address this problem, a
novel network architecture, namely cell-free massive MIMO, was proposed [2,
3]. This is a disruptive emerging technology which has been recognized as a
crucial and core technology for the upcoming beyond 5G and 6G networks[4].
Such technology is expected to bring important benefits, including huge data
throughput, ultra-low latency, ultra-high reliability, a huge increase in
mobile energy efficiency, and ubiquitous uniform coverage. The fundamental
idea is to deploy a large number of distributed access points (APs) connected
to a central processing unit (CPU) to serve all the distributed UEs in a wide
area. In particular, each AP serves all UEs in the same time-frequency
resource block via a time-division duplex (TDD) mode. Compared to conventional
colocated massive MIMO, cell-free networks offer more uniform connectivity for
all UEs thanks to the macro-diversity gain obtained from the distributed
antennas. However, the assumption that each AP serves all UEs renders the
system not scalable and incurs huge power and computational resource
consumption for decoding, especially at UEs with low signal-to-interference-
noise-ratios. To tackle the scalability issue, a _user-centric_ dynamic
cooperation clustering (DCC) scheme [5] was introduced [6], where each user is
only served by a subset of APs.

An efficient data detection algorithm is highly desired in large-scale and
complex networks, such as cell-free massive MIMO. In this aspect, some early
attempts were made on centralized algorithms where the detection is totally
implemented at the CPU with the received pilots and data signals reported from
all APs [2, 7]. However, the computational overhead of such a centralized
detection scheme is prohibitively high when the network size becomes large. To
address this challenge, distributed detectors have been recently investigated.
In [8], one centralized and three distributed receivers with different levels
of cooperation among APs were compared in terms of SE. Unfortunately, the
distributed receivers investigated in [8] are linear receivers and therefore
highly suboptimal in terms of the bit-error rate (BER) performance. Therefore,
it is of great importance to develop a distributed and _non-linear_ receiver
to achieve a better BER performance.

In this paper, we propose a non-linear detector for cell-free massive MIMO
networks, which is derived based upon the EP principle [9] with a distributed
approach [10, 11]. Specifically, by adopting the linear minimum mean-square
error (MMSE) estimator, the APs first detect the symbols with the local
channel state information and transfer the posterior mean and variance
estimates to the CPU. Then, the extrinsic information for each AP is computed
and integrated at the CPU by utilizing maximum-ratio combining (MRC).
Subsequently, the CPU uses the posterior mean estimator to refine the
detection and the extrinsic information is transferred to each AP from the CPU
via the fronthaul. Compared to other distributed linear detectors, the
proposed distributed EP detector improves the detection performance by
introducing more computation overhead mainly at the computationally powerful
CPU. Simulation results will demonstrate that the proposed method outperforms
existing distributed detectors and even the centralized MMSE detector in terms
of BER.

_Notations_ —For any matrix $\mathbf{A}$, $\mathbf{A}^{H}$ and
${\mathrm{tr}}(\mathbf{A})$ denote the conjugate transpose and trace of
$\mathbf{A}$, respectively. In addition, $\mathbf{I}$ is the identity matrix
and $\mathbf{0}$ is the zero matrix. We use $\mathrm{D}z$ to denote the real
Gaussian integration measure. That is,

$\mathrm{D}z=\phi(z)dz,\quad\mathrm{where}\quad\phi(z)\triangleq\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}.$

A complex Gaussian distribution with mean $\boldsymbol{\mu}$ and covariance
$\boldsymbol{\Omega}$ can be described by the probability density function,

$\mathcal{N}_{\mathbb{C}}(\mathbf{z};\boldsymbol{\mu},\boldsymbol{\Omega})=\frac{1}{\mathrm{det}(\pi\boldsymbol{\Omega})}e^{-(\mathbf{z}-\boldsymbol{\mu})^{H}\boldsymbol{\Omega}^{-1}(\mathbf{z}-\boldsymbol{\mu})}.$

The remaining part of this paper is organized as follows. Section II
formulates the cell-free massive MIMO detection problem. The distributed EP
detector is proposed in Section III and an analytical framework is provided in
Section IV. Numerical results are then presented in Section V and Section VI
concludes the paper.

##  II System Model

In this section, we first present the system model and formulate the cell-free
massive MIMO detection problem. Then, four commonly-adopted receivers are
briefly introduced.

###  II-A Cell-Free Massive MIMO

As illustrated in Fig. 1, we consider a cell-free massive MIMO network with
$L$ distributed APs, each equipped with $N$ antennas to serve $K$ single-
antenna UEs. The system can be AP-centric or user-centric. All APs are
connected to a CPU that has a high computational capability. Denote
${\mathbf{h}}_{kl}\sim\mathcal{N}_{{\mathbb{C}}}(\mathbf{0},{\mathbf{R}}_{kl})$
as the channel between the $k$-th user and the $l$-th AP, where
${\mathbf{R}}_{kl}\in{\mathbb{C}}^{N\times N}$ is the spatial correlation
matrix and $\beta_{k,l}=\mathrm{tr}({\mathbf{R}}_{kl})/N$ is the large-scale
fading coefficient involving the geometric path loss and shadowing. In the
uplink data transmission phase, we consider
$\mathcal{M}_{k}\subset\\{1,\ldots,L\\}$ as the subset of APs that serve the
$k$-th UE and define the DCC matrices ${\mathbf{D}}_{kl}$ based on
$\mathcal{M}_{k}$ as

${\mathbf{D}}_{kl}=\left\\{\begin{aligned} &{\mathbf{I}}_{N}\quad\mathrm{if}\
l\in\mathcal{M}_{k}\\\ &\mathbf{0}_{N\times N}\quad\mathrm{if}\
l\notin\mathcal{M}_{k},\end{aligned}\right.$ (1)

which is constructed by the DCC strategy. Furthermore, we define
$\mathcal{D}_{l}$ as the set consisting of the UE indices that are served by
$l$-th AP :

$\mathcal{D}_{l}=\bigg{\\{}k:\mathrm{tr}({\mathbf{D}}_{kl})\geq
1,k\in\\{1,\ldots,K\\}\bigg{\\}},$ (2)

where the cardinality of $\mathcal{D}_{l}$ is denoted as $|\mathcal{D}_{l}|$.
When each AP serves all users, $\mathcal{D}_{l}=K$ accordingly. We assume that
perfect channel state information (CSI) is available at the local APs.
Therefore, the received signal at the $l$-th AP is given by,

${\mathbf{y}}_{l}=\sum_{k=1}^{|\mathcal{D}_{l}|}\sqrt{p_{k}}{\mathbf{h}}_{kl}x_{k}+{\mathbf{n}}_{l},$
(3)

where $x_{k}\in{\mathbb{C}}$ is the transmitted symbol drawn from the $M$-QAM
constellation, and $p_{k}>0$ is the transmit power at the $k$-th UE. The
additive noise at the $l$-th AP is denoted as
${\mathbf{n}}_{l}\sim\mathcal{N}_{{\mathbb{C}}}(\mathbf{0}_{N},\sigma^{2}{\mathbf{I}}_{N})$.
Let ${\mathbf{x}}=[x_{1},\ldots,x_{K}]$ denote the transmitted vector from all
UEs and
${\mathbf{h}}_{l}=[{\mathbf{h}}_{1l},\ldots,{\mathbf{h}}_{Kl}]^{T}\in{\mathbb{C}}^{N\times
K}$ as the channel of the AP $l$ to all UEs. If the $k$-th UE is not
associated with the $l$-th AP, the channel vector
${\mathbf{h}}_{kl}=\mathbf{0}$ accordingly. Furthermore, we denote
${\mathbf{H}}=[{\mathbf{h}}_{k}^{T},\ldots,{\mathbf{h}}_{K}^{T}]^{T}\in{\mathbb{C}}^{LN\times
K}$ as the channel matrix between all UEs and APs. The uplink detection
problem for cell-free massive MIMO is to detect the transmitted data
${\mathbf{x}}$ based on the received signals
${\mathbf{y}}_{l}\,(l=1,2,\ldots,L)$, channel matrix ${\mathbf{H}}$, and noise
power $\sigma^{2}$.

Figure 1: . A typical cell-free massive MIMO system.

###  II-B Linear Receivers

For the cell-free massive MIMO receiver design, four receivers with different
levels of cooperation among APs were introduced and compared in terms of SE
[8].

  * •

Level $4$ is a fully centralized receiver where the pilot and data signals
received at all APs are sent to the CPU for channel estimation and data
detection.

  * •

Level $3$ involves two stages. First, each AP estimates the channels and uses
the linear MMSE detector to detect the received signals. Then, the detected
signals are collected at the CPU for joint detection for all UEs by utilizing
the large-scale fading decoding (LSFD) method. Compared to Level $4$, only the
channel statistics are utilized at the CPU but the pilot signals are not
required to be sent to the CPU.

  * •

Level $2$ is a special case of Level $3$. The CPU performs joint detection for
all UEs by simply taking the average of the local estimates. Thus, no channel
statistics are required to be transmitted to CPU via the fronthaul.

  * •

Level $1$ is a fully distributed approach in which the data detection is
performed at the APs based on the local channel estimates. No information is
required to be transferred to the CPU.

Although the aforementioned four receivers have low complexity, the
performance is far from optimal due to the linear structure. This is because
all of the four receivers are linear receivers. In contrast, non-linear
receivers have shown great advantages in terms of BER while the main
limitation is the high computational complexity [13]. Thanks to the relatively
high computing abilities at the CPU and the large number of APs in cell-free
massive MIMO systems, we can offload parts of the computational-intensive
operations to the CPU and distribute the partial computation tasks to APs.
Next, we will propose a distributed non-linear detector for cell-free massive
MIMO systems, which takes advantage of this idea.

Figure 2: . Block diagram of the proposed distributed EP detector. The “EXT.”
blocks represents the extrinsic information computation.

##  III Proposed Distributed EP Detector

In this section, we apply the EP principle to develop a powerful distributed
MIMO detector for cell-free massive MIMO systems. After introducing the
iterative process of the algorithm, we will analyze the computational
complexity and fronthaul overhead of the proposed detector.

###  III-A Distributed Bayesian MIMO Detector

We first utilize the Bayesian inference to recover the signals ${\mathbf{x}}$
from the received signal ${\mathbf{y}}$ in the data detection stage, which is
related to the following linear model
${\mathbf{y}}={\mathbf{h}}{\mathbf{x}}+{\mathbf{n}}$. Based on Bayes’ theorem,
the posterior probability is given by

$\mathtt{P}({\mathbf{x}}|{\mathbf{y}},{\mathbf{H}})=\frac{\mathtt{P}({\mathbf{y}}|{\mathbf{x}},{\mathbf{H}})\mathtt{P}({\mathbf{x}})}{\mathtt{P}({\mathbf{y}}|{\mathbf{H}})}=\frac{\mathtt{P}({\mathbf{y}}|{\mathbf{x}},{\mathbf{H}})\mathtt{P}({\mathbf{x}})}{\int\mathtt{P}({\mathbf{y}}|{\mathbf{x}},{\mathbf{H}})\mathtt{P}({\mathbf{x}})d{\mathbf{x}}}.$
(4)

Given the posterior probability
$\mathtt{P}({\mathbf{x}}|{\mathbf{y}},{\mathbf{H}})$, the Bayesian MMSE
estimate is obtained by

$\hat{{\mathbf{x}}}=\int{\mathbf{x}}\mathtt{P}({\mathbf{x}}|{\mathbf{y}},{\mathbf{H}})d{\mathbf{x}}.$
(5)

However, the Bayesian MMSE estimator is not computationally tractable because
the marginal posterior probability in (5) involves a high-dimensional
integral. The EP algorithm, proposed in [9], provides an iterative method to
recover the transmitted ${\mathbf{x}}$ from the received signal
${\mathbf{y}}$. It is derived from the factor graph with the messages updated
and passed between different pairs of nodes assumed to follow Gaussian
distributions. As the Gaussian distribution can be fully characterized by its
mean and variance, only mean and variance are required to be calculated and
passed.

Input: Received signal ${\mathbf{y}}_{l}$, channel matrix ${\mathbf{h}}_{l}$,
noise level $\sigma^{2}$.

Output: Recovered signal ${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post},T}$.

Initialize: $\boldsymbol{\gamma}_{l}^{(0)}\leftarrow\mathbf{0}$,
$\lambda_{l}^{(0)}\leftarrow\frac{1}{E_{x}}$

for  _$t=1,\cdots,T$ _ do

Module A in APs: (1) Compute the posterior mean and variance of
${\mathbf{x}}_{\mathrm{A},l}$:
$v_{\mathrm{A},l}^{\mathrm{post}}\leftarrow\boldsymbol{\Sigma}_{l}^{t}=(\sigma^{-2}{\mathbf{h}}^{H}_{l}{\mathbf{h}}_{l}+\lambda_{l}^{(t-1)}{\mathbf{I}})$
(6)
${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{post}}\leftarrow\boldsymbol{\mu}_{l}^{t}=\boldsymbol{\Sigma}_{l}^{t}\bigg{(}\sigma^{-2}{\mathbf{h}}_{l}{\mathbf{y}}_{l}+\boldsymbol{\gamma}_{l}^{(t-1)}\bigg{)}$
(7) Module B in CPU: (2) Compute the extrinsic mean and variance of
${\mathbf{x}}_{\mathrm{A},l}$:
$v_{\mathrm{A},l}^{\mathrm{ext}}\leftarrow\bigg{(}\frac{1}{v_{\mathrm{A},l}^{\mathrm{post}}}-\lambda_{l}^{(t-1)}\bigg{)}^{-1}$
(8) ${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{ext}}\leftarrow
v_{\mathrm{A},l}^{\mathrm{ext}}\bigg{(}\frac{\boldsymbol{\mu}_{l}^{t}}{v_{\mathrm{A},l}^{\mathrm{post}}}-\boldsymbol{\gamma}_{l}^{(t-1)}\bigg{)}^{-1}$
(9) (3) MRC combining of ${\mathbf{x}}_{\mathrm{A},l}$:
$\frac{1}{v_{\mathrm{A}}^{\mathrm{ext}}}=\sum_{l=1}^{L}\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext}}}$
(10)
${\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}=v_{\mathrm{A}}^{\mathrm{ext}}\sum_{l=1}^{L}\frac{{\mathbf{x}}_{\mathrm{A},l}^{\mathrm{ext}}}{v_{\mathrm{A},l}^{\mathrm{ext}}}$
(11) (4) Compute the posterior mean and variance of
${\mathbf{x}}_{\mathrm{B}}$:
${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post}}=\mathtt{E}\\{{\mathbf{x}}|{\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}},v_{\mathrm{A}}^{\mathrm{ext}}\\}$
(12)
${\mathbf{v}}_{\mathrm{B}}^{\mathrm{post}}=\mathtt{var}\\{{\mathbf{x}}|{\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}},v_{\mathrm{A}}^{\mathrm{ext}}\\}$
(13) (5) Compute the extrinsic mean and variance
${\mathbf{x}}_{\mathrm{B},l}$:
$\frac{1}{v_{\mathrm{B},l}^{\mathrm{ext}}}\leftarrow\lambda_{l}^{(t)}=\frac{1}{\mathrm{mean}({\mathbf{v}}_{\mathrm{B}}^{\mathrm{post}})}-\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext}}}$
(14)
$\frac{{\mathbf{x}}_{\mathrm{B},l}^{\mathrm{ext}}}{v_{\mathrm{B},l}^{\mathrm{ext}}}\leftarrow\boldsymbol{\gamma}_{l}^{(t)}=\frac{{\mathbf{x}}_{\mathrm{B}}^{\mathrm{post}}}{\mathrm{mean}({\mathbf{v}}_{\mathrm{B}}^{\mathrm{post}})}-\frac{{\mathbf{x}}_{\mathrm{A},l}^{\mathrm{post}}}{v_{\mathrm{A},l}^{\mathrm{ext}}}$
(15)

Algorithm 1 Distributed EP for cell-free massive MIMO detection

Different from the conventional EP-based [13], the posterior probability in
(4) has to be rewritten in a distributed way as follows,

$\mathtt{P}({\mathbf{x}}|{\mathbf{y}},{\mathbf{H}})\varpropto\mathtt{P}({\mathbf{x}})\prod_{l=1}^{L}\mathrm{exp}(-\|{\mathbf{y}}_{l}-{\mathbf{h}}_{l}{\mathbf{x}}\|^{2}/\sigma^{2}).$
(16)

Leveraging the computational capability in the AP, we can deploy partial
calculation modules of the EP detector at the AP based on the local
information and then send the posterior mean and variance estimates to the CPU
for combining. The distributed EP-based detector is illustrated in Algorithm
1. The input of the algorithm is the received signal ${\mathbf{y}}_{l}$,
channel matrix ${\mathbf{h}}_{l}$, and noise level $\sigma^{2}$ while the
output is the recovered signal ${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post},T}$
in the $T$-th iteration. The initial parameters are
$\boldsymbol{\gamma}_{l}^{(0)}=\mathbf{0}$,
$\lambda_{l}^{(0)}=\frac{1}{E_{x}}$, where

$E_{x}=\mathtt{E}\\{\|{\mathbf{x}}\|^{2}\\}/K.$ (17)

$\boldsymbol{\gamma}_{l}^{(0)}$ and $\lambda_{l}^{(0)}$ are the initialized
extrinsic information and $E_{x}$ is the power of the transmitted symbol
$x_{k}$. Furthermore, the block diagram of the proposed distributed detector
is illustrated in Fig. 2, which is composed of module A and module B. Each
module uses the turbo principle in iterative decoding. That is, each module
passes the extrinsic messages to its next module and this process is executed
until convergence.

To better understand the distributed EP detection algorithm, we elaborate the
details for each iteration in Algorithm 1. Specifically, module A is the
linear MMSE (LMMSE) estimator performed at the APs according to the following
linear model

${\mathbf{y}}_{l}={\mathbf{h}}_{l}{\mathbf{x}}+{\mathbf{n}}_{l}.$ (18)

In the $t$-iteration of the algorithm, the explicit expression for the
posterior covariance matrix $\boldsymbol{\Sigma}_{l}^{t}$ and mean vector
$\boldsymbol{\mu}_{l}^{t}$ are given by (6) and (7), respectively. Note that
each AP only uses the local channel ${\mathbf{h}}_{l}$ to detect the
transmitted signal ${\mathbf{x}}$. For ease of notation, we omit the iteration
index $t$ for all mean and variance estimates. Then, the variance
$v_{\mathrm{A},l}^{\mathrm{post}}=\mathrm{tr}(\boldsymbol{\Sigma}_{l}^{t})/|\mathcal{D}_{l}|$
and mean estimate
${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{post}}=\boldsymbol{\mu}_{l}^{t}$ are
transferred to the CPU to compute the extrinsic information
$v_{\mathrm{A},l}^{\mathrm{ext}}$ (8) and
${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{ext}}$ (9), respectively. The extrinsic
information ${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{\mathrm{ext}}}$ can be
regarded as the AWGN observation given by

${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{ext}}={\mathbf{x}}+{\mathbf{n}}^{\mathrm{eq}}_{l},$
(19)

where
${\mathbf{n}}^{\mathrm{eq}}_{l}\sim\mathcal{N}_{{\mathbb{C}}}(0,v_{\mathrm{A},l}^{\mathrm{ext}}{\mathbf{I}})$
[11]. Therefore, the linear model in (3) is decoupled into $K$ parallel and
independent AWGN channels with equivalent noise
$v_{\mathrm{A},l}^{\mathrm{ext}}$. Subsequently, the CPU collects all
extrinsic means $\\{{\mathbf{x}}_{\mathrm{A},l}^{\mathrm{ext}}\\}_{l=1}^{L}$
and variances $\\{v_{\mathrm{A},l}^{\mathrm{ext}}\\}_{l=1}^{L}$ and performs
MRC. The MRC expressions (10) and (11) are obtained by maximizing the post-
combination signal-to-noise ratio (SNR) of the final AWGN observation
${\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}$ at the CPU, given by

${\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}={\mathbf{x}}+{\mathbf{n}}^{\mathrm{eq}},$
(20)

where
${\mathbf{n}}_{\mathrm{A}}^{\mathrm{ext}}\sim\mathcal{N}_{{\mathbb{C}}}(0,v_{\mathrm{A}}^{\mathrm{ext}}{\mathbf{I}})$.
The CPU uses the posterior mean estimator to detect the signal ${\mathbf{x}}$
from the equivalent AWGN model (20). Then, the posterior mean and variance are
computed by the posterior MMSE estimator for the equivalent AWGN model in (12)
and (13). As the transmitted symbol is assumed drawn from the $M$-QAM set
$\mathcal{S}=\\{s_{1},s_{2},\ldots,s_{M}\\}$, the corresponding expressions
for each element in (12) and (13) are given by

$[x_{\mathrm{B}}^{\mathrm{post}}]_{k}=\frac{\sum_{s_{i}\in\mathcal{S}}s_{i}\mathcal{N}_{{\mathbb{C}}}(s_{i};[x_{\mathrm{A}}^{\mathrm{ext}}]_{k},v_{\mathrm{A}}^{\mathrm{ext}})p(s_{i})}{\sum_{s_{i}\in\mathcal{S}}\mathcal{N}_{{\mathbb{C}}}(s_{i};[x_{\mathrm{A}}^{\mathrm{ext}}]_{k},v_{\mathrm{A}}^{\mathrm{ext}})p(s_{i})}$
(21)

$[v_{\mathrm{B}}^{\mathrm{post}}]_{k}=\frac{\sum_{s_{i}\in\mathcal{S}}|s_{i}|^{2}\mathcal{N}_{{\mathbb{C}}}(s_{i};[x_{\mathrm{A}}^{\mathrm{ext}}]_{k},v_{\mathrm{A}}^{\mathrm{ext}})p(s_{i})}{\sum_{s_{i}\in\mathcal{S}}\mathcal{N}_{{\mathbb{C}}}(s_{i};[x_{\mathrm{A}}^{\mathrm{ext}}]_{k},v_{\mathrm{A}}^{\mathrm{ext}})p(s_{i})}-|[x_{\mathrm{B}}^{\mathrm{post}}]_{k}|^{2}.$
(22)

where $[x_{\mathrm{B}}^{\mathrm{post}}]_{k}$,
$[v_{\mathrm{B}}^{\mathrm{post}}]_{k}$, and
$[x_{\mathrm{A}}^{\mathrm{ext}}]_{k}$ are the $k$-th element in
${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post}}$,
${\mathbf{v}}_{\mathrm{B}}^{\mathrm{post}}$, and
$[{\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}]_{k}$, respectively. The posterior
mean and variance ${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post}}$ and
${\mathbf{v}}_{\mathrm{B}}^{\mathrm{post}}$ are then utilized to compute the
extrinsic information $\lambda_{l}^{(t)}$ and $\boldsymbol{\gamma}_{l}^{(t)}$
for each AP in (14) and (15), where the function $\mathrm{mean}(\cdot)$ is
used to compute the mean. Finally, the extrinsic information
$\lambda_{l}^{(t)}$ and $\boldsymbol{\gamma}_{l}^{(t)}$ are transferred to
each AP in the next iteration. The whole procedure is executed iteratively
until terminated by a certain stopping criterion or a maximum number of
iterations.

TABLE I: . Complexity of different detectors Detectors | Distributed EP | Level 4 | Level 3&2 | Level 1
---|---|---|---|---
AP | $O(T|\mathcal{D}_{l}|N^{2})$ | $0$ | $O(|\mathcal{D}_{l}|N^{2})$ | $O(|\mathcal{D}_{l}|N^{2})$
CPU | $O(T|\mathcal{D}_{l}|^{2})$ | $O(LN)^{3}$ | $O(|\mathcal{D}_{l}|)$ | $0$

###  III-B Computational Complexity and Fronthaul Overhead

In the following, we provide the complexity analysis for different detectors
in Tables I and II. For the proposed distributed EP detector, the
computational complexity at each AP is dominated by the LMMSE estimator for
estimating the signal ${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{post}}$, which is
$O(|\mathcal{D}_{l}|^{3})$ because of the matrix inversion is required in (6)
while the computational complexity at the CPU is $O(|\mathcal{D}_{l}|^{2})$ in
each iteration. Furthermore, if the number of antennas $N$ is less than
$|\mathcal{D}_{l}|$, we can use the matrix inversion lemma to carry out the
matrix inversion in (6) as follows

$\displaystyle(\sigma^{-2}{\mathbf{h}}_{l}^{H}{\mathbf{h}}_{l}+{\mathbf{D}})^{-1}={\mathbf{D}}^{-1}$
(23)
$\displaystyle-\sigma^{-2}{\mathbf{D}}^{-1}{\mathbf{h}}_{l}^{H}({\mathbf{I}}+\sigma^{-2}{\mathbf{h}}_{l}{\mathbf{D}}^{-1}{\mathbf{h}}_{l}^{H})^{-1}{\mathbf{h}}_{l}{\mathbf{D}}^{-1},$

where ${\mathbf{D}}=\lambda_{l}^{(t-1)}{\mathbf{I}}$ and the computational
complexity is reduced to $O(|\mathcal{D}_{l}|N^{2})$. Therefore, the overall
computational complexity at each AP is $O(T|\mathcal{D}_{l}|N^{2})$ while the
overall computational complexity at the CPU is $O(T|\mathcal{D}_{l}|^{2})$ for
$T$ iterations. As can be observed in Table I, the distritbuted EP detector
mainly increases the computational complexity at the CPU when compared with
other distributed linear receivers (Levels 1-3).

TABLE II: . Fronthaul overhead Detectors | Coherence block | Statistical parameters
---|---|---
Distributed EP | $\sum\limits_{l=1}^{L}(\tau_{c}-\tau_{p})2T(|\mathcal{D}_{l}|+1)$ | $0$
Level 4 | $\tau_{c}NL$ | $KLN^{2}/2$
Level 3 | $(\tau_{c}-\tau_{p})KL$ | $KL+(L^{2}K^{2}+KL)/2$
Level 2 | $(\tau_{c}-\tau_{p})KL$ | $0$
Level 1 | $0$ | $0$

We compare the number of complex scalars that need to be transmitted from the
APs to the CPU via the fronthauls in Table II. We assume that $\tau_{c}$ and
$\tau_{p}$ are the coherence time and pilot length, respectively. The results
for the 4 level detectors are from [8]. With the proposed detector,
$\sum\limits_{l=1}^{L}(\tau_{c}-\tau_{p})2T(|\mathcal{D}_{l}|+1)$ scalars need
to be passed from the APs to the CPU and no statistical parameters are
required to be passed, where $T$ denotes the total iterations. The fronthaul
overhead of the proposed distributed EP detector is similar to Level 4 and the
detailed comparison is determined by the value of the system parameters.

##  IV State Evolution Analysis

In this section, we provide an analytical framework to predict the asymptotic
performance of the distributed EP in the large system limit. We consider
$L,K,N\rightarrow\infty$ and fix

$\alpha_{l}=\frac{K}{N},\alpha=\frac{K}{LN}.$ (24)

Then, we have following proposition.

######  Proposition 1.

In the large-system limit, the asymptotic behavior (such as MSE and BER) of
Algorithm 1 can be described by the following equations:

$\displaystyle v_{\mathrm{A},l}^{\mathrm{ext},t}$
$\displaystyle=\frac{\alpha_{l}\sigma^{2}+(\alpha_{l}-1)\lambda_{l}^{(t-1)}}{2}$
(25a)
$\displaystyle+\frac{\sqrt{\bigg{(}\alpha_{l}\sigma^{2}+(\alpha_{l}-1)\lambda_{l}^{(t-1)}\bigg{)}^{2}+4\alpha_{l}\sigma^{2}\lambda_{l}^{(t-1)}}}{2}$
$\displaystyle
v_{\mathrm{A}}^{\mathrm{ext},t}=\bigg{(}\sum_{l=1}^{L}\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext},t}}\bigg{)}^{-1}$
(25b)
$\displaystyle\lambda_{l}^{(t)}=\frac{1}{\mathrm{MSE}(v_{\mathrm{A}}^{\mathrm{ext},t})}-\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext},t}}.$
(25c)

$\blacksquare$

The function $\mathrm{MSE}(\cdot)$ is given by

$\mathrm{MSE}(v_{\mathrm{A},l}^{\mathrm{ext},t})\triangleq\mathtt{E}\bigg{\\{}|x-\mathtt{E}\\{x|x_{\mathrm{A},l}^{\mathrm{ext},t},v_{\mathrm{A},l}^{\mathrm{ext},t}|^{2}\bigg{\\}}$
(26)

and the expectation is with respect to $x$. The state equations can be proved
using the method in [12], and we give an intuitive explanation here, which is
linked to Algorithm 1. By substituting (6) into (8), we have following
expression as

$v_{\mathrm{A},l}^{\mathrm{ext},t}=\frac{1}{\frac{1}{K}\mathrm{tr}(\sigma^{-2}{\mathbf{h}}^{H}_{l}{\mathbf{h}}_{l}+\lambda_{l}^{(t-1)}{\mathbf{I}})}-\lambda_{l}^{(t-1)}.$
(27)

Notably, $v_{\mathrm{A},l}^{\mathrm{ext}}$ can be characterized by the
eigenvalues of ${\mathbf{h}}^{H}_{l}{\mathbf{h}}_{l}$, which converge to a
deterministic distribution in the large system limit. Because of the limited
space, we only show the final result of the Proposition 1 and omit the derived
process. By adopting the result of the $\mathcal{R}$-transform of the average
empirical eigenvalue distribution, we have an asymptotic expression for
$v_{\mathrm{A},l}^{\mathrm{ext},t}$ in (25). Then, the asymptotic expressions
for $v_{\mathrm{A}}^{\mathrm{ext},t}$ and $\lambda_{l}^{(t)}$ can be derived
from (10) and (14). Finally, the asymptotic MSE can be interpreted as the MSE
of the decoupled scalar AWGN channels (20) and is related to distribution of
transmitted signal ${\mathbf{x}}$. Next, we will give a specific example for
Proposition 1.

Example 1: If the data symbol is drawn from a quadrature phase-shift keying
(QPSK) constellation, the $\mathrm{MSE}$ is expressed by

$\mathrm{MSE}=1-\int\mathrm{D}z\tanh(\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext},t}}+\sqrt{\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext},t}}}z).$
(28)

Furthermore, the final BER w.r.t. ${\mathbf{x}}$ can also be evaluated through
the equivalent AWGN channel (20) with an equivalent
$\mathrm{SNR}=1/v_{\mathrm{A},l}^{\mathrm{ext},T}$ and is given by [14]

$\mathrm{BER}=2Q(\sqrt{\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext},T}}})-[Q(\sqrt{\frac{1}{v_{\mathrm{A},l}^{\mathrm{ext},T}}})]^{2},$
(29)

where $Q(x)=\int_{x}^{\infty}\mathrm{D}z$ is the $Q$-function. In fact, the
MSE and BER are determined on the basis of the knowledge of the AWGN channel
(20) with $\mathrm{SNR}=1/v_{\mathrm{A},l}^{\mathrm{ext},T}$, which is known
as the decoupling principle. Thus, if the data symbol is drawn from other
$M$-QAM constellations, then the corresponding BER can be easily obtained
using the closed-form BER expression in [15]. Note that the state evolution
equations are accurate when each AP serves all UEs in the network and the
channel matrix ${\mathbf{H}}$ is a Rayleigh fading channel. This assumption
can be achieved by perfect power control to compensate the large-scaled fading
effect.

##  V Simulation Results

In this section, we provide sample simulation results to demonstrate the
performance of the proposed distributed EP detector for cell-free massive
MIMO. Consider $L=16$, $N=8$, $K=16$, and the other simulation parameters are
the same as those in [8]. The SNR is defined as $\mathrm{SNR}=1/\sigma^{2}$.

Figure 3: . BER performance comparisons of different detectors in the original
cell-free massive MIMO system.

###  V-A Original Cell-Free Massive MIMO

In this subsection, we consider the BER of the proposed detector in the
original cell-free massive MIMO where each AP serves all UEs. We assume that
perfect CSI can be obtained at each AP. Fig. 3 compares the achievable BER of
the proposed distributed EP detector with other detectors investigated in [8].
The results are obtained by the Monte Carlo simulation with $10,000$
independent channel realizations. We denote “deEP” as the distributed EP
detector. It can be observed that the proposed distributed EP detector is
comparable with the Level 1 detector with only one EP iteration and
outperforms the centralized Level 4 detector with $T=5$ iterations.

We next the accuracy of the analytical framework in Fig. 4 with different
modulation schemes. As shown in the figure, the BERs of the proposed detector
match well with the derived analytical results, which demonstrate the accuracy
of the analytical framework. Therefore, instead of performing time-consuming
Monte Carlo simulations to obtain the corresponding performance metrics, we
can predict the theoretical behavior by state equations. Furthermore, the
analytical framework can be further utilized to optimize the system design.

Figure 4: . BER performance comparisons of the analytical and simulation
results of the distributed EP detector with different modulation schemes.

###  V-B Scalable Cell-Free Massive MIMO

As the _user-centric_ approach is more attractive for cell-free massive MIMO,
we then investigate the distributed EP detector with DCC. The accessing UE
first appoints a master AP according to the large-scale fading factor and
assigns a pilot to the appointed AP. Then, other neighboring APs determine
whether they serve the accessing UE according to the assigned pilot. Finally,
the cluster for the $k$-th UE is constructed. Fig. 5 shows that the
distributed EP detector outperforms both the centralized and distributed MMSE
detectors. Furthermore, the performance loss is acceptable when compared to
original cell-free massive MIMO and the computational complexity is
significantly decreased from $O(KN^{2})$ to $O(|\mathcal{D}_{l}|N^{2})$.

Figure 5: . BER performance comparisons of different detectors in cell-free
massive MIMO with dynamic cluster cooperation.

##  VI Conclusion

In this paper, we proposed a distributed EP detector for cell-free massive
MIMO. It is shown that such detector can achieve better performance than other
linear receivers for both original and scalable cell-free massive MIMO
networks. Compared to other distributed detectors, it can achieve a better BER
performance and mainly increase the computational overhead at the CPU. An
analytical framework was also provided to describe the asymptotic performance
of the proposed detector in a large system setting. Simulation results have
also been provided and demonstrated that the proposed method outperforms the
existing distributed detectors for cell-free massive MIMO in terms of BER
performance.

## References

  * [1] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” IEEE Trans. Wireless Commun., vol. 9, no. 11, pp. 3590–3600, Nov. 2010.
  * [2] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, ”Cell-free massive MIMO versus small cells,” IEEE Trans. Wireless Commun., vol. 16, no. 3, pp. 1834–1850, Mar. 2017.
  * [3] J. Zhang, S. Chen, Y. Lin, J. Zheng, B. Ai, and L. Hanzo, “Cell-free massive MIMO: A new next-generation paradigm,” IEEE Access, vol. 7, pp. 99878–99888, Sep. 2019.
  * [4] K. B. Letaief, W. Chen, Y. Shi, J. Zhang, and Y.-J.-A. Zhang, “The roadmap to 6G: AI empowered wireless networks,” IEEE Commun. Mag., vol. 57, no. 8, pp. 84-90, Aug. 2019.
  * [5] E. Björnson, N. Jaldén, M. Bengtsson, and B. Ottersten, “Optimality properties, distributed strategies, and measurement-based evaluation of coordinated multicell OFDMA transmission,” IEEE Trans. Signal Process., vol. 59, no. 12, pp. 6086–6101, Dec. 2011.
  * [6] E. Björnson and L. Sanguinetti, “Scalable cell-free massive MIMO systems,” IEEE Trans. Commun., vol. 68, no. 7, pp. 4247–4261, Jul. 2020.
  * [7] E. Nayebi, A. Ashikhmin, T. L. Marzetta, and B. D. Rao, “Performance of cell-free massive MIMO systems with MMSE and LSFD receivers,” in Proc. 50th Asilomar Conf. Signals, Syst. Comput., Nov. 2016, pp. 203–207.
  * [8] E. Björnson and L. Sanguinetti, “Making cell-free massive MIMO competitive with MMSE processing and centralized implementation,” IEEE Trans. Wireless Commun., vol. 19, no. 1, pp. 77–90, Jan. 2020.
  * [9] T. P. Minka, “A family of algorithms for approximate Bayesian Inference,” Ph.D. dissertation, Dept. Elect. Eng. Comput. Sci., MIT, Cambridge, MA, USA, 2001.
  * [10] C. Jeon, K. Li, J. R. Cavallaro, and C. Studer, “Decentralized equalization with feedforward architectures for massive MU-MIMO,” IEEE Trans. Signal Process., vol. 67, no. 17, pp. 4418–4432, Sep. 2019.
  * [11] H. Wang, A. Kosasih, C. Wen, S. Jin, and W. Hardjawana, “Expectation propagation detector for extra-large scale massive MIMO,” IEEE Trans. Wireless Commun., vol. 19, no. 3, pp. 2036–2051, Mar. 2020.
  * [12] K. Takeuchi, “Rigorous dynamics of expectation-propagation-based signal recovery from unitarily invariant measurements,” IEEE Trans. Inf. Theory., vol. 66, no. 1, 368–386, Oct. 2019.
  * [13] J. Céspedes, P. M. Olmos, M. Sánchez-Fernández, and F. Perez-Cruz, “Expectation propagation detection for high-order high-dimensional MIMO systems,” IEEE Trans. Commun., vol. 62, no. 8, pp. 2840–2849, Aug. 2014.
  * [14] H. He, C.-K. Wen, and S. Jin, “Bayesian optimal data detector for hybrid mmWave MIMO-OFDM systems with low-resolution ADCs,” IEEE J. Sel. Topics Signal Process., vol. 12, no. 3, pp. 469–483, Jun. 2018.
  * [15] J. G. Proakis, _Digital Communications_. Boston, USA: McGraw-Hill Companies, 2007.