# Next-to-eikonal corrections to dijet production in Deep Inelastic Scattering
in the dilute limit of the Color Glass Condensate

Pedro Agostinia,b, Tolga Altinoluka and Néstor Armestob a Theoretical Physics
Division, National Centre for Nuclear Research, Pasteura 7, Warsaw 02-093,
Poland
b Instituto Galego de Física de Altas Enerxías IGFAE, Universidade de Santiago
de Compostela, 15782 Santiago de Compostela, Galicia-Spain

###### Abstract

We analyze the effects of next-to-eikonal corrections on dijet production in
Deep Inelastic Scattering off nuclear targets in the framework of the Color
Glass Condensate. They require the knowledge of correlators of fields in the
target beyond those computed in the standard McLerran-Venugopalan model,
specifically those between transverse and boost-enhanced components, and of
the recoil of the fields. We neglect the latter, while for the former we
develop a linear model valid for large nuclei. We considered the unpolarized
cross sections for dijet production in the approximation of a homogenous
dilute nucleus, obtaining simple analytic expressions for the cross sections
at nex-to-eikonal accuracy, valid in the limit of total dijet momentum and
dijet momentum imbalance larger than the saturation scale of the nucleus. We
perform a numerical study of the results at energies of the Electron Ion
Collider, finding $\mathcal{O}(10\%)$ effects in the cross sections at large
total momentum. We also analyze the azimuthal asymmetries between total
momentum and imbalance, finding that non-eikonal corrections induce odd
azimuthal harmonics for the situation of jets with equal momentum fractions
from the virtual photon, where they are absent in the eikonal approximation.
Finally, in the eikonal approximation we have compared the results of our
analytic expansion valid in the dilute limit of the target, and the full Color
Glass Condensate results in the McLerran-Venugopalan model and their
correlation limit. Our analytic expressions match the correlation limit ones
in the region where both should be simultaneously valid and reproduce very
well the full Color Glass Condensate results in its validity region.

##  I Introduction

Particle production at high energies in collisions of a light projectile with
heavy nuclei is typically calculated within the framework of the Color Glass
Condensate (CGC) [1, 2]. This weak-coupling, non-perturbative effective field
theory treats the projectile as a dilute object with a wave function
computable order-by-order in perturbation theory, often employing light-cone
perturbation theory [3]. The projectile partons subsequently rescatter off the
nuclear target, modeled as an ensemble of classical fields. These
rescatterings are resummed to all orders in the target’s high density, which,
in the eikonal approximation where terms subleading in energy are neglected in
the calculation of observables, translates to describing them as Wilson lines
at fixed transverse positions.

The CGC yields non-linear evolution equations for the target ensembles of
Wilson lines, leading to gluon density saturation for gluons with transverse
momenta smaller than the hadron or nucleus saturation scale, $Q_{s}$.
Unveiling evidence of such novel non-linear Quantum Chromodynamics (QCD)
phenomena and of saturation has been a central focus in Deep Inelastic
Scattering (DIS) at HERA at DESY, and in proton-proton, proton-nucleus,
nucleus-nucleus collisions, and photoproduction in ultraperipheral collisions,
at the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron
Collider (LHC) at CERN, see e.g. [4] and references therein. Despite claims at
HERA [5, 6] and RHIC [7, 8, 9, 10, 11], no conclusive evidence exists to date.

Future studies will continue at RHIC-II and LHC Runs from 3 onwards, with the
dedicated subdetector FoCal [12] proposed for installation in ALICE.
Furthermore, the Electron Ion Collider (EIC) [13, 14], to be built at BNL and
expected to begin data acquisition in the mid-2030s, offers exceptionally high
luminosities, variable center-of-mass energies ranging from $20-140$
GeV/nucleon, and the possibility of utilizing multiple nuclear targets and
both lepton projectile and proton/light ion target polarization. The quest for
the saturation regime constitutes one of the pillars of its physics program.

The eikonal approximation simplifies the treatment of target gluon fields in
three key ways. First, it treats the fields as localized in the longitudinal
direction, akin to a ”shockwave,” neglecting their longitudinal extent.
Second, it focuses solely on the leading longitudinal component of the field,
the one enhanced by the boost, effectively ignoring other components. Finally,
the dynamics of the target itself are disregarded.

As a result, the target gluon field becomes a single component with a delta-
function dependence on the longitudinal coordinate, leaving only the
transverse position as relevant variable. In light-cone coordinates
$(x^{+},x^{-})$, where $x^{\pm}=(x^{0}\pm x^{3})/\sqrt{2}$ for a left-moving
target, this simplified field can be expressed as

$A^{\mu}_{a}(x^{+},x^{-},\mathbf{x})\approx\delta^{\mu-}\delta(x^{+})\alpha_{a}(\mathbf{x}),$
(1)

where $a=1,\dots,N_{c}^{2}-1$ is the color index, boldface denotes transverse
coordinates and the function $\alpha_{a}$ represents the transverse dependence
of the gauge field.

Studies of non-eikonal corrections in the CGC emerged roughly a decade ago
[15, 16]. These studies relaxed the shockwave approximation by considering a
target with a finite longitudinal extent, a central focus in medium-induced
gluon radiation studies for the past three decades (see, e.g., the reviews
[17, 18, 19]). Subsequently, corrections to the other two eikonal
approximations were calculated for scalar, quark, and gluon propagators,
additionally incorporating quark exchanges with the target [20, 21, 22, 23,
24, 25, 26, 27, 28, 29]. These corrections are crucial for spin physics, as
spin exchanges are subleading in energy and therefore subeikonal [30, 31, 32,
33, 34, 35, 36, 37, 38].

Until now, numerical studies of the impact of non-eikonal corrections on
unpolarized observables have exclusively focused on the relaxation of the
shockwave approximation in single and double inclusive particle production in
proton-proton and proton-nucleus collisions [39, 40, 41, 42, 43]. For the
remaining non-eikonal corrections, either the cross sections have not yet been
calculated with sufficient accuracy, or they necessitate further modeling of
the target field averages. The general conclusion from these studies suggests
that the considered non-eikonal effects are significant for
$\sqrt{s_{NN}}\lesssim 100$ GeV, thus below the top RHIC energies. However,
these are precisely the energies of interest for the EIC, which provides the
motivation for examining the size of non-eikonal effects in DIS observables.

Recently, dijet production in DIS has garnered significant attention for its
potential to reveal saturation effects [44, 45, 46, 47, 48, 49]. This includes
the recent computation of the cross section at next-to-leading order in the
coupling constant [50, 51, 52, 53, 54]. Notably, the relevant cross section
incorporating non-eikonal corrections was computed in [55]. However, its
numerical implementation necessitates knowledge of target field averages
beyond the standard $\langle A^{-}A^{-}\rangle$ provided by the McLerran-
Venugopalan (MV) model [56, 57]. Specifically, it requires information about
averages involving transverse components ${\bf A}^{i}$, such as $\langle
A^{-}{\bf A}^{i}\rangle$ and $\langle{\bf A}^{i}{\bf A}^{j}\rangle$, and their
dependence on the longitudinal coordinate $x^{-}$. Unfortunately, such
information is currently unavailable, except for the spin contributions
computed in [36].

Furthermore, dijet production in DIS holds additional interest due to its
potential link between CGC calculations and transverse momentum distributions
(TMDs) in the so-called ”correlation limit” [58, 59, 60, 61, 44, 62, 63, 64,
65]. These TMDs encode three-dimensional information about the hadron and
nuclear structure. Consider two jets with transverse momenta
$\mathbf{k}_{1,2}$ and light-cone momentum fractions $z_{1,2}$ acquired from
the (virtual) photon. We define the momentum imbalance as ${\bf\Delta}={\bf
k}_{1}+{\bf k}_{2}$ and the relative momentum as ${\bf P}=z_{2}{\bf
k}_{1}-z_{1}{\bf k}_{2}$. The correlation limit corresponds to $|{\bf
P}|\gg|{\bf\Delta}|$, under which a connection between target averages of
Wilson lines and target TMDs at momentum fraction $x=0$ is established.

This work investigates the impact of non-eikonal effects on dijet production
in nuclear DIS. Building upon the framework established in [55], which
incorporates next-to-eikonal corrections arising from the transverse component
and $x^{-}$ dependence of the field, we restrict our analysis to the $x^{-}$
independent case. To incorporate the transverse component into the target
average of field correlators, we develop an extension of the MV model,
inspired by the work of [36] for polarized targets.

Subsequently, we investigate the dijet production differential cross section
in the dilute limit of the CGC111This approximation corresponds to high energy
factorization [66, 67] and is also known as the “glasma graph” approximation
in two-particle correlation studies [68]., where partons interact solely
through two-gluon exchange. This approximation holds for large values of
transverse momenta $|{\bf P}|$ and $|{\bf\Delta}|$, exceeding the saturation
scale $Q_{s}$. Leveraging this approximation, we achieve numerical solutions
for the cross section in both the eikonal (see Eqs. (76) and (77)) and next-
to-eikonal (see Eqs. (93) and (94)) scenarios.

Moreover, we conduct a numerical analysis of the analytical results, focusing
on the azimuthal anisotropy expansion of the cross section for $e^{-}{\rm Au}$
collisions at $\sqrt{s}=90$ GeV. Our findings reveal significant non-eikonal
corrections (up to 10%) for moderate values of transverse momenta $P_{\perp}$
and $\Delta_{\perp}$. Notably, non-eikonal effects induce odd azimuthal
harmonics, which are absent in the eikonal approximation.

To finalize our analysis, we compare the dilute limit with the dense and
correlation limit of the CGC. We observe that the eikonal dilute limit of the
DIS differential cross section accurately reproduces both the full eikonal CGC
result and its correlation limit within the kinematic region $|{\bf
P}|,|{\bf\Delta}|\gtrsim Q_{s}$. The validity of the different approximations
is illustrated in Fig. 1, where $P_{\perp}=|{\bf
P}|,\Delta_{\perp}=|{\bf\Delta}|$.

Figure 1: Comparison of our approach with the correlation limit of the CGC. We
see that when $\Delta_{\perp}\ll P_{\perp}$ we get a good agreement between
both approaches as expected.

The manuscript is organized as follows. In Section II we provide the setup for
the calculation. In Section III.3 we derive the field correlators between the
target fields taking into account the $A^{+}$ and ${\bf A}^{i}$ components of
the background field, absent in the standard MV model. In Section IV we review
the differential cross section for $q{\bar{q}}$ production in DIS beyond the
eikonal approximation. In Section V we solve the differential cross section in
the dilute limit of the CGC. In Section VI we present our numerical results.
Section VII contains our conclusions. Finally, Appendix A contains a brief
review of the dense limit of the Wilson line averages as well as its
correlation limit.

##  II Notations and conventions on the nuclear wave function

In this section, we employ the Light-Cone Wave Function (LCWF) formalism,
detailed in references [3, 69, 36] and their citations, to compute the
expectation value of an observable within an ultra-relativistic, large
nucleus. Using this formalism, we describe an unpolarized, left-moving nucleus
consisting of $A$ nucleons with 3-momenta denoted as
$\vec{p}_{i}\equiv(p_{i}^{-},\mathbf{p}_{i})$ via the following wave function:

$\ket{A}=\int\prod_{i=1}^{A}\frac{d^{3}\vec{p}_{i}}{2(2\pi)^{3}p_{i}^{-}}\psi(\\{\vec{p}\\})\ket{\\{\vec{p}\\}}.$
(2)

Here, $\\{\vec{p}\\}\equiv\vec{p}_{1},\dots,\vec{p}_{A}$ represents the
3-momenta of all nucleons, and $\ket{\\{\vec{p}\\}}$ is the Fock space basis
vector for a system with $A$ nucleons. The function
$\psi(\\{\vec{p}\\})=\braket{\\{\vec{p}\\}}{A}$ denotes the amplitude
associated with a specific momentum configuration. We normalize the nuclear
wave function as $\braket{A}{A}=1$, ensuring that the light-cone amplitude
also satisfies the normalization condition:

$\int\prod_{i=1}^{A}\frac{d^{3}\vec{p}_{i}}{2(2\pi)^{3}p_{i}^{-}}|\psi(\\{\vec{p}\\})|^{2}=1.$
(3)

Furthermore, the single-particle state of a nucleon with momentum $\vec{p}$
fulfills

$\braket{\vec{p}^{\prime}}{\vec{p}}=2p^{-}(2\pi)^{3}\delta^{(3)}(\vec{p}-\vec{p}^{\prime}).$
(4)

Referring to Eq. 2, and considering an observable characterized by the
operator $\hat{\mathcal{O}}$, its expectation value within the nuclear
ensemble can be expressed as follows:

$\displaystyle\langle\hat{O}\rangle_{A}\equiv\bra{A}\hat{\mathcal{O}}\ket{A}$
$\displaystyle=\int\prod_{i=1}^{A}\frac{d^{3}\vec{P}_{i}}{2(2\pi)^{3}\left(P_{i}^{-}+\frac{q_{i}^{-}}{2}\right)}\frac{d^{3}\vec{q}_{i}}{2(2\pi)^{3}\left(P_{i}^{-}-\frac{q_{i}^{-}}{2}\right)}\psi^{*}\left(\Big{\\{}\vec{P}-\frac{\vec{q}}{2}\Big{\\}}\right)\psi\left(\Big{\\{}\vec{P}+\frac{\vec{q}}{2}\Big{\\}}\right)$
$\displaystyle\times\Big{\langle}\Big{\\{}\vec{P}-\frac{\vec{q}}{2}\Big{\\}}\Big{|}\hat{\mathcal{O}}\Big{|}\Big{\\{}\vec{P}+\frac{\vec{q}}{2}\Big{\\}}\Big{\rangle}.$
(5)

For convenience, we can express Section II in terms of the Wigner distribution
of the quantum system by rephrasing it as follows:

$\displaystyle\langle\hat{O}\rangle_{A}=\prod_{i=1}^{A}\int\frac{d^{3}\vec{P}_{i}}{(2\pi)^{3}}\int
d^{3}\vec{B}_{i}W(\\{\vec{P},\vec{B}\\})O(\\{\vec{P},\vec{B}\\}),$ (6)

where

$W(\\{\vec{P},\vec{B}\\})=\int\prod_{i=1}^{A}\frac{d^{3}\vec{q}_{i}}{2(2\pi)^{3}}\frac{P_{i}^{-}}{\left(P_{i}^{-}+\frac{q_{i}^{-}}{2}\right)\left(P_{i}^{-}-\frac{q_{i}^{-}}{2}\right)}e^{-i\vec{q}_{i}\cdot\vec{B}_{i}}\psi^{*}\left(\Big{\\{}\vec{P}-\frac{\vec{q}}{2}\Big{\\}}\right)\psi\left(\Big{\\{}\vec{P}+\frac{\vec{q}}{2}\Big{\\}}\right)$
(7)

represents the multi-particle Wigner distribution of the nuclear ensemble,
which provides information on the spatial and momentum distribution of the
nucleons within the nucleus. Furthermore, we have defined

$O(\\{\vec{P},\vec{B}\\})=\int\prod_{i=1}^{A}\frac{d^{3}\vec{q}_{i}}{(2\pi)^{3}2P_{i}^{-}}e^{i\vec{q}_{i}\cdot\vec{B}_{i}}\Big{\langle}\Big{\\{}\vec{P}-\frac{\vec{q}}{2}\Big{\\}}\Big{|}\hat{\mathcal{O}}\Big{|}\Big{\\{}\vec{P}+\frac{\vec{q}}{2}\Big{\\}}\Big{\rangle},$
(8)

which accounts for the expectation value of the observable at a given
configuration of the nucleons inside the nucleus.

While Eq. 6 is entirely general and model-independent, in practical
applications the Wigner distribution is notoriously difficult to calculate
directly and often requires specification through a particular theoretical
framework. Given the vast number of nucleons comprising the nucleus,
correlations between them can often be disregarded as they scale as
$1/\sqrt{A}$. This allows us to factorize the amplitude for the multi-particle
system as a product of single-particle amplitudes, leading to a more tractable
expression for the Wigner distribution:

$W(\\{\vec{P},\vec{B}\\})=\prod_{i=1}^{A}W(\vec{P}_{i},\vec{B}_{i}).$ (9)

Here, it is important to note that we use the same notation, $W$, to represent
both the multi-particle and single-particle Wigner distributions.

Ultimately, the single-particle Wigner function can be described using various
models. Here, we adopt a Glauber-like model where the momentum and spatial
distributions are factorized:

$W(\vec{P},\vec{B})=\zeta(\vec{P})\,\frac{\rho_{A}(\vec{B})}{A},$ (10)

where $\rho_{A}(\vec{x})$ represents the normalized nucleon density, such that
$\int d^{3}\vec{x}\rho_{A}(\vec{x})=A$. Additionally,
$\zeta(\vec{P})\equiv|\phi(\vec{P})|^{2}/2P^{-}$ accounts for the momentum
distribution of the nucleons within the nucleus.

Following the conventions of the McLerran-Venugopalan model as detailed in
[36], we assume that each nucleon possesses a common, large longitudinal
momentum denoted as $P^{-}_{\text{N}}=P^{-}_{A}/A$, where $\vec{P}_{A}$
corresponds to the 3-momentum of the nucleus. Consequently, for a nucleus with
negligible transverse momentum, we can write222In principle, for a nucleus
with negligible transverse momentum, the individual transverse momenta of the
nucleons do not necessarily have to be zero. However, in the context of this
analysis, the expectation value will solely rely on the average transverse
momentum of the nucleons, which will ultimately approach zero. To streamline
our analysis, we make the simplifying assumption that the nucleons carry no
transverse momentum.

$\zeta(\vec{P})=(2\pi)^{3}\delta(P^{-}-P^{-}_{\rm N})\delta^{(2)}({\bf P}).$
(11)

This essentially implies that all nucleons share the same large longitudinal
momentum ($P^{-}_{\rm N}$), while their transverse momenta are distributed
around zero.

By substituting Eqs. 9, 10 and 11 into Eq. 6, we arrive at a simplified
expression for the expectation value:

$\displaystyle\langle\hat{O}\rangle_{A}=\prod_{i=1}^{A}\int
d^{3}\vec{B}_{i}\frac{\rho_{A}(\vec{B}_{i})}{A}O\left(\\{\vec{P}=(P^{-}_{\rm
N},0),\vec{B}\\}\right).$ (12)

This result highlights that the expectation value only depends on the integral
over the nuclear density $\rho_{A}(\vec{B}_{i})$ and the observable evaluated
at the common longitudinal momentum $P^{-}_{\rm N}$ and zero transverse
momentum.

In the CGC framework, the nuclear wave function is characterized by sources
represented as large-$x$ recoilless valence partons, while the small-$x$
gluons are determined by the solutions of the classical Yang-Mills equations
associated with these sources. Consequently, the expectation values of
observables are computed exclusively through the degrees of freedom of the
valence partons. Under the assumption that each nucleon is comprised of
$N_{c}$ quarks333While a complete non-eikonal treatment would necessitate
incorporating “valence” gluons and antiquarks into the nuclear wave function,
as demonstrated in [36], it is crucial to emphasize that in the present study
eikonal corrections solely focus on the classical gluon field generated by the
valence quarks. such that the nucleons are singlets, the single nucleon state
is written as follows [69]:

$\ket{\vec{p}}=\frac{1}{\sqrt{N_{c}!}}\int\prod_{i=1}^{N_{c}}\frac{d^{3}\vec{k}_{i}}{2(2\pi)^{3}k_{i}^{-}}\Psi(\\{\vec{k}\\})\epsilon_{\alpha_{1},\dots,\alpha_{N_{c}}}\ket{\\{\vec{k},\alpha\\}}2(2\pi)^{3}p^{-}\delta^{(3)}\left(\vec{p}-\sum_{n=1}^{N_{c}}\vec{k}_{n}\right),$
(13)

where the basis for the Fock Space of a system with $N_{c}$ quarks is denoted
as
$\ket{\\{\vec{k},\alpha\\}}=\ket{\vec{k}_{1},\alpha_{1}}\otimes\cdots\otimes\ket{\vec{k}_{N_{c}},\alpha_{N_{c}}}$.
Here, $\ket{\vec{k},\alpha}$ represents a single quark state characterized by
momentum $\vec{k}$ and color $\alpha$, while $\Psi(\\{\vec{k}\\})$ represents
the amplitude associated with this particular configuration. The Levi-Civita
tensor, denoted as $\epsilon_{\alpha_{1},\dots,\alpha_{N_{c}}}$, ensures that
the state is a singlet. It is important to mention that, since this study
assumes an unpolarized nucleus, we omit the helicity quantum numbers. The
single quark states are defined by
$\ket{\vec{k},\alpha}=\hat{b}^{\dagger}_{\alpha}(\vec{k})\ket{0}$, with the
creation and annihilation operators adhering to the following normalization:

$\\{\hat{b}_{\alpha^{\prime}}(\vec{k}^{\prime}),\hat{b}^{\dagger}_{\alpha}(\vec{k})\\}=\delta_{\alpha\alpha^{\prime}}2k^{-}(2\pi)^{3}\delta^{(3)}(\vec{k}-\vec{k}^{\prime}),\qquad\\{\hat{b}^{\dagger}_{\alpha^{\prime}}(\vec{k}^{\prime}),\hat{b}^{\dagger}_{\alpha}(\vec{k})\\}=\\{\hat{b}_{\alpha^{\prime}}(\vec{k}^{\prime}),\hat{b}_{\alpha}(\vec{k})\\}=0.$
(14)

The state of $A$ color-neutral and independent nucleons composed by $N_{c}$
quarks is represented as

$\ket{\\{\vec{p}\\}}=\prod_{i=1}^{A}\left[\frac{1}{\sqrt{N_{c}!}}\int\prod_{j=1}^{N_{c}}\frac{d^{3}\vec{k}_{ij}}{2(2\pi)^{3}k_{ij}^{-}}\epsilon_{\\{\alpha_{i}\\}}\Psi(\\{\vec{k}_{i}\\})\ket{\\{\vec{k}_{i},\alpha_{i}\\}}16\pi^{3}p^{-}\delta^{(3)}\left(\vec{p}-\sum_{n=1}^{N_{c}}\vec{k}_{in}\right)\right],$
(15)

where we use the notation
$\\{\alpha_{i}\\}\equiv\alpha_{i1},\dots,\alpha_{iN_{c}}$.

We begin by assuming that quarks within the nucleon wave function behave
independently, separated by distances roughly equivalent to the nucleon’s size
($r_{\rm N}\sim\Lambda_{\rm QCD}^{-1}$, where $\Lambda_{\rm QCD}$ denotes the
QCD non-perturbative scale). This allows us to make similar assumptions as
done at the nucleon level and introduce the quark distribution,
$\rho_{N_{c}}(\vec{x})$, normalized such that $\int
d^{3}\vec{x}\rho_{N_{c}}(\vec{x})=N_{c}$. All in all, we arrive at the
following expression:

$\bra{\\{\vec{p}^{\prime}\\}}\hat{\mathcal{O}}\ket{\\{\vec{p}\\}}=\frac{1}{A!}\Big{(}\braket{\vec{p}_{1}^{\prime}}{\vec{p}_{1}}\cdots\braket{\vec{p}_{A}^{\prime}}{\vec{p}_{A}}+{\rm
permutations}\Big{)}\prod_{i=1}^{A}\prod_{j=1}^{N}\int
d^{3}\vec{b}_{ij}\frac{\rho_{N}(\vec{b}_{ij}-\vec{B}_{i})}{N}\langle\hat{\mathcal{O}}\rangle_{q}(\\{\vec{b}\\}),$
(16)

where $\langle\hat{\mathcal{O}}\rangle_{q}(\\{\vec{b}\\})$ is the expectation
value of the operator at a given quark configuration, defined as follows:

$\langle\hat{\mathcal{O}}\rangle_{q}(\\{\vec{b}\\})=\prod_{i=1}^{A}\left[\frac{1}{N_{c}!}\left(\prod_{j=1}^{N_{c}}\int\frac{d^{3}\vec{q}_{ij}}{2(2\pi)^{3}P_{q}^{-}}e^{i\vec{q}_{ij}\cdot\vec{b}_{ij}}\right)\epsilon_{\\{\alpha_{i}\\}}\epsilon_{\\{\bar{\alpha}_{i}\\}}\right]\Big{\langle}\Big{\\{}P_{q}^{-}-\frac{q^{-}}{2},-\frac{\bf
q}{2},\bar{\alpha}\Big{\\}}\Big{|}\hat{\mathcal{O}}\Big{|}\Big{\\{}P_{q}^{-}+\frac{q^{-}}{2},\frac{\bf
q}{2},\alpha\Big{\\}}\Big{\rangle}.$ (17)

Here, $P_{q}^{-}\equiv P_{A}^{-}/AN_{c}$ defines the longitudinal momentum of
quarks within the nucleus. Since $P_{A}^{-}=M_{A}e^{\omega}$ (where $M_{A}$
and $\omega$ represent the nucleus’s mass and rapidity, respectively),
$P_{q}^{-}$ becomes the dominant momentum scale when the nucleus is highly
boosted.

Finally, by substituting Eq. 16 into Eq. 8, we can express the expectation
value of the observable $\hat{\mathcal{O}}$ in the nuclear ensemble, as
presented in Eq. 12, as follows:

$\langle\hat{O}\rangle_{A}=\prod_{i=1}^{A}\prod_{j=1}^{N_{c}}\int
d^{3}\vec{B}_{i}\frac{\rho_{A}(\vec{B}_{i})}{A}\int
d^{3}\vec{b}_{ij}\frac{\rho_{N_{c}}(\vec{b}_{ij}-\vec{B}_{i})}{N_{c}}\langle\hat{\mathcal{O}}\rangle_{q}(\\{\vec{b}\\}).$
(18)

We can further approximate by recognizing that $\rho_{N_{c}}(\vec{b}-\vec{B})$
is predominantly centered around $\vec{B}=\vec{b}$ and confined within the
nucleon’s size. Meanwhile, $\rho_{A}(\vec{B})$ extends across the entire
nuclear volume, allowing us to treat it as nearly constant over the scale of a
single nucleon. With this, we can write

$\langle\hat{O}\rangle_{A}=\prod_{i=1}^{A}\prod_{j=1}^{N_{c}}\int
d^{3}\vec{b}_{ij}\frac{\rho_{A}(\vec{b}_{ij})}{A}\langle\hat{\mathcal{O}}\rangle_{q}(\\{\vec{b}\\}).$
(19)

##  III Color source and field correlator

We turn our attention to computing the expectation value of the color current
density correlator within the nuclear state discussed previously. As a first
step, we derive the expectation value of the current operator for illustrative
purposes. We then proceed to calculate the correlator itself.

###  III.1 The color current density operator

The color current density operator, at light-cone time $x^{-}=0$, is given in
terms of the fermion fields by the following expression:

$\hat{J}_{a}^{\mu}(x^{-}=0,\vec{x})=g\bar{\psi}_{\alpha}(\vec{x})\gamma^{\mu}t^{a}_{\alpha\beta}\psi_{\beta}(\vec{x}).$
(20)

In this equation, $g$ represents the QCD coupling constant, $t^{a}$ denotes
the SU($N_{c}$) generators in the fundamental representation, $\gamma^{\mu}$
corresponds to the gamma matrices, and the fermion fields with color $\alpha$
are given by

$\displaystyle\psi_{\alpha}(\vec{x})=\int_{\vec{p}}\hat{b}_{\alpha,\sigma}(\vec{p})\frac{u_{\sigma}(\vec{p})}{\sqrt{2p^{-}}}e^{-i\vec{p}\cdot\vec{x}},\qquad\bar{\psi}_{\alpha}(\vec{x})=\int_{\vec{p}}\hat{b}^{\dagger}_{\alpha,\sigma}(\vec{p})\frac{\bar{u}_{\sigma}(\vec{p})}{\sqrt{2p^{-}}}e^{i\vec{p}\cdot\vec{x}}.$
(21)

Here, we have introduced the notation

$\int_{\vec{p}}\equiv\int\frac{d^{3}\vec{p}}{(2\pi)^{3}\sqrt{2p^{-}}}.$ (22)

We exclude the antiquark component of the field due to the exclusively quark
composition of the nucleus, as the antiquark component holds no relevance to
our present analysis. While we introduce the helicity index $\sigma$
temporarily, it is to be noted that the nucleus is unpolarized, necessitating
an averaging over helicity at the end of the calculation.

Therefore, by substituting Eq. 21 into Eq. 20, we can express the color
current density as follows:

$\hat{J}_{a}^{\mu}(\vec{x})=gt_{\alpha\beta}^{a}\int_{\vec{p},\vec{p}^{\prime}}e^{-i\vec{x}\cdot(\vec{p}-\vec{{p}}^{\prime})}\Gamma^{\mu}_{\sigma{\sigma}^{\prime}}(\vec{p},\vec{{p}}^{\prime})\
\hat{b}_{\alpha,{\sigma}^{\prime}}^{\dagger}(\vec{{p}}^{\prime})\hat{b}_{\beta,\sigma}(\vec{p}),$
(23)

where

$\displaystyle\Gamma^{\mu}_{\sigma,{\sigma}^{\prime}}(\vec{p},\vec{{p}}^{\prime})=\frac{\bar{u}_{{\sigma}^{\prime}}(\vec{p}^{\prime})}{\sqrt{2p^{\prime-}}}\gamma^{\mu}\frac{u_{\sigma}(\vec{p})}{\sqrt{2p^{-}}}.$
(24)

To compute this object, we can utilize the Brodsky-Lepage spinors for left-
moving massless fermions [70, 36]:

$\displaystyle
u_{\sigma}(\vec{p})=\frac{1}{\sqrt{\sqrt{2}p^{-}}}\left(\sqrt{2}p^{-}+\gamma^{0}\gamma^{i}{\bf
p}^{i}\right)\chi_{\sigma},$ (25)

where $\chi_{+}=1/\sqrt{2}(1,0,-1,0)^{T}$ and
$\chi_{-}=1/\sqrt{2}(0,1,0,1)^{T}$, such that

$\displaystyle\Gamma^{-}_{\sigma{\sigma}^{\prime}}(\vec{p},\vec{{p}}^{\prime})$
$\displaystyle=\delta_{\sigma\sigma^{\prime}},$ (26a)
$\displaystyle\Gamma^{+}_{\sigma{\sigma}^{\prime}}(\vec{p},\vec{{p}}^{\prime})$
$\displaystyle=\delta_{\sigma\sigma^{\prime}}\frac{{\bf
p}^{i}}{p^{-}}\frac{{\bf p}^{\prime i}-i\sigma\epsilon^{ij}{\bf p}^{\prime
j}}{2p^{\prime-}},$ (26b)
$\displaystyle\Gamma^{i}_{\sigma{\sigma}^{\prime}}(\vec{p},\vec{{p}}^{\prime})$
$\displaystyle=\delta_{\sigma\sigma^{\prime}}\left(\frac{{\bf p}^{\prime
i}-i\sigma\epsilon^{ij}{\bf p}^{\prime j}}{2p^{\prime-}}+\frac{{\bf
p}^{i}+i\sigma\epsilon^{ij}{\bf p}^{j}}{2p^{-}}\right).$ (26c)

We observe that the transverse component of the current, represented by the
helicity tensor defined in Eq. 26c, was previously investigated in [36] for
the specific case of a polarized target within the study of helicity-dependent
processes.

To calculate the expectation value of the density current presented in Eq. 23,
we leverage the following identities derived from the commutation relations in
Eq. 14. Within an $N$-particle Fock basis characterized by momentum and color
quantum numbers denoted as $s\equiv\vec{k},\alpha$, the following identities
apply:

$\displaystyle\hat{b}_{\beta}(\vec{p})\ket{s_{1},\dots,s_{N}}$
$\displaystyle=\sum_{i=1}^{N}(-1)^{i+1}\delta_{\beta\alpha_{i}}16\pi^{3}p^{-}\delta^{(3)}(\vec{p}-\vec{k}_{i})\ket{s_{1},\dots,s_{i-1},s_{i+1},\dots,s_{N}},$
(27a)
$\displaystyle\hat{b}_{\alpha}^{\dagger}(\vec{{p}}^{\,\prime})\hat{b}_{\beta}(\vec{p})\ket{s_{1},\dots,s_{N}}$
$\displaystyle=\sum_{i=1}^{N}\delta_{\beta\alpha_{i}}16\pi^{3}p^{-}\delta^{(3)}(\vec{p}-\vec{k}_{i})\ket{s_{1},\dots,s_{N}\
/\ s_{i}=\vec{p}^{\,\prime},\alpha}.$ (27b)

Consequently, the matrix element of the operator
$\hat{b}_{\alpha}^{\dagger}(\vec{{p}}^{\,\prime})\hat{b}_{\beta}(\vec{p})$ can
be expressed as

$\displaystyle\bra{s_{1}^{\prime},\cdots,s_{N}^{\prime}}\hat{b}_{\alpha}^{\dagger}(\vec{{p}}^{\,\prime})\hat{b}_{\beta}(\vec{p})\ket{s_{1},\dots,s_{N}}$
$\displaystyle=\sum_{i=1}^{N}\delta_{\beta\alpha}16\pi^{3}p^{-}\delta^{(3)}(\vec{p}-\vec{k}_{i})16\pi^{3}p^{\prime-}\delta^{(3)}(\vec{p}^{\,\prime}-\vec{k}^{\prime}_{i})$
$\displaystyle\times\prod_{j\neq
i}^{N}\delta_{\alpha_{j}\alpha^{\prime}_{j}}16\pi^{3}k_{j}^{-}\delta^{(3)}(\vec{k}_{j}-\vec{k}^{\prime}_{j})+\cdots,$
(28)

where we have retained only the diagonal term, which is the component that
dominates when the particles are independent.

Thus, by using Sections III.1, 23 and 17 and performing some algebraic
manipulations, we arrive at the following expression:

$\displaystyle\langle\hat{J}_{a}^{\mu}(\vec{x})\rangle_{q}(\\{\vec{b}\\})$
$\displaystyle=\frac{g{\rm
Tr}[t^{a}]}{N_{c}}\sum_{i=1}^{A}\sum_{j=1}^{N_{c}}\int\frac{d^{3}\vec{q}_{ij}}{(2\pi)^{3}}e^{-i\vec{q}_{ij}\cdot(\vec{x}-\vec{b}_{ij})}\frac{1}{2}\Gamma_{\sigma\sigma}^{\mu}\left({P}_{q}^{-}+\frac{q_{ij}^{-}}{2},\frac{{\bf
q}_{ij}}{2};{P}_{q}^{-}-\frac{q_{ij}^{-}}{2},-\frac{{\bf
q}_{ij}}{2}\right)\sqrt{1-\left(\frac{q_{ij}^{-}}{2P_{q}^{-}}\right)^{2}}$
$\displaystyle\approx\frac{g{\rm
Tr}[t^{a}]}{N_{c}}\sum_{i=1}^{A}\sum_{j=1}^{N_{c}}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}e^{-i\vec{q}\cdot(\vec{x}-\vec{b}_{ij})}\frac{1}{2}\Gamma_{\sigma\sigma}^{\mu}\left({P}_{q}^{-},\frac{{\bf
q}}{2};{P}_{q}^{-},-\frac{{\bf q}}{2}\right).$ (29)

Here, we have taken an average over helicity by performing the trace over
$\sigma$ and dividing by 2. Furthermore, we have neglected terms proportional
to $q_{ij}^{-}$ because they scale as $q_{ij}^{-}\sim 1/A$ after Fourier
transformation, therefore they are subdominant in the large-$A$ limit.
Ultimately, the expectation value of the color current density is
proportionate to ${\rm Tr}[t^{a}]=0$ and, hence, it vanishes, in line with the
expectation for color-neutral nucleons.

###  III.2 Color current density correlator

We proceed to calculate the color current density correlator within the
nucleus. The operator for the helicity-averaged density correlator is
expressed as

$\displaystyle\hat{J}_{a}^{\mu}(\vec{x})\hat{J}_{b}^{\nu}(\vec{y})=g^{2}t_{\alpha_{1}\beta_{1}}^{a}t_{\alpha_{2}\beta_{2}}^{b}\int_{\vec{p}_{1},\vec{p}_{1}^{\,\prime},\vec{p}_{2},\vec{p}_{2}^{\,\prime}}$
$\displaystyle
e^{-i\vec{x}\cdot(\vec{p}_{1}-\vec{p}_{1}^{\,\prime})-i\vec{y}\cdot(\vec{p}_{2}-\vec{p}_{2}^{\,\prime})}\frac{1}{2}\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}(\vec{p}_{1},\vec{p}_{1}^{\,\prime})\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}(\vec{p}_{2},\vec{p}_{2}^{\,\prime})$
$\displaystyle\times\hat{b}_{\alpha_{1}}^{\dagger}(\vec{p}_{1}^{\,\prime})\hat{b}_{\beta_{1}}(\vec{p}_{1})\hat{b}_{\alpha_{2}}^{\dagger}(\vec{p}_{2}^{\,\prime})\hat{b}_{\beta_{2}}(\vec{p}_{2}).$
(30)

For practical applications, we focus on the regime where
$|\vec{x}-\vec{y}|\ll\Lambda_{\rm QCD}^{-1}$, which permits the application of
weak coupling techniques. Based on the assumption that valence quarks have a
typical separation of the order of $\Lambda_{\rm QCD}^{-1}$, both currents
probe the same quark within the perturbative regime. From an operator
perspective, this implies that the one-body component in Section III.2
dominates over the two-body component. Consequently, we can make the following
approximation:

$\displaystyle\hat{J}_{a}^{\mu}(\vec{x})\hat{J}_{b}^{\nu}(\vec{y})$
$\displaystyle\approx
g^{2}t_{\alpha_{1}\beta_{1}}^{a}t_{\alpha_{2}\beta_{2}}^{b}\int_{\vec{p}_{1},\vec{p}_{1}^{\,\prime},\vec{p}_{2},\vec{p}_{2}^{\,\prime}}e^{-i\vec{x}\cdot(\vec{p}_{1}-\vec{p}_{1}^{\,\prime})-i\vec{y}\cdot(\vec{p}_{2}-\vec{p}_{2}^{\,\prime})}\frac{1}{2}\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}(\vec{p}_{1},\vec{p}_{1}^{\,\prime})\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}(\vec{p}_{2},\vec{p}_{2}^{\,\prime})$
$\displaystyle\times\delta_{\alpha_{2}\beta_{1}}16\pi^{3}p_{1}^{-}\delta^{(3)}(\vec{p}_{1}-\vec{p}_{2}^{\,\prime})\hat{b}_{\alpha_{1}}^{\dagger}(\vec{p}_{1}^{\,\prime})\hat{b}_{\beta_{2}}(\vec{p}_{2})$
$\displaystyle=g^{2}t_{\alpha_{1}\beta_{1}}^{a}t_{\beta_{1}\beta_{2}}^{b}\int_{\vec{p}_{1},\vec{p}_{1}^{\,\prime},\vec{p}_{2}}e^{-i\vec{x}\cdot(\vec{p}_{1}-\vec{p}_{1}^{\,\prime})-i\vec{y}\cdot(\vec{p}_{2}-\vec{p}_{1})}\sqrt{2p_{1}^{-}}\frac{1}{2}\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}(\vec{p}_{1},\vec{p}_{1}^{\,\prime})\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}(\vec{p}_{2},\vec{p}_{1})\hat{b}_{\alpha_{1}}^{\dagger}(\vec{p}_{1}^{\,\prime})\hat{b}_{\beta_{2}}(\vec{p}_{2}).$
(31)

By once more employing Sections III.1 and 17, in conjunction with the
aforementioned approximations, we can derive the following expression:

$\displaystyle\Bigg{\langle}\hat{J}_{a}^{\mu}\left(\vec{B}+\frac{\vec{R}}{2}\right)\hat{J}_{b}^{\nu}\left(\vec{B}-\frac{\vec{R}}{2}\right)\Bigg{\rangle}_{q}(\\{\vec{b}\\})$
$\displaystyle=\frac{g^{2}{\rm
Tr}[t^{a}t^{b}]}{N_{c}}\sum_{i=1}^{A}\sum_{j=1}^{N_{c}}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{d^{3}\vec{p}}{(2\pi)^{3}}e^{-i\vec{q}\cdot(\vec{B}-\vec{b}_{ij})-i\vec{p}\cdot\vec{R}}$
$\displaystyle\hskip
0.0pt\times\frac{1}{2}\,\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}\left(P_{q}^{-}+p^{-},{\bf
p};P_{q}^{-},-\frac{\textbf{q}}{2}\right)\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}\left(P_{q}^{-},\frac{\textbf{q}}{2};P_{q}^{-}+p^{-},{\bf
p}\right).$ (32)

The dependence of the integrand on $p^{-}$ results in a non-trivial
correlation function, $F(R^{+};P_{q}^{-})$, after the transformation. This
function introduces a correlation length of $\lambda^{+}\sim 1/P_{q}^{-}$.
Notably, the longitudinal correlations become localized in the eikonal limit,
meaning $F(R^{+};P_{q}^{-}\to\infty)\sim\delta(R^{+})$. For practical
applications, where we integrate $R^{+}$ over the entire nuclear width, the
contributions from these longitudinal correlations scale as
$\lambda^{+}/L^{+}\sim A^{-1/3}$, where $L^{+}$ represents the light-cone
longitudinal width of the medium. Consequently, these contributions become
subdominant for large nuclei. For this reason, we do not take into account
corrections due to longitudinal particle correlations444While neglecting the
longitudinal correlations based on their subdominant scaling, we acknowledge
that their non-eikonal effects could be explicitly incorporated in a more
practical way by modifying the $\delta(R^{+})$ function. As explored in
references [39, 40], this can be achieved by substituting $\delta(R^{+})$ with
$\Theta(\lambda^{+}-|R^{+}|)/2\lambda^{+}$, where $\Theta$ represents the
Heaviside function.  in the present manuscript and we drop dependence of the
integrand on $p^{-}$:

$\displaystyle\Bigg{\langle}\hat{J}_{a}^{\mu}\left(\vec{B}+\frac{\vec{R}}{2}\right)\hat{J}_{b}^{\nu}\left(\vec{B}-\frac{\vec{R}}{2}\right)\Bigg{\rangle}_{q}(\\{\vec{b}\\})$
$\displaystyle=\frac{g^{2}\delta^{ab}}{N_{c}}\sum_{i=1}^{A}\sum_{j=1}^{N_{c}}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{d^{3}\vec{p}}{(2\pi)^{3}}e^{-i\vec{q}\cdot(\vec{B}-\vec{b}_{ij})-i\vec{p}\cdot\vec{R}}$
$\displaystyle\hskip
0.0pt\times\frac{1}{2}\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}\left(P_{q}^{-},{\bf
p};P_{q}^{-},-\frac{\textbf{q}}{2}\right)\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}\left(P_{q}^{-},\frac{\textbf{q}}{2};P_{q}^{-},{\bf
p}\right),$ (33)

where we have written ${\rm Tr}[t^{a}t^{b}]=\delta^{ab}/2$.

Finally, by substituting Section III.2 into Eq. 12, we arrive at the
expression for the color current density correlator in a large nucleus555We
use the following definitions for the transverse integrals:
$\displaystyle\int_{\bf x}\equiv\int d^{2}{\bf x},\qquad\int_{\bf
p}\equiv\frac{d^{2}{\bf p}}{(2\pi)^{2}}.$ (34) :

$\displaystyle\Bigg{\langle}\hat{J}_{a}^{\mu}\left(\vec{b}+\frac{\vec{r}}{2}\right)\hat{J}_{b}^{\nu}\left(\vec{b}-\frac{\vec{r}}{2}\right)\Bigg{\rangle}_{A}$
$\displaystyle=\delta^{ab}\delta(r^{+})\int_{\textbf{B}}\mu^{2}(b^{+},\textbf{b}-\textbf{B})\int_{{\bf
q},{\bf p}}e^{i\textbf{q}\cdot\textbf{B}+i\textbf{p}\cdot\textbf{r}}$
$\displaystyle\times\frac{1}{2}\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}\left(P_{q}^{-},{\bf
p};P_{q}^{-},-\frac{\textbf{q}}{2}\right)\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}\left(P_{q}^{-},\frac{\textbf{q}}{2};P_{q}^{-},{\bf
p}\right),$ (35)

where we have defined

$\displaystyle\mu^{2}(\vec{x})=\frac{g^{2}}{2}\rho_{A}(\vec{x}),$ (36)

which can be interpreted as the mean color charge density (squared) per unit
volume.

###  III.3 Color field correlation

Within the CGC framework, observables depend on the correlators of the color
field generated by nuclear current sources. When considering a highly boosted,
large nucleus, the resultant color field is strong and effectively classical.
Consequently, it can be expressed in terms of the color sources via solutions
to the classical Yang-Mills (YM) equations. However, solving these equations
for a strong, non-eikonal field presents significant analytical and numerical
challenges due to their inherent non-linearity.

To address this complexity, we adopt the dilute limit of the CGC. This allows
us to treat the field as weak yet classical, enabling us to get rid of the
subleading non-linear terms in the YM equations and enabling analytical
solutions even for non-eikonal configurations. Furthermore, we restrict
ourselves to light-cone time-independent fields. Under these assumptions, the
classical YM equation in covariant gauge ($\partial_{\mu}A^{\mu}=0$) for a
time-independent field takes the following form:

$\displaystyle\partial^{2}\hat{A}_{a}^{\mu}(\vec{x})=-\partial_{\perp}^{2}\hat{A}_{a}^{\mu}(\vec{x})=\hat{J}_{a}^{\mu}(\vec{x}),$
(37)

with the solution

$\displaystyle\hat{A}_{a}^{\mu}(\vec{x})=\int_{\bf z}\int_{\bf q}e^{i{\bf
q}\cdot({\bf x}-{\bf z})}\frac{\hat{J}_{a}^{\mu}(x^{+},{\bf z})}{{\bf
q}^{2}}.$ (38)

Consequently, the correlation of the color field within the nuclear ensemble
can be written as

$\displaystyle\big{\langle}\hat{A}_{a}^{\mu}(\vec{x})\hat{A}_{b}^{\nu}(\vec{y})\big{\rangle}_{A}=\int_{{\bf
z}_{1},{\bf z}_{2}}\int_{{\bf q}_{1},{\bf q}_{2}}e^{i{\bf q}_{1}\cdot({\bf
x}-{\bf z}_{1})+i{\bf q}_{2}\cdot({\bf y}-{\bf
z}_{2})}\frac{\big{\langle}\hat{J}_{a}^{\mu}(x^{+},{\bf
z}_{1})\hat{J}_{b}^{\nu}(y^{+},{\bf z}_{2})\big{\rangle}_{A}}{{\bf
q}_{1}^{2}{\bf q}_{2}^{2}}.$ (39)

Using the result of Section III.2, we arrive at

$\displaystyle\Bigg{\langle}A_{a}^{\mu}\left(\vec{b}+\frac{\vec{r}}{2}\right)A_{b}^{\nu}\left(\vec{b}-\frac{\vec{r}}{2}\right)\Bigg{\rangle}_{A}=\delta^{ab}\delta(r^{+})$
$\displaystyle\int_{{\bf P},{\bf Q}}e^{i{\bf P}\cdot{\bf
r}}\frac{\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}\left(P_{q}^{-},{\bf
P};P_{q}^{-},-\frac{\textbf{Q}}{2}\right)\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}\left(P_{q}^{-},\frac{\textbf{Q}}{2};P_{q}^{-},{\bf
P}\right)}{2\left({\bf P}+\frac{\bf Q}{2}\right)^{2}\left({\bf P}-\frac{\bf
Q}{2}\right)^{2}}$ $\displaystyle\times$
$\displaystyle\int_{\textbf{B}}e^{i{\bf Q}\cdot({\bf b}-{\bf
B})}\mu^{2}(b^{+},\textbf{B}).$ (40)

Here, we drop the hat notation for field operator from now on.

Further simplification is possible under the assumption of a homogeneous
nucleus, $\mu(\vec{x})=\mu^{2}$. Introducing the transverse color charge
density $\tilde{\mu}^{2}=\int_{x^{+}}\mu^{2}(x^{+})=L^{+}\mu^{2}$, where
$L^{+}$ represents the light-cone longitudinal width of the nucleus, Section
III.3 can be simplified to give

$\displaystyle\big{\langle}A^{\mu}\left(\vec{x}\right)A^{\nu}\left(\vec{y}\right)\big{\rangle}_{A}=\delta^{ab}\delta(x^{+}-y^{+})\frac{\tilde{\mu}^{2}}{L^{+}}G^{\mu\nu}({\bf
x}-{\bf y}),$ (41)

where we define

$\displaystyle G^{\mu\nu}({\bf r})=\int_{{\bf P}}e^{i{\bf P}\cdot{\bf
r}}\frac{\Gamma^{\mu}_{\sigma_{1}\sigma_{2}}\left(P_{q}^{-},{\bf
P};P_{q}^{-},0\right)\Gamma^{\nu}_{\sigma_{2}\sigma_{1}}\left(P_{q}^{-},0;P_{q}^{-},{\bf
P}\right)}{2{\bf P}^{4}}.$ (42)

Equation 41 is a key result of this work, as it allows us to calculate the
correlation of the classical gluon background field generated by a large
nucleus composed of quasi-independent quarks. In the following sections, we
will leverage this equation to compute the nuclear DIS cross section at sub-
eikonal accuracy.

Prior to proceeding with our analysis, let us examine the scaling behavior of
the correlator with respect to the nucleus rapidity $\omega$,
$\gamma=e^{\omega}$. As established in the preceding section, the partonic
longitudinal momentum scales with the nucleus rapidity as
$P_{q}^{-}\propto\gamma$. Therefore, we can deduce the following scaling with
$\gamma$ from Eq. 26:

$\displaystyle G^{--}\propto 1,\qquad G^{-i}\propto\gamma^{-1},\qquad
G^{-+},G^{ij}\propto\gamma^{-2},\qquad G^{+i}\propto\gamma^{-3},\qquad
G^{++}\propto\gamma^{-4}.$ (43)

As we will see in the next section, at sub-eikonal level, i.e., at
$\mathcal{O}(\gamma^{-1})$, the only components of the field correlator that
contribute are $--$, $-i$, $i-$ and $ij$. In this case, the non-trivial
components of Eq. 42 that we need are the following:

$\displaystyle G^{--}({\bf r})$ $\displaystyle=\int_{{\bf P}}e^{i{\bf
P}\cdot{\bf r}}\frac{1}{{\bf P}^{4}},$ (44) $\displaystyle G^{i-}({\bf r})$
$\displaystyle=\frac{1}{2P_{q}^{-}}\int_{{\bf P}}e^{i{\bf P}\cdot{\bf
r}}\frac{{\bf P}^{i}}{{\bf P}^{4}},$ (45) $\displaystyle G^{ij}({\bf r})$
$\displaystyle=\frac{1}{(2P_{q}^{-})^{2}}\int_{{\bf P}}e^{i{\bf P}\cdot{\bf
r}}\frac{{\bf P}^{i}{\bf P}^{j}+\epsilon^{im}\epsilon^{jn}{\bf P}^{m}{\bf
P}^{n}}{{\bf P}^{4}}.$ (46)

On the other hand, due to Lorentz contraction, $L^{+}\propto\gamma^{-1}$ and
$\delta(x^{+}-y^{+})\propto\gamma$. This leads to the following scaling
behavior for the color field correlators:

$\displaystyle\langle A_{a}^{-}(x^{+},{\bf x})A_{b}^{-}(y^{+},{\bf
y})\rangle_{A}\propto\gamma^{2},\qquad\langle A_{a}^{-}(x^{+},{\bf x}){\bf
A}_{b}^{j}(y^{+},{\bf y})\rangle_{A}\propto\gamma,\qquad\langle{\bf
A}_{a}^{i}(x^{+},{\bf x}){\bf A}_{b}^{j}(y^{+},{\bf
y})\rangle_{A}\propto\gamma^{0}.$ (47)

This observed scaling aligns with our expectations. Under a boost in the left
direction, the gauge field transforms as
$A^{\mu}=\Lambda^{\mu}_{\nu}A_{0}^{\nu}=(\gamma
A_{0}^{-},\gamma^{-1}A_{0}^{+},{\bf A}_{0}^{i})$, with $A^{\mu}_{0}$ the
background field in the rest frame of the nucleus.

##  IV DIS dijet production beyond the eikonal approximation

In this section, we review the calculation of dijet production in nuclear deep
inelastic scattering (DIS) beyond the eikonal approximation within the CGC
framework, presented in [55]. We employ the dipole picture of high-energy
nuclear DIS, where the production of a forward quark-antiquark dijet arises
from the splitting of a virtual photon $\gamma^{*}$ into a quark-antiquark
dipole followed by its scattering with the nucleus. We work in a frame where
the target’s transverse momentum vanishes. The incoming electron has momentum
$k_{e}^{\mu}$, colliding with a nucleus with momentum $P_{A}^{\mu}$. The
interaction is mediated by a virtual photon with momentum
$q^{\mu}=k_{e}^{\prime\mu}-k_{e}^{\mu}$, where $k_{e}^{\prime\mu}$ represents
the final momentum of the electron. The virtual photon subsequently decays
into a quark-antiquark pair with individual momenta $k_{1}^{\mu}$ and
$k_{2}^{\mu}$, respectively. The quark and antiquark interact with the
classical background field, $A^{\mu}(x)$, generated by the color sources
within the nucleus.

Within the single-photon exchange approximation, the DIS process is described
by the product of a leptonic tensor, characterizing the virtual photon
emission by the electron, and a hadronic tensor, characterizing the virtual
photon’s interaction with the target. After integrating over the azimuthal
angle of the scattered electron, the hadronic tensor is projected onto the
cross section of the $\gamma^{*}+A$ interaction:

$\displaystyle\frac{d\sigma^{e^{-}+A\to e+q\bar{q}+X}}{dW^{2}dQ^{2}d^{2}{\bf
k}_{1}d^{2}{\bf k}_{2}d\eta_{1}d\eta_{2}}=\frac{\alpha_{\rm em}}{\pi
Q^{2}sy}\left[\left(1-y+\frac{y^{2}}{2}\right)\frac{d\sigma^{\gamma^{*}_{T}+A\to
q\bar{q}+X}}{d^{2}{\bf k}_{1}d^{2}{\bf
k}_{2}d\eta_{1}d\eta_{2}}+(1-y)\frac{d\sigma^{\gamma^{*}_{L}+A\to
q\bar{q}+X}}{d^{2}{\bf k}_{1}d^{2}{\bf k}_{2}d\eta_{1}d\eta_{2}}\right].$ (48)

Here, $y=\frac{W^{2}+Q^{2}-m_{\rm N}^{2}}{s-m_{\rm N}^{2}}$ is the
inelasticity, $m_{\rm N}$ is the nucleon mass, $Q^{2}=-q^{2}$ is the photon
virtuality, $\sqrt{s}$ is the center-of-mass (CoM) energy per nucleon of the
$e^{-}+A$ system and $W$ is the CoM energy per nucleon of the $\gamma^{*}+A$
system.

The cross-section for $q\bar{q}$ production in a $\gamma^{*}+A$ collision,
extending beyond the eikonal approximation, can be expressed as follows:

$\displaystyle\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d^{2}{\bf
k}_{2}d\eta_{1}d\eta_{2}}=\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to
q\bar{q}+X}}{d^{2}{\bf k}_{1}d^{2}{\bf k}_{2}d\eta_{1}d\eta_{2}}\Bigg{|}_{\rm
Eik.}+\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d^{2}{\bf k}_{2}d\eta_{1}d\eta_{2}}\Bigg{|}_{\rm
NEik.}+\mathcal{O}(\gamma^{-2}),$ (49)

where $\lambda=L,T$ denotes the polarization of the virtual photon, either
longitudinal ($L$) or transverse ($T$). The first term represents the eikonal
approximation, which assumes an infinitely large rapidity difference between
the virtual photon and the target. The second term introduces the
$\gamma^{-1}$ order correction, denoted as NEik., which accounts for the first
order deviations from the eikonal approximation.

Figure 2: Effective eikonal amplitude diagram for $q\bar{q}$ production in
nuclear DIS. A virtual photon with polarization $\lambda$ and momentum
$q^{\mu}$ splits into a quark-antiquark pair. The color dipole interacts with
the nuclear target via light-like Wilson lines. This eikonal interaction is
represented by the blue rectangle.

Figure 2 shows the effective diagram for calculating the scattering amplitude
for $q\bar{q}$ production in $\gamma^{*}+A$ collisions within the eikonal
approximation. The incoming right-moving virtual photon undergoes splitting
into a quark-antiquark dipole. This splitting process is analytically
calculable using perturbative Quantum Electrodynamics (QED). On the other
hand, the subsequent multiple interactions of the dipole with the background
field generated by the target require non-perturbative treatment and are
expressed in terms of light-like Wilson lines in the fundamental
representation of SU($N_{c}$):

$\displaystyle\mathcal{U}_{[x^{+},y^{+}]}({\bf
z})=\mathcal{P}_{+}\exp\left\\{-ig\int_{y^{+}}^{x^{+}}dz^{+}A^{-}(z^{+},{\bf
z})\right\\},$ (50)

where $\mathcal{P}_{+}$ is the path order operator and, due to the infinite
rapidity difference between the projectile and target, only the longitudinal
component of the field contributes to the scattering. Equation 50 represents
the eikonal propagation of a fermion in a gluonic background field from
longitudinal point $y^{+}$ to $x^{+}$.

Henceforth, by combining these elements, the cross section reads [59]

$\displaystyle\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d^{2}{\bf k}_{2}d\eta_{1}d\eta_{2}}\Bigg{|}_{\rm Eik.}$
$\displaystyle=\int_{{\bf v},{\bf v}^{\prime},{\bf w},{\bf
w}^{\prime}}e^{i{\bf k}_{1}\cdot({\bf v}^{\prime}-{\bf v})+i{\bf
k}_{2}\cdot({\bf w}^{\prime}-{\bf w})}\mathcal{C}_{\lambda}({\bf
w}^{\prime}-{\bf v}^{\prime},{\bf w}-{\bf v})$ $\displaystyle\hskip
28.45274pt\times\Big{[}Q({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf
w})-d({\bf w}^{\prime},{\bf v}^{\prime})-d({\bf v},{\bf w})+1\Big{]},$ (51)

where

$\displaystyle\mathcal{C}_{L}({\bf r}_{1},{\bf r}_{2})$
$\displaystyle=\sum_{f}\frac{8N_{c}\alpha_{\rm
em}e_{f}^{2}Q^{2}}{(2\pi)^{6}}\delta_{z}z_{1}^{3}z_{2}^{3}K_{0}(\epsilon_{f}|{\bf
r}_{1}|)K_{0}(\epsilon_{f}|{\bf r}_{2}|),$ (52)
$\displaystyle\mathcal{C}_{T}({\bf r}_{1},{\bf r}_{2})$
$\displaystyle=\sum_{f}\frac{2N_{c}\alpha_{\rm
em}e_{f}^{2}}{(2\pi)^{6}}\delta_{z}z_{1}z_{2}\left[m_{f}^{2}+(z_{1}^{2}+z_{2}^{2})\partial_{{\bf
r}_{1}^{j}}\partial_{{\bf r}_{2}^{j}}\right]K_{0}(\epsilon_{f}|{\bf
r}_{1}|)K_{0}(\epsilon_{f}|{\bf r}_{2}|).$ (53)

Here, we introduce the notation $\delta_{z}\equiv\delta(z_{1}+z_{2}-1)$, where
the variables $z_{1,2}=k_{1,2}^{+}/q^{+}$ represent the fraction of +-momentum
of the photon carried by the (anti)quark. We also define
$\epsilon_{f}^{2}=m_{f}^{2}+z_{1}z_{2}Q^{2}$, where $m_{f}$ denotes the mass
of the quark with flavor $f$ and fractional charge $e_{f}$. Section IV is
factorized into two components: the factor $\mathcal{C}_{\lambda}$, which,
aside from constant factors, represents the probability of splitting, and the
dipole and quadrupole functions

$\displaystyle d({\bf v},{\bf w})$
$\displaystyle=\frac{1}{N_{c}}\Big{\langle}{\rm Tr}\left[\mathcal{U}({\bf
v})\mathcal{U}^{\dagger}({\bf w})\right]\Big{\rangle},\qquad\hskip
11.38092ptQ({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf
w})=\frac{1}{N_{c}}\Big{\langle}{\rm Tr}\left[\mathcal{U}({\bf
w}^{\prime})\mathcal{U}^{\dagger}({\bf v}^{\prime})\mathcal{U}({\bf
v})\mathcal{U}^{\dagger}({\bf w})\right]\Big{\rangle},$ (54)

which account for the eikonal interaction of the quark and antiquark with the
nucleus. In this paper, we adopt the convention that when the particle
propagates through the entire extent $L^{+}$ of the target (as is always
occurs in the eikonal limit due to Lorentz contraction), the Wilson line is
written as $\mathcal{U}({\bf x})\equiv\mathcal{U}_{[L^{+}/2,-L^{+}/2]}({\bf
x})$.

At next-to-eikonal (NEik) accuracy, the classical background field loses its
infinite contraction and acquires a finite light-cone width, denoted by
$L^{+}$. Furthermore, the components $A^{+}$ and ${\bf A}^{i}$ along with the
$x^{-}$ dependence of the field become non-zero. As previously mentioned, this
manuscript excludes non-eikonal corrections stemming from both the $x^{-}$
dependence and the $A^{+}$ component of the field, as the latter only
contributes at next-to-next-to-eikonal accuracy. In this case, the field
correlators derived in Eq. 41 within the covariant gauge
($\partial_{\mu}A^{\mu}=0$) are equally valid in the light-cone gauge
($A^{+}=0$). This equivalence eliminates the need for gauge transformations
when employing the correlators, as the results of [55] were obtained in the
light-cone gauge. Consequently, our study solely concentrates on next-to-
eikonal corrections arising from the finite target width and the non-vanishing
transverse components of the field.

Figure 3: Effective NEik amplitude diagrams for $q\bar{q}$ production in
nuclear deep inelastic scattering (DIS). The diagrams depict scenarios where
the virtual photon splits into a quark-antiquark pair before entering the
nucleus. After splitting, either the quark or antiquark (illustrated only for
the quark here) propagates eikonally in the background field up to a
longitudinal point $v^{+}$. At this point, insertions of the transverse
covariant derivative or strength tensor occurs. Subsequently, the particle
propagates eikonally again.

Due to the medium’s finite width, the photon can split either inside or before
encountering it. In Fig. 3, we show a subset of effective diagrams
contributing to the NEik amplitude when the photon splits before entering the
medium. The blue rectangle represents the eikonal scattering of the partons
(either quark or antiquark) with the medium. The non-eikonal corrections arise
from inserting either the transverse covariant derivative $D_{{\bf
z}^{j}}(v^{+})={\partial_{{\bf z}^{j}}}-ig{\bf A}^{j}(v^{+},{\bf z})$, or the
transverse field strength tensor $F_{ij}(z^{+},{\bf z})=\partial_{{\bf
z}^{i}}{\bf A}_{j}(z^{+},{\bf z})-\partial_{{\bf z}^{j}}{\bf A}_{i}(z^{+},{\bf
z})-ig[{\bf A}_{i}(z^{+},{\bf z}),{\bf A}_{j}(z^{+},{\bf z})]$, at a
longitudinal point $v^{+}$ along the parton’s path. On the other hand, Fig. 4
shows the case where the photon splits inside the target. In this case, the
only type of next-to-eikonal corrections appears through a single transverse
covariant derivative insertion in the longitudinal point where the photon
splits.

Figure 4: Effective NEik amplitude diagram for $q\bar{q}$ production in
nuclear deep inelastic scattering (DIS). The diagram depicts the scenario
where the virtual photon splits into a quark-antiquark pair inside the
nucleus.

All in all, the non-eikonal corrections originating from these diagrams have
been computed in [55]. The final result for the NEik corrections to the cross
section for $q\bar{q}$ production in a $\gamma_{\lambda}^{*}+A$ collision is:

$\displaystyle\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d^{2}{\bf k}_{2}d\eta_{1}d\eta_{2}}\Bigg{|}_{\rm
NEik.}=\frac{1}{q^{+}}{\rm Re}\int_{{\bf v},{\bf w},{\bf v}^{\prime},{\bf
w}^{\prime}}e^{i{\bf k}_{1}\cdot({\bf v}^{\prime}-{\bf v})+i{\bf
k}_{2}\cdot({\bf w}^{\prime}-{\bf w})}\mathcal{C}_{\lambda}({\bf
w}^{\prime}-{\bf v}^{\prime},{\bf w}-{\bf v})$
$\displaystyle\times\Bigg{\\{}\frac{1}{z_{1}}\left[\frac{{\bf k}_{2}^{j}-{\bf
k}_{1}^{j}}{2}+\frac{i}{2}\partial_{{\bf w}^{j}}\right]\left[Q_{j}^{(1)}({\bf
w}^{\prime},{\bf v}^{\prime},{\bf v}_{*},{\bf w})-d_{j}^{(1)}({\bf v}_{*},{\bf
w})\right]-\frac{i}{z_{1}}\left[Q^{(2)}({\bf w}^{\prime},{\bf v}^{\prime},{\bf
v}_{*},{\bf w})-d^{(2)}({\bf v}_{*},{\bf w})\right]$ $\displaystyle\hskip
8.0pt-\frac{1}{z_{2}}\left[\frac{{\bf k}_{2}^{j}-{\bf
k}_{1}^{j}}{2}-\frac{i}{2}\partial_{{\bf v}^{j}}\right]\left[Q_{j}^{(1)}({\bf
v}^{\prime},{\bf w}^{\prime},{\bf w}_{*},{\bf v})^{\dagger}-d_{j}^{(1)}({\bf
w}_{*},{\bf v})^{\dagger}\right]-\frac{i}{z_{2}}\left[Q^{(2)}({\bf
v}^{\prime},{\bf w}^{\prime},{\bf w}_{*},{\bf v})^{\dagger}-d^{(2)}({\bf
w}_{*},{\bf v})^{\dagger}\right]\Bigg{\\}}$ $\displaystyle+\delta^{\lambda
T}\frac{d\sigma^{\rm trans.}}{d\Pi},$ (55)

where we have introduced the factor $\delta^{\lambda T}$ which is 1 when
$\lambda=T$ and 0 when $\lambda=L$, since the function $\frac{d\sigma^{\rm
trans.}}{d\Pi}$, to be defined below, only contributes to the cross section
when the photon is transversely polarized.

Section IV gives the $\mathcal{O}(\gamma^{-1})$ correction to the eikonal
cross section given in Section IV. Indeed, we can see that the only parameter
in Section IV that scales with $\gamma$ under a boost is $q^{+}$, leading to
the expect $\gamma^{-1}$ scaling. Now, the non-eikonal interaction of the
quark and antiquark with the background field is accounted for by the
following modifications of the dipole and quadrupole functions:

$\displaystyle d_{j}^{(1)}({\bf v}_{*},{\bf w})$
$\displaystyle=\frac{1}{N_{c}}\Big{\langle}{\rm
Tr}\left[\mathcal{U}_{j}^{(1)}({\bf v})\mathcal{U}^{\dagger}({\bf
w})\right]\Big{\rangle},\qquad Q_{j}^{(1)}({\bf w}^{\prime},{\bf
v}^{\prime},{\bf v}_{*},{\bf w})=\frac{1}{N_{c}}\Big{\langle}{\rm
Tr}\left[\mathcal{U}({\bf w}^{\prime})\mathcal{U}^{\dagger}({\bf
v}^{\prime})\mathcal{U}_{j}^{(1)}({\bf v})\mathcal{U}^{\dagger}({\bf
w})\right]\Big{\rangle},$ (56) $\displaystyle d^{(2)}({\bf v}_{*},{\bf w})$
$\displaystyle=\frac{1}{N_{c}}\Big{\langle}{\rm
Tr}\left[\mathcal{U}^{(2)}({\bf v})\mathcal{U}^{\dagger}({\bf
w})\right]\Big{\rangle},\qquad Q^{(2)}({\bf w}^{\prime},{\bf v}^{\prime},{\bf
v}_{*},{\bf w})=\frac{1}{N_{c}}\Big{\langle}{\rm Tr}\left[\mathcal{U}({\bf
w}^{\prime})\mathcal{U}^{\dagger}({\bf v}^{\prime})\mathcal{U}^{(2)}({\bf
v})\mathcal{U}^{\dagger}({\bf w})\right]\Big{\rangle},$ (57)

where

$\displaystyle\mathcal{U}_{j}^{(1)}({\bf z})$
$\displaystyle=\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dv^{+}\mathcal{U}_{[\frac{L^{+}}{2},v^{+}]}({\bf
z})\overleftrightarrow{D_{{\bf
z}^{j}}}(v^{+})\mathcal{U}_{[v^{+},-\frac{L^{+}}{2}]}({\bf z}),$ (58)
$\displaystyle\mathcal{U}^{(2)}({\bf z})$
$\displaystyle=\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dv^{+}\mathcal{U}_{[\frac{L^{+}}{2},v^{+}]}({\bf
z})\overleftarrow{D_{{\bf z}^{j}}}(v^{+})\overrightarrow{D_{{\bf
z}^{j}}}(v^{+})\mathcal{U}_{[v^{+},-\frac{L^{+}}{2}]}({\bf z}).$ (59)

In the context of next-to-eikonal corrections, Eqs. 58 and 59 describe the
propagation of a fermion through a medium with finite width, as depicted in
the first and second diagrams of Fig. 3, respectively. Both equations capture
the fermion’s interactions with the background field but differ in the nature
of the non-eikonal interactions at point $v^{+}$.

In Eq. 58, the fermion undergoes eikonal propagation from the medium’s
beginning ($-L^{+}/2$) to $v^{+}$. At $v^{+}$, a non-eikonal interaction
occurs through the insertion of $\overleftrightarrow{D_{{\bf
z}^{j}}}(z^{+})={\partial_{{\bf z}^{j}}}-\overleftarrow{\partial_{{\bf
z}^{j}}}-2ig{\bf A}^{j}(z^{+},{\bf z})$. Thus, the fermion interacts with the
transverse gradient of the longitudinal field and the transverse field itself.
Following this, the fermion experiences another eikonal propagation until
reaching the medium’s end ($L^{+}/2$). Equation 59 describes a similar
scenario with a distinct non-eikonal interaction. Here, the interaction is
represented by $\overleftarrow{D_{{\bf z}^{j}}}(v^{+})\overrightarrow{D_{{\bf
z}^{j}}}(v^{+})$, where $\overleftarrow{D_{{\bf
z}^{j}}}(z^{+})=\overleftarrow{\partial_{{\bf z}^{j}}}+ig{\bf
A}^{j}(z^{+},{\bf z})$. In this case the fermion can interact with two
gradients of the field as well as two transverse field insertions at the same
longitudinal point.

Moreover, the part of the cross section that only contributes for a
transversely polarized photon is given by

$\displaystyle\frac{d\sigma^{\rm trans.}}{d\Pi}$
$\displaystyle=\frac{1}{q^{+}}{\rm Re}\int_{{\bf v},{\bf w},{\bf
v}^{\prime},{\bf w}^{\prime}}e^{i{\bf k}_{1}\cdot({\bf v}^{\prime}-{\bf
v})+i{\bf k}_{2}\cdot({\bf w}^{\prime}-{\bf w})}\mathcal{D}_{T,1}^{ij}({\bf
w}^{\prime}-{\bf v}^{\prime},{\bf w}-{\bf v})$ $\displaystyle\hskip
28.45274pt\times\Bigg{[}\frac{i}{z_{1}}\left[Q_{ij}^{(3)}({\bf
w}^{\prime},{\bf v}^{\prime},{\bf v}_{*},{\bf w})-d_{ij}^{(3)}({\bf
v}_{*},{\bf w})\right]+\frac{i}{z_{2}}\left[Q_{ij}^{(3)}({\bf v}^{\prime},{\bf
w}^{\prime},{\bf w}_{*},{\bf v})^{\dagger}-d_{ij}^{(3)}({\bf w}_{*},{\bf
v})^{\dagger}\right]\Bigg{]}$ $\displaystyle-\frac{1}{q^{+}}{\rm Im}\int_{{\bf
z},{\bf v}^{\prime},{\bf w}^{\prime}}e^{i{\bf k}_{1}\cdot({\bf
v}^{\prime}-{\bf z})+i{\bf k}_{2}\cdot({\bf w}^{\prime}-{\bf
z})}\mathcal{D}_{T,2}^{j}({\bf w}^{\prime}-{\bf v}^{\prime})$
$\displaystyle\hskip
28.45274pt\times\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dz^{+}\Bigg{\langle}\frac{1}{N_{c}}{\rm
Tr}\Big{[}\mathcal{U}({\bf w}^{\prime})\mathcal{U}^{\dagger}({\bf
v}^{\prime})-1\Big{]}\Big{[}\mathcal{U}_{\left[\frac{L^{+}}{2},z^{+}\right]}({\bf
z})\overleftrightarrow{D_{{\bf
z}^{j}}}(z^{+})\mathcal{U}_{\left[\frac{L^{+}}{2},z^{+}\right]}^{\dagger}({\bf
z})-\frac{1}{2}\mathcal{U}({\bf z})\overleftrightarrow{\partial_{{\bf
z}^{j}}}\mathcal{U}^{\dagger}({\bf z})\Big{]}\Bigg{\rangle},$ (60)

where

$\displaystyle\mathcal{D}_{T,1}^{ij}({\bf r}_{1},{\bf
r}_{2})=\sum_{f}\frac{2N_{c}\alpha_{\rm
em}e_{f}^{2}}{(2\pi)^{6}}\delta_{z}z_{1}z_{2}(z_{1}-z_{2})\partial_{{\bf
r}_{1}^{i}}\partial_{{\bf r}_{2}^{j}}K_{0}(\epsilon_{f}|{\bf
r}_{1}|)K_{0}(\epsilon_{f}|{\bf r}_{2}|),$ (61)
$\displaystyle\mathcal{D}_{T,2}^{j}({\bf r})=\sum_{f}\frac{N_{c}\alpha_{\rm
em}e_{f}^{2}}{2(2\pi)^{5}}\delta_{z}[1+(z_{2}-z_{1})^{2}]\partial_{{\bf
r}^{j}}K_{0}(\epsilon_{f}|{\bf r}|).$ (62)

In first term of Section IV, the quark and antiquark pair is generated before
entering the medium and their interaction with the background field is
accounted by the following modifications of the dipole and quadrupole
operators:

$\displaystyle d_{ij}^{(3)}({\bf v}_{*},{\bf w})$
$\displaystyle=\frac{1}{N_{c}}\Big{\langle}{\rm
Tr}\left[\mathcal{U}_{ij}^{(3)}({\bf v})\mathcal{U}^{\dagger}({\bf
w})\right]\Big{\rangle},\qquad Q_{ij}^{(3)}({\bf w}^{\prime},{\bf
v}^{\prime},{\bf v}_{*},{\bf w})=\frac{1}{N_{c}}\Big{\langle}{\rm
Tr}\left[\mathcal{U}({\bf w}^{\prime})\mathcal{U}^{\dagger}({\bf
v}^{\prime})\mathcal{U}_{ij}^{(3)}({\bf v})\mathcal{U}^{\dagger}({\bf
w})\right]\Big{\rangle},$ (63)

where

$\displaystyle\mathcal{U}_{ij}^{(3)}({\bf z})$
$\displaystyle=\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dv^{+}\mathcal{U}_{[\frac{L^{+}}{2},v^{+}]}({\bf
z})gF_{ij}(v^{+},{\bf z})\mathcal{U}_{[v^{+},-\frac{L^{+}}{2}]}({\bf z}).$
(64)

Equation 64 defines a modified Wilson line, capturing the case where the non-
eikonal corrections arise from the insertion of the transverse strength
tensor, as shown in the third diagram of Fig. 3. The second term of Section IV
represents the scenario where the photon splits inside the medium, as depicted
in Fig. 4.

##  V DIS dijet production in the dilute limit beyond the eikonal
approximation

This section aims to calculate the cross section for dijet production in
nuclear DIS beyond the eikonal approximation using the CGC framework given in
Sections IV and IV. Specifically, we will employ the dilute limit
approximation, which offers significant computational advantages while
remaining applicable under certain regions of the kinematic space. In
particular, when the transverse momenta of the $q\bar{q}$ pair exceeds the
saturation scale of the target nucleus, $Q_{s}$, the interaction regime
switches from dense to dilute. In this scenario, the gluon density within the
nucleus becomes sufficiently low that individual gluon scatterings with the
$q\bar{q}$ pair dominate, rather than collective interactions with the entire
Color Glass Condensate.

The main advantage for utilizing the dilute limit approximation is that the
physics become perturbative and we can approximate the Wilson line as

$\displaystyle\mathcal{U}_{[x^{+},y^{+}]}({\bf
z})=1-ig\int_{y^{+}}^{x^{+}}dz^{+}A^{-}(z^{+},{\bf
z})-g^{2}\int_{y^{+}}^{x^{+}}dz_{1}^{+}\int_{y^{+}}^{z_{1}^{+}}dz_{2}^{+}A^{-}(z_{1}^{+},{\bf
z})A^{-}(z_{2}^{+},{\bf z})+\mathcal{O}(g^{3}),$ (65)

so that the interaction of the partons with the medium can be approximated by
the exchange of two gluons instead of a multi-gluon interaction. This allows
us to directly use the non-eikonal field correlator derived in Section III.3.
In contrast, exploring the dense limit would necessitate: (i) develop a method
to resum the non-eikonal correlators, accounting for the multi-gluon exchange
contributions that become essential in the dense regime, and (ii) take into
account the non-linear terms of the classical YM equation.

Therefore, due to the substantial increase in complexity associated with the
dense limit beyond the eikonal approximation, in this paper we study the cross
section in the kinematic regime where ${\bf k}_{1}^{2},{\bf k}_{2}^{2}\gtrsim
Q_{s}^{2}$. In this case the dilute limit approximation is valid and presents
a more computationally efficient and tractable approach while maintaining
sufficient accuracy for our purposes.

###  V.1 DIS dijet production at eikonal accuracy in the dilute limit

Our calculation begins with determining the cross section at eikonal accuracy.
To achieve this, we require a model for the dipole and quadrupole operators
defined in Eq. 54. Leveraging the expansion of the Wilson line from Eq. 65, we
express the dipole function at second order in the coupling constant as
follows:

$\displaystyle d({\bf x},{\bf y})$
$\displaystyle=1-\frac{g^{2}}{N_{c}}\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}z^{+}\int_{-\frac{L^{+}}{2}}^{z^{+}}dz_{1}^{+}$
$\displaystyle\hskip 56.9055pt\times{\rm Tr}\Big{[}\langle A^{-}(z^{+},{\bf
x})A^{-}(z_{1}^{+},{\bf x})\rangle+\langle A^{-}(z_{1}^{+},{\bf
y})A^{-}(z^{+},{\bf y})\rangle-2\langle A^{-}(z^{+},{\bf
x})A^{-}(z_{1}^{+},{\bf y})\rangle\Big{]}.$ (66)

Next, we employ the expression for the field correlator provided in Eq. 41 to
write

$\displaystyle\frac{g^{2}}{N_{c}}\big{\langle}{\rm
Tr}\left[A^{\mu}\left(\vec{x}\right)A^{\nu}\left(\vec{y}\right)\right]\big{\rangle}=2\pi
Q_{s}^{2}\frac{\delta(x^{+}-y^{+})}{L^{+}}G^{\mu\nu}({\bf x}-{\bf y}).$ (67)

Here, we have defined the saturation momentum as

$\displaystyle
Q_{s}^{2}\equiv\frac{g^{2}}{2\pi}C_{F}\tilde{\mu}^{2}\approx\frac{g^{4}AC_{F}}{4\pi^{2}R_{A}^{2}},$
(68)

where in the last step we used the fact the for a homogeneous circular
nucleus, the transverse color charge density squared is
$\tilde{\mu}^{2}\approx g^{2}A/(2\pi R_{A}^{2})$, where $R_{A}\sim A^{1/3}$ fm
is the radius of the nucleus.

By using the above expression, Eq. (V.1) can be written as follows666We note
that by resumming higher order correlators, we arrive at the well known
expression within the MV model:  $\displaystyle d({\bf x},{\bf
y})=\exp\left\\{-\frac{Q_{s}^{2}}{4}({\bf x}-{\bf y})^{2}\ln\frac{1}{|{\bf
x}-{\bf y}|\Lambda_{\rm QCD}}\right\\},$ (69) when $|{\bf x}-{\bf y}|\lesssim
Q_{s}$. :

$\displaystyle d({\bf x},{\bf y})=1-2\pi Q_{s}^{2}\left[G^{--}(0)-G^{--}({\bf
x}-{\bf y})\right].$ (70)

Similarly, in the dilute limit of the CGC, the quadrupole function is
expressed as

$\displaystyle Q({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf w})=1-2\pi
Q_{s}^{2}\Big{[}$ $\displaystyle 2G^{--}(0)-G^{--}({\bf v}^{\prime}-{\bf
v})-G^{--}({\bf v}-{\bf w})+G^{--}({\bf v}^{\prime}-{\bf w})$
$\displaystyle+G^{--}({\bf w}^{\prime}-{\bf v})-G^{--}({\bf w}^{\prime}-{\bf
v}^{\prime})-G^{--}({\bf w}^{\prime}-{\bf w})\Big{]}.$ (71)

All in all, the eikonal term of the cross section, presented in Section IV,
reads

$\displaystyle\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm Eik.}$
$\displaystyle=2\pi Q_{s}^{2}\int_{{\bf v},{\bf v}^{\prime},{\bf w},{\bf
w}^{\prime}}e^{i{\bf k}_{1}\cdot({\bf v}^{\prime}-{\bf v})+i{\bf
k}_{2}\cdot({\bf w}^{\prime}-{\bf w})}\mathcal{C}_{\lambda}({\bf
w}^{\prime}-{\bf v}^{\prime},{\bf w}-{\bf v})$
$\displaystyle\times\Big{[}G^{--}({\bf v}^{\prime}-{\bf v})-G^{--}({\bf
v}^{\prime}-{\bf w})-G^{--}({\bf w}^{\prime}-{\bf v})+G^{--}({\bf
w}^{\prime}-{\bf w})\Big{]}.$ (72)

It is clear that this equation solely accounts for the interaction of the
quark (or antiquark) with a single gluon in both the amplitude and complex
amplitude, as depicted in Fig. 5.

Figure 5: One of the four diagrams contributing to the eikonal differential
cross section for quark-antiquark production in DIS in the dilute limit. The
circle with a cross represents a longitudinal source.

Unlike in the dense limit, Section V.1 can be solved exactly and it does not
depend on the non-perturbative scale $\Lambda_{\rm QCD}$. This can be
explicitly seen by performing the following change of variables:

$\displaystyle{\bf v}={\bf B}-\frac{{\bf r}+{\bf R}}{2},\qquad{\bf
v}^{\prime}={\bf B}-\frac{{\bf b}-{\bf R}}{2},\qquad{\bf w}={\bf B}+\frac{{\bf
r}-{\bf R}}{2},\qquad{\bf w}^{\prime}={\bf B}+\frac{{\bf b}+{\bf
R}}{2},\qquad$ (73)

which simplifies Section V.1 to

$\displaystyle\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm Eik.}=2\pi
Q_{s}^{2}S_{\perp}\int_{{\bf b},{\bf r}}\left(e^{i{\bf b}\cdot{\bf
k}_{2}}-e^{-i{\bf b}\cdot{\bf k}_{1}}\right)\left(e^{-i{\bf r}\cdot{\bf
k}_{2}}-e^{i{\bf r}\cdot{\bf k}_{1}}\right)\mathcal{C}_{\lambda}({\bf r},{\bf
b})\int_{\bf R}e^{i{\bf R}\cdot({\bf k}_{1}+{\bf k}_{2})}G^{--}({\bf R}).$
(74)

Here, the factor $S_{\perp}\equiv\int d^{2}{\bf B}$ arises, since we are
assuming a homogeneous and therefore translationally invariant target.
Equation 74 can be solved analytically using the definition of $G^{--}$ in Eq.
44 as well as the following identities:

$\displaystyle\int_{\bf x}e^{i{\bf x}\cdot{\bf k}}K_{0}(\epsilon_{f}|{\bf
x}|)=2\pi\frac{1}{{\bf k}^{2}+\epsilon_{f}^{2}},\qquad\int_{\bf x}e^{i{\bf
x}\cdot{\bf k}}\partial_{{\bf x}^{j}}K_{0}(\epsilon_{f}|{\bf x}|)=-2\pi
i\frac{{\bf k}^{j}}{{\bf k}^{2}+\epsilon_{f}^{2}}.$ (75)

By solving the transverse integrals using Eq. 75, we obtain the expression for
the differential cross-section for quark and antiquark production at eikonal
accuracy:

$\displaystyle\frac{d\sigma^{\gamma^{*}_{L}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm Eik.}$
$\displaystyle=\sum_{f}C^{(f)}_{L}\frac{Q^{2}({\bf k}_{1}^{2}-{\bf
k}_{2}^{2})^{2}}{({\bf k}_{1}+{\bf k}_{2})^{4}({\bf
k}_{1}^{2}+\epsilon_{f}^{2})^{2}({\bf k}_{1}^{2}+\epsilon_{f}^{2})^{2}},$ (76)
$\displaystyle\frac{d\sigma^{\gamma^{*}_{T}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm Eik.}$
$\displaystyle=\sum_{f}C^{(f)}_{T}\frac{1}{({\bf k}_{1}+{\bf
k}_{2})^{4}}\left[m_{f}^{2}\frac{({\bf k}_{1}^{2}-{\bf k}_{2}^{2})^{2}}{({\bf
k}_{2}^{2}+\epsilon_{f}^{2})^{2}({\bf
k}_{1}^{2}+\epsilon_{f}^{2})^{2}}+(z_{1}^{2}+z_{2}^{2})\left(\frac{{\bf
k}_{2}^{i}}{{\bf k}_{2}^{2}+\epsilon_{f}^{2}}+\frac{{\bf k}_{1}^{i}}{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}\right)^{2}\right],$ (77)

where we have defined the following dimensionless parameters:

$\displaystyle C^{(f)}_{L}=\frac{8N_{c}\alpha_{\rm
em}e_{f}^{2}Q_{s}^{2}S_{\perp}}{(2\pi)^{3}}\delta_{z}z_{1}^{3}z_{2}^{3},\qquad
C^{(f)}_{T}=\frac{2N_{c}\alpha_{\rm
em}e_{f}^{2}Q_{s}^{2}S_{\perp}}{(2\pi)^{3}}\delta_{z}z_{1}z_{2}.$ (78)

It follows, then, that in the dilute regime of CGC, the eikonal differential
cross-section for dijet production in DIS is analytical. Moreover, unlike in
the dense limit within the MV model, the infrared regulator $\Lambda_{\rm
QCD}$ is not required. While the results presented in Eqs. 76 and 77 thus do
include non-perturbative modes, this does not significantly affect the results
provided that ${\bf k}_{1}^{2},{\bf k}_{2}^{2}\gtrsim Q_{s}^{2}\gg\Lambda_{\rm
QCD}^{2}$.

###  V.2 DIS dijet production beyond eikonal accuracy in the dilute limit

We now turn to calculating the next-to-eikonal corrections to the differential
cross-section using the same approach as before. As outlined in Section IV,
non-eikonal corrections can be categorized into four distinct types based on
how the parton interacts with the medium. If the splitting occurs before the
photon enters the medium, these corrections include: (i) a single covariant
derivative insertion, (ii) a double covariant derivative insertion, and (iii)
a transverse strength field insertion. We also need to consider the
contribution when (iv) the splitting occurs inside the medium. In the
following, we compute all these contributions in the dilute limit of the CGC.

####  V.2.1 Single covariant derivative contribution

The case in which the parton interacts with the background field with a single
transverse covariant derivative insertion is accounted for by the modification
of the Wilson line given in Eq. 58. In the dilute limit, it can be written as
follows:

$\displaystyle\mathcal{U}_{j}^{(1)}({\bf v})$
$\displaystyle=\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dv^{+}\Bigg{[}-ig\int_{-\frac{L^{+}}{2}}^{v^{+}}dz_{1}^{+}\left(\partial_{{\bf
v}^{j}}A^{-}(z_{1}^{+},{\bf v})-\partial_{{\bf v}^{j}}A^{-}(v^{+},{\bf
v})\right)-2ig{\bf A}^{j}(v^{+},{\bf v})$
$\displaystyle+g^{2}\int_{-\frac{L^{+}}{2}}^{v^{+}}dz_{1}^{+}\int_{-\frac{L^{+}}{2}}^{z_{1}^{+}}dz_{2}^{+}\partial_{{\bf
v}^{j}}\left(A^{-}(v^{+},{\bf v})A^{-}(z_{1}^{+},{\bf v})-A^{-}(z_{1}^{+},{\bf
v})A^{-}(z_{2}^{+},{\bf v})\right)$
$\displaystyle-2g^{2}\int_{-\frac{L^{+}}{2}}^{v^{+}}dz_{1}^{+}\Big{(}{\bf
A}^{j}(v^{+},{\bf v})A^{-}(z_{1}^{+},{\bf v})+A^{-}(v^{+},{\bf v}){\bf
A}^{j}(z_{1}^{+},{\bf v})\Big{)}\Bigg{]}+\mathcal{O}(g^{3}).$ (79)

Using this expansion, it is straightforward to compute in the dilute limit the
modifications of the dipole and quadrupole operators given in Eq. 56:

$\displaystyle d_{j}^{(1)}({\bf v}_{*},{\bf w})$ $\displaystyle=4\pi
Q_{s}^{2}\left[G^{-j}({\bf w}-{\bf v})-G^{j-}(0)\right],$ (80) $\displaystyle
Q_{j}^{(1)}({\bf w}^{\prime},{\bf v}^{\prime},{\bf v}_{*},{\bf w})$
$\displaystyle=4\pi Q_{s}^{2}\left[G^{-j}({\bf w}-{\bf v})+G^{-j}({\bf
v}^{\prime}-{\bf v})-G^{-j}({\bf w}^{\prime}-{\bf v})-G^{j-}(0)\right].$ (81)

In conclusion, we note that the next-to-eikonal corrections in Section IV
solely arise from the difference between the modified quadrupole and dipole
functions:

$\displaystyle Q_{j}^{(1)}({\bf w}^{\prime},{\bf v}^{\prime},{\bf v}_{*},{\bf
w})-d_{j}^{(1)}({\bf v}_{*},{\bf w})=4\pi Q_{s}^{2}\left[G^{-j}({\bf
v}^{\prime}-{\bf v})-G^{-j}({\bf w}^{\prime}-{\bf v})\right].$ (82)

This type of non-eikonal correction, independent of the gradient of the
longitudinal field component, describes the scenario where a parton interacts
with the transverse field in the amplitude, while the partons in the complex
amplitude interact with the longitudinal component. Figure 6 illustrates this,
where the red curly line represents the transverse gluon and the black line
represents the longitudinal gluon. We only show the case where the quark
interacts with the transverse component, but there are two additional diagrams
where the antiquark interacts with the transverse field.

Figure 6: Two of the four diagrams contributing to the next-to-eikonal
differential cross section for quark-antiquark production in DIS with a single
insertion. The red curly line represents the exchange of a transverse gluon,
while the black line represents a longitudinal gluon. The red square
represents a transverse source of color charge, and the circle with a cross
represents a longitudinal source. The diagrams with the antiquark interacting
with the transverse gluon are not shown but are analogous to the ones
presented here.

####  V.2.2 Double covariant derivative contribution

The case where the parton interacts with the background field through a double
transverse covariant derivative insertion is accounted for by the modification
of the Wilson line given in Eq. 59. In the dilute limit, this modified Wilson
line can be written as follows:

$\displaystyle\mathcal{U}^{(2)}({\bf v})$
$\displaystyle=-g^{2}\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dv^{+}\Bigg{[}\int_{-\frac{L^{+}}{2}}^{v^{+}}dz_{1}^{+}\int_{-\frac{L^{+}}{2}}^{z_{1}^{+}}dz_{2}^{+}\partial_{{\bf
v}^{j}}A^{-}(z_{1}^{+},{\bf v})\partial_{{\bf v}^{j}}A^{-}(z_{2}^{+},{\bf v})$
$\displaystyle+\int_{-\frac{L^{+}}{2}}^{v^{+}}dz_{1}^{+}\Big{(}\partial_{{\bf
v}^{j}}A^{-}(v^{+},{\bf v}){\bf A}^{j}(z_{1}^{+},{\bf v})-{\bf
A}^{j}(v^{+},{\bf v})\partial_{{\bf v}^{j}}A^{-}(z_{1}^{+},{\bf
v})\Big{)}-{\bf A}^{j}(v^{+},{\bf v}){\bf A}^{j}(v^{+},{\bf
v})\Bigg{]}+\mathcal{O}(g^{3}).$ (83)

Since $\mathcal{U}^{(2)}\propto g^{2}$, the modified dipole and quadrupole
functions defined in Eq. 57 are given by

$\displaystyle Q^{(2)}({\bf w}^{\prime},{\bf v}^{\prime},{\bf v}_{*},{\bf
w})=d^{(2)}({\bf v}_{*},{\bf
w})+\mathcal{O}(g^{3})=\langle\mathcal{U}^{(2)}({\bf
v})\rangle+\mathcal{O}(g^{3}).$ (84)

This implies that these operators do not contribute to the NEik differential
cross-section in the dilute limit:

$\displaystyle Q^{(2)}({\bf w}^{\prime},{\bf v}^{\prime},{\bf v}_{*},{\bf
w})-d^{(2)}({\bf v}_{*},{\bf w})=0.$ (85)

####  V.2.3 Transverse strength tensor contribution

The case where the parton interacts with the background field through a
transverse strength tensor insertion is accounted for by the modification of
the Wilson line given in Eq. 64. In the dilute limit, this modified Wilson
line can be written as follows:

$\displaystyle\mathcal{U}_{ij}^{(3)}({\bf v})$
$\displaystyle=\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dv^{+}\Bigg{[}g\Big{(}\partial_{{\bf
v}^{j}}{\bf A}^{i}(v^{+},{\bf v})-\partial_{{\bf v}^{i}}{\bf A}^{j}(v^{+},{\bf
v})-ig[{\bf A}^{i}(v^{+},{\bf v}),{\bf A}^{j}(v^{+},{\bf v})]\Big{)}$ (86)
$\displaystyle\hskip-28.45274pt-
ig^{2}\int_{-\frac{L^{+}}{2}}^{v^{+}}dz_{1}^{+}\Big{(}A^{-}(v^{+},{\bf
v})[\partial_{{\bf v}^{j}}{\bf A}^{i}(z_{1}^{+},{\bf v})-\partial_{{\bf
v}^{i}}{\bf A}^{j}(z_{1}^{+},{\bf v})]+[\partial_{{\bf v}^{j}}{\bf
A}^{i}(v^{+},{\bf v})-\partial_{{\bf v}^{i}}{\bf A}^{j}(v^{+},{\bf
v})]A^{-}(z_{1}^{+},{\bf v})\Big{)}\Bigg{]}+\mathcal{O}(g^{3}).$

Up to order $g^{2}$, all terms contributing to the modification of the dipole
and quadrupole functions (Eq. (63)) involve correlators like

$\displaystyle\langle[\partial_{{\bf v}^{j}}{\bf A}^{i}(v^{+},{\bf
v})-\partial_{{\bf v}^{i}}{\bf A}^{j}(v^{+},{\bf v})]A^{-}(w^{+},{\bf
w})\rangle$ $\displaystyle\propto\partial_{{\bf v}^{j}}G^{i-}({\bf v}-{\bf
w})-\partial_{{\bf v}^{i}}G^{j-}({\bf v}-{\bf w})$ $\displaystyle=i\int_{\bf
P}e^{i{\bf P}\cdot({\bf x}-{\bf y})}\frac{{\bf P}^{i}{\bf P}^{j}-{\bf
P}^{j}{\bf P}^{i}}{{\bf P}^{4}}=0,$ (87)

where in the last step we have used the definition of $G^{i-}$ given in Eq.
45, and

$\displaystyle\langle[{\bf A}^{i}(v^{+},{\bf v}),{\bf A}^{j}(v^{+},{\bf
v})]\rangle\propto G^{ij}(0)-G^{ji}(0)=0,$ (88)

where we have used the fact that $G^{ij}$, given by Eq. 46, is symmetric in
the transverse indices.

Therefore, in the dilute limit and for a homogeneous target, corrections
stemming from the field strength tensor insertion do not contribute to the
cross section because

$\displaystyle Q_{ij}^{(3)}({\bf w}^{\prime},{\bf v}^{\prime},{\bf v}_{*},{\bf
w})=d_{ij}^{(3)}({\bf v}_{*},{\bf w})=0.$ (89)

We note that the vanishing result in Section V.2.3 relies on the assumption of
a homogeneous nucleus. For an inhomogeneous target with a non-trivial density
distribution, these corrections are expected to become non-negligible. In
fact, the magnitude of their contribution likely increases with the degree of
inhomogeneity in the density distribution. Further investigation of such
inhomogeneous scenarios, using more complex models for $\mu^{2}(b^{+},{\bf
b})$, could be crucial for achieving a more accurate understanding of these
corrections and their impact on the overall calculations.

####  V.2.4 Splitting inside contribution

Finally, we now analyze the corrections stemming from the case in which the
photon splits in the $q{\bar{q}}$ dipole in a coordinate $z^{+},{\bf z}$
inside the target. In this case we need to evaluate the following correlator:

$\displaystyle Q_{\rm in}({\bf w}^{\prime},{\bf v}^{\prime},{\bf
z})\equiv\int_{-\frac{L^{+}}{2}}^{\frac{L^{+}}{2}}dz^{+}\Bigg{\langle}\frac{1}{N_{c}}\Big{[}\mathcal{U}({\bf
w}^{\prime})\mathcal{U}^{\dagger}({\bf
v}^{\prime})-1\Big{]}\Big{[}\mathcal{U}_{\left[\frac{L^{+}}{2},z^{+}\right]}({\bf
z})\overleftrightarrow{D_{{\bf
z}^{j}}}(z^{+})\mathcal{U}_{\left[\frac{L^{+}}{2},z^{+}\right]}^{\dagger}({\bf
z})-\frac{1}{2}\mathcal{U}({\bf z})\overleftrightarrow{\partial_{{\bf
z}^{j}}}\mathcal{U}^{\dagger}({\bf z})\Big{]}\Bigg{\rangle}$ (90)

Expanding the Wilson lines, it is possible to see that the contributions
coming from the longitudinal field gradient vanish. Thus, after performing the
longitudinal integrals this operator can be written as

$\displaystyle Q_{\rm in}({\bf w}^{\prime},{\bf v}^{\prime},{\bf z})=4\pi
Q_{s}^{2}\left[G^{-j}({\bf v}^{\prime}-{\bf z})-G^{-j}({\bf w}^{\prime}-{\bf
z})\right].$ (91)

We see that this is the same result obtained in Eq. 82 but with ${\bf v}={\bf
w}={\bf z}$.

Figure 7: Diagrams contributing to the next-to-eikonal differential cross
section for quark-antiquark production in DIS when the photon splits inside
the target. The red curly line represents the exchange of a transverse gluon,
while the black line represents a longitudinal gluon. The red square
represents a transverse source of color charge, and the circle with a cross
represents a longitudinal source.

####  V.2.5 Summing up the contributions

In summary, for a homogeneous target within the dilute limit, the next-to-
eikonal corrections to the differential cross section are solely determined by
the single covariant derivative contribution (Eq. 82) and the contribution
from in-medium splitting (Eq. 106). These contributions are depicted in the
diagrams presented in Figs. 6 and 7.

By substituting Eqs. 82 and 106 into the NEik differential cross section
definition (Sections IV and IV) and utilizing the change of variables in Eq.
73, we arrive at the following expression for the cross section:

$\displaystyle\frac{d\sigma^{\gamma^{*}_{\lambda}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm NEik.}=\frac{2\pi
Q_{s}^{2}S_{\perp}}{q^{+}}\int_{{\bf b},{\bf r},{\bf R}}e^{i{\bf R}\cdot({\bf
k}_{1}+{\bf k}_{2})}\mathcal{C}_{\lambda}({\bf r},{\bf b})$
$\displaystyle\times\Bigg{[}e^{-i{\bf r}\cdot{\bf k}_{2}}\left(e^{i{\bf
b}\cdot{\bf k}_{2}}-e^{-i{\bf b}\cdot{\bf k}_{1}}\right)\frac{({\bf
k}_{2}-{\bf k}_{1})^{j}}{z_{1}}G^{-j}({\bf R})+e^{-i{\bf r}\cdot{\bf
k}_{1}}\left(e^{i{\bf b}\cdot{\bf k}_{1}}-e^{-i{\bf b}\cdot{\bf
k}_{2}}\right)\frac{({\bf k}_{1}-{\bf k}_{2})^{j}}{z_{2}}G^{-j}({\bf
R})^{*}\Bigg{]}$ $\displaystyle\hskip 85.35826pt-\delta^{\lambda T}\frac{4\pi
Q_{s}^{2}S_{\perp}}{q^{+}}{\rm Im}\int_{{\bf b},{\bf R}}e^{i{\bf R}\cdot({\bf
k}_{1}+{\bf k}_{2})}\mathcal{D}_{T,2}^{j}({\bf b})\left(e^{i{\bf b}\cdot{\bf
k}_{2}}-e^{-i{\bf b}\cdot{\bf k}_{1}}\right)G^{-j}({\bf R}).$ (92)

Finally, solving this integral is straightforward, by using the definition of
$G^{-i}$ given in Eq. 45 together with the Bessel integral given in Eq. 75.
Analogous to Eqs. 76 and 77, the result is also analytical and given by

$\displaystyle\frac{d\sigma^{\gamma^{*}_{L}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm NEik.}$
$\displaystyle=\frac{C_{T}}{2q^{+}P_{q}^{-}}\frac{Q^{2}({\bf k}_{1}^{2}-{\bf
k}_{2}^{2})^{2}}{({\bf k}_{1}+{\bf k}_{2})^{4}({\bf
k}_{1}^{2}+\epsilon_{f}^{2})^{2}({\bf
k}_{1}^{2}+\epsilon_{f}^{2})^{2}}\left[\frac{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}{z_{1}}-\frac{{\bf
k}_{2}^{2}+\epsilon_{f}^{2}}{z_{2}}\right],$ (93)
$\displaystyle\frac{d\sigma^{\gamma^{*}_{T}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm NEik.}$
$\displaystyle=\frac{C_{T}}{2q^{+}P_{q}^{-}}\frac{{\bf k}_{1}^{2}-{\bf
k}_{2}^{2}}{({\bf k}_{1}+{\bf k}_{2})^{4}}\Bigg{\\{}m_{f}^{2}\frac{{\bf
k}_{1}^{2}-{\bf k}_{2}^{2}}{({\bf k}_{1}^{2}+\epsilon_{f}^{2})^{2}({\bf
k}_{1}^{2}+\epsilon_{f}^{2})^{2}}\Bigg{[}\frac{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}{z_{1}}-\frac{{\bf
k}_{2}^{2}+\epsilon_{f}^{2}}{z_{2}}\Bigg{]}$ $\displaystyle\hskip
0.0pt+(z_{1}^{2}+z_{2}^{2})\left(\frac{{\bf k}_{2}^{i}}{{\bf
k}_{2}^{2}+\epsilon_{f}^{2}}+\frac{{\bf k}_{1}^{i}}{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}\right)\Bigg{[}\frac{1}{z_{1}}\frac{{\bf
k}_{2}^{i}}{{\bf k}_{2}^{2}+\epsilon_{f}^{2}}+\frac{1}{z_{2}}\frac{{\bf
k}_{1}^{i}}{{\bf k}_{1}^{2}+\epsilon_{f}^{2}}+\frac{1}{z_{1}z_{2}}\frac{{\bf
k}_{1}^{i}+{\bf k}_{2}^{i}}{{\bf k}_{1}^{2}-{\bf
k}_{2}^{2}}\Bigg{]}\Bigg{\\}}.$ (94)

These equations represent the next-to-eikonal, $\mathcal{O}(\gamma^{-1})$,
corrections to the differential cross section for dijet production in
$\gamma^{*}+A$ collisions in the dilute limit of the nucleus. Indeed, these
corrections are suppressed by $q^{+}P_{q}^{-}\propto e^{\Delta Y}$, where
$\Delta Y$ is the rapidity difference between the quarks inside the nucleus
and the photon. Moreover, it is important to note that since Eqs. 93 and 94
were computed in the dilute approximation, they are only valid in the regime
where ${\bf k}_{1}^{2},{\bf k}_{2}^{2}\gtrsim Q_{s}^{2}$.

###  V.3 The differential cross section at subeikonal order

To finalize the section, we sum the differential cross section at eikonal
accuracy given in Eqs. 76 and 77 with their NEik corrections given in Eqs. 93
and 94. This sum represents the subeikonal differential cross section, which
incorporates the NEik corrections to the eikonal result. Moreover, we express
the suppression in the NEik contribution in terms of the COM energy per
nucleon of the $\gamma^{*}+A$ system. In order to do so, we note that for a
energetic nucleus we can write $W^{2}=(q+P_{\rm N})^{2}\approx
2N_{c}q^{+}P_{q}^{-}$, where we have used the fact that all the quarks in the
nucleon have the same longitudinal momentum, $P_{\rm N}^{-}=N_{c}P_{q}^{-}$.
Therefore, the subeikonal differential cross section for a longitudinally and
transversely polarized photon can be written as follows:

$\displaystyle\frac{d\sigma^{\gamma^{*}_{L}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}$
$\displaystyle=\frac{d\sigma^{\gamma^{*}_{L}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm
Eik.}\Bigg{[}1+\frac{N_{c}}{W^{2}}\left(\frac{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}{z_{1}}-\frac{{\bf
k}_{2}^{2}+\epsilon_{f}^{2}}{z_{2}}\right)\Bigg{]}+\mathcal{O}(\gamma^{-2}),$
(95) $\displaystyle\frac{d\sigma^{\gamma^{*}_{T}+A\to q\bar{q}+X}}{d^{2}{\bf
k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}$
$\displaystyle=\sum_{f}\frac{C_{T}^{(f)}}{({\bf k}_{1}+{\bf
k}_{2})^{4}}\Bigg{\\{}m_{f}^{2}\frac{({\bf k}_{1}^{2}-{\bf
k}_{2}^{2})^{2}}{({\bf k}_{2}^{2}+\epsilon_{f}^{2})^{2}({\bf
k}_{1}^{2}+\epsilon_{f}^{2})^{2}}\Bigg{[}1+\frac{N_{c}}{W^{2}}\left(\frac{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}{z_{1}}-\frac{{\bf
k}_{2}^{2}+\epsilon_{f}^{2}}{z_{2}}\right)\Bigg{]}$
$\displaystyle\hskip-56.9055pt+(z_{1}^{2}+z_{2}^{2})\left(\frac{{\bf
k}_{2}^{i}}{{\bf k}_{2}^{2}+\epsilon_{f}^{2}}+\frac{{\bf k}_{1}^{i}}{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}\right)\Bigg{(}\frac{{\bf k}_{2}^{i}}{{\bf
k}_{2}^{2}+\epsilon_{f}^{2}}\Bigg{[}1+\frac{N_{c}}{W^{2}}\frac{{\bf
k}_{1}^{2}-{\bf k}_{2}^{2}}{z_{1}}\Bigg{]}+\frac{{\bf k}_{1}^{i}}{{\bf
k}_{1}^{2}+\epsilon_{f}^{2}}\Bigg{[}1+\frac{N_{c}}{W^{2}}\frac{{\bf
k}_{1}^{2}-{\bf k}_{2}^{2}}{z_{2}}\Bigg{]}$ $\displaystyle\hskip
113.81102pt+\frac{N_{c}}{W^{2}}\frac{{\bf k}_{1}^{i}+{\bf
k}_{2}^{i}}{z_{1}z_{2}}\Bigg{)}\Bigg{\\}}+\mathcal{O}(\gamma^{-2}).$ (96)

These equations are the main result of this section and describe the
differential cross section for $q{\bar{q}}$ production in DIS at subeikonal
accuracy when the final transverse momentum of the partons are relatively
large, ${\bf k}_{1}^{2},{\bf k}_{2}^{2}\gtrsim Q_{s}^{2}$. In the next
section, we will analyze these solutions numerically to gain further insights
into the dynamics of DIS in the subeikonal regime.

##  VI Numerical results

In this section, we perform a numerical analysis of the subeikonal
differential cross sections, Eqs. 95 and 96, derived in the previous section,
and we contrast them with the eikonal approximation. Within this analysis, the
target nucleus is chosen as Au, with $A=197$ nucleons and characterized by a
nuclear radius $R_{A}=6.37$ fm which translates into a transverse area of
$S_{\perp}=\pi R_{A}^{2}$. A saturation momentum $Q_{s}=0.6$ GeV is taken.
Given the present study’s focus on qualitative insights on the effect of non-
eikonal corrections, the small-$x$ evolution of the dipole and quadrupole
operators is not considered.

We focus on massless, $m_{f}=0$, light quarks. This choice implies
$\epsilon_{f}^{2}=z_{1}z_{2}Q^{2}$, thereby eliminating any quark flavor
dependence in the cross section except for the factor
$\sum_{f}e_{f}^{2}=\left(\frac{2}{3}\right)^{2}+2\left(\frac{1}{3}\right)^{2}$.
Additionally, the CoM energy per nucleon of the collision is fixed to be
$\sqrt{s}=90$ GeV, which corresponds to the top EIC energy for collisions on
nuclear targets. We utilize the standard values for the number of colors
($N_{c}=3$) and the fine structure constant ($\alpha_{\rm em}=1/137$).

The pseudorapidity ($\eta_{1,2}$) and longitudinal momentum fraction
($z_{1,2}$) of the dijet partons are linked by $z_{1,2}=\frac{2E_{\rm
N}}{W^{2}}|\textbf{k}_{1,2}|e^{\eta_{1,2}},$ where $E_{\rm N}$ is the beam
energy per nucleon. Given that $z_{2}=1-z_{1}$, we investigate the asymmetry
between $z_{1}$ and $z_{2}$ defining

$\xi=\ln\frac{1-z_{1}}{z_{1}}=\ln\frac{|\textbf{k}_{2}|}{|\textbf{k}_{1}|}+\Delta\eta,$
(97)

where $\Delta\eta=\eta_{2}-\eta_{1}$ is the pseudorapidity difference. This
variable takes the value $\xi=0$ for symmetric longitudinal momentum fractions
($z_{1}=z_{2}=0.5$), and approaches $\pm\infty$ when $z_{1,2}=0$.

Furthermore, due to rotational symmetry, the cross section depends only on the
azimuthal angle between the quark and antiquark defined as
$\cos\phi_{q\bar{q}}={\bf k}_{1}\cdot{\bf k}_{2}/k_{1\perp}k_{2\perp}$, with
$k_{1,2\perp}\equiv|{\bf k}_{1,2}|$. Hence, integrating over $z_{2}$ using the
Dirac delta function, the phase space measure adopts the form $d\Pi=2\pi
k_{1\perp}dk_{1\perp}k_{2\perp}dk_{2\perp}d\phi_{q\bar{q}}d\xi$.

Figure 8: Differential cross section as a function of the angle between the
quark and antiquark, $\phi_{q\bar{q}}$, at fixed values of $k_{1\perp}=2$ GeV,
$k_{2\perp}=3$ GeV, and $\xi=0$. Left: Transverse (red) and longitudinal
(black) differential cross sections for the $\gamma^{*}+\text{Au}$ system.
Right: Differential cross section per unit of $W^{2}$ and $Q^{2}$ for
$e^{-}+\text{Au}$ collisions. Results are shown for different values of
$Q^{2}$ and $W$, with the solid lines representing the eikonal result and the
dashed lines representing the eikonal+subeikonal one.

Figure 8 showcases the impact of subeikonal corrections on the dijet angular
distribution. We quantitatively compare the eikonal and (eikonal+)subeikonal
differential cross sections as a function of the dijet angle
($\phi_{q\bar{q}}$) at fixed values of transverse momenta ($k_{1\perp}=2$ GeV,
$k_{2\perp}=3$ GeV) and longitudinal momentum asymmetry ($\xi=0$). This
analysis is conducted for various combinations of center-of-mass energy
($W^{2}$) and momentum transfer ($Q^{2}$). Additionally, we quantify the
subeikonal corrections in the cross section for the $e^{-}+A$ system defined
in Eq. 48.

Our observations reveal that the back-to-back peak centered around
$\phi_{q\bar{q}}=\pi$ exhibits a broader profile in the very inelastic
scattering regime ($W\approx\sqrt{s}$) compared to a less inelastic case
($W\ll\sqrt{s}$). Furthermore, for larger $Q^{2}$ values (corresponding to a
higher inelasticity $y$), a double-peak structure emerges. This shows that the
dominance of the back-to-back peak diminishes as momentum transfer increases.
Notably, the non-eikonal corrections become more pronounced for smaller $W$
values, substantially modifying the back-to-back peak of the DIS cross
section. Given the magnitude of these modifications, we anticipate that
subeikonal corrections will hold sizable importance in future experimental
measurements at the Electron-Ion Collider.

###  VI.1 Azimuthal anisotropy

Motivated by its relevance in the correlation limit of the Color Glass
Condensate and to facilitate comparison with existing literature, we describe
quark-antiquark momenta using the momentum imbalance,
$\mathbb{\Delta}=\mathbb{k}_{1}+\mathbb{k}_{2}$, and the relative momentum,
$\mathbb{P}=z_{2}\mathbb{k}_{1}-z_{1}\mathbb{k}_{2}$. $\mathbb{P}$ relates to
the dijet pair’s invariant mass as
$M^{2}_{q\bar{q}}=P_{\perp}^{2}/(z_{1}z_{2})$, where $P_{\perp}=|{\bf P}|$,
$\Delta_{\perp}=|{\bf\Delta}|$, and $\phi$ is the azimuthal angle difference
between $\mathbb{P}$ and $\mathbb{\Delta}$. In these variables, the phase
space measure simplifies to $d\Pi=2\pi
P_{\perp}dP_{\perp}\Delta_{\perp}d\Delta_{\perp}d\phi d\xi$.

For the $\gamma^{*}+A$ system777For the $e^{-}+A$ system, the harmonics
expansion is analogous, with the key difference that the angle-averaged cross
section is given in units of $Q^{2}$ and $W^{2}$. We denote this by
$\tilde{D}_{0}(P_{\perp},\Delta_{\perp})$., the differential cross section can
be expressed as a harmonics expansion

$\displaystyle\frac{d\sigma_{\lambda}^{\gamma^{*}A\to q\bar{q}X}}{d\Pi}$
$\displaystyle=D_{0,\lambda}(P_{\perp},\Delta_{\perp})\left[1+2\sum_{n=1}^{\infty}v_{n,\lambda}(P_{\perp},\Delta_{\perp})\cos\phi\right],$
(98)

where

$\displaystyle D_{n,\lambda}(P_{\perp},\Delta_{\perp})$
$\displaystyle=\int\frac{d\phi}{2\pi}e^{in\phi}\frac{d\sigma^{\gamma_{\lambda}^{*}A\to
q\bar{q}X}}{d\Pi},\qquad v_{n,\lambda}=\frac{D_{n,\lambda}}{D_{0,\lambda}}.$
(99)

Leveraging the expressions derived in Section V, we analyze the angle-averaged
differential cross section, $D_{0,\lambda}$, and the azimuthal harmonics
$v_{n,\lambda}$ as functions of $P_{\perp}$, $\Delta_{\perp}$ and $\xi$.

Figure 9: Angle-averaged dijet cross section as a function of the fraction of
longitudinal momentum asymmetry, $\xi$, at fixed value of $P_{\perp}=3$ GeV
and $\Delta_{\perp}=2$ GeV. Left: Transverse (red) and longitudinal (black)
differential cross sections for the $\gamma^{*}+\text{Au}$ system; Right:
Cross section for $e^{-}+{\rm Au}$ collisions. Results are shown for different
values of $Q^{2}$ and $W$, with the solid line representing the eikonal result
and the dashed lines including the non-eikonal corrections.

Figure 9 shows the dependence of the angle-averaged cross section
$D_{0,\lambda}$ with respect to $\xi$ at both eikonal and subeikonal accuracy.
Fixed values of $P_{\perp}=3$ GeV, $\Delta_{\perp}=2$ GeV, and $Q^{2}=10$ GeV
are employed. The first row displays results for relatively low inelasticity
($y=0.11$), while the second row shows the case for high inelasticity
($y=0.98$). We present the results for both $\gamma^{*}+\text{Au}$ and
$e^{-}+\text{Au}$ collisions in separate panels.

We observe that the eikonal cross section is symmetric with respect to $\xi$,
while the subeikonal corrections break this symmetry. As expected, non-eikonal
corrections are larger for smaller $W$ values but vanish when $\xi=0$. This
behavior can be analyzed analytically through Eq. 93, where we see that the
non-eikonal terms, when $z_{1}=z_{2}$, are proportional to ${\bf
k}_{1}^{2}-{\bf k}_{2}^{2}\propto\cos\phi$ and therefore only modify odd
harmonics.

Figure 10 shows the dependence of the angle-averaged reduced cross section per
unit of $W^{2}$ and $Q^{2}$, $\tilde{D}_{0,\lambda}$, for $e^{-}+{\rm Au}$
collisions, on $P_{\perp}$ and $\Delta_{\perp}$, at eikonal and subeikonal
accuracy. The non-eikonal corrections exhibit a strong dependence on
$P_{\perp}$, but only a weak variation with respect to $\Delta_{\perp}$. For
$P_{\perp}\gtrsim 2.5$ GeV, $\Delta_{\perp}=1$ GeV, and $\xi=1$, the non-
eikonal corrections reach approximately 10 %, while for $P_{\perp}=1$ GeV and
$\xi=1$, they amount to $\sim 5\%$, irrespective of the $\Delta_{\perp}$
value.

Figure 10: Left: Angle averaged dijet cross section as a function of
$P_{\perp}$ at a fixed value of $\Delta_{\perp}=1$ GeV and $\xi=1$. Right:
Angle averaged dijet cross section as a function of $\Delta_{\perp}$ at a
fixed value of $P_{\perp}=2$ GeV and $\xi=1$. The solid line represents the
eikonal result while the dashed line includes the non-eikonal corrections.

Figure 11 depicts the dependence of the elliptic ($v_{2}$) and triangular
($v_{3}$) harmonics on the longitudinal momentum asymmetry, $\xi$, at fixed
$P_{\perp}=3$ GeV and $\Delta_{\perp}=2$ GeV. Similar to Fig. 9, the non-
eikonal corrections break the (anti)symmetry of the even (odd) azimuthal
harmonics with respect to $\xi$. Notably, the subeikonal cross section
exhibits non-zero odd harmonics even when the longitudinal momentum fractions
are symmetric ($\xi=0$), unlike the eikonal case. Importantly, in contrast to
what is found in the correlation limit (see Section A.1), this approach
reveals a dependence of both odd and even harmonics on the longitudinal
momentum asymmetry.

Figure 11: Left: Dependence of the elliptic harmonic, $v_{2}$, with respect to
the longitudinal momenta asymmetry, $\xi$. Right: Dependence of the triangular
harmonic, $v_{3}$, with respect to $\xi$. Both figures are for a fixed
$P_{\perp}=3$ GeV and $\Delta_{\perp}=2$ GeV. We denote the harmonics for the
$e^{-}+{\rm Au}$ collision as ”total”.

Figure 12 illustrates the dependence of the elliptic ($v_{2}$) and triangular
($v_{3}$) harmonics on the transverse momentum transfer, $P_{\perp}$, at fixed
$\Delta_{\perp}=2$ GeV and $\xi=1$. Consistent with Fig. 10, significant
corrections arise from non-eikonal effects at high $P_{\perp}$ values. As
expected, both $v_{2}$ and $v_{3}$ exhibit a saturation behavior at large
momenta.

Furthermore, Fig. 13 depicts the odd azimuthal harmonic as a function of
$P_{\perp}$ and $\Delta_{\perp}$ for symmetric longitudinal momentum fractions
($\xi=0$). Unlike the eikonal case, non-eikonal corrections generate non-zero
odd harmonics, reaching approximately 5%. Measuring these odd harmonics in
nuclear DIS experiments therefore provides a direct probe of the non-eikonal
contributions.

Figure 12: Left: Dependence of the elliptic harmonic, $v_{2}$, with respect to
$P_{\perp}$. Right: Dependence of the triangular harmonic, $v_{3}$, with
respect to $P_{\perp}$. Both figures are for fixed $\Delta_{\perp}=2$ GeV and
$\xi=1$. We denote the harmonics for the $e^{-}+{\rm Au}$ collision as
”total”. Figure 13: Left: Dependence of the triangular harmonic, $v_{3}$, with
respect to $P_{\perp}$ at $\Delta_{\perp}=2$ GeV. Right: Dependence of the
triangular harmonic, $v_{3}$, with respect to $\Delta_{\perp}$ at
$P_{\perp}=3$ GeV. Both plots are for the case in which the quark-antiquark
longitudinal momentum are the same, $\xi=0$. We denote the harmonics for the
$e^{-}+{\rm Au}$ collision as ”total”.

###  VI.2 Comparison with the dense limit and the correlation limit

To finalize, this section presents a comprehensive comparison of the eikonal
angle-averaged cross section. We analyze the expression derived within the
dilute limit of the CGC framework (Eqs. 76 and 77) from Section V. These
expressions are compared to the dense limit introduced in Appendix A and,
owing to its analytical simplicity, to the correlation limit
($\Delta_{\perp}\ll P_{\perp}$) described in Section A.1.

This comparative analysis serves two primary objectives: (i) it offers a
stringent validation of the expressions derived within this manuscript, and
(ii) it establishes quantitatively the kinematic regime of $P_{\perp}$ and
$\Delta_{\perp}$ wherein the methodology employed in this study remains
applicable.

Figure 14 visually depicts the aforementioned comparison of the angle-averaged
cross section in the dilute limit with the dense and correlation limits. In
order to make a better comparison with the approach in Section V, where the
inclusion of a non-perturbative parameter was not required, we use a small
value of $\Lambda_{\rm QCD}=0.1$ GeV. As expected, all approaches converge in
the limit $P_{\perp}\gg\Delta_{\perp}\gg Q_{s}$. Notably, the dilute limit
exhibits a closer agreement with the dense limit than the correlation limit
across a wider region of the kinematic space. This observation is particularly
noteworthy considering the inherent analytical simplicity of the dilute limit
expression. This finding suggests that the proposed approach offers a
computationally efficient tool for calculations within a broader kinematic
range compared to other limits.

Figure 14: Left: Eikonal angle-averaged dijet cross section as a function of
$\Delta_{\perp}$ and $P_{\perp}=1,2,3$ GeV. Right: Eikonal angle-averaged
dijet cross section as a function of $P_{\perp}$ and $\Delta_{\perp}=1,2,3$
GeV. In all the figures we take $\xi=0$ and compare the dilute limit,
represented by a continuous line, with the dense limit (dotted line) and its
correlation limit (dashed line).

##  VII Conclusions

In this work we have started the study of the effects of non-eikonal
corrections in the CGC on observables relevant for nuclear DIS, particularly
at the energies of the future EIC. We have analyzed the effect of non-eikonal
corrections on dijet production in DIS off nuclear targets. We have considered
the next-to-eikonal corrections for this observable derived in [55]. They
require the knowledge of correlators of fields in the target beyond those
computed in the MV model, specifically those between transverse and boost-
enhanced (minus) components, and of the recoil ($x^{-}$-dependence) of the
field. We have ignored the latter, while for the former we have developed a
linear model valid for large nuclei, which reproduces the MV model results for
the correlator of boost-enhanced components and generalizes the one in [36]
going beyond spin exchange contributions.

We have considered the unpolarized cross sections for dijet production in the
approximation of a homogenous dilute nucleus, obtaining simple analytic
expressions for the cross sections in the eikonal and next-to-eikonal
approximations valid in the limit of total dijet momentum $P_{\perp}$ and
dijet momentum imbalance $\Delta_{\perp}$ larger than the saturation scale
$Q_{s}$, $P_{\perp},\Delta_{\perp}\gtrsim Q_{s}$.

We have performed a numerical study of the results applied to
$e^{-}+\text{Au}$ collisions at top EIC energies. We have found
$\mathcal{O}(10\%)$ effects in the cross sections, particularly at large
$P_{\perp}$, thus making the inclusion of non-eikonal corrections necessary
for precision studies at the EIC. We have analyzed the azimuthal asymmetries
between $\mathbf{P}$ and $\mathbf{\Delta}$, $P_{\perp}=|{\bf P}|$,
$\Delta_{\perp}=|{\bf\Delta}|$. We have found that non-eikonal corrections
induce odd azimuthal harmonics for the symmetric situation (i.e., jets with
equal momentum fractions from the virtual photon) where they are absent in the
eikonal approximation, therefore offering a clean signature of non-eikonal
effects.

Finally, in the eikonal approximation we have compared the results of our
analytic expansion valid in the dilute limit of the target, and the full CGC
results in the MV model and their correlation limit [59]. We have shown that
our analytic expressions match the correlation limit ones in the region where
both should be simultaneously valid, $P_{\perp}\gg\Delta_{\perp}\gtrsim
Q_{s}$. We have also observed that our simple analytic formulae reproduce very
well the full CGC results in its validity region
$P_{\perp},\Delta_{\perp}\gtrsim Q_{s}$.

Given the focus of this paper solely on nuclear DIS, one crucial future
direction is the extension to the case of DIS on light nuclei or protons. This
extension necessitates the inclusion of longitudinal correlations and two-body
operators, neglected in Section II. Furthermore, when
$P_{\perp},\Delta_{\perp}<Q_{s}$ multi-gluon exchanges become more
significant, requiring the application of resummation techniques.
Consequently, exploring the dense limit becomes paramount for a comprehensive
understanding of DIS at next-to-eikonal accuracy across the full kinematic
regime.

For a complete picture of next-to-eikonal corrections, it would be interesting
to incorporate the effects of the $x^{-}$ dependence of the field, as well as
quark or antiquark exchanges with the background field. Additionally,
incorporating next-to-next-to-eikonal corrections would enhance the
calculation’s precision, as these corrections include contributions from the
transverse gradient of the $A^{-}$ field component as well as the $A^{+}$
field component. Finally, investigating the relationship between target
averages of Wilson lines and TMD distributions, and the interplay between
next-to-leading and next-to-eikonal effects, could offer valuable insights
into non-eikonal effects.

## Acknowledgments

We thank Guillaume Beuf, Florian Cougoulic and Fabio Domínguez for useful
discussions. P.A. and N.A. have received financial support from Xunta de
Galicia (Centro singular de investigación de Galicia accreditation 2019-2022,
ref. ED421G-2019/05), by the European Union ERDF, by the “María de Maeztu”
Units of Excellence program MDM2016- 0692, and by the Spanish Research State
Agency under project PID2020-119632GBI00. This work was performed in the
framework of the European Research Council project ERC-2018-ADG-835105
YoctoLHC and the MSCA RISE 823947 “Heavy ion collisions: collectivity and
precision in saturation physics” (HIEIC) and has received funding from the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No. 824093. P.A. has been supported by Consellería de Cultura,
Educación e Universidade of Xunta de Galicia under the grant ED481B-2022-050.

##  Appendix A Inclusive dijet production in the dense limit

This section briefly reviews the evaluation of the differential cross section
introduced in Section IV within the dense limit. Here, multiple gluon
exchanges between the quark-antiquark pair and the dense nucleus requires
resummation techniques. For a Gaussian (MV) correlation of gluon fields, this
resummation can be performed analytically, resulting in the following
nonlinear expressions, as derived in [59], for the dipole operator

$\displaystyle d({\bf v},{\bf w})=e^{-\Sigma({\bf v}-{\bf
w})},\quad\Sigma({\bf r})=\frac{Q_{s}^{2}}{4}{\bf
r}^{2}\ln\left(\frac{1}{|{\bf r}|\Lambda_{\rm QCD}}+e\right)$ (100)

and the quadrupole operator:

$\displaystyle Q({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf w})$
$\displaystyle=d({\bf w},{\bf v})d({\bf w}^{\prime},{\bf
v}^{\prime})\left[\left(\frac{F({\bf w}^{\prime},{\bf v},{\bf v}^{\prime},{\bf
w})+\sqrt{\gamma}}{2\sqrt{\gamma}}-\frac{F({\bf w}^{\prime},{\bf
v}^{\prime},{\bf v},{\bf
w})}{\sqrt{\gamma}}\right)e^{\frac{N_{c}}{4}\sqrt{\gamma}}\right.$
$\displaystyle-\left.\left(\frac{F({\bf w}^{\prime},{\bf v},{\bf
v}^{\prime},{\bf w})-\sqrt{\gamma}}{2\sqrt{\gamma}}-\frac{F({\bf
w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf
w})}{\sqrt{\gamma}}\right)e^{-\frac{N_{c}}{4}\sqrt{\gamma}}\right]e^{-\frac{N_{c}}{4}F({\bf
w}^{\prime},{\bf v},{\bf v}^{\prime},{\bf w})+\frac{1}{2N_{c}}F({\bf
w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf w})}.$ (101)

This expression involves the following auxiliary functions:

$\displaystyle F({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf w})$
$\displaystyle=\frac{1}{C_{F}}\ln\frac{d({\bf w}^{\prime}-{\bf v})d({\bf
v}^{\prime}-{\bf w})}{d({\bf w}^{\prime}-{\bf w})d({\bf v}^{\prime}-{\bf
v})},$ (102) $\displaystyle\gamma({\bf w}^{\prime},{\bf v}^{\prime},{\bf
v},{\bf w})$ $\displaystyle=F^{2}({\bf w}^{\prime},{\bf v},{\bf
v}^{\prime},{\bf w})+\frac{4}{N_{c}^{2}}F({\bf w}^{\prime},{\bf
v}^{\prime},{\bf v},{\bf w})F({\bf w}^{\prime},{\bf 2},{\bf v},{\bf
v}^{\prime}),$ (103)

where $\gamma\equiv\gamma({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf w})$.

This approach not only reproduces the dilute approximation from Section V but
also captures higher-order field correlators, extending its validity to lower
momentum regimes. However, its application in this study was limited due to
the current inability to resum non-eikonal field correlators.

The differential cross section can be calculated using the aforementioned
results in Section IV. Furthermore, the Fourier harmonics can be efficiently
computed using the following formula [47]:

$\displaystyle D_{n,\lambda}(P_{\perp},\Delta_{\perp})$
$\displaystyle=(-1)^{n}\int_{{\bf B},{\bf r},{\bf r}^{\prime}}J_{n}(|{\bf
B}|\Delta_{\perp})J_{n}(|{\bf r}-{\bf
r}^{\prime}|P_{\perp})\mathcal{C}_{\lambda}({\bf r},{\bf r}^{\prime})\
\Xi({\bf B}+z_{2}{\bf r},{\bf B}-z_{1}{\bf r},-z_{1}{\bf r}^{\prime},z_{2}{\bf
r}^{\prime}),$ (104)

where

$\displaystyle\Xi({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf w})$
$\displaystyle=Q({\bf w}^{\prime},{\bf v}^{\prime},{\bf v},{\bf w})-d({\bf
w}^{\prime},{\bf v}^{\prime})-d({\bf w},{\bf v})+1.$ (105)

###  A.1 The correlation limit

In the correlation limit, where $\Delta_{\perp}\ll P_{\perp}$ (or
equivalently, $|{\bf r}|,|{\bf r}^{\prime}|\ll|{\bf B}|$), we can expand the
dipole and quadrupole operators in powers of the dipole sizes. This expansion
allows to obtain analytical expressions for the differential cross sections in
the eikonal approximation [59]:

$\displaystyle\frac{d\sigma^{\gamma^{*}_{L}+A\to q\bar{q}+X}}{d^{2}{\bf
P}d^{2}{\bf\Delta}d\eta_{1}d\eta_{2}}\Bigg{|}_{\rm Eik.}=8\alpha_{\rm
em}e_{f}^{2}Q^{2}S_{\perp}\delta_{z}z_{1}^{3}z_{2}^{3}\frac{P_{\perp}^{2}}{(P_{\perp}^{2}+\epsilon_{f}^{2})^{4}}\left[\alpha_{s}xG(\Delta_{\perp})+\alpha_{s}xh(\Delta_{\perp})\cos
2\phi\right],$ (106) $\displaystyle\frac{d\sigma^{\gamma^{*}_{T}+A\to
q\bar{q}+X}}{d^{2}{\bf P}d^{2}{\bf\Delta}d\eta_{1}d\eta_{2}}\Bigg{|}_{\rm
Eik.}=\alpha_{\rm
em}e_{f}^{2}S_{\perp}\delta_{z}z_{1}z_{2}(z_{1}^{2}+z_{2}^{2})\frac{P_{\perp}^{4}+\epsilon_{f}^{4}}{(P_{\perp}^{2}+\epsilon_{f}^{2})^{4}}\left[\alpha_{s}xG(\Delta_{\perp})-\frac{2P_{\perp}^{2}\epsilon_{f}^{2}}{P_{\perp}^{4}+\epsilon_{f}^{4}}\alpha_{s}xh(\Delta_{\perp})\cos
2\phi\right].$ (107)

Here, $xG(\Delta_{\perp})$ and $xh(\Delta_{\perp})$ represent the Weizsäcker-
Williams gluon distribution and its linearly polarized partner, respectively.
In the MV model and within the correlation limit, these functions read

$\displaystyle\alpha_{s}xG(\Delta_{\perp})$
$\displaystyle=\frac{N_{c}^{2}-1}{(2\pi)^{3}N_{c}}\int_{0}^{\infty}dBBJ_{0}(\Delta_{\perp}B)\left[1-e^{-\frac{2N_{c}^{2}}{N_{c}^{2}-1}\Sigma(B)}\right]\frac{1}{\Sigma(B)}\left[\frac{d^{2}}{dB^{2}}+\frac{1}{B}\frac{d}{dB}\right]\Sigma(B),$
(108) $\displaystyle\alpha_{s}xh(\Delta_{\perp})$
$\displaystyle=-\frac{N_{c}^{2}-1}{(2\pi)^{3}N_{c}}\int_{0}^{\infty}dBBJ_{2}(\Delta_{\perp}B)\left[1-e^{-\frac{2N_{c}^{2}}{N_{c}^{2}-1}\Sigma(B)}\right]\frac{1}{\Sigma(B)}\left[\frac{d^{2}}{dB^{2}}-\frac{1}{B}\frac{d}{dB}\right]\Sigma(B).$
(109)

The expressions presented here hold validity only in the correlation limit,
$P_{\perp}\gg\Delta_{\perp}$. Conversely, the expressions derived in Eqs. 76
and 77 are valid within the dilute limit, where both $P_{\perp}$ and
$\Delta_{\perp}$ are significantly larger than the saturation scale $Q_{s}$.
Therefore, we anticipate substantial overlap between the correlation and
dilute limits in the region where $P_{\perp}\gg\Delta_{\perp}\gtrsim Q_{s}$.

Indeed, Fig. 15 visually depicts the relative difference between the dilute
and correlation limits as a function of $P_{\perp}$ and $\Delta_{\perp}$. For
better comparison, we employ a smaller value of $\Lambda_{\rm QCD}=0.1$ GeV.
As the figure demonstrates, good agreement between both approaches is observed
within the expected kinematic regime of $P_{\perp}\gg\Delta_{\perp}\gtrsim
Q_{s}$.

Figure 15: Relative difference between the dilute and correlation limits of
the longitudinal (left) and transverse (right) dijet cross sections in the CGC
as a function of $P_{\perp}$ and $\Delta_{\perp}$.

Finally, Eqs. 76 and 77 can be reduced to the following expressions in the
correlation limit ($P_{\perp}\gg\Delta_{\perp}$):

$\displaystyle\frac{d\sigma^{\gamma^{*}_{L}+A\to q\bar{q}+X}}{d^{2}{\bf
P}d^{2}{\bf\Delta}d\eta_{1}d\eta_{2}}\Bigg{|}_{\rm Eik.}$
$\displaystyle=\frac{2C_{L}}{\Delta_{\perp}^{2}}\frac{Q^{2}P_{\perp}^{2}}{(P_{\perp}^{2}+\epsilon_{f}^{2})^{4}}(1+\cos
2\phi),$ (110) $\displaystyle\frac{d\sigma^{\gamma^{*}_{T}+A\to
q\bar{q}+X}}{d^{2}{\bf k}_{1}d\eta_{1}d^{2}{\bf k}_{2}d\eta_{2}}\Bigg{|}_{\rm
Eik.}$
$\displaystyle=\frac{C_{T}}{\Delta_{\perp}^{2}}\frac{P_{\perp}^{4}+\epsilon_{f}^{4}}{(P_{\perp}^{2}+\epsilon_{f}^{2})^{4}}\Bigg{[}1-\frac{2P_{\perp}\epsilon_{f}}{P_{\perp}^{4}+\epsilon_{f}^{4}}\cos
2\phi\Bigg{]}.$ (111)

Using the saturation momentum definition in Eq. 68 and comparing Eqs. 106 and
107 with Eqs. 110 and 111, we can analytically derive the Weizsäcker-Williams
gluon distributions in the limit of $\Delta_{\perp}\gg Q_{s}$:

$\displaystyle
xG(\Delta_{\perp})=xh(\Delta_{\perp})=\frac{2N_{c}C_{F}\tilde{\mu}^{2}S_{\perp}}{\pi^{2}\Delta_{\perp}^{2}}.$
(112)

Interestingly, this result aligns perfectly with the one obtained from
perturbative QCD [59].

## References

  * Gelis  _et al._ [2010] F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, The Color Glass Condensate, Ann. Rev. Nucl. Part. Sci. 60, 463 (2010), arXiv:1002.0333 [hep-ph] .
  * Kovchegov and Levin [2012] Y. V. Kovchegov and E. Levin,  _Quantum chromodynamics at high energy_, Vol. 33 (Cambridge University Press, 2012).
  * Brodsky  _et al._ [1998] S. J. Brodsky, H.-C. Pauli, and S. S. Pinsky, Quantum chromodynamics and other field theories on the light cone, Phys. Rept. 301, 299 (1998), arXiv:hep-ph/9705477 .
  * Hentschinski  _et al._ [2023] M. Hentschinski _et al._ , White Paper on Forward Physics, BFKL, Saturation Physics and Diffraction, Acta Phys. Polon. B 54, 3 (2023), arXiv:2203.08129 [hep-ph] .
  * Ball  _et al._ [2018] R. D. Ball, V. Bertone, M. Bonvini, S. Marzani, J. Rojo, and L. Rottoli, Parton distributions with small-x resummation: evidence for BFKL dynamics in HERA data, Eur. Phys. J. C 78, 321 (2018), arXiv:1710.05935 [hep-ph] .
  * Abdolmaleki  _et al._ [2018] H. Abdolmaleki _et al._ (xFitter Developers’ Team), Impact of low-$x$ resummation on QCD analysis of HERA data, Eur. Phys. J. C 78, 621 (2018), arXiv:1802.00064 [hep-ph] .
  * Arsene  _et al._ [2004] I. Arsene _et al._ (BRAHMS), On the evolution of the nuclear modification factors with rapidity and centrality in d + Au collisions at s(NN)**(1/2) = 200-GeV, Phys. Rev. Lett. 93, 242303 (2004), arXiv:nucl-ex/0403005 .
  * Adler  _et al._ [2005] S. S. Adler _et al._ (PHENIX), Nuclear modification factors for hadrons at forward and backward rapidities in deuteron-gold collisions at s(NN)**(1/2) = 200-GeV, Phys. Rev. Lett. 94, 082302 (2005), arXiv:nucl-ex/0411054 .
  * Adams  _et al._ [2006] J. Adams _et al._ (STAR), Forward neutral pion production in p+p and d+Au collisions at s(NN)**(1/2) = 200-GeV, Phys. Rev. Lett. 97, 152302 (2006), arXiv:nucl-ex/0602011 .
  * Adare  _et al._ [2011] A. Adare _et al._ (PHENIX), Suppression of back-to-back hadron pairs at forward rapidity in $d+$Au Collisions at $\sqrt{s_{NN}}=200$ GeV, Phys. Rev. Lett. 107, 172301 (2011), arXiv:1105.5112 [nucl-ex] .
  * Abdallah  _et al._ [2022] M. S. Abdallah _et al._ (STAR), Evidence for Nonlinear Gluon Effects in QCD and Their Mass Number Dependence at STAR, Phys. Rev. Lett. 129, 092501 (2022), arXiv:2111.10396 [nucl-ex] .
  * ALI [2020] Letter of Intent: A Forward Calorimeter (FoCal) in the ALICE experiment, (2020).
  * Accardi  _et al._ [2016] A. Accardi _et al._ , Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all, Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Abdul Khalek  _et al._ [2022] R. Abdul Khalek _et al._ , Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report, Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419 [physics.ins-det] .
  * Altinoluk  _et al._ [2014] T. Altinoluk, N. Armesto, G. Beuf, M. Martínez, and C. A. Salgado, Next-to-eikonal corrections in the CGC: gluon production and spin asymmetries in pA collisions, JHEP 07, 068, arXiv:1404.2219 [hep-ph] .
  * Altinoluk  _et al._ [2016] T. Altinoluk, N. Armesto, G. Beuf, and A. Moscoso, Next-to-next-to-eikonal corrections in the CGC, JHEP 01, 114, arXiv:1505.01400 [hep-ph] .
  * Casalderrey-Solana and Salgado [2007] J. Casalderrey-Solana and C. A. Salgado, Introductory lectures on jet quenching in heavy ion collisions, Acta Phys. Polon. B 38, 3731 (2007), arXiv:0712.3443 [hep-ph] .
  * Mehtar-Tani  _et al._ [2013] Y. Mehtar-Tani, J. G. Milhano, and K. Tywoniuk, Jet physics in heavy-ion collisions, Int. J. Mod. Phys. A 28, 1340013 (2013), arXiv:1302.2579 [hep-ph] .
  * Blaizot and Mehtar-Tani [2015] J.-P. Blaizot and Y. Mehtar-Tani, Jet Structure in Heavy Ion Collisions, Int. J. Mod. Phys. E 24, 1530012 (2015), arXiv:1503.05958 [hep-ph] .
  * Altinoluk  _et al._ [2021a] T. Altinoluk, G. Beuf, A. Czajka, and A. Tymowska, Quarks at next-to-eikonal accuracy in the CGC: Forward quark-nucleus scattering, Phys. Rev. D 104, 014019 (2021a), arXiv:2012.03886 [hep-ph] .
  * Altinoluk and Beuf [2022] T. Altinoluk and G. Beuf, Quark and scalar propagators at next-to-eikonal accuracy in the CGC through a dynamical background gluon field, Phys. Rev. D 105, 074026 (2022), arXiv:2109.01620 [hep-ph] .
  * Agostini [2023] P. Agostini, Scalar propagator in a background gluon field beyond the eikonal approximation, JHEP 11, 099, arXiv:2307.13573 [hep-ph] .
  * Balitsky and Tarasov [2015] I. Balitsky and A. Tarasov, Rapidity evolution of gluon TMD from low to moderate x, JHEP 10, 017, arXiv:1505.02151 [hep-ph] .
  * Balitsky and Tarasov [2016] I. Balitsky and A. Tarasov, Gluon TMD in particle production from low to moderate x, JHEP 06, 164, arXiv:1603.06548 [hep-ph] .
  * Chirilli [2019] G. A. Chirilli, Sub-eikonal corrections to scattering amplitudes at high energy, JHEP 01, 118, arXiv:1807.11435 [hep-ph] .
  * Chirilli [2021] G. A. Chirilli, High-energy operator product expansion at sub-eikonal level, JHEP 06, 096, arXiv:2101.12744 [hep-ph] .
  * Altinoluk  _et al._ [2023a] T. Altinoluk, N. Armesto, and G. Beuf, Probing quark transverse momentum distributions in the color glass condensate: Quark-gluon dijets in deep inelastic scattering at next-to-eikonal accuracy, Phys. Rev. D 108, 074023 (2023a), arXiv:2303.12691 [hep-ph] .
  * Jalilian-Marian [2017] J. Jalilian-Marian, Elastic scattering of a quark from a color field: longitudinal momentum exchange, Phys. Rev. D 96, 074020 (2017), arXiv:1708.07533 [hep-ph] .
  * Jalilian-Marian [2019] J. Jalilian-Marian, Quark jets scattering from a gluon field: from saturation to high $p_{t}$, Phys. Rev. D 99, 014043 (2019), arXiv:1809.04625 [hep-ph] .
  * Kovchegov  _et al._ [2016] Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, Helicity Evolution at Small-x, JHEP 01, 072, [Erratum: JHEP 10, 148 (2016)], arXiv:1511.06737 [hep-ph] .
  * Kovchegov  _et al._ [2017a] Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, Small-$x$ asymptotics of the quark helicity distribution, Phys. Rev. Lett. 118, 052001 (2017a), arXiv:1610.06188 [hep-ph] .
  * Kovchegov  _et al._ [2017b] Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, Helicity Evolution at Small $x$: Flavor Singlet and Non-Singlet Observables, Phys. Rev. D 95, 014033 (2017b), arXiv:1610.06197 [hep-ph] .
  * Kovchegov  _et al._ [2017c] Y. V. Kovchegov, D. Pitonyak, and M. D. Sievert, Small-$x$ Asymptotics of the Quark Helicity Distribution: Analytic Results, Phys. Lett. B 772, 136 (2017c), arXiv:1703.05809 [hep-ph] .
  * Kovchegov and Sievert [2019] Y. V. Kovchegov and M. D. Sievert, Small-$x$ Helicity Evolution: an Operator Treatment, Phys. Rev. D 99, 054032 (2019), arXiv:1808.09010 [hep-ph] .
  * Cougoulic and Kovchegov [2019] F. Cougoulic and Y. V. Kovchegov, Helicity-dependent generalization of the JIMWLK evolution, Phys. Rev. D 100, 114020 (2019), arXiv:1910.04268 [hep-ph] .
  * Cougoulic and Kovchegov [2020] F. Cougoulic and Y. V. Kovchegov, Helicity-dependent extension of the McLerran–Venugopalan model, Nucl. Phys. A 1004, 122051 (2020), arXiv:2005.14688 [hep-ph] .
  * Kovchegov and Santiago [2021] Y. V. Kovchegov and M. G. Santiago, Quark sivers function at small $x$: spin-dependent odderon and the sub-eikonal evolution, JHEP 11, 200, [Erratum: JHEP 09, 186 (2022)], arXiv:2108.03667 [hep-ph] .
  * Jalilian-Marian [2020] J. Jalilian-Marian, Rapidity loss, spin, and angular asymmetries in the scattering of a quark from the color field of a proton or nucleus, Phys. Rev. D 102, 014008 (2020), arXiv:1912.08878 [hep-ph] .
  * Altinoluk and Dumitru [2016] T. Altinoluk and A. Dumitru, Particle production in high-energy collisions beyond the shockwave limit, Phys. Rev. D 94, 074032 (2016), arXiv:1512.00279 [hep-ph] .
  * Agostini  _et al._ [2019a] P. Agostini, T. Altinoluk, and N. Armesto, Non-eikonal corrections to multi-particle production in the Color Glass Condensate, Eur. Phys. J. C 79, 600 (2019a), arXiv:1902.04483 [hep-ph] .
  * Agostini  _et al._ [2019b] P. Agostini, T. Altinoluk, and N. Armesto, Effect of non-eikonal corrections on azimuthal asymmetries in the Color Glass Condensate, Eur. Phys. J. C 79, 790 (2019b), arXiv:1907.03668 [hep-ph] .
  * Agostini  _et al._ [2022] P. Agostini, T. Altinoluk, N. Armesto, F. Dominguez, and J. G. Milhano, Multiparticle production in proton–nucleus collisions beyond eikonal accuracy, Eur. Phys. J. C 82, 1001 (2022), arXiv:2207.10472 [hep-ph] .
  * Agostini  _et al._ [2023] P. Agostini, T. Altinoluk, and N. Armesto, Finite width effects on the azimuthal asymmetry in proton-nucleus collisions in the Color Glass Condensate, Phys. Lett. B 840, 137892 (2023), arXiv:2212.03633 [hep-ph] .
  * Dumitru  _et al._ [2019] A. Dumitru, V. Skokov, and T. Ullrich, Measuring the Weizsäcker-Williams distribution of linearly polarized gluons at an electron-ion collider through dijet azimuthal asymmetries, Phys. Rev. C 99, 015204 (2019), arXiv:1809.02615 [hep-ph] .
  * Mäntysaari  _et al._ [2020] H. Mäntysaari, N. Mueller, F. Salazar, and B. Schenke, Multigluon Correlations and Evidence of Saturation from Dijet Measurements at an Electron-Ion Collider, Phys. Rev. Lett. 124, 112301 (2020), arXiv:1912.05586 [nucl-th] .
  * Zhao  _et al._ [2021] Y.-Y. Zhao, M.-M. Xu, L.-Z. Chen, D.-H. Zhang, and Y.-F. Wu, Suppressions of dijet azimuthal correlations in the future EIC, Phys. Rev. D 104, 114032 (2021), arXiv:2105.08818 [hep-ph] .
  * Boussarie  _et al._ [2021] R. Boussarie, H. Mäntysaari, F. Salazar, and B. Schenke, The importance of kinematic twists and genuine saturation effects in dijet production at the Electron-Ion Collider, JHEP 09, 178, arXiv:2106.11301 [hep-ph] .
  * Altinoluk  _et al._ [2021b] T. Altinoluk, C. Marquet, and P. Taels, Low-x improved TMD approach to the lepto- and hadroproduction of a heavy-quark pair, JHEP 06, 085, arXiv:2103.14495 [hep-ph] .
  * van Hameren  _et al._ [2021] A. van Hameren, P. Kotko, K. Kutak, S. Sapeta, and E. Żarów, Probing gluon number density with electron-dijet correlations at EIC, Eur. Phys. J. C 81, 741 (2021), arXiv:2106.13964 [hep-ph] .
  * Caucal  _et al._ [2021] P. Caucal, F. Salazar, and R. Venugopalan, Dijet impact factor in DIS at next-to-leading order in the Color Glass Condensate, JHEP 11, 222, arXiv:2108.06347 [hep-ph] .
  * Taels  _et al._ [2022] P. Taels, T. Altinoluk, G. Beuf, and C. Marquet, Dijet photoproduction at low x at next-to-leading order and its back-to-back limit, JHEP 10, 184, arXiv:2204.11650 [hep-ph] .
  * Caucal  _et al._ [2022] P. Caucal, F. Salazar, B. Schenke, and R. Venugopalan, Back-to-back inclusive dijets in DIS at small x: Sudakov suppression and gluon saturation at NLO, JHEP 11, 169, arXiv:2208.13872 [hep-ph] .
  * Caucal  _et al._ [2023] P. Caucal, F. Salazar, B. Schenke, T. Stebel, and R. Venugopalan, Back-to-back inclusive dijets in DIS at small x: gluon Weizsäcker-Williams distribution at NLO, JHEP 08, 062, arXiv:2304.03304 [hep-ph] .
  * Caucal  _et al._ [2024] P. Caucal, F. Salazar, B. Schenke, T. Stebel, and R. Venugopalan, Back-to-Back Inclusive Dijets in Deep Inelastic Scattering at Small x: Complete NLO Results and Predictions, Phys. Rev. Lett. 132, 081902 (2024), arXiv:2308.00022 [hep-ph] .
  * Altinoluk  _et al._ [2023b] T. Altinoluk, G. Beuf, A. Czajka, and A. Tymowska, DIS dijet production at next-to-eikonal accuracy in the CGC, Phys. Rev. D 107, 074016 (2023b), arXiv:2212.10484 [hep-ph] .
  * McLerran and Venugopalan [1994a] L. D. McLerran and R. Venugopalan, Computing quark and gluon distribution functions for very large nuclei, Phys. Rev. D 49, 2233 (1994a), arXiv:hep-ph/9309289 .
  * McLerran and Venugopalan [1994b] L. D. McLerran and R. Venugopalan, Gluon distribution functions for very large nuclei at small transverse momentum, Phys. Rev. D 49, 3352 (1994b), arXiv:hep-ph/9311205 .
  * Boussarie  _et al._ [2023] R. Boussarie _et al._ , TMD Handbook, (2023), arXiv:2304.03302 [hep-ph] .
  * Dominguez  _et al._ [2011] F. Dominguez, C. Marquet, B.-W. Xiao, and F. Yuan, Universality of Unintegrated Gluon Distributions at small x, Phys. Rev. D 83, 105005 (2011), arXiv:1101.0715 [hep-ph] .
  * Kotko  _et al._ [2015] P. Kotko, K. Kutak, C. Marquet, E. Petreska, S. Sapeta, and A. van Hameren, Improved TMD factorization for forward dijet production in dilute-dense hadronic collisions, JHEP 09, 106, arXiv:1503.03421 [hep-ph] .
  * Dumitru  _et al._ [2015] A. Dumitru, T. Lappi, and V. Skokov, Distribution of Linearly Polarized Gluons and Elliptic Azimuthal Anisotropy in Deep Inelastic Scattering Dijet Production at High Energy, Phys. Rev. Lett. 115, 252301 (2015), arXiv:1508.04438 [hep-ph] .
  * Marquet  _et al._ [2016] C. Marquet, E. Petreska, and C. Roiesnel, Transverse-momentum-dependent gluon distributions from JIMWLK evolution, JHEP 10, 065, arXiv:1608.02577 [hep-ph] .
  * Marquet  _et al._ [2018] C. Marquet, C. Roiesnel, and P. Taels, Linearly polarized small-$x$ gluons in forward heavy-quark pair production, Phys. Rev. D 97, 014004 (2018), arXiv:1710.05698 [hep-ph] .
  * Altinoluk  _et al._ [2019] T. Altinoluk, R. Boussarie, and P. Kotko, Interplay of the CGC and TMD frameworks to all orders in kinematic twist, JHEP 05, 156, arXiv:1901.01175 [hep-ph] .
  * Altinoluk and Boussarie [2019] T. Altinoluk and R. Boussarie, Low $x$ physics as an infinite twist (G)TMD framework: unravelling the origins of saturation, JHEP 10, 208, arXiv:1902.07930 [hep-ph] .
  * Catani  _et al._ [1991] S. Catani, M. Ciafaloni, and F. Hautmann, High-energy factorization and small x heavy flavor production, Nucl. Phys. B 366, 135 (1991).
  * Collins and Ellis [1991] J. C. Collins and R. K. Ellis, Heavy quark production in very high-energy hadron collisions, Nucl. Phys. B 360, 3 (1991).
  * Dumitru  _et al._ [2011] A. Dumitru, K. Dusling, F. Gelis, J. Jalilian-Marian, T. Lappi, and R. Venugopalan, The Ridge in proton-proton collisions at the LHC, Phys. Lett. B 697, 21 (2011), arXiv:1009.5295 [hep-ph] .
  * Dumitru  _et al._ [2018] A. Dumitru, G. A. Miller, and R. Venugopalan, Extracting many-body color charge correlators in the proton from exclusive DIS at large Bjorken x, Phys. Rev. D 98, 094004 (2018), arXiv:1808.02501 [hep-ph] .
  * Brodsky and Lepage [1989] S. J. Brodsky and G. P. Lepage, Exclusive Processes in Quantum Chromodynamics, Adv. Ser. Direct. High Energy Phys. 5, 93 (1989).