Forward trijet production in p-p and p-Pb collisions at LHC

We calculate various azimuthal angle distributions for three jets produced in the forward rapidity region with transverse momenta $p_T>20\,\mathrm{GeV}$ in proton-proton (p-p) and proton-lead (p-Pb) collisions at center of mass energy $13\,\,\mathrm{TeV}$. We use the multi-parton extension of the so-called small-$x$ Improved Transverse Momentum Dependent factorization (ITMD). We study effects related to change from the standard $k_T$-factorization to ITMD factorization as well as changes as one goes from p-p collision to p-Pb. We observe rather large differences in the distribution when we change the factorization approach, which allows to both improve the small-$x$ TMD gluon distributions as well as validate and improve the factorization approach. We also see significant depletion of the nuclear modification ratio, indicating a possibility of searches for saturation effects using trijet final states in a more exclusive way than for dijets.


Introduction
So far, the complete trijet production in high energy hadron-hadron collisons with initial off-shell partons was discussed only within the k T -factorization approach [1,2]. Here we wish to discuss for the first time the production of three jets within the so-called small-x Improved Transverse Momentum Dependent (ITMD) approach [3] that accounts for gluon saturation effects and involves several transverse momentum dependent (TMD) gluon distributions.
The trijet production process at LHC kinematics is of great interest, since, as follows from a recent study by three of us [2], it has great potential to uncover details of dynamics related to transverse momentum dependence of proton constituents and to test properties of parton showers. Furthermore it allows to study the features of ITMD factorization, and, last but not least, it constitutes the real emission contribution to dijet production at NLO.
Before we continue with three jets, let us first summarize the theoretical formalisms used to calculate forward jet production in proton-proton and proton-nucleus collisions. This will also allow to set up the terminology to avoid possible confusion.
There are two quite distinct factorization approaches that make use of the parton distribution functions depending on the transverse momentum. First are the leading power factorization theorems of QCD, usually called the TMD factorization theorems [4]. Because they hold to leading power in k T /µ at fixed energy (k T being the transverse momentum of incoming partons and µ the hard scale), the partonic processes entering the factorization formulae are calculated fully on-mass-shell, i.e. the transverse momentum k T of incoming partons is neglected in the hard part 1 . It is not neglected, however, in the hadronic soft part. Thus the only dependence on incoming parton momenta is in the PDFs. The evolution equations for the TMD PDFs are the Collins-Soper-Sterman equations and they resum the large logarithms of k T /µ [7]. On the other hand the k T -factorization also called the high energy factorization (HEF) focuses on the small Bjorken x limit, not neglecting the powers k T /µ [8][9][10]. Thus, the incoming partons (usually gluons, as they dominate at high energy) carry the transverse momentum and the partonic processes must be calculated off-shell. In this context, the TMD PDFs are also called unintegrated PDFs. The corresponding evolution equations typically resum the logarithms of energy (or equivalently 1/x) by means of the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [11,12] or its extensions, like for instance the Catani-Ciafaloni-Fiorani-Marchesini (CCFM) [13][14][15] or approaches/models combining BFKL and DGLAP like Kwieciński-Martin-Staśto (KMS) equation [16] or Kimber-Martin-Ryskin-Watt [17] (KMRW).
The BFKL-type approach does not account for very large gluon densities inside hadrons -the evolution equations are linear and the gluon densities can grow power-like with energy, eventually violating the unitarity bound. The QCD theory predicts however a nonlinear generalization of the BFKL equationthe Balitsky-Kovchegov (BK) equation [18,19], which exhibits gluon saturation, i.e. a state where almost all gluons have momenta k T ∼ Q s , where Q s is the saturation scale. The BK equation is the mean field approximation to the more general system of equations known as the B-JIMWLK equations [18,[20][21][22][23][24][25][26][27], which describe evolution of various gluon operators supplemented with Wilson lines. These operators have very different behavior for small k T but coincide (or vanish very fast) at large k T Q s , i.e. in the linear regime.
The modern effective theory incorporating saturation is the Color Glass Condensate (CGC) theory (see e.g. [28]; for a comprehensive pedagogical review of high energy QCD see [29]). In this theory, the gluon operators coupling to various particle production processes are averaged over random color sources of a dense target. It is important to underline two basic aspects of this approach: i) the eikonal approximation is assumed (recently also the non-eikonal corrections have been studied [30]), ii) the gluon operators contain both the leading power contribution (leading twist) and subleading power corrections, containing in principle the genuine multi-parton operators. The latter lead to resummation of the multiple parton interactions (MPI).
In jet phenomenology at LHC we typically assume the hard scale µ is set up by the average transverse momentum p T of jets. We consider here not so hard jets, with p T above, say, 15-20 GeV so that we are still sensitive to saturation, at sufficiently large energy and forward rapidity. In that regime, many simplifications occur. First, if one is interested in the back-to-back jet region only (k T µ assuming the collinear projectile), the leading power extraction [31,32] leads to an effective TMD factorization with onshell partonic amplitudes and several small-x leading power TMD gluon distributions containing various Wilson line operators ensuring gauge invariance and resumming collinear gluons [33]. In case we are interested in full description of jet imbalance (thus any k T between Q s and µ), we can extend the above effective TMD formalism to incorporate the kinematic twist corrections. This is done via keeping the incoming gluon off-shell in the amplitude and assuring the gauge invariance [34][35][36], which is equivalent to Lipatov's vertices in quasi-multi-Regge kinematics [37]. Such formalism has been first developed for dijets in [3] and is referred to as small-x improved TMD factorization (ITMD). It is equivalent to CGC expressions for dilute-dense collisions [38] with all kinematic twist corrections isolated and resummed (while neglecting genuine MPI contributions) [39].
The extension of the ITMD formalism for three jets used in the present work can be constructed analogously, provided the TMD gluon distribution definitions are known. Among others, for that purpose the operator structures for three-and four-jet processes have been explicitly calculated in [40]. The precise factorization formula shall be given in the next section. In the end, let us note that at present, the ITMD formalism does not account for the linearly-polarized gluons in unpolarized target. In CGC theory, such a contribution is absent for massless two-particle production, but appears in heavy quark production [41] and will appear in higher multiplicity processes as it has been already observed in the correlation limit [42]. Therefore in what follows we shall call our extension of the ITMD formalism to multipartonic processes ITMD*, to indicate that it does not take linearly polarized gluons into account yet. The construction of the full ITMD framework is left for the future.

ITMD* for three jets
The generic formula for multiparticle production within the ITMD* approach has been given in [40] in terms of color-ordered amplitudes (see e.g. [43] for a review of the color decomposition technique). For a specific case of forward particle production, where a dilute proton p (probed at large x) collides with a dense target A (probed at small x) the generic formula reads: where dΓ n is the n-particle phase space, f a/p the collinear PDF for parton a, depending on longitudinal momentum fraction x 1 and factorization scale µ, b 1 , . . . , b n are various final state partons contributing to ag → b 1 . . . b n partonic sub-process. Further A is a vector of tree-level color-ordered amplitudes for given partonic sub-process, C is the color matrix and the symbol • is the Hadamard (element-wise) Finally the Φ is the matrix of unpolarized TMD gluon distributions in the color-ordered basis. Entries of this matrix consist in linear combinations of the basis TMD gluon distributions given below. For three and four jets they have been explicitly calculated in [40]. It has to be noted that the formalism is not restricted to a specific color representation, and explicit formulas for a particular one are given in Appendix A. The operator definitions of the basic unpolarized TMD gluon distributions are 2 : where F.T. stands for the Fourier transform The angle brackets represent the hadronic matrix element . . . = P | . . . |P , where P is the momentum of the hadron. FurtherF µν = F µν a T a , where T a are the color generators. We employ standard light-cone basis, with hadron traveling along the 'plus' direction. The fields are separated in the light-cone 'minus' and transverse directions: The two staple-like fundamental representation Wilson lines connecting the fields are where the square brackets are the straight segments of the Wilson link. The Wilson loop is just two staples glued together: TMD matrices for the trijet case for two different color representations can be found in [40] and in Appendix A. The gauge invariant color ordered amplitudes with off-shell inital state gluon can be calculated automatically at tree-level using the methods of [34][35][36]. In the following work we use two latter methods independently to cross check the results. They are independently implemented in two different Monte Carlo programs generating weighted or unweighted events according to the formula 1: KaTıe [44] and LxJet [45].
Let us note that the ITMD formalism uses the small-x limit of the basis TMD gluon distributions. In that limit they can be rewritten in terms of matrix elements of CGC-style infinite Wilson lines with fixed transverse positions. Indeed, in the limit x → 0 only transverse position survives in the definitions (2)- (11) and the x dependence emerges from the evolution in energy. Whilst the full evolution equation adequate for moderate and small x for all correlators is not known (see [46,47] for initial attempts), the high energy limit is well controlled by the B-JIMWLK equations. Since B-JIMWLK is an evolution of the functional representing random color configurations in target it can be used for any operator. Assuming common initial distribution for operators contributing to dijet production, the proof of principle was given in [41,48]. This type of calculation can be carried keeping the subleading 1/N c corrections, but so far no distributions have been produced that incorporate data driven input.
In the following work we follow another path, first employed in [49]. We are going to use the TMD distribution appearing in the inclusive DIS processes, the so-called dipole distribution (2). In particular, we shall use the TMD coming from the BK equation supplemented with subleading corrections following the KMS framework [50] and fitted to F 2 data [51]. This equation is actually more suitable for harder jets because of including the DGLAP and kinematical constraint contributions. Having the dipole gluon distribution, all other distributions appearing in the dijet production: F qg can be calculated in the mean field approximation often used in CGC theory and to leading number of colors, see [49] for details.
In the above setup, i.e. to leading number of colors and in the mean field approximation for the distribution of color sources in the target, the cross section for trijet production can be calculated using the same basis TMD distributions as for the dijet case.

Numerical results
Before we present our results for the cross section, let us first discuss in more detail the basic TMD gluon distributions F qg , as well as the Weizsäcker-Williams gluon distribution F (3) gg calculated in [52]. This will be necessary to properly interpret the results for the cross section, as the trijet topology probes the kinematic range so far unexplored in inclusive and dijet calculations.
The KS dipole gluon distribution F qg , has two trends: the dependence on k T below ∼ 1 GeV was approximated by a power-like falloff, whereas above that threshold the TMD is given by the solution KS gluon TMDs in proton F qg (1) F qg (2) F gg (1) F gg (2) F gg (3) F gg (6) 10 KS gluon TMDs in lead F qg (1) F qg (2) F gg (1) F gg (2) F gg (3) F gg (6) 10 F qg (2) F gg (1) F gg (2) F gg (3) F gg (6) 10 F qg (2) F gg (1) F gg (2) F gg (3) F gg (6)  of BK equation with subleading corrections. The saturation of the distribution is visible as the clear smooth maximum developing for k T 1 GeV and moving towards larger k T with decrease of x (solid red line in Fig. 1). The remaining gluon distributions obtained from F (1) qg have universal behavior at large k T -they decay like ∼ 1/k T . The exception is the F (2) gg which decays much faster, so that it does not contribute to the perturbative tail 1/k T . Further, it becomes negative for some values of k T and approaches zero from below. Therefore, we plot its absolute value on the logarithmic scale in Fig. 1. To see in greater detail the differences as we go from proton to lead we plot also the ratios of gluon densities, Fig. 2. Interestingly, whereas at very low x all the distributions in lead are suppressed as compared to the proton, for moderate x and large k T the ratios exceed one. The interpretation of this result is not obvious. One possible explanation is that for a given x, the gluon distribution in lead can be significantly broader compared to proton, and thus create an enhancement for large k T , being still suppressed for small k T values. Now we are ready to discuss our present results for the trijet cross section. The processes that we are after are the p-p and p-Pb collisions at √ s = 13 TeV per nucleon and with demand that the three    jets are produced in the forward rapidity window 3.2 < |η 1 , η 2 , η 3 | < 4.9. The jets are defined as the outgoing on-shell partons satisfying the ∆φ − ∆η cut with the radius R = 0.5. We demand the jets with the transverse momentum of at least p T > 20 GeV. We order the jets according to their p T , so that we can distinguish leading, sub-leading, and the soft jet.
We perform calculation with the ITMD* framework described in the preceding section. For the collinear parton distributions we use the CTEQ10NLO set [53] obtained from LHAPDF6 [54], and the Kutak-Sapeta (KS) dipole gluon distribution [51] to get the five TMD gluon distributions [49] needed for the ITMD* at leading number of colors and in the mean field approximation, as discussed in details in the previous Section and above. We consider the following partonic channels, for 5 flavors of quarks: where q and q are quarks with necessarily different flavor. We do not include the sub-process with incoming anti-quark as it gives negligible contribution in the forward jet production case. The off-shell gauge invariant amplitudes were obtained numerically for fixed helicity and summed over on the eventby-event basis.
In the following we are interested in the azimuthal angle distributions of the jets. Thus, we consider the differential cross sections as a function of: (i) the azimuthal angle between the leading and sub-leading jets ∆φ 12 , (ii) the azimuthal angle between the leading and the soft jet ∆φ 13 , (iii) the azimuthal angle between the plane spanned by the two leading jets and the soft jet ∆φ (12)3 .
In addition to the absolute differential cross sections, a very useful observable is the nuclear modification ratio that quantifies saturation effects. It is defined generically as where the numerator corresponds to an observable in pA collision and the denominator to an observable in pp collisions, scaled by number of nucleons A. The deviation from unity suggests emergence of novel effects as one goes from pp to pA. In our case it reflects the emergence of nonlinearities leading to the gluon saturation (for R < 1) and possibly to anti-shadowing (for R > 1) due to the broadening of the TMD gluon distributions used in the calculations. In Fig. 3 we plot nuclear modification ratio as a function of ∆φ 12 , ∆φ 13 and ∆φ (12)3 . We consider the following scenarios: • the ITMD* case with the KS TMD gluon distributions with and without the x-dependent nuclear target area 3 , • the HEF case with nonlinear and linear KS dipole gluon distribution. The latter is the exact dilute limit of the ITMD* factorization and the CGC formalism.
Comparison of the above cases allows us to quantify the role of gauge links by comparing ITMD* with HEF when the gauge links are neglected but the nonlinearity is kept, as well as to quantify the combined effect of gauge links and nonlinearities by comparing ITMD* to completely dilute limit. From the panels in Fig. 3 we see that in all considered scenarios (except the fully dilute limit) a deviation from unity is clearly visible. Especially sensitive is ∆φ (12)3 , which shows a significant suppression in the back-to-back region, indicating strong saturation effects. Furthermore, we see that for some scenarios R pPb exceeds unity towards smaller values of azimuthal angles. We link this behavior to already discussed properties of ITMD* gluons where the ratio F P b /F p exceeds unity. In Figs. 4, 5, 6 we plot the absolute cross section, differential in ∆φ 12 , ∆φ 13 , and ∆φ (12)3 . We see that in addition to the already discussed suppression and enhancement of the p-Pb cross section (per nucleon), the normalisation of ITMD* is significantly larger than for HEF in the correlation region. We attribute this feature to a visibly different shape and larger normalization of the TMD gluons not present in the HEF formalism. Indeed, as seen from Figs. 1, 2 they start to dominate over the dipole gluon density as one enters the saturation region. It is clearly visible in Fig. 1 Figure 3: Nuclear modification ratio in azimuthal differences ∆φ 12 , ∆φ 13 and ∆φ (12)3 . The ITMD* calculation predicts less suppression in the corresponding back-to-back regions comparing to HEF and displays up to 10% of an enhancement of the p-Pb cross section away from the correlation limit (for the possible origin of this effect see the discussion in the main text).

Conclusions
In this work, we have presented calculations within an extension of small-x Improved TMD factorization, which was originally proposed in [3] for two final-state partons, to the case of three final-state partons, and designated it ITMD*. Our extension takes into account gauge invariant hard matrix elements involving off-shell eikonally-coupled initial state gluons, split into several nonequivalent color flows and corresponding TMD gluon distributions. At leading order, the formalism can be used to calculate trijet production in forward rapidity region in p-p and p-Pb collisions. The construction relies on an application of correlators obtained in [40] together with the sets of gauge invariant hard coefficient that match the color correlators.
Using the ITMD* factorization we calculated various azimuthal-angle-related observables, both for pp and p-Pb, as well as related nuclear modification ratios R pPb . We observe significant saturation effects, visible especially in R pPb as a function of the azimuthal angle difference between the plane spanned by two leading jets and the third jet. In addition, our results show that there is a significant difference between results obtained using the ITMD* and the standard k T -factorization/high energy factorization (HEF). First of all, the ITMD* results give higher cross sections in the correlation region compared to HEF, which is visible in all absolute cross section plots. Since the differences are rather large we expect them to be a good discriminator of theoretical frameworks. Secondly, the ITMD* result is not bound to unity, which for the present set of TMD gluon distributions is a consequence of broadening of the k T distribution for large k T and moderate x.
In the end, let us stress that ITMD factorization is a consequence of saturation effects. If there is no saturation, then the ITMD cross section formula reduces to HEF and all the TMD distributions reduce to a single gluon density being a solution to the linear evolution equation. In that sense, the differences between the ITMD* and HEF frameworks we observe reflect the proper account for saturation effects, and thus provide an excellent discrimination tool.

Acknowledgements
The authors are grateful to C. Marquet for useful comments and H. Van Heavermaet, P. Van Mechelen, M. Pieters, P. Taels for discussions. KK and AvH thank for support by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 824093. PK acknowledges partial support by Narodowe Centrum Nauki grant DEC-2017/27/B/ST2/01985.

A ITMD factorization in the color connection representation
The factorization formula for hybrid k T -factorization involves integrals over kinematical variables and an integrand including the factors where F is the k T -dependent pdf, and |M| 2 is the matrix element of the hard scattering process. The essential difference in the ITMD factorization formula involves these two factors. In order to illustrate this, it is useful to represent the matrix element in the color connection representation of [55,56], because it leads to particularly transparent formulas for ITMD factorization. This is also the color representation employed in KaTıe. The matrix element involves a sum over all color degrees of freedom of the external particles in the hard process. For a hard process with n g external gluons and n q external quark-antiquark pairs, the squared amplitude implies where a 1 , . . . , a ng are the adjoint color indices of the gluons, i 1 , . . . , i nq are the fundamental color indices of the quarks, and j 1 , . . . , j nq are those of the antiquarks. The color connection representation is obtained by introducing fundamental color indices for the gluons through contracting every adjoint index a with ( √ 2 T a ) k l , where the T a are the SU (N c ) generators, so Because of the identity the matrix element can now be written as ..,kn g l1,...,ln g i1,...,in q j1,...,jn q In [56] it is explained howM can be calculated directly with adjusted Feynman rules regarding color. Notice that the unavoidable "1/N c -correction" is in the quark-gluon vertex rather than in the gluon propagator like in [57], avoiding the need for projectors in Eq. (23). From now on we will use the same symbol i for color-indices of gluons and quarks, and j for anticolor indices for gluons and anti-quarks, and write n = n g + n q . The scattering amplitudeM can be decomposed into color factors and partial amplitudes following were S n is the group of all permutations of (1, 2, . . . , n), and where the partial amplitudes A σ do not depend on color, but may include factors of 1/N c . If the scattering process does not involve quarkantiquark pairs, then A σ actually vanishes for many permutations, which can be expressed explicitly with formula (3) in [57]. The formula above, however, holds for any process. Inserting Eq. (24) into Eq. (23), the matrix element can be expressed in terms of partial amplitudes via a color matrix with C στ = i1,...,in j1,...,jn It is not difficult to see that each entry of the matrix C στ consists of a single power of N c . Despite the fact that one can simply ignore vanishing partial amplitudes, the formula is still not optimal from the point of view of computational efficiency, since the partial amplitudes are linearly dependent. This is exploited in the color representation of [58], leading to the smallest possible color matrices, but with more complicated entries consisting of polynomials in N c . Below we will stick to the color representation with the big matrices with simple entries.
Let us assume that the off-shell gluon is the one carrying label number 1. In ITMD factorization, Eq. (19) is replaced with where the symbols U [λ] denote the Wilson lines of Eq. (14), with λ = ± depending on whether the parton whose color it connects is incoming or outgoing, andF + is the field strength. The double brackets represent both the hadronic matrix element and the Fourier transform. Notice that the formula does not distinguish between gluons and quark-antiquark pairs, and that compared to the formulas in Table 2  Table 1: TMD-valued color matrix and vector of partial amplitudes for the processes g * 1 q 2 → q 4q 3 q 5 and g * 1 q 2 → q 4q3 q 5 . The 6 partial amplitudes are explicitly labeled with their associated permutation. The logic in the enumation of the partons is that gluons come first, and then anti-quarks, where initial-state quarks count as negative-energy antiquarks. For clarity, all 6 partial amplitudes are included, also the non-contributing ones.  Table 2: Representation of the TMD-valued color matrix and vector of partial amplitudes for the processes g * 1 q 4 → g 2 g 3 q 5 . Each pair ij of intergers represents N i c F (j) qg , and a single 0 means that the entry vanishes. The last column gives the permutation associated with the partial amplitudes.
of [40], the "1/N c -terms" are absent. Now, we can insert Eq. (24) again, and find that the right-hand side of Eq. (27) can be written as where the entries of the "TMD-valued" color matrixĈ στ consist exactly of a single power of N c times one of the 10 F-functions from Eq. (2) to Eq. (11). The relevant matrices for 3-jet production are given in Table 1, Table 2, Table 3, and Table 4.  Table 3: Representation of the TMD-valued color matrix and vector of partial amplitudes for the processes g * 1 g 2 → g 3q4 q 5 . Each pair ij of intergers represents N i c F (j) gg , and a single 0 means that the entry vanishes. The last column gives the permutation associated with the partial amplitudes.  Table 4: Representation of the TMD-valued color matrix and vector of partial amplitudes for the processes g * 1 g 2 → g 3 g 4 g 5 . Each pair ij of intergers represents N i c F (j) gg , and a single 0 means that the entry vanishes. The last column gives the permutation associated with the partial amplitudes.