# Single Inclusive Hadron Production in DIS at Small $x$: Next to Leading
Order Corrections

Filip Bergabo<EMAIL_ADDRESS>Department of Natural Sciences,
Baruch College, CUNY, 17 Lexington Avenue, New York, NY 10010, USA  City
University of New York Graduate Center, 365 Fifth Avenue, New York, NY 10016,
USA  Jamal Jalilian-Marian<EMAIL_ADDRESS>Department
of Natural Sciences, Baruch College, CUNY, 17 Lexington Avenue, New York, NY
10010, USA  City University of New York Graduate Center, 365 Fifth Avenue, New
York, NY 10016, USA

###### Abstract

We calculate the one-loop corrections to single inclusive hadron production in
Deep Inelastic Scattering (DIS) at small $x$ in the forward rapidity region
using the Color Glass Condensate formalism. We show that the divergent parts
of the next to leading order (NLO) corrections either cancel among each other
or lead to $x$ (rapidity) evolution of the leading order (LO) dipole cross
section according to the JIMWLK evolution equation and DGLAP evolution of the
parton-hadron fragmentation function. The remaining finite parts constitute
the NLO ($\alpha_{s}$) corrections to the LO single inclusive hadron
production cross section in DIS at small $x$.

##  I Introduction

Gluon saturation Gribov et al. (1983); Mueller and Qiu (1986) at small $x$ as
encoded in the Color Glass Condensate (CGC) formalism Iancu and Venugopalan
(2003); Iancu et al. (2002); Jalilian-Marian and Kovchegov (2006); Weigert
(2005); Morreale and Salazar (2021) has been the subject of intense
theoretical studies and experimental searches. Theoretical work based on
leading order (LO) or leading log (LL) approximations to gluon saturation have
successfully described structure functions, suppression of the single
inclusive hadron transverse momentum spectrum and disappearance of the away
side peak in dihadron angular correlations in high energy
proton(deuteron)-gold/lead collisions at RHIC and the LHC Kovner and Wiedemann
(2001); Jalilian-Marian and Kovchegov (2004); Dumitru et al. (2006); Jalilian-
Marian (2006); Marquet (2007); Albacete and Marquet (2010); Stasto et al.
(2012); Lappi and Mantysaari (2013); Jalilian-Marian and Rezaeian (2012a, b);
Zheng et al. (2014); Stasto et al. (2018); Albacete et al. (2019); Mäntysaari
et al. (2020); Hatta et al. (2021a); Jia et al. (2020); Gelis and Jalilian-
Marian (2002); Dominguez et al. (2011); Metz and Zhou (2011); Dominguez et al.
(2012); Iancu and Laidet (2013); Altinoluk et al. (2016a); Hatta et al.
(2016); Dumitru et al. (2015); Kotko et al. (2015); Marquet et al. (2016); van
Hameren et al. (2016); Marquet et al. (2018); Dumitru et al. (2019); Dumitru
and Jalilian-Marian (2002a, b); Mäntysaari et al. (2019); Salazar and Schenke
(2019); Boussarie et al. (2021); Ayala et al. (1996); Jalilian-Marian et al.
(1997a); Kotko et al. (2017); Hagiwara et al. (2017); Henley and Jalilian-
Marian (2006); Klein and Mäntysaari (2019); Hatta et al. (2021b); Kolbé et al.
(2021); Gelis and Jalilian-Marian (2002); Altinoluk et al. (2019); Boussarie
et al. (2019, 2016, 2014); Dumitru and Jalilian-Marian (2010). Nevertheless
firmly establishing gluon saturation as the QCD dynamics responsible for these
experimental observations requires more precise theoretical calculations. The
ongoing work on improving the accuracy of leading order CGC calculations can
be broadly put into three categories; higher order in $\alpha_{s}$ corrections
to leading order results Fadin and Lipatov (1998); Chirilli et al. (2012a, b);
Balitsky and Chirilli (2013a, b); Grabovsky (2013); Caron-Huot (2015); Kovner
et al. (2014); Lublinsky and Mulian (2017); Caron-Huot and Herranen (2018);
Boussarie et al. (2018); Beuf et al. (2022a, b, 2021); Mäntysaari and Penttala
(2021); Mäntysaari and Penttala (2022a, b); Lappi et al. (2022); Iancu and
Mulian (2021); Roy and Venugopalan (2020); Hatta et al. (2022); Iancu et al.
(2022); Taels et al. (2022); Caucal et al. (2021); Bergabo and Jalilian-Marian
(2022a, b), sub-eikonal corrections which aim to relax the infinite energy
assumption inherent to eikonal approximation Kovchegov et al. (2017a);
Cougoulic and Kovchegov (2019); Kovchegov and Sievert (2019); Kovchegov et al.
(2017b, c, d); Kovchegov et al. (2016); Agostini et al. (2019a, b); Altinoluk
and Dumitru (2016); Altinoluk et al. (2016b); Altinoluk et al. (2014), and
inclusion of intermediate/large $x$ dynamics into CGC in order to generalize
CGC to include DGLAP evolution and collinear factorization and high $p_{t}$
physics Jalilian-Marian (2021, 2020, 2019, 2017); Hentschinski et al. (2018,
2016); Gituliar et al. (2016); Balitsky and Tarasov (2016, 2015). Here we will
focus on next to leading order corrections to single inclusive hadron
production in Deep Inelastic Scattering (DIS) at small $x$ in the forward
rapidity region Marquet et al. (2009) (virtual photon going direction) for the
case when the virtual photon is longitudinal. We note that leading order
results for single inclusive hadron production in DIS in the midrapidity
region were obtained in Kovchegov and Tuchin (2002).

The ideal environment in which to investigate gluon saturation and CGC is DIS
experiments at high energy as the incoming virtual photon does not interact
strongly. Single inclusive hadron production in DIS (SIDIS) at small $x$ is
one of the most attractive channels for gluon saturation studies as it is not
sensitive to Sudakov effects which can obscure saturation dynamics in dihadron
production and angular correlations. Furthermore, it is more discriminatory
than the total cross section (structure functions) so that it contains more
information about the QCD dynamics of the target. While there exists leading
order calculations of single inclusive hadron production in DIS at small $x$
in the CGC framework Marquet et al. (2009); Kovchegov and Tuchin (2002) it is
highly desirable and in fact urgently needed to perform a next to leading
order calculation which can then be used for quantitative studies of the
transverse momentum spectra of produced hadrons in DIS with proton and nuclear
targets at the proposed Electron Ion Collider (EIC).

Here we calculate the next to leading order corrections to single inclusive
hadron production in DIS at small $x$ in the forward rapidity region using the
Color Glass Condensate formalism. To do so we use our recent results for next
to leading order corrections to dihadron production Bergabo and Jalilian-
Marian (2022b) in DIS and integrate out one of the final state partons. As
expected we encounter various divergences which appear when we integrate over
the phase space of the final state parton. We show that UV and soft
divergences cancel among each other while the collinear divergences associated
with radiation of a massless parton are absorbed into the parton-hadron
fragmentation function. We show that all quadrupole terms appearing in the
intermediate steps of the calculation cancel among various terms and one is
left with dipoles (and squared dipoles) only. The rapidity divergences arising
from integrating over longitudinal phase space of the final state parton are
absorbed into evolution of the dipoles describing the target dynamics and lead
to JIMWLK evolution of the leading order cross section. The remaining terms
are finite and constitute the $O(\alpha_{s})$ corrections to leading order
single inclusive hadron production in DIS at small $x$.

##  II Leading Order Cross Section

To get the leading order single inclusive hadron production in DIS at small
$x$ we start with the quark antiquark production cross section in DIS given by

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to
q\bar{q}X}}{\differential^{2}\mathbf{p}\,\differential^{2}\mathbf{q}\,\differential
y_{1}\,\differential y_{2}}$ $\displaystyle=$
$\displaystyle\frac{e^{2}Q^{2}(z_{1}z_{2})^{2}N_{c}}{(2\pi)^{7}}\delta(1-z_{1}-z_{2})\int\differential^{8}\mathbf{x}\left[S_{122^{\prime}1^{\prime}}-S_{12}-S_{1^{\prime}2^{\prime}}+1\right]$
(1) $\displaystyle
e^{i\mathbf{p}\cdot(\mathbf{x}_{1}^{\prime}-\mathbf{x}_{1})}e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}\bigg{[}4z_{1}z_{2}K_{0}(|\mathbf{x}_{12}|Q_{1})K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{1})+$
$\displaystyle(z_{1}^{2}+z_{2}^{2})\,\frac{\mathbf{x}_{12}\cdot\mathbf{x}_{1^{\prime}2^{\prime}}}{|\mathbf{x}_{12}||\mathbf{x}_{1^{\prime}2^{\prime}}|}\,K_{1}(|\mathbf{x}_{12}|Q_{1})K_{1}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{1})\bigg{]}.$

where ($\mathbf{p},y_{1}$) and ($\mathbf{q},y_{2}$) are the transverse
momentum and rapidity of the produced quark and antiquark, respectively, and
$Q^{2}$ is the virtuality of the incoming photon. We have made the following
definitions and short hand notations,

$\displaystyle
Q_{i}=Q\sqrt{z_{i}(1-z_{i})},\,\,\,\,\,\,\mathbf{x}_{ij}=\mathbf{x}_{i}-\mathbf{x}_{j},\,\,\,\,\,\,\differential^{8}\mathbf{x}=\differential^{2}\mathbf{x}_{1}\,\differential^{2}\mathbf{x}_{2}\,\differential^{2}\mathbf{x}_{1^{\prime}}\,\differential^{2}\mathbf{x}_{2^{\prime}}.$
(2)

We have also defined
$z_{1}\equiv\frac{p^{+}}{l^{+}}\,\,,\,\,z_{2}\equiv\frac{q^{+}}{l^{+}}$ as the
momentum fractions carried by the final state quark and antiquark relative to
the photon’s longitudinal momentum $l^{+}$. In terms of these momentum
fractions the rapidity is related via $\differential y_{i}=\frac{\differential
z_{i}}{z_{i}}$. All the dynamics of the strong interactions and gluon
saturation are contained in the dipoles $S_{ij}$ and quadrupoles $S_{ijkl}$,
normalized correlation functions of two and four Wilson lines

$\displaystyle S_{ij}=\frac{1}{N_{c}}{\rm tr}\left\langle
V_{i}V_{j}^{\dagger}\right\rangle,\,\,\,\,\,\,\,\,\,S_{ijkl}=\frac{1}{N_{c}}{\rm
tr}\left\langle V_{i}V_{j}^{\dagger}V_{k}V_{l}^{\dagger}\right\rangle,$ (3)

where the index $i$ refers to the transverse coordinate $\mathbf{x}_{i}$ and
the following notation is used for Wilson lines,

$\displaystyle V_{i}$ $\displaystyle=\hat{P}\exp\left(ig\int\differential
x^{+}A^{-}(x^{+},\mathbf{x}_{i})\right).$ (4)

The Wilson lines efficiently resum the multiple scatterings of the quark and
antiquark from the target hadron or nucleus. The angle brackets in Eq. 3
signify color averaging 111Throughout the paper we assume that these dipoles
and quadrupoles are real, nevertheless both can have imaginary parts which
however do not contribute here.. It is important to keep in mind that as this
is a classical result the cross section has no non-trivial $x$ (or
rapidity/energy) dependence. It is also easy to check that if one integrates
over the phase space of the quark and antiquark one recovers the standard
expressions for the virtual photon-target total cross section at small $x$.

Integrating over the quark’s momentum then sets $z_{1}=1-z_{2}$ and
$\mathbf{x}^{\prime}_{1}=\mathbf{x}_{1}$ and gives

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to\bar{q}X}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}$ $\displaystyle=$
$\displaystyle\frac{e^{2}Q^{2}z_{2}^{2}(1-z_{2})N_{c}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}\left[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1\right]$
(5) $\displaystyle
e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}\bigg{[}4z_{2}(1-z_{2})K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})+$
$\displaystyle\left[z_{2}^{2}+(1-z_{2})^{2}\right]\,\frac{\mathbf{x}_{12}\cdot\mathbf{x}_{12^{\prime}}}{|\mathbf{x}_{12}||\mathbf{x}_{12^{\prime}}|}\,K_{1}(|\mathbf{x}_{12}|Q_{2})K_{1}(|\mathbf{x}_{12^{\prime}}|Q_{2})\bigg{]}$

where the first (second) term inside the big square bracket corresponds to
contribution of longitudinal (transverse) photons. To get the full single
inclusive production cross section one must also consider the case when one
integrates out the antiquark. It can however be shown that the two results are
identical so that we will only integrate out the quark and multiple our final
results by a factor of $2$. This can also be shown to be true when we
calculate the next to leading order corrections. Therefore we will consider
only the case when the quark is integrated out. Furthermore and as before we
will consider only the case of longitudinal photons in this paper.

##  III One-loop corrections

Figure 1: The real corrections
$i\mathcal{A}_{1}^{a},...,i\mathcal{A}_{4}^{a}$. The arrows on Fermion lines
indicate Fermion number flow, all momenta flow to the right. The thick solid
line indicates interaction with the target.

Figure 2: The ten virtual NLO diagrams
$i\mathcal{A}_{5},...,i\mathcal{A}_{14}$. All momenta flow to the right,
except for gluon momenta.

To calculate the next to leading order corrections to single inclusive hadron
production we start with our next to leading order results for dihadron
production computed in Bergabo and Jalilian-Marian (2022b). The real
corrections labeled $\differential\sigma_{i\times j}$ come from squaring the
diagrams in Fig. 1 (these were first calculated in Ayala et al. (2016, 2017))
and the virtual corrections $\differential\sigma_{i}$ from multiplying the
diagrams in Fig. 2 with the leading order amplitude. These must then be
multiplied by their corresponding phase space differentials
$\differential\Phi^{(n)}$. The explicit details are shown in Eq. 6 where we
have also defined $i\mathcal{M}_{i}^{a}$ via
$i\mathcal{A}_{i}^{a}=2\pi\delta(l^{+}-p^{+}-q^{+}-k^{+})i\mathcal{M}_{i}^{a}$
and $i\mathcal{M}_{i}$ via
$i\mathcal{A}_{i}=2\pi\delta(l^{+}-p^{+}-q^{+})i\mathcal{M}_{i}$. We take the
flux factor $\mathcal{F}$ to be $2l^{+}$. This gives the NLO corrections to
quark antiquark production, and so we then integrate out the quark (i.e.:
perform the integral over $\mathbf{p}$ and $y_{1}$ in each expression) as
before to obtain the single inclusive results in Eq. 7 \- 22. We note that
obtaining the complete result for single inclusive hadron production starting
from our original expressions for one loop corrections to dihadron production
would require going back and integrating out any two of the three partons in
the final state (in real corrections). As our main goal here is to demonstrate
factorization of the cross section and cancellation/absorption of all
divergences it is enough to focus on the case when the radiated gluon was
first integrated out to get the NLO corrections to dihadron production.

$\displaystyle\differential\sigma^{L}_{\text{NLO}}=$
$\displaystyle\sum_{i,j=1}^{4}\differential\sigma^{L}_{i\times
j}+2\,\text{Re}\sum_{i=5}^{14}\differential\sigma_{i}^{L},$
$\displaystyle\differential\sigma^{L}_{i\times j}=$
$\displaystyle\frac{1}{\mathcal{F}}\int_{z,\mathbf{k}}\left[(i\mathcal{M}_{i}^{a})(i\mathcal{M}_{j}^{a})^{*L}\right]\differential\Phi^{(3)},\,\,\,\,\,\,\differential\sigma^{L}_{i}=\frac{1}{\mathcal{F}}\left[(i\mathcal{M}_{i})(i\mathcal{M})^{*,L}\right]\differential\Phi^{(2)},$
$\displaystyle\differential\Phi^{(3)}=$ $\displaystyle
2l^{+}\frac{\differential^{2}\mathbf{p}\,\differential^{2}\mathbf{q}\,\differential^{2}\mathbf{k}\,\differential
y_{1}\,\differential y_{2}\,\differential
z}{(2\pi)^{8}(4l^{+})^{2}z}\delta(1-z_{1}-z_{2}-z),\,\,\,\,\,\,\,\differential\Phi^{(2)}=2l^{+}\frac{\differential^{2}\mathbf{p}\,\differential^{2}\mathbf{q}\,\differential
y_{1}\,\differential y_{2}}{2(2\pi)^{5}(2l^{+})^{2}}\delta(1-z_{1}-z_{2}).$
(6)

$\displaystyle\frac{\differential\sigma^{L}_{2\times
2}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}}{(2\pi)^{8}z_{2}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}(1-z_{2}-z)^{2}(z+z_{2})^{2}\left[z_{2}^{2}+(z+z_{2})^{2}\right]\int\differential^{8}\mathbf{x}\,K_{0}(|\mathbf{x}_{12}|Q_{1})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{1})\Delta^{(3)}_{22^{\prime}}$
$\displaystyle[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}e^{i\frac{z}{z_{2}}\mathbf{q}\cdot(\mathbf{x}^{\prime}_{2}-\mathbf{x}_{2})}$
(7) $\displaystyle\frac{\differential\sigma^{L}_{1\times
2}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}(1-z_{2}-z)(z_{2}+z)\left[z_{2}(1-z_{2})+(1-z_{2}-z)(z_{2}+z)\right]\int\differential^{8}\mathbf{x}\,K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{2})$
$\displaystyle
K_{0}(|\mathbf{x}_{12}|Q_{1})\Delta^{(3)}_{12}[S_{12}S_{1^{\prime}2^{\prime}}-S_{12}-S_{1^{\prime}2^{\prime}}+1]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}e^{-i\frac{z}{z_{2}}\mathbf{q}\cdot(\mathbf{x}_{2}-\mathbf{x}_{3})}$
(8) $\displaystyle\frac{\differential\sigma^{L}_{3\times
3}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{3}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}\left[(1-z_{2}-z)^{2}+(1-z_{2})^{2}\right]\int\differential^{8}\mathbf{x}\,K_{0}(QX)K_{0}(QX_{1}^{\prime})\frac{1}{\mathbf{x}_{31}^{2}}$
$\displaystyle[S_{22^{\prime}}-S_{13}S_{23}-S_{13}S_{2^{\prime}3}+1]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(9) $\displaystyle\frac{\differential\sigma^{L}_{4\times
4}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}(1-z_{2}-z)^{2}\left[z_{2}^{2}+(z+z_{2})^{2}\right]\int\differential^{8}\mathbf{x}\,K_{0}(QX)K_{0}(QX_{1}^{\prime})\Delta^{(3)}_{22^{\prime}}$
$\displaystyle[S_{22^{\prime}}-S_{13}S_{23}-S_{13}S_{2^{\prime}3}+1]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(10) $\displaystyle\frac{\differential\sigma^{L}_{3\times
4}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{2}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}(1-z_{2}-z)\left[z_{2}(1-z_{2})+(1-z_{2}-z)(z_{2}+z)\right]\int\differential^{8}\mathbf{x}\,K_{0}(QX)K_{0}(QX_{1}^{\prime})\Delta^{(3)}_{12^{\prime}}$
$\displaystyle[S_{22^{\prime}}-S_{13}S_{23}-S_{13}S_{2^{\prime}3}+1]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(11) $\displaystyle\frac{\differential\sigma^{L}_{2\times
3}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}(1-z_{2}-z)(z_{2}+z)\left[(1-z_{2}-z)(z_{2}+z)+z_{2}(1-z_{2})\right]\int\differential^{8}\mathbf{x}K_{0}(|\mathbf{x}_{12}|Q_{1})K_{0}(QX_{1}^{\prime})\Delta^{(3)}_{21}$
$\displaystyle[S_{23}S_{2^{\prime}3}-S_{13}S_{2^{\prime}3}-S_{12}+1]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}e^{i\frac{z}{z_{2}}\mathbf{q}\cdot(\mathbf{x}_{3}-\mathbf{x}_{2})}$
(12) $\displaystyle\frac{\differential\sigma^{L}_{2\times
4}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}(1-z_{2}-z)^{2}(z_{2}+z)\left[z_{2}^{2}+(z_{2}+z)^{2}\right]\int\differential^{8}\mathbf{x}K_{0}(|\mathbf{x}_{12}|Q_{1})K_{0}(QX_{1}^{\prime})\Delta^{(3)}_{22^{\prime}}$
$\displaystyle[S_{23}S_{2^{\prime}3}-S_{13}S_{2^{\prime}3}-S_{12}+1]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}e^{i\frac{z}{z_{2}}\mathbf{q}\cdot(\mathbf{x}_{3}-\mathbf{x}_{2})}$
(13)

$\displaystyle\frac{\differential\sigma^{L}_{6}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}(1-z_{2})^{2}}{(2\pi)^{8}}\int_{0}^{z_{2}}\frac{\differential
z}{z}\int\differential^{8}\mathbf{x}[S_{32^{\prime}}S_{23}-S_{13}S_{23}-S_{12^{\prime}}+1]\left[z_{2}^{2}+(z_{2}-z)^{2}\right]$
$\displaystyle\frac{K_{0}(QX_{6})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})}{\mathbf{x}_{32}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}e^{-i\frac{z}{z_{2}}\mathbf{q}\cdot(\mathbf{x}_{3}-\mathbf{x}_{2})}$
(14)
$\displaystyle\frac{\differential\sigma^{L}_{8}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}(1-z_{2})}{(2\pi)^{8}}\int_{0}^{z_{2}}\frac{\differential
z\,(z_{2}-z)}{z}\int\differential^{8}\mathbf{x}\,[S_{32^{\prime}}S_{23}-S_{13}S_{23}-S_{12^{\prime}}+1]\left[z_{2}(1-z_{2})+(z_{2}-z)(1-z_{2}+z)\right]$
$\displaystyle
K_{0}(QX_{6})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})\Delta^{(3)}_{12}e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}e^{-i\frac{z}{z_{2}}\mathbf{q}\cdot(\mathbf{x}_{3}-\mathbf{x}_{2})}$
(15)
$\displaystyle\frac{\differential\sigma_{10}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{-e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}(1-z_{2})^{2}}{(2\pi)^{6}}\int\differential^{6}\mathbf{x}\big{[}S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1\big{]}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
$\displaystyle\times\int_{0}^{z_{2}}\frac{\differential
z}{z}\,\left[z_{2}^{2}+(z_{2}-z)^{2}\right]\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{1}{\left(\mathbf{k}-\frac{z}{z_{2}}\mathbf{q}\right)^{2}}$
(16)
$\displaystyle\frac{\differential\sigma_{11}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{3}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}\big{[}S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1\big{]}K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
$\displaystyle\int_{0}^{1-z_{2}}\frac{\differential
z}{z}[(1-z_{2})^{2}+(1-z_{2}-z)^{2}]\int\frac{\differential^{2}\mathbf{k}_{2}}{(2\pi)^{2}}\int\frac{\differential^{2}\mathbf{k}_{1}}{(2\pi)^{2}}$
$\displaystyle\frac{e^{i\mathbf{k}_{1}\cdot(\mathbf{x}_{1}-\mathbf{x}_{2})}}{\big{[}\mathbf{k}_{1}^{2}+Q_{2}^{2}\big{]}\left[\left(\mathbf{k}_{2}-\frac{z}{1-z_{2}}\mathbf{k}_{1}\right)^{2}+\frac{z(1-z_{2}-z)}{(1-z_{2})^{2}z_{2}}\mathbf{k}_{1}^{2}+\frac{z}{1-z_{2}}(1-z_{2}-z)Q^{2}\right]}$
(17)
$\displaystyle\frac{\differential\sigma_{12}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}(1-z_{2})^{2}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}\big{[}S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1\big{]}K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})\,e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
$\displaystyle\int_{0}^{z_{2}}\frac{\differential
z}{z}\left[z_{2}^{2}+(z_{2}-z)^{2}\right]\int\frac{\differential^{2}\mathbf{k}_{2}}{(2\pi)^{2}}\int\frac{\differential^{2}\mathbf{k}_{1}}{(2\pi)^{2}}\frac{e^{i\mathbf{k}_{1}\cdot(\mathbf{x}_{1}-\mathbf{x}_{2})}}{\big{[}\mathbf{k}_{1}^{2}+Q_{2}^{2}\big{]}\left[\left(\mathbf{k}_{2}-\frac{z}{z_{2}}\mathbf{k}_{1}\right)^{2}+\frac{z(z_{2}-z)}{(1-z_{2})z_{2}^{2}}\mathbf{k}_{1}^{2}+\frac{z}{z_{2}}(z_{2}-z)Q^{2}\right]}$
(18)
$\displaystyle\frac{\differential\sigma_{13(1)}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{2}(1-z_{2})}{(2\pi)^{6}}\int_{0}^{z_{2}}\differential
z(1-z_{2}+z)(z_{2}-z)\int\differential^{8}\mathbf{x}[S_{12}S_{1^{\prime}2^{\prime}}-S_{12}-S_{1^{\prime}2^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}$
$\displaystyle
K_{0}\left(|\mathbf{x}_{12}|Q\sqrt{(1-z_{2}+z)(z_{2}-z)}\right)K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{2})\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}e^{i\mathbf{k}\cdot\mathbf{x}_{21}}$
$\displaystyle\int\frac{\differential^{2}\mathbf{p}}{(2\pi)^{2}}\Bigg{[}\frac{\frac{z_{2}(z_{2}-z)}{z}}{(z_{2}\mathbf{k}-z\mathbf{q})^{2}}+z_{2}z\frac{\left(\mathbf{k}-\mathbf{q}+\frac{(z_{2}-z)}{1-z_{2}}\mathbf{p}\right)\cdot\left(\mathbf{k}+\mathbf{p}-\frac{(1-z_{2}+z)}{z_{2}}\mathbf{q}\right)}{(1-z_{2}+z)(z_{2}\mathbf{k}-(1-z_{2})\mathbf{q})^{2})\left[\frac{(z_{2}\mathbf{k}-(1-z_{2})\mathbf{q})^{2}}{z_{2}(z_{2}-z)}-\frac{((1-z_{2})\mathbf{k}-z\mathbf{p})^{2}}{(1-z_{2})(1-z_{2}+z)}\right]}\Bigg{]}e^{i\mathbf{p}\cdot\mathbf{x}_{1^{\prime}1}}$
(19)
$\displaystyle\frac{\differential\sigma_{13(2)}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{2}(1-z_{2})}{(2\pi)^{6}}\int_{0}^{1-z_{2}}\differential
z(1-z_{2}-z)(z_{2}+z)\int\differential^{8}\mathbf{x}[S_{12}S_{1^{\prime}2^{\prime}}-S_{12}-S_{1^{\prime}2^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}$
$\displaystyle
K_{0}\left(|\mathbf{x}_{12}|Q_{1}\right)K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{2})\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}e^{i\mathbf{k}\cdot\mathbf{x}_{12}}$
$\displaystyle\int\frac{\differential^{2}\mathbf{p}}{(2\pi)^{2}}\Bigg{[}\frac{\frac{(1-z_{2})(1-z_{2}-z)}{z}}{((1-z_{2})\mathbf{k}-z\mathbf{p})^{2}}+(1-z_{2})z\frac{\left(\mathbf{k}-\mathbf{p}+\frac{(1-z_{2}-z)}{z_{2}}\mathbf{q}\right)\cdot\left(\mathbf{k}+\mathbf{q}-\frac{(z_{2}+z)}{1-z_{2}}\mathbf{p}\right)}{(z_{2}+z)((1-z_{2})\mathbf{k}-z\mathbf{p})^{2}\left[\frac{((1-z_{2})\mathbf{k}-z\mathbf{p})^{2}}{(1-z_{2})(1-z_{2}-z)}-\frac{(z_{2}\mathbf{k}-z\mathbf{q})^{2}}{z_{2}(z_{2}+z)}\right]}\Bigg{]}e^{i\mathbf{p}\cdot\mathbf{x}_{1^{\prime}1}}.$
(20)

$\displaystyle\frac{\differential\sigma_{14(1)}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{2}(1-z_{2})}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
$\displaystyle\int_{0}^{1-z_{2}}\frac{\differential
z}{z}\frac{\differential^{2}\mathbf{k}_{1}}{(2\pi)^{2}}\frac{\differential^{2}\mathbf{k}_{2}}{(2\pi)^{2}}e^{i\mathbf{k}_{2}\cdot\mathbf{x}_{12}}\Bigg{[}\frac{z_{2}(1-z_{2})+(1-z_{2}-z)(z_{2}+z)}{\Big{[}\mathbf{k}_{2}^{2}+Q_{2}^{2}\Big{]}\left[\left(\mathbf{k}_{1}-\frac{1-z_{2}-z}{1-z_{2}}\mathbf{k}_{2}\right)^{2}+\frac{z(1-z_{2}-z)}{z_{2}(1-z_{2})^{2}}\mathbf{k}_{2}^{2}+\frac{z}{1-z_{2}}(1-z_{2}-z)Q^{2}\right]}$
$\displaystyle+\frac{\frac{(1-z_{2}-z)(z_{2}+z)}{z_{2}(1-z_{2})}((z_{2}+z)(1-z_{2}-z)+z_{2}(1-z_{2})))}{\Big{[}\mathbf{k}_{1}^{2}+(1-z_{2}-z)(z_{2}+z)Q^{2}\Big{]}\left[\left(\mathbf{k}_{1}-\frac{1-z_{2}-z}{1-z_{2}}\mathbf{k}_{2}\right)^{2}+\frac{z(1-z_{2}-z)}{z_{2}(1-z_{2})^{2}}\mathbf{k}_{2}^{2}+\frac{z}{1-z_{2}}(1-z_{2}-z)Q^{2}\right]}$
$\displaystyle-\frac{\frac{z(1-z_{2}-z)}{1-z_{2}}\left[z_{2}(1-z_{2})+(1-z_{2}-z)(z_{2}+z)\right]Q^{2}}{\Big{[}\mathbf{k}_{1}^{2}+(1-z_{2}-z)(z_{2}+z)Q^{2}\Big{]}\Big{[}\mathbf{k}_{2}^{2}+Q_{2}^{2}\Big{]}\left[\left(\mathbf{k}_{1}-\frac{1-z_{2}-z}{1-z_{2}}\mathbf{k}_{2}\right)^{2}+\frac{z(1-z_{2}-z)}{z_{2}(1-z_{2})^{2}}\mathbf{k}_{2}^{2}+\frac{z}{1-z_{2}}(1-z_{2}-z)Q^{2}\right]}$
$\displaystyle-\frac{z\left[(1-z_{2})(1-z_{2}-z)+z_{2}(z_{2}+z)\right]}{\Big{[}\mathbf{k}_{1}^{2}+(1-z_{2}-z)(z_{2}+z)Q^{2}\Big{]}\Big{[}\mathbf{k}_{2}^{2}+Q_{2}^{2}\Big{]}}\Bigg{]}$
(21)
$\displaystyle\frac{\differential\sigma_{14(2)}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{2}(1-z_{2})}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
$\displaystyle\int_{0}^{z_{2}}\frac{\differential
z}{z}\frac{\differential^{2}\mathbf{k}_{1}}{(2\pi)^{2}}\frac{\differential^{2}\mathbf{k}_{2}}{(2\pi)^{2}}e^{i\mathbf{k}_{2}\cdot\mathbf{x}_{12}}\Bigg{[}\frac{z_{2}(1-z_{2})+(1-z_{2}+z)(z_{2}-z)}{\Big{[}\mathbf{k}_{2}^{2}+Q_{2}^{2}\Big{]}\left[\left(\mathbf{k}_{1}-\frac{z_{2}-z}{z_{2}}\mathbf{k}_{2}\right)^{2}+\frac{z(z_{2}-z)}{(1-z_{2})z_{2}^{2}}\mathbf{k}_{2}^{2}+\frac{z}{z_{2}}(z_{2}-z)Q^{2}\right]}$
$\displaystyle+\frac{\frac{(1-z_{2}+z)(z_{2}-z)}{z_{2}(1-z_{2})}\left[z_{2}(1-z_{2})+(1-z_{2}+z)(z_{2}-z)\right]}{\Big{[}\mathbf{k}_{1}^{2}+(1-z_{2}+z)(z_{2}-z)Q^{2}\Big{]}\left[\left(\mathbf{k}_{1}-\frac{z_{2}-z}{z_{2}}\mathbf{k}_{2}\right)^{2}+\frac{z(z_{2}-z)}{(1-z_{2})z_{2}^{2}}\mathbf{k}_{2}^{2}+\frac{z}{z_{2}}(z_{2}-z)Q^{2}\right]}$
$\displaystyle-\frac{\frac{z(z_{2}-z)}{z_{2}}\left[z_{2}(1-z_{2})+(1-z_{2}+z)(z_{2}-z)\right]Q^{2}}{\Big{[}\mathbf{k}_{1}^{2}+(1-z_{2}+z)(z_{2}-z)Q^{2}\Big{]}\Big{[}\mathbf{k}_{2}^{2}+Q_{2}^{2}\Big{]}\left[\left(\mathbf{k}_{1}-\frac{z_{2}-z}{z_{2}}\mathbf{k}_{2}\right)^{2}+\frac{z(z_{2}-z)}{(1-z_{2})z_{2}^{2}}\mathbf{k}_{2}^{2}+\frac{z}{z_{2}}(z_{2}-z)Q^{2}\right]}$
$\displaystyle-\frac{z\left[z_{2}(z_{2}-z)+(1-z_{2})(1-z_{2}+z)\right]}{\Big{[}\mathbf{k}_{1}^{2}+(1-z_{2}+z)(z_{2}-z)Q^{2}\Big{]}\Big{[}\mathbf{k}_{2}^{2}+Q_{2}^{2}\Big{]}}\Bigg{]}$
(22)

where we have
$\mathbf{x}_{1}^{\prime}\equiv\mathbf{x}_{1}+\frac{z}{1-z_{2}}(\mathbf{x}_{3}-\mathbf{x}_{1})$
in $\sigma_{1\times 2}$. Here we have also defined

$\displaystyle\Delta^{(3)}_{ij}=\frac{\mathbf{x}_{3i}\cdot\mathbf{x}_{3j}}{\mathbf{x}_{3i}^{2}\mathbf{x}_{3j}^{2}}.$
(23)

We also note that since we have integrated over $z_{1}$ the definition of
$Q_{2}$ remains the same but $Q_{1}$ (which is still used in some expressions)
has now changed.

$\displaystyle
Q_{2}=Q\sqrt{z_{2}(1-z_{2})},\,\,\,\,\,\,\,\,Q_{1}=Q\sqrt{(1-z_{2}-z)(z_{2}+z)}.$
(24)

We have also used a shorthand notation for the coordinate dependence in some
of the Bessel functions.

$\displaystyle X$
$\displaystyle=\sqrt{(1-z_{2}-z)z_{2}\mathbf{x}_{12}^{2}+(1-z_{2}-z)z\mathbf{x}_{13}^{2}+z_{2}z\mathbf{x}_{23}^{2}},$
$\displaystyle X_{6}$
$\displaystyle=\sqrt{(1-z_{2})(z_{2}-z)\mathbf{x}_{12}^{2}+(1-z_{2})z\,\mathbf{x}_{13}^{2}+z(z_{2}-z)\mathbf{x}_{23}^{2}}.$
(25)

$X^{\prime}$ is the same as $X$ but with primed coordinates (except for
$\mathbf{x}_{3}$ which is never primed). $X^{\prime}_{1}$ is the same as
$X^{\prime}$ but $\mathbf{x}_{1}^{\prime}$ has become unprimed.

These expressions provide the formal results for the one-loop corrections to
single inclusive hadron production. Looking at the results (Eq. 7 \- 22) one
can see that some corrections appear to be missing. In particular, we have not
written $\sigma_{1\times 1},\sigma_{1\times 3},\sigma_{1\times
4},\sigma_{5},\sigma_{7},$ and $\sigma_{9}$. This is because $\sigma_{1\times
1}$ exactly cancels $\sigma_{9}$ ($\sigma_{9}$ gets an extra factor of $2$ due
to it being a cross term). Similarly, $\sigma_{1\times 3}$ cancels
$\sigma_{5}$, and also $\sigma_{1\times 4}$ cancels $\sigma_{7}$. We include
the expressions for these in appendix A for completeness. To understand this
cancellation, one can draw these corrections in cut diagram notation. In Fig.
3 we show that $\sigma_{1\times 1}$ and $\sigma_{9}$ become diagrammatically
identical when one integrates the final state quark. Similarly, one finds that
$\sigma_{1\times 3}$ becomes the same diagram as $\sigma_{5}$, and
$\sigma_{1\times 4}$ becomes the same as $\sigma_{7}$ (see Fig. 4).

Figure 3: When one integrates out the final state quark, the real correction
$\sigma_{1\times 1}$ and the virtual correction $\sigma_{9}$ become the same
diagram. Here the ‘x’ at the end of a solid line indicates produced quarks and
antiquarks.

Figure 4: Here the left diagram is $\sigma_{1\times 3}$ and $\sigma_{5}$, and
the right diagram is $\sigma_{1\times 4}$ and $\sigma_{7}$. The ‘x’ indicates
the produced antiquark.

Therefore since the expressions corresponding to these diagrams differ
mathematically only by a sign, they all cancel each other completely.

$\displaystyle\differential\sigma_{1\times 1}+2\differential\sigma_{9}=0,$
$\displaystyle\differential\sigma_{1\times 3}+\differential\sigma_{5}=0,$
$\displaystyle\differential\sigma_{1\times 4}+\differential\sigma_{7}=0.$ (26)

Note that the factor of 2 on $\differential\sigma_{9}$ comes due to the fact
that it’s a cross term and therefore gets double counted relative to
$\differential\sigma_{1\times 1}$. So it is justified to ignore these
corrections (see appendix A for the expressions). In the remaining terms, one
finds divergences which must either be canceled or absorbed into the
renormalization of physical parameters. This is the topic of the next section.

##  IV Divergences

As in the case of dihadron production in Bergabo and Jalilian-Marian (2022b)
there are four categories of divergences; UV, soft, rapidity and collinear.
The cancellation of UV and soft divergences proceeds as follows. UV
divergences appear as $\mathbf{x}_{3}\to\mathbf{x}_{2}$ in $\sigma_{6}$, as
$\mathbf{x}_{3}\to\mathbf{x}_{1}$ in $\sigma_{3\times 3}$, as
$\mathbf{k}\to\infty$ in $\sigma_{10}$, as $\mathbf{k}_{2}\to\infty$ in
$\sigma_{11}$ and $\sigma_{12}$, and as $\mathbf{k}_{1}\to\infty$ in
$\sigma_{14(1)}$ and $\sigma_{14(2)}$. We find that these UV divergences all
cancel according to Eq. 27. Note that all other UV divergent terms were
already canceled in the first two lines of Eq. 26.

$\displaystyle\left[\differential\sigma_{6}+\differential\sigma_{12}\right]_{\text{UV}}=0,$
$\displaystyle\left[\differential\sigma_{3\times
3}+\differential\sigma_{11}+\differential\sigma_{10}+\differential\sigma_{14(1)}+\differential\sigma_{14(2)}\right]_{\text{UV}}=0.$
(27)

Soft divergences occur when $\mathbf{k}$ and $z$ both go to zero (in
coordinate space this becomes $\mathbf{x}_{3}\to\infty$ and $z\to 0$). The
cancellation of soft divergences in all remaining terms proceeds identically
to that in dihadron production Bergabo and Jalilian-Marian (2022b) with

$\displaystyle\left[\differential\sigma_{2\times
2}+2\,\differential\sigma_{10}\right]_{\text{soft}}=0,$
$\displaystyle\left[\differential\sigma_{1\times
2}+\differential\sigma_{13(1)}+\differential\sigma_{13(2)}\right]_{\text{soft}}=0,$
$\displaystyle\left[\differential\sigma_{3\times
3}+\differential\sigma_{4\times 4}+2\,\differential\sigma_{3\times
4}\right]_{\text{soft}}=0,$ $\displaystyle\left[\differential\sigma_{2\times
3}+\differential\sigma_{2\times 4}\right]_{\text{soft}}=0,$
$\displaystyle\left[\differential\sigma_{6}+\differential\sigma_{8}\right]_{\text{soft}}=0,$
$\displaystyle\left[\differential\sigma_{11}+\differential\sigma_{14(1)}\right]_{\text{soft}}=0,$
$\displaystyle\left[\differential\sigma_{12}+\differential\sigma_{14(2)}\right]_{\text{soft}}=0.$
(28)

Rapidity divergences appear when $z\rightarrow 0$ at finite (non-zero)
transverse momentum. To isolate those it is customary to introduce a rapidity
factorization scale $z_{f}$ and write the $z$ integral as

$\displaystyle\int_{0}^{1}\frac{\differential
z}{z}f(z)=\left\\{\int_{0}^{z_{f}}\frac{\differential
z}{z}+\int_{z_{f}}^{1}\frac{\differential z}{z}\,\right\\}\,f(z).$ (29)

so that the rapidity divergences will come from the first integral while the
cross section in the second integral will contain no rapidity divergence.
Therefore we first focus on the first integration region containing the
rapidity divergence. As all soft divergences have already been canceled only
terms of the form
$\frac{\mathbf{x}_{ij}^{2}}{\mathbf{x}_{3i}^{2}\mathbf{x}_{3j}^{2}}$ with
$i\neq j$ remain (see Bergabo and Jalilian-Marian (2022b) for the explicit
expressions). These can be added to give

$\displaystyle\frac{\differential\sigma^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{4e^{2}Q^{2}N_{c}z_{2}^{3}(1-z_{2})^{2}}{(2\pi)^{7}}\,\int\differential^{6}\mathbf{x}K_{0}(|\mathbf{x}_{12}|Q_{1})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{1})$
$\displaystyle\bigg{\\{}\frac{N_{c}\,\alpha_{s}}{2\pi^{2}}\,\int_{0}^{z_{f}}\frac{\differential
z}{z}\int\differential^{2}\mathbf{x}_{3}\bigg{[}\frac{\mathbf{x}_{22^{\prime}}^{2}}{\mathbf{x}_{32}^{2}\mathbf{x}_{32^{\prime}}^{2}}\left(S_{23}S_{2^{\prime}3}-S_{22^{\prime}}\right)-\frac{\mathbf{x}_{12}^{2}}{\mathbf{x}_{31}^{2}\mathbf{x}_{32}^{2}}\left(S_{13}S_{32}-S_{12}\right)-\frac{\mathbf{x}_{12^{\prime}}^{2}}{\mathbf{x}_{31}^{2}\mathbf{x}_{32^{\prime}}^{2}}\left(S_{13}S_{2^{\prime}3}-S_{12^{\prime}}\right)\bigg{]}\bigg{\\}}.$
(30)

Comparing this to the LO result in Eq. 5 it is clear that the terms inside the
curly bracket correspond to the BK/JIMWLK evolution Balitsky (1996); Kovchegov
(2000); Jalilian-Marian et al. (1997b, 1998a, 1998b) of the dipoles that
appear in the LO cross section. As the contribution of the second term in the
$z$ integral contains no rapidity divergence this shows that all rapidity
divergences can be absorbed into BK/JIMWLK evolution of the LO cross section.
The contribution of the $\int_{z_{f}}^{1}\differential z$ region is now free
of divergences (after absorbing the collinear divergences into scale dependent
fragmentation functions done in the following pages) and constitute the NLO
correction to the LO result.

Collinear divergences are identical to the case of dihadrons and were treated
in full detail in Bergabo and Jalilian-Marian (2022b), therefore here we will
skip some details. When we integrate out quarks there remains collinear
divergences involving the antiquark when its transverse momentum $\mathbf{q}$
becomes parallel to the loop transverse momentum $\mathbf{k}$
($\theta\rightarrow 0$ with
$\cos\theta\equiv\frac{\mathbf{q}\cdot\mathbf{k}}{|\mathbf{q}||\mathbf{k}|}$
at finite transverse momenta). One can write the hadronic cross section in
terms of the partonic cross section as follows,

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to
hX}}{\differential^{2}\mathbf{q}_{h}\,\differential
y_{2}}=2\int_{0}^{1}\frac{\differential
z_{h}}{z_{h}^{2}}\frac{\differential\sigma^{\gamma^{*}A\to\bar{q}X}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}D_{h/\bar{q}}(z_{h}).$ (32)

Here we have included only the case where the quark is integrated out and the
antiquark with momentum $q$ fragments into a hadron $h$ with with momentum
$z_{h}q$. As mentioned earlier, the opposite case (antiquark integrated out
and the quark fragments into a hadron) is mathematically the same after a
relabeling of some variables. Therefore we account for the opposite case with
the overall factor of 2 above. The renormalized fragmentation function
$D_{h/\bar{q}}(z_{h},\mu^{2})$ can be written in terms of the _bare_
fragmentation function $D^{0}_{h/\bar{q}}$ and higher order corrections.

$\displaystyle
D_{h/\bar{q}}(z_{h},\mu^{2})=D^{0}_{h/\bar{q}}(z_{h})+\mathcal{O}(\alpha_{s})$
(33)

The $\mathcal{O}(\alpha_{s})$ corrections are expected to come from
$\sigma_{2\times 2}$ (Eq. 7) and $\sigma_{10}$ (Eq. 16). Using Eq. 32 with the
leading order partonic cross section (longitudinal part of Eq. 5) and the two
relevant corrections (in momentum space), we can write

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to hX}_{\text{LO}+2\times
2+10}}{\differential^{2}\mathbf{q}_{h}\,\differential y_{2}}=$
$\displaystyle\int_{0}^{1}\frac{\differential
z_{h}}{z_{h}^{2}}D^{0}_{h/\bar{q}}(z_{h})\frac{2e^{2}Q^{2}N_{c}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})$
$\displaystyle\Bigg{[}4z_{2}^{3}(1-z_{2})^{2}+\frac{2g^{2}N_{c}}{(2\pi)}\int\frac{\differential
z}{z}\frac{(1-z_{2}-z)^{2}(z_{2}+z)^{2}[z_{2}^{2}+(z_{2}+z)^{2}]}{z_{2}}\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{e^{i\mathbf{k}\cdot\mathbf{x}_{2^{\prime}2}}}{\left(\mathbf{k}-\frac{z}{z_{2}}\mathbf{q}\right)^{2}}$
$\displaystyle-\frac{2g^{2}N_{c}}{(2\pi)}\int_{0}^{z_{2}}\frac{\differential
z}{z}z_{2}(1-z_{2})^{2}[z_{2}^{2}+(z_{2}-z)^{2}]\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{1}{\left(\mathbf{k}-\frac{z}{z_{2}}\mathbf{q}\right)^{2}}\Bigg{]}.$
(34)

Here the first term in the square brackets is the leading order contribution,
the second term is the real correction $\sigma_{2\times 2}$, and the last term
is the virtual correction $\sigma_{10}$ (doubled here since it’s a cross
term). Next we follow identical steps to what was done for dihadrons and
rewrite this expression using $g^{2}=4\pi\alpha_{s}$ and relax the
large-$N_{c}$ approximation inside the square brackets taking $N_{c}\to
2C_{F}$. We also define a new variable $\xi$ for the $z$ integration in both
the real and virtual corrections. For the real correction
$\xi=z_{2}/(z_{2}+z)$ and for the virtual correction $\xi=(z_{2}-z)/z_{2}$.

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to hX}_{\text{LO}+2\times
2+10}}{\differential^{2}\mathbf{q}_{h}\,\differential y_{2}}=$
$\displaystyle\int_{0}^{1}\frac{\differential
z_{h}}{z_{h}^{2}}D^{0}_{h/\bar{q}}(z_{h})\frac{8e^{2}Q^{2}N_{c}z_{2}^{3}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})$
$\displaystyle\Bigg{[}(1-z_{2})^{2}+2\alpha_{s}C_{F}\int\frac{\differential\xi}{\xi^{5}}\frac{(1-z_{2}/\xi)^{2}(1+\xi^{2})}{(1-\xi)}\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{e^{i\mathbf{k}\cdot\mathbf{x}_{2^{\prime}2}}}{\left(\mathbf{k}-\frac{(1-\xi)}{\xi}\mathbf{q}\right)^{2}}$
$\displaystyle-2(1-z_{2})^{2}\alpha_{s}C_{F}\int_{0}^{1}\differential\xi\frac{(1+\xi^{2})}{(1-\xi)}\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{1}{\left(\mathbf{k}-(1-\xi)\mathbf{q}\right)^{2}}\Bigg{]}.$
(35)

Now, we want to regulate the divergences in the two integrals over
$\mathbf{k}$. To do so we follow Bergabo and Jalilian-Marian (2022b) and write

$\displaystyle\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{e^{i\mathbf{k}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}}{\left(\mathbf{k}-\frac{(1-\xi)}{\xi}\mathbf{q}\right)^{2}}\to\frac{e^{i\frac{(1-\xi)}{\xi}\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}}{2\pi}\left[\frac{1}{\epsilon}-\log\left(\pi
e^{\gamma_{E}}\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]+\mathcal{O}(\epsilon),\,\,\,\,\,\epsilon=d-2>0.$
(36)

$\displaystyle\left[\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{1}{\mathbf{k}^{2}}\right]_{IR}=\frac{1}{2\pi}\left[-\frac{1}{\epsilon_{IR}}-\log\left(e^{\gamma_{E}}\pi\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]+\mathcal{O}(\epsilon).$
(37)

$\displaystyle\left[\int\frac{\differential^{2}\mathbf{k}}{(2\pi)^{2}}\frac{1}{\mathbf{k}^{2}}\right]_{UV}=\frac{1}{2\pi}\left[\frac{1}{\epsilon_{UV}}+\log\left(e^{\gamma_{E}}\pi\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]+\mathcal{O}(\epsilon).$
(38)

Here Eq. 36 is used for the real integral. We also make the approximation
$\xi\approx 1$ inside the exponential on the right side to ignore it, this we
motivate by noting that there is a $1-\xi$ in the denominator so the integral
is dominated by the region $\xi\approx 1$. Eq. 37 and 38 are used for the
virtual correction, after a shift on the $\mathbf{k}$ integral. This shift
makes the collinear divergence look infrared, but this is distinct from the
soft divergences discussed earlier and is in fact the collinear divergence.
The UV divergence is canceled against other virtual corrections (Eq. 27).
Therefore what remains here is the finite part of Eq.38 which remains as part
of the finite NLO corrections, and Eq. 37 which can now be added to the real
correction by setting $\epsilon_{IR}=-\epsilon$.

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to hX}_{\text{LO}+2\times
2+10}}{\differential^{2}\mathbf{q}_{h}\,\differential
y_{2}}=\int_{0}^{1}\frac{\differential
z_{h}}{z_{h}^{2}}D^{0}_{h/\bar{q}}(z_{h})\frac{8e^{2}Q^{2}N_{c}z_{2}^{3}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})$
$\displaystyle\Bigg{[}(1-z_{2})^{2}+\left\\{\frac{\alpha_{s}C_{F}}{\pi}\int\frac{\differential\xi}{\xi^{5}}\frac{(1-z_{2}/\xi)^{2}(1+\xi^{2})}{(1-\xi)}-\frac{(1-z_{2})^{2}\alpha_{s}C_{F}}{\pi}\int_{0}^{1}\differential\xi\frac{(1+\xi^{2})}{(1-\xi)}\right\\}\left[\frac{1}{\epsilon}-\log\left(\pi
e^{\gamma_{E}}\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]\Bigg{]}$
(39)

Next, let’s note that $z_{2}$ is no longer an independent variable, it can be
written in terms of $z_{h}$ using

$\displaystyle z_{2}=\frac{q_{h}^{+}}{z_{h}l^{+}}.$ (40)

Writing all the $z_{2}$’s this way, and rearranging the result we have

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to hX}_{\text{LO}+2\times
2+10}}{\differential^{2}\mathbf{q}_{h}\,\differential
y_{2}}=\int_{0}^{1}\differential
z_{h}\frac{8e^{2}Q^{2}N_{c}(q_{h}^{+})^{3}}{(2\pi)^{5}(l^{+})^{3}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})$
$\displaystyle\Bigg{[}\frac{\left(1-\frac{q_{h}^{+}}{z_{h}l^{+}}\right)^{2}}{z_{h}^{5}}D^{0}_{h/\bar{q}}(z_{h})+\left\\{\frac{\alpha_{s}C_{F}}{\pi}\int\frac{\differential\xi}{\xi^{5}}\frac{(1-\frac{q_{h}^{+}}{z_{h}l^{+}\xi})^{2}(1+\xi^{2})}{(1-\xi)}\frac{D^{0}_{h/\bar{q}}(z_{h})}{z_{h}^{5}}-\frac{(1-\frac{q_{h}^{+}}{z_{h}l^{+}})^{2}\alpha_{s}C_{F}}{\pi}\int_{0}^{1}\differential\xi\frac{(1+\xi^{2})}{(1-\xi)}\frac{D^{0}_{h/\bar{q}}(z_{h})}{z_{h}^{5}}\right\\}$
$\displaystyle\times\left[\frac{1}{\epsilon}-\log\left(\pi
e^{\gamma_{E}}\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]\Bigg{]}.$
(41)

Now that all the $z_{h}$ dependence is explicit, let’s perform a substitution
on the $z_{h}$ integral only in the real correction term (first term in the
curly brackets). We’ll define $z_{h}^{\prime}=\xi z_{h}$, and once the
substitution is complete we’ll remove the prime. This yields

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to hX}_{\text{LO}+2\times
2+10}}{\differential^{2}\mathbf{q}_{h}\,\differential
y_{2}}=\int_{0}^{1}\frac{\differential
z_{h}}{z_{h}^{2}}\frac{8e^{2}Q^{2}N_{c}z_{2}^{3}(1-z_{2})^{2}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})$
$\displaystyle\Bigg{[}D^{0}_{h/\bar{q}}(z_{h})+\left\\{\frac{\alpha_{s}C_{F}}{\pi}\int_{z_{h}}^{1}\frac{\differential\xi}{\xi}\frac{(1+\xi^{2})}{(1-\xi)}D^{0}_{h/\bar{q}}(z_{h}/\xi)-\frac{\alpha_{s}C_{F}}{\pi}\int_{0}^{1}\differential\xi\frac{(1+\xi^{2})}{(1-\xi)}D^{0}_{h/\bar{q}}(z_{h})\right\\}\left[\frac{1}{\epsilon}-\log\left(\pi
e^{\gamma_{E}}\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]\Bigg{]}.$
(42)

Here we have written things back in terms of $z_{2}$ for aesthetics. Finally,
the real and virtual corrections can be combined using the antiquark-antiquark
splitting function $\mathcal{P}_{\bar{q}\bar{q}}$ defined as

$\displaystyle\mathcal{P}_{\bar{q}\bar{q}}(\xi)=C_{F}\left[\frac{(1+\xi^{2})}{(1-\xi)_{+}}+\frac{3}{2}\delta(1-\xi)\right]$
(43)

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to hX}_{\text{LO}+2\times
2+10}}{\differential^{2}\mathbf{q}_{h}\,\differential
y_{2}}=\int_{0}^{1}\frac{\differential
z_{h}}{z_{h}^{2}}\frac{8e^{2}Q^{2}N_{c}z_{2}^{3}(1-z_{2})^{2}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})$
$\displaystyle\int_{z_{h}}^{1}\frac{\differential\xi}{\xi}D^{0}_{h/\bar{q}}(z_{h}/\xi)\Bigg{[}\delta(1-\xi)+\frac{\alpha_{s}}{\pi}P_{\bar{q}\bar{q}}(\xi)\left[\frac{1}{\epsilon}-\log\left(\pi
e^{\gamma_{E}}\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]\Bigg{]}.$
(44)

For more details on these calculations, see Bergabo and Jalilian-Marian
(2022b). So, we are now able to define the DGLAP evolved fragmentation
function

$\displaystyle
D_{h/\bar{q}}(z_{h},\mu^{2})=\int_{z_{h}}^{1}\frac{\differential\xi}{\xi}D^{0}_{h/\bar{q}}(z_{h}/\xi)\Bigg{[}\delta(1-\xi)+\frac{\alpha_{s}}{\pi}P_{\bar{q}\bar{q}}(\xi)\left[\frac{1}{\epsilon}-\log\left(\pi
e^{\gamma_{E}}\mu|\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2}|\right)\right]\Bigg{]},$
(45)

in terms of which our result becomes

$\displaystyle\frac{\differential\sigma^{\gamma^{*}A\to hX}_{\text{LO}+2\times
2+10}}{\differential^{2}\mathbf{q}_{h}\,\differential
y_{2}}=\int_{0}^{1}\frac{\differential
z_{h}}{z_{h}^{2}}\frac{8e^{2}Q^{2}N_{c}z_{2}^{3}(1-z_{2})^{2}}{(2\pi)^{5}}\int\differential^{6}\mathbf{x}[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1]e^{i\mathbf{q}\cdot\mathbf{x}_{2^{\prime}2}}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})D_{h/\bar{q}}(z_{h}).$
(46)

Thus we have shown that all divergences appearing in our next-to-leading order
results are either canceled or absorbed into evolution of dipoles and
fragmentation functions. The full result for the single inclusive cross
section at next-to-leading order can be written schematically as

$\displaystyle\differential\sigma^{\gamma^{*}A\to
hX}=\differential\sigma_{LO}\otimes\text{JIMWLK}+\differential\sigma_{LO}\otimes
D_{h/\bar{q}}(z_{h},\mu^{2})+\differential\sigma_{NLO}^{\text{finite}}.$ (47)

Here we imply the presence of a bare fragmentation function
$D^{0}_{h/\bar{q}}(z_{h})$ in the first and last terms.

In summary we have derived the next to leading order corrections to single
inclusive hadron production in the forward rapidity region in DIS at small
$x$. We have shown that all divergences either cancel among various terms or
can be absorbed into BK/JIMWLK evolution of dipoles and DGLAP evolution of the
parton-hadron fragmentation function.

### Acknowledgement:

We gratefully acknowledge support from the DOE Office of Nuclear Physics
through Grant No. DE-SC0002307 and by PSC-CUNY through grant No. 63158-0051.
We would like to thank Y. Kovchegov, C. Marquet, B. Xiao and F. Yuan for
helpful discussions.

##  Appendix A Cancellations

Here we list the corrections that all cancel each other exactly according to
Eq. 26 and can therefore be neglected from the results.

$\displaystyle\frac{\differential\sigma^{L}_{1\times
1}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{3}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}\left[(1-z_{2}-z)^{2}+(1-z_{2})^{2}\right]\int\differential^{8}\mathbf{x}\,K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})\frac{1}{\mathbf{x}_{31}^{2}}$
$\displaystyle\left[S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1\right]e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(48) $\displaystyle\frac{\differential\sigma^{L}_{1\times
3}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{3}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}\int\differential^{8}\mathbf{x}\,[S_{322^{\prime}1^{\prime}}S_{13}-S_{13}S_{23}-S_{1^{\prime}2^{\prime}}+1]\left[(1-z_{2})^{2}+(1-z_{2}-z)^{2}\right]$
$\displaystyle\frac{K_{0}(QX)K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{2})}{\mathbf{x}_{31}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(49) $\displaystyle\frac{\differential\sigma^{L}_{1\times
4}}{\differential^{2}\mathbf{q}\,\differential y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{2}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z\,(1-z_{2}-z)}{z}\int\differential^{8}\mathbf{x}\,[S_{322^{\prime}1^{\prime}}S_{13}-S_{13}S_{23}-S_{1^{\prime}2^{\prime}}+1]\left[(1-z_{2})z_{2}+(1-z_{2}-z)(z_{2}+z)\right]$
$\displaystyle
K_{0}(QX)K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{2})\Delta^{(3)}_{12}e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(50)
$\displaystyle\frac{\differential\sigma^{L}_{5}}{\differential^{2}\mathbf{q}\differential
y_{2}}=$
$\displaystyle\frac{2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{3}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z}{z}\int\differential^{8}\mathbf{x}\,[S_{322^{\prime}1^{\prime}}S_{13}-S_{13}S_{23}-S_{1^{\prime}2^{\prime}}+1]\left[(1-z_{2})^{2}+(1-z_{2}-z)^{2}\right]$
$\displaystyle\frac{K_{0}(QX)K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{2})}{\mathbf{x}_{31}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(51)
$\displaystyle\frac{\differential\sigma^{L}_{7}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{-2e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{2}}{(2\pi)^{8}}\int_{0}^{1-z_{2}}\frac{\differential
z\,(1-z_{2}-z)}{z}\int\differential^{8}\mathbf{x}\,[S_{322^{\prime}1^{\prime}}S_{13}-S_{13}S_{23}-S_{1^{\prime}2^{\prime}}+1]\left[(1-z_{2})z_{2}+(1-z_{2}-z)(z_{2}+z)\right]$
$\displaystyle
K_{0}(QX)K_{0}(|\mathbf{x}_{1^{\prime}2^{\prime}}|Q_{2})\Delta^{(3)}_{12}e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
(52)
$\displaystyle\frac{\differential\sigma_{9}^{L}}{\differential^{2}\mathbf{q}\,\differential
y_{2}}=$
$\displaystyle\frac{-e^{2}g^{2}Q^{2}N_{c}^{2}z_{2}^{3}}{(2\pi)^{8}}\int\differential^{8}\mathbf{x}\big{[}S_{22^{\prime}}-S_{12}-S_{12^{\prime}}+1\big{]}K_{0}(|\mathbf{x}_{12}|Q_{2})K_{0}(|\mathbf{x}_{12^{\prime}}|Q_{2})e^{i\mathbf{q}\cdot(\mathbf{x}_{2}^{\prime}-\mathbf{x}_{2})}$
$\displaystyle\times\int_{0}^{1-z_{2}}\frac{\differential
z}{z}[(1-z_{2})^{2}+(1-z_{2}-z)^{2}]\frac{1}{\mathbf{x}_{31}^{2}}$ (53)

Here
$\mathbf{x}_{1}^{\prime}\equiv\mathbf{x}_{1}+\frac{z}{1-z_{2}}(\mathbf{x}_{3}-\mathbf{x}_{1})$
in these expressions.

## References

  * Gribov et al. (1983) L. V. Gribov, E. M. Levin, and M. G. Ryskin, Phys. Rept. 100, 1 (1983).
  * Mueller and Qiu (1986) A. H. Mueller and J.-w. Qiu, Nucl. Phys. B 268, 427 (1986).
  * Iancu and Venugopalan (2003) E. Iancu and R. Venugopalan, _The Color glass condensate and high-energy scattering in QCD_ (2003), pp. 249–3363, eprint hep-ph/0303204.
  * Iancu et al. (2002) E. Iancu, A. Leonidov, and L. McLerran, in _Cargese Summer School on QCD Perspectives on Hot and Dense Matter_ (2002), pp. 73–145, eprint hep-ph/0202270.
  * Jalilian-Marian and Kovchegov (2006) J. Jalilian-Marian and Y. V. Kovchegov, Prog. Part. Nucl. Phys. 56, 104 (2006), eprint hep-ph/0505052.
  * Weigert (2005) H. Weigert, Prog. Part. Nucl. Phys. 55, 461 (2005), eprint hep-ph/0501087.
  * Morreale and Salazar (2021) A. Morreale and F. Salazar, Universe 7, 312 (2021), eprint 2108.08254.
  * Kovner and Wiedemann (2001) A. Kovner and U. A. Wiedemann, Phys. Rev. D 64, 114002 (2001), eprint hep-ph/0106240.
  * Jalilian-Marian and Kovchegov (2004) J. Jalilian-Marian and Y. V. Kovchegov, Phys. Rev. D 70, 114017 (2004), [Erratum: Phys.Rev.D 71, 079901 (2005)], eprint hep-ph/0405266.
  * Dumitru et al. (2006) A. Dumitru, A. Hayashigaki, and J. Jalilian-Marian, Nucl. Phys. A 765, 464 (2006), eprint hep-ph/0506308.
  * Jalilian-Marian (2006) J. Jalilian-Marian, Nucl. Phys. A 770, 210 (2006), eprint hep-ph/0509338.
  * Marquet (2007) C. Marquet, Nucl. Phys. A 796, 41 (2007), eprint 0708.0231.
  * Albacete and Marquet (2010) J. L. Albacete and C. Marquet, Phys. Rev. Lett. 105, 162301 (2010), eprint 1005.4065.
  * Stasto et al. (2012) A. Stasto, B.-W. Xiao, and F. Yuan, Phys. Lett. B 716, 430 (2012), eprint 1109.1817.
  * Lappi and Mantysaari (2013) T. Lappi and H. Mantysaari, Nucl. Phys. A 908, 51 (2013), eprint 1209.2853.
  * Jalilian-Marian and Rezaeian (2012a) J. Jalilian-Marian and A. H. Rezaeian, Phys. Rev. D 86, 034016 (2012a), eprint 1204.1319.
  * Jalilian-Marian and Rezaeian (2012b) J. Jalilian-Marian and A. H. Rezaeian, Phys. Rev. D 85, 014017 (2012b), eprint 1110.2810.
  * Zheng et al. (2014) L. Zheng, E. C. Aschenauer, J. H. Lee, and B.-W. Xiao, Phys. Rev. D 89, 074037 (2014), eprint 1403.2413.
  * Stasto et al. (2018) A. Stasto, S.-Y. Wei, B.-W. Xiao, and F. Yuan, Phys. Lett. B 784, 301 (2018), eprint 1805.05712.
  * Albacete et al. (2019) J. L. Albacete, G. Giacalone, C. Marquet, and M. Matas, Phys. Rev. D 99, 014002 (2019), eprint 1805.05711.
  * Mäntysaari et al. (2020) H. Mäntysaari, N. Mueller, F. Salazar, and B. Schenke, Phys. Rev. Lett. 124, 112301 (2020), eprint 1912.05586.
  * Hatta et al. (2021a) Y. Hatta, B.-W. Xiao, F. Yuan, and J. Zhou, Phys. Rev. Lett. 126, 142001 (2021a), eprint 2010.10774.
  * Jia et al. (2020) J. Jia, S.-Y. Wei, B.-W. Xiao, and F. Yuan, Phys. Rev. D 101, 094008 (2020), eprint 1910.05290.
  * Gelis and Jalilian-Marian (2002) F. Gelis and J. Jalilian-Marian, Phys. Rev. D 66, 094014 (2002), eprint hep-ph/0208141.
  * Dominguez et al. (2011) F. Dominguez, C. Marquet, B.-W. Xiao, and F. Yuan, Phys. Rev. D 83, 105005 (2011), eprint 1101.0715.
  * Metz and Zhou (2011) A. Metz and J. Zhou, Phys. Rev. D 84, 051503 (2011), eprint 1105.1991.
  * Dominguez et al. (2012) F. Dominguez, J.-W. Qiu, B.-W. Xiao, and F. Yuan, Phys. Rev. D 85, 045003 (2012), eprint 1109.6293.
  * Iancu and Laidet (2013) E. Iancu and J. Laidet, Nucl. Phys. A 916, 48 (2013), eprint 1305.5926.
  * Altinoluk et al. (2016a) T. Altinoluk, N. Armesto, G. Beuf, and A. H. Rezaeian, Phys. Lett. B 758, 373 (2016a), eprint 1511.07452.
  * Hatta et al. (2016) Y. Hatta, B.-W. Xiao, and F. Yuan, Phys. Rev. Lett. 116, 202301 (2016), eprint 1601.01585.
  * Dumitru et al. (2015) A. Dumitru, T. Lappi, and V. Skokov, Phys. Rev. Lett. 115, 252301 (2015), eprint 1508.04438.
  * Kotko et al. (2015) P. Kotko, K. Kutak, C. Marquet, E. Petreska, S. Sapeta, and A. van Hameren, JHEP 09, 106 (2015), eprint 1503.03421.
  * Marquet et al. (2016) C. Marquet, E. Petreska, and C. Roiesnel, JHEP 10, 065 (2016), eprint 1608.02577.
  * van Hameren et al. (2016) A. van Hameren, P. Kotko, K. Kutak, C. Marquet, E. Petreska, and S. Sapeta, JHEP 12, 034 (2016), [Erratum: JHEP 02, 158 (2019)], eprint 1607.03121.
  * Marquet et al. (2018) C. Marquet, C. Roiesnel, and P. Taels, Phys. Rev. D 97, 014004 (2018), eprint 1710.05698.
  * Dumitru et al. (2019) A. Dumitru, V. Skokov, and T. Ullrich, Phys. Rev. C 99, 015204 (2019), eprint 1809.02615.
  * Dumitru and Jalilian-Marian (2002a) A. Dumitru and J. Jalilian-Marian, Phys. Lett. B 547, 15 (2002a), eprint hep-ph/0111357.
  * Dumitru and Jalilian-Marian (2002b) A. Dumitru and J. Jalilian-Marian, Phys. Rev. Lett. 89, 022301 (2002b), eprint hep-ph/0204028.
  * Mäntysaari et al. (2019) H. Mäntysaari, N. Mueller, and B. Schenke, Phys. Rev. D 99, 074004 (2019), eprint 1902.05087.
  * Salazar and Schenke (2019) F. Salazar and B. Schenke, Phys. Rev. D 100, 034007 (2019), eprint 1905.03763.
  * Boussarie et al. (2021) R. Boussarie, H. Mäntysaari, F. Salazar, and B. Schenke, JHEP 09, 178 (2021), eprint 2106.11301.
  * Ayala et al. (1996) A. Ayala, J. Jalilian-Marian, L. D. McLerran, and R. Venugopalan, Phys. Rev. D 53, 458 (1996), eprint hep-ph/9508302.
  * Jalilian-Marian et al. (1997a) J. Jalilian-Marian, A. Kovner, L. D. McLerran, and H. Weigert, Phys. Rev. D 55, 5414 (1997a), eprint hep-ph/9606337.
  * Kotko et al. (2017) P. Kotko, K. Kutak, S. Sapeta, A. M. Stasto, and M. Strikman, Eur. Phys. J. C 77, 353 (2017), eprint 1702.03063.
  * Hagiwara et al. (2017) Y. Hagiwara, Y. Hatta, R. Pasechnik, M. Tasevsky, and O. Teryaev, Phys. Rev. D 96, 034009 (2017), eprint 1706.01765.
  * Henley and Jalilian-Marian (2006) E. M. Henley and J. Jalilian-Marian, Phys. Rev. D 73, 094004 (2006), eprint hep-ph/0512220.
  * Klein and Mäntysaari (2019) S. R. Klein and H. Mäntysaari, Nature Rev. Phys. 1, 662 (2019), eprint 1910.10858.
  * Hatta et al. (2021b) Y. Hatta, B.-W. Xiao, F. Yuan, and J. Zhou, Phys. Rev. D 104, 054037 (2021b), eprint 2106.05307.
  * Kolbé et al. (2021) I. Kolbé, K. Roy, F. Salazar, B. Schenke, and R. Venugopalan, JHEP 01, 052 (2021), eprint 2008.04372.
  * Altinoluk et al. (2019) T. Altinoluk, R. Boussarie, and P. Kotko, JHEP 05, 156 (2019), eprint 1901.01175.
  * Boussarie et al. (2019) R. Boussarie, A. V. Grabovsky, L. Szymanowski, and S. Wallon, Phys. Rev. D 100, 074020 (2019), eprint 1905.07371.
  * Boussarie et al. (2016) R. Boussarie, A. V. Grabovsky, L. Szymanowski, and S. Wallon, JHEP 11, 149 (2016), eprint 1606.00419.
  * Boussarie et al. (2014) R. Boussarie, A. V. Grabovsky, L. Szymanowski, and S. Wallon, JHEP 09, 026 (2014), eprint 1405.7676.
  * Dumitru and Jalilian-Marian (2010) A. Dumitru and J. Jalilian-Marian, Phys. Rev. D 82, 074023 (2010), eprint 1008.0480.
  * Fadin and Lipatov (1998) V. S. Fadin and L. N. Lipatov, Phys. Lett. B 429, 127 (1998), eprint hep-ph/9802290.
  * Chirilli et al. (2012a) G. A. Chirilli, B.-W. Xiao, and F. Yuan, Phys. Rev. Lett. 108, 122301 (2012a), eprint 1112.1061.
  * Chirilli et al. (2012b) G. A. Chirilli, B.-W. Xiao, and F. Yuan, Phys. Rev. D 86, 054005 (2012b), eprint 1203.6139.
  * Balitsky and Chirilli (2013a) I. Balitsky and G. A. Chirilli, Phys. Rev. D 87, 014013 (2013a), eprint 1207.3844.
  * Balitsky and Chirilli (2013b) I. Balitsky and G. A. Chirilli, Phys. Rev. D 88, 111501 (2013b), eprint 1309.7644.
  * Grabovsky (2013) A. V. Grabovsky, JHEP 09, 141 (2013), eprint 1307.5414.
  * Caron-Huot (2015) S. Caron-Huot, JHEP 05, 093 (2015), eprint 1309.6521.
  * Kovner et al. (2014) A. Kovner, M. Lublinsky, and Y. Mulian, Phys. Rev. D 89, 061704 (2014), eprint 1310.0378.
  * Lublinsky and Mulian (2017) M. Lublinsky and Y. Mulian, JHEP 05, 097 (2017), eprint 1610.03453.
  * Caron-Huot and Herranen (2018) S. Caron-Huot and M. Herranen, JHEP 02, 058 (2018), eprint 1604.07417.
  * Boussarie et al. (2018) R. Boussarie, A. V. Grabovsky, D. Y. Ivanov, L. Szymanowski, and S. Wallon, PoS DIS2017, 062 (2018), eprint 1709.04422.
  * Beuf et al. (2022a) G. Beuf, T. Lappi, and R. Paatelainen, Phys. Rev. D 106, 034013 (2022a), eprint 2204.02486.
  * Beuf et al. (2022b) G. Beuf, T. Lappi, and R. Paatelainen, Phys. Rev. Lett. 129, 072001 (2022b), eprint 2112.03158.
  * Beuf et al. (2021) G. Beuf, T. Lappi, and R. Paatelainen, Phys. Rev. D 104, 056032 (2021), eprint 2103.14549.
  * Mäntysaari and Penttala (2021) H. Mäntysaari and J. Penttala, Phys. Lett. B 823, 136723 (2021), eprint 2104.02349.
  * Mäntysaari and Penttala (2022a) H. Mäntysaari and J. Penttala, Phys. Rev. D 105, 114038 (2022a), eprint 2203.16911.
  * Mäntysaari and Penttala (2022b) H. Mäntysaari and J. Penttala, JHEP 08, 247 (2022b), eprint 2204.14031.
  * Lappi et al. (2022) T. Lappi, H. Mäntysaari, and J. Penttala, SciPost Phys. Proc. 8, 133 (2022), eprint 2106.12825.
  * Iancu and Mulian (2021) E. Iancu and Y. Mulian, JHEP 03, 005 (2021), eprint 2009.11930.
  * Roy and Venugopalan (2020) K. Roy and R. Venugopalan, Phys. Rev. D 101, 034028 (2020), eprint 1911.04530.
  * Hatta et al. (2022) Y. Hatta, B.-W. Xiao, and F. Yuan (2022), eprint 2205.08060.
  * Iancu et al. (2022) E. Iancu, A. H. Mueller, and D. N. Triantafyllopoulos, Phys. Rev. Lett. 128, 202001 (2022), eprint 2112.06353.
  * Taels et al. (2022) P. Taels, T. Altinoluk, G. Beuf, and C. Marquet (2022), eprint 2204.11650.
  * Caucal et al. (2021) P. Caucal, F. Salazar, and R. Venugopalan, JHEP 11, 222 (2021), eprint 2108.06347.
  * Bergabo and Jalilian-Marian (2022a) F. Bergabo and J. Jalilian-Marian, Nucl. Phys. A 1018, 122358 (2022a), eprint 2108.10428.
  * Bergabo and Jalilian-Marian (2022b) F. Bergabo and J. Jalilian-Marian, Phys. Rev. D 106, 054035 (2022b), eprint 2207.03606.
  * Kovchegov et al. (2017a) Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, JHEP 10, 198 (2017a), eprint 1706.04236.
  * Cougoulic and Kovchegov (2019) F. Cougoulic and Y. V. Kovchegov, Phys. Rev. D 100, 114020 (2019), eprint 1910.04268.
  * Kovchegov and Sievert (2019) Y. V. Kovchegov and M. D. Sievert, Phys. Rev. D 99, 054032 (2019), eprint 1808.09010.
  * Kovchegov et al. (2017b) Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, Phys. Lett. B 772, 136 (2017b), eprint 1703.05809.
  * Kovchegov et al. (2017c) Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, Phys. Rev. D 95, 014033 (2017c), eprint 1610.06197.
  * Kovchegov et al. (2017d) Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, Phys. Rev. Lett. 118, 052001 (2017d), eprint 1610.06188.
  * Kovchegov et al. (2016) Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, JHEP 01, 072 (2016), [Erratum: JHEP 10, 148 (2016)], eprint 1511.06737.
  * Agostini et al. (2019a) P. Agostini, T. Altinoluk, and N. Armesto, Eur. Phys. J. C 79, 790 (2019a), eprint 1907.03668.
  * Agostini et al. (2019b) P. Agostini, T. Altinoluk, and N. Armesto, Eur. Phys. J. C 79, 600 (2019b), eprint 1902.04483.
  * Altinoluk and Dumitru (2016) T. Altinoluk and A. Dumitru, Phys. Rev. D 94, 074032 (2016), eprint 1512.00279.
  * Altinoluk et al. (2016b) T. Altinoluk, N. Armesto, G. Beuf, and A. Moscoso, JHEP 01, 114 (2016b), eprint 1505.01400.
  * Altinoluk et al. (2014) T. Altinoluk, N. Armesto, G. Beuf, M. Martínez, and C. A. Salgado, JHEP 07, 068 (2014), eprint 1404.2219.
  * Jalilian-Marian (2021) J. Jalilian-Marian, Nucl. Phys. A 1005, 121943 (2021).
  * Jalilian-Marian (2020) J. Jalilian-Marian, Phys. Rev. D 102, 014008 (2020), eprint 1912.08878.
  * Jalilian-Marian (2019) J. Jalilian-Marian, Phys. Rev. D 99, 014043 (2019), eprint 1809.04625.
  * Jalilian-Marian (2017) J. Jalilian-Marian, Phys. Rev. D 96, 074020 (2017), eprint 1708.07533.
  * Hentschinski et al. (2018) M. Hentschinski, A. Kusina, K. Kutak, and M. Serino, Eur. Phys. J. C 78, 174 (2018), eprint 1711.04587.
  * Hentschinski et al. (2016) M. Hentschinski, A. Kusina, and K. Kutak, Phys. Rev. D 94, 114013 (2016), eprint 1607.01507.
  * Gituliar et al. (2016) O. Gituliar, M. Hentschinski, and K. Kutak, JHEP 01, 181 (2016), eprint 1511.08439.
  * Balitsky and Tarasov (2016) I. Balitsky and A. Tarasov, JHEP 06, 164 (2016), eprint 1603.06548.
  * Balitsky and Tarasov (2015) I. Balitsky and A. Tarasov, JHEP 10, 017 (2015), eprint 1505.02151.
  * Marquet et al. (2009) C. Marquet, B.-W. Xiao, and F. Yuan, Phys. Lett. B 682, 207 (2009), eprint 0906.1454.
  * Kovchegov and Tuchin (2002) Y. V. Kovchegov and K. Tuchin, Phys. Rev. D 65, 074026 (2002), eprint hep-ph/0111362.
  * Ayala et al. (2016) A. Ayala, M. Hentschinski, J. Jalilian-Marian, and M. E. Tejeda-Yeomans, Phys. Lett. B 761, 229 (2016), eprint 1604.08526.
  * Ayala et al. (2017) A. Ayala, M. Hentschinski, J. Jalilian-Marian, and M. E. Tejeda-Yeomans, Nucl. Phys. B 920, 232 (2017), eprint 1701.07143.
  * Balitsky (1996) I. Balitsky, Nucl. Phys. B 463, 99 (1996), eprint hep-ph/9509348.
  * Kovchegov (2000) Y. V. Kovchegov, Phys. Rev. D 61, 074018 (2000), eprint hep-ph/9905214.
  * Jalilian-Marian et al. (1997b) J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Nucl. Phys. B 504, 415 (1997b), eprint hep-ph/9701284.
  * Jalilian-Marian et al. (1998a) J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Phys. Rev. D 59, 014014 (1998a), eprint hep-ph/9706377.
  * Jalilian-Marian et al. (1998b) J. Jalilian-Marian, A. Kovner, and H. Weigert, Phys. Rev. D 59, 014015 (1998b), eprint hep-ph/9709432.


# GEVD-based Low-Rank Channel Covariance Matrix Estimation and MMSE Channel
Estimation for Uplink Cellular Massive MIMO Systems

###### Abstract

Uplink channel estimation is a crucial component for the performance of
cellular massive MIMO systems. However, when the number of user equipments
(UEs) grows, the sharing of the available resources causes interference
between UEs in neighboring cells. Minimum mean squared error (MMSE) channel
estimators have been proposed to mitigate this interference, but these require
estimates of the channel covariance matrices. Therefore, a new channel
covariance matrix estimator for low-rank channel covariance matrices is
presented in this paper, using a generalized eigenvalue decomposition (GEVD)
of two covariance matrices that can be estimated from the available uplink
data. The requirements for the systems are minimal and, except for
synchronization, there is no need for communication between the different
cells and no prior knowledge on the background noise is required. Approximate
MMSE estimators are also derived based on the newly proposed channel
covariance matrix estimator. The effectiveness of the proposed methods is
demonstrated in numerical simulations.

Index Terms—  Covariance matrix estimation, random pilot allocation, low-rank
channel model, generalized eigenvalue decomposition (GEVD), massive MIMO

##  1 Introduction

Massive MIMO (multiple input multiple output) is a key technology for current
generation cellular networks [1, 2, 3]. Base stations (BSs) are deploying
hundreds of antennas, enabling spatial multiplexing of several user equipments
(UEs) per cell, both in uplink and downlink, with simple linear signal
processing methods [4]. The increasing number of UEs, and as a consequence the
sharing of the available pilots whenever the coherence blocks are limited in
size, causes interference between UEs in neighboring cells. To make it
possible for a BS to separate the UEs that it serves from the interference,
second-order statistical information present in the channel covariance
matrices has to be exploited during the channel estimation.

The channel covariance matrices are commonly assumed to be perfectly known in
massive MIMO literature, which is a strong assumption since the matrix
dimension grows with the number of antennas and furthermore the statistics are
changing over time. Existing covariance matrix estimators either use
additional pilot overhead, sacrifice samples for data transmission [5], use
strong separability conditions that might not always hold [6] or require
knowledge of the used pilots of all the UEs in the network [7, 8]. Therefore,
a new channel covariance matrix estimator for low-rank channel covariance
matrices is presented in this paper, using a generalized eigenvalue
decomposition (GEVD) of two covariance matrices that can be estimated from the
available uplink data. The requirements for the system are minimal and, except
for synchronization, there is no need for communication between the different
cells and no prior knowledge on the background noise is required. Approximate
MMSE estimators are also derived based on the newly proposed channel
covariance matrix estimator. The effectiveness of the proposed methods is
demonstrated in numerical simulations.

##  2 Channel and Signal Model

###  2.1 Channel Model

A cellular massive MIMO system with $L$ cells is considered where each cell
contains a BS with $N$ antennas and $K$ single-antenna UEs. The channel from
UE $k$ in cell $l$ to the BS in cell $j$ is denoted by
$\mathbf{h}_{jlk}\in\mathbb{C}^{N}$ and is assumed to remain constant during a
coherence block of $\tau_{c}$ samples. Superscript $.^{t}$ will be used to
denote the quantities in coherence block $t$, e.g. $\mathbf{h}^{t}_{jkl}$. A
common channel model for $\mathbf{h}^{t}_{jkl}$ is that it is drawn from a
correlated Rayleigh fading realization
$\mathcal{NC}(\mathbf{0},\overline{\mathbf{R}}_{jlk})$, where
$\overline{\mathbf{R}}_{jlk}\in\mathbb{C}^{N\times N}$ is the positive semi-
definite channel covariance matrix describing the large-scale fading,
including geometric pathloss, shadowing, antenna gains, and spatial channel
correlation [9]. The complex Gaussian distribution around
$\overline{\mathbf{R}}_{jlk}$ models the small-scale fading. Due to the
different positions of the UEs, it is assumed that the channels for different
UEs are uncorrelated, i.e.
$E\\{\mathbf{h}^{t}_{jlk}\mathbf{h}^{t,H}_{jmi}\\}=\mathbf{0}$ if
$\\{l,k\\}\neq\\{m,i\\}$, where $E\\{.\\}$ denotes the expected value
operator, with respect to different channel realizations and $.^{H}$ is the
Hermitian transpose operator. Measurements in [10] show that the matrices
$\\{\overline{\mathbf{R}}_{jlk}\\}$ remain constant over several thousands of
coherence blocks, where this large number is denoted by $T$.

The scattering is mostly localized around the UE, since the BSs are elevated
and hence have limited scattering in their near-field [8]. The multipath
components defining the total channel $\mathbf{h}_{jlk}$ are therefor arriving
from a specific localized region in the network, so that the matrices
$\\{\overline{\mathbf{R}}_{jlk}\\}$ can be approximated by low-rank matrices,
where the rank is denoted by $R$. $\overline{\mathbf{R}}_{jlk}$ can thus be
approximated as a sum of $R$ rank-1 terms:

$\overline{\mathbf{R}}_{jlk}=\sum_{r=1}^{R}\overline{\mathbf{q}}_{jlk,r}\overline{\mathbf{q}}_{jlk,r}^{H}$
(1)

where $\overline{\mathbf{q}}_{jlk,r}$ will be defined later. This implies that
the channel $\mathbf{h}^{t}_{jlk}$ in coherence block $t$ can be modeled as

$\mathbf{h}^{t}_{jlk}=\sum_{r=1}^{R}\overline{\mathbf{q}}_{jlk,r}\overline{z}^{t}_{jlk,r}=\overline{\bm{\mathds{Q}}}_{jlk}\overline{\mathbf{z}}_{jlk}^{t}$
(2)

where the normal complex variables
$\overline{\mathbf{z}}_{jlk}^{t}=[\overline{z}^{t}_{jlk,1}\ ...\
\overline{z}^{t}_{jlk,R}]^{T}$ with $.^{T}$ the transpose operator, are
independently drawn from $\mathcal{NC}(0,1)$ and
$\overline{\bm{\mathds{Q}}}_{jlk}=[\overline{\mathbf{q}}_{jlk,1}\ ...\
\overline{\mathbf{q}}_{jlk,R}]$.

###  2.2 Signal Model

The samples in one coherence block of $\tau_{c}$ samples are divided in
$\tau_{p}$ samples for uplink channel estimation and $\tau_{u}$ samples for
uplink data transmission with $\tau_{c}=\tau_{p}+\tau_{u}$. With $\tau_{p}$
samples available, the network can predefine $\tau_{p}$ different orthogonal
and unitary pilot sequences $\\{\mathbf{s}_{b}\\}$ where

$\mathbf{s}_{b}=[s_{b}(1)\ ...\ s_{b}(\tau_{p})]^{T},\ b=1...\tau_{p}$ (3)

with $\mathbf{s}^{H}_{b}\mathbf{s}_{b}=\tau_{p}$ and
$\mathbf{s}^{H}_{b}\mathbf{s}_{c}=0$ if $b\neq c$. In the uplink channel
estimation phase, each UE $k$ in each cell $l$ selects randomly one of the
$\tau_{p}$ pilot sequences in each coherence block $t$ , the index of this
pilot sequence is denoted by $b^{t}_{lk}$. An example is provided in Table 1.
The UEs transmits the chosen pilot sequence with powers $\\{p_{lk}\\}$. The
signal received at BS $j$ is then given as

$\mathbf{y}^{t}_{j}(p)=\sum_{l=1}^{L}\sum_{k=1}^{K}\sqrt{p_{lk}}\mathbf{h}^{t}_{jlk}s_{b_{lk}^{t}}(p)+\mathbf{n}^{t}_{j}(p),\
p=1...\tau_{p}$ (4)

where $p$ is the sample index and $\mathbf{n}^{t}_{j}(p)$ is the background
noise at BS $j$.

|  | Coherence block
---|---|---
|  |  |  |  |  |
|  | 1 | 2 | 3 | … | $T$
|  | $\mathbf{s}_{1}$ | $\mathbf{s}_{2}$ | $\mathbf{s}_{3}$ | $\mathbf{s}_{4}$ | $\mathbf{s}_{1}$ | $\mathbf{s}_{2}$ | $\mathbf{s}_{3}$ | $\mathbf{s}_{4}$ | $\mathbf{s}_{1}$ | $\mathbf{s}_{2}$ | $\mathbf{s}_{3}$ | $\mathbf{s}_{4}$ | … | $\mathbf{s}_{1}$ | $\mathbf{s}_{2}$ | $\mathbf{s}_{3}$ | $\mathbf{s}_{4}$
UE | 1 |  |  |  |  |  |  |  |  |  |  |  |  |  | … |  |  |  |  |
2 |  |  |  |  |  |  |  |  |  |  |  |  | … |  |  |  |
3 |  |  |  |  |  |  |  |  |  |  |  |  | … |  |  |  |
4 |  |  |  |  |  |  |  |  |  |  |  |  | … |  |  |  |
… |  | … |  |  | … |  |  | … |  | … |  | … |
| $K$ |  |  |  |  |  |  |  |  |  |  |  |  | … |  |  |  |

Table 1: Example of random allocation of 4 pilot sequences over $K$ UEs in a
certain cell.

In the uplink data transmission phase, each UE $k$ in each cell $l$ transmits
the unitary signal $s_{lk}(u)$ with zero mean and power $p_{lk}$. The received
signal at BS $j$ is then given as

$\mathbf{y}^{t}_{j}(u)=\sum_{l=1}^{L}\sum_{k=1}^{K}\sqrt{p_{lk}}\mathbf{h}^{t}_{jlk}s_{lk}(u)+\mathbf{n}^{t}_{j}(u),\
u=1...\tau_{u}.$ (5)

Note that the indices $p$ and $u$ are used to distinguish between the
different phases. BS $j$ has to detect the uplink data from the UEs in its
cell, i.e. $s_{jk}(u)$ for $k=1...K$. To perform receive combining strategies
like maximum-ratio (MR), (regularized) zero-forcing (RZF) or minimum mean
squared error (MMSE) combining [9, 11], good estimates of the channels of the
$K$ UEs in cell $j$ are required. The background noise $\mathbf{n}_{j}^{t}$ is
assumed to be zero-mean complex Gaussian, but with an unknown spatial
covariance matrix
$\overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}}=E\\{\mathbf{n}_{j}^{t}\mathbf{n}_{j}^{t,H}\\}$.
It is often assumed to be white noise with a diagonal covariance matrix, but
here a more general noise model is used, where the noise can be correlated
between the antennas of a BS. This happens for example when there are located
jamming signals nearby the BS.

##  3 Optimal MMSE Channel Estimation

The signal obtained from (4) after despreading with pilot sequence
$b^{t}_{jk}$ can be represented as

$\displaystyle\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]=$
$\displaystyle\sum_{u=1}^{\tau_{p}}\mathbf{y}^{t}_{j}(p)s^{*}_{b_{jk}^{t}}(p)$
(6) $\displaystyle=$
$\displaystyle\sum_{l=1}^{L}\sum_{i=1}^{K}\delta^{t}_{jk,li}\sqrt{p_{li}}\tau_{p}\mathbf{h}^{t}_{jli}+\underbrace{\sum_{u=1}^{\tau_{p}}\mathbf{n}^{t}_{j}(p)s^{*}_{b_{jk}^{t}}(p)}_{\mathbf{n}_{j}^{t,\text{pilot}}[b^{t}_{jk}]}$

where the random variable $\delta^{t}_{jk,li}$ is 1 with probability
$1/\tau_{p}$ and 0 otherwise since all other UEs randomly pick 1 of the
$\tau_{p}$ pilot sequences and $\delta^{t}_{jk,jk}$ is deterministic and equal
to 1.

Based on the despread signal in (6) a linear MMSE channel estimation
criterion, similar to [12, Sect. 4.2], can be defined to obtain an estimate
$\overline{\mathbf{h}}^{t}_{jjk}=\overline{\mathbf{W}}^{H}_{jk}\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]$
of $\mathbf{h}_{jjk}^{t}$ in the current coherence block $t$:

$\overline{\mathbf{W}}_{jk}=\operatorname*{arg\,min}_{\bm{W}_{jk}}E\\{||\mathbf{h}^{t}_{jjk}-\bm{W}^{H}_{jk}\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]||^{2}\\}.$
(7)

The optimal solution of (7) is given as

$\overline{\mathbf{W}}_{jk}=\sqrt{p_{jk}}\left(\overline{\mathbf{R}}^{\text{pilot}}_{jjk}\right)^{-1}\overline{\mathbf{R}}_{jjk}$
(8)

where

$\displaystyle\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$
$\displaystyle=\frac{1}{\tau_{p}}E\\{\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]\mathbf{y}_{j}^{t,\text{pilot},H}[b^{t}_{jk}]\\}$
(9) $\displaystyle=p_{jk}\tau_{p}\overline{\mathbf{R}}_{jjk}+\sum_{i\neq
k}p_{ji}\overline{\mathbf{R}}_{jji}+\sum_{l\neq
j}\sum_{i=1}^{K}p_{li}\overline{\mathbf{R}}_{jli}+\overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}}$

since the channels between different UEs are uncorrelated and are also
uncorrelated with the background noise. The fact that
$E\\{{\delta^{t^{2}}_{jk,li}}\\}=\frac{1}{\tau_{p}}$ is also used. For the
background noise term, the following derivation is used:

$\displaystyle E\\{$
$\displaystyle\mathbf{n}_{j}^{t,\text{pilot}}[b^{t}_{jk}]\mathbf{n}_{j}^{t,\text{pilot},H}[b^{t}_{jk}]\\}$
(10)
$\displaystyle=E\left\\{\left(\sum_{u=1}^{\tau_{p}}\mathbf{n}_{j}^{t}(p)s^{*}_{b_{jk}^{t}}(p)\right)\left(\sum_{u=1}^{\tau_{p}}\mathbf{n}_{j}^{t,H}(p)s_{b_{jk}^{t}}(p)\right)^{H}\right\\}$
$\displaystyle=\sum_{u=1}^{\tau_{p}}E\\{|s_{b_{jk}^{t}}(p)|^{2}\\}E\\{\mathbf{n}_{j}^{t}\mathbf{n}_{j}^{t,H}\\}=\tau_{p}\overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}}$

where the cross terms in the second step are zero since the noise is
uncorrelated between different samples.

An important difference with the MMSE channel estimation method proposed in
[12, Sect. 4.1,] is that now all UEs appear as interfering UEs in (9), while
in [12] it is assumed that the pilot sequences for all the UEs are
predetermined and fixed, so that only UEs sharing the same pilot sequence
appear as interfering UEs. However, this assumption requires prior knowledge
of the channel statistics of all the UEs and the background noise, while the
proposed method allows for a more efficient estimation method which will be
discussed next.

##  4 Low-rank Channel Estimation

###  4.1 Low-rank covariance matrix estimator

In practice the covariance matrices
$\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ and $\overline{\mathbf{R}}_{jjk}$
have to be estimated from the available data.
$\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ can be estimated using time
averaging over the last $T$ coherence blocks, possibly using a diagonal
regularizer as in [13, 5], since perfect synchronization and knowledge of
$b_{jk}^{t}$ at the BS $j$ is assumed. The estimation will be denoted as
$\mathbf{R}^{\text{pilot}}_{jjk}$ (without the overline).
$\overline{\mathbf{R}}_{jjk}$ as in (1) can not immediately be estimated,
since the channel is unknown. However, the combined covariance matrix
$\mathbf{R}^{\text{all}}_{j}$ can be estimated from all available antenna
signals during both the channel estimation phase and uplink data transmission
phase, which is an estimate of

$\displaystyle\overline{\mathbf{R}}^{\text{all}}_{j}$
$\displaystyle=E\left\\{\frac{1}{\tau_{p}}\sum_{p=1}^{\tau_{p}}\mathbf{y}_{j}^{t}(p)\mathbf{y}_{j}^{t,H}(p)+\frac{1}{\tau_{u}}\sum_{u=1}^{\tau_{u}}\mathbf{y}_{j}^{t}(u)\mathbf{y}_{j}^{t,H}(u)\right\\}$
$\displaystyle=\sum_{l=1}^{L}\sum_{i=1}^{K}p_{li}\overline{\mathbf{R}}_{jli}+\overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}}.$
(11)

When comparing (9) and (4.1), it can be stated that
$\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ is the covariance matrix
estimated after despreading and $\overline{\mathbf{R}}^{\text{all}}_{j}$ is
the covariance matrix estimated before despreading. They differ only in a
contribution depending on $\overline{\mathbf{R}}_{jjk}$. The following
relation can thus be derived:

$p_{jk}\overline{\mathbf{R}}_{jjk}=\frac{1}{\tau_{p}-1}\left(\overline{\mathbf{R}}^{\text{pilot}}_{jjk}-\overline{\mathbf{R}}^{\text{all}}_{j}\right).$
(12)

However, simply replacing $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ and
$\overline{\mathbf{R}}^{\text{all}}_{j}$ with the estimated quantities
$\mathbf{R}^{\text{pilot}}_{jjk}$ and $\mathbf{R}^{\text{all}}_{j}$ in (12)
will not result in a good estimate $\mathbf{R}_{jjk}$ of
$\overline{\mathbf{R}}_{jjk}$ since it will generally not satisfy the assumed
rank assumption in (1) and might even result in an indefinite estimate
$\mathbf{R}_{jjk}$. One known remedy consists in using a rank-approximation
based on the GEVD [14, 15] of the matrix pencil
$\\{{\mathbf{R}}^{\text{pilot}}_{jjk},{\mathbf{R}}^{\text{all}}_{j}\\}$ given
as

$\displaystyle{\mathbf{R}}^{\text{pilot}}_{jjk}$
$\displaystyle=\mathbf{Q}_{jjk}\bm{\Sigma}_{jjk}\mathbf{Q}^{H}_{jjk}=\sum_{r=1}^{N}\sigma_{jjk,r}\mathbf{q}_{jjk,r}\mathbf{q}^{H}_{jjk,r}$
(13) $\displaystyle{\mathbf{R}}^{\text{all}}_{j}$
$\displaystyle=\mathbf{Q}_{jjk}\mathbf{Q}^{H}_{jjk}=\sum_{r=1}^{N}\mathbf{q}_{jjk,r}\mathbf{q}^{H}_{jjk,r}$
(14)

where $\bm{\Sigma}_{jjk}=\text{diag}\\{\sigma_{jjk,1},...,\sigma_{jjk,N}\\}$
is a diagonal matrix with the generalized eigenvalues sorted from large to
small and the columns of $\mathbf{X}_{jjk}=\mathbf{Q}^{-H}_{jjk}$ contain the
generalized eigenvectors. The eigenvectors are normalized such that
$\mathbf{X}^{H}_{jjk}{\mathbf{R}}^{\text{all}}_{j}\mathbf{X}_{jjk}=\mathbf{I}$.
The optimal rank $R$ estimate $\mathbf{R}_{jjk}$ is given by keeping only the
rank-1 terms belonging to the $R$ largest eigenvalues greater than one in (12)

$p_{jk}\mathbf{R}_{jjk}=\sum_{r=1}^{R}\frac{\sigma_{jjk,r}-1}{\tau_{p}-1}\mathbf{q}_{jjk,r}\mathbf{q}^{H}_{jjk,r}=\mathds{Q}_{jjk}\mathbb{\Lambda}_{jjk}\mathds{Q}^{H}_{jjk}$
(15)

where
$\mathbb{\Lambda}_{jjk}=\text{diag}\\{\frac{\sigma_{jjk,1}-1}{\tau_{p}-1},...,$
$\frac{\sigma_{jjk,R}-1}{\tau_{p}-1}\\}$.

As stated in [14], the GEVD-based approximation effectively selects the modes
with the highest SINR, since approximation errors are weighted relative to the
interference and noise aggregated in ${\mathbf{R}}^{\text{all}}_{j}$. Also
note that unlike other rank-approximations of (12), the GEVD-based
approximation is immune to scaling and linear combining of the signals, i.e.
the output signal and output SINR is independent of such scaling and
combining, which is a desirable property.

The proposed covariance matrix estimator is data-driven and requires only
synchronization and knowledge of the chosen pilot sequence $b^{t}_{jk}$ of the
UEs in cell $j$ in each coherence block. There is no need for prior knowledge
of $\overline{\mathbf{R}}_{jjk}$ (except for the rank $R$ of
$\overline{\mathbf{R}}_{jjk}$),
$\overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}}$, or any prior knowledge
on the statistics or pilot sequences of the UEs in other cells $l\neq j$.

###  4.2 Low-rank approximate MMSE estimator

When the estimated quantities are used to calculate the filter in (7), the
following filter is obtained:

$\displaystyle\mathbf{W}_{jk}$
$\displaystyle=\sqrt{p_{jk}}\left(\mathbf{Q}_{jjk}\bm{\Sigma}_{jjk}\mathbf{Q}^{H}_{jjk}\right)^{-1}\frac{1}{p_{jk}}\mathds{Q}_{jjk}\mathbb{\Lambda}_{jjk}\mathds{Q}^{H}_{jjk}$
(16)
$\displaystyle=\frac{1}{\sqrt{p_{jk}}}\mathds{X}_{jjk}\mathbb{\Upsilon}_{jjk}\mathds{Q}^{H}_{jjk}$

with
$\mathbb{\Upsilon}_{jjk}=\text{diag}\\{\upsilon_{jjk,1},...,\upsilon_{jjk,R}\\}=\frac{1}{\tau_{p}-1}\text{diag}\\{\frac{\sigma_{jjk,1}-1}{\sigma_{jjk,1}},...,$
$\frac{\sigma_{jjk,R}-1}{\sigma_{jjk,R}}\\}$ and
$\mathds{X}_{jjk}=[\mathbf{x}_{jjk,1}\ ...\ \mathbf{x}_{jjk,R}]$ contains the
first $R$ columns of $\mathbf{X}_{jjk}$. The final low-rank approximate MMSE
estimate for the channel of UE $k$ in coherence block $t$ is then given as

$\displaystyle\hat{\mathbf{h}}^{t}_{jjk}$
$\displaystyle=\mathbf{W}^{H}_{jk}\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]=\frac{1}{\sqrt{p_{jk}}}\mathds{Q}_{jjk}\mathbb{\Upsilon}_{jjk}\mathds{X}^{H}_{jjk}\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]$
$\displaystyle=\frac{1}{\sqrt{p_{jk}}}\sum_{r=1}^{R}\mathbf{q}_{jjk,r}\underbrace{\upsilon_{jjk,r}\left(\mathbf{x}^{H}_{jjk,r}\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]\right)}_{z^{t}_{jjk,r}}.$
(17)

Comparing this expression for $\hat{\mathbf{h}}^{t}_{jjk}$ with the model used
for ${\mathbf{h}}^{t}_{jjk}$ in (2), it can be observed that the estimation
method produces an estimate that exactly fits in this model. The random
variable $\bar{\mathbf{z}}^{t}_{jjk}$ of (2) is estimated from the despread
signal $\mathbf{y}_{j}^{t,\text{pilot}}[b^{t}_{jk}]$ using a simple inner-
product with the corresponding eigenvector $\mathbf{x}_{jjk,r}$ and some
scalars depending on the generalized eigenvalues.

###  4.3 Improved low-rank approximate MMSE estimator

An improved low-rank approximate MMSE estimator can be derived based on (6).
When computing the expectation in (7), the activity of all the UEs is
represented by the random variable $\delta^{t}_{jk,li}$. Since BS $j$ is aware
of the used pilot sequences of the UEs in its cell, this information can be
used to improve the MMSE estimator:

$\displaystyle\mathbf{y}_{j}^{t,\text{pilot,impr}}[b^{t}_{jk}]=$
$\displaystyle\sum_{i=1}^{K}\bar{\delta}^{t}_{jk,ji}\sqrt{p_{ji}}\tau_{p}\mathbf{h}^{t}_{jji}$
(18) $\displaystyle+\sum_{l\neq
j}\sum_{i=1}^{K}\delta^{t}_{jk,li}\sqrt{p_{li}}\tau_{p}\mathbf{h}^{t}_{jli}+\mathbf{n}_{j}^{t,\text{pilot}}[b^{t}_{jk}]$

where $\bar{\delta}^{t}_{jk,ji}$ is now deterministic and either equal to 0 or
1 and $\bar{\delta}^{t}_{jk,jk}=1$. The corresponding cell-specific covariance
matrix that now depends on the coherence block index $t$ is given as
(alternative to (9))

$\displaystyle\overline{\mathbf{R}}^{\text{pilot,impr,t}}_{jjk}$
$\displaystyle=\frac{1}{\tau_{p}}E\\{\mathbf{y}_{j}^{t,\text{pilot,impr}}[b^{t}_{jk}]\mathbf{y}_{j}^{t,\text{pilot,impr},H}[b^{t}_{jk}]\\}$
(19)
$\displaystyle=\sum_{i=1}^{K}\bar{\delta}^{t^{2}}_{jk,ji}p_{ji}\tau_{p}\overline{\mathbf{R}}_{jji}+\sum_{l\neq
j}\sum_{i=1}^{K}p_{li}\overline{\mathbf{R}}_{jli}+\overline{\mathbf{R}}_{\mathbf{n}_{j}\mathbf{n}_{j}}$
$\displaystyle=\overline{\mathbf{R}}^{\text{pilot}}_{jjk}+\sum_{i\neq
k}(\bar{\delta}^{t^{2}}_{jk,ji}\tau_{p}-1)p_{ji}\overline{\mathbf{R}}_{jji}.$

where $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$ is defined as in (9). By
replacing the covariance matrices $\overline{\mathbf{R}}^{\text{pilot}}_{jjk}$
and $\overline{\mathbf{R}}_{jji}$ with their estimated quantities as explained
in Subsection 4.1, an improved low-rank approximate MMSE estimate
$\hat{\mathbf{h}}^{t,\text{impr}}_{jjk}=\mathbf{W}^{\text{impr},t,H}_{jk}\mathbf{y}_{j}^{t,\text{pilot,impr}}[b^{t}_{jk}]$
can be obtained where

$\mathbf{W}^{\text{impr},t}_{jk}=\Big{(}\mathbf{R}^{\text{pilot}}_{jjk}+\sum_{i\neq
k}(\bar{\delta}^{t^{2}}_{jk,ji}\tau_{p}-1)p_{ji}\mathbf{R}_{jji}\Big{)}^{-1}\mathbf{R}_{jjk}.$
(20)

Although this estimator uses more of the information available at BS $j$, it
is important to note that the provided improved low-rank approximate MMSE
estimator requires an extra matrix inversion for each UE in its cell, while
this is not required for the low-rank approximate MMSE estimator of Subsection
4.2.

##  5 Numerical Simulations

The same multicell setup with $L=7$ hexagonal cells as in [5, 7] is
considered. There are $K=10$ UEs per cell in a ring around the BS (there are
70 UEs in total) and each BS has a uniform linear array with $N=100$ antennas.
The multipath components from a UE arrive uniformly distributed in an angular
interval ($10^{\circ}$) centered around the geographical angle to the UE [16].
This results in channel covariance matrices with approximately 25 dominant
eigenvalues.

The normalized MSE (NMSE) (averaged over 20 Monte-Carlo runs),
$E\\{||\mathbf{h}_{jjk}-\hat{\mathbf{h}}_{jjk}||^{2}\\}/\text{tr}\\{\overline{\mathbf{R}}_{jjk}\\}$,
for an average UE in the center cell is used as a performance metric. The
proposed low-rank approximate MMSE estimator is calculated for different
ranks: $R=30$ and $R=60$, denoted with GEVD 30 and GEVD 60 respectively, and
also when the subtraction in (12) is used instead of the GEVD to estimate the
channel covariance matrices, denoted with SUBT. The performance of the MMSE
estimator with perfect channel covariance matrix knowledge (MMSE random) in
(8) is also shown. Similar results are shown for the proposed improved
approximate MMSE estimator, denoted with the word impr. As baseline, the
simple least-square estimation (LS fixed) that does not require covariance
matrix information is shown in a scenario where each UE is allocated with a
fixed pilot sequence for all coherence blocks [5]. Also the MMSE estimate with
perfect channel covariance matrix knowledge (MMSE fixed) for the fixed pilot
allocation regime is provided as absolute lower bound.

Figure 1 shows that a relative small number of coherence blocks (300) suffices
to obtain an accurate channel estimate for the proposed methods. The
estimators with a rank constraint on the channel covariance matrix outperform
the estimators without the rank constraint (SUBT) and $R=30$ works better for
the current setup than $R=60$, since it can already capture the dominant
eigenspace. The proposed GEVD-based procedure focuses on the modes with the
highest power in the estimated covariance matrices and neglect the other
modes, that typically contain more estimation error noise than useful signal
power. The improved approximate MMSE estimator is also outperforming the
approximate MMSE estimator when there are enough coherence blocks available
for the estimation procedure. Since all the estimates of the channel
covariance matrices are used to calculate the filter in (20), this will only
perform well when they are estimated using enough samples. Although not shown
in Figure 1, it should be noted that the SUBT-curve will converge to the MMSE
fixed-curve when $T\rightarrow\infty$ while the GEVD-curves perform worse,
since the rank assumption on the channel covariance matrices is only
approximately valid for this setup.

In Figure 2, the performance with respect to the number of available pilot
sequences $\tau_{p}$ is shown for $T=1500$ coherence blocks. The seemingly
arbitrary shapes of the curves with a fixed pilot allocation emerge from the
fact that the pilots are allocated in a cyclic way, starting from the leftmost
UE in a cell and the cycle is not restarted in the neighboring cell. It is
again observed that the proposed low-rank estimators outperform the estimates
obtained with subtraction of the covariance matrices and the baseline LS-
method. GEVD 30 provides the best performance, independently of the number of
available pilot sequences. Increasing $\tau_{p}$, increases the dominance of
$\overline{\mathbf{R}}_{jjk}$ in (9), and has thus a positive impact on the
channel estimation performance.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0$$400$$800$$1{,}200$$1{,}600$$2{,}000$$2{,}400$$2{,}800$$T$NMSEMMSE
fixedMMSE randomLS fixedSUBTSUBT imprGEVD 30GEVD 30 imprGEVD 60GEVD 60 impr

Fig. 1: NMSE for an average intra-cell UE for different estimators with
$\tau_{p}=10$ when the number of coherence block is varied.

$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$\tau_{p}$NMSEMMSE
fixedMMSE randomLS fixedSUBTSUBT imprGEVD 30GEVD 30 imprGEVD 60GEVD 60 impr

Fig. 2: NMSE for an average intra-cell UE for different estimators with
$T=1500$ when the number of pilot sequences is varied.

##  6 Conclusion

In this paper, a new uplink channel covariance matrix estimator has been
presented for low-rank channel covariance matrices, using a GEVD of two
covariance matrices that are estimated from the available uplink data. The
requirements for the systems are minimal, i.e. UEs in a cell can randomly
choose a pilot sequence in each coherence block and only the BS in its cell
should be aware of its choice. Except for synchronization, there is no need
for communication between the different BSs and no prior knowledge on the
background noise is required. The derived approximate MMSE and improved
approximate MMSE estimators have been shown to provide a good estimate of the
true channel, even when only a small number of coherence blocks is available,
while a simple subtraction of the two covariance matrices provides inferior
performance due to estimation errors.

## References

  * [1] Jakob Hoydis, Stephan Ten Brink, and Merouane Debbah,  “Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?,”  IEEE Journal on Selected Areas in Communications, vol. 31, no. 2, pp. 160–171, 2013.
  * [2] Emil Björnson, Jakob Hoydis, and Luca Sanguinetti,  “Massive MIMO Has Unlimited Capacity,”  IEEE Transactions on Wireless Communications, vol. 17, no. 1, pp. 574–590, 2018.
  * [3] Emil Björnson, Luca Sanguinetti, Henk Wymeersch, Jakob Hoydis, and Thomas L. Marzetta,  “Massive MIMO is a reality—What is next?: Five promising research directions for antenna arrays,”  Digital Signal Processing, vol. 94, pp. 3–20, 2019.
  * [4] Fredrik Rusek, Daniel Persson, Buon Kiong Lau, Erik G. Larsson, Thomas L. Marzetta, Ove Edfors, and Fredrik Tufvesson,  “Scaling up MIMO : Opportunities and challenges with very large arrays,”  IEEE Signal Processing Magazine, vol. 30, no. 1, pp. 40–60, 2013\.
  * [5] Emil Björnson, Luca Sanguinetti, and Merouane Debbah,  “Massive MIMO with imperfect channel covariance information,”  2016 50th Asilomar Conference on Signals, Systems and Computers, pp. 974–978, 2016.
  * [6] Saeid Haghighatshoar and Giuseppe Caire,  “Massive MIMO pilot decontamination and channel interpolation via wideband sparse channel estimation,”  IEEE Transactions on Wireless Communications, vol. 16, no. 12, pp. 8316–8332, 2017.
  * [7] David Neumann, Michael Joham, and Wolfgang Utschick,  “Covariance Matrix Estimation in Massive MIMO,”  IEEE Signal Processing Letters, vol. 25, no. 6, pp. 863–867, 2018\.
  * [8] Hongxiang Xie, Feifei Gao, and Shi Jin,  “An Overview of Low-Rank Channel Estimation for Massive MIMO Systems,”  IEEE Access, vol. 4, pp. 7313–7321, 2016.
  * [9] Emil Björnson, Jakob Hoydis, and Luca Sanguinetti,  “Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency,”  Foundations and Trends® in Signal Processing, vol. 11, no. 3-4, pp. 154–655, 2017.
  * [10] Ingo Viering, Helmut Hofstetter, and Wolfgang Utschick,  “Spatial Long-term Variantions in Urban, Rurla and Indoor Enviroments,”  in in the 5th Meeting of COST273, Lisbon, Portugal, 2002.
  * [11] Özgecan Özdogan, Emil Björnson, and Erik G. Larsson,  “Massive MIMO with Spatially Correlated Rician Fading Channels,”  IEEE Transactions on Communications, vol. 67, no. 5, pp. 3234–3250, 2019.
  * [12] Özlem Tugfe Demir, Emil Björnson, and Luca Sanguinetti,  “Foundations of User-Centric Cell-Free Massive MIMO,”  Foundations and Trends® in Signal Processing, vol. 14, no. 3-4, pp. 162–472, 2021.
  * [13] Olivier Ledoit and Michael Wolf,  “A well-conditioned estimator for large-dimensional covariance matrices,”  Journal of Multivariate Analysis, vol. 88, no. 2, pp. 365–411, 2004\.
  * [14] Romain Serizel, Marc Moonen, Bas Van Dijk, and Jan Wouters,  “Low-rank Approximation Based Multichannel Wiener Filter Algorithms for Noise Reduction with Application in Cochlear Implants,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 22, no. 4, pp. 785–799, apr 2014.
  * [15] M. Dendrinos, S. Bakamidis, and G. Carayannis,  “Speech enhancement from noise: A regenerative approach,”  Speech Commun., vol. 10, no. 1, pp. 45–57, feb 1991.
  * [16] Ansuman Adhikary, Junyoung Nam, Jae Young Ahn, and Giuseppe Caire,  “Joint spatial division and multiplexing-The large-scale array regime,”  IEEE Transactions on Information Theory, vol. 59, no. 10, pp. 6441–6463, 2013.