Forward dijet production at the LHC within an impact parameter dependent TMD approach

We investigate possible signatures of gluon saturation using forward $p+A \to j+j+X$ di-jet production processes at the Large Hadron Collider. In the forward rapidity region, this is a highly asymmetric process where partons with large longitudinal momentum fraction \(x\) in the dilute projectile are used as a probe to resolve the small \(x\) partonic content of the dense target. Such dilute-dense processes can be described in the factorization framework of Improved Transverse Momentum Distributions (ITMDs). We present a new model for ITMDs where we explicitly introduce the impact parameter (\(b\)) dependence in the ITMDs, to properly account for the nuclear enhancement of gluon saturation effects, and discuss the phenomenological consequences for $p-Pb$, $p-Xe$ and $p-O$ collisions at the LHC. While the case of $p-p$ and $e-p$ collisions is used to fix the model parameters, we find that, on average, the nuclear enhancement of the saturation scale is noticeably weaker than expected from naive scaling with a simple dependence on the atomic number. Since our model explicitly accounts for event-by-event fluctuations of the nuclear geometry, it can also be applied to study forward central correlations in $p-A$ collisions.


Introduction
QCD predictions in the high-energy regime as well as the precise determination of the parton distribution functions (PDFs) have been one of the main goals at the HERA, Tevatron and now LHC colliders.The observation of a steep rise of the gluon density in the proton towards small Bjorken-x at HERA [1] led to the concept of saturation where some recombination (non-linear) mechanisms that tames the (linear) growth in gluon multiplicitly should set in at some values of x and Q 2 to preserve unitarity.The gluon recombination mechanism is best described in the framework of Color-Glass Condensate (CGC) [2], where the linear and non-linear corrections are formalized into the JIMWLK/BK [3][4][5][6][7][8][9] evolution equations.As a consequence of the interplay of these two competing mechanisms of linear emission and non-linear absorption, a dynamical transverse momentum scale arises (see for instance [10]).The understanding of the mechanisms that dynamically give rise to the energy scale of saturation Q s (x) is related to dynamics of partons inside hadrons.
Forward di-jet production is one of the preferred processes where the effects of the dynamics of saturation can be investigated at hadron colliders, and has been the subject of numerous theoretical [11][12][13][14] and phenmenological studies [15][16][17][18][19][20][21] (see e.g.[22] for recent reviews).The detection of jets in the very forward region in the same hemisphere of a typical detector at the LHC imposes a dilute-dense asymmetry where the interacting parton in one of the (dilute) hadrons carries a large (to intermediate) fraction x of the momentum of the hadron and is used to probe the (dense) gluon distribution of the other hadron at much smaller values of x (see Fig. 3.1 in Sec.3.1 for an illustration of such a configuration).The process is sensitive to saturation effects when the transverse energy of the small-x gluon, which is controlled by the transverse momentum imbalance between the two jets k r , is comparable to the saturation scale.The interested reader can also find di-jet production analyses dedicated to saturation physics in the context of DIS at e.g.[14,23,24].
In this paper, we will compute forward di-jet production cross-section using the unified formalism given in Ref. [25] that allows to interpolate between the high-energy factorization (HEF) framework [26] valid in the regime where the momentum imbalance k t of the jets is of the order of the transverse momentum p j of the jets k t ≃ p j and the Transverse-Momentum-Distribution (TMD) framework [27,28] valid when the momentum imbalance of the two jets is on the order of the saturation scale k t ≲ Q s .This Improved-TMD (ITMD) framework has been applied to obtain predictions for forward dijet production at high rapidity in p − p and p − A collisions [17].The authors argued that the larger nuclear density in large nuclei results in an earlier onset of the nonlinear regime in p − P b versus p − p collisions as the jet geometry approaches the "back-to-back" configuration (k t ≪ p j ).A significant reduction in jet azimuthal correlation appears in p − A versus p − p collisions.The nucleus distribution was extrapolated from the proton distribution by simply imposing the expected mass number A scaling, thus treating the nucleus as a uniform object resembling a large proton.
In this paper, we discuss a more accurate approach to describe the nucleus as an incoherent collection of nucleons.In order to account for these effects on the saturation dynamics, we developed a simple model in which the nuclear TMDs can be computed from the proton TMDs, by promoting the distributions to take into account the impact parameter dependence, given the spatial distribution of nucleons inside the nuclear target.
The starting point of our model is an accurate description of the nuclear matter distribution at small-x.DIS data at HERA are used as the main source of information for the proton structure because of the high precision of the measurement [1,29,30] and of the simplicity of the theoretical interpretation in terms of the dipole picture [31][32][33][34]34,35].In the dipole picture, DIS cross sections are written in a factorized form that separates the transition of a (virtual) photon to a quarkantiquark (q q) dipole from the quark-antiquark dipole scattering off the dense nuclear matter in the proton target.The nuclear matter distribution cannot follow from a pure perturbative description.The small-x evolution can be computed using perturbative techniques but the nuclear configuration at the starting scale of the evolution must be extracted from experimental data.
In our model we include the small-x evolution following the procedure of Ref. [36].The authors (AAMQS) provide a good fit of the most recent HERA data [1] including c and b heavy quark contents, where the small-x evolution is driven by the BK equation supplemented with a subset of the next-to-leading corrections to account for the running of the coupling (rcBK).Such a choice represents a good compromise between the simplicity of its implementation and its accuracy.A more rigorous treatment would entail a renormalization in terms of the BK-JIMWLK equations which, unfortunately, are notorious for being hard to implement numerically.The performance of the fit gives further confirmation of the validity of this approach.
The AAMQS fit involves only quantities averaged over impact parameter since it is based on the HERA measurements that were not dependent on the centrality of the scattering.A rigorous treatment of the impact parameter dependence of the dipole amplitude requires the solution of the rcBK equation beyond the translational invariant approximation [37][38][39].Such study is complicated by non-perturbative effects that appear at large impact parameters where the onset of the physics of confinement emerges.Even though a satisfactory treatment of the impact parameter dependent evolution remains elusive, several attempts [40][41][42] have been made to introduce the impact parameter dependence into a consistent physical description without having to meddle with the intricacies of the b-dependent evolution equations.
One prominent example is the Impact Parameter Saturation (IPSat) model [40], which is based on the particularly simple Golec-Biernat and Wusthoff (GBW) [43] dipole saturation model, but contains on one hand, a DGLAP extension towards higher Q [44] and, on the other hand, an impact parameter dependence in the dipole dynamics.A key advantage of the IPSat model is the fact that it is immediately generalizable from protons to nuclei, such that the model has frequently been employed in studies of e − A, p − A and even A − A collisons [45][46][47][48].On the other hand, as originally noted by the authors, an obvious shortcoming is the model inability to describe the lower end of the x evolution since the gluon structure function is only evolved through the (lowest order) DGLAP evolution kernel.
Our model which was inspired by the success of the IPSat model and of the AAMQS fit with rcBK evolution can be seen as an hybrid description combining the two.In short, we propose to keep the impact parameter dependence of the dipole amplitude N in the eikonal exponentiated form as done in the IPSat model but we substitute the DGLAP evolution of IPSat with the rcBK evolution extracted from the AAMQS fit.In this way, the dipole amplitude can then be expected to be more reliable toward small-x while at the same time, capturing some of the effects that the b-averaged AAMQS description lacks.Since in our study the (transverse) geometry of the nuclei is modeled using the rather sophisticated MCGlauber model [49] that is frequently used in heavyion physics, one can readily anticipate that the model predictions for nuclei deviate significantly from the "naive" scaling with the atomic number.Notably, this more detailed treatment will also allow for future extensions of the model to include e.g.centrality dependence or study correlations between forward di-jet production and event activity in the central rapidity region.
The structure of the paper is divided in four sections.Section II starts with a summary of the AAMQS fit model and it is also used to set a common physical framework and notations used in the following sections.Then, what follows is dedicated to a description of our model.It presents the impact parameter dependence and the values of the free parameters of the model.To that end, we compare the dipole TMDs computed from the dipole AAMQS fits and the b-dependent TMDs calculated within our model, averaged over impact parameter.Finally, in section III, we use our model to make predictions for di-jet production in p − p and p − A collisions as a probe of saturation effects.Concluding remarks are outlined in Sec.IV.

Model Description of Nuclear TMDs
Since our model was designed with the goal of computing the di-jet production cross section in the asymmetric dilute-dense regime in p − p and p − A collisions, the task at hand is to infer the various nuclear TMDs that enter the calculation of the ITMD factorized cross-section of Eq. 3.3 (see also [25]).The main interest in these processes derives from the fact that they represent an almost ideal probe to study the onset of saturation at hadron colliders.Since the saturation scale is expected to increase with the atomic number A as , it has often been argued that the measurement of the di-jet distributions in p − A and p − p collisions will be sensitive to the differences in saturation scales [50].However, in order to make precise predictions, one needs to properly model the distribution of nuclear matter inside the nucleus that induces the rise of the saturation scale in nuclei.The challenge is thus to find a consistent way to generalize the available knowledge for protons to arbitrary nuclei.
Our general strategy for tackling this problem can be briefly summarized as follows.Starting point of discussion are the HERA measurements of the proton structure functions F 2 , F L , which can be used [36,47,48] to constrain the (average) dipole amplitude of the proton down to very low values of x.We obtain this information about the average dipole amplitude N (x, r) from the AAMQS fit and subsequently generalize the model to account for impact parameter b dependence, by introducing a dependence on the local thickness T p (b).Based on this procedure, we obtain an expression for the dipole amplitude N (x, r, b), which can simply be generalized from protons to nuclei, by replacing the local thickness of the proton T p (b) with the nuclear thickness T A (b) obtained as the sum of the thickness of all nucleons inside a nucleus T A (b) = i∈A T i (b).By assuming Gaussian fluctuations of color fields inside the nuclei, we then derive a set of relations, whereby the various TMDs relevant for forward di-jet production can all be expressed in terms of the dipole amplitude N (x, r, b).Since this calculation is technically rather involved, we defer the details to Appendix C, and simply note that our results generalize earlier results of Ref. [51] for some of the TMDs.Based on these relations, we can then perform a sampling of the spatial distributions of nucleons inside nuclei using the MC Glauber model, and calculate the x and k t dependence of the various TMDs by integrating over the impact parameter b and averaging over the spatial distribution of nucleons.Details on the individual steps are provided in the following subsections followed by a brief discussion of our results for the relevant TMDs for proton and nuclear targets.

Deep-inelastic scattering (DIS) in the Dipole Picture & AAMQS fit to HERA data
In the dipole picture used to describe DIS physics at small-x, the interactions between the highenergy photon and the target nuclear matter (the proton) are described in terms of the light-cone wave function Ψ f T,L , which quantifies the probability for a virtual photon to fluctuate into a quarkantiquark dipole pair of flavor f , and the dipole amplitude N , which describes the interaction of a q q color dipole with the hadronic target.The longitudinal and transverse components of the cross section can be given in the factorized form where m f , e f are the mass and charge of the quark flavors in the dipole, Q2 the photon virtuality, z the longitudinal momentum fraction carried by the quark (see Fig. Generally speaking, models of saturation physics aim to extract the behavior of the dipole amplitude N (b, r, x) from measurements of the inclusive cross-section (and in some cases more exclusive processes [41,45]) at HERA.Since our calculation of small x TMDs will be based on the extraction of the dipole amplitude of Ref. [36] (AAMQS), we will now briefly discuss the simplifications made in the AAMQS analysis.Since the AAMQS analysis was based on inclusive DIS, only impact parameter averaged quantities are considered.Introducing (half) the average of the quark distribution in the transverse plane σ p , the cross section in Eq. 2.1 is expressed in terms of the averaged dipole amplitude N (r, x) through the substitution (2.2) where σ p was left free in the AAMQS fit analysis and led to the value 3.54fm2 for the light quark data set.We keep it fixed to the same value in our analysis.For simplicity, the heavy quark (charm and bottom) data set has been excluded from our analysis.Within the analysis, the averaged dipole amplitude N (x, r) was evolved starting from a nonperturbative initial condition modeled following the semi-classical Mueller-Venugopalan model where Q s0 , Λ, γ are additional fitted parameters and e is the Euler constant.Eq. 2.3 gives the "exponentiated" form of the dipole amplitude originating from the multiple rescatterings of the dipole passing through a dense target shock-wave [40].Starting from the initial scale x = 0.01, the evolution towards smaller x was performed using the Balintsky-Kovchegov (BK) equation [3,4] with running coupling corrections (rcBK) [37,55].This evolution represents the dynamical input based on the QCD perturbative predictions.The explicit form of the rcBK integral-differential equation can be found in Eqs.(2.7, 2.8) of Ref. [36].
We use the FORTRAN routine provided in Ref. [36] to compute the dipole amplitude N (x, |r|) and generate a data set mapping for any chosen value of x ∈ (10 −12 , 0.008 of the dipole sizes r = |r| with the corresponding dipole distributions N (r, x).The available dipole sizes span several order of magnitudes r i ∈ 10 −8 , 10 .
Since the calculation of the TMDs involves taking derivatives over r, a function N * (r, x n ), where the subscript n(= 1 . . .8) enumerates the x sampled values was fitted to each data set corresponding to a fixed value of x = x n .The choice of the function N * used for the fit was inspired by the M V γ form of the dipole amplitude 2 We note that in the above expressions, the dipole size r is rescaled as r = r/r s , where r s is the inverse of the saturation scale, i.e.
which is computed for each data set, by following the standard procedure in solving the implicit equation By interpolating the results based on this parametrization, the exponent of the dipole amplitude G(r) = G{q n , λ n , γ n , e cn }(r) is then given as a smooth function of r that depends on a family of dimensionless parameters 3 .We note that the quality of the fit in Eqs.2.4 remains within 10% of data throughout the full range of dipole sizes.

Impact Parameter Dependence
Evidently, in the proton case, the fit to HERA data is of good quality which means that the physics of inclusive DIS can be well described without taking into account any impact parameter dependence.However, we expect the impact parameter dependence to play an important role in the context of p − A collisions.In nuclei, the nuclear densities can vary greatly according to the nucleon positions resulting in an irregular shape with lumps and valleys whose effects are hard to capture with a simple scaling correction factor.Hence, we will now discuss how to introduce the b-dependence into the description starting from the b-averaged rcBK fit.
Inspired by the simple but successful IPSat Model [40], we include a proportionality to the nuclear matter density in the exponential, as it is done for the Glauber-Mueller formula [33] for the dipole cross section.Starting from Eq. 2.4, we thus impose the following substitution where the nuclear thickness T (b) is the two-dimensional projection of the nuclear density ρ(z, b) over the transverse impact parameter space (c.f.Appendix A), S p ⊥ = σ p /2 is the proton transverse area fixed to the AAMQS value.σ, η and T cut are model parameters to be fixed according to the procedure detailed in Sec.2.4.Note that having set an overall S p ⊥ factor in Eq. 2.6 the parameter η is dimensionless while [σ] = fm 2 .Specifically, for the proton we follow the IPSat model and employ a Gaussian profile with a nucleon (square) width B G = 4GeV −2 [41].When considering nuclei, we simply replace the T p → T A without modifying the parameters of our model, where the nuclear thickness T A is obtained as the sum of the thickness functions T p/n of all nucleons 4 .Due to the fluctuating positions of the nucleons in the nucleus, which we sample according to a MC Glauber procedure, as explained in Appendix A, the nuclear thickness fluctuates on an event-by-event basis.The final results for the TMDs are obtained taking the statistical averages of the TMDs computed from typically ≈ 10 nucleon configurations.We note that a lower bound T cut > 0 on the nuclear density is necessary in Eq. 2.6 to recover a finite result once the d 2 b integral is computed.Physically, a finite cut-off can be justified since a description in terms of saturation physics is valid only when the dipole interacts with a substantial amount of nuclear matter, and thus has a high probability of subsequent multiple rescatterings.Due to this cut-off, it is necessary to introduce a further dimensionless parameter η, so that the bdependent model correctly reproduces the average proton TMDs upon integrating over the impact parameter.Since the behavior of G(r, x) is readily fixed by the AAMQS fit, the parameters σ, η and T cut are then the only free parameters of our model, which can be fixed by requiring that the proton TMDs calculated in our model match the proton TMDs computed from the average dipole amplitude in Eqs.2.4.

Calculation of small x TMDs
Based on the knowledge of the parametrization of the dipole amplitude in Eq. 2.6, we can proceed with the calculation of the small x TMDs that enter the di-jet production cross section.Since, in the small x limit, the TMDs can be formulated in terms of correlation functions of light-like Wilson lines [56][57][58], it is possible to express the various TMDs in terms of the exponent G(x, r, b) of the fundamental dipole amplitude, if one invokes the additional assumption of Gaussian statistics for the distribution of color charges inside the nucleus, that is common to almost all saturation models (see e.g.[59] for a discussion).Denoting the two end points of the dipole as x, y, such that x + y 2 one finds for example, that the first TMD F qg in the quark-gluon channel -which corresponds to the usual fundamental dipole gluon distribution -can be expressed, in (transverse) coordinate space, in terms of G x,y and its partial derivatives G (i,j) x,y (x) 5 as F (1)  qg (x, y) = e G (0,0) x,y up to a factor of c = Nc αs(2π) 2 π , which we omitted from the definition of the TMD in (transverse) coordinate space and also exclude from all plots of TMDs in the following.Similar expressions can also be derived for all the other TMDs, following the procedure outlined in Appendix C. Since the calculation is rather cumbersome, and the resulting formulae are not particularly enlightning, we refrain from reproducing all the expressions in the main part of this paper, and instead refer the interested reader directly to Appendix C, where we provide the explicit expressions of all TMDs relevant to our analysis.
Based on these expressions, we follow the same logic as in Eq. 2.6 and employ to calculate the various TMDs F in coordinate space, by means of the substitution in the respective expression.Subsequently, the TMDs in momentum space are obtained by a simple Fourier transform, i.e.
where the transverse momentum k is the Fourier conjugated variable to the dipole size r.Since, for unpolarized nuclei, the distributions are azimuthaly symmetric, this amounts in practice to calculate the Bessel integral which we obtain by means of numerical integration.

Fixing the Model Parameters
Now that we have established the procedure to calculate the various TMDs, we will discuss how the parameters of our model are fixed by requiring that the b-dependent proton dipole TMD F (1) qg reproduces the rcBK fit once it is averaged over impact parameter space.In the following, the former and the latter will be indicated as p and p respectively, such that the matching condition takes the form F (1) where ) Based on this requirement, the parameters σ and η can be chosen so that the perturbative large k t "tails" of the proton dipole TMDs match as shown in Fig. 2.2.By looking at the explicit expressions for the TMDs (2.13a) one realizes that in the large k limit, the exponentials can be approximated as unity, giving rise to the following matching conditions We note that T p is the average thickness value, while the parameter A p ≲ 1 because of the finite cut-off T cut .Based on our discussion, a reasonable value for T cut is provided by Since this choice is to some extent ad-hoc, the model uncertainty was estimated by varying this cut-off value by ±50%.We note that -as can be easily understood from their definitions (see Eq. 2.14a) -the model parameters σ and η must be recomputed when T cut varies.Their explicit values are collected in Table 1.
We include this model uncertainty in all figures, where the colored bands indicate the uncertainties due to the cut-off value.We verify this procedure in Fig. 2.2, where we show the dipole TMDs for both the b-dependent and b-averaged approaches for three different values of x.We find that, as can be expected from the matching procedure, the p and p TMDs show small differences at low values of k t ∼ Q s (x), but are in excellent agreement in the perturbative large k t ≫ Q s (x) region.
Similar behavior can also be observed for the other TMDs that enter the calculation of the di-jet production cross-section, which are compactly summarized in Fig. 2.3, where we present the results for the seven relevant proton TMDs. 6Notably -with the exception of F (2) gg (c.f.discussion in [25]) -all the TMDs exhibit the same perturbative large k t /Q s (x) behavior, which is in excellent agreement between the average p and impact parameter dependent p model.Conversely, the low momentum behavior is different between the different TMDs, and it turns out to be important to include the full finite N c expressions in App.C for the various TMDs, to properly capture the low momentum behavior.While at low momentum k t /Q s (x) ≲ 1 there are visible differences between the results for the average proton and the impact parameter dependent model, the overall behavior of the two model TMDs is still quite similar.

Model predictions for the proton F 2 structure function measured at HERA
Before we proceed to the calculation of nuclear TMDs, we perform a closure test on the b-dependent model by comparing its predictions for the DIS proton structure function F 2 to the HERA data used in the AAMQS fit.The structure function F 2 is given in terms of the dipole amplitude N as where Ψ L,T (Q 2 , r, z) are the longitudinal (L) and transverse (T) photon wave functions.Their explicit expressions can be found at Ref. [60].Thus, we compare x = 0.003 x = 0.003 gg F (4)   gg F (5)   gg where N and N * are coming from the AAMQS fits (see Eq. 2.4) and the model dipole amplitudes (see Eq. 2.6) respectively.
The fit quality depicted in Fig. 2.4, which is comparable using both parametrizations (N and N * ) supports the validity of our model for the b-dependence.n values starting from Q 2 n=0 = 0.85 GeV 2 , to avoid the overlapping between the different curves.We note a good agreement with data using both approaches, which validates our model for the b dependence.

Model predictions for nuclear TMDs
Now that we have established the procedure for calculating small-x TMDs, we proceed to the calculation of the nuclear TMDs.We note that in order to perform an efficient computation of the TMDs, which as shown in Eq. 2.10a involve numerical integrals of Bessel functions that demand a large number of integrand evaluations, it is convenient to tabulate the TMDs as a function of the nuclear thickness value T .One example of these unintegrated (w.r.t.b) TMDs is shown in the left panel of Fig. 2.5, where we present the unintegrated dipole TMD F (x, T = T (b)) as a function of the rescaled momentum √ σT denotes the saturation scale for different nuclear thicknesses, which is depicted in the right panel of Fig. 2.5 as a function of x for three values of T .While at larger values of x, the dependence on the nuclear thickness T can almost entirely be accounted for by use of the scaling variable k t /Q s (x, T ), a more non-trivial T dependence emerges at smaller values x following the rcBK evolution.
Knowing the un-integrated TMDs F (x, T = T (b)), the integral over impact parameters in Eq. 2.10a can then simply be carried out by mapping the value of b to the corresponding unin-tegrated TMD F (x, T = T (b)).The TMDs shown in all other figures that will be discussed in the following have already been integrated over b.Once again, we first present results for the dipole TMD, which is shown in the top left panel of Fig. 2.6 for proton (p), Oxygen (O 16 ), Xenon (Xe 129 ) and Lead (P b 208 ) nuclei.Solid bands in Fig. 2.6 show the results obtained from the impact parameter dependent model, while dashes lines correspond to the naive scalings of F ∝ A 2/3 and k 2 t ∝ A 1/3 with the atomic mass number A. Generally, we find that the uncertainties of the impact parameter dependent model appear to be well under control, and the results show the expected behavior that the dipole TMD of a nucleus is larger due to the increase of the transverse area and peaks at larger values of the transverse momentum k t due to the increase of the saturation scale.However, the naive scaling with atomic number does not properly account for this dependence, as can be seen from the sizeable deviations from the impact parameter dependent predictions for all nuclei.We further elaborate on this behavior in the top right and bottom panel Fig. 2.6, where we present the A dependence of the saturation scale at three different values of x and of the nuclear TMD rescaling ratio where estimates the effective transverse areas of protons (S p ⊥ ) or nuclei (S A ⊥ ) in the model.The A A and T p values for each nucleus are computed analogously to A p and T p of Eq. 2.13 and are reported in Table 2.
When considering the results in Fig. 2.6, we find that as a function of the atomic number the saturation momentum Q A s (x) 2 increase more slowly than the naive A 1/3 scaling, whereas the effective transverse area increases more rapidly than the naive A 2/3 power law 7 .Intuitively, this can be understood from the fact that a large part of the cross section is dominated by scattering off the relatively large edges of the nucleus, where the nuclear matter densities do not significantly exceed those inside a proton, such that the nuclear enhancement of saturation effects is rather mild.We thus conclude that a proper treatment of the impact parameter dependence is important in order to make accurate saturation physics predictions for nuclei.Despite the relatively weak nuclear enhancement of saturation effects, it is important to point out, that the TMDs in our impact parameter dependent saturation model exhibit geometric scaling to a rather good accuracy.This phenomenon can be observed in the plot in the bottom left panel of Fig. 2.7 which shows the dipole TMDs as absolute values and as ratios to the proton case.Similar scaling behavior holds also for all other TMDs as shown in the top panels of Fig. 2.7, where we present results for the seven TMDs relevant to forward di-jet production for various nuclei.
Importantly, in all of the plots the nucleus TMDs are rescaled along the x and y axes according to such that the result for the various TMDs in Fig. 2.6 nearly collapse onto a single curve for all the different nuclei.While this clearly indicates the property of geometric scaling, it is important to point out once again, that the nuclear enhancement of saturation effects is rather mild, and merely gives rise to about a factor of Q P b s /Q p s ∼ 2 for the Lead, as can be seen in the bottom panel of Fig. 2.7, where we show the x dependence of the nuclear saturation scale (  x = 0.01 [GeV] 3 Di-jet production at the LHC This section is dedicated to study how the model predictions for the impact parameter dependent TMDs, detailed in the previous section, may affect the search for saturation effects at hadron x = 0.003 . By construction, the model prediction for the proton saturation scale matches exactly the "average" proton case colliders.To that end, we first introduce the selected experimental observables and then we discuss the results with emphasis on the distinctive features that separate the model predictions from what a simpler naive scaling approximation would have suggested. The forward di-jet production in proton-proton and proton-nucleus collisions is one of the preferred processes [61] through which the target (proton or nucleus) can be probed at very small longitudinal fractions, i.e. when the gluon density is large and recombination mechanisms related to saturation are expected to appear especially in heavy nuclei.We consider the special kinematics in hadronic (p − p or p − A) collisions when both partons (jets) are measured at large rapidities in the same forward hemisphere.The strongly asymmetric configuration in jet emission results into a dilute-dense asymmetry in the partonic distributions of target and probe that interact in the scattering, which is due to the strong imbalance between the longitudinal momentum factions (x 2 ≪ x 1 ) of the respective interacting partons.
At high energies, the longitudinal momentum of the jets is entirely provided by the co-moving parent parton carrying sizable part of the proton momentum 8 which can, thus, be described in terms of the well known PDFs.On the other hand, the partonic content in the target (the counter-moving hadron) is probed at much smaller x values where effects of gluon saturation can be probed.
Besides for the longitudinal fractions x also the initial transverse momentum of the scattered gluon in the target k t plays a crucial role in the sensitivity to saturation effects.Its value can discriminate between different physical regimes that require descriptions withing different frameworks.On one hand, when it is of the order of the jet transverse momenta k t ≃ p j the usual high energy factorization applies [46].On the other hand, when k t approaches the saturation scale Q s in the target for which p j ≫ Q s ∼ k t , the process can be described in the framework of Transverse Momentum Dependent gluon distributions (TMDs) supplied with small-x corrections.
Quite recently, the Improved-TMD (ITMD) factorization scheme -a convenient unifying framework that encompasses the two regimes and provides a smooth interpolation between the two descriptions -was proposed in Ref. [25].In this scheme, cross sections in the dilute-dense regimes are given as convolution of the probe (proton) PDF, the target (proton or nucleus) TMDs and several hard factors regardless of the relative size of k t and p j (provided that p j ≫ Q s ).Each TMD couples with a specific hard factor which represents the scattering off-shell matrix element.The universality of the partonic distribution on the target side is restored at large k t where all the TMDs collapse to the same perturbative power law behavior, which appears in the high-energy factorization framework 9 In Ref. [17] this framework was applied to describe forward di-jet production in proton-proton (p-p) and proton-lead (p-Pb) collisions at center of mass energies √ s N N = 8.16 TeV at the LHC.We follow this path and proceed to employ the ITMD framework in our analysis.The cross section (Eq.3.3) and the hard factor (App. B) formulae that we use are identical to those of Ref. [17].However, in contrast to previous works 10 , we employ the TMDs calculated in Sec.2.3 and App.C within our impact parameter dependent description of the nucleus, rather than assuming "naive" scaling for the saturation scale Q 2 s ∼ A 1/3 and transverse area S ⊥ ∼ A 2/3 .Consequently, important differences can be expected when considering the cross-section in p − A collisions.

Determination of the forward jet production cross section
We will compute the cross section of inclusive di-jet production in proton-hadron collisions when both jets are emitted in the forward region 11p(p 1 ) + A(p 2 ) → j 1 (p j 1 ) + j 2 (p j 2 ) + X , where X stands for anything (or possibly nothing) produced in addition to the two jets.In this process, a projectile -a proton with momentum p 1 -scatters against a target -another proton (A = p) or a nucleus (A = {O, Xe, P b}) with momentum p 2 -to produce two jets in the forward direction of the projectile (y j 1 , y j 2 > 0).Fig. 3.1 displays a schematic of the considered process.We note that both jets are required to be "hard", i.e. |p j | ≫ Λ QCD .Within our analysis, we will consider two sets of generic kinematic ranges within the acceptances of a typical detector at the LHC such as the CMS or ATLAS central detector (CTR) or the CMS CASTOR [62] and the ALICE FOCAL [63] (FWD) detectors (that allow to reach lower values of x) at a typical center-of-mass energy of √ s = 8.16 TeV.These are shown in Table 3. Experimental data are already available for a forward jet measurement using the CMS CASTOR detector.The acceptance of the forward detector in rapidity and transverse energy makes it the ideal settings to probe the lowest x tail of the target nuclear distribution.
The longitudinal momentum fraction (with respect to the parent) of the interacting parton in the projectile x 1 and in the target x 2 are defined as ) where y j 1 /2 and p j 1/2 are the jet rapidities and transverse momenta respectively.With the acceptances given in Table 3, the lowest possible values of x 2 allow probing the kinematical domain where saturation effects can appear .2: Schematic of the forward jet process In the TMD formalism, the interacting gluon in the target can be connected to the upper parton from the projectile as a gluon-parton scattering everywhere within the gray area.
The "master-formula" for the calculation of the cross section in the I-TMD factorization framework reads where J = (p j , y j ), d 3 J = d 2 p j dy j is a compact notation to indicate the jet kinematics.P t = |P| is the jet hard scale given in Eq.B.3. k t = |k| is the jet transverse momentum imbalance k = p j 1 +p j 2 that can be traced back to the transverse momentum of the small x gluon inside the target (as the other partons are considered to be collinear to the beam in this approximation) 12 .Also, we impose p j 1 > p j 2 in the phase-space integrals.The Φ of the target TMDs that couple with the hard factors ag * →cd computed in the ITMD formalism.Their explicit expressions are given in Table 5 of Appendix B. The large x partons in the projectile are described by the usual collinear proton PDFs f a (x 1 , µ F ).For the collinear PDFs, we use the CTEQ18 [64] parametrization.The number of active flavors is set to N f = 3.We exclude the heavy flavors to match that setting in the TMD calculation (see.Sec.2.1).The same CTEQ18 package is used to compute the strong coupling α s .The factorization and renormalization scales are defined as µ R = µ F = (p j 1 + p j 2 )/2.Note that hadronization effects and jet recombination mechanisms are entirely neglected in this analysis.In this approximation, a one-to-one correspondence between a given jet and its partonic parent can be drawn.

Results
In this section, we present the model predictions for the forward di-jet production in p−p and p−A collisions at √ s = 8. 16 TeV.The goal is to observe the effects on the cross section of the dynamics of emission and recombination shaping the nuclear distribution in the target, thus probing the onset of the saturation regime.A simple observable that is less sensitive to the precise measurement of the jet transverse momentum, difficult in the very forward region, is the jet azimuthal decorrelation where the right hand side of this equation was given in Eq. 3.3.Another observable is the nuclear modification factor which weights the saturation effects (normalized by number of nucleons) in nuclei versus proton targets.The larger saturation scale in nuclei is expected to result in a suppression towards ∆φ j 1,2 ∼ π, where the transverse momentum imbalance k t can be on the order of the saturation energy Differential cross sections as a function of the azimuthal angle between the two jets defined in Eq.  3).Different solid curves in Figs.3.3 and 3.4, show the results obtained using our impact parameter dependent model TMDs, whereas dashed lines show the results for the same nuclei, obtained using in Eq. 3.4 the average proton TMDs rescalend as In the case of the impact parameter dependent model, the colored bands in Figs.3.3 and 3.4 indicate the uncertainties of the cross section due to the variation of the large distance cut-off T cut (see Sec. 2.4).Before we comment on the comparison between the two, we first explain the general features of the cross-sections in Figs.3.3 and 3.4.Indeed, the interpretation of the observed distribution behavior in the context of saturation dynamics is straightforward if one keeps in mind the connection between the investigated jet configurations and the two variables x 2 and k t that describe the target gluon distributions.Firstly, the typical value of x at which the gluons in the target are resolved in the scattering depends on the jet rapidity and momentum ranges (see Eq. 3.1b) and dictates the magnitude of the saturation scale Q s .In particular, looking at Figs. 3.3 and 3.4, in this order, and at the plots within each figure from right to left we observe a decrease in the typical x value and, thus, an increase in Q s (x) value at which the target is probed.Secondly, looking at these plots in the same sequential order, we expect that, due to decreasing jet energies, the azimuthal region sensitive to saturation -where the jet momentum imbalance is of the order of the saturation scale k t ∼ Q s -to slowly expand from the vicinity of the back-to-back region ∆φ j 1,2 = π towards smaller angles.The consequence of these two reinforcing effects are the observed flatter distributions with shallower peaks at ∆φ j 1,2 = π when more forward and softer jets are measured.Similarly, the proton-proton cross sections are more peaked around the back-to-back configuration compared to the proton-nucleus case, as expected since the saturation scale is smaller for the proton as compared to the heavier ions.Even though these general features of the cross section can be observed for both the nuclear model TMDs (solid) and the naive scaling ansatz (dashed), there are discernible differences in the predictions for the cross section in proton-nucleus collisions.While for the case of the proton, the results for average proton TMDs fall within the uncertainty band obtained using the b-dependent proton model TMDs, differences between the "naive" and impact parameter dependent description, regarding both the normalization and shape of the cross section increase with increasing atomic number from Oxygen to Xenon and Lead.Notably, the largest discrepancies arise for low momentum jets and in the forward rapidity region, where the di-jet production is sensitive to small x in the nucleus, and thus to the dynamics around the nuclear saturation scale Q s .However, important unaccounted corrections are expected to substantially affect the corresponding predictions since the underlying hypothesis p j ≫ Q s , laying the foundation of the TMD description, is no longer valid.In this regime, a complete description within the full-complexity of the CGC framework would be required to take these corrections into account, and we refer the interested reader to [19] for recent progress in this direction.Conversely, for high momentum jets that mostly probe the perturbative tail of the TMDs, the agreement between the impact parameter dependent model and the naive scaling formula is rather good.We will continue on the comparison between naive versus model predictions further down after commenting on the model predictions for the nuclear modification factors.
Evidently, the precise modeling of nuclear TMDs also has immediate consequences for the behavior of the nuclear modification factor R pA .In Figs.3.5 and 3.6 we show the nuclear modification factors for Oxygen, Xenon and Lead, in the two different kinematical ranges accessible within present and future LHC experiments.Once again, the colored bands represent the model uncertainties, associated with the ±50% variation of the large distance cut-off in the nuclear TMD calculation 13 .Interestingly, we find that the modification factors R pXe and R pP b are almost identical.On the other hand, R pO remains around 1 within uncertainties except in the immediate vicinity of ∆φ = π as it is best seen in Fig. 3.6.Generally, the R pA shows a modest suppression of back-to-back jets and enhancement of nearly back-to-back jets, which can be attributed to the azimuthal de-correlation of the jet pair due to the transverse momentum k t provided by the nuclear small x TMD gluon.Since the typical transverse momentum k t is of the order of the saturation scale Q s (x), we use arrows (indicated by matching colors) to mark the azimuthal angle below which the respective nuclear TMDs are sampled at values k t /Q A s ≤ 3 in the scattering. 14Naturally, we find that the location of the maxima coincides rather well with the position of the arrows, indicating that at the respective angular scale the measurement is indeed sensitive to saturation effects.Now that we have established the basic features of the nuclear modification of di-jet production, we will contrast our results from the impact parameter dependent TMD model with those obtained using the "naive" scaling approximation, where the nuclear TMDs are obtained by scaling the saturation scale of the nucleus ( We exemplify this for the Lead nuclear modification factor shown in Figs.3.7 and 3.8, where the nuclear modification factor in the "naive" scaling approximation (brown line) is calculated as in Eq. 3.5 using the TMDs defined in Eq. 3.6.By comparing the results, one observes that in the saturation region, i.e. at low momentum, forward rapidity or close to the back to back limit, the model predictions differ significantly from the "naive" scaling expectation.Conversely, for large momentum imbalance and at large jet momentum, where one probes the perturbative tail of the gluon distributions, the model predicition is in reasonable agreement with the naive scaling prediction.Clearly, the different behavior of the two models can be directly understood from the behavior of the underlying nuclear TMDs in the top left panel of Fig. 2.6.
While the perturbative tails match between the model and the naive scaling predicition, the impact parameter dependent model predicts broader and higher peaks of the nuclear TMDs (solid lines) around the saturation scale as compared to the naive scaling (dashed lines) expectations.Noteably, the larger ratio of the nuclear TMDs to proton TMDs, also shown in Fig. 2.6, is directly reflected in the ratio of the cross-sections R pA , and gives rise to a "knee" shaped enhancement of the nuclear modification factor at intermediate angles, which is a distinctive feature of the model predictions that is not present with the naive scaling of the saturation scale of the nucleus.Generally, we find that the suppression (R pA < 1) close to back-to-back kinematics (∆φ ≃ pi) appears to be a robust model independent feature, which can be attributed to the larger transverse momentum acquired from multiple scattering in the nucleus.However, our analysis also shows that the behaviour at intermediate angles can in fact be quite sensitive to the behaviour of nuclear TMDs, and can vary from enhancement (R pA > 1) obtained in the presumably more realistic impact parameter dependent model, to suppression (R pA < 1) obtained with naive scaling of the saturation scale. 13The upper edge (+) of the error bands is defined as A dσ p (T + cut ) d∆φj 12 .
The lower edge (−) is defined likewise. 14To be more precise, the azimuthal angles marked by the arrows are found by solving the implicit equation min(kt/Q A s ) = 3, where the left hand side is a function of ∆φj 1,2 .

Conclusions
We performed a phenomenological analysis of forward dijet production in p − p and p − A (Oxygen, Xenon and Lead nuclei) collisions in dilute-dense kinematics which can be understood within the ITMD framework.The goal is to analyze the effects on the dynamics of saturation due to the complex transverse distribution of color charges in nuclei.
We developed a simple model to compute the TMDs for any nucleus given its nuclear matter distribution.By assuming Gaussian statistics for the nuclear color charge distributions, all the small-x TMDs can be computed from the knowledge of the fundamental nuclear dipole amplitude N (x, r), where r is the dipole size, which can be inferred with great precision and down to very small values of x from DIS measurements at HERA.One typical example, which we make use of, is the AAMQS fit using rcBK evolution equation.
Since saturation model fits to HERA data are often performed at the level of inclusive crosssection, and thus not sensitive to the impact parameter dependence of saturation effects, extending the calculations to the complex geometry of heavy and light nuclei is non-trivial and can result in a significant uncertainty of the predictions.Inspired by the IPSat model, we therfore propose to re-introduce the impact parameter dependence at the level of the dipole amplitude, to consistently derive TMDs of heavy and light nuclei from the underlying proton distributions.This procedure introduces a few free parameters that we omit here, which can be fixed by matching (the transverse integral of) the b-dependent TMDs with the b averaged TMDs computed from the inclusive dipole amplitude provided by the underlying AAMQS fit.As cross-check, we verified that this procedure can still reproduce the proton structure function F 2 well.Clearly the advantage of the impact parameter dependent distributions is that they can be straightforwardly generalized to nuclei, and the respective TMDs can be easily obtained by finally integrating over the impact parameter.Based on this impact parameter dependent treatment of saturation effects, the newly found TMDs show a dependence upon the nuclear mass number that deviate significantly from a naive scaling expectation.More precisely, the dependence on the atomic number A is stronger then A 2 /3 for the overall magnitude of the TMD but it is weaker then A 1 /6 for the characteristic momentum Q A s .These effects also have important consequences on the p − A scattering cross sections.
The forward di-jet cross section in p − p and p − A collisions are constructed convoluting (the impact parameter average of) the nuclear TMDs with the appropriate hard factors in the usual ITMD scheme.We presented predictions corresponding to the kinematics acceptance of the forward portion of the CMS detector and the Castor (CMS)/Focal (ALICE) very forward calorimeters.
As expected the p − A and even more so p − p cross sections are strongly peaked at ∆φ j 12 = π corresponding to small gluon transverse momenta k t in the target.The peaks are more pronounced for the highest jet momentum and more central rapidity acceptances since in this case the azimuthal angular region sensitive to saturation is squeezed towards the back-to-back configuration.Furthermore, we observe that the p − A cross sections and the nuclear modification factors deviate significantly from the naive scaling approximation, especially when the jets are softer and more forward.We confirm that, the earlier onset of saturation in nuclei versus protons induces a suppression on the nuclear modification factors R pA towards ∆φ ∼ π.However, the model predictions deviate significantly from what the naive scaling would suggest.A "knee" region develops while moving from small to large azimuthal angles, before the saturation regime is reached and the nuclear modification factors are suppressed.At intermediate angles R pA grows and can exceed unity (for Xe and P b, less so for O) before starting to drop.This effect is due to the broader and higher peaks that the impact parameter dependent TMDs show with respect to naive expectations.We therefore conclude, that a proper treatment of the impact parameter dependence of the nuclear matter distribution inside nuclei is in fact rather important in order to provide accurate predictions for saturation effects in nuclei.
Beyond the inclusive predictions presented here, it is also important to stress, that a proper treatment of the impact parameter dependence is crucial to also provide predictions for forwardcentral correlations, where e.g. the de-correlation of forward di-jets in p − A collisions can be investigated as a function of the central event activity to enhance saturation signals.Even though such calculations will require further modeling of the event activity in the central detector region, our calculations clearly provide a first important step in this direction.
We finally remark, that in addition to a proper treatment of the nuclear geometry, there are other important phenomena that were not take into account in our analysis but can be expected to affect the predictions appreciably.First of all, an important unaccounted contribution is the already mentioned correction due to Sudakov resummation of log(p j /k t ) enhanced diagrams.This resummation can be obtained in parallel with the BK log(1/x) resummation [65], and its effects in this context have recently been investigated in Refs.[20,66].We also did not include the heavy flavor (c and b) quarks in the present TMD calculations, and we intend to return to both of these points in a forthcoming work.
We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) under grant CRC-TR 211 "Strong-interaction matter under extreme conditions" project no.315477589-TRR 211.We also thank Cyrille Marquet and Martin Hentschinski for useful discussions at the beginning of this study.

A Nuclear Thickness in the MC Glauber Model
Below we explain how the nuclear thickness T A ,which describes the nuclear matter density projected in the transverse impact parameter space, is determined within the MC Glauber Model [49].Starting point is the distribution of nucleons in side a nucleus, which is generated stochastically, event-by-event, by a Glauber Monte-Carlo sampling where, for each nucleon, random values in the range −5R, 5R are assigned to each of its 3-dimensional coordinates (x = {x, y, z}) in a system of reference with its origin placed at the nucleus center of mass.Subsequently, the probability of having a nucleon at such x position is weighted by the Wood-Saxon (WS) distribution where R is the nuclear radius and a is the "surface diffusiveness" which characterizes the steepness of the distribution drop out in the proximity of its edges.The normalization constant in the Monte Carlo sampling is given by ρ 0 = 1 + exp(−R/a), and corresponds to the nucleon density in the center of the nucleus.Numerical values of the Wood-Saxon parameters for all nuclei are collected in Tab. 4.
Based on the nucleon positions, the nuclear thickness T A , is then found by convoluting the nucleon distribution in nucleus with the transverse matter density within each nucleon, i.e. given the nucleon distribution at the n-th MC iteration where x i with i = 1 . . .A specify all the nucleon positions, the nuclear thickness is given by the superposition of the corresponding proton (or neutron) distributions where b i are the nucleon positions underlinex i projected over the transverse space.We employ a Gaussian profile with a width of √ B G = 4GeV −1 .for the transverse profile T p/n of the nucleon, with no distinction between protons and neutrons (T p = T n ), with

B Hard Factors
All the expression collected in this section were taken from Ref. [25].They are reported here for convenience of consultation.The gluon TMDs Φ (i) are given by a linear combination of the actual TMDs xx introduced in Sec. C. They are collected in Table 5 together with the corresponding hard factors K (i) which are given in terms of the Maldestan invariants Eqs. B.1a sum up to ŝ + t + û = k 2 t and s + t + ū = 0. i 1 2 gg + F (6)   gg − 2F (3)  gg + F (4)  gg + F (5)   gg gg − F gg + F qg − F ag * →cd and the gluon TMDs Φ ag→cd in the finite-N c limit.

C TMD calculation
In this appendix we describe the calculation techniques that have been employed in order to compute the gluon TMDs needed for the analysis which were briefly introduced in Sec.2.3.Starting point of our discussion is the operator definition of the TMDs in the small-x limit, where the relevant TMDs for di-jet production take the form Ref. [51] for the qg channel and for the gg channel where S ⊥ denotes the transverse area of the target.
We have calculated the TMDs using an ensemble of lattice Wilson lines (see Sec. C.2) as well as in the Gaussian approximation (see Sec. C.1).Note that in the large N c limit, expressions for the various TMDs have been provided in the literature, and the structure is simple especially in coordinate space (Eqs.C.9a,C.10a).In the qg channel, they are given by Gluon TMD: g Figure C.1: Gluon TMDs in the MV model.Data points correspond to numerical lattice simulation.Solid (dashed) lines show the result obtained from the Gaussian appproximation at finite N c (in the large N c limit).
while, in the gg channel, they are while in the gg channel the various TMDs are given by Even though in the large N c limit, the momentum space expressions can be compactly expressed in terms of convolution of various distributions (see e.g.[51]), this is no longer the case at finite N c and the corresponding momentum space expressions have to be obtained by performing the Fourier transformation numerically according to15 In Fig. C.1 we present a comparison of numerical lattice results and the Gaussian approximation at large and finite N c = 3 for the MV γ model, where the Gaussian approximation is exact.Indeed we observe excellent agreement between the different calculations, providing an important cross-check for our calculation.Interestingly, one also observes from this comparison that finite N c corrections appear to be important for |k| < Q s .

C.1.1 Expressions for gluon TMDs in terms of derivatives of local multipole operators
We relate the position space expressions of the various TMDs to the fundamental Dipole, Quadrupole, Sextupole and Octupole operators such that at the level of a single dipoles (n = 1) one obtains the dipole gluon distribution F (1) Similarly, at the level of two dipoles or a quadrupole (n = 2) one obtains F gg , F gg , F gg , F gg as 1 2) (x 1 , y 1 )S (2) (x 2 , y 2 ) x 1 =x,y 1 =y x 2 =y,y 2 =x , (C.17) qg is related to the product of a quadrupole and a dipole (n = 3) 1 and F (5 gg are related to the octupol and quadrupole-dipole-dipole operators involving eight Wilson lines (n = 4) where in some of the above expressions we used the symmetry V x i ↔ V † x i to unify the operator structures.

C.1.2 Calculating local multipole operators
We follow the standard procedure for evaluating the expectation values of the correlation functions of n Wilson lines, by reformulating the problem to a coupled set of evolution equations of various multipole operators.Since the color fields are assumed to follow local Gaussian statistics, the averaging procedure is equivalent to solving the stochastic differential equation from z = 0 to z = 1 with V x (z = 0) = I and Gaussian noise ξ a x (z) such that Evolution equations for arbitrary products of 2n Wilson lines can then be obtained in a straightforward fashion by applying the rules of stochastic calculus.Considering e.g. the simplest possible example of a single dipole (n = 1) one finds such that upon evaluating the noise expectation value according to Eq. (C.23) and using t a t a = C F to simplify the color structure, one obtains where we defined the correlation function which frequently appears in color singlet operators.Since we will neglect the impact parameter dependence in the target, the function G xy is assumed to be a function of only |x − y|, such that while the second derivative is given by which will be relevant for evaluating the TMDs.We also note that by virtue of the definitions the correlation function satisfies the following symmetry properties We illustrate the procedure at the example of F (1) qg (x, y), which according to Eq. C.15 is related to the fundamental dipole S (2) (x 1 , y 1 ).By using the prescription in Eq. (C.22) and using the rules of stochastic calculus, it is straightforward to derive the evolution equation for S (2) (x 1 , y 1 ) as which is trivially solved by S (2) (x 1 , y 1 |z) = e zGx 1 y 1 .By evaluating the expression at z = 1 and taking the appropriate derivatives, one then immediately obtains However, from the point of view of evaluating higher order correlation functions n ≥ 2 the exponentiation of the operator on the r.h.s of the evolution equation can be difficult and it is more convenient to evaluate the derivatives prior to exponentiation, by formally expressing the gluon distribution as which obviously gives the same results upon integration.

C.1.4 Evolution of n = 2 operators
Clearly, this strategy becomes advantageous already in the case of n = 2 relevant for the calculation of F gg , F gg , F gg , F gg .By following the procedure outlined above, one obtains a coupled set of evolution equations for the dipole-dipole and quadrupole operators where the transition matrix element takes the form Even though it is still possible in this case to perform the exponentiation of the operator on the r.h.s.(see e.g.[XXX]) it is much easier to first take the derivatives and set the coordinates to the desired values which is easily exponentiated.By taking the various derivatives of the evolution operator on the r.h.s of Eq. (C.36), the gluon TMDs F gg , F gg can then be evaluated as in Eq. (C.35) yielding Gxy . (C.40) Similarly one also finds the well known result for the Weiszaecker-Williams distribution F gg , F [colorredXXX], which in our notation takes the form one has a number of different possibilities for the insertion of the Gaussian noise ξ a t a , which upon performing the Gaussian averaging schematically correspond to a) the one possible insertion of generators between neighboring Wilson lines at position q of operator number p which leaves the operator structure intact.Or b) the four possible insertions of generators between non-neighboring Wilson lines at position q and l > q of the same operator number p By using the Fierz identity t a ij t a kl = 1 2 δ il δ kj − 1 2Nc δ ij δ kl , the expression can be further simplified and separated into two distinct contributions.While the first contribution generates an additional trace, the second term proportional to the identity does not change the original operator structure.

By adding the two contributions one obtains
where we defined the transition matrix element Of course, there are also c) the four possible insertions of generators between Wilson lines at position q and l of the different operators numbered p and k > p which following the same steps as indicated above yields Even though the general expressions are rather lengthy, and the number of distinct operators grows rapidly as a function of n the procedure can be automated in a straightforward way and is easily tractable for the n-values of interest.One first constructs a basis of all possible configurations.Starting from a configuration of n fundamental dipoles n i=1 S (2) (x i , y i ), one considers all possible generator insertions a),b) and c) and add the new operator structures to the basis.By repeating the second step until no new operators are generated (n − 1 times) one obtains a complete basis for the evolution of the 2n Wilson line correlators.Next one constructs the elements of the evolution matrix, by considering all possible generator insertions a),b) and c) on each basis state.By identifying the relevant operators for the TMD calculation, one can then take the relevant derivatives and subsequently set the coordinates (x 1 , y 1 , • • • , x n , y n ) to x, y as required for the evaluation of the operator.By performing the exponentiation as in Eq. (C.35) and projecting the result on the appropriate operator configuration one obtains the expression for the TMD.Since the expressions involved can be rather lengthy, we refrain from providing further details on the intermediate steps of the calculation and simply quote the final results.

C.2 Calculation of Gluon TMD's on the lattice
Below we explain the calculation of the various gluon TMDs from the Wilson line configurations V x discretized on a 2-D transverse lattice, as it is standard for studying e.g.JIMWLK evolution.We discretize the light-cone fields E (+/−) µ,x .By virtue of the SU (N c ) Fierz identity t a ij t a kl = 1 2 δ il δ kj − 1 2Nc δ ij δ kl one then has (C.54b) We note that based on the unitarity relation we can then express all derivative of Wilson lines according to such that all relevant correlation functions can be expressed solely in terms of products of light-like Wilson lines V x and light-cone fields E (+/−) µ,x . Based on the above defintions of lattice operators, the small-x TMDs in the qg channel can be compactly expressed as  to express Based on this procedure, the small-x gluon TMDs in the qg channel are then given by g 2 (2π) 3 Based on this procedure, the small-x gluon TMDs in the gg channel are then given by such that all the TMDs can be calculated as local operators in Fourier (k) space.

Figure 2 . 1 :
Figure 2.1: Schematic of DIS in the dipole picture.

Figure 2 . 3 :
Figure 2.3: Different TMDs (up to a constant c = Nc αsπ(2π) 3 ) entering the di-jet production calculation for the proton (p) for three different values of x.Dashed lines show the b-averaged proton cases (p), while solid lines correspond to the impact parameter dependent model.The colored bands indicate the model uncertainty range.

Figure 2 . 4 :
Figure 2.4: F 2 (x, Q 2 ) (dashed line) and ⟨F 2 (x, Q 2 , b)⟩ (full line) compared to HERA data.The plots corresponding to different photon virtualities Q 2 = Q 2 n are shifted along the y-axis by a factor of 2 n , where n = 0 to 16 enumerates the Q 2n values starting from Q 2 n=0 = 0.85 GeV 2 , to avoid the overlapping between the different curves.We note a good agreement with data using both approaches, which validates our model for the b dependence.

A=
A A (T −=+ cut ) and average thickness value T −=+ A = T A (T −=+ cut ) where the superscripts −, =, + indicate the lower, central and upper T cut values respectively.

Figure 2 . 5 :
Figure 2.5: Left: (Unintegrated) F (1) qg TMDs plotted as a function of the (rescaled) transverse momentum kt = k t /Q s (x, T ) for small, intermediate and large values of the nuclear thicknesses for four different values of x.The line thickness indicates the uncertainty due to the choice of the cut-off value.Right: Saturation scale for three different values of T as a function of x.

Figure 2 . 7 :
Figure 2.7: Top: All 7 independent (rescaled) TMDs (F (3) gg = F (4) gg ) (up to a constant c = Nc αsπ(2π) 3 ) at x = 0.001 for b-dependent (p) and b-independent (p) proton, Oxygen (O), Xenon (Xe) and Lead (P b).The colored bands represent the uncertainty due to the cut-off variations.Bottom Left: Model predictions for proton, Oxygen, Xenon and Lead dipole TMDs F (1) qg (up to a constant c = Nc αsπ(2π) 3 ).The bottom plot shows the result for different ions relative to the proton.At intermediate to large transverse momenta the (rescaled) TMDs match exactly whereas at smaller values the ratio can (slightly) exceeds values of 1.5 at about k t /Q A s ≲ 0.1.Bottom Right: Nucleus squared saturation scale as a function of x.The saturation scale (squared) from our model (shown in full line) grows slower than the "naive" Q A s (x)2 ∝ Q 2 s (x)A 1/3 expectation (shown in dashed line).By construction, the model prediction for the proton saturation scale matches exactly the "average" proton case Q p s ≡ σT p Q s = Q s .

Figure 3 . 1 :
Figure 3.1: Schematic of central and forward di-jet hadro-production.Only the forward jet configuration can be sensitive to the saturation region x ≪ 1, k t ∼ Q s .
3.4 for p − p, p − O, p − Xe and p − P b are shown in Figs.3.3 and 3.4 in the CTR and FWD detector kinematics respectively (see Table

Figure 3 . 3 :
Figure 3.3: Azimuthal angle distributions for p, p, Xe, O and P b for three bins in jet p j in the acceptance of the central detectors as predicted by the model (solid bands) and in the "naive" scaling approximation (dashed lines).The p (no b-dependence) lays within the (b-dependent) p error band.

Figure 3 . 4 :Figure 3 . 5 :Figure 3 . 6 :Figure 3 . 7 :
Figure 3.4: Azimuthal angle distributions for p, p, Xe, O and P b for two bins in jet p j in the acceptance of the forward detectors as predicted by the model (solid bands) and in the "naive" scaling approximation (dashed lines).

Figure 3 . 8 :
Figure 3.8: Nuclear modification factor R pP b as a function of jet ∆Φ in two bins in p j from our model compared to the "naive" scaling prediction in the acceptance of the forward detectors.The vertical arrow gives an indication on the ∆φ * value where saturation effects become important defined as k t /Q A s < 3 on the right from the arrow.It is not visible on the left plot where k t /Q A s < 3 is everywhere.

2 c N 2 c − 1 Gxy − 1 . (C. 41 )C. 1 . 5
Evolution of general operatorsBy using the prescription in Eq. (C.22) it is straightforward to obtain the evolution equations for arbitrary products of 2n Wilson lines.Considering the most general operator consisting of 2n Wilson lines

4a 2 sF 2 (2π) 3 4a 2 sF ( 1 ) 1 N 1 N 2 sF ( 5 )
V y Tr V † x V y = 1 N c x,y e −ik(x−y) Tr E (−) µ,x E (−) µ,y Tr V † x V y , (C.57b)and similary one finds the following expressions for the small-x TMD's in the gg channelg gg (k) = 1 N c x,y e −ik(x−y) Tr V † x E V y Tr V x V † y = c x,y e −ik(x−y) Tr E c x,y e −ik(x−y) Tr V x E e −ik(x−y) Tr V x E (−) µ,x V † x V y E (−) µ,y V † y (C.58c) g 2 (2π) 3 4a gg (k) =x,y e −ik(x−y) V y Tr[V x V † y ]Tr[V y V † x ] (C.58e)Notably, the definitions in terms of E (−) µ,x are simpler to implement, as they require only a single operator structure for each channel.Specifically, by defining operators with open color indices according to XXX for the qg channelE (−) µ V † ij,klk use of the identities V kl = V † lk * and E

Table 1 :
σ and η values for different values of T cut .

Table 3 :
Kinematic cuts corresponding to the acceptances of the CMS central detector (CTR) and ALICE/FOCAL and CMS/CASTOR forward detectors (FWD).