Searching for saturation in forward dijet production at the LHC

We review recent results for forward jests at the LHC and EIC as obtained within small-x Improved Transverse Momentum Dependent factorization (ITMD). In addition to elementary overview of various approaches to perturbative QCD at high energy, including High Energy Factorization, Color Glass Condensate and ITMD, we describe the Monte Carlo implementation and discuss the existing and unpublished phenomenological results for forward dijets.


Introduction
Quantum Chromodynamics (QCD) is a well established theory that describes interactions of quarks and gluons.However, it still has its challenges.In the high energy domain, one of the long standing problems is finding clear experimental signals of gluon saturation, which is a signature of quasi equilibrium between gluon splitting and gluon fusion in dense nuclear systems.Gluon saturation has been predicted from QCD long time ago [1,2] and has been extensively studied using various approaches, most recently the Color Glass Condensate (CGC) -i.e.effective theory obtained within QCD (for a reviews see [3][4][5][6][7] and the textbook [8]).
While there is no doubt that gluon distributions must saturate at some point due to the unitarity constraints on the cross section, and there are strong indications of saturation in the data [9][10][11][12][13][14][15][16][17][18][19][20], a complete consensus on reaching it is still to be achieved.This is mainly due to the demanding kinematics.It requires the final states to be measured in the forward region of the detectors and dealing with nuclear targets, which posses additional challenges due to the rich collision environment (see e.g.[21]).Moreover, theoretical predictions in high energy QCD are at present not as precise as those requiring ordinary collinear factorization.
Among various final states that can be measured in the context of saturation searches, the system consisting of two identified jets in the forward region (and everything else) plays a special role [22][23][24].In addition to the possibility of studying correlations, it has an important advantage over the single-inclusive jet production.Namely, the jet transverse momenta can be quite large, of the order of twenty GeV or so, and still be sensitive to the saturation effects in the back-toback region.In single inclusive jet production, the jet transverse momentum is demanded to be of the order of the saturation scale Q s and determines the average transverse momenta of gluons being "saturated".The latter is of the order of a few GeV, for sufficiently high collision energy and sufficiently forward.Jets with larger transverse momenta are easier to reconstruct and are "cleaner".Moreover, on theoretical grounds, the description of harder dijets becomes possible in terms of "generalized" transverse momentum dependent (TMD) factorization, which makes it easier to implement in Monte Carlo generators than the full CGC.
Dijet yields at relatively forward rapidity have already been measured at the Large Hadron Collider (LHC) by ATLAS collaboration for both proton-proton and proton-lead collisions [25].Since there was no cross sections measurement, the conclusions regarding saturation signals do not give yet a definite answer to whether saturation has been observed or not.Inclusive single forward jet energy spectra for proton-lead collisions measured at CMS CASTOR detector [26] does not provide convincing proof either.Further research and analysis are necessary to gain a better understanding of gluon saturation and its effects on nuclear systems.In particular, more measurements of both proton-proton and proton-lead collisions in the same kinematics are needed.In addition to the ATLAS forward detector, the ALICE collaboration plans to build a more forward detector FoCal [27] that hopefully will shed more light on the saturation phenomenon.The forward jet physics at LHC is complementary to the physics of the Electron Ion Collider (EIC) [28,29] one of the primarily goals of which is to study gluon saturation physics.
In this work, we review the theoretical framework suitable for a description of forward dijet production in hadro and lepto-production in the full azimuthal angle range.The formalism lies in the intersection of the CGC theory and more traditional factorization approach, utilizing TMD gluon distributions.Since the formalism accounts for power corrections it has been dubbed as Improved TMD Factorization (ITMD) [30].A review of some essential aspects of high energy QCD and the ITMD framework itself is contained in Section 2. We shall also review Monte Carlo implementation of the ITMD framework KaTie (Section 3), construction of the TMD gluon distributions, as well as recent phenomenological predictions for azimuthal dijet correlations, both for proton-proton and proton-lead collisions, in the kinematics of the forward calorimeters of ATLAS detector and planned FoCal of ALICE (Section 4).We also include new unpublished earlier computations of the rapidity distributions for different kinematic cuts.
2 Forward dijet production in high energy QCD

Kinematics
The process in question is the inclusive production of two jet, that is, there are at least two jets with transverse momenta above a certain threshold, and both are produced in the forward rapidity interval in the center of mass frame: where A is either the proton state p or a nucleus having the four momentum P A that is hit by the proton state with four momentum P B (see Fig. 1).The produced jets have four momenta p 1 and p 2 and are defined by a proper jet algorithm.X corresponds to other particles produced in the process; there are no kinematic cuts imposed on those additional states.We shall also discuss the complementary DIS process A(P A ) + e − (P B ) → J(p 1 ) + J(p 2 ) + e − (P ′ B ) + X , which will be studied at EIC and will be essential in testing the theoretical formalism and constraining the nonperturbative input.
In what follows, we define the P A and P B four momenta as being the two light cone vectors, defining the"plus" and "minus" light cone components: with the center of mass energy s = 2P A • P B .Using these, the Sudakov decomposition of any four vector reads: where ( The central assumption in our investigations is that the nucleus (or a proton) is probed at small longitudinal momentum transfers compared to the hadron longitudinal momentum.That is, we define longitudinal momentum fractions where k 1 is the momentum transferred to the target A, whereas k 2 is the momentum transferred to the proton target.We assume These fractions can be also expressed in terms of the final state rapidities y i and transverse momenta Therefore, restricting the final states to a large rapidity window and keeping the jet transverse momenta relatively low guarantees the smallness of x A .
In phenomenological applications, we shall follow the realistic setup of the LHC experiments that measured (or plan to measure) forward dijets in proton-proton and proton-nucleus collisions: ATLAS and ALICE with their planned forward calorimeter FoCal.We shall discuss the kinematic cuts more precisely in Section 4 devoted to phenomenology.

QCD at high energy
One of the biggest achievements in theoretical developments of perturbative QCD are hard factorization theorems (see [31] for a review).These include the collinear factorization and the TMD factorization, both types proved to all orders in coupling constant α s and the leading power in the hard scale Q 2 , for few sufficiently inclusive processes.In practice, not only the truncation of perturbation series in α s is necessary (at present at very low orders), but also the limit Q 2 → ∞ is rarely adequate.Indeed, when comparing with experimental data, for example for large transverse momenta dijet production [32,33], both a normalization factor (so-called K-factor) is needed, as well as a resummation (parton shower) with the addition of semi-nonperturbative effects (like multiple interactions).
In the collinear factorization or TMD factorization, one encounters powers of logarithms of both x A and x B at every perturbative order.Therefore, at very large energies, small transverse momenta, forward production, or a mixture of all of these conditions, the perturbative expansion becomes less reliable.Although within the collinear factorization, it is in principle possible to resum the logarithms of small x (see [34] and references therein), this does not address the two important issues.First, the evolution of the PDFs is always linear, whereas at high energy QCD, one must include the nonlinear effects that will tame the growth of the cross section as required by the unitarity of the scattering matrix [35].Second, the collinear factorization neglects higher twist effects (that contribute to the power corrections).For a discussions of higher twists in inclusive DIS and Drell-Yan see e.g.[36][37][38][39][40][41].
Below, we briefly recall the existing approaches that attempt to address the above issues (for an extensive and pedagogical review of high energy QCD see [8]).We focus on the central ideas and main references, skipping the technical details.

Pomeron and The Reggeon Field Theory
Historically, Pomeron -the Regge trajectory with intercept greater than one and vacuum quantum numbers, was introduced in order to explain the growth of the total hadronic cross section with energy.The nature of this quasi-particle is essentially non-perturbative.Together with its parityodd partner, the Odderon, they constitute a possible effective theory of strong interactions with color singlet degrees of freedom.Typically, their interactions are described by the Euclidean field theory called the Regge field theory [42] (see also [43,43] and [44,45]).In its most general form, it includes the multiple Pomeron and Odderon interaction, but in the context of unitarity, one typically studies a truncation to the triple Pomeron vertex.
In perturbative QCD the parity even color singlet state exchanged in the t channel can be made out of two gluons.This leads to the perturbative Pomeron of Balitsky, Fadin, Kuraev and Lipatov (BFKL) [46][47][48][49] (for a pedagogical review see eg. [50]).The corresponding energy evolution equation of the Pomeron Green function gives the power like behavior of the cross section with energy.The triple Pomeron vertex needed for nonlinear taming of the growth can also be constructed in QCD [51][52][53].The difference in the Pomeron intercept derived in perturbative QCD and the one needed to describe for example elastic hadron-hadron scattering leads to a distinction between the hard and the soft Pomeron.

Lipatov effective action
The effective building blocks of the diagrams in the BFKL approach are the "reggeized" gluons R and effective interaction vertices RR → g.The former is characterized by the eikonal coupling to other vertices (due to the high energy approximation) and the "reggeization" factor (−s) ω(t) appearing due to radiative corrections, where ω(t) is the perturbative Regge trajectory.The RRg effective vertices are separated in rapidity (the quasi-multi-Regge kinematics).
Stripping off the reggeization factors, the building blocks can be formalized into an effective gauge invariant action [62,63] that contains both the QCD degrees of freedom, and the reggeized gluon fields A + , A − , that, roughly, by virtue of equations of motions are straight semi-infinite Wilson lines along two light-cone directions P A and P B .The vertices are naturally generalized into RRg . . .g vertices, as well as multiple reggeon vertices.For a review see [64] (chapter 11 in [65]).Using Lipatov's high energy effective action beyond tree level is rather cumbersome due to double counting between gluon field and composite reggeons, but several results have been obtained (see for example [66][67][68][69][70]).The important feature of that action is that it includes all the necessary ingredients to unitarize the cross section.In particular in [71][72][73][74] a relation to the Color Glass Condensate was established and investigated.

High Energy Operator Product Expansion
Consider a projectile, a photon for simplicity, scattering off a hadron in a frame where it is moving slowly.The photon fluctuates into a pair of quarks well before it hits the target and -at high energies -interacts with the target eikonally.In the considered high energy kinematics, the target is made of gluons, that are treated as background field, which shrinks to a shock wave.The quark-antiquark dipole interacts with that shockwave becoming a straight infinite Wilson lines of the background field.Since the scattering cross section is related to the hadronic matrix element of the time ordered product of quark current, effectively, the above picture provides a means of decomposing the hadronic tensor into products of Wilson lines and the so-called impact factors.This is the essence of the high energy operator product expansion [75].More precisely where is the straight infinite Wilson line along the plus light cone direction and fixed transverse position, whereas I µν 2 is the leading order impact factor.The dots represent the subleading corrections.Similar to the treatment of the ordinary factorization, where the hadronic matrix element is traded to the partonic one when deriving the evolution equation, here it is traded to the matrix element of the background gluon field The distinction between the background field and the projectile is related to the rapidity; formally the matrix elements depend on the rapidity cutoff.This dependence is explicitly introduced when regulating the rapidity divergence appearing when considering perturbative correction.However, a perturbative gluon interacting with the shock wave generates an adjoint Wilson line that gives a product of four fundamental Wilson lines.The evolution of this four-point correlator generates in turn further higher point correlators.Therefore the resulting evolution equation is not closed and consists of a tower of products of an increasing number of Wilson lines, known as the Balitsky hierarchy.However, at the large N c limit, the hierarchy is terminated and one obtains a single evolution equation for the product of two Wilson lines, the so-called dipole amplitude where η is the implicit rapidity cutoff.The evolution equation [75] reads It was also obtained independently in [76] using the Mueller dipole approach [77] and is therefore called the Balitsky-Kovchegov (BK) equation.

Color Glass Condensate
The High Energy Operator Product expansion does provide the evolution equation but does not allow for the calculation of the non-perturbative correlators itself.T is tamed at kT ∼ Qs, with Qs increasing with energy.
On the other hand, the effective theory of high energy QCD -the Color Glass Condensate (CGC) is a convenient framework to approach both the evolution and averages of Wilson lines in the background color field (for a review see for example [5]).
It is convenient to review now a famous model of large nucleus -the McLerran-Venugopalan (MV) model [107,108].Assume we view a large nucleus with a large number of nucleons in a frame where it is moving with a very large velocity.In that boosted frame, large x partons are localized in one of the light-cone directions and can be considered as color sources for the classical Yang-Mills fields.These sources are some, apriori unknown, distributions in the transverse plane ρ a (x T ) and are considered independent within the large nucleus because there are large numbers of nucleons inside.One can solve the Yang-Mills equation for such configuration to obtain the gluon "wee" fields.In order to compute observables, one needs to average over the color sources.For example, a basic quantity to consider is the particle number density with the generic definition where Ãi b are partially Fourier transformed gauge fields, k = (k + , ⃗ k T ) is a three-vector conjugate to y = (y − , ⃗ y T ).We have assumed here that the large boost is in the light cone "plus" direction.In the MV model, the hadronic matrix element of the gluon field correlator is traded to where the averaging is over the color sources in the functional sense.W x is the functional weight for random color sources, which in the MV model is a Gaussian.Direct computation of the above gluon distribution leads to the so-called Weizsacker-Williams (WW) gluon distribution [109,110].Interestingly, despite it being the most basic quantity in the MV model it is not the easiest one to compute nor can it be probed in the simplest scattering processes.Rather, as we discuss later, the quantities that appear for example in inclusive DIS are correlators of Wilson lines.The WW gluon distribution can be approximated as where µ A is the average color charge per unit transverse area per color.The saturation phenomenon is visible in the above distribution as follows.If the two-point correlator ( 16) was an ordinary perturbative correlator it would behave like ∼ 1/k 2 T ; this behavior leads to the power-like growth of the gluon distribution.It turns out that (18) behaves like ∼ 1/k 2 T in the perturbative domain (large k T ), but at small k T it behaves like 1 αs log , where Q s is the scale at which the suppression happens, called the saturation scale (see Fig. 2).This dynamical scale depends on x (through the evolution) and increases with decreasing x.
The subscript x in the functional weight W x denotes a longitudinal cutoff between the sources and the "wee" partons.Decreasing the cutoff produces new color sources and thus the Gaussian distribution gets distorted.The CGC high energy effective theory predicts the evolution equation of W x in x.It is the so-called Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation [111][112][113][114][115][116] and has the following general form where H JIMWL is the so-called JIMWLK Hamiltonian (we skip its exact form here).It gives a nonlinear evolution of W x and in turn also a nonlinear evolution of any gauge field-dependent quantity due to the averaging procedure in ρ.These equations are consistent with the Balitsky hierarchy mentioned earlier, therefore it is often called B-JIMWLK equation.

Small-x Improved TMD (ITMD) factorization
In this section, we shall review the formalism that is suitable for phenomenological studies of forward jet production processes at moderate transverse momenta.The formalism lies at the intersection of the two aspects of QCD that in the past were rarely overlapping and, essentially, were exercised by distinct communities.One is the high energy scattering reviewed in the previous Section and the other is the hard factorization theorems, in particular the TMD factorization (see eg. [31] for a review).

Dijet production in dilute-dense collisions
We start with a CGC description of a scattering of a proton or a photon off a nuclear target to produce two partons.The basic assumption here is that, in the case of the proton projectile, it is dilute, that is the partons extracted from the projectile are moderate x partons.With that assumption, we can treat the scattering of color dipoles, such as γ → q q, q → qg, etc., off a color field of a nucleus.This is the essence of the dilute-dense (or hybrid) approach which, as explained in Section 2.1, is kinematically well suited to study forward jets [117].Consider, as an example, the scattering of a q → qg color system off a nuclear target (Fig. 3).In the shock wave approximation, the interaction with the color field of the nucleus can take place either before the quark splits or after.For example, in Fig. 3 the scattering takes place before the splitting in the amplitude and after the splitting in the amplitude conjugate.Each colored particle gains a nonabelian phase when scattering, given by the Wilson line where the generator T a R is in the fundamental representation (R = F ) for a quark or in the adjoint representation (R = A) for a gluon.The cross section for the scattering of that system can be written as [22] where ψ q (z, ⃗ x T ) is the wave function of the dipole, with z being a fraction of the incoming quark momentum taken by a gluon.This wave function can be computed in perturbative QCD.The quantities in the curly bracket correspond to the color averages of the Wilson lines: The above correlators depend on x through the B-JIMWLK equations.In practice, one often uses the Gaussian approximation [118], which states that the color source distribution stays Gaussian throughout the evolution.This allows us to compute the evolution of correlators in closed form.
In order to compute the differential cross section for dijet production in proton-nucleus scattering in the above setup, one needs to find the remaining color dipole contributions [22,119] and convolute the parton-nucleus cross section with the ordinary collinear PDFs.
The CGC formulation of the dilute-dense scattering provides the high-energy description of the jet production at both small transverse momenta p T of jets and moderate (at very large p T one should switch the approach from leading energy to leading power in the hard scale).The drawback of the above approach, in addition to being hard to generalize and implement in a Monte Carlo code, is that the perturbative information is contained not only in the color dipole wave function but also in the correlators.As we discuss later, in certain limits this information can be extracted and combined with the dipole wave function to obtain the hard matrix elements, similar to those known from hard factorization.

High Energy (or k T ) factorization
Usually, the high energy factorization (HEF) (or k T -factorization) refers to a description of particle production at asymptotically high energies in terms of "unintegrated" PDFs that undergo the BFKL [141,142] or CCFM [143,144] evolution.More precisely, within HEF the cross section for a multi-jet production in a collision of two hadrons can be written as where F g/H are unintegrated gluon distributions defined as where Φ H is the non-perturbative impact factor of a hadron and G is the BFKL Green's function.
At leading order, the partonic cross section in (25) reads where V RR→P ...P (the bar denotes the usual spin/color summation and averaging) is the suitably normalized tree-level RR → P . . .P vertex to produce n partons (jets), with R being reggeon states understood as fields from the high energy Lipatov's action (i.e.without the "reggeization" factors), while P being on-shell quarks or gluons.Above, dΓ n is ordinary on-shell phase space and F is a jet function that implements the cuts.The momenta of reggeons are Notice that they are off-shell and lack one of the light-cone components, due to the high energy approximations.
In the original work on HEF in heavy quark production [142] the partonic cross section (27) was expressed in terms of off-shell amplitude to produce heavy quark pair, with initial gluons having momenta (28) and projected onto P A and P B (eikonal coupling).In that simple case, it turns out that the amplitude calculated in terms of ordinary diagrams is gauge invariant, despite the initial gluons being off-shell.In practical applications, it proves to be very convenient to stick to that logic also for more complicated processes, instead of using Lipatov's effective action.Indeed, it is possible to define off-shell gauge invariant amplitude for arbitrary processes.In addition to ordinary Feynman diagrams (with off-shell space-like gluons), one needs "gauge restoring" contributions.Those additional diagrams can be reconstructed for example by the "embedding" method [145], where an off-shell process is embedded in an on-shell process with special kinematics of the auxiliary quarks or gluons.This method was automatized for arbitrary Standard Model tree level process in the KaTie Monte Carlo discussed in Section 31 .Also, other methods suitable for automatization were also developed [147][148][149][150].In Fig. 4 we illustrate one of the methods.
In the context of forward jets, one can apply similar considerations to the hybrid approach explained before.Namely, when, say, x B is moderate and not small, one should extract partons from collinear PDF rather than unintegrated PDF undergoing BFKL evolution.This was discussed in [151] and the corresponding HEF formula reads where we have explicitly denoted the fact that the partonic cross section is constructed from amplitude with one off-shell (gauge invariant) gluon.Above f a/B is the collinear PDFs for parton a (quark or gluon) in hadron B.
Let us now refocus our attention to the dijet case.It turns out that the k T -factorization formula with proper gauge invariant off-shell amplitudes can be retrieved from the CGC expressions discussed before [30,119].In the dilute limit, which corresponds to | ⃗ k T | ≫ Q s , one can neglect the multiple scattering off the target.This means that the triple and quadrupole operators S (3) , S (4)  can be expressed only in terms of the dipole S (2) .In that limit, one obtains for the q → qg dipole scattering where and the unintegrated gluon distribution is related to the average of the weak field limit of the CGC dipole operator as follows The hard factor obtained above turns out to be exactly the off-shell gauge invariant amplitude.

TMD factorization
Formally, the TMD factorization is the leading power (in the hard scale µ 2 ) factorization of a cross section into TMD dependent PDFs and hard factors, that to leading power are on-shell (for a review see [31] and [152,153]).This factorization does not resum the large small-x logarithms, therefore becomes unreliable at very high energies.Moreover, formally it does not hold for processes we are interested in.However, the formalism provides sturdy theoretical definitions of the TMD gluon distributions that, as it turns out, can be matched to the CGC correlators.
In TMD factorization, gluon distribution is given by the Fourier transform of the bilocal matrix element of the gluon field strength tensor where F j+ = F j+ a t a is the gluon field strength tensor; the two operators are displaced in both the light cone and transverse direction (unlike in the collinear PDF, where the displacement is only along the light cone).The bilocal operator would not be gauge invariant, therefore the proper general definition requires gauge links U C1 and U C2 (here everything is in the fundamental representation) that connect the two space-time points.There is also a possibility of the double-trace over the fundamental representation (see below).We consider here only the unpolarized case, therefore the transverse index j is summed over.The TMD gluon distributions are connected to a partonic process by virtue of factorization.Since we consider gauge theory, there are multiple soft and collinear gluons that can be connected to various places in the diagrams.The definition given above corresponds to the bare operator.In perturbation theory, it contains UV and rapidity divergences that, for some processes, can be removed order-by-order by the operator renormalization and by (part of) the soft factor accumulating soft gluons.This gives the hard scale and rapidity evolution of the TMD PDFs.The collinear gluons (collinear to the target hadron) can be resummed into the Wilson lines U C1 and U C2 .As can be easily understood, the form of these Wilson lines depends on the actual hard process (its color flow).In [154] a general procedure of determining these Wilson lines via resummation of collinear gluons was given.It turns out, that they can become quite complicated for colored partonic processes, see for example [155].In Table 1 we collect the operators relevant for the dijet production.The notation F (i) gg and F (i) gg corresponds to TMD PDFs that appear for incoming gluons (the "gg" subscript) or quark-gluon system (the "qg" subscript), and various color flows (the (i) superscript).The Wilson lines U [±] are defined as which is the past-pointing staple-like gauge link and is the future pointing "staple", see Fig. 5.The square brackets above [x, y] is a standard notation for segments of straight gauge link between the points x and y.Out of two staples it is possible to make a Wilson loop The relation of the field theoretical definitions of TMD PDFs and small-x QCD has been a subject of intense work, see for example [156][157][158][159][160][161][162][163][164][165][166][167].In the context of forward dijet production, we are interested in the small-x limit of the TMD gluon distributions.This is achieved by literally taking the limit x A → 0 in the definitions.On the other hand, one can consider the leading power limit of the CGC expressions (21).This leads to the identification of the TMD gluon distributions and the leading terms in the gradient expansions of the CGC correlators [157].For example that is, the F qg is identified with the expansion of the dipole operator.For the complete list of similar relations for other TMDs see [157,168].Ultimately, in the leading power limit, the CGC formula for the qg → qg contribution can be written as where s is the hadronic center of mass energy and qg→qg are on-shell hard factors corresponding to the two independent color flows.They read where ŝ, t, û are the Mandelstam variables.
In collinear LO factorization, both TMD gluon distributions would be replaced by the collinear gluon PDFs and the above hard factors would be added.It is easy to check that the sum corresponds to the known on-shell hard factor for qg → qg process.The list of all hard factors for all subprocesses beyond the large N c limit is given in [148].
In the end, let us mention progress in establishing the CGC-leading power TMD factorization correspondence.In [155] the Authors explicitly obtained the structure of all TMD operators corresponding to five and six colored parton processes.In [140] the TMD factorization limit was studied for three jet production in CGC.Finally, in [137] the TMD factorization was studied in the context of dijet production in photoproduction at NLO.

ITMD and resummation of kinematic twists
It is important to stress, that the TMD factorization in dilute-dense collisions, despite being leading power, does take into account gluon saturation.First, notice that the leading power means here gg F (3)   gg F (4)   gg F (5)   gg F (6)   gg qg Table 1: Gauge links UC 1 and UC 2 in terms of the "staple-like" Wilson lines contributing to TMD gluon distributions that are coupled to independent color flows of gg → gg and gg → q q (upper table) and qg → qg (lower table) processes.The operators in the column marked with a star * should be traced with the gluon field strength tensor independently.
k T ≪ P T .The saturation scale is not neglected and is of the order of k T .Second, the leading twist TMD gluon distribution is considered in the strict high energy limit and they still can undergo nonlinear evolution, even at leading twist (for an interesting discussion of two types of saturation see [169]).To summarize, the dilute-dense TMD factorization is suitable when the transverse momenta of dijets are rather large and we are interested in the back-to-back dijet region.For very large transverse momenta of the jets, the highest scale is not given by the energy but the hard scale and one should switch the framework to the collinear factorization.
The main limitation of the dilute-dense TMD factorization is that it works only in the backto-back region.In dijet studies, especially at small x, dijet azimuthal correlations are the most important observables and can be measured over a wide range of the azimuthal angle.Therefore, it is essential to account not only for jet correlations but also decorrelations.A related issue is that of the transverse momentum conservation.In the TMD factorization, the transverse momentum k T of the gluon scattered (or extracted) from the target nucleus enters only TMD gluon distributions, and not the hard process.This complicates for example Monte Carlo realization of such an approach.
It turns out, that, in practical terms, it is actually quite easy to improve the dilute-dense TMD factorization.The hard factors have to take into account power corrections due to the transverse momentum of the incoming gluon.This is uniquely realized in the high energy limit by promoting the on-shell hard factors to the gauge invariant off-shell ones, following the rules described in the part devoted to HEF.This procedure was described in detail in [148], where also off-shell gauge invariant hard factors are calculated for all channels and color flows.A proof and proper interpretation within the CGC was later given in [170].First, a formal distinction between kinematic and genuine twist is made.The genuine twist counting can be simply understood as counting gluon operators in the TMD matrix element (but not the Wilson lines), for example, an operator with two strength field tensors corresponds to twist two and two-body contribution.The kinematic twists are power corrections in k T to the n-body process and come from the hard matrix element.In [170] a resummation of the kinematic twists was performed for the two-body contributions, showing that the resulting hard factors indeed correspond to the off-shell gauge invariant hard factors.
The ITMD factorization (pure gluon channel) of the p+A dijet cross section into TMD PDF (upper blob), the collinear PDF (lower blob) and the off-shell hard part (center blow).Since the exchanged momentum between the upper blob and the center is off-shell, multiple eikonal gluon exchanges are required to maintain the gauge invariance.In ITMD factorization, these gluons do not increase the genuine twist of the TMD operator.The TMD distribution is given by a certain linear combination of the operators listed in Table 1.
Before we review the dilute-dense improved TMD (ITMD) factorization framework, let us first comment, that when determining the TMD gluon distribution functions it is important to work with gauge invariant subset of diagrams and not the individual ones.The systematic way of doing that in small-x physics was suggested in [148], where a decomposition of an amplitude into color-ordered amplitudes was used.The color-ordered amplitudes [171] are gauge invariant, but contain only planar diagrams and correspond to the different ordering of the external partons.Each color-ordered amplitude comes with a color structure, for which the TMD operator can be determined.It is then usually given as some combination of the basic operators build from U [±]  Wilson lines, but corresponds to the gauge invariant hard factor.The ITMD framework given below was formulated with that feature in mind.
The ITMD factorization formula, accounting for all channels, reads  2. They are expressed in terms of the ordinary Mandelstam variables, as well as "modified" Mandelstam variables.The former read explicitly where the incoming momenta are and They sum up to ŝ + t + û = −| ⃗ k T | 2 .The "modified" Mandelstam variables take into account only the longitudinal component of the off-shell initial state k 1 and read which are related via the equation s + t + ū = 0 .
The ITMD factorization was also investigated for other processes then dijets in proton-nucleus collisions.In [172] heavy quark pair production was studied and the problem of longitudinal gluons was discussed in depth.In [173] the effect of the genuine twists vs kinematic twists was studied in detail for dijets in DIS.Finally, neglecting the longitudinal gluon contribution, the ITMD framework was formulated and applied to the trijet production [155,174].

ITMD with the Monte Carlo tool KaTie
KaTie [175] is a parton-level Monte Carlo event generator that can deal with space-like initialstate partons for arbitrary tree-level processes within the Standard Model.This means that it can be provided with k T -dependent PDFs and that it will automatically calculate the necessary matrix elements with space-like initial-state partons to create parton-level event files.These event files can be chosen to be in the LHEF-format [176].The k T -dependent PDFs can be provided via TMDlib [177], or via independent grid files, which however must be in a specific format.Furthermore, KaTie can perform calculations within ITMD factorization, and automatically calculates the necessary gauge invariant matrix elements corresponding to the color structures associated with the several gluon distributions that appear in this factorization scheme.While KaTie primarily creates event files, it also provide tools to create histograms of differential distributions.
KaTie can be downloaded from https://bitbucket.org/hameren/katie/downloads/, and a description of use is provided in the manual there.Here, we provide some background on the Monte Carlo method it employs, elucidating the procedures to be followed when using KaTie.
In hybrid factorization like ITMD factorization, the cross section can be written as a 3ndimensional integral where n is the number of final-state momenta.There are 4 initial-state variables, and 3n − 4 final-state variables.For hadro-production of dijets at tree-level the number n = 2.The function F includes the collinear PDF, the TMDS, the flux factor, and the hard matrix elements.Given a probability density G(ω) in the integration space (or phase space) that is non-zero whenever F (ω) is non-zero, we can write The Monte Carlo method is based on the Central Limit Theorem, which dictates that if {ω 1 , ω 2 , . . ., ω N } is a sequence of random phase space points, or events, independently drawn from density G(ω), then Clearly, the approximation converges to the correct result faster if the fluctuation over the terms in the sum is smaller, with the optimum when G(ω) = F (ω)/σ.Successfully applying the Monte Carlo method means that one solves the problem of bringing G(ω) to this optimum to satisfactory degree.Reaching the optimum implies having solved the integration problem and eliminates the need for the Monte Carlo procedure.Satisfaction means that the result can be obtained within an i 1 2 gg − 2F (3)   gg + F (4)  gg + F (5)  gg + N 2 c F (6)   gg acceptable time, and in practice means finding a compromise between the number of terms needed, and the complexity of the algorithm to produce the sequence of events.An advantage of the Monte Carlo method is that, given a "satisfactory" sequence generator to calculate the cross section σ, it can also be used to estimate differential distributions.Let φ(ω) be an observable, and b j (φ) = θ(φ − φ j )θ(φ j+1 − φ) be a bin for this observable between values [φ j , φ j+1 ], then The idea of an event file is to store events ω i and their weight W i = F (ω i )/G(ω i ) to produce arbitrary distributions without having to redo the generation of events.
In the context of the foregoing, KaTie performs roughly speaking two tasks: it calculates the hard matrix elements, automatically and efficiently, as part of the evaluation of F (ω), and, provided with PDFs and TMDs, it creates reasonably-sized event files.In order to perform the first task, Ka Tie employs the Dyson-Schwinger approach to calculate tree-level helicity amplitudes numercially, as first proposed in [178] and also utilized in other tree-level programs, notably Alpgen [179].The feature of space-like initial-state momenta is dealt with following the auxiliary-parton method outlined in [145].The second task is achieved through the application of many well-known optimization techniques, a few of which are worth to mention in order to understand the procedures in KaTie's operation.
One is adaptive importance sampling, which optimizes adaptable probability densities iteratively, using the information of events generated so far in each iterative step.While being very effective, it can cause bias in those events, and they must not be used for the actual event file.Consequently, there must be a separate optimization stage before the generation of the event file starts.This happens for each partonic subprocess separately.Once this has been performed, the information of the optimized densities is stored, and is from then on used to generate an arbitrary number of event files.
Secondly, KaTie employs rejection.Phase space cuts are required to mimic detector acceptance and to avoid singularities in tree-level matrix elements.These typically are formulated in terms of variables that are non-trivial functions of the variables in ω that are actually generated, and cannot be implemented as exact integration bounds.Instead, the integrand F (ω) is imagined to vanish outside those phase space cuts, and a bigger enveloping space is generated.This leads to many events with W i = 0, which are however not stored but are included in the eventual normalization of the weights.In order to achieve the unbiased normalization of the event weights in case several event files are generated, so called raw files are stored instead of actual event files, containing non-normalized weights and more statistical information.Event files in the LHF format can be extracted from these.
The third method that needs to be mentioned is unweighting.This is a statistical procedure to reduce the fluctuations of the weights W i while keeping the event file sound.In practice this procedure "weeds out" low-weight events and reweighs remaining ones.The cross section estimate is not affected (is also not getting better) but it allows to reduce the required size of the event file.While it causes many events to get the same constant weight, there may occur events with higher weights.
The phase space is only 3n-dimensional, but it is much more convenient to write the events in terms of the 4n variables being the parton-level initial-state and final-state momenta.We will still denote events with the symbol ω.While all components of F (ω i ), like the matrix element, the PDFs, the strong coupling etc., can be re-evaluated using these, it is more convenient to also store their values besides W i .This is useful if one wants to employ the reweighting procedure.For example, let an event file be created within hybrid k T -factorization with a TMD F(ω) (we simply imply that this function selects the appropriate variables from ω. Realize that this in practice includes the final-state ones due to the factorization scale dependence).Suppose a user has their own TMD F ′ but only in a numerical form that cannot (yet) be married with KaTie or TMDlib.The value of F i = F(ω i ) is stored for each event, and the event file can be transformed to be valid for F ′ by multiplying for each event.It will increase the fluctuation of the weighs and reduce the quality of the event file, but most likely not in a drastic way.As mentioned earlier, the use of KaTie is described in the manual, but below we present the complete single input card that allows to create event files for dijets within ITMD.

TMD gluon distributions
From the BK equation discussed in the earlier sections one can obtain TMD dipole gluon density.In this section we will use k 2 = | ⃗ k T | 2 as argument of gluon density.This is the standard notation used in discussion of angular averaged distributions.Let's define the following Fourier transforms [180,181] ) as well as explicit relation between the two functions F(x, k 2 ) and Φ(x, k 2 ).
Since we are working in the large target approximation as well as in the forward limit i.e. we neglect any momentum transfer therefore we introduce the following normalization conditions for integration over the impact parameter where S(b) is the profile function The momentum space equation for the Φ(x, k 2 ) assumes the form: where α = N c α s /π.The the equation for F(x, k 2 ) is obtained from eq. ( 56) by inserting in the nonlinear part of it relation expressing Φ(x, k 2 ) in terms of F(x, k 2 ) and acting on the whole equation with the operator that transforms Φ(x, k 2 ) to F(x, k 2 ) see [181] for the details of this transformations.In the end one obtains [53,181,182] The linear part of the equation above is the well known BFKL kernel while the nonlinear part is the triple pomeron vertex.The triple pomeron vertex has such property that it is dominated by the anticollinear pole.We see that as one evaluates the gluon at lower and lower values of k 2 the integration in the nonlinear term is over larger domain giving larger contribution and suppressing the gluon density.One can also see that in the collinear limit the nonlinear term completely drops and one is left with linear equation.

BK equation within Kwiecinski Martin Staśto model
The ( 57) is a leading order BK equation for dipole gluon density.In order to make it applicable to phenomenology it has been extended to account for higher order corrections following Kwiecinski-Martin-Stasto (KMS) prescription [183] that was originally applied to the BFKL equation.Those include • kinematical constraints which enforce that the virtuality of exchanged gluon is dominated by its transverse component This constrain suppresses anticollinear pole and therefore suppresses the diffusion into infrared • nonsingular at small z pieces of the splitting function • running coupling constant • sea quark singlet contribution to match the DGLAP limit at large z With these corrections the equation assumes the form [181] where Σ(x, k 2 ) is a sea quark distribution obeying essentially DGLAP equation in unintegrated form (further details can be found in [181]).Please note that in the equation above lower cut in k 2 = k 2 0 was introduced.The origin of the cut is related to the method of solving of the equation and in principle can be set to arbitrarily small value.
In Ref. [12], the following initial condition for Eq. ( 59) was fitted to the F 2 proton structure data from HERA [184] with For k 2 ≤ 1 GeV 2 , the gluon distribution was taken as F(x, k 2 ) = k 2 F(x, 1), which is motivated by the shape obtained from the solution of the LO BK equation in the saturation regime [185].
The fitting procedure gave the following numerical values for the parameters: N = 0.94, β = 18.6, D = −82.1 and R = 2.40 GeV −1 .The overall quality of the fit was good, with χ 2 /ndof = 1.73.We shall refer to this gluon distribution as the Kutak-Sapeta or KS gluon.
For completeness, the fit of linearized version of Eq. ( 59), i.e. with the last term dropped, was performed as well, and the following parameters were obtained: N = 0.004, β = 26.7,D = −51102 and χ 2 /ndof = 3.86.The presence of the parameter R, characterizing the target, allows one to obtain dipole gluon distribution of nuclei.In order to do that, one uses relation R A = d A 1/3 R, which in the end gives the enhancement by A 1/3 of the nonlinear term for gluon density normalized to the number of nucleons [12].The parameter d is a phenomenological factor that was was varied between d = 0.5 and d = 1.0.In the following computations we used the "least saturation" scenario with d = 0.5.

The Sudakov Resummation
An important class of perturbative corrections to scattering process, appearing when the emitted partons are both collinear and soft, is resummed in terms of the Sudakov form factor [186][187][188].It appears due to not exact cancellation of virtual corrections and real corrections, as a consequence of certain exclusivity of the final state.In the processes considered here the largest effect of the Sudakov is expected to affect the cross section when the jets are in nearly back-to-back configuration.The large logarithm comes since the transverse momenta of jets can be sizable while the imbalance k T of incoming space-like parton is small.The effect of Sudakov form factor is to suppress the back-to-back configuration and enhance the moderate angle part of the distribution, leading to so-called broadening.Within the small-x formalism the Sudakov form factor can be factorized, in the coordinate space, from the hard process [186][187][188].The complete cross-section accounting for the Sudakov form factor reads: where ag→cd is the Fourier transform of the TMD gluon distributions and S ag→cd are the Sudakov factors written for each each channel where S i p and S i np are the perturbative and non-perturbative contributions.The perturbative Sudakov factors, including double and single logarithms, are given by [186,189] where β 0 = (11 − 2n f /N c )/12.The gg → q q channel is negligible in kinematic domain considered here.In the above the scale µ b is the inverse of the impact parameter: Given this selection, the scale µ b becomes constant at the point of high b T , where its value is determined to be 2e −γ E /b max , which is considerably greater than Λ QCD .Following Ref. [190], we will adopt b max = 0.5, GeV −1 .
For completeness we comment briefly on DGLAP based prescriptions for incorporating the Sudakov form factor.The method relies on constructing the Sudakov form factor from DGLAP splitting function and using it to reshuffle events according to relation between hard scale and transverse momentum of the gluon.In such constructions one chooses some inclusive quantity to be unmodified while allowing for modification of unintegrated quantity.Two such methods were presented in [24,191].In the former the total cross section was preserved while in the later the integrated gluon density was unmodified.

Kutak-Sapeta (KS) gluon distribution
We shall now discuss the KS gluon, introduced above, in variants with and without the Sudakov resummation.The parameters were fixed by the original fit [12] with no Sudakov factors and the gluon was later used without any modifications.Hence, combining it with the Sudakov does not introduce new parameters.This is true because the perturbative part ( 64) and ( 65) is parameterfree while the non-perturbative terms are universal in the kinematic domain of our study [192].We introduce the Sudakov effects into the KS gluon distribution following the formalism described above.In addition, for reference, we use two methods employed in our earlier studies [24,191].Those calculations used the Sudakov form factor, understood as the DGLAP evolution kernel, that has been applied on the top of the gluon TMD, together with constrains such as unitarity.Those methods should therefore be considered as models, in contrast to the proper resummation of Sudakov logarithms described in the preceding section.Nevertheless, the approaches used in Refs.[24,191] were phenomenologically successful (see also [14]), and it is therefore useful and interesting to compare the predictions of those simplistic models with the proper way of including the Sudakov effects into the small-x gluon.
The reference models are: • Model 1: The survival probability model [191], where the Sudakov factor of the form [193] T is imposed at the level of the cross section.This procedure corresponds to performing a DGLAP-type evolution from the scale µ 0 ∼ k T to µ, decoupled from the small-x evolution.
• Model 2: The model with a hard scale introduced in Ref. [24].The Sudakov form factor of the same form as in Eq. ( 68) is imposed on top of the KS gluon distribution in such a way that, after integration of the resulting hard scale dependent gluon TMD, one obtains the same result as by integrating the KS gluon distribution.
In Fig. 7, we show the KS gluon distributions, with and without Sudakov form factors, as functions of the transverse momentum k T and the hard scale µ.Three columns correspond to three different x values.The first row shows the original KS gluon distribution, which, as expected, does not depend on the value of µ.In the second row, we show the KS hardscale gluon distribution of Ref. [24] (the other model [191] does not allow one to plot gluon distribution, as it applies Sudakov effects at the cross section level via a reweighting procedure).Here, the dependence on µ is nontrivial and we see that the gluon develops a maximum in that variable.As shown in the figure, this maximum is rather broad.In the third and the fourth row of Fig. 7, we present the KS gluon distribution with the Sudakov form factor from Eqs. ( 63)- (65).As explained earlier, this gluon exists in two versions, one for the qg and the other for the gg channel.The dependence on k T and µ is qualitatively similar between the these gluons and the naive KS hardscale gluon distribution.In the former case, however, the peak is significantly narrower in µ as compared to the naive model of Ref. [24].It is interesting to note that the qg gluon is broader than the gg gluon.This can be understood by comparing the colour factors in the Sudakov functions (64) and (65).Since the colour factor for the gg channel is larger than for the qg, the Sudakov suppression is stronger along the µ direction in the former case.
We have as well computed linear versions of the KS gluon distributions with the Sudakov, using the KS linear gluon distribution of Ref. [12].The gluons discussed in this section are available publicly as part of the KS package and can be downloaded from http://nz42.ifj.edu.pl/~sapeta/ KSgluon-2.0.tar.gz.

ITMD distributions from KS gluon
All the gluon TMDs ( 69)-( 73) can be calculated from a single xG (2) x, | ⃗ k T | distrubution in the above Gaussian approximation.However, because the KS gluon provides directly an impactparameter-integrated distribution, it is not straightforward to identify S ⊥ and obtain F x, | ⃗ k T | from Eq. (74).To address this issue, we aplied the following procedure [194].We first computed the dipole cross section σ dipole (x, r = |r|) = 2 d 2 b N F (x, r) from xG (2) x, | ⃗ k T | by inverse Fourier transformation of Eq. ( 75), and then defined S ⊥ as its value at large r, i.e. when it saturates (since in that limit N F → 1): We can now obtain F x, | ⃗ k T | and calculate all the needed gluon TMDs.Their behavior as a function of k t = | ⃗ k T | is plotted in Fig. 8, both for the proton and the lead nucleus.The small mismatch between their high-k t behavior, expected due to the initial condition for the x evolution, can be observed.Similarly, we computed [195] the ITMD distributions from the KS gluon with the proper QCD Sudakov, shown in Fig. 7.
Within the Gaussian approximation, one can also derive the following formula for the WW gluon density [194] where F qg is the dipole gluon density and S ⊥ is the target's transverse area.Using the procedure described above, we computed also the WW gluon, which we show in Fig. 9 in proton (left) and lead (right), with and without Sudakov form factors, as functions of the transverse momentum | ⃗ k T | and the hard scale µ, for one particular x = 10 −3 .(The gluon density is available from the TMDlib [177].) First of all, let us notice that the WW gluon distribution has no maximum, contrary to the dipole gluon [12,194].Secondly, we see that the Sudakov factor suppresses the gluon distribution at low | ⃗ k T | and enhances it at higher | ⃗ k T |.Because the Sudakov form factor is derived in the regime µ ∝ |⃗ p T | ≫ | ⃗ k T |, we apply it only to that part of the gluon density where µ > | ⃗ k T |.In the remaining domain, we use the gluon without Sudakov, given in Eq. ( 80).This is visible in Fig. 9 as a kink of the curve corresponding to µ = 17 GeV.(A similar kink exits also for the µ = 67 GeV curve but it is located at larger values of log All the ITMD gluons discussed in this above are available publicly and can be downloaded from http://nz42.ifj.edu.pl/~sapeta/itmd-KS.tar.gz.

Forward dijets at the LHC
In the following section we review some existing predictions for forward dijet production in protonproton and proton-lead collisions at LHC energies, obtained within the ITMD framework.We also include some new computations, not published before.
Before we proceed, it is important to mention, that the first predictions for forward dijets in saturation formalism using the ITMD, that -as discussed in Section 2.3 -approximates the CGC for sufficiently large transverse momenta, appeared in [194].This computation significantly improved the predictions of [23] obtained with k T -factorization, which is not entirely accurate at small dijet imbalance.Later, in [196] a full computation within CGC was compared to ITMD, confirming, that at larger transverse momenta ITMD is adequate.None of the above computations included the Sudakov resummation, however.In [14] a calculation has been performed that included both the full ITMD and the Sudakov resummation "model" [191], based on the reweighting the events with the DGLAP Sudakov form factor.The Authors compared the shape of the obtained azimuthal dijet correlations for proton-proton and proton-lead with those obtained by ATLAS collaboration [25] (no dijet cross section was actually measured, only the conditional yields).Interestingly, when the shapes of the azimuthal correlations for p-p and p-Pb are stacked together so that they match in the first bin, one can clearly see broadening of the p-Pb cross section.It turns out, that similar broadening is obtained within the ITMD, if both the saturation as well as the Sudakov resummation are present.We show this result in Fig. 10.
In the following, we shall report on further advances in dijet computations within ITMD, with the full Sudakov resummation [197].In particular, we will present the numerical results for the differential cross sections in terms of the azimuthal angle ∆ϕ between the leading and the subleading jets for the proton-proton and the proton-lead collisions at LHC energies.We shall also discuss the nuclear modification ratio R p−Pb where O is a differential related to an observable.Finally, we shall investigate differential cross sections in the rapidity of the leading and sub-leading jet.The cross sections were computed using the KaTie Monte Carlo program [175] within the ITMD factorization, both described in the preceding Sections.We considered the proton-proton and proton-lead collisions at √ s = 5.02 TeV, 8.16 TeV and 8.8 TeV per nucleon.For protonproton collisions, we also computed the proton-proton cross section for √ s = 14 TeV.In order to define the leading and the sub-leading jets, mentioned above, we used the anti-k T jet clustering algorithm [198] with a radius of R = 0.4.Since our computation is leading order, the jet algorithm is actually equivalent to a simple cut in rapidity-azimuthal plane.Motivated by the current and planned LHC experiments, we applied the following cuts to the transverse momentum of these jets: Specifically, we used the first three cuts i) − iii) for the transverse momentum of the jets in the rapidity range 2.7 < y ⋆ 1 , y ⋆ 2 < 4.0, both for proton-proton and proton-lead collisions at √ s = 8.16 TeV.These cuts correspond to the FCal calorimeter of the ATLAS detector and are motivated by the measurement [25].The last two cuts iv)−v) were applied in the rapidity range 3.8 < y ⋆ 1 , y ⋆ 2 < 5.1, both for proton-proton and proton-lead collisions at √ s = 5.02 TeV, and 8.8 TeV energies per nucleon.These correspond to the planned FoCal extension of the ALICE detector [27].For the same kinematic domain (rapidity and transverse momentum cuts for the jets), we also considered protons collisions at √ s = 14 TeV (proton-lead collisions are not feasible at this energy).The result for 5.02 TeV and 14 TeV were not published in [197] and are thus new.
The factorization and renormalization scales were set using the transverse momentum of the leading and the sub-leading jets µ = (p T 1 + p T 2 )/2.For the TMD gluon distributions, we used the ones calculated in [194] based on the Kutak-Sapeta (KS) fit of the dipole gluon density [12], as described in the preceding subsections.For the collinear PDFs we used the CTEQ10NLO PDF set [199] from LHAPDF6 [200].For the computation of the cross sections, we included the channels qg * → qg, with five quark flavours, and gg * → gg.The channel gg * → qq was neglected because its contribution, for the considered kinematic domain, is quite small [12,194].
Let us now discuss the results.In Fig. 11 we show the results for the azimuthal correlations for p-p and p-Pb collisions in the ATLAS kinematics at √ s = 8.16 TeV.We compare the Sudakov resummation obtained via two methods, the full-b space resummation (dotted) and the simplified approach, where the collinear PDF is not affected by the Sudakov form factor.As can be seen, they differ a bit, but as we shall see below, the difference cancels in the nuclear modification ratio, thus does not affect the saturation signal.Moreover, as we also show below, there are very large scale uncertainties, within which both methods are qualitatively similar.In our computations we used the ITMD factorization alone, without parton shower or hadronisation corrections.In order to asses this effect, we computed cross section with PYTHIA [201,202] with all correction turned on, and then just with initial-state parton shower.The latter roughly corresponds to the TMD framework (as discussed eg. in [203]), therefore both calculations allow for extraction of a "correction factor".We repeated the same procedure using the nucler PDFs in PYTHIA.In Fig. 12 we applied that correction to KaTie results and compared with PYTHIA computations.Although the correction factor is rather large, as we shall see below it does not affect the saturation signal in the nuclear modification ratio.Next, in Fig. 13 we show the error bands due to both scale variation.Finally, in Fig. 14 we compute the nuclear modification ratio.For all three p T cuts, we show the ITMD result with simiplified Sudakov resummation, with the full-b space resummation, and the result where the non-perturbative correction factor from PYTHIA was applied.We also compute the error band due to scale variation and due to the correction factor.We see, that the saturation effect due to the nonlinear evolution is visible, but not very significant for the considered kinematics.For the lower p T cut we observe suppression of about 15%, but the uncertainty due to the correction factor is not much less.Interestingly (but understandably) this error decreases for larger p T bins.
Next, we move on to the more forward region with FoCal kinematics.Here we focus first on the p T > 10 GeV cut.The results for azimuthal correlations are shown in Fig. 15.In one figure we show ITMD result with Sudakov obtained via both models, as well as PYTHIA corrections and errorbands, as discussed before.In Fig. 16 we compute the nucler modification ration for that case.As can be seen from the plot, the suppression due to saturation is very large, about 25-30%, and is not destroyed by the errors, both due to the scale variationn and due to the non-perturbative effects.In Fig. 17 we show latest results, where we compare p-p and p-Pb cross sections at different energies, computed also with larger p T > 20 GeV cut.As expected, the difference between p-p and p-Pb is smaller in that case.Next, in Fig. 18-19 we present latest results for rapidity distributions for the two p T cuts, for the leading jet, and the sub-leading jet.As the shape of the rapidity distribution is correlated with the x dependence, and thus with the evolution in the energy, measurement of the rapidity distribution may provide a valuable discriminatory tool for the evolution equations.Finally, in Fig. 20 we show the nuclear modification ratio for the two energies and both p T cuts.As can be seen, with forward FoCal kinematic, even with the p T cut of 20 GeV the suppression is about 20%.42) (solid lines), the full b-space resummation Eq. ( 62) (dotted lines).The plots were taken from [197].TeV in the FCal ATLAS kinematics.The solid lines represent the results from KaTie computed within the ITMD approach, and the dotted lines represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The bands represent the error due to the variation of the factorization/renormalization scales from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The plots were taken from [197].62) (blue line).The band represents the error due to the variation of the factorization/renormalization scales from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The points ▽ represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The error bars associated with these points account for both errors due to variation of scale in KaTie and the statistical uncertainties in the correction factor from Pythia.The plots were taken from [197].62) (dotted lines).In the second plot, the solid lines represent the cross sections from KaTie computed within the ITMD approach, the points represent the results from Pythia with different components, and the dotted lines represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.In the third plot, the solid lines represent the cross sections from KaTie computed within the ITMD approach, and the dotted lines represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The bands represent the error due to the variation of the factorization/renormalization scales in KaTie from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The plots were taken from [197].62) (blue line).The band represents the error due to the variation of the factorization/renormalization scales from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The points ▽ represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The error bars associated with these points account for both errors due to variation of scale in KaTie and the statistical uncertainties in the correction factor from Pythia.The plots were taken from [197].

Forward dijets at Electron Ion Collider
While the focus of this review is on forward jets at LHC, the new accelerator Electron Ion Collider (EIC) is just around the corner [29].A vital part of its research program is to look for gluon saturation.However, since the involved hard scales are lower than at the LHC, in fact, hadrons constitute better probes of the saturation physics than jets.Nevertheless, there has been a large interest in jets at EIC, see for example [173,204,205].Moreover, at the partonic level, both hadron and jet production need classes of the same diagrams.Therefore, the upcoming EIC experiment triggered tremendous progress in computing NLO corrections for the DIS process within the CGC formalism, see for example [94,99,132,137,206,207].
In the context of the ITMD factorization, dijet production at EIC is actually a very interesting and important process.As discussed eg. in [157] the only TMD gluon distribution that is relevant for dijets in DIS, in the back-to-back limit, is the Weizsacker-Williams (WW) gluon distribution F (3) gg (called also xG (1) )2 .For the ITMD factorization, this is only true for sufficiently small photon virtuality Q 2 ≪ p T , otherwise longitudinaly polarized gluons contribute, and they are not formally accounted for in the ITMD approach.Due to this fact, for processes where such contributions formally enter, the ITMD factorization is called ITMD * (see for example [174] for the trijet production case).
The dijet production in the DIS process gives actually more interesting observables than dijets in hadron-hadron scattering.This is because there is also an identified electron in the final state.Thus, in addition to dijet azimuthal correlations, one can also study correlations with the outgoing electron.
The ITMD * factorization formula for the process Figure 17: The differential cross sections in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 5.02 TeV (blue) and 8.8 TeV (red) in the ALICE FoCal kinematics.For the same kinematics, each plot also represents the differential cross section for the proton-proton collision at √ s = 5.02 TeV (black).Since proton-lead collisions are not feasible at this energy, we didn't compute those.All of these cross sections were computed using KaTie within the ITMD factorization formula with the simplified Sudakov resummation Eq. ( 42) The solid lines represent the results for proton-lead collisions and the dotted lines represent the results for the protonproton collisions.reads As mentioned above, there is only WW gluon distribution in the above formula.
In [209] we computed several observables using the above framework, not only to quantify possible saturation effects in dijet observable at EIC but also to demonstrate the relevance of interplay of the saturation, the Sudakov effect, and the exact kinematics.
First, we found that in the context of saturation physics, it is best to consider forward dijets.This provides a very good focusing of the x A distribution around the sufficiently small values so that the use of the small-x formalism is well justified.For the WW gluon distribution, we used the Kutak-Sapeta fit of the dipole TMD to HERA data, and then the WW TMD distribution was calculated in the mean field approximation (see Section 4.5).The calculation was supplemented with the Sudakov form factor, see Fig. 9 and the corresponding discussion in the text.In Fig. 21 we show prediction for the cross section dependence on azimuthal angle between the dijet system and the scattered electron for CM energy per nucleon √ s = 90 GeV.The calculations are done both in the lab frame and the Breit frame.One clearly sees effects coming from the Sudakov form factor while the saturation is rather mild, despite the cut on the jet transverse momentum (in the Breit frame) is very low.Primarily, it is a consequence of the saturation pattern visible in the WW gluon distribution.Unlike the dipole gluon distribution, it does not have a peak as a function of the transverse momentum, c.f. Fig. 9 and 8.In our computation, we put quite a low cut on the photon virtuality Q 2 > 1 GeV 2 in order to suppress the longitudinal gluons.For a complete discussion see [209] and [210].

Summary and outlook
The text discusses the challenges and selected ongoing research related to searches of gluon saturation in QCD using dijet observables.Although there are indications of saturation in the experimental data, achieving a complete consensus on its existence is still a challenge due to demanding kinematics and the complex collision environment.While there are simpler observables than dijets, the latter have the advantage of providing the azimuthal correlations, which in the back-to-back kinematics are sensitive to saturation even for relatively hard jets.
The study of forward dijet production, using the Improved TMD Factorization (ITMD) framework, is a relatively novel proposal in the search for gluon saturation.The framework has been already used in addressing many dijet related observables that were reviewed in this work.The text highlights the need for more measurements and the construction of the FoCal detector by the ALICE collaboration to shed more light on the saturation phenomenon.Furthermore, it discusses complementary dijet measurements in photon nucleus scattering that can be done at the Electron Ion Collider.
While in the review we focused our attention on dijets, there are other forward physics final states that can be used to search for saturation effects [211].Those include di-hadrons [13,212,213], J/ψ production [214], and trijets [140,174].The processes are complementary to dijets, however, in some aspects more challenging.The di-hadrons require knowledge of fragmentation function and are sensitive to longitudinal polarization effects while trijets are just more difficult to measure and interpret.
The future research will be focused on increasing the precision of the calculations.The ITMD framework is well suited to extrapolate methods known in collinear factorization to compute higher order corrections to off-shell gauge invariant matrix elements, and in principle to automatize NLO computations.First steps toward that goal have already been undertaken [215][216][217][218]. Also, it is important to increase the accuracy of the evolution equations and there is also a substantial progress in this direction, see for example [105,167,[219][220][221][222][223][224][225].Last but not least, the ITMD framework should account for the longitudinally polarized gluons, which contribute to some processes (see eg. [226]).The corresponding amplitudes should be obtainable automatically, similarly to to off-shell gauge invariant amplitudes coupled eikonally that are used currently in the ITMD.Research in that direction is also ongoing.

Figure 1 :
Figure 1: The momentum assignment in the inclusive dijet production in p+A collision.The multiple double lines represent the nuclear target A. J are the jets with momenta p1 and p2.X denotes arbitrary final states.

Figure 3 :
Figure3: Scattering of a color dipole q → qg off a nuclear target in CGC.The shaded ellipses symbolize the color field of a nucleus.This particular diagram corresponds to a contribution, where in the amplitude the incoming quark scatters off the color field of the nucleus before the splitting, while in the amplitude conjugate both gluon and a quark scatter after splitting.Particles scattering off the color field of the nucleus gain infinite straight Wilson line along the light cone, in their respective representations.

Figure 4 :
Figure 4: Off-shell gauge invariant amplitude (upper blob -the zigzag lines represent off-shell gauge invariant gluons) can be constructed by promoting a single off-shell gluon coupled eikonally to a straight infinite Wilson line along corresponding hadron momenta PA or PB.Here we show diagrams for the production of two gluons at tree level.The double line represents the momentum space Wilson line along PA (for the top line) and PB (for the bottom line).Gluons couple to these Wilson lines via igt a p µ A(B) and the double-line propagators have form i/(k • p A(B) + iϵ).Only planar diagrams are shown.The blob represents all possible connections of gluons via standard vertices.

Figure 5 :
Figure 5:The shape of the "staple-like" gauge links U[+] and U[−] .The horizontal axis represents the light-cone minus direction, the vertical axis symbolizes the transverse displacement.The transverse pieces are placed at +∞ and −∞, respectively.
(i) ag * →cd ⃗ k T are off-shell gauge invariant hard factors and Φ (i) ag→cd x A , | ⃗ k T | are the corresponding TMD gluon distributions.We collect them in Table

Figure 7 :
Figure 7: KS gluon distribution, without and with the Sudakov form factors.The second row corresponds to the simple model-Sudakov given in Eq. (68), while the third and the fourth rows show results obtained with the Sudakov factors derived from QCD and given in Eqs.(63)-(65).

Figure 8 :
Figure 8: The KS gluon TMDs as a function of log(| ⃗ k T | 2 /GeV 2 ) at x = 1.6 10 −4 for the proton (left) and the lead nucleus (right).Since F

Figure 9 :
Figure 9: The WW gluon density in the proton (left) and lead (right), with and without Sudakov resummation, as a function of the transverse momentum of the gluon, for various values of the hard scale.

Figure 10 :
Figure 10: Broadening of azimuthal correlations in p-Pb collisions vs p-p collisions for different sets of cuts imposed on the jets' transverse momenta.The blue and red bands show the normalized differential cross sections in azimuthal angle ∆ϕ, respectively for for p-Pb and p-p, shifted so that they match in the first bin.The points show the experimental data [25] for p-p and p-Pb, where the p-Pb data were shifted by a pedestal, so that the values in the bin ∆ϕ ∼ π are the same.Theoretical calculations are represented by the histograms with uncertainty bands coming from varying the scale by factors 1/2 and 2.

Figure 11 :
Figure 11: The differential cross sections in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 8.16 TeV in the FCal ATLAS kinematics.These were computed using KaTie within the ITMD factorization formula with: the simplified Sudakov resummation Eq. (42) (solid lines), the full b-space resummation Eq. (62) (dotted lines).The plots were taken from[197].

Figure 12 :
Figure 12:The differential cross sections in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 8.16 TeV in the FCal ATLAS kinematics.The solid lines represent the results from KaTie computed within the ITMD approach, the points represent the results from Pythia with different components, and the dotted lines represent the Ka Tie results corrected with the non-perturbative correction factor extracted from Pythia.The plots were taken from[197].

Figure 13 :
Figure 13: The differential cross sections in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 8.16TeV in the FCal ATLAS kinematics.The solid lines represent the results from KaTie computed within the ITMD approach, and the dotted lines represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The bands represent the error due to the variation of the factorization/renormalization scales from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The plots were taken from[197].

Figure 14 :
Figure 14: Nuclear modification ration R p−Pb in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 8.16 TeV in the FCal ATLAS kinematics.The solid lines represent the results from KaTie computed within the ITMD approach with: the simplified Sudakov resummation Eq. (42) (red line), the full b-space resummation Eq. (62) (blue line).The band represents the error due to the variation of the factorization/renormalization scales from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The points ▽ represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The error bars associated with these points account for both errors due to variation of scale in KaTie and the statistical uncertainties in the correction factor from Pythia.The plots were taken from[197].

Figure 15 :
Figure 15: The differential cross sections in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 8.16 TeV in the ALICE FoCal kinematics.The first plot represents the corss sections computed using KaTie within the ITMD factorization formula with: the simplified Sudakov resummation Eq. (42) (solid lines), the full b-space resummation Eq. (62) (dotted lines).In the second plot, the solid lines represent the cross sections from KaTie computed within the ITMD approach, the points represent the results from Pythia with different components, and the dotted lines represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.In the third plot, the solid lines represent the cross sections from KaTie computed within the ITMD approach, and the dotted lines represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The bands represent the error due to the variation of the factorization/renormalization scales in KaTie from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The plots were taken from[197].

Figure 16 :
Figure 16: Nuclear modification ration R p−Pb in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 8.16 TeV in the ALICE FoCal kinematics.The solid lines represent the results from KaTie computed within the ITMD approach with: the simplified Sudakov resummation Eq. (42) (red line), the full b-space resummation Eq. (62) (blue line).The band represents the error due to the variation of the factorization/renormalization scales from a value of (pT 1 + pT 2)/2 by a factor of 1/2 and 2. The points ▽ represent the KaTie results corrected with the non-perturbative correction factor extracted from Pythia.The error bars associated with these points account for both errors due to variation of scale in KaTie and the statistical uncertainties in the correction factor from Pythia.The plots were taken from[197].

Figure 18 :
Figure 18: The differential cross sections in terms of the rapidity of the leading jet y ⋆ 1 for the proton-proton and the proton-lead collisions at √ s = 5.02 TeV (blue) and 8.8 TeV (red) in the ALICE FoCal kinematics.For the same kinematics, each plot also represents the differential cross section for the proton-proton collision at √ s = 5.02 TeV (black).Since proton-lead collisions are not feasible at this energy, we didn't compute those.All of these cross sections were computed using KaTie within the ITMD factorization formula with the simplified Sudakov resummation Eq. (42) The solid lines represent the results for protonlead collisions and the dotted lines represent the results for the proton-proton collisions.

Figure 19 :
Figure 19: The differential cross sections in terms of the rapidity of the sub-leading jet y ⋆ 2 for the protonproton and the proton-lead collisions at √ s = 5.02 TeV (blue) and 8.8 TeV (red) in the ALICE FoCal kinematics.For the same kinematics, each plot also represents the differential cross section for the protonproton collision at √ s = 5.02 TeV (black).Since proton-lead collisions are not feasible at this energy, we didn't compute those.All of these cross sections were computed using KaTie within the ITMD factorization formula with the simplified Sudakov resummation Eq. (42) The solid lines represent the results for proton-lead collisions and the dotted lines represent the results for the proton-proton collisions.

Figure 20 :
Figure 20: Nuclear modification ration R p−Pb in terms of the azimuthal angle ∆ϕ between the leading and the sub-leading jets for the proton-proton and the proton-lead collisions at √ s = 5.02 TeV (blue) and 8.8 TeV (red) in the ALICE FoCal kinematics.The solid lines represent the results from KaTie computed within the ITMD approach with the simplified Sudakov resummation Eq. (42) for pT 1, pT 2 > 10 GeV.And the dotted lines represent the results for pT 1, pT 2 > 20 GeV.

Figure 21 :
Figure 21:  Azimuthal correlations between the total transverse momentum of the dijet system and the transverse momentum of the scattered electron at EIC in two frames: the LAB frame (left), the Breit frame (right).The calculation has been done within the ITMD* framework with the Weizsäcker-Williams gluon distribution obtained from the Kutak-Sapeta fit to HERA data.

Table 2 :
The top table lists the combinations of the TMD gluon distributions listed in Table1that correspond to gauge invariant off-shell hard factors of the ITMD factorization formula.The bottom table gives explicit formulae for the LO off-shell gauge invariant hard factors of the ITMD factorization.The ŝ, t and û are ordinary Mandelstam variables, whereas s, t, ū are invariants where instead of the incoming off-shell gluon momentum its longitudinal component is used.