Saturation effects in SIDIS at very forward rapidities

Using the dipole picture for electron-nucleus deep inelastic scattering at small Bjorken $x$, we study the effects of gluon saturation in the nuclear target on the cross-section for SIDIS (single inclusive hadron, or jet, production). We argue that the sensitivity of this process to gluon saturation can be enhanced by tagging on a hadron (or jet) which carries a large fraction $z \simeq 1$ of the longitudinal momentum of the virtual photon. This opens the possibility to study gluon saturation in relatively hard processes, where the virtuality $Q^2$ is (much) larger than the target saturation momentum $Q_s^2$, but such that $z(1-z)Q^2\lesssim Q_s^2$. Working in the limit $z(1-z)Q^2\ll Q_s^2$, we predict new phenomena which would signal saturation in the SIDIS cross-section. For sufficiently low transverse momenta $k_\perp\ll Q_s$ of the produced particle, the dominant contribution comes from elastic scattering in the black disk limit, which exposes the unintegrated quark distribution in the virtual photon. For larger momenta $k_\perp\gtrsim Q_s$, inelastic collisions take the leading role. They explore gluon saturation via multiple scattering, leading to a Gaussian distribution in $k_\perp$ centred around $Q_s$. When $z(1-z)Q^2\ll Q^2$, this results in a Cronin peak in the nuclear modification factor (the $R_{pA}$ ratio) at moderate values of $x$. With decreasing $x$, this peak is washed out by the high-energy evolution and replaced by nuclear suppression ($R_{pA}<1$) up to large momenta $k_\perp\gg Q_s$. Still for $z(1-z)Q^2\ll Q_s^2$, we also compute SIDIS cross-sections integrated over $k_\perp$. We find that both elastic and inelastic scattering are controlled by the black disk limit, so they yield similar contributions, of zeroth order in the QCD coupling.


Contents
1 Introduction One of the main objectives of the future experimental program at the Electron-Ion Collider (EIC) [1,2] 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 1 of nucleons in the structure of the nucleus (with A 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 Q 2 /s 1 when s Q 2 . Here, Q 2 is the virtuality of the space-like photon γ * exchanged between the electron and one nucleon from the nuclear target and s is the center-of-mass energy squared for the γ * -nucleon scattering.
One expects this rise to be eventually tamed by gluon saturation, i.e. by non-linear effects associated with gluon self-interactions, which limit their occupation numbers to values of order 1/α 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 Q 2 s (A, x). The saturation momentum Q s (A, x) is the typical transverse momentum of the saturated gluons. It increases with the energy, Q 2 s ∝ x −λs with λ s 0.2, and also with the nucleon number A (Q 2 s ∝ A 1/3 when A 1). Yet, for the values of x that should be accessible at the EIC and for A 200, this scale Q 2 s (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 Q 2 s can be hindered by non-perturbative contamination: the perturbative quark-antiquark (qq) 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 1 ÷ 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 Q 2 s [3][4][5][6]. The "golden probe" which attracted most attention over the last years is the production of a pair of hadrons (or jets) in "dilutedense" collisions (eA or pA) at forward rapidities (in the fragmentation region of the dilute projectile) [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Even when the final particles/jets have relatively large transverse momenta (k 2 ⊥ ∼ Q 2 Q 2 s ), 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") [24]. In order to reduce the latter, one must again limit oneself to semi-hard values for Q 2 [25].
In this paper, we make the observation that saturation effects in relatively hard eA collisions with Q 2 Q 2 s could be observed at the EIC in an even simpler process, the semi-inclusive production of a single hadron (or jet), a.k.a. SIDIS, provided the measured particle/jet is very forward -meaning that it carries a large fraction z 1 of the longitudinal momentum of the virtual photon. The precise condition is that Q 2 (1 − z) Q 2 s , which for large virtualities Q 2 Q 2 s forces z to be very close to one: 1 − z 1. When this happens, saturation effects become important because the qq fluctuation which mediates the interactions between the virtual photon and the nuclear target has a large transverse size, r 1/Q s , hence it typically scatters off saturated gluons with transverse momenta p ⊥ ∼ Q s . As implicit in the above discussion, throughout this paper we use the "colour-dipole" picture for DIS at small Bjorken x 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 γ * -nucleon collision -in which the virtual photon has a large longitudinal momentum, but zero transverse momentum (q ⊥ = 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 qq fluctuation of the virtual photon (a "colour dipole") and the small-x gluons from the target. The transverse size of this dipole 1 can be estimated as r ∼ 1/Q, whereQ 2 ≡ z(1 − z)Q 2 , with z and 1 − z the longitudinal fractions of the two quarks.
So long asQ 2 Q 2 s , 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 crosssection is proportional to the gluon distribution. However, if z(1 − z) is sufficiently small, one can haveQ 2 Q 2 s even at high virtualities Q 2 Q 2 s , 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 ⊥ of the produced hadron, the polarisation of the virtual photon, and the type of process under consideration -elastic 2 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 [26][27][28], 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 [29][30][31][32][33][34][35], or their BK truncation at large N c [29,36]). The leading-order predictions of this theory for both single particle (quark) and two particle (quark-antiquark) production 3 are well known [4,5,37] 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 effective virtualityQ 2 ≡ z(1 − z)Q 2 . We shall immediately specialise to the interesting regime atQ 2 Q 2 s , 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 [38,39] 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 caseQ 2 Q 2 s ; 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 ⊥ -distribution of the produced quark is found to be a Gaussian with dispersion k 2 ⊥ = Q 2 s , which describes transverse momentum broadening via multiple scattering. This 1 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 ⊥ Q . 2 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-Nc approximation that we shall generally use. 3 In applications to phenomenology, such partonic cross-sections should be convoluted with parton-to-hadron fragmentation functions.
Gaussian is surrounded by power-low distributions at both lower momenta and higher momenta. At lower momenta k 2 ⊥ Q 2 s , we find a 1/k 2 ⊥ -distribution which is in fact the distribution generated by the decay γ * T → qq of the transverse virtual photon and which is preserved by the elastic scattering in the black disk limit. At larger momenta k 2 ⊥ Q 2 s , we find a power tail ∝ 1/k 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 [40,41] and theoretically studied in several publications [42][43][44][45][46][47]: that is, it shows nuclear suppression at low momenta, k ⊥ < Q s , and nuclear enhancement at larger momenta, k ⊥ Q s , with the possible appearance of a (wide) Cronin peak near k ⊥ = 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 ∼ 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 [48,49] (see also the related works [50][51][52][53][54][55][56]). 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 ⊥ 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 [40,41,[44][45][46].
Since the identification of very forward jets, or hadrons, with semi-hard k ⊥ ∼ Q s may be experimentally challenging, we also provide estimates for the SIDIS cross-sections integrated over the transverse momentum k ⊥ -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 atQ 2 Q 2 s . 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 ⊥ -dependent physics which is quite different in the two cases: the γ * T -wavefunction in the case of the elastic scattering and, respectively, the k ⊥ -broadening of the measured quark for the inelastic one. We also find it interesting to study the weak scattering regime at large effective virtualityQ 2 Q 2 s , where saturation is still visible, albeit rather indirectly, via the phenomenon of geometric scaling [57][58][59].
A striking property of the integrated cross-section is its rapid increase with z when z is large enough (say, z ≥ 0.9). Such a rise should be easily seen experimentally, although it could be more difficult to distinguish between the power-law increase expected in the dilute regime atQ 2 Q 2 s and, respectively, the logarithmic increase in the black-disk regime atQ 2 Q 2 s . Once again, this distinction can become sharper by constructing appropriate R pA ratios, that we shall introduce and study in Sect. 3.
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 1 1 − z) dipole configurations of interest for us here represent only a small fraction of the total DIS cross-section when Q 2 Q 2 s . Of course, asymmetric configurations per se are not that rare: the 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 qq pairs with z(1 − z) 1). Yet, the configurations that we consider in this paper are even more asymmetric: they must obey z(1 − z) Q 2 s /Q 2 1 (in order to be sensitive to gluon saturation), unlike the traditional aligned-jet configurations, for which Q 2 s /Q 2 z(1 − z) 1 (corresponding to weak scattering and dilute parton distributions; see the discussion in Sect. 5).
But albeit very asymmetric and rare when Q 2 Q 2 s , the configurations relevant for saturation become more and more symmetric, and also typical, when decreasing Q 2 towards Q 2 s . 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 Q 2 s 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 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 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 ⊥ ∼ Q s and at very forward rapidities.
The paper is organised as follows. Sect. 2 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. 3 and 4 contain our main new results. In Sect. 3 we consider the double-differential (in z and k ⊥ ) SIDIS cross-sections in the saturation regime atQ 2 = z(1 − z)Q 2 < Q 2 s . We present both analytic results (valid in the strictly-ordered limitQ 2 Q 2 s and for different ranges in k ⊥ ) and numerical results, which include the effects of the high-energy evolution. In Sect. 4, we discuss the SIDIS cross-section differential in z but integrated over k ⊥ . Once again, we combine analytic approximations (in the limiting casesQ 2 Q 2 s andQ 2 Q 2 s ) and numerical results (for the whole kinematics). In Sect. 5 we discuss the inclusive DIS cross-section (as obtained by integrating the SIDIS cross-section over both z and k ⊥ ), 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 Q 2 s . We notably explain the difference between these configurations and the better known aligned jets. Sect. 6 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.
2 DIS in the dipole picture: from dijets to SIDIS 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 dipole frame, by which one means a Lorentz frame in which the virtual photon (γ * ) has a large longitudinal momentum, q 2 z Q 2 , and zero transverse momentum, q ⊥ = 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 4 . That said, both for the purposes of the physical discussion and for the description of the high energy evolution 5 , it is convenient to work in a special frame where both γ * 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 (qq) 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 whereas the nucleus is a left-mover with 6 4-momentum P µ = δ µ− P − per nucleon. In the high-energy or small Bjorken x regime, with the coherence time ∆x + 2q + /Q 2 of the virtual photon, i.e. the typical lifetime of its qq fluctuation, is much larger than the longitudinal extent ∼ A 1/3 /P − of the nuclear target. This justifies the use of the eikonal approximation when computing the dipole-target scattering.

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 qq pair. This depends upon the polarisation (transverse or longitudinal) state of the virtual photon, with the following results [3,5,6]: and respectively: In writing these expressions, we used longitudinal momentum conservation, q + = k + 1 + k + 2 , to write k + 1 = zq + and k In applications to the phenomenology, these results must of course be translated to the laboratory frame at the EIC. 5 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 [48,49] (see also Appendix D for a summary). 6 We neglect the proton mass M which is much smaller than all the other scales in the problem, M 2 Q 2 2P · q.
r ≡ x 1 − x 2 are the transverse sizes of the dipole in the direct amplitude (DA) and in the complexconjugate amplitude (CCA), respectively. The modified Bessel functions K 0 (Qr) and K 1 (Qr) belong to the light-cone wavefunctions describing the qq fluctuation of the virtual photon for the two possible polarisations. These functions exponentially vanish whenQr 1, showing that the transverse size of the dipole cannot be much larger than 1/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 prod i⊥ ∼Q. Of course, these production momenta can be subsequently modified by the scattering, so their measured values k 1⊥ and k 2⊥ can be very different fromQ (typically, larger than it).
The QCD scattering between the dipole and the gluons from the nuclear target is encoded in the following S-matrix structure, where the two-point (dipole) and four-point (quadrupole) functions are defined as and, respectively, Here, V (x) and V † (x) are Wilson lines in the fundamental representation, which describe multiple scattering for the quark and the anti-quark, in the eikonal approximation; explicitly, 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 qq pair) be on mass-shell: the total light-cone energy in the final state, that is, 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 g P − transferred by the gluons from the target, via the collision. This condition yields Note that this effective value x g can be significantly larger than the Bjorken x variable (2.1) if the final transverse momenta k i⊥ are considerably larger thanQ. For what follows, it is useful to keep in mind that the typical transverse momentum that can be transferred by the collision is of the order Q s ≡ Q s (A, x g ) -the nuclear saturation momentum at the longitudinal scale x g . So, we typically expect k 2 i⊥ ∼ max Q 2 , Q 2 s . The total transverse momentum ∆ ≡ k 1 + k 2 of the produced pair is generated via multiple scattering, hence its magnitude is of order Q s . 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 (2.2)-(2.3). Within the time-ordered LC perturbation theory, there are two contributing amplitudes, illustrated in Fig. 1: the virtual photon can split into the qq 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 γ * → qq in the same colour state. For the process in Fig. 1. 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 α and β, with transition amplitudes measured by the matrix elements of the Wilson lines. For the other process, cf. Fig. 1. right, one clearly has α = β. This discussion motivates the following expression for the total amplitude [3], where Ψ(r, z; q + , Q 2 ) generically denotes the qq component of the γ * wavefunction, for either transverse or longitudinal polarisations (see Eqs. (19)- (20) in [5] for explicit expressions). The negative sign in front of the unit matrix in Eq. (2.9) 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 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. (2.4). This discussion also shows that the two dipole S-matrices which enter Eq. (2.4) with negative signs are interference terms between the two processes shown in Fig. 1. It will be also interesting to consider dijet production via elastic (or "coherent") scattering. The "elastic" version of Eq. (2.9) is obtained by demanding the qq 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 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 2.2 Single inclusive jet production in DIS at small x 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. (2.2)-(2.3). The integral over k + 2 simply removes the δ-functions for longitudinal momentum conservation. That over k 2 identifies the transverse coordinates of the antiquark in the DA and in the CCA: x 2 = x 2 . In turn, this introduces an important simplification in the colour structure: the original quadrupole in Eq. (2.6) reduces to an effective dipole, S (2) xg (x 1 , 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 πR 2 A . Assembling everything together, we find the following expressions for the longitudinal and transverse SIDIS cross sections: (2.14) We have introduced here the reduced cross-sections, which encode all the non-trivial dynamics. As suggested by their notations, J L and J T depend upon the virtuality of the initial photon Q 2 and upon the splitting fraction z for the γ * → qq decay only via the productQ 2 = z(1 − z)Q 2 (the "effective virtuality").
In writing Eqs. (2.15)-(2.16), we found it useful to replace the dipole S-matrix by the corresponding scattering amplitude, T (r, x g ) ≡ 1 − S (2) xg (x 1 , x 2 ), with r = x 1 − 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. (2.11)-(2.12). For the transverse photon, one has 18) 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 Eq. (2.8), one sees that, in the case of two-particle production, x g depends upon the transverse momenta k 1 and k 2 of both final partons. So, our above treatment of the integral over k 2 , which ignored this additional dependence, may look as just an approximation. This is indeed true in general, but not also for the elastic scattering: one can easily check that, for the case of a homogeneous target, the r.h.s of Eq. (2.17) is proportional to the deltafunction δ (2) (k 1 + 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 k 2 simply identifies k 2 = −k 1 , and the version of Eq. (2.8) for SIDIS in the elastic channel reads without any approximation. As for the inelastic channel, it is easy to identify what should be the right approximation to Eq. (2.8): the integral over k 2 is controlled by the typical transverse momenta of the unmeasured parton, which as discussed after Eq. (2.8) are of the order of the largest amongQ 2 and Q 2 s . So in practice we shall use The integrations over r and r in Eqs. (2.15)-(2.18) are controlled by a product of 3 factors: the modified Bessel functions K 0 (Qr) and K 1 (Qr) describing the wavefunction of the virtual photon, the amplitude T (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 virtualityQ 2 = z(1 − z)Q 2 , the target saturation momentum Q s (A, x g ), and, respectively, the transverse momentum k ⊥ ≡ |k| of the produced quark. As already mentioned, the Bessel functions effectively restrict the integrations to values r 1/Q (larger values are exponentially suppressed). Also, the Fourier phase is rapidly oscillating for dipole sizes r 1/k ⊥ . Finally, the dipole amplitude is small, but rapidly increasing with r, for small dipoles, T (x g , r) ∼ r 2 Q 2 s 1 when r 1/Q s ("colour transparency"), but this growth is tamed by multiple scattering when r 1/Q s . For larger dipole sizes, r 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 ∼ 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 [60].) Hence, we would like to choose the external scales k ⊥ , z, and Q 2 in such a way that the integrations favour dipole sizes r, r 1/Q s . For this to be true, we clearly need the two following conditions: (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 rQ s 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 2 ⊥ Q 2 s . 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.

Saturation effects in SIDIS at
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 Q 2 s , for the perturbative treatment of the γ * decay to be justified, while at the same time z is sufficiently close to 1 forQ 2 = z(1 − z)Q 2 to be (much) smaller than Q 2 s -meaning that the dipole scattering is sensitive to gluon saturation in the nuclear target. In this section we shall analytically study the strongly ordered caseQ 2 Q 2 s , 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:Q 2 Q 2 s . In the subsequent analysis, it will be convenient to distinguish between three regimes for the transverse momentum k ⊥ of the produced quark: probing the virtual photon wavefunction When bothQ 2 and k 2 ⊥ are much smaller than Q 2 s , 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 Eq. (2.17) for J T,el are dominated by relatively large dipole sizes, which obey 1/Q s r, r 1/Q, for which the scattering is as strong as possible: T 1. (The respective contributions from smaller dipoles are suppressed by colour transparency: T (r) ∝ r 2 for r 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 Eq. (2.18) mutually cancel in their sum. Accordingly, the net contribution to J T,in comes from smaller dipoles r, r 1/Q s and thus is parametrically suppressed compared to J T,el . Clearly, the presence of the last term in the r.h.s. of Eq. (2.18), which is quadratic in T and becomes important at strong scattering, was essential for the above argument.

Elastic scattering near the unitarity limit
So let us focus on the elastic cross-section in this section. As visible in Eq. (2.17), 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 r. As just mentioned, the dominant contribution whenQ 2 , k 2 ⊥ Q 2 s can be computed by replacing T → 1 in Eq. (2.17). After this replacement, the elastic amplitude reduces to the Fourier transform (FT) of the γ * wavefunction (in particular, it becomes independent of the QCD coupling α s ). By also using the integrals (A.1) and (A.7) from Appendix A, one finds and 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 ⊥ ∼Q (in particular, it vanishes in the photo-production limitQ 2 → 0), whereas for a transverse photon, it is logarithmically distributed within the rangeQ 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 1/Q s . In terms of Feynman graphs, the contributions shown in Eqs. (3.1)-(3.2) fully correspond to the second graph in Fig. 1 -that is, the decay of γ * occurs after crossing the nuclear target. The contribution of the first graph to the elastic scattering is proportional to the dipole Smatrix (recall Eq. (2.10)), 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 atQ 2 k 2 ⊥ Q 2 s , the cross-section for SIDIS is controlled by the elastic scattering of the qq fluctuation of a virtual photon with transverse polarisation and can be approximated as Although not manifest in this result for J 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 2 s ≡ Q 2 s (A, x g ) on k 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 Eq. (2.19), which for the present kinematics reduces to 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 ∼ 10 −2 , where the MV model applies.

Elastic scattering in the MV model
The McLerran-Venugopalan (MV) model [38,39] 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 Eq. (2.7)), 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. [26,61]) where the quantity in the exponent is the dipole amplitude in the single scattering approximation (and in the MV model). This quantity Q 2 A , 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 α s like Q 2 A ∝ α 2 s A 1/3 . (One factor of α s is inherent in the colour charge density, the other one comes from the coupling to the dipole.) Furthermore, Λ is the QCD confinement scale and the logarithm has been generated by integrating over Coulomb exchanges with transverse momenta q ⊥ within the range Λ 2 q 2 ⊥ 1/r 2 . Eq. (3.5) is derived under the assumption that both Q 2 A and 1/r 2 are much larger than Λ 2 and is valid to leading logarithmic accuracy w.r.t. to the Coulomb log.
Our present interest in Eq. (3.5) 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, Eq. (3.5) 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 J T,el , which extends the previous result in Eq. (3.3) to larger transverse momenta, k ⊥ ∼ Q s . The corresponding approximation for the inelastic cross-section J T,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 8 : this is the value Q s (A) = 2/r for which the exponent in Eq. (3.5) becomes equal to one; this condition yields We now return to Eq. (2.17) for J T,el and construct an analytic approximation valid within the MV model and within the extended rangeQ 2 k 2 ⊥ Q 2 s . Since the interesting dipole sizes obey r 1/Q, one can keep only the leading order term, K 1 (Qr) 1/(Qr), in the small-argument expansion of the Bessel function; then Eq. (2.17) reduces to 9 Since independent ofQ, this result is formally the same as the limit of Eq. (2.17) whenQ → 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 Eq. (3.5), enables us to almost exactly compute the integral in Eq. (3.7). 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 → 2/Q s within the slowly varying logarithmic factor in the exponent of Eq. (3.5).
That is, one can effectively replace the dipole S-matrix by the Gaussian S 0 (r) exp(−r 2 Q 2 s /4), where we have also used the definition (3.6) for Q 2 s . The ensuing integral over r can be exactly evaluated, as explained in Appendix B; one thus finds This result looks indeed as a natural extension of our previous estimate in Eq. (3.3) to larger transverse momenta k ⊥ ∼ Q s . It has a natural physical interpretation. The power 1/k 2 ⊥ is the transverse momentum distribution produced by the decay of the transverse virtual photon, as already discussed. The Gaussian exp(−2k 2 ⊥ /Q 2 s ) describes the subsequent broadening of this distribution via elastic multiple scattering between the qq dipole and the gluons from the nucleus. Eq. (3.8) shows that the deviations from the pure power law should become important for k 2 ⊥ ∼ Q 2 s /2, where the factor 1/2 reflects the elastic nature of the scattering process.

Numerical results for
As promised, we shall now compare the analytic approximations that we previously developed for the kinematical range atQ 2 Q 2 s and k 2 ⊥ Q 2 s with numerical calculations. We start with the MV model, for which we shall use a slightly modified version of Eq. (3.5), which is well defined for any value of r, including r > 1/Λ, and thus allows for unrestricted numerical integrations; this reads (below, e = 2.718... is the Euler number) 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 (3.6), but the difference is irrelevant for our purposes. In practice, we have chosen Q 2 A and Λ 2 such as Q 2 s0 = 2 GeV 2 for the MV model S-matrix in Eq.
JT,el JT,inel respectively. The (reduced) cross-sections multiplied by k 2 ⊥ are shown as functions of k ⊥ , for k ⊥ ≤ Q s0 and several valuesQ 2 Q 2 s0 . Here, Q 2 s0 = 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 ⊥ , shows a maximum around k 2 ⊥ ∼QQ s0 , and then decreases. Inelastic scattering keeps increasing until k ⊥ ∼ Q s0 ; it becomes the dominant contribution already below Q s0 (see also the discussion in Sect. 3.2). clarity, in all the plots and the associated discussions we use the specific notation Q 2 s0 for the saturation momentum in the MV model.) Our numerical results corresponding to the MV model are shown in the two plots in Fig. 2. For convenience, we display results for the reduced cross-sections multiplied by k 2 ⊥ . That is, we plot dimensionless quantities like k 2 ⊥ J T,el (k ⊥ ) as functions of k ⊥ , in the saturation region at k ⊥ ≤ Q s0 . We use Q 2 s0 = 2 GeV 2 together with 3 values for the effective virtualityQ 2 , namelyQ 2 = 0.02, 0.1, and 0.2 GeV 2 , which are well separated from each other and also from Q 2 s0 , and thus are useful for indicating theoretical trends. These values forQ 2 are quite low and should be taken with a grain of salt. The largest among them,Q 2 = 0.2 GeV 2 , corresponds to a transverse momentum k ⊥ =Q 450 MeV, which is a reasonable scale for the transition from non-perturbative to perturbative regimes. The even lower values, likeQ 2 = 0.02 GeV 2 , are strictly speaking non-perturbative, but they are only used to illustrate the interesting phenomena, outside any phenomenological context.
Specifically, in Fig. 2a we compare the elastic cross-sections in the transverse and in the longitudinal sector, respectively, whereas in Fig. 2b we compare the elastic and the inelastic cross-sections, both in the transverse sector. The curves in Fig. 2a have the expected shapes in light of the previous discussion in this paper: k 2 ⊥ J L,el is peaked around k ⊥ =Q, whereas k 2 ⊥ J T,el shows a wide maximum ("plateau") at, roughly,Q < k ⊥ < Q s0 / √ 2. When increasingQ, the height of this maximum is decreasing and its position gets displaced towards larger values of k ⊥ . This can be understood as follows: whenQ → 0, Eq. (3.8) shows that the maximum of k 2 ⊥ J T,el lies at k ⊥ = 0 and its height is equal to 1/(4π) 0.08. The effect of increasingQ can be appreciated by using the following approximation, which interpolates between Eq. (3.2) at low k ⊥ ∼Q Q s0 and Eq. (3.8) at larger momenta Q k ⊥ ∼ Q s0 : It is not difficult to check that the function k 2 ⊥ J T,el develops a maximum at k 2 ⊥ ∼QQ s0 and that the height of this maximum is proportional to exp(−2Q/Q s0 ) -in qualitative agreement with Fig. 2a. Fig. 2a also shows that the longitudinal piece J L,el dominates at very low momenta k ⊥ Q , but it is much smaller than the transverse piece J T,el at any k ⊥ Q . Of course, one should be careful when comparing numerical results for J T and J L , as these quantities are weighted by different factors in the respective physical cross-sections, cf. Eqs. (2.13)-(2.14). Notably, the longitudinal cross-section has an additional suppression when z → 1, due to the prefactor z(1 − z) in (2.13). Altogether, we conclude that the longitudinal contribution is tiny and can be safely neglected in the physical regime under consideration. Concerning Fig. 2b, we notice that the elastic cross-section J T,el is indeed larger than the inelastic one J T,in at sufficiently small k ⊥ Q s , as expected, but this hierarchy is changing quite fast when increasing k ⊥ . The reasons for this behaviour should become clear in the next section, where the inelastic contribution will be explicitly computed.
Furthermore, Fig. 4 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 [48,49]. (For convenience, this equation is briefly summarised in Appendix D.) The evolution "time" is the target rapidity 10 η = 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, Eq. (3.9), as an initial condition at η = 0 and evolve the collinearly-improved rcBK equation up to η = 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 (η ≤ 6). These are numerically extracted from the solution to the (improved) BK equation, with the results shown in Fig. 3.
In Fig. 4a we show the elastic cross-section k 2 ⊥ J T,el for η ≤ 3 and the same 3 values for the ratiō Q 2 /Q 2 s0 as in Fig. 2. The results are plotted as functions of k ⊥ , up to k ⊥ = 2 GeV. (For comparison, we note that Q s (η = 3) 1.77 GeV, cf. Fig. 3.) In Fig. 4b, we fixQ 2 /Q 2 s0 = 0.05 and display the functions k 2 ⊥ J T,el and k 2 ⊥ J T,in for the same values of η and in the same range in k ⊥ 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 η 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 ⊥ , as shown in Eq. (3.4) 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. 4, and the logarithmic dependence of η upon x g , it should be quite clear that such a k ⊥ -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.

k 2
⊥ ∼ Q 2 s : transverse momentum broadening In the previous section we have seen that, when increasing k 2 ⊥ aboveQ 2 , the elastic contribution to SIDIS is rapidly decreasing, according to Eq. (3.3). On the other hand, the inelastic contribution is expected to be most important when k ⊥ ∼ 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. 2 and 4, which show that k 2 ⊥ J T,in keeps increasing when k ⊥ increases towards Q s . In this section, we shall construct an analytic approximation for J T,in for transverse momenta in the rangeQ 2 k 2 ⊥ ∼ Q 2 s (by which we mean that k ⊥ can be both smaller, or larger, than Q s , but not much larger; the tail of the cross-section at k ⊥ Q s will be studied in Sect. 3.3). This will allow us to physically understand the respective numerical results in Figs. 2 and 4, as well as their extension to larger values k ⊥ Q s .

Inelastic scattering off the saturated gluons
To study the inelastic cross-section (2.18) for k ⊥ ∼ Q s , it is convenient to return to its representation in terms of dipole S-matrices (recall Eq. (2.12)), rather than scattering amplitudes. One has where we recall thatQ 2 k 2 ⊥ ∼ Q 2 s . In this regime, the double integral is dominated by dipole sizes within the range 1/Q s r, r 1/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 → 0 (or T → 1), since the net result would vanish in that limit. Yet, when S 1, the second term within the square brackets in Eq. (3.11), 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 e −ik·(r−r ) S(r − r ) restricts the double integration to values r and r such that |r − r | 1/k ⊥ ∼ 1/Q s 1/Q. This makes it convenient to change one integration variable, say r → ρ ≡ r − r . Then, after also using the small-argument expansion of the Bessel function, K 1 (z) 1/z for z 1, one sees that the integral over r at fixed ρ shows a logarithmic enhancement over the range ρ r 1/Q, with ρ = |ρ|. Indeed, one can write r·r r r K 1 (Qr)K 1 (Qr ) 1 Q 2 r 2 for 1/Q r ρ , (3.12) and the ensuing integral over r yields the logarithm ln(1/ρ 2Q2 ), which is large when 1/ρ ∼ Q s . To leading logarithmic accuracy with respect to this transverse log, one can further write where the second line follows from the fact that the integral over ρ is controlled by values ρ ∼ 1/Q s , as we shall shortly check; hence, to the leading-log accuracy of interest, it is legitimate to replace ρ → 1/Q s within the argument of the logarithm from the first line.
One can verify that the inelastic contribution J L,in of the longitudinal photon has no logarithmic enhancement 11 and hence it is parametrically smaller than the transverse contribution in Eq. (3.15).
The final result in Eq. (3.13) has a simple physical interpretation. The logarithm ln(Q 2 s /Q 2 ) comes from integrating over the initial transverse momentum distribution of the quark, as generated by the decay of the virtual photon γ * T and shown in Eq. (3.3). The second factor in the r.h.s. of Eq. (3.13) is recognised as the Fourier transform of the dipole S-matrix. It describes the final transverse momentum distribution of the quark, as generated via (multiple) scattering off the nuclear target.
The dipole S-matrix appearing in Eq. (3.3) should be evaluated at a rapidity η = ln(x 0 /x g ) with x g given by Eq. (2.20); within the present kinematics, that is, k 2 ⊥ ∼ Q 2 s Q 2 , this yields with x as defined in Eq. (2.1). In the regime of interest here, the rapidity η can be considerably smaller than the value ln(x 0 /x) that would be normally used for studies of inclusive DIS.

Inelastic scattering in the MV model
To obtain a more explicit expression for J T,in , let us evaluate Eq. (3.13) within the MV model, that is, by using the dipole S-matrix in Eq. (3.5). This calculation can be simplified via the same argument as discussed in relation with Eq. (3.7): in the interesting range at k ⊥ Q s , the FT in Eq. (3.13) is dominated by the largest value of ρ that is allowed by the dipole S-matrix, namely ρ 2/Q s . To the accuracy of interest, one can replace the S-matrix in Eq. (3.5) by the Gaussian S 0 (r) exp{−r 2 Q 2 s /4}, with Q 2 s defined in Eq. (3.6). The one easily finds This Gaussian distribution in k ⊥ is the hallmark of a diffusive process: it describes transverse momentum broadening via soft (Λ 2 q 2 ⊥ Q 2 s ) 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 (J T,el ) and the inelastic (J T,in ) cross-sections, and for transverse momenta within the rangē Q 2 k 2 ⊥ ∼ Q 2 s . By comparing the respective results in Eqs. (3.8) and (3.15), one finds that the elastic scattering dominates over the inelastic one so long as k ⊥ is small enough for Q 2 s /k 2 ⊥ > ln(Q 2 s /Q 2 ). Using Q 2 s Q 2 , it is not difficult to check that the "critical" value k 2 ⊥ Q 2 s / ln(Q 2 s /Q 2 ) at which the elastic and the inelastic contributions are roughly equal to each other lies in between the respective maxima 12 , namely k 2 ⊥ Q sQ for the elastic piece and k 2 ⊥ Q 2 s for the inelastic one. This is in agreement with the numerical results shown in Figs. 2b and 4b. More results will be presented in the next subsection.

Numerical results for
The two plots shown in Fig. 5 extend the previous ones in Fig. 2b and, respectively, Fig. 4b to larger values for the transverse momenta, namely up to k 2 ⊥ = 3Q 2 s . Notice also that, as compared to Figs. 2 and 4, the rescaled cross-sections k 2 ⊥ J T,in and k 2 ⊥ J T,el are now plotted as functions of k 2 ⊥ , and not of k ⊥ . Fig. 5a shows the respective predictions of the MV model for Q 2 s0 = 2 GeV 2 and 3 values different for the effective virtualityQ 2 , namelyQ 2 /Q 2 s0 = 0.01, 0.05, 0.1. Fig. 5a illustrates the effects of the high-energy evolution up to η = 3, for k 2 ⊥ 3Q 2 s (η = 3) and a fixed ratioQ 2 /Q 2 s0 = 0.05. We have already discussed the results for the elastic cross-section in some detail, in relation with Figs. 2 and 4. So, let us now focus on the inelastic contributions. The respective MV-model curves in Fig. 5a show the same qualitative behaviour as the analytic estimate in Eq. (3.15) (multiplied by k 2 ⊥ ): a rapid increase at low k 2 ⊥ Q 2 s , followed by a wide peak around Q 2 s , and finally a slow decrease at larger values k 2 ⊥ > Q 2 s . Also, when increasingQ 2 , the position of the peak remains unchanged -as it should, since this is located at k 2 ⊥ = Q 2 s -but its height is decreasing, due to the reduction in the logarithmic factor ln(Q 2 s /Q 2 ) in Eq. (3.15). 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 withQ, as discussed in relation with Fig. 2 (and also visible in Fig. 5a).
This general trend is preserved by the BK evolution up to η = 3, as shown in Fig. 5b. Due to the increase in the saturation momentum with η, cf. Fig. 3, the position of the peak is (slightly) moving towards larger values of k 2 ⊥ , for both the elastic and the inelastic cross-sections. Also, the magnitude of the cross-sections is growing with η, albeit quite slowly and only at sufficiently large transverse momenta k 2 ⊥ Q 2 s (η). This growth is more important for the inelastic channel, where it 11 Using K0(z) ln(1/z) for z 1, one can check that, in the kinematical regime at study, the integrals in (2.18) are dominated by very large dipoles, with r, r ∼ 1/Q, but such that ρ = |r − r | 1/k ⊥ . Accordingly, there is no logarithmic phase-space for the integral over r at fixed ρ. 12 More precisely, these maxima refer to the reduced cross-sections multiplied by k 2 ⊥ , as in Figs. 2 and 4. can be attributed to the corresponding increase in the logarithm ln(Q 2 s (η)/Q 2 ), cf. Eq. (3.13).

k 2
⊥ Q 2 s : hard scattering and BK evolution When the transverse momentum k ⊥ of the produced quark is much larger than the target saturation momentum, k 2 ⊥ Q 2 s , the physical picture of SIDIS is quite different: the (elastic and inelastic) crosssections are controlled by relatively small dipoles, with transverse sizes r, r 1/Q s , which scatter only weakly off the hadronic target: T (r) 1. An immediate consequence of that is that the elastic cross-section, which is quadratic in T (cf. Eq. (2.17)), is even more strongly suppressed relative to the inelastic one, which in turn is dominated by the first 3 terms, linear in T , in Eq. (2.18). That said, gluon saturation may still influence SIDIS even for such high transverse momenta, via its consequences (like geometric scaling [57][58][59]) for the high-energy evolution of the dipole scattering amplitude.
In this section, we shall consider the situation where k ⊥ is much larger than both Q s andQ, without any strong ordering between the two latter scales (this ordering becomes irrelevant when k 2 ⊥ Q 2 s ,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.

Hard scattering in the MV model
To compute the total (or inelastic) cross-section for k 2 ⊥ Q 2 s ,Q 2 within the MV model, one can use the single scattering approximation to the dipole S-matrix (3.5), that is, together with the small-argument expansion of the Bessel function: K 1 (z) 1/z for z 1. Under these assumptions, the dominant contribution to J T (the only one to be enhanced by a large transverse logarithm) is that generated by the last term in Eq. (2.16), i.e. T (r−r ), which refers to the scattering of the produced quark. (This term enters Eq. (2.16) 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/Q, since it is only their difference ρ = r −r which is constrained by the Fourier phase to small values ρ 1/k ⊥ . Changing integration variables from r, r to r, ρ, the integral over r is logarithmic within the interval ρ r 1/Q (recall Eq. (3.12)), and one finds The expression in the first line is recognised as the single-scattering approximation to the first line of Eq. (3.13). The obtention of the final result in the second line is a bit subtle and requires the more elaborated calculations presented in Appendix C. Here, we shall rather focus on its physical origin. Notice first that, for the kinematics at hand, the integral over ρ is controlled by the Fourier phase, which selects ρ ∼ 1/k ⊥ . One may naively think that the dominant contribution can be obtained by replacing ρ 2 → 1/k 2 ⊥ within both logarithms inside the integrand: such a contribution would be quadratic in the transverse logarithms and hence parametrically larger than the linear combination shown in the second line of Eq. (3.13). However, if one does so, then the ensuing double logarithm gets multiplied by an integral over ρ which vanishes when k ⊥ > 0 (since its result is proportional to ∇ 2 k δ (2) (k)). So, it is natural to find a result which is only linear in the transverse logs, and also symmetric under the exchangeQ 2 ↔ Λ 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 Eq. (3.13), the ln(k 2 ⊥ /Q 2 ) comes from integrating over the initial k ⊥ -distribution of the quark produced by the decay of γ * T , Eq. (3.3). As for the associated power tail ∝ 1/k 4 ⊥ , this is generated by a hard single scattering between the quark and a nucleon from the target, with a transferred momentum q ⊥ k ⊥ .
Concerning the first term in Eq. (3.17), it is quite clear that the logarithm ln(k 2 ⊥ /Λ 2 ) is directly inherited from the dipole amplitude (3.16): this is a Coulomb logarithm, generated via a relatively soft exchange, within the range Λ 2 q 2 ⊥ k 2 ⊥ . Accordingly, the power-like tail ∝ 1/k 4 ⊥ for this term should somehow be associated with the decay of the virtual photon. This may look surprising since, as discussed in Sect. 3.1, the γ * T decay is expected to produce a much softer tail ∝ 1/k 2 ⊥ , cf. Eq. (3.3). Yet, as shown by the calculation in Appendix C, a 1/k 4 ⊥ tail is in fact generated by the difference between the two such distributions, with and without soft rescattering in the final state.
In principle, Eq. (3.17) should be used in the case whereQ 2 Λ 2 , to justify the use of a perturbative (partonic) description for the γ * wavefunction; in such a case, the first logarithm in the final result would dominate over the second one. In practice though, the difference betweenQ 2 and Λ 2 may not be that important, which explains why we have kept both logarithms in Eq. (3.17).
For completeness, let us also observe that the elastic cross-section J T,el in the MV model and in the single scattering approximation can be easily estimated from Eqs. (2.17) and (3.16). One thus finds J T,el ∼ Q 4 A /k 6 ⊥ (see Appendix B for an explicit calculation and Eq. (B.10) for the final result). This is suppressed by a power of Q 2 A /k 2 ⊥ 1 compared to the inelastic piece (a "higher twist effect"), as expected.

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 rQ s 1 is known as [59] T (r) (r 2 Q 2 s ) γ ln 1 r 2 Q 2 s for rQ s 1, (3.18) where the "anomalous dimension" γ 0.63 is significantly smaller than 1. It is here understood that Q s ≡ Q s (A, η) depends upon the rapidity η = ln(x 0 /x g ), assumed to be large:ᾱ s η 1. Eq. (3.18) shows geometric scaling [57][58][59]: it depends upon the 2 independent kinematical variables r and η (and also upon the nucleon number A) only via the product r 2 Q 2 s (A, η). This property is an imprint of the physics of saturation on the behaviour of the amplitude in the weak scattering regime at rQ s 1. Clearly, the asymptotic behaviour at high energies is not precisely the regime to be explored via DIS at EIC, where only moderate values η 3 should be accessible. Here, we shall rather use Eq.
However, since the power γ is strictly smaller than one, the mathematical treatment of Eq. (3.19) is considerably simpler. Namely, it is now legitimate to replace ρ 2 → 1/k 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 (E.1) in Appendix E). One thus finds As anticipated, this is strongly suppressed as compared to the inelastic part in Eq. (3.20), due to the larger power of Q 2 s /k 2 ⊥ : 2γ instead of γ. One interesting conclusion which emerges from the calculations in Sects. 3.2 and 3.3 is that, for all transverse momenta k ⊥ Q s and for any energy (from η 0, where the MV model applies, up the very large η 1/ᾱ s , where one can use the asymptotic solution to the BK equation), the dominant contribution to SIDIS at low virtuality (Q 2 Q 2 s ) is generated by the scattering of the measured quark alone (as opposed to the scattering of the quark-antiquark dipole). Of course, the scattering of the whole dipole does matter in general -this is represented by the terms T (r) and T (r ) in Eq. (2.16) -, but its contribution turns out to be (parametrically) less important than that from the measured quark alone -the piece involving T (r−r ) in Eq. (2.16).

Numerical results for k
In Fig. 6 we show numerical results for the quantity k 2 ⊥ J T,in at different rapidities η ≤ 6 in terms of the dimensionless variable k 2 ⊥ /Q 2 s (η), up to relatively large values k 2 ⊥ Q 2 s (η). In the left plot, we use a linear scale and show all the transverse momenta up to k 2 ⊥ /Q 2 s (η) = 6, in order to exhibit the global shape of the function k 2 ⊥ J T,in . In the right plot, we use a logarithmic scale and show only the high momentum regime at 1 < k 2 ⊥ /Q 2 s (η) < 40. Note that after multiplication by k 2 ⊥ , the inelastic (a)  cross-section J T,in is very close to a scaling function (see e.g. Eqs. (3.15), (3.17) and (3.20)) -that is, a function of the ratio k 2 ⊥ /Q 2 s (η). This scaling behaviour is however not exact (not even in our analytic approximations): it is broken by logarithmic factors. For instance, Eq. (3.13), which we recall is valid when k ⊥ ∼ Q s , involves an overall factor ln(Q 2 s (η)/Q 2 ), which increases with η, due to the respective behaviour of the saturation momentum. Similarly, at larger momenta k 2 ⊥ Q 2 s , one can rely on Eq. (3.20) to deduce that the function k 2 ⊥ J T,in can be written as the sum of a scaling function and a scaling-violating term proportional to ln(Q 2 s (η)/Q 2 ). (To that aim, one should rewrite one of the logarithms in Eq. (3.20) as ln(k 2 ⊥ /Q 2 ) = ln(k 2 ⊥ /Q 2 s ) + ln(Q 2 s /Q 2 ).) These analytic arguments suggest that, when plotted in terms of the scaling variable k 2 ⊥ /Q 2 s (η), the function k 2 ⊥ J T,in should still increase with η, albeit only slightly. This is in agreement with the numerical findings in Fig. 6a.
Another interesting feature of the numerical results, better seen in Fig. 6b, is the emergence of an anomalous dimension γ < 1 as a consequence of the high-energy evolution: indeed, the high-k ⊥ tails corresponding to η 4 are clearly seen to be softer (i.e. to show a slower decrease with increasing k 2 ⊥ ) then the MV model curve, for which γ = 1. That said, the rapidities η ≤ 6 considered here are still too small to feature the asymptotic expectation γ 0.63. Besides, the numerical extraction of this power γ is further hindered by the additional logarithmic dependences, as visible in Eq. (3.20).

Cronin peak
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 ⊥ ) 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 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 ⊥ ∼ Q s (A).
For an analytic estimate of Eq. (3.22) at relatively low energies (say, x g ∼ 0.01), we can use the MV model for both types of targets, with Q 2 A = A 1/3 Q 2 p . Specifically, for k ⊥ ∼ Q s (A), the saturation effects are important for the nuclear target, for which one should use the full MV result in Eq. (3.15), but they are not important for the proton target, for which one should rather use the high-k ⊥ approximation in Eq. (3.17). To make the argument clearer, let us focus on the physically interesting regime atQ 2 Λ 2 ; in that case, one can neglect the second logarithm in the final result in Eq. (3.17), which allows us to write where the last estimate holds for k ⊥ ∼ Q s (A). (We have also used Q 2 A = A 1/3 Q 2 p together with Eq. (3.6) for the nuclear saturation momentum.) The first equality in Eq. (3.23) also shows that R pA (k ⊥ ) increases with k ⊥ so long k ⊥ < Q s (A), but it decreases at larger momenta k ⊥ Q s (A). Hence, Eq. (3.23) shows a peak at k ⊥ Q s (A) and the height of this peak is expected to be larger than one in the limit where Q 2 s (A) Q 2 . This is similar to the Cronin peak observed in hadron production at midrapidities in d+Au collisions at RHIC [40,41]. In that context too, this enhancement can be explained by gluon saturation in the nuclear wavefunction and its consequences in terms of multiple  scattering [42][43][44][45][46][47].
The existence of such a peak at the level of the MV model is numerically studied in Fig. 7a. This figure shows the R pA ratio in Eq. (3.22) as a function of k 2 ⊥ (up to 5Q 2 s = 10 GeV 2 ) and for 3 values, Q 2 /Q 2 s = 0.01, 0.05, 0.1, for the ratio between the relevant scales. A mild peak, centered at k 2 ⊥ slightly above Q 2 s , is indeed seen in the caseQ 2 /Q 2 s = 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 ⊥ ) increases quite fast at low momenta until it reaches a value of order 1 for k ⊥ ∼ Q s . For larger k ⊥ > Q s , it shows either a mild peak, or a plateau with height R pA > 1. The rapid increase at low k ⊥ Q s is a consequence of gluon saturation in the nucleus, as quite clear by inspection of the calculation in Eq. (3.23).
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. 7b shows our respective predictions for the most favourable casē Q 2 /Q 2 s = 0.01, as obtained from numerical solutions to the BK equation up to η = 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 2 ⊥ ≤ 10 GeV 2 for η ≥ 2. A similar behaviour was experimentally observed in d+Au collisions at RHIC [40,41] and theoretically explained (within the CGC framework) in Refs. [44][45][46]. Clearly, the same physical mechanism is working in the present set-up: for the values of k ⊥ considered in Fig. 7b, 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 η. This difference transmits to the respective dipole amplitudes and hence to the two cross-sections in the numerator and the denominator of Eq. (3.22), respectively.
So far in this section, we have focused on the well-separated regime atQ 2 Q 2 s , where the saturation effects are more striking and which allowed for analytic approximations. In view of the phenomenology, it is however important to notice that the saturation effects should remain important when the two scales are comparable to each other,Q 2 ∼ Q 2 s , since the corresponding dipole sizes are large (r ∼ 1/Q s when k 2 ⊥ Q 2 s ), hence the unitarity corrections are of O(1). To illustrate this point, we have repeated some of our calculations for larger values ofQ 2 , namelyQ 2 /Q 2 s = 0.25, 0.5 and 1, with the results shown in Fig. 8 (for the MV model, for definiteness). As visible in these plots, the saturation effects -a maximum in the rescaled cross-section k 2 ⊥ J T near k 2 ⊥ = Q 2 s (cf. Fig. 8a) and a modest enhancement (R pA 1) in the R pA ratio at k 2 ⊥ > Q 2 s (cf. Fig. 8b) -are clearly seen forQ 2 = 0.25Q 2 s , but smoothly disappear when further increasingQ 2 . Such a behaviour is indeed consistent with our previous results, notably those appearing in Eqs. (3.15) and (3.23).

k ⊥ -integrated SIDIS and saturation
In the previous section, we have studied the differential cross-section dσ γ * A→qX dz d 2 k for producing a quark with a given longitudinal momentum fraction z and a given transverse momentum k in DIS processes with relatively low effective virtualitiesQ 2 = z(1 − z)Q 2 Q 2 s . Since this double differential quantity may be difficult to measure in the actual experiments, in what follows we shall discuss its more inclusive version dσ/dz 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: s , which is the most interesting regime for a study of saturation, and (ii)Q 2 Q 2 s , where saturation can manifest itself indirectly, via the phenomenon of geometric scaling. But before that, let us briefly list the general formulae valid for any kinematics.

General expressions
Within our current formalism, the fact of integrating out k simply amounts to identifying the transverse coordinates, r and 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. (2.15)-(2.16)) for the elastic piece, for the inelastic one, and, finally, for the total cross-section. We have also used the fact that T (−r) = T * (r), as it can be easily checked by using T = 1−S together with the definition (2.5) of the dipole S-matrix 13 . The T -matrix structure 13 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 [62,63].
in Eq. (4.2) can be equivalently written as where the r.h.s. is recognised as the probability for the dipole with size r to suffer an inelastic scattering (recall that |S(r)| 2 is the respective survival probability). Similarly, Eq. (4.3) involves T (r) + T (−r) = 2Re T (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 Eq. (4.3) 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): with σ dipole (η, r) = 2πR 2 A T (η, r) when T is real (as generally in this paper).

4.2
Low virtuality z(1 − z)Q 2 Q 2 s : 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 atQ 2 Q 2 s . In this regime, all the cross-sections shown in Eqs. (4.1)-(4.3) are dominated by large dipoles such that 1/Q s r 1/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 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 α s , except for the weak (logarithmic) dependence introduced by the scale Q 2 s in the argument of the logarithm (recall that Q 2 s ∼ α 2 s ). To better understand the physical origin of these results, it is instructive to re-derive them by integrating over k the double-differential cross-sections computed in Sect. 3.
Concerning the elastic cross-section, it is quite clear that the respective logarithm in Eq. (4.6) is generated by integrating the k ⊥ -distribution in Eq. (3.3) over the rangeQ 2 k 2 ⊥ Q 2 s . More precisely, the MV estimate in Eq. (3.8) shows that the logarithmic integration over k 2 ⊥ is cut off by multiple scattering already at k 2 ⊥ ∼ Q 2 s /2. Hence the corresponding estimate in Eq. (4.6), 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 2 s /Q 2 ) is already visible in the differential distribution in the second line of Eq. (3.13), or in its MV model estimate (3.15). By integrating any of these two expressions over all k ⊥ , one finds exactly the result displayed in Eq. (4.6). Of course, the Gaussian broadening in Eq. (3.15) only holds within a limited range of momenta around Q s , but this is precisely the range which dominates the k ⊥ -integration to the accuracy of interest 14 . These considerations are corroborated by the numerical results in Fig. 9a, which show the integrated elastic, inelastic, and total cross-sections in terms of ln(Q 2 /Q 2 s ) in the saturation regime at Q 2 Q 2 s . Together with the predictions of the MV model (the curves with η = 0), we also exhibit the results of the BK evolution up to η = 3. ForQ 2 /Q 2 s 0.1, all the curves are roughly straight lines, with slopes which are consistent with Eq. (4.6). 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.

High virtuality z(1 − z)Q 2
Q 2 s : geometric scaling Even though not so directly relevant for studies of saturation (since controlled by weak scattering; see below), the caseQ 2 Q 2 s 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. 4.4. In this regime, the dominant contributions to the integrals in Eqs. (4.7) In principle, one should cut the above integration at r ∼ 1/Q s , but since the Bessel function is exponentially decaying forQr Q s r 1, we can safely extend the integration back to infinity. The angular integration is trivial and keeping only the logarithmically enhanced term, we find It is easy to check that this dominant contribution arises from dipole sizes r ∼ 1/Q 1/Q s , thus justifying our above use of the weak scattering approximation. Up to a numerical factor, the quantity (Q 2 A /α s ) ln(Q 2 /Λ 2 ) represents the nuclear gluon distribution per unit transverse area generated by the valence quarks, as probed on a transverse resolution scaleQ 2 . So, Eq. (4.8) has a simple physical interpretation: the small projectile dipole with transverse size r ∼ 1/Q 1/Q s scatters with a probability ∼ α 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/Q 2 .
The same exercise can be repeated for the amplitude given the BK solution at asymptotically high energies, shown in Eq. (3.18). One finds where we omitted the overall proportionality constant since this is anyway not under control when writing the BK amplitude as in Eq. (3.18). We have here anticipated that the elastic cross section is more strongly suppressed in this limit. One can indeed check that The high-energy and large-virtuality results in Eqs. (4.9)-(4.10) show two remarkable features: (i) the "anomalous" dimension γ < 1 specific to the high-energy evolution (BK/JIMWLK) in the presence of saturation, and (ii) geometric scaling -the fact that these expressions depend upon the (effective) virtualityQ 2 = z(1 − z)Q 2 and upon the COM energy of the scattering only via the ratiō Q 2 /Q 2 s (η). Notice that geometric scaling also holds at low virtualitiesQ 2 Q 2 s (η), at least within the black disk approximation in Eq. (4.6). So, for sufficiently high energies (or large η), this is expected to be a robust feature of the integrated SIDIS cross-sections at anyQ 2 .  Figure 11: The normalized k-integrated transverse cross section-defined in Eq. (4.11) plotted as a function of z for various values of η at (a) Q 2 = Q 2 s0 ("lowQ 2 ") and (b) Q 2 = 5Q 2 s0 ("highQ 2 "). When z approaches 1, one sees a sharper increase in Fig. 11b than in Fig. 11a.
We have checked these properties against numerical solutions to the BK equation, with the results shown in Fig. 9b. With increasing η, one clearly sees a softening of the distribution at highQ 2 as compared to the low-energy (MV model) result Eq. (4.8). This softening reflects the emergence of the anomalous dimension, although the values of η ≤ 3 considered there are still to small to observe the asymptotic value γ 0.63 expected at very high energies.
To also check geometric scaling, we have plotted in Fig. 10 the results for the integrated crosssection d 2 k J T in terms of the scaling variableQ 2 /Q 2 s (η), up to the larger rapidity η = 6. In the low virtuality regime atQ 2 Q 2 s (η), one finds a quasi-perfect scaling (see Fig. 10a), in agreement with the leading-order estimate in Eq. (4.6). At higher virtualities,Q 2 Q 2 s (η), the scaling becomes better with increasing rapidity (see Fig. 10b) -again, as expected.

The z-dependence of the integrated cross-sections
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 2 s (η). In Fig. 11 we show our respective numerical results for two choices for the virtuality, Q 2 = Q 2 s0 (left plot) and respectively Q 2 = 5Q 2 s0 (right plot), and for η ≤ 3. The plotted quantity is the ratio where σ T (z) is merely a short-hand notation for the respective differential cross-section per unit of z, that is dσ T /dz. As compared to Fig. 9, we now implicitly concentrate on a more limited range of values for the ratioQ 2 /Q 2 s0 , 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.  Figure 12: (a) The ratio R pA for the SIDIS cross-section integrated over k, as defined in Eq. (4.12), and its rapidity evolution of up to η = 3. The ratio is always smaller than 1 in contrast to the unintegrated case (cf. Fig. 7). (b) The nuclear modification factor for the normalised cross-section introduced in Eq. (4.11) plotted as a function of z, for fixed Q 2 and for various values of η ≤ 3. As z approaches 1 the ratio drops significantly below 1, due to the different z-dependences observed in Figs. 11a and 11b, respectively.
This strong z-dependence can be understood from our analytic approximations in Eqs. (4.6) and (4.9). These expressions show that, when increasing z towards 1, σ T (z) increases like ln 1 1−z so long as Q 2 Q 2 s (η) and (roughly) like the power 1/(1 − z) γ whenQ 2 Q 2 s (η). For z 0.8, the two plots in Fig. 11 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. 11b ("highQ 2 ") than in Fig. 11a ("lowQ 2 "), and it becomes less pronounced when increasing η, thus reflecting the onset of an anomalous dimension γ < 1.

R pA ratios for the integrated cross-section
In Sect. 3.4 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 ⊥ , but a similar ratio can of course be defined for the cross-sections integrated over k ⊥ .
Specifically, in Fig. 12a we display the ratio as a function ofQ 2 = z(1 − z)Q 2 up toQ 2 = 5Q 2 s0 (A) = 10 GeV 2 and for different rapidities up to η = 3. This ratio would be equal to one in the absence of nuclear effects, but by inspection of Fig. 12a one sees that this is clearly not the case: the ratio is strictly smaller than one and monotonously increasing withQ 2 (albeit almost flat forQ 2 Q 2 s0 ) within the whole range under consideration. Moreover, the nuclear suppression is significantly increasing with η at any value forQ 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 [45,46]: the low momentum modes at k ⊥ < Q s get depleted and the gluons which are missing at low k ⊥ are redistributed in the modes with larger momenta k ⊥ Q s . Hence, if one integrates the k 2 ⊥ -distribution up to some valueQ 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 increasingQ 2 . With increasing energy (or η), 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. 12a. 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 Eq. (3.13), 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") [4,5].
Furthermore, Fig. 12b shows the nuclear modification factor for the quantity depicted in Fig. 11, which is itself a ratio, cf. Eq. (4.11). That is, the quantity depicted in Fig. 11a for the case of a nucleus is divided by the same quantity computed for a proton 15 . In both calculations, one uses Q 2 = Q 2 s0 (A) = 2 GeV 2 . Recalling that Q 2 s0 (A) = A 1/3 Q 2 s0 (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. 11b. In other terms, each of the curves shown in Fig. 12b can be thought of as the ratio of the 2 corresponding curves in Fig. 11a and Fig. 11b, respectively. From this perspective, the behaviour visible in Fig. 12b is easy to understand: when Q 2 = Q 2 s0 (A) = 6Q 2 s0 (p), one probes a dense gluon distribution in the nuclear, but a dilute one in the proton. Hence, the cross-section σ 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. 11). Accordingly, the nucleus-to-proton ratio is decreasing, as shown in Fig. 12b. Interestingly, this suppression is quite strong already for not too large values of z (say z 0.8) and thus it might be a robust signal of saturation in the EIC data.

Inclusive DIS: symmetric vs. aligned-jet configurations
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 thatQ 2 ≡ z(1 − z)Q 2 Q 2 s . Asymmetric dipole configurations, such that z is either very close to one, 1 − z 1, or very small z 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 σ T ≡ σ γ * T A (x, Q 2 ), cf. Eq. (4.5). 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.
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 2 s : (i) Low virtuality Q 2 Q 2 s : in this case, one clearly hasQ 2 Q 2 s 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 Λ 2 .
15 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 (4.11).
Consider first the "low-energy" case, as described by the MV model; using Eq. (4.8), one finds 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 2 s ). (A similar contribution from z ∼ 1 has been included via a factor of 2.) Clearly, this dominant contribution comes from the aligned-jet configurations, as previously defined -that is, z is small, but not too small (and similarly for 1 − z), in such a way that 1 z(1 − z) Q 2 s /Q 2 . Consider now the high-energy regimeᾱ s η 1, where the BK evolution becomes important. Using Eq. (4.9) and neglecting proportionality constants, we obtaiñ Since γ < 1, the above integral is dominated by large values 16 z ∼ 1/2 (in particular, the lower limit at z = Q 2 s /Q 2 1 played no role for the final result), i.e. by symmetric dipole configurations, with very small sizes r ∼ 1/Q. This change of behaviour w.r.t. the MV model calculation in Eq. (5.3) is due to the fact that, when γ < 1, the dipole amplitude for a single scattering rises more slowly with the dipole size r than in the regime of "colour transparency" (γ = 1).
The results in Eqs. (5.3) and (5.4) confirm our earlier statement in this section, namely the fact that the aligned-jet configurations dominate the DIS cross-section σ 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 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 η even faster than Q 2 s ) [46,58,59]). The above results also show that, irrespective of the precise value of γ ≤ 1, inclusive DIS at high Q 2 Q 2 s is controlled by relatively small dipoles undergoing weak scattering (r 1/Q s ). The corresponding contribution of the large dipoles with r 1/Q s , as estimated in Eq. (5.2), is parametrically smaller, in the sense that it lacks the logarithm enhancement visible in Eqs. (5.3) and (5.4). We thus conclude that the SIDIS configurations which are most interesting for studies of saturation yield a relatively small contribution to the inclusive cross-section, which is however sizeable -it is only logarithmically suppressed compared to the dominant contribution.
Before we conclude this section, let us notice that, when γ < 1, the parametric estimate shown in Eq. (5.4) for σ T also holds for the longitudinal cross-section σ L ; in that case too, this result is generated by symmetric dipole configurations with z ∼ 1/2 and r ∼ 1/Q. But unlike σ T , the longitudinal crosssection is controlled by symmetric configurations even when γ = 1, so in that case σ L is lacking the 16  large logarithm ln(Q 2 /Q 2 s ) visible in Eq. (5.3) and thus is suppressed relative to σ T .

Conclusions
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 1 (meaning 1 − z 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 ∼ 1/Q 1/Q, with Q 2 = z(1 − z)Q 2 . When moreover one has also hasQ 2 Q 2 s , 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 ⊥ of the final jet (or hadron), we have argued that the most interesting kinematics for a study of saturation is atQ 2 k 2 ⊥ Q 2 s . We identified two interesting phenomena in this regime. For sufficiently low transverse momenta k 2 ⊥ Q 2 s , 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 qq pair occurs after the photon has crossed the nucleus without interacting with it.
On the other hand, for larger transverse momenta k 2 ⊥ ∼ Q 2 s , 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 ⊥ ∼ 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 ratioQ 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 collinearlyimproved version of the BK equation (with running coupling) recently proposed in Refs. [48,49]. 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 ⊥ . Once again, the interesting regime for saturation is atQ 2 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 /Q 2 ), with the logarithm coming from integrating over the wavefunction of the virtual photon, that is, over the k ⊥ -distribution of the qq pair produced by the decay of the virtual photon. For the elastic scattering, this is also the final distribution. For the inelastic channel, the final distribution is generated via multiple scattering, so the integral over k ⊥ is controlled by values k ⊥ ∼ Q s .
To better emphasise the distinguished nature of this behaviour at lowQ 2 Q 2 s , we have also computed the integrated SIDIS cross-section in the weak scattering regime atQ 2 Q 2 s . We thus found the expected "leading-twist" suppression by the factor Q 2 s /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 Q 2 s (cf. Fig. 11): in the weak scattering regime atQ 2 Q 2 s one sees a rapid, power-like, growth when increasing z towards 1, which is eventually tamed (whenQ 2 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. 12b.
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 Q 2 s . Although both types of configurations are associated with dipole fluctuations which are very asymmetric in their z sharing (i.e. z(1−z) 1) and hence relatively large (r ∼ 1/Q 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 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 1/Q s ) even at high Q 2 Q 2 s . Their contribution to the inclusive cross-section is parametrically suppressed, but only logarithmically, so it should represent a sizeable fraction of the total cross-section measured in the experiments.
Throughout this analysis, we restricted ourselves to the leading-order pQCD approximation, where the particle "measured" in SIDIS at large z is a bare, on-shell, quark. In general, the quark emerging from the scattering with the hadronic target will have some time-like virtuality, so it will evolve via final-state radiation, thus giving rise to a jet. In an actual experiment, one could measure either a jet with a large value of z (defined as the total longitudinal momentum fraction carried by the jet constituents), or a single, "leading", hadron, which itself has a large value of z.
Leaving aside the non-perturbative effects, like hadronisation or higher-twist corrections, it is important to observe that already the higher-order pQCD corrections are expected to be quite different in the two cases (jet production and single hadron production, respectively). Indeed, for an exclusive process like hadron production, the limit z → 1 introduces a strong constraint on the radiation in the final state: the total energy fraction carried by the unmeasured emissions cannot exceed the small value 1 − z. Such a constraint is well known to introduce special higher-order corrections, where each power of α s is accompanied by the double logarithm ln 2 1 1−z . (See App. F for an explicit, albeit admittedly sketchy, calculation at next-to-leading order.) Physically, these large corrections express the fact that it is very unlikely to observe a final state in which most of the total energy of the incoming photon is carried by a hadron alone. The existence of these double-logarithmic corrections spoils the validity of any fixed-order approximationa fortiori, of our present estimates at leading order.
However, the situation becomes much simpler if in the final state one measures a jet. Unmeasured emissions which are relatively energetic (with energy fractions larger than 1 − z) are now allowed, so long as they belong to the forward jet. The final state is now composed with more likely configurations and, as a result, there are no double-logarithmic corrections anymore. (Technically, this is ensured by an efficient compensation between "real" and "virtual" emissions; see App. F for an example at NLO.) In view of this, our present results can be seen as a bone fidae leading-order approximation to the SIDIS cross-section for forward jets with large z 1. On the other hand, much more workincluding all-order resummations of the double-logarithmic corrections -is needed to give a reliable estimate for the single hadron production in SIDIS at large z.
Still for the case of forward jet production, one can also worry about the importance of the highertwist corrections -the non-perturbative corrections suppressed by inverse powers of the jet virtuality. As already stressed, both final jets have relatively small transverse momenta k 2 ⊥ (1 − z)Q 2 Q 2 (when z 1). If this transverse scale were to be representative for the jet virtuality, then indeed the twist-expansion would be questionnable. However, this is not the case for the measured jet, with a large energy fraction z ∼ 1. By inspection of the kinematics, one can deduce that this jet should be time-like and relatively hard, with an invariant mass squared M 2 z ≡ k µ k µ of the order of the initial photon virtuality: M 2 z ∼ Q 2 . This follows from the appropriate generalisation of the kinematical constraint (2.8), as obtained by replacing 17 k 2 1⊥ → k 2 1⊥ + M 2 z : after also using k 2 1⊥ k 2

2⊥
(1 − z)Q 2 when z 1, one sees that the limit z → 1 puts no special constraint on the jet invariant mass M 2 z , which therefore is of order Q 2 .
The fact that the produced jet is hard when Q 2 Q 2 s also ensures that the SIDIS process at large z is leading-twist in the sense of the QCD parton picture -it measures a genuine quark distribution in the hadronic target -even whenQ 2 = z(1 − z)Q 2 is low enough,Q 2 Q 2 s , for the saturation effects to be important. In other terms, the saturation effects that we discussed in this paper can be characterised as leading-twist quark shadowing. This may appear as a surprise since in our calculations the saturation effects manifest themselves as multiple scattering -a phenomenon generally associated with higher-twist corrections. Recall however that the standard twist expansion for DIS is formulated in a special frame (a target infinite momentum frame, such as the Bjorken frame) and a special gauge (the target light-cone gauge, that is, A − = 0 in the present conventions), which are not the frame and gauge employed by the dipole picture that we have used throughout. It would be conceptually interesting to repeat some of the previous calculations in the traditional QCD parton picture, to verify our above statements about the twist expansion and, especially, to understand the role of gluon saturation and of the unitarity corrections in that framework.
The previous remarks clearly show that our present study is only preliminary. Further theoretical studies are needed, especially for the case of forward hadron production. On the experimental side, one still need to ascertain the 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.
We have performed the angular integration using the first equality in Eq. (A.7). In the above T = 1−S is the MV model amplitude, with S given in (3.5). We shall follow the procedure developed in [46] (cf. Sect. 2.2 and Appendix A there), which is in fact applicable to a wide class of "observables" in the MV model. Using Eq. (3.6) to replace Q A in terms of Q s and after doing simple algebraic manipulations we can rewrite S as In writing the second, approximate, equality we have done an expansion to first order in 1/(ln Q 2 s /Λ 2 ) which is a small number, roughly equal to the QCD running coupling evaluated at the saturation momentum. Such an expansion is clearly valid for r 1/Q s and the only potential dangerous regime could be when r 1/Q s , that is, very close to the unitarity limit. However in that regime the Smatrix is anyway close to zero and thus the error in our expansion is innocuous. It is important to emphasize that this is not a twist expansion, since the first factor clearly contains multiple scattering to all orders. Rather it should be viewed as an expansion in r around 1/Q s in a well-defined and controlled approximation scheme. The two terms in Eq. (B.2) give rise to two contributions to ∇ k W, more precisely These saturation and twist pieces, with the names to be shortly justified, correspondingly read and The integration in the saturation contribution is easily done and we get while the twist contribution can be given in terms of the exponential integral function Ei(x) = PV x −∞ dt e t /t, where PV stands for the Principal Value, and reads in Eq. (3.16), rather we shall take a step back and write [61] T (Expanding the above to second order for small ρ indeed leads to Eq. (3.16).) Inserting the above into Eq. (C.1), reversing the order of integrations and using Eq. (A.1) we find where in writing the second approximate equality we have just combined the two fractions and used k 2 ⊥ Q 2 . The above contains two logarithmic regimes of integrations. The "harder" one is for Q p ⊥ k ⊥ , with p ≡ q − k and from only the first term in the first equality above we find Note that the power 1/k 4 ⊥ has been generated by the exchange (a hard scattering with transfer momentum q ⊥ k ⊥ ), whereas the logarithm came from the integration (d 2 p/p 2 ) over the wavefunction of the virtual photon.
The second, "softer", regime is for Λ q ⊥ k ⊥ and one must be very careful in expanding the fraction in the integrand after the second equality in Eq. (C.3) to keep all the terms that eventually lead to O(1/k 4 ⊥ ) contributions in J T . Notice that the leading 1/k 2 ⊥ contributions from the two terms in the first expression in Eq. (C.3) cancel each other in the limit of interest, and this is why we prefer to work with the second, equivalent, expression. One has After doing the angular integration the last term vanishes while the first one acquires an averaging factor of 1/2, so overall Eq. (C.5) is equal to q 2 ⊥ /k 2 ⊥ . The ensuing integral over q generates the Coulomb logarithm characteristic of soft scattering. We finally deduce the following contribution to J T in Eq. (C.3): As it should be clear from the above manipulations, the power 1/k 4 ⊥ now originates from the original transverse momentum distribution 1/(k 2 ⊥ +Q 2 ) generated by the decay of the virtual photon -more precisely, it arises as the difference between two such distributions, with one of them shifted with the momentum q transferred by the scattering.

D Collinearly improved BK equation
The BK evolution equation we have used for the numerical calculations in this work reads [48,49] (see also [50,51]) S xz (η−δ xz;r )S zy (η−δ zy;r ) − S xy (η) (D.1) and we briefly explain below all the various elements in the above. The gluon is emitted at the transverse position z by either the quark or the antiquark of the parent dipole located at (x, y) and the indices in the S-matrices stand for its dependence on the corresponding dipole coordinates. When compared to the LO BK equation we note various differences which all arise from the resummation of higher order corrections enhanced by large logarithms.
(i) The rapidity shifts in the arguments of the S-matrices for the daughter dipoles resum double logarithms associated with the time ordering in the successive emissions and they are given by with r ≡ |x − y|, and similarly for δ zy;r . They become significant for emissions in which one of the daughter dipoles is much smaller than the parent one.
(ii) The factor r 2 /z 2 ±A 1 , wherez ≡ min{|x − z|, |z − y|} and with A 1 = 11/12, includes the resummation of the first set of single DGLAP logarithms. The sign in the exponent ±A 1 is taken to be plus for r 2 <z 2 and minus r 2 >z 2 , so it is always suppressing the evolution.
(iii) The use of a running coupling resums the corrections related to the one-loop QCD β-function. When one of the three dipoles is much smaller than the other two, pQCD requires that the scale should be the smallest dipole size, i.e. one should useᾱ s (r min ) in Eq. (D.1), where r min ≡ min{|x − y|, |x − z|, |z − y|}. Here, we shall use what is close to a BLM prescription (called "fast apparent convergence" in [52]), which minimizes the NLO correction to the BK equation proportional to the one-loop β-function. It is defined as withb 0 = 0.75 (corresponding to 3 flavors) and Λ = 0.2 GeV. The Landau pole is avoided by the introduction of r * according to r * = r/ 1 + r 2 /r 2 max with r max = 4 GeV −1 , which freezes the coupling in the deep IR to the valueᾱ s (r r max ) 0.728. Finally, we point out that Eq. (D.1) is solved as an initial value problem, but one must pay particular attention to the non-locality in η. The shift in Eq. (D.2) introduces a dependence to rapidities smaller than η and thus when we start the evolution at a rapidity η 0 (which is taken to vanish throughout the current paper) we need to specify the initial condition for η ≤ η 0 . We shall assume a constant behaviour in η (cf. the discussion in Sect. 9 of Ref. [48]), that is S xy (η < η 0 ) = S xy (η 0 ) (D.5) and the latter is simply given by the MV model as written in Eq. (3.9).

E Useful integrals
First, we list some integrals which are useful in calculating the high-momentum tail of the k-dependent cross sections. Generally, for −1 < γ < −1/2, one has By analytic continuation we shall take the above to be true for values of γ outside the aforementioned interval. Differentiating with γ we can also calculate integrals involving a power of ln ρ 2 . For example we have which holds for any γ > 0. In particular for γ = 1 and γ = 2 we respectively get 2/3 and 8/5.

F Double-logarithmic corrections: hadron measurement versus jet measurement
In the concluding section, we have mentioned an important difference between measuring a quark ("hadron") and, respectively, a jet, in a SIDIS experiment where the longitudinal fraction z of the measured system is close to one. In the first case -hadron measurement -, one expects large next-to-leading order (NLO) corrections enhanced by a double-logarithm ln 2 1 1−z , which make the leading-order predictions unreliable. In the second case -jet measurement -, such large corrections are however absent. In this Appendix, we would like to make this statement more precise, by presenting an admittedly sketchy calculation of the relevant NLO corrections, to the double-logarithmic accuracy of interest. We plan to present this calculation in full detail in a subsequent publication.
In the case of a quark measurement, that we shall address first, the large NLO corrections of order α s ln 2 1 1−z are generated by late gluon emissions, which are a part of the familiar DGLAP evolution of the final quark. More precisely, they are the result of a kinematical mismatch between "real" and "virtual" gluon emissions, as represented by the graphs in Fig. 13a and Fig. 13b, respectively. Namely, the momentum of the measured quark is equal to k for the "virtual" graph 13b and to k − for the "real" graph 13a. Here, is the gluon 4-momentum and is integrated over (since the gluon is not measured in the final state). In both cases, the longitudinal momentum of the final quark -that is, k + for Fig. 13b and respectively k + − + for Fig. 13a -is measured to be a fraction z of the initial momentum q + of the virtual photon.
As generally throughout this paper, we use the (time-orderd) light-cone perturbation theory together with the projectile LC gauge A + = 0. In particular, all the (real or virtual) quanta are on-shell and the lack of energy conservation at the emission vertices translates into energy denominators which encode the duration of the quantum emission process. Also, the LC longitudinal momentum of each particle is positive and bound by the respective momentum of its parent. For the case of the "real" emissions, it is furthermore constrained by the kinematics of the final state.
Accordingly, the longitudinal momentum of the "virtual" gluon in Fig. 13b is constrained by + ≤ k + = zq + , whereas that of the "real" gluon in Fig. 13a rather obeys + ≤ (1 − z)q + . (Indeed, once the quark is measured with zq + , the total longitudinal momentum carried by the unmeasured antiquark and gluon is known to be (1 − z)q + , and this represents the upper bound on + .) This kinematical mismatch lies at the origin of the double-logarithmic NLO corrections, as we now explain.
The calculation of the DGLAP-like NLO corrections in the double logarithmic approximation (DLA) being a standard exercice, here we shall merely exhibit the final result. These corrections factorise, so they can be fully encoded in a factor, F real or F virt , multiplying the leading-order SIDIS cross-section. This factor has a similar structure for the real and virtual graphs: in both cases, it features the integral over the phase-space of the emitted gluon, in logarithmic units for both the longitudinal ( + ) and the transverse ( ) momentum. However, F real and F virt differ from each other in their integration limits over + (as already discussed) and also in their global sign: the virtual corrections are negative, as they describe a reduction in the probability to find a "bare" quark in the final state.
Specifically, the factor F virt encoding the effect of the "virtual" graph in Fig. 13b is found as The longitudinal integration being by now quite obvious, let us give some explanations about the structure of the transverse integration. First, we have introduced a non-zero quark mass m, to screen the would-be collinear divergence in the limit where the gluon emission angle relative to its parent quark approaches to zero: θ k → 0 with θ 2 k ≡ + − k k + 2 . A non-zero mass m limits this angle to a minimal value θ k ≥ m/k + , via the well-known "dead cone effect": gluon radiation within an angle smaller than m/k + around the direction of the emitter is forbidden by the kinematics. The other important ingredient in Eq. (F.1) is the step function enforcing the formation time τ = 2 + / 2 of the gluon to be larger than the coherence time τ q = 2q + /Q 2 of the virtual photon. This is the condition that the gluon emissions under consideration correspond indeed to final-state radiation. One can check that early radiation, with τ < τ q , does not generate collinear divergences, hence it does not contribute to DLA. After a change in the integration variable, → + + k + k, and to the DLA accuracy of interest, Eq. (F.1) can be further simplified to As compared to Eq. (F.1), we have also approximated k + = zq + q + in the integration limits, in view of the fact that 1 − z 1. The corresponding "real" contribution is given by a similar expression, except for a change of sign and for the much lower value of the upper limit on + : Clearly, the net contribution F quark ≡ F real + F virt reads (the lower script "quark" is meant to recall that we measure a quark in the final state) The final result in Eq. (F.4) exhibits the anticipated double-logarithmic correction, which becomes arbitrarily large in the limit z → 1. Notice that this correction is independent of the transverse momentum k of the measured quark, so it would be the same for both the double differential crosssection dσ dz d 2 k and the cross-section dσ dz integrated over k. To summarise, this large NLO correction to the SIDIS cross-section for quark (more generally, hadron) production occurs because of the large disparity between the longitudinal phase-spaces for "real" and respectively "virtual" gluon emissions. In turn, this disparity is a consequence of the fact that the "real" gluon is not measured in the final state, so its longitudinal momentum must be very low, + ≤ (1 − z)q + , in order for the measured quark to carry a large longitudinal fraction z 1 by itself. However, if instead a quark, we measure a jet, then the "real" gluon can be a part of that jet and its longitudinal momentum is not constrained anymore. This will be explained in what follows.
When measuring jets in the final state, we would like to be able to distinguish between a hard jet with longitudinal fraction z 1 and a soft one, with fraction 1 − z 1. To that aim, the angular opening of the hard jet must be limited by the angle δ between the 2 jets, that can be estimated as . (F.5) Concerning the "real" gluon emissions, it is convenient to distinguish between soft emissions, with + < (1 − z)q + , and hard ones, with + > (1 − z)q + . (The intermediate values with + ∼ (1 − z)q + do not contribute at DLA.) The soft emissions can have relatively large emission angles, hence for them the jet cone condition θ k < δ may be quite restrictive. However, such soft emissions do not count for the jet energy balance, for none of the 2 jets, so in practice they are not observable and can be neglected. The hard emissions on the other hand must belong to the hard jet. Recalling that θ k = ⊥ / + (we use the ⊥ variable occurring in Eq. (F.2)), we have the condition ⊥ + < δ or 2 ⊥ < ( + ) 2 What we would like to show is that this constraint is in fact automatically satisfied by the "real" emissions with + > (1 − z)q + . Indeed, the upper limit on 2 ⊥ in the second inequality above is in fact larger than the upper limit on the respective integral in Eq. (F.3), as we now verify: the inequality that we need, namely, is equivalent to which is obviously satisfied since k 2 ⊥ ≥ Q 2 z(1 − z) Q 2 (1 − z), and + (1 − z)q + . Accordingly, when going from a quark to a jet measurement, the integral over "real" gluon emissions extends up to k + q + , so the virtual contributions (which are still given by Eq. (F.2)) are exactly cancelled by real emissions of gluons lying within the jet cone (F.5) -at least to the double-logarithmic accuracy of interest. In other terms, they double-logarithmic corrections found in Eq. (F.4) for the case of quark production are absent in the case of jet production.