text:
\section{Introduction}
\label{sect:intro}

One of the main objectives of the future experimental program at the Electron-Ion
Collider (EIC) \cite{Accardi:2012qut,Aschenauer:2017jsk}
 is an in-depth study of the regime of high parton densities in the wavefunction of
a large nucleus accelerated to ultrarelativistic energies. The parton densities, and notably the
gluon distribution, are amplified both by the quantum evolution with increasing energy
and by the presence of a large number $A\gg 1$ of nucleons in the structure of the
nucleus (with $A\simeq 200$ for lead or gold). As a result of this evolution,
the gluon distribution, as measured (at least, indirectly) in deep inelastic scattering (DIS),
rises rapidly with the energy --- roughly like an inverse power of the Bjorken variable $x$,
which is small at high energy:  $x\simeq Q^2/s\ll 1$ when $s\gg Q^2$.
Here, $Q^2$ is the  virtuality of the space-like photon $\gamma^*$
exchanged between the electron and one nucleon from the nuclear target and $s$
is the center-of-mass energy squared for the $\gamma^*$--nucleon scattering.

One expects this rise to be eventually tamed by {\it gluon saturation}, i.e. by non-linear
effects associated with gluon self-interactions, which limit their occupation numbers
to values of order $1/\alpha_s$. The saturation effects should become manifest in the DIS
structure functions at sufficiently small values of $x$ and for relatively low virtualities
$Q^2\lesssim Q_s^2(A,x)$. The saturation momentum $Q_s(A,x)$ is the
typical transverse momentum of the saturated gluons. It increases with the energy,
$Q_s^2\propto x^{-\lambda_s}$ with $\lambda_s\simeq 0.2$,
and also with the nucleon number $A$ ($Q_s^2\propto A^{1/3}$ when $A\gg 1$).
Yet, for the values of $x$ that should be accessible at the EIC and for $A\simeq
200$, this scale $Q_s^2(A,x)$ is still quite  low, in the ballpark of 1 to 2~GeV$^2$. Because of that,
the studies of gluon saturation in inclusive DIS at $Q^2\lesssim Q_s^2$ can be hindered
by non-perturbative contamination: the perturbative quark-antiquark ($q\bar q$)
fluctuation of the virtual photon that one would like to use as a probe of the saturated gluon
distribution can mix with non-perturbative, hadronic (vector mesons), fluctuations. When
$Q^2\simeq 1\div 2$~GeV$^2$, the respective contributions to the inclusive DIS cross-section
may be difficult to separate from each other.

One common strategy to overcome such a difficulty is to consider less inclusive observables,
like multi-particle production, which involve several scales, some of which can be ``semi-hard''
(i.e. of order $Q_s$, or smaller) even for relatively large virtualities $Q^2\gg Q_s^2$
\cite{Gelis:2002nn,Marquet:2009ca,Dominguez:2011wm,Mantysaari:2019hkq}.
The ``golden probe'' which attracted most attention over the last years is the
production of a pair of hadrons (or jets)  in ``dilute-dense'' collisions ($eA$ or $pA$) at forward
rapidities (in the fragmentation region of the dilute projectile)
\cite{Marquet:2007vb,Albacete:2010pg,Dominguez:2011wm,Stasto:2011ru,Metz:2011wb,Lappi:2012nh,Iancu:2013dta,Kotko:2015ura,Dumitru:2015gaa,Altinoluk:2015dpi,Hatta:2016dxp,Marquet:2016cgx,vanHameren:2016ftb,Albacete:2018ruq,Dumitru:2018kuw,Salazar:2019ncp,Mantysaari:2019hkq,Kolbe:2020tlq,Mantysaari:2020lhf}.
Even when the final particles/jets have relatively large transverse momenta ($k_{\perp}^2\sim Q^2\gg Q_s^2$),
the physics of saturation is still visible in the
broadening of the back-to-back peak in their azimuthal angle distribution.
That said, this effect is not exempt from ambiguities either: it can hardly be distinguished
from the broadening due to radiation in the final state (the ``Sudakov effect'') \cite{Mueller:2013wwa}.
In order to reduce the latter, one must again
limit oneself to semi-hard values for $Q^2$ \cite{Zheng:2014vka}.

In this paper, we make the observation that saturation effects in relatively hard $eA$ collisions with
$Q^2\gg Q_s^2$ could be observed at the EIC in an even simpler process, the semi-inclusive
production of a {\it single} hadron (or jet), a.k.a. SIDIS,
provided the measured particle/jet is {\it very forward} --- meaning that it
carries a large fraction $z\simeq 1$ of the longitudinal momentum of the virtual photon.
The precise condition is that
$Q^2(1-z)\lesssim Q_s^2$, which for large virtualities  $Q^2\gg Q_s^2$
forces $z$ to be very close to one: $1-z\ll 1$.
When this happens, saturation effects become important because the
 $q\bar q$ fluctuation which mediates the interactions between the virtual
photon and the nuclear target has a large transverse size, $r\gtrsim 1/Q_s$, hence it typically
scatters off saturated gluons with transverse momenta $p_\perp \sim Q_s$.

As implicit in the above discussion, throughout this paper we use the  ``colour-dipole'' picture for DIS
at small Bjorken $x\ll 1$. This is the natural description in any Lorentz frame  --- such as
the rest frame of the hadronic target, or the center-of-mass frame for the $\gamma^*$--nucleon collision
 --- in which the virtual photon has a large longitudinal momentum, but zero
 transverse momentum ($q_\perp=0$).
In  such a frame and to leading order in pQCD
(but including the high-energy evolution and the non-linear effects associated with gluon saturation), the
DIS process can be described as the interaction between a long-lived $q\bar q$ fluctuation of the virtual
photon (a ``colour dipole'') and the small-$x$ gluons from the target.
The transverse size of this dipole\footnote{This is more precisely the dipole size
in the absence of any scattering; the dipole could be localised on a smaller scale by a
hard scattering with transferred momentum $k_\perp\gg \bar Q$.}
can be estimated as $r\sim 1/\bar Q$, where $\bar Q^2\equiv
z(1-z)Q^2$, with $z$ and $1-z$ the longitudinal fractions of the two quarks.

So long as $\bar Q^2\gg Q_s^2$, we are in the standard ``leading-twist'' regime, where the
dipole-hadron scattering is weak and essentially ``counts'' the number of gluons in the target:
the dipole cross-section is proportional to the gluon distribution.
However, if $z(1-z)$ is sufficiently small, one can have $\bar Q^2\lesssim Q_s^2$ even
at high virtualities $Q^2\gg  Q_s^2$, and then the scattering is strong ---
that is, multiple scattering becomes  important and the dipole amplitude approaches the unitarity
(or ``black disk'') limit, which is the dual description of gluon saturation in a dipole frame.
The impact of multiple scattering/gluon saturation on the final state also depends upon other
factors, like the transverse momentum $k_\perp$ of the produced hadron, the polarisation of
the virtual photon, and the type of process under consideration --- elastic\footnote{By ``elastic'',
we mean processes where also the hadronic target scatters coherently, that is,
where the CGC average is performed as the level of the amplitude. This is in line with the large-$N_c$
approximation that we shall generally use.} or inelastic.

In our subsequent analysis, we shall identify several physical regimes where saturation effects are expected
to be important for SIDIS and will provide leading order calculations for the relevant observables.
Our general framework is the Colour Glass Condensate (CGC)
effective theory \cite{Iancu:2003xm,Gelis:2010nm,Kovchegov:2012mbw},
where the ``leading order approximation''
already includes all-order resummations of multiple scattering (via Wilson lines) and of soft gluon
emissions (via the Balitsky--JIMWLK evolution equations \cite{Balitsky:1995ub,JalilianMarian:1997jx,JalilianMarian:1997gr,Kovner:2000pt,Iancu:2000hn,Iancu:2001ad,Ferreiro:2001qy},
or their BK truncation at large $N_c$  \cite{Balitsky:1995ub,Kovchegov:1999yj}).
The leading-order predictions of this theory for both single particle (quark) and two particle
(quark-antiquark) production\footnote{In applications to phenomenology, such partonic
cross-sections should be convoluted with parton-to-hadron fragmentation functions.}
are well known \cite{Mueller:1999wm,Marquet:2009ca,Dominguez:2011wm}
and will represent the starting point of our analysis.
Up to trivial factors, these cross-sections depend upon $Q^2$ and $z$ only via the {\it effective}
virtuality $\bar Q^2\equiv z(1-z)Q^2$. We shall immediately specialise to the
interesting regime at $\bar Q^2\lesssim Q_s^2$, that we shall study via both analytic and
numerical calculations. As customary in the CGC framework, we shall describe the dipole-nucleus
scattering by using the McLerran-Venugopalan (MV)  model \cite{McLerran:1993ni,McLerran:1994vd}
at low energies and the BK equation for the evolution with increasing energy.

In our analytic studies, we shall mostly focus on the strongly ordered case $\bar Q^2 \ll Q_s^2$;
besides allowing for controlled approximations, this also has the merit to be
conceptually clear: the interesting physical regimes are well separated and the physical interpretation
of the results is transparent. In this regime, the dipole-hadron scattering is typically strong
and implies a transverse momentum transfer of order $Q_s$.
Accordingly, the $k_\perp$-distribution of the produced quark
is found to be a Gaussian with dispersion
$\langle k_\perp^2\rangle =Q_s^2$, which describes transverse momentum broadening via multiple scattering.
This Gaussian is surrounded by power-low distributions at both lower momenta and higher momenta.
At lower momenta $ k_\perp^2 \ll Q_s^2$, we find a $1/k_\perp^2$--distribution which is in fact
the distribution generated by the decay $\gamma^*_{\rm T}\to q\bar q$ of the transverse virtual photon
and which is preserved by the elastic scattering in the black disk limit. At larger momenta
$ k_\perp^2 \gg Q_s^2$, we find a power tail $\propto 1/k_\perp^4$ (up to logarithms), which is
produced via a single hard scattering with one of  the valence quarks in the nucleus.

A particularly suggestive way to exhibit the effects of saturation is by constructing the $R_{pA}$ ratio
for SIDIS, that is, the ratio of the SIDIS cross-sections in $eA$ and respectively $ep$ collisions,
with the latter scaled up by the number of participants. We find that this ratio behaves very
similarly to the corresponding ratio  in d+Au collisions, as measured at
RHIC \cite{Arsene:2004ux,Adams:2006uz} and theoretically studied in several publications
\cite{Kharzeev:2002pc,Baier:2003hr,
Albacete:2003iq,Kharzeev:2003wz,Iancu:2004bx,Blaizot:2004wu}: that is, it shows nuclear suppression
at low momenta, $k_\perp <  Q_s$, and nuclear enhancement at larger momenta, $k_\perp \gtrsim  Q_s$,
with the possible appearance of a (wide) Cronin peak near $k_\perp =  Q_s$.
This similarity is in fact natural: in both problems, the cross-section for
particle production is driven by the ``dipole'' unintegrated gluon distribution in the hadronic
target (a proton or a nucleus).

The results above mentioned strictly hold in the MV model, i.e. at relatively low energies (say, $x\sim 10^{-2}$).
We will numerically study their high-energy evolution using the collinearly improved version of
the BK equation with running coupling recently proposed in \cite{Ducloue:2019ezk,Ducloue:2019jmy}
(see also the related works \cite{Beuf:2014uia,Iancu:2015vea,Iancu:2015joa,Lappi:2015fma,Lappi:2016fmu,Albacete:2015xza,Beuf:2020dxl}).
We shall thus find that the effects of the evolution are quite mild up to moderately large
rapidities (as relevant for the phenomenology at the EIC). The only qualitative changes are the emergence
of an ``anomalous'' dimension in the power-law distribution at high $k_\perp$
and, especially,
a rapid suppression of the Cronin peak, which disappears after just one unit in rapidity --- once again, very similarly to what happens in
d+Au collisions at RHIC \cite{Arsene:2004ux,Adams:2006uz,Albacete:2003iq,Kharzeev:2003wz,Iancu:2004bx}.

Since the identification of very forward jets, or hadrons, with semi-hard $k_\perp\sim Q_s$
may be experimentally challenging, we also provide estimates for the SIDIS cross-sections integrated
over the transverse momentum $k_\perp$ --- that is, the cross-sections for producing a ``quark''
(physically, a hadron, or jet) with a given longitudinal fraction $z$, but any transverse momentum.
Once again, the most interesting regime for studies of saturation lies at $\bar Q^2\ll Q_s^2$. In this regime,
we find that the elastic and inelastic cross-sections are both dominated by the black disk limit and
hence they give identical contributions (at least, at leading order), as expected from the optical theorem.
This is especially interesting since the respective contributions come from integrating $k_\perp$--dependent
physics which is quite different in the two cases: the $\gamma^*_{\rm T}$--wavefunction in the case of the elastic
scattering and, respectively, the $k_\perp$--broadening of the measured quark for the inelastic one.
We also find it interesting to study the weak scattering regime at large effective virtuality $\bar Q^2\gg
Q_s^2$, where saturation is still visible, albeit rather indirectly, via the phenomenon of {\it geometric
scaling}  \cite{Stasto:2000er,Iancu:2002tr,Mueller:2002zm}.

A striking property of the integrated cross-section is its rapid increase with $z$ when
 $z$  is large enough (say, $z\ge 0.9$). Such a rise should be easily
seen experimentally, although it could be more difficult to distinguish between the {\it power-law} increase
expected in the dilute regime at $\bar Q^2\gg Q_s^2$ and, respectively, the {\it logarithmic} increase
in the black-disk regime at  $\bar Q^2\lesssim Q_s^2$. Once again, this distinction can become
sharper by constructing appropriate $R_{pA}$ ratios, that we shall introduce and study
in Sect.~\ref{sect:qllqs}.

Let us conclude this discussion with some remarks on the experimental
challenges to be faced when searching for saturation effects in SIDIS at forward rapidities.
First, these effects are associated with relatively rare events:
the very asymmetric ($z\simeq 1\gg 1-z$) dipole
configurations of interest for us here represent only a small
fraction of the total DIS cross-section when $Q^2\gg Q_s^2$.
Of course, asymmetric configurations {\it per se} are {\it not} that rare:
the {\it aligned jet} configurations, which dominate the DIS structure
functions at small $x$ and very large $Q^2$, are asymmetric too
(when viewed in the dipole frame, they correspond to $q\bar q$ pairs with $z(1-z)\ll 1$).
Yet, the configurations that we consider in this paper are even {\it more} asymmetric:
they must obey $z(1-z)\lesssim Q_s^2/Q^2\ll 1$ (in order to be sensitive to gluon saturation),
unlike the traditional aligned-jet configurations, for which $Q_s^2/Q^2\ll z(1-z) \ll 1$
(corresponding to weak scattering and dilute parton distributions; see the discussion
in Sect.~\ref{sect:sym-al}).

But albeit very asymmetric and rare when $Q^2\gg Q_s^2$, the configurations relevant for saturation
become more and more symmetric, and also typical, when decreasing $Q^2$ towards $Q_s^2$.
So, in practice, one can use SIDIS to study the saturation physics by simultaneously
varying $Q^2$ and $z$, provided the constraint $z(1-z)Q^2\lesssim Q_s^2$ remains satisfied.
We believe that this leaves a sizeable phase-space to be explored by the experiments.

Another tricky issue, that we shall not explicitly address in this work, is the final-state evolution
of the leading parton and, related to that, the possibility to reconstruct its $z$ fraction from
the measured final state (a hadron or a jet). Clearly, this fraction can only be degraded
by the final state radiation. In order to keep a good control on it,
we foresee two possible strategies, depending upon the type of measurement.
If one measures a hadron in the final state, then one should focus on hadrons
with a large longitudinal momentum fraction $z_h\simeq 1$ w.r.t. the virtual photon. Since
the longitudinal fraction of the leading parton is even larger, $z_h< z< 1$, we would already
have some control on the physical regime at work.
To have a more precise estimate for $z$, the data should be unfolded with the parton-to-hadron
fragmentation functions near $z=1$. Measuring jets instead of hadrons may look like a better
option, since a jet captures the quasi-totality of the longitudinal momentum of the leading
parton: $z_{\rm jet}$ is a good proxy for the $z$ of the quark. But this rises the question
about the possibility to reconstruct jets with semi-hard (total) transverse momenta $k_\perp
\sim Q_s$ and at very forward rapidities.

The paper is organised as follows. Sect.~\ref{sect:dijet} presents a brief summary of known results
for particle production in DIS in the dipole picture and in the CGC formalism at leading order.
We start with the dijet cross-section, for which the physical and diagrammatic interpretations are
more transparent, and deduce the SIDIS cross-section by ``integrating out'' the kinematics of the
unmeasured quark.  Sects.~\ref{sect:qllqs} and \ref{sect:intk} contain our main new results.
In Sect.~\ref{sect:qllqs} we consider the double-differential (in $z$ and $k_\perp$) SIDIS
cross-sections in the saturation regime at $\bar Q^2=z(1-z)Q^2< Q_s^2$. We present both
analytic results (valid in the strictly-ordered limit $\bar Q^2\ll Q_s^2$ and for
different ranges in $k_\perp$) and numerical results, which include the effects of the
high-energy evolution. In Sect.~\ref{sect:intk}, we discuss the SIDIS cross-section differential in $z$
but integrated over $k_\perp$. Once again, we combine analytic approximations (in the limiting
cases $\bar Q^2\ll Q_s^2$ and $\bar Q^2\gg Q_s^2$) and numerical results (for the whole kinematics).
In Sect.~\ref{sect:sym-al} we discuss the inclusive DIS cross-section (as obtained by integrating
the SIDIS cross-section over both $z$ and $k_\perp$), with the purpose to expose the signatures
of saturation and to clarify the relative importance of the large-$z$ SIDIS configurations at high
$Q^2\gg Q_s^2$. We notably explain the difference between these configurations and the
better known aligned jets. Sect.~\ref{sec:conc} contains our summary and conclusions.
We have relegated some of the calculations to five appendices, out of which one is devoted
to a brief review of the collinearly-improved BK equation.

\section{DIS in the dipole picture: from dijets to SIDIS}
\label{sect:dijet}

As mentioned in the Introduction, both the physical picture and the mathematical description of DIS
at small Bjorken $x$ are greatly simplified if one describes this process in a {\it dipole frame},
by which one means a Lorentz frame in which the virtual photon ($\gamma^*$) has a large longitudinal
momentum, $q_z^2\gg Q^2$, and zero transverse momentum, $q_\perp=0$.
This choice is not unique: two such frames can be related
to each other by a boost  along the longitudinal direction.
As a matter of facts, most of the subsequent developments
in this paper are boost-invariant, so they remain valid in any dipole frame\footnote{In
applications to the phenomenology, these results must of course be translated to the
laboratory frame at the EIC.}.
That said, both for the purposes of the physical discussion and for the description of
the high energy evolution\footnote{At leading order, the predictions of the BK equation are
boost-invariant, but this symmetry is spoilt by the resummation of the higher-order
corrections amplified by large collinear logarithms. The collinear resummation  that
we shall employ here is adapted to the evolution of the target
 \cite{Ducloue:2019ezk,Ducloue:2019jmy} (see also Appendix~\ref{sec:BK} for a summary).},
 it is convenient to work in a special
frame where both $\gamma^*$ and the target are ultrarelativistic,
but such that most of the total energy is still carried by the target. In this frame, the
high-energy evolution is fully encoded in the gluon distribution of the target,
whereas the virtual photon has just enough energy to fluctuate
into a quark-antiquark ($q\bar q$) pair  --- a colour dipole ---, which then scatters off the target.
As a consequence of this scattering, the quark and the antiquark from the dipole can
lose their initial coherence and evolve into independent jets, or hadrons, in the final state.

Specifically, we chose the virtual photon to be the right mover, with light-cone (LC) 4-momentum
$q^\mu\equiv (q^+, q^-, \bm{q}_\perp)= (q^+, -\frac{Q^2}{2q^+},
\bm{0}_\perp)$, whereas the nucleus is a left-mover with\footnote{We neglect the
proton mass $M$ which is much smaller than all the other scales in the
problem, $M^2\ll Q^2\ll 2P\cdot q$.} 4-momentum $P^\mu=\delta^{\mu-}P^-$ per nucleon.
In the high-energy or small Bjorken $x$ regime, with
\begin{eqnarray}\label{xBj}
x\equiv\,\frac{Q^2}{2P\cdot q}=\,\frac{Q^2}{2P^- q^+}\,\ll\,1\,,
\end{eqnarray}
the coherence time $\Delta x^+\simeq 2q^+/Q^2$ of the virtual photon,
i.e.\ the  typical lifetime of its $q\bar q$
fluctuation, is much larger than the longitudinal extent $\sim A^{1/3}/P^-$ of
the nuclear target. This justifies the use of the eikonal approximation when computing the
dipole-target scattering.

\subsection{Dijet production in DIS at small $x$}

To leading order in perturbative QCD, the most detailed information about the final state
which is still inclusive with respect to the nucleus is contained in the cross-section for producing
an on-shell $q\bar q$ pair. This depends upon the polarisation (transverse or longitudinal) state
of the virtual photon, with the following  results \cite{Gelis:2002nn,Dominguez:2011wm,Mantysaari:2019hkq}:
\begin{align}
\label{2j-sigmaL}
\frac{{\rm d}\sigma^{\gamma^*_{{\rm L}} A \to q\bar q X}}{ {\rm d}^3k_1\, {\rm d}^3k_2} = &\, 8(2\pi)^2
 \frac{N_c \alpha_{\rm em} e_q^2}{q^+}\,\delta(q^+-k_1^+-k_2^+)\,
 z (1-z)
\bar{Q}^2
	\int \frac{{\rm d}^2 \bm{x}_1}{(2\pi)^2}\frac{{\rm d}^2 \bm{x}_1^\prime}{(2\pi)^2}
	\frac{{\rm d}^2 \bm{x}_2}{(2\pi)^2}\frac{{\rm d}^2 \bm{x}_2^\prime}{(2\pi)^2}
	 \nonumber\\*[0.2cm]
     &\times{\rm e}^{-{\rm i} \bm{k}_1 \cdot (\bm{x}_1 - \bm{x}_1^\prime)} \,
	{\rm e}^{-{\rm i} \bm{k}_2 \cdot (\bm{x}_2 - \bm{x}_2^\prime)}\,	{\rm K}_0(\bar{Q}r) {\rm K}_0(\bar{Q}r')\,
	\mathcal{O}^{(4)}_{x_g}(\bm{x}_1, \bm{x}_2, \bm{x}_1^\prime, \bm{x}_2^\prime),
\end{align}
and respectively:
\begin{align}
\label{2j-sigmaT}
\frac{{\rm d}\sigma^{\gamma^*_{{\rm T}} A \to q\bar q X}}{ {\rm d}^3k_1\, {\rm d}^3k_2} = &\, 2(2\pi)^2
 \frac{N_c \alpha_{\rm em} e_q^2}{q^+}\,\delta(q^+-k_1^+-k_2^+)\,
 [z^2+ (1-z)^2]
\bar{Q}^2
	\int \frac{{\rm d}^2 \bm{x}_1}{(2\pi)^2}\frac{{\rm d}^2 \bm{x}_1^\prime}{(2\pi)^2}
	\frac{{\rm d}^2 \bm{x}_2}{(2\pi)^2}\frac{{\rm d}^2 \bm{x}_2^\prime}{(2\pi)^2}
	 \nonumber\\*[0.2cm]
	&\times	{\rm e}^{-{\rm i} \bm{k}_1 \cdot (\bm{x}_1 - \bm{x}_1^\prime)} \,
	{\rm e}^{-{\rm i} \bm{k}_2 \cdot (\bm{x}_2 - \bm{x}_2^\prime)}     \,
	\frac{\bm{r} \!\cdot \!\bm{r}'}{r\, r'}\,{\rm K}_1(\bar{Q}r) {\rm K}_1(\bar{Q}r')\,
     \mathcal{O}^{(4)}_{x_g}(\bm{x}_1, \bm{x}_2, \bm{x}_1^\prime, \bm{x}_1^\prime).
\end{align}
In writing these expressions,
we used longitudinal momentum conservation, $q^+=k_1^++k_2^+$, to write
$k_1^+=zq^+$ and  $k_2^+=(1-z)q^+$, with $0< z<1$.
Also, $\bar{Q}^2 \equiv z(1-z) Q^2$, while $\bm{r}\equiv\bm{x}_1-\bm{x}_2 $
and $\bm{r}'\equiv \bm{x}_1^\prime - \bm{x}_2^\prime$ are the transverse sizes of the dipole in the direct amplitude (DA) and in the complex-conjugate amplitude (CCA), respectively.
The modified Bessel functions  ${\rm K}_0(\bar{Q}r) $ and ${\rm K}_1(\bar{Q}r) $
belong to the light-cone wavefunctions describing the $q\bar q$ fluctuation of the virtual photon
for the two possible polarisations. These functions exponentially vanish
when $\bar{Q}r\gtrsim 1$, showing that the transverse size of the dipole cannot be much larger than
$1/\bar Q$.
Via the uncertainty principle, this in turn implies that the quark and the antiquark are produced by the
decay of the photon with typical transverse momenta $k_{i\perp}^{\rm prod}\sim \bar Q$. Of course, these
{\it production} momenta can be subsequently modified by the scattering,
so their {\it measured} values $k_{1\perp}$ and $k_{2\perp}$ can be very
different from $\bar Q$.

The QCD scattering between the dipole and the gluons from the nuclear target
is encoded in the following $S$-matrix structure,
\begin{eqnarray}\label{O4}
\mathcal{O}^{(4)}_{x_g}(\bm{x}_1, \bm{x}_2, \bm{x}_1^\prime, \bm{x}_2^\prime)
\equiv     1 + S^{(4)}_{x_g}(\bm{x}_1, \bm{x}_2 ; \bm{x}_2^\prime, \bm{x}_1^\prime)-
     S^{(2)}_{x_g}(\bm{x}_1, \bm{x}_2) - S^{(2)}_{x_g}(\bm{x}_2^\prime, \bm{x}_1^\prime),
     \end{eqnarray}
where the two-point (dipole) and four-point (quadrupole) functions are defined as
\begin{equation}\label{dipole}
S^{(2)}_{x_g}(\bm{x}_1, \bm{x}_2)
\,\equiv\,\frac{1}{N_{c}}\left\langle \mathrm{tr} \,\big(
V(\bm{x}_1)V^{\dagger}(\bm{x}_2)\big)\right\rangle_{x_g},
\end{equation}
and, respectively,
\begin{equation}
\label{quadrupole}
S^{(4)}_{x_g}(\bm{x}_1,\bm{x}_2;\bm{x}_2^\prime,\bm{x}_1^\prime)\,\equiv\,\frac{1}{N_{c}}\,\left\langle
\mathrm{tr}\big(V(\bm{x}_1)V^{\dagger}(\bm{x}_2)V(\bm{x}_2^\prime)
V^{\dagger}(\bm{x}_1^\prime)\big)\right\rangle_{x_g}.
\end{equation}
Here, $V(\bm{x})$ and $V^{\dagger}(\bm{x})$
are Wilson lines in the fundamental representation, which describe multiple scattering for
the quark and the anti-quark, in the eikonal approximation; explicitly,
\begin{eqnarray}\label{Wilson}
  V(\bm{x})\,=\,{\rm T}\exp\left\{ ig\int dx^{+}\, t^{a}A^{-}_a(x^{+},\, \bm{x})\right\}.
   \end{eqnarray}
The colour field  $A^{-}_a$ is a random quantity whose correlations (computed
within the CGC effective theory) encode the information about the gluon distribution in the target.
Via the (non-linear) B-JIMWLK evolution, these correlations depend upon the value $x_g$ of the longitudinal momentum fraction of the gluons which participate in the collision. To the leading order approximation at hand, this value is determined by energy-momentum conservation
together with the condition that the produced ``jets'' (the $q\bar q$ pair) be on mass-shell:
the total light-cone energy in the final state, that is, $k_1^-+k_2^-$ with $k_i^-=k^2_{i\perp}/2k_i^+$,
must be equal to the sum of the (negative) LC energy of the virtual photon, $q^-=-Q^2/2q^+$,
 and the (positive) LC energy $p^-=x_gP^-$ transferred by the gluons from the target, via the collision. This
 condition yields
 \begin{eqnarray}\label{xg}
 x_g = \frac{1}{2P\cdot q}\left(\frac{k^2_{1\perp}}{z} + \frac{k^2_{2\perp}}{1-z} + Q^2\right).
 \end{eqnarray}
Note that this effective value $x_g$ can be significantly larger than the Bjorken $x$ variable \eqref{xBj}
if the final transverse momenta $k_{i\perp}$ are considerably larger than $\bar Q$. For what follows,
it is useful to keep in mind that the {\it typical} transverse momentum that can be transferred
by the collision is of the order $Q_s\equiv Q_s(A, x_g)$ --- the nuclear saturation momentum at the longitudinal
scale $x_g$. So, we typically expect  $k_{i\perp}^2\sim {\rm max}\big(\bar Q^2, \,Q_s^2\big)$.
 The total transverse momentum ${\bm \Delta}
\equiv \bm{k}_1 + \bm{k}_2$ of the produced pair is generated via multiple scattering, hence its
magnitude is of order $Q_s$.

\begin{figure}
  \centerline{
    \includegraphics[width=0.85\textwidth]{Figs/DIS-2jets-amplitude}}
  \caption{\small Left: The virtual photon splits before crossing the  nuclear shockwave;
  the originally colourless $q \bar q$ pair acquires colour via the scattering off the gluons
in the target. Right: The virtual photon splits
  after crossing the target, thus giving rise to a colourless  $q \bar q$ pair in the final state.}
  \label{fig:2jets-ampli}
\end{figure}

Notice that the non-linear effects associated with the high gluon density in the target
enter the cross-sections \eqref{2j-sigmaL}--\eqref{2j-sigmaT} at two levels: on one hand, via
the multiple scattering of the dipole projectile, as encoded in the Wilson lines, cf. \eqn{Wilson};
on the other hand, via the non-linear effects in the B-JIMWLK evolution of the
multi-partonic $S$-matrices which appear in Eqs.~\eqref{O4}--\eqref{quadrupole}.

For what follows, it is useful to have a more detailed microscopic understanding of this $S$-matrix
structure. To that aim, let us move one step
backwards and study the amplitudes underlying the cross-sections \eqref{2j-sigmaL}--\eqref{2j-sigmaT}.
Within the time-ordered LC perturbation theory, there are two contributing amplitudes,
illustrated in Fig.~\ref{fig:2jets-ampli}: the virtual photon
can split into the  $q\bar q$ pair either long before, or long after, it crosses the nuclear ``shockwave''.
In both cases, the quark and the antiquark are originally produced by the decay $\gamma^*\to
q\bar q$ in the same colour state. For the process in Fig.~\ref{fig:2jets-ampli}.~left, this colour state
can be rotated by the scattering off the nucleus, separately for the quark and the antiquark;
hence, the 2 fermions will generally end up in different colours states, say $\alpha$ and
$\beta$, with transition amplitudes measured by the matrix elements of the Wilson lines.
For the other process, cf. Fig.~\ref{fig:2jets-ampli}.~right,  one clearly has $\alpha=\beta$.
 This discussion motivates the following expression for the total amplitude \cite{Gelis:2002nn},
\begin{eqnarray}\label{dijet-ampl}
\mathcal{A}_{\alpha\beta}(z, \bm{k}_1, \bm{k}_2; q^+, Q^2)
=
\int_{\bm{x}_1,\bm{x}_2}
{\rm e}^{-{\rm i}(\bm{x}_1\cdot\bm{k}_1+\bm{x}_2\cdot\bm{k}_2)}\,\Psi(r, z; q^+, Q^2)
\, \big(V({\bm{x}_1}) V^\dagger({\bm{x}_2})-1\big)_{\alpha\beta},
\end{eqnarray}
where $\Psi(r, z; q^+, Q^2)$ generically denotes the  $q\bar q$ component of the $\gamma^*$
wavefunction, for either transverse or longitudinal polarisations (see Eqs.~(19)--(20) in
 \cite{Dominguez:2011wm} for explicit expressions). The negative sign in front of the unit matrix in \eqn{dijet-ampl} can be viewed as a consequence of probability conservation. The total amplitude
vanishes in the absence of scattering, as expected: the space-like photon cannot decay into a
pair of {\it on-shell} fermions without additional interactions.
After taking the modulus squared of this amplitude and performing the CGC average
over the target field, one indeed generates the $S$-matrix  structure in Eq.~\eqref{O4}.
This discussion also shows that the two dipole $S$-matrices which enter
 \eqn{O4} with negative signs are interference terms between the two processes shown
 in Fig.~\ref{fig:2jets-ampli}.

 It will be also interesting to consider dijet production via {\it elastic} (or ``coherent'') scattering.
 The ``elastic'' version of   \eqn{dijet-ampl} is
 obtained by demanding the $q\bar q$ pair to remain a colour singlet  in the final state and by
 performing the CGC averaging already at the level of the amplitude (thus enforcing elastic
 scattering also for the target); this yields
  \begin{eqnarray}\label{dijet-elas}
  \hspace*{-.4cm}
\mathcal{A}_{\rm el}(z, \bm{k}_1, \bm{k}_2; q^+, Q^2) =
\int_{\bm{x}_1,\bm{x}_2}
{\rm e}^{-{\rm i}(\bm{x}_1\cdot\bm{k}_1+\bm{x}_2\cdot\bm{k}_2)}\,\Psi(r, z; q^+, Q^2)
\, \frac{1}{N_{c}}\left\langle \mathrm{tr} \,\big(
V(\bm{x}_1)V^{\dagger}(\bm{x}_2)-1\big)\right\rangle_{x_g}.
\end{eqnarray}
Clearly the respective cross-sections take the general form in Eqs.~\eqref{2j-sigmaL}--\eqref{2j-sigmaT},
but with a new $S$-matrix structure: $\mathcal{O}^{(4)}_{x_g}\to \mathcal{O}^{\,\rm el}_{x_g}$, with\footnote{One can also introduce a generalisation of \eqn{Oel} in which the CGC average
is performed only at the level of the cross-section and which describes a part of the inelastic diffraction
\cite{Mantysaari:2019hkq}. However the difference w.r.t. our ``fully elastic'' cross-section vanishes
in the large-$N_c$ limit, which is the limit that we shall use in practice anyway.}:
\begin{eqnarray}\label{Oel}
\mathcal{O}^{\,\rm el}_{x_g}(\bm{x}_1, \bm{x}_2, \bm{x}_1^\prime, \bm{x}_2^\prime)
\equiv     1 + S^{(2)}_{x_g}(\bm{x}_1, \bm{x}_2)  S^{(2)}_{x_g}(\bm{x}_2^\prime, \bm{x}_1^\prime)-
     S^{(2)}_{x_g}(\bm{x}_1, \bm{x}_2) - S^{(2)}_{x_g}(\bm{x}_2^\prime, \bm{x}_1^\prime).
     \end{eqnarray}
Finally, the cross-section for inelastic dijet production is  the difference between
the total cross-section and its elastic component, so the associated $S$-matrix structure reads
\begin{eqnarray}\label{Oinel}
\mathcal{O}^{\,\rm inel}_{x_g}=\mathcal{O}^{(4)}_{x_g}- \mathcal{O}^{\,\rm el}_{x_g}=
S^{(4)}_{x_g}(\bm{x}_1, \bm{x}_2 ; \bm{x}_2^\prime, \bm{x}_1^\prime)-
S^{(2)}_{x_g}(\bm{x}_1, \bm{x}_2)  S^{(2)}_{x_g}(\bm{x}_2^\prime, \bm{x}_1^\prime)\,.\end{eqnarray}

\subsection{Single inclusive jet production in DIS at small $x$}
\label{sect:general}

Given the dijet results in the previous section, it is easy to deduce the corresponding
results for the quantity in which we are primarily interested in this work: the cross-section
for the inclusive production of a single ``jet'' (a quark or an antiquark), or SIDIS.
Assume that we measure the quark, for definiteness. Then we have to integrate out
the kinematics of the unmeasured antiquark in Eqs.~\eqref{2j-sigmaL}--\eqref{2j-sigmaT}.
The integral over  $k_2^+$ simply removes the $\delta$-functions for longitudinal
momentum conservation. That over  $\bm{k}_2$ identifies the transverse coordinates
of the antiquark in the DA and in the CCA: $\bm{x}_2=\bm{x}_2^\prime$. In turn, this introduces an
important simplification in the colour structure: the original quadrupole in
\eqn{quadrupole} reduces to an effective dipole, $S^{(2)}_{x_g}(\bm{x}_1, \bm{x}'_1)$, built with the
transverse coordinates of the measured quark  in the DA and in the CCA, respectively.
By also assuming that the target is homogeneous in the transverse plane (a disk with radius $R_A$), a further 2-dimensional integration can be performed to give a factor $\pi R^2_A$. Assembling everything together,
we find the following expressions for the longitudinal and transverse SIDIS cross sections:
\begin{align}
\label{sigmaL}
	\frac{{\rm d}\sigma^{\gamma^*_{{\rm L}} A \to q X}}{{\rm d} z\, {\rm d}^2\bm{k}} =
	8 R^2_A N_c \alpha_{\rm em} e_q^2 z (1-z) \mathcal{J}_{\rm L},
\end{align}
\begin{align}
\label{sigmaT}
	\frac{{\rm d}\sigma^{\gamma^*_{{\rm T}} A \to q X}}{{\rm d} z\, {\rm d}^2\bm{k}} =
	2 R^2_A N_c \alpha_{\rm em} e_q^2 \big[z^2 +(1-z)^2\big]
	\mathcal{J}_{\rm T}.
\end{align}
We have introduced here the reduced cross-sections,
\begin{align}
\label{JL}
	\mathcal{J}_{\rm L}(\bm{k},\bar Q) \equiv
	\pi \bar{Q}^2
	\int \frac{{\rm d}^2 \bm{r}}{(2\pi)^2}\frac{{\rm d}^2 \bm{r}'}{(2\pi)^2}\,
	{\rm e}^{-{\rm i} \bm{k} \cdot (\bm{r} - \bm{r}')}
	{\rm K}_0(\bar{Q}r) {\rm K}_0(\bar{Q}r')
	[T(\bm{r}) + T(-\bm{r}') - T(\bm{r} \!-\! \bm{r}')] ,
\end{align}
\begin{align}
\label{JT}
	\mathcal{J}_{\rm T}(\bm{k},\bar Q)\equiv
	\pi \bar{Q}^2
	\int \frac{{\rm d}^2 \bm{r}}{(2\pi)^2}\frac{{\rm d}^2 \bm{r}'}{(2\pi)^2}
	{\rm e}^{-{\rm i} \bm{k} \cdot (\bm{r} - \bm{r}')}
	\frac{\bm{r} \!\cdot \!\bm{r}'}{r\, r'}\, {\rm K}_1(\bar{Q}r) {\rm K}_1(\bar{Q}r')
	[T(\bm{r}) + T(-\bm{r}') - T(\bm{r} \!-\! \bm{r}')],
\end{align}
which encode all the non-trivial dynamics. As suggested by their notations, $\mathcal{J}_{\rm L}$
and $\mathcal{J}_{\rm T}$ depend upon the virtuality of the
initial photon $Q^2$ and upon the splitting fraction $z$ for the $\gamma^*\to q\bar q$
decay only via the product $\bar Q^2=z(1-z)Q^2$ (the ``effective virtuality'').

In writing Eqs.~\eqref{JL}--\eqref{JT}, we found it useful
to replace the dipole $S$-matrix by the corresponding scattering amplitude,
$T(\bm{r},x_g)\equiv  1-S^{(2)}_{x_g}(\bm{x}_1, \bm{x}_2)$, with $\bm{r}=\bm{x}_1- \bm{x}_2$. Also,
the dependence of $T$ upon the longitudinal fraction $x_g$ of the gluons from the
target was kept implicit,  to simplify the notations.
The relevant values for $x_g$ will be discussed in a moment. But before that,
it is useful to separate the above cross sections into elastic and inelastic pieces, following
Eqs.~\eqref{Oel}--\eqref{Oinel}. For the transverse photon, one has
\begin{align}
\label{JTel}
	\mathcal{J}_{{\rm T},\rm{el}} =
	\pi \bar{Q}^2
	\int \frac{{\rm d}^2 \bm{r}}{(2\pi)^2}\frac{{\rm d}^2 \bm{r}'}{(2\pi)^2}
	{\rm e}^{-{\rm i} \bm{k} \cdot (\bm{r} - \bm{r}')}
	\frac{\bm{r} \!\cdot \!\bm{r}'}{r\, r'}\, {\rm K}_1(\bar{Q}r) {\rm K}_1(\bar{Q}r') T(\bm{r}) T(-\bm{r}'),
\end{align}
\begin{align}
\label{JTinel}
\hspace*{-0.8cm}
	\mathcal{J}_{{\rm T},\rm{in}}=
	\pi \bar{Q}^2\!
	\int \frac{{\rm d}^2 \bm{r}}{(2\pi)^2}\frac{{\rm d}^2 \bm{r}'}{(2\pi)^2}
	{\rm e}^{-{\rm i} \bm{k} \cdot (\bm{r} - \bm{r}')}
	\frac{\bm{r} \!\cdot \!\bm{r}'}{r\, r'}\, {\rm K}_1(\bar{Q}r) {\rm K}_1(\bar{Q}r') [T(\bm{r}) + T(-\bm{r}') - T(\bm{r} \!-\! \bm{r}')-T(\bm{r}) T(-\bm{r}')],
\end{align}
and similarly for the longitudinal one.

As promised, let us now specify the values for $x_g$ (the longitudinal momentum fraction of the gluons
from the target which are involved in the scattering) that should be used in evaluating the above SIDIS
cross-sections. By inspection of \eqn{xg}, one sees that, in the case of two-particle production,
 $x_g$ depends upon the transverse momenta $\bm{k}_1$ and $\bm{k}_2$ of both final partons.
 So, our above treatment of the integral over $\bm{k}_2$, which ignored this additional dependence,
 may look as just an approximation. This is indeed true in general, but not also for the {\it elastic}
 scattering: one can easily check that, for the case of a homogeneous target, the r.h.s of \eqn{JTel}
 is proportional to the delta-function $\delta^{(2)}(\bm{k}_1+\bm{k}_2)$. (In the elastic scattering, there is
 no net transfer of transverse momentum from the target to the produced partons.)
 Accordingly, the integral over $\bm{k}_2$ simply identifies $\bm{k}_2=-\bm{k}_1$,
 and the version of  \eqn{xg} for SIDIS in the elastic channel reads
  \begin{eqnarray}\label{xgel}
 x_g = \frac{1}{2P\cdot q}\left(\frac{k^2_{\perp}}{z(1-z)} + Q^2\right)\qquad\mbox{for elastic SIDIS},
 \end{eqnarray}
 without any approximation. As for the inelastic channel, it is easy to identify what should be
 the right approximation to \eqn{xg}: the integral over $\bm{k}_2$ is controlled
 by the {\it typical} transverse momenta of the unmeasured parton, which as discussed after
  \eqn{xg} are of the order of the largest among $\bar Q^2$ and $Q_s^2$. So in practice
  we shall use  \begin{eqnarray}\label{xginel}
 x_g = \frac{1}{2P\cdot q}\left(\frac{k^2_{\perp}}{z} +
 \frac{{\rm max}\big(\bar Q^2, \,Q_s^2(A,x_g)\big)}{1-z} +
 Q^2\right)\qquad\mbox{for inelastic SIDIS}.
 \end{eqnarray}

 The integrations over $\bm{r}$ and $\bm{r}'$ in Eqs.~\eqref{JL}--\eqref{JTinel} are controlled by a product
of 3 factors: the modified Bessel functions ${\rm K}_0(\bar{Q}r) $ and ${\rm K}_1(\bar{Q}r) $ describing
the wavefunction of the virtual photon, the amplitude $T(\bm{r})$ for the elastic dipole-nucleus scattering, and the exponential phase for the Fourier transform. Each of these factors introduces its
own characteristic scale: the effective virtuality $\bar Q^2=z(1-z)Q^2$, the target
saturation momentum $Q_s(A,x_g)$, and, respectively, the transverse momentum $k_\perp\equiv
|\bm{k}|$ of the produced quark. As already mentioned, the Bessel functions effectively restrict the
integrations to values $r\lesssim 1/\bar Q$ (larger values are exponentially suppressed). Also,
the Fourier phase is rapidly oscillating for dipole sizes $r \gg 1/k_\perp$. Finally, the dipole
amplitude is small, but rapidly increasing with $r$, for small dipoles,
$T(x_g,\bm{r})\sim r^2Q_s^2\ll 1$ when $r\ll 1/Q_s$ (``colour transparency''),
but this growth is tamed by multiple scattering when $r\gtrsim 1/Q_s$.
For larger dipole sizes, $r\gg  1/Q_s$, the amplitude saturates the unitarity limit $T=1$.
So, clearly, the result of these integrations depends upon the competition between the 3 scales
aforementioned.

In this work, we would like to explore the regime of strong scattering where $T\sim 1$. This is interesting
since multiple scattering it the reflexion of gluon saturation when DIS is viewed in the dipole frame.
(The partonic picture of saturation in terms of gluon occupation numbers becomes manifest only in an infinite momentum frame for the target, such as the Bjorken frame \cite{Mueller:2001fv}.)
Hence, we would like to choose the external scales $k_\perp$, $z$, and $Q^2$ in such a way that the integrations
favour dipole sizes $r,\,r'\gtrsim 1/Q_s$. For this to be true, we clearly need the two following conditions:
\begin{equation}
	\label{cons}
	\bar{Q}^2 \lesssim Q_s^2 \quad \textrm{and} \quad k_{\perp}^2 \lesssim Q_s^2.
\end{equation}
(Otherwise the integrations will be heavily suppressed either due to the exponential fall-off of the Bessel function, or due to the oscillations introduced by the Fourier transform.)
Notice that these conditions may be relaxed a bit, since the weak scattering regime at $r Q_s \ll 1$ is affected
too by saturation, via the phenomenon of geometric scaling. Accordingly, in what follows we shall also
study the high-momentum tails at $k_{\perp}^2 \gg Q_s^2$. Yet, as we shall see, the most interesting
signatures of multiple scattering and gluon saturation come indeed from the regime
where the scattering is strong.

\section{Saturation effects in SIDIS at $z(1-z)Q^2\ll Q_s^2$}
\label{sect:qllqs}

The fact that the longitudinal momentum fraction $z$ of the produced particle can be measured
in SIDIS, at least in principle, opens a very interesting possibility for studies of saturation: one can
select processes such that $Q^2$ is sufficiently high,  $Q^2\gtrsim Q_s^2$, for the perturbative treatment
of the $\gamma^*$ decay to be justified, while at the same time $z$ is sufficiently close to 1
for $\bar Q^2=z(1-z)Q^2$ to be (much) smaller than $Q_s^2$ --- meaning that the dipole scattering
is sensitive to gluon saturation in the nuclear target. In this section we shall analytically
study the strongly ordered case $\bar Q^2 \ll Q_s^2$, via suitable approximations.
The analytical results will be supported by numerical calculations which also cover
more general situations, where the two scales are comparable to each other,
with $\bar Q^2 < Q_s^2$ though.
In this analysis, it will be convenient to distinguish between tthree regimes for  the
transverse momentum $k_{\perp}$ of the produced quark:
{\ttfamily (i)} $k_{\perp}^2 \ll Q_s^2$, {\ttfamily (ii)}
$k_{\perp}^2 \sim Q_s^2$, and  {\ttfamily (iii)}  $k_{\perp}^2 \gg Q_s^2$.

\subsection{$k_{\perp}^2 \ll Q_s^2$:  probing the virtual photon wavefunction}
\label{sect:LCWV}

When both $\bar Q^2$ and $k_{\perp}^2$ are much smaller than $Q_s^2$,
the scattering is predominantly elastic. This will be confirmed by the
subsequent calculations, but can  also be understood via the following argument.
In this regime, the integrals in \eqn{JTel} for $\mathcal{J}_{{\rm T},\rm{el}}$
are dominated by relatively large dipole sizes, which obey $1/Q_s \lesssim
r,\,r'\ll 1/\bar Q$, for which the scattering is as strong as possible:  $T\simeq 1$.
(The respective contributions from smaller dipoles are suppressed by colour transparency:
  $T(r)\!\propto\! r^2$ for $r\ll 1/Q_s$.)  On the contrary, for the inelastic cross-section,
 the would-be dominant,  ``black disk'' ($T=1$), contributions to the
 individual terms in  \eqn{JTinel} mutually cancel in their sum. Accordingly, the net
 contribution to $\mathcal{J}_{{\rm T},\rm{in}}$
 comes from smaller dipoles $r,\,r'\lesssim 1/Q_s$ and thus is parametrically suppressed
 compared to $\mathcal{J}_{{\rm T},\rm{el}}$. Clearly, the presence of the last term
 in the r.h.s. of  \eqn{JTinel},  which is quadratic in $T$ and becomes
 important at strong scattering, was essential for the above argument.

 \subsubsection{Elastic scattering near the unitarity limit}

 So let us focus on the elastic cross-section in this section.
As visible in \eqn{JTel}, this has a factorised structure
--- the cross-section is the modulus squared of the elastic amplitude ---, hence it is
enough to discuss one of the two integrations, say the one over $\bm{r}$. As just mentioned,
the dominant contribution when $\bar{Q}^2, \,k_{\perp}^2 \ll Q_s^2$ can be computed
by replacing $T\to 1$ in  \eqn{JTel}. After this replacement, the elastic
amplitude reduces to the Fourier transform (FT) of the $\gamma^*$ wavefunction
(in particular, it becomes independent of the QCD coupling $\alpha_s$). By also using the
integrals \eqref{int1} and \eqref{int2} from Appendix~\ref{sec:app1}, one finds
\begin{align}
\label{JLlow}
\mathcal{J}_{{\rm L},\rm{el}} (\bm{k}) = \frac{1}{4\pi} \frac{\bar{Q}^2}{(k_{\perp}^2 +\bar{Q}^2)^2}
\quad \textrm{for} \quad \bar{Q}^2, k_{\perp}^2 \ll Q_s^2,
\end{align}
and
\begin{align}
\label{JTlow}
\mathcal{J}_{{\rm T},\rm{el}}(\bm{k})  = \frac{1}{4\pi} \frac{k_{\perp}^2}{(k_{\perp}^2 +\bar{Q}^2)^2}
\quad \textrm{for} \quad \bar{Q}^2, k_{\perp}^2 \ll Q_s^2,
\end{align}
for longitudinal and transverse polarisations, respectively.
By definition, these expressions represent the  transverse momentum distributions of the
quark and the antiquark generated by the decay of the virtual photon.  For a virtual photon
with longitudinal polarisation, this distribution is peaked near $k_\perp\sim \bar Q$ (in particular,
it vanishes in the photo-production limit $\bar Q^2\to 0$), whereas for a transverse photon,
it is logarithmically distributed within the range $ \bar{Q}^2 \ll k_{\perp}^2 \ll Q_s^2$, where
$\mathcal{J}_{{\rm T},\rm{el}} \simeq 1/(4\pi k_{\perp}^2)$.

It is quite
remarkable that a DIS measurement in the regime where the QCD scattering is as strong as possible
($T=1$) carries no information about the target, but merely unveils the structure
of the virtual photon. This is due to the fact that, in this strong scattering regime, the internal
structure of the target is not visible for the dipole projectile: the target looks ``black'',
i.e. totally absorptive, for any $r\gg 1/Q_s$.
In terms of Feynman graphs, the contributions shown in Eqs.~\eqref{JLlow}--\eqref{JTlow}
fully correspond to the second graph in Fig.~\ref{fig:2jets-ampli} --- that is, the decay of $\gamma^*$
occurs {\it after} crossing the nuclear target.
The contribution of the first graph to the elastic scattering is proportional to the dipole $S$-matrix
(recall \eqn{dijet-elas}), which vanishes in the  ``black disk'' limit.
But in the second graph, the dipole does not ``see'' the target at all,
so in that case $S(r)=1$ independently of the value of $r$.

To summarise, in the black disk/deeply saturated regime at $\bar{Q}^2\ll k_{\perp}^2 \ll Q_s^2$,
the cross-section for SIDIS is controlled by the elastic scattering of the $q\bar q$ fluctuation of
a virtual photon with transverse polarisation and can be approximated as
\begin{align}
\label{JTellow}
\mathcal{J}_{{\rm T}} \simeq \mathcal{J}_{{\rm T},{\rm el}} \simeq
\frac{1}{4\pi k_{\perp}^2}
\quad \textrm{for} \quad \bar{Q}^2 \ll k_{\perp}^2 \ll Q_s^2.
\end{align}

Although not manifest in this result for $\mathcal{J}_{{\rm T}}$,
the details of the QCD scattering (in particular, the nature
of the target) are of course important, as they determine its validity range (notably, the upper
limit $Q_s^2\equiv Q_s^2(A,x_g)$ on $k_{\perp}^2$), as well as the corrections beyond this universal
form. We shall shortly evaluate these corrections, both analytically (within the context of the
McLerran-Venugopalan model)
and numerically (from solutions to the BK equation).
To that aim, we also need the relevant value of $x_g$ ---
the longitudinal momentum fraction of the participating gluons from the target  ---, which fixes
the rapidity phase-space for the high-energy evolution. For elastic processes, this is
generally given by \eqn{xgel}, which for the present kinematics reduces to
  \begin{eqnarray}\label{xgel2}
 x_g \simeq \frac{1}{2P\cdot q}\,\frac{k^2_{\perp}}{z(1-z)}\,=x\,
 \frac{k^2_{\perp}}{\bar Q^2}\,\gg\,x\,.
 \end{eqnarray}
As indicated by the last inequality, this value $x_g$ is typically much larger than the Bjorken $x$
variable relevant for inclusive DIS, meaning that the effects of the high energy evolution are less
important than one might naively expect. But before we study this evolution in more detail,
let us focus on the situation at lower energies, say $x_g\sim 10^{-2}$, where the MV model applies.

\subsubsection{Elastic scattering in the MV model}

The McLerran-Venugopalan (MV) model \cite{McLerran:1993ni,McLerran:1994vd} is
a semi-classical approximation for the gluon distribution in a large nucleus, valid at weak coupling and for
moderately high energies --- namely,  so long as one can neglect the effects of the high-energy
evolution. In this model, the nuclear gluon distribution is assumed to be
generated via incoherent, classical, radiation by the valence quarks. This leads to a Gaussian
distribution for the target colour fields $A^{-}_a$ (recall \eqn{Wilson}), for which it is possible to analytically
compute the relevant Wilson-lines correlators. In particular, the dipole $S$-matrix is found as
(see e.g. \cite{Iancu:2002xk,Iancu:2003xm})
\begin{eqnarray}\label{SMV}
S_0(r)=\exp\left\{-\frac{r^2Q_{A}^2}{4}\ln\frac{1}{r^2\Lambda^2}\right\}\,,
\end{eqnarray}
where the quantity in the exponent is the dipole amplitude in the single scattering approximation
(and in the MV model). This quantity $Q_{A}^2$, which has the dimension of a transverse momentum
squared, is proportional to the colour charge density squared of the valence quarks per unit transverse
area. It scales with the nucleon number $A$
and with the QCD coupling $\alpha_s$ like $Q_{A}^2\propto \alpha_s^2 A^{1/3}$. (One factor
of $\alpha_s$ is inherent in the colour charge density, the other one comes from the coupling
to the dipole.) Furthermore,  $\Lambda$ is the QCD confinement scale and the
logarithm has been generated by integrating over Coulomb exchanges with transverse
momenta $q_\perp$ within the range $\Lambda^2 \ll q_\perp^2 \ll 1/r^2$.  \eqn{SMV} is derived
under the assumption that both $Q_A^2$ and  $1/r^2$ are much larger than $\Lambda^2$
and is valid to leading logarithmic accuracy w.r.t. to the Coulomb log.

Our present interest in \eqn{SMV} is twofold. First, as already mentioned,
this will be used as the initial condition for the BK evolution with decreasing $x_g$.
Second, due to its relative simplicity, \eqn{SMV} allows us to deduce relatively accurate, analytic,
approximations for the SIDIS cross-sections of interest for us here. In particular, in this section
we shall construct an approximation to the elastic cross-section $\mathcal{J}_{{\rm T},{\rm el}}$,
which extends the previous result in \eqn{JTellow} to larger transverse momenta,  $k_\perp\sim Q_s$.
The corresponding approximation for the inelastic cross-section $\mathcal{J}_{{\rm T},{\rm in}}$ will
be discussed in the next section.

Before we proceed, let us specify the definition of the saturation momentum $Q_s(A)$ that we shall
use in this context\footnote{Notice that the definition of $Q_s(A)$ to be used in the analytic calculations
based on the MV model is slightly different from the one that we shall generally use, for practical
convenience, in the numerical calculations.}: this is the value $Q_s(A)=2/r$ for which the exponent
in \eqn{SMV} becomes equal to one; this condition yields
\begin{eqnarray}\label{QsMV}
Q_s^2(A)\,\simeq\,Q_{A}^2 \ln\frac{Q_s^2(A)}{\Lambda^2}\,.\end{eqnarray}

We now return to \eqn{JTel} for $\mathcal{J}_{{\rm T},{\rm el}}$ and construct an analytic approximation
valid within the MV model and within the extended range
$\bar{Q}^2 \ll k_{\perp}^2 \lesssim Q_s^2$.  Since the interesting dipole sizes
obey $r \ll 1/\bar Q$, one can keep only the leading order term,
$ {\rm K}_1(\bar{Q}r) \simeq 1/(\bar{Q}r)$, in the small-argument
expansion of the Bessel function; then \eqn{JTel} reduces to\footnote{Incidentally,
the quantity $ \mathcal{W}(\bm{k})$ is formally the same as the so-called Weisz\"acker-Williams
unintegrated gluon distribution as computed in the MV model \cite{Iancu:2003xm}.}
\begin{align}
\label{JTelMV}
	\mathcal{J}_{{\rm T},\rm{el}} \simeq
	\pi \big |\nabla^i_{\bm{k}} \mathcal{W}\big |^2,\qquad\mbox{with}\qquad
	 \mathcal{W}(\bm{k})\equiv
	\int \frac{{\rm d}^2 \bm{r}}{(2\pi)^2}\,
	{\rm e}^{-{\rm i} \bm{k} \cdot \bm{r}}\,
	\frac{T(\bm{r})}{r^2}\,.
\end{align}
Since independent of $\bar Q$, this result is formally the same as the limit of \eqn{JTel}
when $\bar Q\to 0$. As such, it holds independently of the precise form of
the dipole amplitude. But the use of the MV model, i.e.  $T(r)= 1-S_0(r)$ with $S_0(r)$
shown in \eqn{SMV},  enables us to almost exactly compute the integral in \eqn{JTelMV}.
Indeed, one can easily see that the piece of the integral involving $S_0$ is
dominated by $r=2/Q_s$ (the largest value of $r$ which is allowed by the decay of $S_0$).
Hence, in evaluating that piece, one can replace $r\to 2/Q_s$ within the slowly varying
logarithmic factor in the exponent of \eqn{SMV}. That is, one can effectively replace
the dipole $S$-matrix by the Gaussian $S_0(r)\simeq \exp(-{r^2Q_{s}^2}/{4})$, where we have also
used the definition \eqref{QsMV} for  $Q_{s}^2$.
The ensuing integral over $\bm{r}$ can be exactly evaluated, as explained
in Appendix~\ref{sec:elastic}; one thus finds
\begin{align}
\label{JTelGBW}
 \mathcal{J}_{{\rm T},{\rm el}} \simeq
\frac{1}{4\pi k_{\perp}^2}\,{\rm e}^{-\frac{2k_\perp^2}{Q_s^2}}
\qquad \textrm{for} \qquad \bar{Q}^2 \ll k_{\perp}^2 \lesssim Q_s^2.
\end{align}
This result looks indeed as a natural extension of our previous estimate  in \eqn{JTellow}
 to larger transverse momenta  $k_\perp\sim Q_s$.  It has a natural physical interpretation.
 The power $1/ k_{\perp}^2$ is the transverse momentum distribution produced by the decay
 of the transverse virtual photon, as already discussed. The Gaussian $\exp(-{2k_\perp^2}/{Q_s^2})$
 describes the subsequent  broadening of this distribution via {\it elastic}
multiple scattering between the $q\bar q$ dipole and the gluons from the nucleus.
 \eqn{JTelGBW}  shows that the deviations from the pure power law should become
important  for $k_\perp^2 \sim Q_s^2/2$, where the factor 1/2 reflects the elastic nature
of the scattering process.

\subsubsection{Numerical results for $k_{\perp}^2 \lesssim Q_s^2$}

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/JTel_vs_JLel_MV}
    \caption{}\label{fig:lowkT_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/Jel_vs_Jin_MV}
    \caption{}\label{fig:lowkT_b}
  \end{subfigure}
  \caption{\small (a)
  Comparison between the elastic cross-sections in the longitudinal and the transverse sector,
  respectively.
 The (reduced) cross-sections multiplied by $k_{\perp}^2$ are shown as functions of $k_\perp$, for
 $k_\perp\le  Q_{s0}$ and several
 values $\bar{Q}^2 \ll Q_{s0}^2$. Here, $Q_{s0}^2=2$~GeV$^2$ is the saturation momentum
 in the MV model. (b)  Comparison between the elastic and inelastic cross-sections in the transverse sector.
 Elastic scattering dominates for low enough $k_\perp$,
shows a maximum around $k_{\perp}^2 \sim \bar{Q} Q_{s0}$, and then decreases. Inelastic scattering keeps increasing until $k_{\perp} \sim Q_{s0}$; it becomes the dominant contribution already below $Q_{s0}$
   (see also the discussion in Sect.~\ref{sec:medkT}).}
    \label{fig:lowkT}
    \end{figure}

As promised, we shall now compare the  analytic approximations that we previously developed
for the kinematical range at $\bar{Q}^2 \ll Q_s^2$ and $k_{\perp}^2 \lesssim Q_s^2$
with numerical calculations.

We start with the MV model, for which we shall use a slightly modified version of \eqn{SMV},
which is well defined for any value of $r$, including $r> 1/\Lambda$, and thus allows
for unrestricted numerical integrations; this reads (below, $e=2.718...$ is the Euler number)
\begin{eqnarray}\label{SMVnum}
S_0(r)=\exp\left\{-\frac{r^2Q_{A}^2}{2}\ln\left(\frac{1}{r\Lambda}+e\right)\right\}\,.
\end{eqnarray}
Also, in all the numerical simulations, we shall universally define the saturation momentum
as the value $Q_s=2/r$ for which the dipole amplitude is equal to 1/2: $T(r=2/Q_s)=1/2$. This is
slightly different from its previous definition \eqref{QsMV}, but the difference is irrelevant for our purposes.
In practice, we have chosen $Q_{A}^2$ and $\Lambda^2$ such as
$Q_{s0}^2=2$~GeV$^2$ for the MV model $S$-matrix in \eqn{SMVnum}.
(For more clarity, in all the plots and the associated discussions
we use the specific notation $Q_{s0}^2$ for the saturation momentum in the MV model.)

Our numerical results corresponding to the MV model are shown
in the two plots in Fig.~\ref{fig:lowkT}. For convenience, we display results for
the reduced cross-sections multiplied by $k_\perp^2$. That is, we plot
dimensionless quantities like  $k_\perp^2 \mathcal{J}_{{\rm T},{\rm el}}(k_\perp)$
in the saturation region at $k_\perp\le Q_{s0}\simeq 1.4$~GeV and
for 3 different values for the ratio $\bar Q^2/Q_{s0}^2$,
namely, $\bar Q^2/Q_{s0}^2=0.01,\, 0.05,\,0.1$.  In Fig.~\ref{fig:lowkT_a}
we compare the elastic cross-sections in the transverse and in
the longitudinal sector, respectively, whereas
in Fig.~\ref{fig:lowkT_b} we compare
the elastic and the inelastic cross-sections, both in the transverse sector.

The curves in Fig.~\ref{fig:lowkT_a} have the expected shapes in light of the previous
discussion in this paper: $k_\perp^2 \mathcal{J}_{{\rm L},{\rm el}}$ is peaked around
$k_\perp=\bar Q$, whereas  $k_\perp^2 \mathcal{J}_{{\rm T},{\rm el}}$ shows a wide
maximum (``plateau'') at, roughly, $\bar Q< k_\perp < Q_{s0}/\sqrt{2}$. When
increasing $\bar Q$, the height of this maximum is decreasing and its position
gets displaced towards larger values of $k_\perp$. This can be understood
as follows: when $\bar Q\to 0$,  \eqn{JTelGBW} shows that the maximum of
$k_\perp^2 \mathcal{J}_{{\rm T},{\rm el}}$ lies at $k_\perp=0$ and its height is equal to
$1/(4\pi)\simeq 0.08$. The effect of increasing $\bar Q$ can be appreciated by using
the following approximation, which interpolates between  \eqn{JTlow} at low $k_\perp\sim \bar Q\ll Q_{s0}$
and  \eqn{JTelGBW} at larger momenta $\bar Q\ll k_\perp\sim Q_{s0}$:
 \begin{align}
\label{JTlow-inter}
\mathcal{J}_{{\rm T},\rm{el}}(\bm{k})  \simeq \frac{1}{4\pi} \frac{k_{\perp}^2}{(k_{\perp}^2 +\bar{Q}^2)^2}
\,{\rm e}^{-\frac{2k_\perp^2}{Q_s^2}}
\qquad \textrm{for} \quad \bar{Q}^2 \ll Q_s^2 \quad\mbox{and}\quad  k_{\perp}^2\lesssim Q_s^2.
\end{align}
It is not difficult to check that the function $k_\perp^2 \mathcal{J}_{{\rm T},{\rm el}}$
develops a maximum at $k_\perp^2\sim \bar Q Q_{s0}$ and that the height of this maximum
is proportional to $\exp(-2\bar Q/{Q_{s0}})$ --- in qualitative agreement with Fig.~\ref{fig:lowkT_a}.

Fig.~\ref{fig:lowkT_a} also shows that
the longitudinal piece $\mathcal{J}_{{\rm L},{\rm el}}$ dominates at very low momenta
$k_\perp\lesssim \bar Q$, but it is much smaller than the transverse piece
$\mathcal{J}_{{\rm T},{\rm el}}$ at any $k_\perp\gg \bar Q$. Of course, one should be
careful when comparing numerical results for  $\mathcal{J}_{{\rm T}}$ and  $\mathcal{J}_{{\rm L}}$,
as these quantities are weighted by different factors in the respective {\it physical}
cross-sections, cf. Eqs.~\eqref{sigmaL}--\eqref{sigmaT}.  Notably, the longitudinal cross-section
has an additional suppression when $z\to 1$, due to the prefactor $z(1-z)$ in \eqref{sigmaL}.
Altogether, we conclude that the longitudinal contribution is tiny and can
be safely neglected in the physical regime under consideration.

Concerning Fig.~\ref{fig:lowkT_b}, we notice that the elastic cross-section $\mathcal{J}_{{\rm T},{\rm el}}$
is indeed larger than the inelastic one  $\mathcal{J}_{{\rm T},{\rm in}}$ at sufficiently small
$k_\perp\ll Q_s$, as expected, but this hierarchy is changing quite fast when increasing $k_\perp$.
The reasons for this behaviour
should become clear in the next section, where the inelastic contribution will be explicitly
computed.

Furthermore, Fig.~\ref{fig:BKlowkT} shows the results of the high-energy evolution
for the elastic and inelastic (reduced) cross-sections, in the transverse sector alone.
These results have been obtained by solving an improved version of the
BK equation with running coupling (rcBK), which also includes all-order resummations
of large radiative corrections enhanced by double and single collinear logarithms
 \cite{Ducloue:2019ezk,Ducloue:2019jmy}. (For convenience,
this equation is briefly summarised in Appendix~\ref{sec:BK}.)
The evolution ``time'' is the target rapidity\footnote{The use of the target rapidity as opposed
 to the projectile rapidity not only is natural in
the context of DIS, but also has the advantage to automatically ensure the proper
time-ordering (with ``time'' = the LC time $x^+$ of the projectile) of the successive, soft,
gluon emissions, thus avoiding large higher-order corrections --- the {\it anti-collinear}
double logarithms --- which would lead to an unstable evolution \cite{Lappi:2015fma}; see the
discussion in  \cite{Ducloue:2019ezk}.}
 $\eta=\ln(x_0/x_g)$, related to the ``minus''  longitudinal momentum ($p^-=x_g P^-$)
 of the gluons from the target which are involved in the scattering.
Specifically, we use the $S$-matrix in the MV model, \eqn{SMVnum}, as an initial condition at $\eta=0$
and evolve the collinearly-improved rcBK equation up to $\eta=3$.
The value $x_0$ of $x_g$ at which one starts this evolution is irrelevant for our purposes;
in applications to phenomenology, one often takes $x_0=10^{-2}$.
For the subsequent discussion, it is also interesting to know what are the values of the target
saturation momentum for the rapidities to be explored in this paper ($\eta\le 6$). These are
numerically extracted from the solution to the (improved) BK equation, with the results
 shown in Fig.~\ref{fig:Qseta}.

 \begin{figure}
 \centering
  \includegraphics[width=0.5\textwidth,page=1]{Figs/saturation_scales}
 \caption{\small Rapidity dependence of the saturation momentum squared,
 as obtained by solving the collinearly resummed BK equation described in Appendix \ref{sec:BK} and with MV model initial condition with $Q_{s0}^2=2$~GeV$^2$.}\label{fig:Qseta}
\end{figure}

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/Jel_vs_eta}
    \caption{} \label{fig:BKlowkT_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/Jel_vs_Jin_eta}
    \caption{}
    \label{fig:BKlowkT_b}
  \end{subfigure}
  \caption{\small (a) Rapidity evolution of the  transverse elastic cross-section multiplied by $k_{\perp}^2$ for $Q_{s0}^2=2$~GeV$^2$, for 3 different values for $\bar{Q}^2$ and 4 different rapidities ($\eta=0,\,1,\,2,\,3$);
  see the main text for details.
    (b) Rapidity evolution of the transverse elastic and inelastic cross-sections multiplied by $k_{\perp}^2$, for
    fixed $\bar{Q}^2=0.1$~GeV$^2$ and the same 4 values for $\eta$.}
  \label{fig:BKlowkT}
    \end{figure}

In Fig.~\ref{fig:BKlowkT_a} we show the elastic cross-section
$k_\perp^2 \mathcal{J}_{{\rm T},{\rm el}}$ for $\eta\le 3$ and the
same 3 values for the ratio $\bar Q^2/Q_{s0}^2$ as in  Fig.~\ref{fig:lowkT}. The results are plotted
as functions of $k_\perp$, up to $k_\perp= 2$~GeV. (For comparison, we note that
 $Q_s(\eta=3)\simeq 1.77$~GeV, cf.  Fig.~\ref{fig:Qseta}.)
In Fig.~\ref{fig:BKlowkT_b}, we fix $\bar Q^2/Q_{s0}^2=0.05$ and display the functions
$k_\perp^2 \mathcal{J}_{{\rm T},{\rm el}}$ and $k_\perp^2 \mathcal{J}_{{\rm T},{\rm in}}$ for the same
values of $\eta$ and in the same range in $k_\perp$ as in the left plot.  By inspection
of these curves, one sees that the effects of the evolution are rather mild (at least,
for the moderate values of $\eta$ under consideration). This is a consequence of the
aforementioned resummations (the running coupling corrections and the collinear logarithms),
which considerably slow down the evolution as compared to the leading-order BK equation.

Note finally a subtle point that we have glossed over in the above discussion:
strictly speaking, the relevant values for $x_g$ may
depend upon the measured transverse momentum $k_\perp$, as shown in \eqn{xgel2} for the case
of elastic scattering. In practice, we have neglected this dependence, since it would greatly complicate the numerical solutions to the BK equation. However, given the smallness of the evolution effects,
as visible in Fig.~\ref{fig:BKlowkT}, and the logarithmic dependence of $\eta$ upon $x_g$,
 it should be quite clear that such a $k_\perp$--dependent change in the value of $x_g$ cannot
 change our qualitative discussion. As for the inelastic sector, we shall shortly discover (in the next section)
 that this ambiguity is even less important.

\subsection{$k_{\perp}^2 \sim Q_s^2$:  transverse momentum broadening}
\label{sec:medkT}

In the previous section we have seen that, when increasing $k_\perp^2$ above $\bar Q^2$,
the elastic contribution to SIDIS is rapidly decreasing, according to \eqn{JTellow}.
On the other hand, the inelastic contribution is expected to be most important when
$k_\perp \sim Q_s$, since this is the typical transverse momentum transferred by the target to the produced quark via incoherent multiple scattering. This expectation is confirmed by the numerical results in
Figs.~\ref{fig:lowkT} and \ref{fig:BKlowkT}, which show that $k_\perp^2\mathcal{J}_{{\rm T},{\rm in}}$ keeps
increasing when $k_\perp$ increases towards $Q_s$. In this section, we shall construct an analytic
approximation for $\mathcal{J}_{{\rm T},{\rm in}}$ for transverse momenta in the range
$\bar{Q}^2\ll k_{\perp}^2 \sim Q_s^2$ (by which we mean that $k_{\perp}$ can be both smaller,
or larger, than $Q_s$, but not {\it much} larger; the tail of the cross-section at
$k_{\perp} \gg Q_s$ will be studied in Sect.~\ref{sec:highkT}).
This will allow us to physically understand the respective numerical
results in Figs.~\ref{fig:lowkT} and \ref{fig:BKlowkT}, as well as their extension to larger values
$k_{\perp} \gtrsim Q_s$.

\subsubsection{Inelastic scattering off the saturated gluons}

To study the inelastic cross-section \eqref{JTinel} for $k_\perp \sim Q_s$, it is convenient to return
to its representation in terms of dipole $S$-matrices (recall \eqn{Oinel}),
rather than scattering amplitudes. One has
\begin{align}
\label{JTinel2}
\hspace*{-0.1cm}
	\mathcal{J}_{{\rm T},\rm{in}}=
	\pi \bar{Q}^2
	\int \frac{{\rm d}^2 \bm{r}}{(2\pi)^2}\frac{{\rm d}^2 \bm{r}'}{(2\pi)^2}
	{\rm e}^{-{\rm i} \bm{k} \cdot (\bm{r} - \bm{r}')}
	\frac{\bm{r} \!\cdot \!\bm{r}'}{r\, r'}\, {\rm K}_1(\bar{Q}r) {\rm K}_1(\bar{Q}r') [S(\bm{r} \!-\! \bm{r}')-S(\bm{r}) S(-\bm{r}')],
\end{align}
where we recall that $\bar Q^2\ll k_\perp^2\sim Q_s^2$. In this regime, the double
integral is dominated by dipole sizes within the range $1/Q_s\lesssim  r,\,r'\, \ll 1/\bar Q$, which
 (at least,  marginally)  correspond to
strong scattering: the dipole $S$-matrices are significantly smaller than one, albeit
their actual values cannot be neglected, i.e. one cannot work in the unitarity limit $S\to 0$
(or $T\to 1$), since the net result would vanish in that limit. Yet, when $S\ll 1$, the second term
within the square brackets in \eqn{JTinel2}, which is quadratic in $S$, can be neglected next
to the first term there, which is linear in $S$.

 To evaluate the dominant contribution, we first observe that the product ${\rm e}^{-{\rm i} \bm{k} \cdot (\bm{r} - \bm{r}')}
 S(\bm{r} \!-\! \bm{r}')$ restricts the double integration to values $\bm{r}$ and $\bm{r}'$ such that $|\bm{r}- \bm{r}'|
 \lesssim 1/k_\perp \sim 1/Q_s\ll 1/\bar Q$.  This makes it convenient to change
one integration variable, say $\bm{r}'\to \bm{\rho} \equiv \bm{r} - \bm{r}'$. Then, after also using the small-argument
 expansion of the Bessel function,  $ {\rm K}_1(z)\simeq 1/z$ for $z\ll 1$, one sees
 that the integral over $\bm{r}$ at fixed $\bm{\rho}$ shows a logarithmic enhancement over the range
 $\rho \ll r\ll 1/\bar Q$, with $\rho=|\bm{\rho}|$. Indeed, one can write
 \begin{eqnarray}\label{logint}
\frac{\bm{r} \!\cdot \!\bm{r}'}{r\, r'}\, {\rm K}_1(\bar{Q}r) {\rm K}_1(\bar{Q}r')
\,\simeq\,\frac{1}{\bar{Q}^2r^2}\quad\mbox{for}\quad 1/\bar Q \gg r\gg\rho\,,\end{eqnarray}
and the ensuing integral over $\bm{r}$ yields the logarithm $\ln(1/\rho^2 \bar Q^2)$, which is
large when $1/\rho\sim Q_s$. To leading logarithmic accuracy with respect to this transverse
log, one can further write
\begin{align}
\label{JTinlow}
\mathcal{J}_{{\rm T},{\rm in}}
\simeq  &\,
\frac{1}{16\pi^2}
\int {\rm d}^2\bm{\rho}\, {\rm e}^{-{\rm i} \bm{k} \cdot \bm{\rho}} S(\bm{\rho}) \,\ln \frac{1}{\rho^2\bar{Q}^2}
\nonumber\\*[0.2cm]
\simeq  &\,
\frac{1}{16\pi^2}\ln \frac{Q_s^2}{\bar{Q}^2}
\int {\rm d}^2\bm{\rho}\, {\rm e}^{-{\rm i} \bm{k} \cdot \bm{\rho}} S(\bm{\rho})
\quad \textrm{for} \quad \bar{Q}^2 \ll k_{\perp}^2 \sim Q_s^2,
\end{align}
where the second line follows from the fact that the integral over $\rho$
is controlled by values $\rho\sim 1/Q_s$, as we shall shortly check; hence,
to the leading-log accuracy of interest, it is legitimate to replace  $\rho\to 1/Q_s$
within the argument of the logarithm from the first line.

One can verify that the inelastic contribution $\mathcal{J}_{{\rm L},{\rm in}}$ of
the longitudinal photon has no logarithmic enhancement\footnote{Using
$ {\rm K}_0(z)\simeq \ln(1/z)$ for $z\ll 1$, one can check that, in the kinematical
regime at study, the integrals in \eqref{JTinel} are dominated by very large dipoles,
with $r,\,r'\sim 1/\bar Q$, but such that $\rho=|\bm{r} - \bm{r}'|\lesssim 1/k_\perp$. Accordingly,
there is no logarithmic phase-space for the integral over $r$ at fixed $\rho$.}
and hence it is parametrically smaller than the transverse contribution in \eqn{JTbroad}.

The final result in \eqn{JTinlow} has a simple physical interpretation. The
logarithm $\ln({Q_s^2}/{\bar{Q}^2})$
comes from integrating over the {\it initial} transverse momentum distribution
of the quark, as generated by the decay of the virtual photon $\gamma^*_{\rm T}$
and shown in \eqn{JTellow}. The second factor in
the r.h.s. of \eqn{JTinlow} is recognised as the Fourier transform of the dipole $S$-matrix.
It describes the {\it final} transverse momentum distribution of the quark, as generated via
(multiple) scattering off the nuclear target.

The dipole $S$-matrix appearing in \eqn{JTellow} should be evaluated at a rapidity
$\eta=\ln(x_0/x_g)$ with $x_g$ given by \eqn{xginel}; within the present kinematics,
that is,  $k_{\perp}^2 \sim Q_s^2\gg \bar{Q}^2$,  this yields
\begin{eqnarray}\label{xgQs}
 x_g = \frac{1}{2P\cdot q}
\, \frac{Q_s^2}{z(1-z)} \ \Longrightarrow \ \eta
  = \ln \frac{x_0}{x} - \ln  \frac{Q_s^2}{\bar Q^2}\,,
  \end{eqnarray}
with $x$ as defined in \eqn{xBj}.
In the regime of interest here, the rapidity  $\eta$ can be considerably smaller than the
value $\ln(x_0/x)$ that would be normally used for studies of inclusive  DIS.

\subsubsection{Inelastic scattering in the MV model}

To obtain a more explicit expression for $\mathcal{J}_{{\rm T},{\rm in}} $, let us evaluate
 \eqn{JTinlow} within the MV model, that is, by using the dipole $S$-matrix in \eqn{SMV}.
This calculation can be simplified via the same argument as discussed in relation
with \eqn{JTelMV}:
in the interesting range at $k_\perp\lesssim Q_s$, the FT in \eqn{JTinlow} is
dominated by the largest value of $\rho$ that is allowed by the dipole $S$-matrix,
namely $\rho\simeq 2/Q_s$. To the accuracy of interest, one can replace the $S$-matrix
in \eqn{SMV} by the Gaussian $S_0(r)\simeq \exp\{-{r^2Q_{s}^2}/{4}\}$, with $Q_{s}^2$
defined in \eqn{QsMV}. The one easily finds
\begin{align}
\label{JTbroad}
\mathcal{J}_{{\rm T},{\rm in}}
\simeq
\frac{1}{4\pi}\ln \frac{Q_s^2}{\bar{Q}^2}\,\frac{1}{Q_s^2}\,{\rm e}^{-{k_\perp^2}/{Q_s^2}}
\quad \textrm{for} \quad \bar{Q}^2 \ll k_{\perp}^2 \sim Q_s^2.
\end{align}
This Gaussian distribution in $k_\perp$ is the hallmark of a diffusive process:
it describes transverse momentum broadening via soft ($\Lambda^2 \ll q_\perp^2 \lesssim Q_s^2$)
multiple scattering between the produced quark and the gluons in the nucleus.

At this point, we have analytic estimates for the predictions of the MV model
for both the elastic ($\mathcal{J}_{{\rm T},{\rm el}}$)
and the inelastic ($\mathcal{J}_{{\rm T},{\rm in}}$) cross-sections, and for transverse
momenta within the range $\bar Q^2\ll k_\perp^2\sim Q_s^2$. By comparing
the respective results in Eqs.~\eqref{JTelGBW} and \eqref{JTbroad}, one finds that the
elastic scattering dominates over the inelastic one so long as $k_\perp$ is small enough for
$Q_s^2/ k_\perp^2 > \ln ({Q_s^2}/{\bar{Q}^2})$. Using $Q_s^2\gg \bar Q^2$, it is not difficult
to check that the ``critical'' value $k_\perp^2\simeq Q_s^2/\ln ({Q_s^2}/{\bar{Q}^2})$ at which
the elastic and the inelastic contributions are roughly equal to each other lies in between the
respective maxima\footnote{More precisely, these maxima refer to the reduced cross-sections
multiplied by $k_\perp^2$, as in Figs.~\ref{fig:lowkT} and \ref{fig:BKlowkT}.},
namely $k_\perp^2\simeq Q_s\bar Q$ for
the elastic piece and $k_\perp^2\simeq Q_s^2$ for the inelastic one.
This is in agreement with the numerical results shown in
Figs.~\ref{fig:lowkT_b} and \ref{fig:BKlowkT_b}. More results
will be presented in the next subsection.

 \subsubsection{Numerical results for $k_{\perp}^2 \sim Q_s^2$}

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/Jel_vs_Jin_MV_Qs}
    \caption{}\label{fig:el-inel_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/Jel_vs_Jin_eta_Qs}
    \caption{}\label{fig:el-inel_b}
  \end{subfigure}
  \caption{\small (a) The transverse elastic and inelastic cross-sections multiplied by $k_{\perp}^2$,
  plotted as functions of $k_{\perp}^2$ up to $k_\perp^2= 3Q_{s0}^2=6$~GeV$^2$ and
  for 3 different values of $\bar{Q}^2$.
  (b) Rapidity evolution up to $\eta=3$ of the results from the left plot corresponding to
   $\bar{Q}^2=0.1$~GeV$^2$. We note that $Q_s^2(\eta=3)\simeq 3.13$~GeV$^2$.}\label{fig:el-inel}
\end{figure}

The two plots shown in Fig.~\ref{fig:el-inel} extend the previous ones in
Fig.~\ref{fig:lowkT_b} and, respectively, Fig.~\ref{fig:BKlowkT_b}
to larger values for the transverse momenta, namely up to  $k_\perp^2= 3Q_{s}^2$.
Notice also that, as compared to Figs.~\ref{fig:lowkT} and \ref{fig:BKlowkT}, the
rescaled cross-sections  $k_\perp^2 \mathcal{J}_{{\rm T},{\rm in}}$ and
$k_\perp^2 \mathcal{J}_{{\rm T},{\rm el}}$ are now plotted as functions of $k_\perp^2$,
and not of $k_\perp$.
Fig.~\ref{fig:el-inel_a} shows the respective predictions of the MV model
for $Q_{s0}^2=2$~GeV$^2$ and 3 values different for the effective virtuality $\bar Q^2$, namely
$\bar Q^2/Q_{s0}^2=0.01,\, 0.05,\,0.1$.  Fig.~\ref{fig:el-inel_a} illustrates
the effects of the high-energy evolution up to $\eta =3$, for
$k_\perp^2\lesssim 3Q_s^2(\eta=3)$ and  a fixed ratio $\bar Q^2/Q_{s0}^2=0.05$.

We have already discussed the results for the elastic cross-section in some detail,
in relation with Figs.~\ref{fig:lowkT} and \ref{fig:BKlowkT}. So, let us now focus on the
inelastic contributions. The respective MV-model curves in Fig.~\ref{fig:el-inel_a} show the
same qualitative behaviour as the analytic estimate in \eqn{JTbroad}: a rapid
increase at low $k_{\perp}^2\ll Q_s^2$, followed by a wide peak around $Q_s^2$, and finally a
rather slowly decrease at larger values $k_{\perp}^2>Q_s^2$.   Also, when increasing
$\bar Q^2$, the position of the peak remains unchanged --- as it should, since this is located at
$k_\perp^2=Q_s^2$ --- but its height is decreasing, due to the reduction in the logarithmic
factor $\ln({Q_s^2}/{\bar{Q}^2})$ in \eqn{JTbroad}. Notice the difference between
this behaviour and that of the elastic cross-section, for which both the height and the position
of the peak are changing with $\bar{Q}$, as discussed in relation with
Fig.~\ref{fig:lowkT} (and also visible in Fig.~\ref{fig:el-inel_a}).

This general trend is preserved by the BK evolution up to $\eta=3$, as shown in  Fig.~\ref{fig:el-inel_b}.
Due to the increase in the saturation momentum with $\eta$, cf. Fig.~\ref{fig:Qseta}, the
position of the peak is (slightly) moving towards larger values of $k_\perp^2$, for both the
elastic and the inelastic cross-sections. Also, the magnitude of the cross-sections is
growing with $\eta$, albeit quite slowly and only at sufficiently large transverse momenta
 $k_\perp^2\gtrsim Q_s^2(\eta)$. This growth is more important
 for the inelastic channel, where it can be attributed
to the corresponding increase in the logarithm $\ln({Q_s^2(\eta)}/{\bar{Q}^2})$, cf. \eqn{JTinlow}.

\subsection{$k_{\perp}^2 \gg Q_s^2$: hard scattering and BK evolution}
\label{sec:highkT}

When the transverse momentum $k_\perp$ of the produced quark is much larger than
the target saturation momentum,  $k_\perp^2\gg Q_s^2$, the physical picture of SIDIS
is quite different: the (elastic and inelastic) cross-sections are controlled by
relatively small dipoles, with transverse sizes $r,\,r'\ll 1/Q_s$, which scatter only weakly off the
hadronic target: $T(r)\ll 1$. An immediate consequence of that is that the elastic cross-section,
which is quadratic in $T$ (cf. \eqn{JTel}), is even more strongly suppressed relative to the
inelastic one, which in turn is dominated by the first 3 terms, linear in $T$, in \eqn{JTinel}.
That said, gluon saturation may still influence SIDIS even for such high transverse momenta,
via its consequences
(like {\it geometric scaling}  \cite{Stasto:2000er,Iancu:2002tr,Mueller:2002zm})
for the high-energy evolution of the dipole scattering amplitude.

In this section, we shall consider the situation where $k_\perp$ is much larger than both
$Q_s$ and $\bar Q$, without any strong ordering between the two latter scales
(this ordering becomes irrelevant when $k_\perp^2\gg Q_s^2,\,\bar Q^2$).
We shall provide analytic estimates in two limiting situations: for relatively small
energies, where we shall rely on the MV model, and for very high energies,
where we shall use an analytic approximation to the solution to the BK equation.
These estimates will be completed with
numerical solutions to the BK equation, which also apply to intermediate energies.

\subsubsection{Hard scattering in the MV model}

To compute the total (or inelastic) cross-section for $k_\perp^2\gg Q_s^2,\,\bar Q^2$ within
the MV model, one can use the single scattering approximation to the dipole $S$-matrix \eqref{SMV}, that is,
 \begin{eqnarray}\label{MVT0}
 T_0\equiv 1-S_0\simeq \frac{r^2Q_{A}^2}{4}\ln\frac{1}{r^2\Lambda^2}\,,
 \end{eqnarray}
together with the small-argument expansion of
the Bessel function: $ {\rm K}_1(z)\simeq 1/z$ for $z\ll1$.  Under these assumptions,
the dominant contribution to  $\mathcal{J}_{{\rm T}}$
(the only one to be enhanced by a large transverse logarithm)
is that generated by the last term in
\eqn{JT}, i.e. $T(\bm{r} \!-\! \bm{r}')$, which refers to the scattering of the
produced quark. (This term enters \eqn{JT} with a negative sign, but it actually gives a positive
contribution to the FT; see below.) Indeed, for this term, the dipole sizes $r$ and $r'$
can be as large as $1/\bar Q$, since it is only their difference $\bm{\rho}=\bm{r} \!-\! \bm{r}'$
which is constrained by the Fourier phase to small values $\rho\lesssim 1/k_\perp$.
Changing integration variables from $\bm{r},\,\bm{r}'$ to  $\bm{r},\,\bm{\rho}$, the integral
over $r$ is logarithmic within the interval $\rho\ll r\ll 1/\bar Q$ (recall \eqn{logint}),
and one finds
\begin{align}
\label{JTsingle}
\mathcal{J}_{{\rm T}}
& \simeq
\frac{1}{16\pi^2}
\int {\rm d}^2\bm{\rho}\, {\rm e}^{-{\rm i} \bm{k} \cdot \bm{\rho}} \,\ln \frac{1}{\rho^2\bar{Q}^2}
\left(- \frac{\rho^2Q_{A}^2}{4}\ln\frac{1}{\rho^2\Lambda^2}
\right)\nonumber\\*[.2cm]
&\simeq \,\frac{1}{4\pi}\,\frac{Q_{A}^2}{k_\perp^4}\left(\ln \frac{k_\perp^2}{\Lambda^2}
+\ln \frac{k_\perp^2}{\bar{Q}^2}\right)\,
\quad \textrm{for} \quad \bar{Q}^2 \ll Q_s^2 \ll k_{\perp}^2.
\end{align}
The expression in the first line is recognised as the single-scattering
approximation to the first line of  \eqn{JTinlow}. The obtention of the final result
in the second line
is a bit subtle and requires the more elaborated calculations presented in Appendix~\ref{sec:mv}.
Here, we shall rather focus on its physical origin.

Notice first that, for the kinematics
at hand, the integral over $\rho$ is controlled by the Fourier phase, which selects
$\rho\sim 1/k_\perp$. One may naively think that the dominant contribution can be
obtained by replacing $\rho^2\to 1/k_\perp^2$ within both logarithms inside the integrand:
such a contribution would be {\it quadratic} in the transverse logarithms and hence
parametrically larger than the {\it linear} combination shown
in the second line of  \eqn{JTinlow}.  However, if one does so, then the ensuing double
logarithm gets multiplied by an integral over $\bm{\rho}$ which vanishes
when $k_\perp > 0$ (since its result is proportional to $\nabla^2_{\bm{k}}\delta^{(2)}(\bm{k})$).
So, it is natural to find a result which is only {\it linear} in the transverse logs, and also
{\it symmetric} under the exchange $\bar Q^2\leftrightarrow \Lambda^2$, like the integrand
in the first line.

Consider now the physical interpretation of the two terms in the final result. We start
with the second one, whose origin is quite transparent.
As in the case of \eqn{JTinlow}, the $\ln({k_\perp^2}/{\bar{Q}^2})$
comes from integrating over the initial $k_\perp$--distribution
of the quark produced by the decay of $\gamma_{\rm T}^*$,
 \eqn{JTellow}. As for the associated power tail $\propto 1/k_\perp^4$, this
is generated by a hard single scattering between the quark and a nucleon from the
target, with a transferred momentum $q_\perp\simeq k_\perp$.

Concerning the second term,
it is quite clear that the logarithm $\ln({k_\perp^2}/{\Lambda^2})$ is directly inherited from
the dipole amplitude \eqref{MVT0}: this is a {\it Coulomb} logarithm, generated via a relatively
soft exchange, within the range $\Lambda^2\ll q_\perp^2\ll k_\perp^2$. Accordingly, the power-like
tail $\propto 1/k_\perp^4$ for this term should somehow be associated with the decay of the virtual
photon. This may look surprising since, as discussed in Sect.~\ref{sect:LCWV},
the $\gamma_{\rm T}^*$ decay is expected to produce a much softer tail $\propto 1/k_\perp^2$, cf. \eqn{JTellow}.
Yet, as shown by the calculation in Appendix~\ref{sec:mv}, a $1/k_\perp^4$
tail is in fact generated by the {\it difference}
$\frac{1}{(\bm{k}\!-\!\bm{q})^2}-\frac{1}{k_\perp^2}$ between the two such distributions,
 with and without soft rescattering in the final state.

 In principle, \eqn{JTsingle} should be used in the case where
 $\bar Q^2\gg \Lambda^2$, to justify the use of a perturbative (partonic) description
for the $\gamma^*$ wavefunction; in such a case, the first logarithm in the final result would dominate
over the second one. In practice though, the difference between $\bar Q^2$ and  $\Lambda^2$ may
not be that important, which explains why we have kept both logarithms in  \eqn{JTsingle}.

For completeness, let us also observe that the elastic cross-section $\mathcal{J}_{{\rm T}, {\rm el}}$
in the MV model and in the single scattering approximation can be easily estimated from
Eqs.~\eqref{JTel} and \eqref{MVT0}. One thus finds $\mathcal{J}_{{\rm T}, {\rm el}}\sim Q_{A}^4/{k_\perp^6}$
(see Appendix~\ref{sec:elastic} for an explicit calculation and \eqn{Wlim1} for the final result).
This is suppressed  by a power of $ Q_{A}^2/k_\perp^2\ll 1$
compared to the inelastic piece (a ``higher twist effect''), as expected.

\subsubsection{Hard scattering at very high energies}

We now turn to the high-energy limit, where the asymptotic solution to the BK equation in the regime
$r Q_s \ll 1$ is known as \cite{Mueller:2002zm}
\begin{align}
	\label{tbkdil}
	T(r) \simeq (r^2 Q_s^2)^\gamma \ln
	\frac{1}{r^2 Q_s^2}
	\quad \textrm{for} \quad r Q_s \ll 1,
\end{align}
where the ``anomalous dimension''  $\gamma \simeq 0.63$ is significantly smaller than 1.  It is here understood
that $Q_s\equiv Q_s(A,\eta)$ depends upon the rapidity $\eta=\ln(x_0/x_g)$, assumed to be large:
$\bar{\alpha}_s\eta \gg 1$. \eqn{tbkdil} shows {\it geometric scaling} \cite{Stasto:2000er,Iancu:2002tr,Mueller:2002zm}:
it depends upon the 2 independent kinematical
variables $r$ and $\eta$ (and also upon the nucleon number $A$) only via the product $r^2 Q_s^2(A,\eta)$.
This property is an imprint of the physics of saturation on the behaviour of the amplitude in the weak
scattering regime at $rQ_s\ll 1$.

Clearly, the asymptotic behaviour at high energies is not precisely the regime to be explored via
DIS at EIC, where only moderate values $\eta\lesssim 3$ should be accessible. Here, we shall
rather use  \eqn{tbkdil} as a convenient {\it Ansatz} which interpolates
(when decreasing $\gamma$ from 1 to 0.63) between the MV model at low
energy, cf.  \eqn{MVT0}, and the prediction of the BK equation at very high energy.

Given this formal similarity between the BK amplitude in \eqn{tbkdil} and the MV amplitude in \eqn{MVT0},
it is clear that the argument previously used in relation with \eqn{JTsingle}  can be taken
over to the present situation.  The proper generalisation of the first line in  \eqn{JTsingle}  reads
\begin{align}
\label{JTasymp}
\mathcal{J}_{{\rm T}}
\simeq
\frac{1}{16\pi^2}
\int {\rm d}^2\bm{\rho}\, {\rm e}^{-{\rm i} \bm{k} \cdot \bm{\rho}} \,\ln \frac{1}{\rho^2\bar{Q}^2}
\left(- (\rho^2Q_s^2)^\gamma\,
\ln\frac{1}{\rho^2Q_s^2}
\right).
\end{align}
However, since the power $\gamma$ is strictly smaller than one, the mathematical
treatment of \eqn{JTasymp} is considerably simpler. Namely, it is now legitimate to replace
$\rho^2\to 1/k_\perp^2$ within the arguments of both logarithms inside the integrand,
since the large contribution thus obtained
has a non-zero coefficient (given by the integral \eqref{int4} in Appendix~\ref{sec:int}).
One thus finds
\begin{align}
\label{JTbkhigh}
\mathcal{J}_{{\rm T}} \simeq
\mathcal{J}_{{\rm T},{\rm in}} \simeq
\frac{4^{\gamma} \gamma  \Gamma(1+\gamma)}{\Gamma(1-\gamma)}
\frac{1} {8\pi k_{\perp}^2}
\left(\frac{Q_s^2}{k_{\perp}^2}\right)^\gamma
\ln \frac{k_{\perp}^2}{\bar{Q}^2}
\ln \frac{k_{\perp}^2}{Q_s^2}
\quad \textrm{for} \quad \bar{Q}^2 \ll Q_s^2 \ll k_{\perp}^2.
\end{align}
The longitudinal cross section $\mathcal{J}_{{\rm L},{\rm in}}$ can be similarly estimated: one finds the same power-like behaviour as in Eq.~\eqref{JTbkhigh}, but with a single logarithmic enhancement.
Hence,  Eq.~\eqref{JTbkhigh} provides the dominant contribution ---
in the sense of a leading logarithmic accuracy ---  to SIDIS in this high-$k_\perp$ regime,
Comparing this result to the low energy behaviour shown in \eqn{JTsingle}, we observe that the
the high-energy evolution in presence of saturation produces a softer tail at large transverse momenta.
Also,  Eq.~\eqref{JTbkhigh} shows no special scaling, like geometrical scaling, due to the presence
of {\it two} external scales, $\bar Q^2$ and $Q_s^2$, in the final result
(rather than only one such a scale in the dipole amplitude \eqref{tbkdil}). We shall
return to this point in the numerical analysis to follow.

For completeness, let us also show the corresponding result for the
elastic cross section, as obtained via similar manipulations starting with Eq.~\eqref{JTel}:
\begin{align}
\label{JTelbkhigh}
\mathcal{J}_{{\rm T},{\rm el}} \simeq
\frac{4^{2\gamma} \Gamma^2(1+\gamma)}{\Gamma^2(1-\gamma)}
\frac{1} {4\pi k_{\perp}^2}
\left(\frac{Q_s^2}{k_{\perp}^2}\right)^{2\gamma}
\ln^2 \frac{k_{\perp}^2}{Q_s^2}
\quad \textrm{for} \quad \bar{Q}^2 \ll Q_s^2 \ll k_{\perp}^2\,.
\end{align}
As anticipated, this is strongly suppressed as  compared to the inelastic part in Eq.~\eqref{JTbkhigh},
due to the larger power of ${Q_s^2}/{k_{\perp}^2}$: $2\gamma$ instead of $\gamma$.

One interesting conclusion which emerges from the calculations in Sects.~\ref{sec:medkT} and
\ref{sec:highkT} is that, for all transverse momenta $k_\perp\gtrsim
Q_s$ and for any energy (from $\eta\simeq 0$, where the MV model applies,
up the very large $\eta\gg 1/\bar{\alpha}_s$, where one can use the asymptotic solution to the BK equation),
the dominant contribution to SIDIS at low virtuality ($\bar Q^2\ll Q_s^2$) is generated
by the scattering of the measured quark {\it alone} (as opposed to the scattering of the quark-antiquark
dipole).  Of course, the scattering of the whole dipole does matter {\it in general} --- this is represented
by the terms $T(\bm{r})$ and $T(\bm{r}')$ in \eqn{JT} ---, but its contribution turns out to be (parametrically)
less important than that from the measured quark alone --- the piece involving $T(\bm{r} \!-\! \bm{r}')$ in \eqn{JT}.

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/high_kt_tail}
    \caption{} \label{fig:highkt_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/high_kt_tail_log}
    \caption{}\label{fig:highkt_b}
  \end{subfigure}
  \caption{\small Approximate geometric scaling behaviour of the transverse inelastic cross section multiplied by $k_\perp^2$ at (a) low and intermediate transverse momenta and (b) high transverse momenta. With increasing $\eta$ the high-$k_{\perp}$ tail becomes slightly softer.}\label{fig:highkt}
\end{figure}

 \subsubsection{Numerical results for $k_{\perp}^2 \gtrsim Q_s^2$}

In Fig.~\ref{fig:highkt} we show numerical results for the quantity $k_{\perp}^2\mathcal{J}_{{\rm T},{\rm in}}$
at different rapidities $\eta\le 6$ in terms of the dimensionless variable $k_{\perp}^2/Q_s^2(\eta)$, up
to relatively large values $k_{\perp}^2\gg Q_s^2(\eta)$. In the left plot, we use a linear scale and
show all the transverse momenta up to $k_{\perp}^2/Q_s^2(\eta) = 6$, in order to exhibit
the global shape of the
function $k_{\perp}^2\mathcal{J}_{{\rm T},{\rm in}}$. In the right plot, we use a logarithmic
scale and show only the high momentum regime at $1 < k_{\perp}^2/Q_s^2(\eta) < 40$.
Note that after multiplication by $k_{\perp}^2$, the inelastic cross-section $\mathcal{J}_{{\rm T},{\rm in}}$
is very close to a {\it scaling function} (see e.g. Eqs.~\eqref{JTbroad}, \eqref{JTsingle} and
\eqref{JTbkhigh}) --- that is, a function of the ratio  $k_{\perp}^2/Q_s^2(\eta)$. This scaling behaviour
is however not exact (not even in our analytic approximations): it is broken by logarithmic factors.
For instance, \eqn{JTinlow}, which we recall is valid when $k_\perp\sim Q_s$,
 involves an overall factor $\ln ({Q_s^2(\eta)}/{\bar{Q}^2})$, which increases with $\eta$,
 due to the respective behaviour of the saturation momentum. Similarly, at larger momenta
 $k_\perp^2\gg Q_s^2$, one can rely on \eqn{JTbkhigh} to deduce that the function
$k_{\perp}^2\mathcal{J}_{{\rm T},{\rm in}}$ can be written as the sum of a scaling
function and a scaling-violating term proportional to $\ln ({Q_s^2(\eta)}/{\bar{Q}^2})$.
(To that aim, one should rewrite one of the logarithms in \eqn{JTbkhigh} as
$\ln ({k_{\perp}^2}/{\bar{Q}^2}) = \ln ({k_{\perp}^2}/{{Q}_s^2}) + \ln({Q_s^2}/{\bar{Q}^2})$.)
These analytic arguments suggest that, when plotted in terms of the scaling
variable $k_{\perp}^2/Q_s^2(\eta)$, the function
$k_{\perp}^2\mathcal{J}_{{\rm T},{\rm in}}$ should still increase with $\eta$,
albeit only slightly. This is in agreement with the numerical findings in Fig.~\ref{fig:highkt_a}.

Another interesting feature of the numerical results, better seen in Fig.~\ref{fig:highkt_b}, is the
emergence of an anomalous dimension $\gamma < 1$ as a consequence of the high-energy
evolution: indeed, the high-$k_\perp$ tails corresponding to $\eta\gtrsim 4$ are clearly seen
to be softer (i.e. to show a slower decrease with increasing $k_{\perp}^2$) then the MV model
curve, for which $\gamma=1$. That said, the rapidities $\eta\le 6$ considered here are still too small
to feature the asymptotic expectation $\gamma\simeq 0.63$. Besides, the numerical extraction
of this power $\gamma$ is further hindered by the additional logarithmic dependences, as visible
in \eqn{JTbkhigh}.

\subsection{Cronin peak}
\label{sec:cronin}

 As an application of the previous developments in this section, let us discuss
 an interesting physical consequence of the multiple
 scattering and the associated transverse momentum broadening, that might be measured at the
 EIC. Namely, let us compute the nuclear modification factor $R_{pA}$ for the SIDIS process,
 that is, the ratio between particle production (with a given $z$ and $k_\perp$) in $eA$ and
 respectively $ep$ collisions, with the latter scaled up by the number of participants.
 In the approximations at hand, this is obtained as
 \begin{eqnarray}\label{RpA}
 R_{pA}(z,\bm{k})\,=\,\frac{\frac{{\rm d}\sigma^{\gamma^*\!A \to q X}}{{\rm d} z\, {\rm d}^2\bm{k}}}
 {A \frac{{\rm d}\sigma^{\gamma^*p \to q X}}{{\rm d} z\, {\rm d}^2\bm{k}}}\,
 \simeq\,\frac{\mathcal{J}_{\rm T}^A(z,\bm{k})}{A^{1/3}\mathcal{J}_{\rm T}^p(z,\bm{k})}\,,
 \end{eqnarray}
 where in the second equality we have restricted ourselves to the transverse sector, since this
 dominates in the interesting kinematical regime. Namely we are interested
 in transverse momenta of the order of the saturation momentum in the nucleus: $k_\perp\sim Q_s(A)$.

 For an analytic estimate of \eqn{RpA} at relatively low energies (say, $x_g\sim 0.01$), we can use
 the MV model for both types of targets, with $Q_A^2= A^{1/3}Q_p^2$. Specifically, for
  $k_\perp\sim Q_s(A)$, the saturation effects are important for the nuclear target, for which one
  should use the full MV result in \eqn{JTbroad}, but they are not important for the proton target,
for which one should rather use the high-$k_\perp$ approximation in \eqn{JTsingle}.
To make the argument clearer, let us focus on the physically interesting regime at $\bar Q^2\gg
\Lambda^2$; in that case, one can neglect the second logarithm in the final result in \eqn{JTsingle},
which allows us to write
 \begin{eqnarray}\label{RpAMV}
 R_{pA}(k_\perp)\,\simeq\,\frac{\ln ({Q_s^2(A)}/{\bar{Q}^2})}{\ln ({k_\perp^2}/\Lambda^2)}\,
 \frac{ k_\perp^4 \,{\rm e}^{-{k_\perp^2}/{Q_s^2(A)}}}{A^{1/3} Q_{p}^2\, Q_s^2(A)}\
  \underset{k_\perp\sim Q_s(A)}{\longrightarrow}\ \ln\frac{Q_s^2(A)}{\bar Q^2}  \,>\,1,
 \end{eqnarray}
 where the last estimate holds for $k_\perp\sim Q_s(A)$. (We have also used $Q_A^2= A^{1/3}Q_p^2$
 together with \eqn{QsMV} for the nuclear saturation momentum.)
 The first equality in \eqn{RpAMV} also shows that $R_{pA}(k_\perp)$ increases with  $k_\perp$
 so long $k_\perp< Q_s(A)$, but it decreases at larger momenta $k_\perp\gg Q_s(A)$.
 Hence, \eqn{RpAMV} shows a peak at $k_\perp\simeq Q_s(A)$ and the height of this peak
 is expected to be larger than one in the limit where $Q_s^2(A)\gg  \bar Q^2$. This is similar
to the Cronin peak observed in hadron production at midrapidities in d+Au collisions at RHIC
\cite{Arsene:2004ux,Adams:2006uz}. In that context too, this enhancement
can be explained by gluon saturation in the nuclear wavefunction and its consequences in terms
of multiple scattering \cite{Kharzeev:2002pc,Baier:2003hr,
Albacete:2003iq,Kharzeev:2003wz,Iancu:2004bx,Blaizot:2004wu}.

The existence of such a peak at the level of the MV model is numerically studied in
Fig.~\ref{fig:Cronin_a}. This figure shows the $R_{pA}$ ratio in \eqn{RpA} as a function of $k_\perp^2$
(up to $5Q_s^2=10$~GeV$^2$) and for 3 values,  $\bar Q^2/Q_s^2=0.01,\, 0.05,\,0.1$,
for the ratio between the relevant scales. A mild peak, centered at $k_\perp^2$ slightly above
 $Q_s^2$, is indeed seen in the case $\bar Q^2/Q_s^2=0.01$, but not also in the 2 other cases.
A perhaps more robust prediction, valid for all the 3 cases under investigation, is that the ratio
$R_{pA}(k_\perp)$ increases quite fast at low momenta until it reaches a value of order 1 for
$k_\perp\sim Q_s$. For larger $k_\perp> Q_s$, it shows either a mild peak, or a plateau with height
$R_{pA}\gtrsim 1$. The rapid increase at low $k_\perp\ll Q_s$ is a consequence of gluon saturation in
the nucleus, as quite clear by inspection of the calculation in \eqn{RpAMV}.

From the lessons at RHIC, we also know that it is interesting to study the evolution of  the $R_{pA}$ ratio
with increasing rapidity. Fig.~\ref{fig:Cronin_b} shows our respective predictions for the most
favourable case $\bar Q^2/Q_s^2=0.01$, as obtained from
numerical solutions to the BK equation up to $\eta=3$. Remarkably, one can see that the peak
flattens out already after one step in rapidity and that the ratio  $R_{pA}$ remains strictly below 1
 at any  $k_\perp^2\le 10$~GeV$^2$ for $\eta\ge 2$. A similar behaviour was experimentally observed
 in d+Au collisions at RHIC \cite{Arsene:2004ux,Adams:2006uz}
 and theoretically explained (within the CGC framework)
 in Refs.~\cite{Albacete:2003iq,Kharzeev:2003wz,Iancu:2004bx}. Clearly, the same
 physical mechanism is working in the present set-up: for the values of $k_\perp$ considered
 in Fig.~\ref{fig:Cronin_b}, the nucleus is still close to saturation, hence it evolves only slowly, unlike
 the proton, which is in a dilute regime, so its gluon distribution rises much faster with increasing $\eta$.
 This difference transmits to the respective dipole amplitudes and hence to the two cross-sections
 in the numerator and the denominator of  \eqn{RpA}, respectively.

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/RpA_total_cs}
    \caption{}\label{fig:Cronin_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/RpA_total_cs_rapidity}
    \caption{}\label{fig:Cronin_b}
  \end{subfigure}
  \caption{\small The ratio $R_{pA}$ defined in Eq.~\eqref{RpA}, plotted as a function of $k_{\perp}^2$.
  (a) The respective MV model predictions
   for 3 different values of $\bar{Q}^2$. For sufficiently small $\bar{Q}^2$ the ratio goes well above 1 at mid transverse momenta. (b) Rapidity evolution of the results from the left plot corresponding to
   $\bar{Q}^2=0.02$~GeV$^2$.
   After a couple of rapidity steps the ratio stays below 1 for all values of $k_\perp^2$ under
   consideration.}\label{fig:Cronin}
\end{figure}

\section{$k_{\perp}$-integrated SIDIS and saturation}
\label{sect:intk}

So far we have studied the double differential cross-section
$\frac{{\rm d}\sigma^{\gamma^* \!A \to q X}}{{\rm d} z\, {\rm d}^2\bm{k}}$ for producing a quark with a given
longitudinal momentum fraction $z$ and a given transverse momentum $\bm{k}$ in DIS processes with
relatively low effective virtualities $\bar Q^2=z(1-z)Q^2\ll Q_s^2$. Since this differential quantity may be difficult to
measure in the actual experiments, in what follows we shall discuss its more inclusive version
${{\rm d}\sigma}/{{\rm d} z}$ where the transverse momentum distribution of the
produced quark is integrated out --- only its longitudinal fraction $z$ is measured. We shall study in more
detail two kinematical limits: \texttt{(i)} $\bar Q^2\ll Q_s^2$, which is the most interesting regime
for a study of saturation, and  \texttt{(ii)} $\bar Q^2\gg Q_s^2$, where saturation can manifest itself
{\it indirectly}, via the phenomenon of geometric scaling. But before that, let us briefly list the
general formulae valid for any kinematics.

\subsection{General expressions}

Within our current formalism, the fact of integrating out $\bm{k}$ simply amounts to identifying the transverse coordinates, $\bm{r}$ and $\bm{r}'$, of the measured quark in the DA and in the CCA, respectively. So,
it is very easy to deduce the general respective expressions for the elastic, inelastic and total cross-sections.
To avoid a proliferation of formulae, we consider the transverse sector alone, for which we find
(for the reduced cross-sections defined in Eqs.~\eqref{JL}--\eqref{JT})
\begin{align}
\label{intJel}
	\int {\rm d}^2 \bm{k}\, \mathcal{J}_{{\rm T},{\rm el}} =
	\bar{Q}^2
	\int \frac{{\rm d}^2 \bm{r}}{4\pi} \,
	{\rm K}_1^2(\bar{Q}r) \,|T(\bm{r})|^2,
\end{align}
for the elastic piece,
\begin{align}
\label{intJin}
	\int {\rm d}^2 \bm{k}\, \mathcal{J}_{{\rm T},{\rm in}} =
	\bar{Q}^2
	\int \frac{{\rm d}^2 \bm{r}}{4\pi} \,
	{\rm K}_1^2(\bar{Q}r) \left[T(\bm{r}) + T(-\bm{r})- |T(\bm{r})|^2\right],
\end{align}
for the inelastic one, and, finally,
\begin{align}
\label{intJT}
	\int {\rm d}^2 \bm{k}\, \mathcal{J}_{\rm T} =
	\bar{Q}^2
	\int \frac{{\rm d}^2 \bm{r}}{4\pi} \,
	{\rm K}_1^2(\bar{Q}r)  \left[T(\bm{r}) + T(-\bm{r})\right],
\end{align}
for the total cross-section. We have also used the fact that $T(-\bm{r})=T^*(\bm{r})$, as it can be
easily checked by using $T=1-S$ together with the definition \eqref{dipole}
of the dipole $S$-matrix\footnote{$T$ and $S$ are real for all the explicit calculations in this paper,
but in general they can include an imaginary part describing the exchange of an object
with negative charge parity (``odderon''). In pQCD, such $C$-odd exchanges require al least
three gluons; see e.g. the discussions in \cite{Kovchegov:2003dm,Hatta:2005as}.}
. The $T$-matrix structure in \eqn{intJin}  can be equivalently written as
\begin{eqnarray}
T(\bm{r}) + T(-\bm{r})- |T(\bm{r})|^2 \,=\,1 - |S(\bm{r})|^2\,,
\end{eqnarray}
where the r.h.s. is recognised as the probability for the dipole with size $\bm{r}$ to suffer
an inelastic scattering (recall that $|S(\bm{r})|^2$ is the respective survival probability).
Similarly, \eqn{intJT} involves $T(\bm{r}) + T(-\bm{r})=2 {\rm Re}\,T(\bm{r})$, which is the expected
prediction of the optical theorem for the dipole total cross-section per unit impact parameter.

As  a further check, let us observe that by integrating \eqn{intJT} over $z$ one
recovers the expected, leading-order, result for the fully inclusive DIS cross-section in the dipole
factorisation  (here, for a transverse photon and for one flavour of massless quark):
\begin{equation}
\label{sigmat}
	\sigma_{\gamma_{{\rm T}}^*A} (x,Q^2)= e_q^2 \frac{\alpha_{\rm em} N_c}{2\pi^2}\,
			\int_0^1 {\rm d} z\int {\rm d}^2\bm{r}\,
	 \left[z^2+(1-z)^2\right] \bar{Q}^2
{\rm K}_1^2\big(\bar{Q}r\big)\,\sigma_{\rm dipole}(\eta, r),
	\end{equation}
with $\sigma_{\rm dipole}(\eta, r)=2\pi R_A^2T(\eta, r)$ when $T$ is real (as generally in this paper).

\subsection{Low virtuality $z(1-z)Q^2\ll Q_s^2$: the black disk limit}

It should be clear by now that the physics of strong scattering or gluon saturation can be most
directly investigated by focusing on the low-virtuality regime at $\bar Q^2\ll Q_s^2$. In this regime,
all the cross-sections shown in Eqs.~\eqref{intJel}--\eqref{intJT} are dominated by large dipoles
such that $1/ Q_s \ll r \ll 1/ \bar{Q}$, for which one can use the black disk limit $T(r)=1$ together
with the expansion of the Bessel function for small values of its argument. Indeed, this intermediate
range in $r$ yields a logarithmically enhanced contribution, which is easily found as
\begin{align}
\label{intJlow}
	\int {\rm d}^2 \bm{k}\, \mathcal{J}_{\rm T} \simeq
	2 \int {\rm d}^2 \bm{k}\, \mathcal{J}_{{\rm T},{\rm el}} \simeq
	2 \int {\rm d}^2 \bm{k}\, \mathcal{J}_{{\rm T},{\rm in}} \simeq
	\frac{1}{2} \ln \frac{Q_s^2}{\bar{Q}^2}
	\quad \textrm{for} \quad \bar{Q}^2 \ll Q_s^2.
\end{align}
Since obtained by invoking the black disk limit, these results are almost
insensitive to the details of the QCD scattering. In particular, they are independent
of the QCD coupling $\alpha_s$, except for the weak (logarithmic) dependence introduced
by the scale $Q_s^2$ in the argument of the logarithm (recall that $Q_s^2\sim \alpha_s^2$).
To better understand the physical origin of these results,
it is instructive to re-derive them by integrating over $\bm{k}$ the double-differential
cross-sections computed in Sect.~\ref{sect:qllqs}.

Concerning the elastic cross-section, it is quite clear that the respective logarithm in \eqn{intJlow}
is generated by integrating the $k_\perp$--distribution in \eqn{JTellow} over
the range $\bar{Q}^2 \ll k_{\perp}^2 \ll Q_s^2$. More precisely, the MV estimate in \eqn{JTelGBW}
shows that the logarithmic integration over $k_\perp^2$ is cut off by multiple scattering already
at $k_\perp^2\sim Q_s^2/2$. Hence the corresponding  estimate in  \eqn{intJlow}, although correct
to leading logarithmic accuracy, is expected to be larger than the actual result.

Consider now the inelastic cross-section: in this case, the logarithm  $\ln({Q_s^2}/{\bar{Q}^2})$
is already visible in the differential distribution in the second line of  \eqn{JTinlow}, or in its MV
model estimate \eqref{JTbroad}. By integrating any of these two expressions over all $k_\perp$,
one finds exactly the result displayed in \eqn{intJlow}. Of course, the Gaussian broadening
in  \eqn{JTbroad} only holds within a limited range of momenta around $Q_s$, but this is precisely
the range which dominates the $k_\perp$--integration to the accuracy of interest\footnote{One can
easily check that the power-like tail of the distribution at $k_\perp\gg Q_s$, e.g. $\mathcal{J}_{{\rm T}} \propto
Q_s^2/k_\perp^4$ in the MV model, merely yields a contribution of order 1 to the integrated cross-section.}.

These considerations are corroborated by the numerical results in Fig.~\ref{fig:int_a}, which show
the integrated elastic, inelastic, and total cross-sections in terms of  $\ln(\bar{Q}^2/ Q_s^2)$
in the saturation regime at $\bar{Q}^2 \ll Q_s^2$. Together with the predictions
of the MV model (the curves with $\eta=0$),
we also exhibit the results of the BK evolution up to $\eta=3$. For $\bar{Q}^2/Q_s^2 \lesssim 0.1$,
all the curves are roughly straight lines, with slopes which are consistent with \eqn{intJlow}.
As expected, the curves corresponding to elastic scattering lie below the inelastic ones,
but they are still comparable in so far as orders of magnitude are concerned.

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/integrated_lowkT}
    \caption{}\label{fig:int_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/integrated_highkT}
    \caption{}\label{fig:int_b}
  \end{subfigure}
  \caption{\small SIDIS cross-sections integrated over $\bm{k}$ and their high-energy evolution up to $\eta=3$:
    the elastic, inelastic and total cross-sections are plotted in terms of $\bar{Q}^2/Q_s^2$ and
  for 4 values of the rapidity $\eta$. (a) Logarithmic increase at low $\bar{Q}^2$: elastic and inelastic
  cross-sections are still comparable with each other. (b) Power law behaviour at high $\bar{Q}^2$: inelastic scattering strongly dominates the total cross-section.}\label{fig:int}
\end{figure}

\subsection{High virtuality $z(1-z)Q^2\gg Q_s^2$: geometric scaling}
\label{sec:highbarQ}

Even though not so directly relevant for studies of saturation (since controlled by weak scattering; see below),
the case $\bar{Q}^2 \gg Q_s^2$ is still interesting, notably for studies of the high-energy evolution
in the presence of saturation. Also, it features a strong $z$ -dependence of the cross-section, to be
discussed in Sect.~\ref{sec:zdep}.

In this regime, the dominant contributions to the integrals in Eqs.~\eqref{intJel}--\eqref{intJT} come
from very small dipoles with $r \ll 1/Q_s$, for which one can use the single
scattering approximation. Consider first the low energy case, where one can use
the MV model, cf. Eq.~\eqref{MVT0}, to deduce
\begin{align}
\label{intJhighMV}
        \int {\rm d}^2 \bm{k}\ \mathcal{J}_{{\rm T}} \simeq
       \bar{Q}^2
      \int
	\frac{{\rm d}^2 \bm{r}}{2\pi} \,
	{\rm K}_1^2(\bar{Q}r)\,
	 \frac{r^2 Q_{{\rm A}}^2}{4}
	 \ln \frac{1}{r^2\Lambda^2}\,.
\end{align}
In principle, one should cut the above integration at $r \sim 1/Q_s$, but since the Bessel function
is exponentially decaying for $\bar{Q}r\gg Q_sr\gg 1$, we can safely extend the integration back to infinity. The angular integration is trivial and  keeping only the logarithmically enhanced term, we find
\begin{align}
\label{intJhighMV2}
        \int {\rm d}^2 \bm{k}\ \mathcal{J}_{{\rm T}} \simeq
\frac{1}{6}\,
        \frac{Q_{\rm A}^2}{\bar{Q}^2}\,
        \ln \frac{\bar{Q}^2}{\Lambda^2}
        \qquad \textrm{for} \quad \bar{Q}^2 \gg Q_s^2.
\end{align}
It is easy to check that this dominant contribution arises from dipole sizes $r\sim 1/\bar Q\ll 1/Q_s$,
thus justifying our above use of the weak scattering approximation. Up to a numerical factor,
the quantity $(Q_{A}^2/\alpha_s)\ln({\bar{Q}^2}/{\Lambda^2})$ represents the nuclear gluon distribution
per unit transverse area
generated by the valence quarks, as probed on a transverse resolution scale $\bar{Q}^2$.
So, \eqn{intJhighMV2} has a simple physical interpretation: the small projectile dipole with transverse size
$r\sim 1/\bar Q\ll 1/Q_s$  scatters with a probability $\sim\alpha_s$ off all the gluons with
transverse sizes similar to its own and that it encounters while crossing the nucleus --- those
within a longitudinal tube with transverse area $1/\bar Q^2$.

The same exercise can be repeated for the amplitude given the BK solution at  asymptotically high energies,
shown in Eq.~\eqref{tbkdil}. One finds
\begin{align}
\label{intJhigh}
	\int {\rm d}^2 \bm{k}\, \mathcal{J}_{\rm T} \simeq
	\int {\rm d}^2 \bm{k}\, \mathcal{J}_{{\rm T},{\rm in}} \approx

	\left(\frac{Q_s^2}{\bar{Q}^2}\right)^\gamma
	\ln \frac{\bar{Q}^2}{Q_s^2}
	\qquad \textrm{for} \quad \bar{Q}^2 \gg Q_s^2,
\end{align}
where we omitted the overall  proportionality constant since this is anyway not under control when writing
the BK amplitude as in \eqn{tbkdil}.  We have here anticipated that the elastic cross section is more
strongly suppressed in this limit. One can indeed check that
\begin{align}
\label{intJelhigh}
	\int {\rm d}^2 \bm{k}\, \mathcal{J}_{{\rm T},{\rm el}} \approx

	\left(\frac{Q_s^2}{\bar{Q}^2}\right)^{2\gamma}
	\ln^2 \frac{\bar{Q}^2}{Q_s^2}
	\qquad \textrm{for} \quad \bar{Q}^2 \gg Q_s^2.
\end{align}

The high-energy and large-virtuality
results in Eqs.~\eqref{intJhigh}--\eqref{intJelhigh}  show two remarkable features:
\texttt{(i)}  the  ``anomalous'' dimension $\gamma <1$
specific to the high-energy evolution (BK/JIMWLK)
in the presence of saturation, and \texttt{(ii)} geometric scaling --- the fact that these expressions
depend upon the (effective) virtuality $\bar Q^2=z(1-z) Q^2$ and upon the COM energy of the scattering
only via the ratio ${\bar{Q}^2}/{Q_s^2(\eta)}$. Notice that geometric scaling also holds at low virtualities
${\bar{Q}^2}\ll {Q_s^2(\eta)}$, at least within the black disk approximation in \eqn{intJlow}.  So, for
sufficiently high energies (or large $\eta$), this is expected to be a robust feature of the
integrated SIDIS cross-sections at any $\bar{Q}^2$.

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/int_scaling_lowkT}
    \caption{}\label{fig:scaling_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/int_scaling_highkT}
    \caption{}\label{fig:scaling_b}
  \end{subfigure}
  \caption{\small Testing geometric scaling for the SIDIS transverse cross-section integrated over $\bm{k}$,
  from solutions to the collinearly improved rcBK up to $\eta=6$.
  (a) Low $\bar{Q}^2$: the scaling is quasi-perfect. (b) High $\bar{Q}^2$:
   with increasing $\eta$, the power law tail softens and the scaling becomes more accurate.}\label{fig:intscaling}
\end{figure}

We have checked these properties against numerical solutions to the BK equation,
 with the results shown in  Fig.~\ref{fig:int_b}. With increasing $\eta$, one clearly sees
a softening of the distribution at high $\bar Q^2$ as compared to the low-energy (MV model) result
\eqn{intJhighMV2}. This softening reflects the emergence of the  anomalous dimension, although
 the values of $\eta\le 3$ considered there are still to small to observe the asymptotic value
 $\gamma \simeq 0.63$ expected at very high energies.

 To also check geometric scaling, we have plotted in
 Fig.~\ref{fig:intscaling} the results for the integrated cross-section
 $\int {\rm d}^2 \bm{k}\ \mathcal{J}_{{\rm T}} $ in terms of the scaling
variable $\bar{Q}^2/ Q_s^2(\eta)$, up to the larger rapidity $\eta=6$.
In the low virtuality regime at ${\bar{Q}^2}\ll {Q_s^2(\eta)}$, one finds a quasi-perfect scaling
(see Fig.~\ref{fig:scaling_a}),
in agreement with the leading-order estimate in  Eq.~\eqref{intJlow}.
At higher virtualities, ${\bar{Q}^2}\gtrsim {Q_s^2(\eta)}$, the scaling becomes better
with increasing rapidity (see Fig.~\ref{fig:scaling_b}) --- again, as expected.

\subsection{The $z$-dependence of the integrated cross-sections}
\label{sec:zdep}

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/z-dependence}
    \caption{}\label{fig:zdep_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[width=\textwidth]{Figs/z-dependence-5}
    \caption{}\label{fig:zdep_b}
  \end{subfigure}
  \caption{\small The normalized $\bm{k}$-integrated transverse cross section-defined in \eqn{sigmaz} plotted
   as a function of $z$ for various values of $\eta$ at (a) $Q^2=Q_{s0}^2$ (``low $\bar{Q}^2$'')
   and (b) $Q^2=5Q_{s0}^2$ (``high $\bar{Q}^2$''). When $z$ approaches 1, one sees
   a sharper increase in Fig.~\ref{fig:zdep_b} than in Fig.~\ref{fig:zdep_a}.}\label{fig:zdep}
\end{figure}

Yet a different way of looking at the results for the integrated cross-sections
is to plot them in terms of the longitudinal momentum fraction $z$ of the measured quark, for
fixed values of the virtuality $Q^2$ and of the target saturation momentum ${Q_s^2(\eta)}$.
In  Fig.~\ref{fig:zdep} we show our respective numerical results for two choices for the virtuality,
 $Q^2=Q_{s0}^2$ (left plot) and respectively  $Q^2=5Q_{s0}^2$ (right plot), and for $\eta\le 3$.
 The plotted quantity is the ratio
 \begin{eqnarray}\label{sigmaz}
 \frac{\sigma_{{\rm T}}(z)}{\sigma_{{\rm T}}(z=0.5)}\,=\,\frac{\int {\rm d}^2 \bm{k}\, \mathcal{J}_{\rm T}(\bm{k}, z)}
 {\int {\rm d}^2 \bm{k}\, \mathcal{J}_{\rm T}(\bm{k}, z=0.5)}\,,\end{eqnarray}
 where $\sigma_{{\rm T}}(z)$ is merely a short-hand notation for the respective
 differential cross-section per unit of $z$, that is ${{\rm d}\sigma}_{{\rm T}}/{{\rm d} z}$.
As compared to  Fig.~\ref{fig:int},  we now implicitly concentrate on a more limited range of values for
 the ratio $\bar Q^2/Q_{s0}^2$, which are more realistic from the viewpoint of the phenomenology.
This new way of plotting the results is interesting in that it emphasises their strong dependence
upon $z$ in the vicinity of $z=1$, a feature that could be observed in the experiments.

This strong $z$--dependence can be understood from our analytic approximations
in Eqs.~\eqref{intJlow} and \eqref{intJhigh}. These expressions show that,
when increasing $z$ towards 1, $\sigma_{{\rm T}}(z)$ increases like $\ln\frac{1}{1-z}$ so long as
${\bar{Q}^2}\ll {Q_s^2(\eta)}$ and (roughly) like the power $1/(1-z)^\gamma$ when
${\bar{Q}^2}\gg {Q_s^2(\eta)}$. For $z\gtrsim 0.8$, the two plots in Fig.~\ref{fig:zdep}
are representative for these two limiting kinematical regimes. And indeed, all the curves there show a rather
pronounced increase with increasing $z$ above 0.8. This increase is stronger in Fig.~\ref{fig:zdep_b}
(``high $\bar Q^2$'') than in Fig.~\ref{fig:zdep_a}  (``low $\bar Q^2$''), and it becomes less pronounced
when increasing $\eta$, thus reflecting the onset of an anomalous dimension $\gamma < 1$.

\subsection{$R_{pA}$ ratios for the integrated cross-section}

In Sect.~\ref{sec:cronin} we observed that a suggestive way to visualise the nuclear modifications
is by plotting the ratio $R_{pA}$ between the SIDIS cross-sections in $eA$ collisions and in $ep$ collisions
(times the nucleon number $A$), respectively. Previously, we studied this ratio for
the case of the differential cross-section in $k_\perp$, but a similar ratio can of course be defined
for the cross-sections integrated over $k_\perp$.

\begin{figure}
  \centering
  \begin{subfigure}[t]{0.47\textwidth}
    \includegraphics[width=\textwidth]{Figs/integrated_x_rpa_qbar}
    \caption{}
    \label{fig:Rpaint_a}
  \end{subfigure}
  \hfill
  \begin{subfigure}[t]{0.49\textwidth}
    \includegraphics[width=\textwidth]{Figs/rpa_z-dependence}
    \caption{}
    \label{fig:Rpaint_b}
  \end{subfigure}
  \caption{\small (a) The ratio $R_{pA}$ for the  SIDIS cross-section integrated over $\bm{k}$,
  as defined in Eq.~\eqref{RpAz}, and its rapidity evolution of up to $\eta=3$.
   The ratio is always smaller than 1 in contrast to the unintegrated case (cf.~Fig.~\ref{fig:Cronin}). (b)
   The nuclear modification factor for the normalised cross-section introduced  in \eqn{sigmaz}
   plotted as a function of $z$, for fixed $Q^2$ and for various values of $\eta\le 3$. As $z$ approaches
   1 the ratio drops significantly below 1, due to the different $z$-dependences observed in
   Figs.~\ref{fig:zdep_a} and \ref{fig:zdep_b}, respectively.}\label{fig:Rpaint}
\end{figure}

Specifically, in  Fig.~\ref{fig:Rpaint_a} we display the ratio
\begin{eqnarray}\label{RpAz}
 R_{pA}(z)\,\equiv\,
 \frac{\int {\rm d}^2 \bm{k}\,\mathcal{J}_{\rm T}^A(z,\bm{k})}{A^{1/3}\int {\rm d}^2 \bm{k}\,\mathcal{J}_{\rm T}^p(z,\bm{k})}\,,
 \end{eqnarray}
 as a function of $\bar Q^2=z(1-z)Q^2$ up to $\bar Q^2=5Q_{s0}^2(A)=10$~GeV$^2$ and
 for different rapidities up to $\eta=3$. This ratio would be
 equal to one in the absence of nuclear effects, but by inspection of  Fig.~\ref{fig:Rpaint_a}
 one sees that this is clearly not the case: the ratio is strictly smaller than one and monotonously
 increasing with $\bar Q^2$ (albeit almost flat for  $\bar Q^2 \gtrsim Q_{s0}^2$) within the whole
  range under consideration.
   Moreover, the nuclear suppression is significantly increasing with $\eta$ at any value for $\bar Q^2$.

 This behaviour can be easily understood, at least qualitatively. At the level of the MV model,
 it is well known that the effect of gluon saturation in a large nucleus is to push the gluon
 distribution towards large transverse momenta
 \cite{Kharzeev:2003wz,Iancu:2004bx}: the low momentum modes at $k_\perp< Q_s$
 get depleted and the gluons which are missing at low $k_\perp$ are redistributed in the modes
 with larger momenta  $k_\perp\gtrsim Q_s$.
Hence, if one integrates the  $k_\perp^2$-distribution up to some value $\bar Q^2$, the total
number of gluons thus counted is smaller than it would have been in the absence of saturation.
Since the saturation effects are much smaller in a proton, the
$R_{pA}$ ratio for the integrated gluon distributions is smaller than one, but it rises towards
one when increasing $\bar Q^2$. With increasing energy (or $\eta$), the gluon distribution rises
faster in the proton than in the nucleus, hence the $R_{pA}$ ratio  becomes even smaller.

This is precisely the kind of behaviour one sees in Fig.~\ref{fig:Rpaint_a}. Strictly speaking, this plot does
not apply to the integrated gluon distribution in the target, but rather to the integrated
cross-section for quark production in DIS. Yet, these two quantities are closely related to each other.
This is visible
e.g. at the level of  \eqn{JTinlow}, where the Fourier transform of the dipole $S$-matrix is
recognised as the ``dipole'' version of the unintegrated gluon distribution in the target
(the ``dipole TMD'') \cite{Marquet:2009ca,Dominguez:2011wm}.

Furthermore, Fig.~\ref{fig:Rpaint_b} shows the nuclear modification factor for the quantity depicted in
Fig.~\ref{fig:zdep}, which is itself a ratio, cf. \eqn{sigmaz}. That is, the quantity depicted
in Fig.~\ref{fig:zdep_a} for the case of a nucleus is divided by the same quantity computed for
a proton\footnote{Note that, in forming this ratio, there is no  factor of  $A$
in the denominator. Indeed, the global proportionality with $A$ of the cross-section
corresponding to $eA$ collisions has already cancelled in the ratio \eqref{sigmaz}.}.
In both calculations, one uses $Q^2=Q_{s0}^2(A)=2$~GeV$^2$. Recalling that $Q_{s0}^2(A)=
A^{1/3}Q_{s0}^2(p)$ and $A^{1/3}=6$, it is easy
to see that the quantity appearing in the denominator of this ratio
 (the one corresponding to the proton) should behave similarly
to that depicted  Fig.~\ref{fig:zdep_b}. In other terms, each of the curves
shown in Fig.~\ref{fig:Rpaint_b} can be thought of as the ratio of the 2 corresponding curves
in Fig.~\ref{fig:zdep_a} and
Fig.~\ref{fig:zdep_b},  respectively. From this perspective, the behaviour visible in  Fig.~\ref{fig:Rpaint_b}
is easy to understand: when $Q^2=Q_{s0}^2(A)=6Q_{s0}^2(p)$,
one probes a dense gluon distribution in the nuclear, but a dilute one in the proton. Hence,
the cross-section $\sigma_{{\rm T}}(z)$ increases much slower when increasing $z$ towards 1
for the nucleus than for the proton (as also visible by comparing the two plots in Fig.~\ref{fig:zdep}).
Accordingly, the nucleus-to-proton ratio is decreasing, as shown in Fig.~\ref{fig:Rpaint_b}.
Interestingly, this suppression is quite strong already for not too large values
of $z$ (say $z\simeq 0.8$) and thus it might be a robust signal of saturation in the EIC data.

 \section{Inclusive DIS: symmetric vs. aligned-jet configurations}
\label{sect:sym-al}

In most of the previous developments in this paper, we assumed that the longitudinal momentum fraction
$z$ of the measured quark is very close to one, in such a way that $\bar{Q}^2 \equiv z(1-z)Q^2 \ll Q_s^2$.
Asymmetric dipole configurations, such that
$z$ is either very close to one, $1-z\ll 1$, or very small $z\ll 1$, have often been
discussed in a different context, namely in relation with the so-called
``aligned jet'' contributions to the inclusive DIS cross-section
$\sigma_{\rm T}\equiv
\sigma_{\gamma_{{\rm T}}^*A} (x,Q^2)$, cf. \eqn{sigmat}. So, it is natural to wonder what is the
relation between these two physical regimes, if any. In this section, we would like to explain that
these situations are in fact different and do not overlap with each other. The subsequent discussion will
also make clear that the large-$z$ kinematics relevant for this paper corresponds to {\it rare}
configurations, which brings only small contributions to the overall DIS cross-section.

Before we discuss the inclusive DIS cross-section, in which the longitudinal fraction
$z$ is integrated over, let us first summarise
the interesting kinematical regimes for fixed values of $Q^2$, $z$ and $Q_s^2$:

 \texttt{(i)} Low  virtuality $Q^2\ll Q_s^2$: in this case, one clearly has $\bar{Q}^2\ll Q_s^2$ for any value of $z$.
 In this regime, saturation effects are certainly important, but in practice they may mix with
 non-perturbative effects like vector meson fluctuations of the virtual photon, due to the fact that
 such low values for $Q^2$ are close to the QCD confinement scale $\Lambda^2$.

 \texttt{(ii)} High virtuality $Q^2\gg Q_s^2$, but such that $z$ is sufficiently small, or sufficiently close to one,
 in order to have  $\bar{Q}^2 \ll Q_s^2$. This is the kinematical regime that we almost exclusively considered
 throughout this paper, in the particular context of SIDIS and for $z\simeq 1$.

 \texttt{(iii)} High virtuality  $Q^2\gg Q_s^2$ and $z$ values such that $\bar Q^2 \gg Q_s^2$ as well.
 Clearly, this regime is not very interesting for studies of gluon saturation (except marginally, via effects
 like geometric scaling and the BK anomalous dimension discussed in Sect.~\ref{sec:highbarQ}).
In turn, this regime includes two different situations, depending upon the values of $z$: (a) aligned-jet configurations, which are such that $1\gg z(1-z)\gg Q_s^2/Q^2$,  and (b) symmetric configurations,
for which $z$
is neither very close to zero, nor very close to 1; we shall succinctly refer to them as $z\sim 1/2$.

It is generally stated that the aligned-jet configurations
give the dominant contribution to $\sigma_{\rm T}$ at high $Q^2$,
but as we shall shortly see this conclusion can be modified by the high-evolution evolution
and the associated anomalous dimension, cf. \eqn{tbkdil}.

 We now present leading-order, parametric,  estimates for the inclusive DIS
 cross-section $\sigma_{\rm T}$. For convenience, we shall use the simple notation
 $ \tilde\sigma_{\rm T}$ for the cross-section in \eqn{sigmaT} {\it without} the factor $R^2_A N_c \alpha_{\rm em} e_q^2$, that is, without the constant geometrical, colour and electromagnetic factors. Still, we shall take into account the numerical constant and the $z$-dependent factor. We successively consider the low-virtuality and the high-virtuality regimes and
distinguish between contributions from various regions in $z$.

\texttt{(i)} When $Q^2 \ll Q_s^2$, one can use the ``black-disk'' approximation in \eqn{intJlow} for any value of $z$,
to deduce
\begin{align}
        \label{sigmalow}
         \tilde\sigma_{\rm T} \simeq
        \int\limits_0^1 {\rm d} z \left[z^2 + (1-z)^2 \right]
        \left[ \ln \frac{Q_s^2}{Q^2} + \ln \frac{1}{z(1-z)}
        \right]
        \simeq \frac{2}{3}\ln \frac{Q_s^2}{Q^2}
        \qquad \textrm{for} \quad Q^2 \ll Q_s^2,
\end{align}
where only the first term in the square brackets is important to logarithmic accuracy. Clearly, this
dominant contribution comes from the symmetric configurations with  $z\sim 1/2$.

 \texttt{(ii)} When $Q^2\gg Q_s^2$, but $z$ is restricted to the corners in the phase-space where
 $z(1-z)\lesssim Q_s^2/Q^2$, we can still rely on \eqn{intJlow}, but with appropriate limits on the
integral over $z$; namely, instead of \eqn{sigmalow}, one should write
 \begin{align}
        \label{sigmahighsat}
         \tilde\sigma_{\rm T} \simeq
        2 \int\limits_0^{Q_s^2/Q^2} {\rm d} z
        \left[ \ln \frac{1}{z} - \ln \frac{Q^2}{Q^2_s}
        \right]
        =\, \frac{2 Q_s^2}{Q^2}
        \qquad \textrm{for} \quad Q^2 \gg Q_s^2 \gtrsim \bar{Q}^2,
\end{align}
where we integrated over small values $z \le Q_s^2/Q^2$ and multiplied the result by a factor of 2
to account for the region at $1-z\le Q_s^2/Q^2$. (Notice however that the upper limit in Eq.~\eqref{sigmahighsat} should involve a fudge factor of order $\mathcal{O}(1)$, thus the proportionality constant is not really under control.) This is the contribution to the inclusive
cross-section from the SIDIS configurations which have been our main focus in this paper --- those
where the dipole fluctuations have very large sizes $r\sim 1/\bar Q\gtrsim 1/Q_s$ and therefore
suffer strong scattering.

\texttt{(iii)} When both $Q^2 \gg Q_s^2$ and $\bar Q^2 \gg Q_s^2$, the dipole scattering is weak
($r\sim 1/\bar Q\ll 1/Q_s$) and we can use the results for ${{\rm d}\sigma_{\rm T}}/{{\rm d} z}$
from Sect.~\ref{sec:highbarQ},  but with the subsequent integration over $z$ restricted to
$z(1-z) \gg Q_s^2/Q^2$.

Consider first the  ``low-energy'' case, as described by the MV model; using
 \eqn{intJhighMV2}, one finds
 \begin{align}
        \label{sigmahighMV}
         \tilde\sigma_{\rm T} \simeq
        \frac{2}{3}\, \frac{Q_{\rm A}^2}{Q^2}
        \int\limits_{Q_s^2/Q^2}
        \frac{{\rm d} z}{z} \left(\ln \frac{Q^2}{\Lambda^2} + \ln z \right)
        \simeq \frac{2}{3}\,
        \frac{Q_{\rm A}^2}{Q^2}\, \ln\frac{Q^2}{\Lambda^2}\,
        \ln \frac{Q^2}{Q_s^2}\,
        \quad \textrm{for} \quad Q^2 \gg Q_s^2 \quad \textrm{(MV model)},
\end{align}
 where in the final result we kept only the term
generated by the logarithmic integration over $z$, which
is enhanced by the large logarithm $\ln ({Q^2}/{Q_s^2})$. (A similar contribution from $z\sim 1$
has been included via a factor of 2.)
Clearly, this dominant contribution comes from the {\it aligned-jet configurations},
as previously defined --- that is, $z$ is small, but not {\it too} small (and similarly for $1-z$), in such
a way that $1\gg z(1-z)\gg Q_s^2/Q^2$.

Consider now the high-energy regime $\bar{\alpha}_s\eta \gtrsim 1$,
where  the BK evolution becomes important. Using Eq.~\eqref{intJhigh}
and neglecting proportionality constants, we obtain
\begin{align}
        \label{sigmahigh}
        \!\!\!\tilde \sigma_{\rm T} \approx
        2 \left(\frac{Q_s^2}{Q^2}\right)^\gamma\!\!\!
        \int\limits_{Q_s^2/Q^2}^{1}\!\! {\rm d} z\,
        \frac{z^2+(1-z)^2}{[z(1-z)]^\gamma}
        \left[ \ln \frac{Q^2}{Q_s^2} + \ln z(1-z)\right]
        \approx
        \left(\frac{Q_s^2}{Q^2}\right)^\gamma \!\ln \frac{Q^2}{Q_s^2}
        \quad \textrm{for} \quad Q^2 \gg Q_s^2.
\end{align}
Since $\gamma < 1$, the above integral is dominated by large values\footnote{Unlike in
\eqn{sigmahighMV}, the logarithm $ \ln({Q^2}/{Q_s^2}) $ which appears in
 the final result in \eqn{sigmahigh}
is not generated by the integral over $z$, rather it comes from the transverse logarithm
originally present in the BK solution for the dipole amplitude at high energy, cf.
 \eqn{tbkdil}.} $z\sim 1/2$ (in particular, the
lower limit at $z= Q_s^2/Q^2\ll 1$ played no role for the final result), i.e.~by {\it symmetric} dipole
configurations, with very small sizes $r\sim 1/Q$. This change of behaviour w.r.t. the MV model
calculation in Eq.~\eqref{sigmahighMV} is due to the
fact that, when $\gamma < 1$, the dipole amplitude for a single
scattering rises more slowly with the dipole size $r$
than in the regime of ``colour transparency'' ($\gamma=1$).

The results in Eqs.~\eqref{sigmahighMV} and \eqref{sigmahigh} confirm our earlier statement
in this section, namely the fact that the aligned-jet configurations dominate the DIS
cross-section $\sigma_{\rm T}$ at high $Q^2$ only so long as one can neglect the onset
of the BK anomalous dimension.
This is the actual physical situation at not too high energies (in particular, in the MV model and,
more generally, in the regime where the DGLAP evolution is the dominant QCD dynamics).
It is furthermore the case at {\it any} energy, provided $Q^2$ is sufficiently high $Q^2$  ---
 larger than the upper limit of the geometric scaling region (which, via
the high-energy evolution, grows with $\eta$ even faster than $Q_s^2$)
 \cite{Iancu:2002tr,Mueller:2002zm,Iancu:2004bx}).

The above results also show that, irrespective of the precise value of $\gamma\le 1$,
inclusive DIS at high $Q^2\gg Q_s^2$ is controlled by relatively small dipoles undergoing weak scattering
($r\ll 1/Q_s$). The corresponding contribution of the large dipoles with $r\gtrsim 1/Q_s$,
as estimated in \eqn{sigmahighsat}, is parametrically smaller. This confirms another statement that
we made earlier, namely, the fact that the SIDIS configurations which are most interesting for
studies of saturation yield only a negligible contribution to the inclusive cross-section.

Before we conclude this section, let us notice that, when $\gamma < 1$,
 the parametric estimate shown in \eqn{sigmahigh} for $\sigma_{\rm T}$ also holds for the
 longitudinal cross-section $\sigma_{\rm L}$; in that case too, this result is generated by
 symmetric dipole configurations with $z\sim 1/2$ and $r\sim 1/Q$.  But unlike $\sigma_{\rm T}$,
 the longitudinal cross-section is controlled by symmetric  configurations even
when $\gamma=1$, so in that case $\sigma_{\rm L}$ is lacking the large
logarithm $ \ln({Q^2}/{Q_s^2}) $ visible in \eqn{sigmahighMV} and thus is suppressed relative
to $\sigma_{\rm T}$.

\comment{
\noindent \hrulefill

\noindent \hrulefill

Here I rewrite \eqref{sigmahighsat} including the $\gamma$-dependent prefactor to show the singularity in the $\gamma \to 1$ limit.
\begin{align}
        \label{testt1}
        \tilde \sigma_{\rm T}
        \approx
        \frac{4 \Gamma(1-\gamma) \Gamma(3-\gamma)}{\Gamma(4-2\gamma)}
        \left(\frac{Q_s^2}{Q^2}\right)^\gamma \!\ln \frac{Q^2}{Q_s^2}
        \quad \textrm{for} \quad Q^2 \gg Q_s^2.
\end{align}
I also add the respective formulae for  $\sigma_{{\rm L}}$ (including the prefactor $8 z(1-z)$ in Eq.~\eqref{sigmaL}).
\begin{align}
        \label{testl1}
         \tilde\sigma_{\rm L} \simeq
        8 \int\limits_0^1 {\rm d} z\, \frac{ z  (1-z)}{2}
        \simeq \frac{2}{3}
        \qquad \textrm{for} \quad Q^2 \ll Q_s^2,
\end{align}
so as expected and in contrast to $\tilde{\sigma}_{\rm L}$, there is no large logarithm (no logarithmic integration in $k_\perp^2$).
\begin{align}
        \label{testl2}
         \tilde\sigma_{\rm L} \simeq
        8 \!\!\int\limits_{Q_s^2/Q^2}^1\!\!
        {\rm d} z\, z (1-z)\,\frac{Q_A^2}{12 \bar{Q}^2}
        \ln \frac{\bar{Q}^2}{\Lambda^2}\,
        \simeq   \frac{2}{3}\,
        \frac{Q_A^2}{Q^2}\,
        \ln \frac{Q^2}{\Lambda^2}\,
        \qquad \textrm{for} \quad Q^2 \gg Q_s^2 \quad \textrm{(MV model)},
\end{align}
and again one of the logarithms is missing compared to the transverse case (no logarithmic integration in $z$).
\begin{align}
        \label{testl3}
        \!\!\!\tilde \sigma_{\rm L} \approx
        8 \left(\frac{Q_s^2}{Q^2}\right)^\gamma\!\!\!
        \int\limits_{Q_s^2/Q^2}^{1}\!\!
        \frac{{\rm d} z}{[z(1-z)]^{\gamma-1}}
        \ln \frac{Q^2z(1-z)}{Q_s^2}
        \approx
        \frac{8 \Gamma(2-\gamma)^2}{\Gamma(4-2\gamma)}
        \left(\frac{Q_s^2}{Q^2}\right)^\gamma \!\ln \frac{Q^2}{Q_s^2}
        \quad \textrm{for} \quad Q^2 \gg Q_s^2.
\end{align}
It is of the same order as \eqref{testt1} when $\gamma>1$, but there is no singularity as $\gamma \to 1$.

\noindent \hrulefill

\noindent \hrulefill
}

\section{Conclusions}
\label{sec:conc}

In this paper, we have identified a new physical regime that should be interesting for studies
of unitarity corrections and gluon saturation in electron-nucleus deep inelastic scattering at the
future Electron-Ion Collider. In the dipole picture for DIS at small Bjorken $x$, this regime
corresponds to very asymmetric quark-antiquark fluctuations, such that one of the two quarks
carries a large fraction $z\simeq 1$ (meaning $1-z\ll 1$) of the longitudinal momentum of its
parent virtual photon. This longitudinal asymmetry implies that the dipole fluctuation has a
very large transverse size $r\sim 1/\bar Q\gg 1/Q$, with $\bar Q^2=z(1-z)Q^2$. When moreover
one has also has $\bar Q^2\lesssim Q_s^2$, with $Q_s$ the saturation momentum in the
nuclear target, then the dipole-nucleus scattering is strong and thus probes high gluon density
effects, like gluon saturation, in the target wavefunction.

We have proposed several observables in view of experimental studies of this regime.
They all refer to SIDIS experiments in which one measures the jet (or hadron) produced by the
final-state evolution of the quark carrying the large longitudinal momentum fraction $z$.

For experiments which measure both the longitudinal fraction $z$ and the transverse momentum
$k_\perp$ of the final jet (or hadron), we have argued that the most interesting kinematics
for a study of saturation is at $\bar Q^2\ll k_\perp^2\lesssim Q_s^2$.  We identified two
interesting phenomena in this regime.

For sufficiently low transverse momenta $k_\perp^2\ll Q_s^2$, we found that
the scattering probes the ``black disk limit'', i.e. the regime
where the dipole amplitude reaches the unitarity limit $T=1$. In this case,
the inelastic scattering is suppressed by unitarity, hence the total SIDIS cross-section
is controlled by its  elastic (or coherent)  component, as given by the ``shadow'' of the
strong scattering. In our dipole picture,
this corresponds to configurations where the splitting of the virtual photon
into a $q\bar q$ pair occurs {\it after} the photon has crossed the nucleus without
interacting with it.

On the other hand, for larger transverse momenta $k_\perp^2\sim Q_s^2$,
the dominant component is the inelastic scattering of the measured quark, leading to
transverse momentum broadening, that is, to a Gaussian distribution
which peaks at $k_\perp\sim Q_s$. A similar, ``Cronin'', peak could also be seen in the $R_{pA}$ ratio
between the SIDIS cross-sections in $eA$ and in $ep$ collisions, respectively, albeit only for
extremely small values of the ratio $\bar Q^2/Q^2_s$, that might be difficult to reach in practice.
Yet, the deviation of $R_{pA}$ from unity is found to be significantly large including
for more realistic kinematical conditions, so the nuclear effects should be easy to observe.
From the viewpoint of the target, the multiple scattering of the measured quark
probes the nuclear gluon distribution at, or near, saturation, as measured by a special unintegrated
gluon distribution known as the ``dipole TMD".

We have also studied the effects of the high energy evolution, as described by the collinearly-improved
version of the BK equation (with running coupling) recently proposed in
Refs.~\cite{Ducloue:2019ezk,Ducloue:2019jmy}. We found that, in general, these effects are quite small,
due to the slowing down of the evolution by the various (collinear and running-coupling) resummations,
and also to the fact that the we have only considered relatively small rapidity evolutions, in line
with the kinematics expected at the EIC. The only observable for which even a small rapidity evolution
may have a qualitative effect is the $R_{pA}$ ratio: the Cronin peak disappears after just one unit
of rapidity and is replaced by uniform nuclear suppression
up to relatively large transverse momenta, well above $Q_s$.

Since the transverse momentum of the tagged jet (or hadron) may be difficult to measure
in practice, especially at very forward rapidities, we also study the saturation effects on the
SIDIS cross-section integrated over  $k_\perp$. Once again, the interesting regime for saturation is at
$\bar Q^2\ll Q^2_s$. In this regime, we find that the (integrated) elastic and inelastic cross-sections
are comparable with each other --- they are even identical to leading logarithmic accuracy.
Moreover, they are both of order $\ln(Q^2_s/\bar Q^2)$, with the logarithm coming from
integrating over the wavefunction of the virtual photon, that is, over the $k_\perp$-distribution
of the $q\bar q$ pair produced by the decay of the virtual photon. For the elastic scattering, this is
also the {\it final} distribution. For the inelastic channel, the final distribution is generated via
multiple scattering, so the integral over $k_\perp$ is controlled by values $k_\perp\sim Q_s$.

To better emphasise the distinguished nature of this behaviour at low $\bar Q^2\ll Q^2_s$, we
have also computed the integrated SIDIS cross-section in the weak scattering regime at
$\bar Q^2\gg Q^2_s$. We thus found the expected ``leading-twist'' suppression by the factor
$Q^2_s/\bar Q^2$ (possibly modified by an anomalous dimension introduced by the high
energy evolution). A particularly suggestive way to present our results is by plotting the
integrated cross-section as a function of $z$ for a fixed and relatively large value of the virtuality
 $Q^2\gg Q_s^2$ (cf. Fig.~\ref{fig:zdep}):  in the weak scattering
regime at $\bar Q^2\gg Q_s^2$ one sees a rapid, power-like, growth when increasing $z$ towards 1,
which is eventually tamed (when $\bar Q^2\lesssim Q^2_s$) by saturation. This change of behaviour
can perhaps be better visualised by plotting the corresponding $R_{pA}$ ratio,
as shown in Fig.~\ref{fig:Rpaint_b}.

We have finally considered the inclusive DIS cross-section, with the purpose
of clarifying the relation between the special, large-$z$, SIDIS configurations that we have
investigated in this paper and the ``aligned jet'' configurations traditionally discussed in relation
with DIS at large $Q^2\gg Q_s^2$.
Although both types of configurations are associated with dipole fluctuations
which are very asymmetric in their $z$ sharing (i.e. $z(1-z)\ll 1$) and hence relatively large
($r\sim 1/\bar Q\gg 1/Q$), they nevertheless  correspond to
different physical regimes, which do not overlap with each other.
The ``aligned jets'' correspond to dipoles which are still small enough ($r\ll 1/Q_s$) to scatter
only weekly; these are the typical configurations which control inclusive DIS at high $Q^2$.
On the contrary, the ``large-$z$'' SIDIS configurations which are interesting for
saturation are sufficiently large to suffer strong scattering ($r\gtrsim 1/Q_s$)
even at high $Q^2\gg Q_s^2$. But they are
rare events, with only a small contribution to the inclusive cross-section.

Our present study is mostly preliminary. For applications to the phenomenology, it should be completed with an analysis of the final-state parton evolution into hadrons or jets --- an issue which is especially important for such large values of $z$ --- and also with a study of the experimental feasibility of hadron or jet measurements at very forward rapidities and semihard transverse momenta. We hope that our present analysis will give an incentive for further developments and will inspire new studies of gluon saturation in deep inelastic scattering at the EIC.

\section*{Acknowledgements}
The work of E.I.~is supported in part by the Agence Nationale de la Recherche project ANR-16-CE31-0019-01. The work of A.H.M.~is supported in part by the U.S.~Department of Energy Grant \# DE-FG02-92ER40699.

redpajama_set_name: RedPajamaArXiv