QCD resummation of dijet azimuthal decorrelations in pp and pA collisions

We study the azimuthal angular decorrelations of dijet production in both proton-proton (pp) and proton-nucleus (pA) collisions. By utilizing soft-collinear effective theory, we establish the factorization and resummation formalism at the next-to-leading logarithmic accuracy for the azimuthal angular decorrelations in the back-to-back limit in pp collisions. We propose an approach where the nuclear modifications to dijet production in pA collisions are accounted for in the nuclear modified transverse momentum dependent parton distribution functions (nTMDPDFs), which contain both collinear and transverse dynamics. This approach naturally generalizes the well-established formalism related to the nuclear modified collinear parton distribution functions (nPDFs). We demonstrate strong consistency between our methodology and the CMS measurements in both pp and pA collisions, and make predictions for dijet production in the forward rapidity region in pA collisions at LHC kinematics and for mid-rapidity kinematics at sPHENIX. Throughout this paper, we focus on the application of this formalism to a simultaneous fit to both collinear and transverse momentum dependent contributions to the transverse momentum dependent distributions.


Introduction
The investigation of high-energy proton-proton (pp) and proton-nucleus (pA) collisions is a crucial area of study in particle and nuclear physics, as it provides valuable insight into the fundamental structure of matter and the strong interaction among its constituents [1][2][3][4].Jet production is a crucial observable in these collisions, where collimated sprays of particles produced by the strong force, described by quantum chromodynamics (QCD), are observed.One of the key features of jet production in proton-proton and proton-nucleus collisions is the azimuthal angular distribution, or the difference in the azimuthal angle between the two jets.In the perturbative region, this decorrelation is a result of emissions from both the initial and final states that can alter the direction of the jets.The study of azimuthal decorrelation is critical for a deeper understanding of QCD jets and for testing QCD predictions and searching for new physics.
When one studies the dijet pseudorapidity spectrum while integrating over the full range of the azimuthal angle, the observable can be studied within the usual collinear factorization [5] and such a pseudorapidity spectrum is directly sensitive to the collinear parton distribution functions (PDFs), allowing us to constrain longitudinal motion of partons inside a free nucleon [6][7][8].When going from pp to pA collisions, there have been two approaches to deal with the nuclear modification [1], especially at the kinematic region where one probes the small-x parton physics.One is a DGLAP-based approach, while the other one is the saturation-based or color glass condensate (CGC) approach.In the DGLAP-based approach, one replaces the usual proton PDFs with the nuclear modified PDFs (nPDFs) [9][10][11][12][13][14] and follows the exact same collinear factorization.In this approach, the nuclear modification is included in the parameterization of the initial conditions for the DGLAP evolution of the nPDFs.On the other hand, in the saturation/CGC approach, gluon mergers and interactions dynamically lead to the nonlinear BK-JIMWLK evolution equations [15][16][17][18][19][20][21].For the theoretical formalism of the dijet production in the CGC framework, see for example Refs.[22,23].See also other work [24][25][26][27][28] along this direction.
Alternatively, when one studies more differential dijet observables, e.g.dijet azimuthal decorrelation, the conventional pQCD collinear factorization could be impaired.In the nearly back-to-back region where δϕ = π − ∆ϕ → 0, the perturbative expansion of the azimuthal angle decorrelation diverges due to logarithmic singularities at δϕ → 0 [29,30].The pioneering work in this field has highlighted the necessity of all-order resummation for accurately describing hadronic radiation, leading to a TMD-like factorization as shown below.This conclusion has been supported by numerous studies that have performed all-order resummation for various processes .In Fig. 1 we depict this back-to-back configuration for a narrow jet radius (R ≪ 1), where R is the radius of the jet.Fortunately, the azimuthal decorrelation of QCD jets in the nearly back-to-back region is sensitive to the intrinsic motion of the bound partons, allowing us to perform three-dimensional (3D) quantum imaging of the proton at high-energy facilities such as the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC).This three dimensional structure is encoded in the transverse momentum dependent parton distribution functions (TMD-PDFs), which contain both collinear and transverse momentum degrees of freedom.
While studying the nuclear modification to the inclusive dijet pseudorapidity spectrum in pA collisions, in the DGLAP-based approach, one encodes nuclear modification inside the nPDFs within the collinear factorization formalism.The natural question is how one handles the nuclear modification of the dijet production in the nearly back-to-back region when going from pp to pA collisions.As a natural generalization, we could encode nuclear modification of back-to-back dijet production inside nuclear modified TMDPDFs (nTMD-PDFs) within the TMD-like factorization formalism.Following such an approach, a recent global extraction of nuclear-modified TMDPDFs has successfully described world data for semi-inclusive electron-nucleus deep inelastic scattering and Drell-Yan processes in protonnucleus collisions in Ref. [55].Furthermore, an independent cross check of this analysis was performed in Ref. [56], verifying the results of Ref. [55].However, the applicability of nT-MDPDFs to other processes, such as dijet production, is yet to be determined.Finally, the study of QCD jet production in forward rapidity regions where one probes small-x parton dynamics is crucial for investigating the phenomenon of gluon saturation or CGC.Just like nPDF vs CGC approaches, to confirm saturation effects, it is important to have a proper understanding of the impact of nTMDPDFs vs CGC approaches in the back-to-back dijet production.For recent studies that deal with the back-to-back dijet production within the CGC formalism, see for example Refs.[57][58][59].
Experimental measurements of the azimuthal angular decorrelations in proton-proton and proton-lead (pPb) collisions at the LHC were performed in [60,61], respectively; while in [61,62] the integrated dijet azimuthal angle decorrelation in the region ∆ϕ > 2π/3 was measured.The first phenomenological studies of these data have been used to further constrain the nuclear modified collinear PDFs, see for instance in [12,63,64], by ap- proximating the integrated azimuthal angular decorrelations with the dijet pseudorapidity spectrum within a next-to-leading order (NLO) collinear factorization formalism.However, in the back-to-back region, which is encapsulated by ∆ϕ > 2π/3, the TMD effects, such as non-perturbative corrections and resummation can also be explored.Due to the sensitivity of these data to both collinear and transverse momentum contributions, these data can serve as a window into a simultaneous extraction of both collinear and transverse momentum effects in bound nucleons inside the heavy nucleus, which has so far not been performed.
In this study, we investigate the azimuthal angular decorrelation of dijet production in proton-proton collisions using the soft-collinear effective theory (SCET) framework [65][66][67][68][69].The utilization of the SCET framework enables us to perform QCD resummation of the large logarithmic terms in the azimuthal angle and jet radius at next-to-leading logarithmic (NLL) accuracy.Additionally, we examine the effects of nuclear modification on the azimuthal angular distribution in proton-nucleus collisions through the incorporation of nTMDPDFs and comment on the implications of our formalism to measuring nTMDPDFs as well as understanding nuclear modification of both collinear and transverse motions of the partons inside the nucleus.
Two predominant approaches are typically utilized for calculating the resummation formula in azimuthal decorrelation, known as the indirect [32] and the direct [29] methods.The indirect strategy focuses on the extraction of an all-order factorization and resummation formula for the two-dimensional transverse momentum imbalance q T of dijet pairs and the subsequent development of the azimuthal decorrelation ∆ϕ distribution originating from the q T distribution.In contrast, the direct method underpins the derivation of a factorization formula for the azimuthal angular distribution in the back-to-back limit, followed by the direct computation of all-order resummation results.While the association between these two methods is explicit for Drell-Yan-like procedures, it becomes increasingly intricate for processes implicating jet production, necessitating the resummation of sizable logarithms from final-state QCD radiation.Historically, it has been demonstrated that the indirect method could induce divergences in the azimuthal integral for a narrow jet radius [37,38].To mitigate these issues, various regularization schemes have been recommended [37,38,46].In this study, to evade such complexities, we have opted for the application of the direct method.
The rest of this paper is organized as follows.In section 2 we first discuss the factorization and resummation formula for nearly back-to-back dijet production in proton-proton collisions.Then we present the nuclear modified resummation formula in proton-nucleus collisions.In sub-section 3.1, we provide information for the numerical parameterization of the non-perturbative physics as well as the non-global logarithms (NGLs).We present the numerical results using the theoretical formula, enumerate all theoretical uncertainties and compare our predictions with the LHC experimental data in sub-section 3.2.We also make predictions for the azimuthal decorrelation of dijet production at the LHC, as well as for the sPHENIX kinematics region at the RHIC.We summarize our paper in section 4. The details of anomalous dimensions are provided in the appendix.

Factorization and resummation formula
In this section, we present our factorization and resummation formalism for the azimuthal decorrelation of dijet production in pp and pA collisions in the back-to-back limit.

Factorization in SCET for pp collisions
In the back-to-back limit and with the narrow jet approximation, the QCD modes which contribute to the dijet cross section are given by hard : p µ h ∼ p T (1, 1, 1), (2.1) ) ) where the momentum p µ is expressed in light-cone coordinates as and n µ i are light-like vectors associated with the initial-state proton beams (n a,b ) or finalstate jets (n c,d ).The n a,b -collinear, n c,d -collinear-soft and soft modes all have the same invariant mass and will result in rapidity divergences in the factorization formula.We address these divergences using the standard Collins-Soper-Sterman (CSS) treatment [70,71] and collinear anomaly [72,73] method, as explained in the next subsection.The contribution from the Glauber modes, which would result in the breaking of TMD factorization [74][75][76][77], is neglected in this study.The magnitude of factorization breaking effects from the Glauber mode can be explored by comparing theoretical predictions with future highprecision experimental data.
Based on the assumption of the above kinematic modes, we follow the standard steps in SCET [78][79][80] to obtain the following factorization formula1 where we have taken the short-hand The cross section is differential with respect to: the x component of the transverse momentum imbalance of the jet pair (|q x | = p T δϕ), the outgoing rapidities of jets c and d (y c,d ), the jet transverse momentum (p T ).In this expression, a, b, c, d represent parton flavors which are summed over in the cross section.The Kronecker delta symbol δ cd in the prefactor on the right side of this expression arises from the symmetry factor due to identical partons in the final state.Additionally, in this expression we introduced the partonic center-of-mass energy reads ŝ = x a x b s, and t = −x a p T √ se −yc , where s is the hadronic CM energy and x a and x b represent the Bjorken variables which are defined in terms of our phase space variables through the relations where E p is the energy of the incoming protons in the lab frame.The functions f unsub a,b/p represent the one-dimensional unsubtracted TMDPDFs for the incoming parton of flavor a, b [81].For these distributions, µ and ν are standard renormalization scale and rapidity scales, while ζ a,b represent the Collins-Soper parameters [82,83].
The function H ab→cd and S ab→cd are the hard and soft functions.In our formalism, we follow the work of [84,85] to organize the hard and soft functions into matrices, denoted by the bold characters.In this formalism, the IR divergent, UV finite scattering amplitudes for the 2 → 2 process can be written as vectors in color space where |C I ⟩ denote basis vectors in the color space while I is an index that runs over the dimensionality of the color space, which is determined purely through the species and the number of the external particles in the hard partonic process.The prefactors of the color basis vectors contain the kinematic contributions and the IR divergences of the amplitudes.
Following the work of [85], the basis vectors are absorbed into the soft sector.We now note that the integration of the virtual partons in the amplitudes of Eq. (2.9) contain interactions at the hard scale as well as interactions at scales associated with the IR modes in Eqs.(2.2), (2.3), (2.4), and (2.5).To define a purely hard scattering amplitude, one needs to subtract off the virtual loop contributions from these IR modes.As the virtual loop integrals of the IR modes are scaleless, this subtraction scheme swaps the IR divergences in the scattering amplitudes of Eq. (2.9) to UV ones.Thus we can define the purely hard scattering amplitudes through the subtraction where i runs over the IR modes in (2.2), (2.3), (2.4) and (2.5).The divergences entering into the hard scattering amplitude are now UV and can therefore be handled in a multiplicative renormalization procedure.Thus we can define UV subtracted amplitudes as where Z H is the hard multiplicative renormalization factor and is a matrix in color space.From this expression, the evolution of the subtracted scattering amplitudes is given by the expression where the hard anomalous dimension is defined as (2.13) In the following section, we will provide the hard anomalous dimension matrix, while we will further summarize the formalism in this section.
In SCET, the soft contributions enter as vacuum matrix elements.In our formalism, we define a b-space unsubtracted global soft function as with In this expression, b µ = (0, b, 0, 0), n µ i are the light-like vectors defined below Eq. (2.5), and T ( T) represents (anti-) time ordering.The soft Wilson line is given by where P denotes path ordering.We stress that, since we derive the factorization formalism in the direct method, the transverse vector b µ points along the x-direction, which is perpendicular to all vectors n a,b,c,d .This differs from the TMD soft function which was derived in [42], where the TMD factorization was derived for the two-dimensional transverse momentum imbalance of dijet pairs.As a result, the operator definition of the TMD soft function in this paper is different from that in [42].The soft function in Eq. (2.14) also enters into the factorization in the transverse energy-energy correlator event shape in [86].
To define the color matrix, we follow the work of Ref. [87] to absorb the color vectors into the soft function as where the SU (3) generators in the Wilson lines beyond tree level modify the color structure of the soft color matrices.
Aside from these complications associated with the hard and soft color matrices, to describe this observable, we must account for two final-state radiative effects.Firstly, in the narrow jet approximation (R ≪ 1), radiative corrections of the final-state partons are encoded in the jet and collinear-soft functions, J i and Scs i .The one loop exclusive jet function is well-known, see for instance [88], while the one-loop calculation of the collinear-soft function can be found in the appendix of [53].In addition to the standard ϵ divergences in dimensional regularization, the collinear-soft function that enters into our factorization also contains rapidity poles.We stress that these rapidity poles enter into the direct computation of the azimuthal angle decorrelation.However these poles do not enter into the collinear-soft function for the two dimension dijet transverse momentum imbalance in [42].Secondly, as the observable is insensitive to radiative emissions within the jet, this observable is non-global and is thus sensitive to NGLs [89].Such NGLs modify the factorization structure of the jet and collinear-soft function at two loops.The full factorization formula can be obtained by introducing the multi-Wilson structure in SCET [90,91].For simplicity, we do not write down the full formula in this paper, and in the resummation calculation, we use the fitting function [89,92] to include their contribution at the NLL accuracy.
After taking these effects into account, we note that the convolution in the cross section, C x , can be simplified by working in b-space, the conjugate space to q x .After performing the Fourier transform, the convolutional integral can be written as where the b-space functions are defined as After taking into consideration the simplification when working in b-space, the expression for the factorized cross section is given by the expression In the following sections, we will summarize the expressions for the evolution and resummation of each contribution in this cross section.

RG evolution and resummation formula
In the above subsection, we have obtained a factorization formula for azimuthal angular distribution in the joint back-to-back and small jet radius region.To achieve the resummation formula, one solves the RG equations for each of the ingredients in (2.21).In this section, we begin by performing resummation for pp scattering and then discuss our treatment for the pA scattering.
The hard functions for all 2 → 2 processes in massless QCD are given up to next-tonext-to-leading order (NNLO) in Ref. [93].To ensure consistency in the expressions for the hard anomalous dimensions between this study and our work, we choose to use the same color basis as this reference.Using these bases, the hard function satisfies the RG equation as where the anomalous dimension takes the form with C H = n q C F + n g C A and γ H = n q γ q + n g γ g .Here n q and n g indicate the number of quark and gluon, respectively.The matrix M reads where the dimensionless parameter r is defined as r = − t/ŝ.The expressions for M 1,2 can be found in Ref. [93].In this work, we consider QCD resummation at NLL accuracy, thus, we include the double logarithms anomalous dimension up to two-loop order and the single logarithms anomalous dimension up to one-loop order.The coefficients of all anomalous dimensions used in our calculation are given in the appendix A and we remark that the anomalous dimensions for quadrupole color and kinematic entanglement have been ignored in (2.23), since they contribute at three-loop order and beyond [94,95].Lastly, we remark that information associated with solving the RG equations in color space is provided in [85].
The jet functions in Eq. (2.21) fulfill the RG equation where the anomalous dimension of the jet is given by In this expression, C i = C F or C A is the Casimir of the parton i.It is worth noting that our analysis here does not account for the non-global structures in the factorization formula (2.21).As shown in [38], to obtain a complete description, the contribution of non-global structures must also be incorporated.In our current study, we do not take into consideration these structures in the factorization formula.However, the leading logarithmic (LL) NGLs are resummed by a fitting function, which is explained later in the paper.
In addition to the hard and jet function, all other terms in (2.21) also depend on the rapidity scale ν.For the TMDPDFs, we resum the large logarithms using the Collins-Soper equation.Specifically, in the Collins-11 treatment [82,83], the properly-defined TMDPDFs are obtained by absorbing the standard TMD soft function in the Dell-Yan process, Sab (b, µ, ν), and we have where C a denote the color of the incoming parton.In comparison to the two-dimensional transverse momentum resummation formula [42], the presence of rapidity divergence in the collinear-soft functions represents a new property.This divergence arises from the small jet approximation [53] and requires resummation of the corresponding rapidity logarithms.Two commonly used approaches to achieve this resummation are the rapidity RG [96,97] and collinear anomaly [72,73] framework.In this study, we choose to use the collinear anomaly framework.
In our study, we re-factorize the product of the global soft function and two collinearsoft functions.Using the collinear anomaly framework, we define a novel soft function W as where the rapidity logarithms arising from the narrow jet approximation in the collinearsoft functions are refactorized through the collinear anomaly exponent Notice that in this expression, we have subtracted the back-to-back soft function, which has already been included in the properly-defined TMDPDFs as in Eq. (2.26).This subtraction is required to avoid double counting of the soft modes in the final factorization formalism.Their renormalization group equations have the form as where Γ W is expressed as A rigorous test of our formalism is that we can obtain the RG invariance of the cross section as At NLL accuracy, the TMDPDF matches onto the collinear PDF through the relation where we have used the fact that the rapidity anomalous dimension vanishes at the scale µ b , κa (b, µ b ) = 0 at NLL accuracy and the f on the right hand side denotes the collinear PDF.Additionally, to circumvent the issue of the Landau pole in the large b limit, we have introduced the b * prescription that will be discussed in more detail in Sec. 3. Lastly in Eq. (2.35), we have introduced the non-perturbative Sudakov, which parameterizes the intrinsic motion of the bound partons and depends on the initial TMD scale Q 0 .
Combining the results for the hard, jet, TMDPDFs, and soft functions at NLL accuracy, our final resummed expression for azimuthal angular distribution is In this expression, the quantity λ K represents the eigenvalue of the matrices M 1,2 .In the small jet radius regime, the resummation of NGLs is achieved through a non-linear RG evolution between the jet and collinear-soft functions [38] that is contained in the U i NG functions.Lastly, µ h and µ j are the hard and jet scales which will be discussed in Sec.3.1.

Nuclear modified formalism for pA collisions
Having established the factorization and resummation for dijet production in pp collisions in the previous section, in this section we extend this formalism to incorporate the nuclear modifications in pA collisions.
As we mentioned in the Introduction, for observables that can be described by the collinear factorization formalism, a DGLAP-based approach can be used to deal with the nuclear modification when going from pp to pA collisions.In this approach, one assumes the same collinear factorization while replacing the proton PDFs with the nuclear modified PDFs [9][10][11][12][13].Now for the azimuthal decorrelation of dijet production in the nearly backto-back region, a TMD factorization and resummation in Eq. (2.6) is derived.Thus, as a natural generalization of the idea implemented in nPDFs, when going from pp to pA collisions, we assume that the same factorization and resummation formalism in Eq. (2.6) holds for pA collisions, while replacing the proton TMDPDF fb/p with the nuclear modified TMDPDF fb/A for the target nucleus.The nTMDPDF fb/A (x b , b, µ, ζ b,f ) contains the nuclear modification of both collinear (associated with x) and transverse (associated with b) motions for the partons inside the nucleus.Follow the assumptions made in Ref. [55], these nuclear modification will be absorbed into the non-perturbative parameterizations for the collinear PDF and the non-perturbative Sudakov.Thus under this assumption the nTMDPDF fb/A (x b , b, µ, ζ b,f ) can be matched onto the nPDF through the NLL relation Here, besides the collinear nPDF f b/A (x b , µ b * ), we have introduced the medium modified non-perturbative Sudakov S b,A NP (b, Q 0 , √ ŝ), whose parameterization will be discussed in the next section.Note that we only keep the leading power term in the OPE matching in Eq. (2.37) where the nTMDPDF is matched onto the collinear nPDF.In principle, there could be power corrections O b 2 Q 2 s (A) in the expansion which are associated with highertwist nuclear matrix elements [98].Here Q s (A) is a dynamical scale, often referred to as the saturation scale [99], associated with multiple scattering in the nuclear medium.We do not consider the effect of such power corrections in this paper.
With this replacement for nTMDPDF and following the same resummation procedure, the factorization and resummation formalism for back-to-back dijet production at NLL accuracy in pA collisions is given by (2.38) In the next section, we will discuss our parameterization for both the collinear nPDFs fb/A (x b , µ b * ) and our nuclear modified Sudakov factor S b,A NP (b, Q 0 , √ ŝ).
3 Numerical parameterization and results

Parameterization
To capture the resummation of the NGLs, we follow the prescription of Ref. [89] to parameterize the U function as where u = ln[α s (µ b * )/α s (µ j )]/β 0 , a = 0.85 C A , b = 0.86 C A and c = 1.33 [89].Since the factorized formula (2.6) involves two jet functions, the square of U NG is required to incorporate the NGL resummation associated with each jet.
In previous work on CSS resummation, the b * -prescription was introduced along with non-perturbative Sudakov factors, which were modeled through various functional forms and obtained by fitting to experimental data [100][101][102][103][104][105][106].In this work, we follow the standard b * -prescription where as in [71].Since we also need to study the impact of the nuclear modification on the azimuthal angular distribution in proton-nucleus collisions, we adopt the same functional form used in Refs.[102,107] which was employed in the extraction of nTMDPDFs [55].Specifically, the non-perturbative Sudakov factors in the last line of Eq. (2.36) are given by with g f 1 = 0.106 GeV −2 , g 2 = 0.84 and Q 2 0 = 2.4 GeV 2 .Finally, in our numerical calculations the intrinsic scales in the resummation formula (2.36) are chosen as To obtain numerical results for the pA collisions, we need a parameterization for nT-MDPDF in Eq. (2.37), which contains both the collinear and transverse motion for partons inside the nucleus.To describe the medium modifications to the collinear PDF, specifically f b/A (x b , µ b * ), we follow the parameterization in [55] to use the EPPS16 parameterization given in [63] while describing the collinear PDF for the proton, we use CT14nlo parameterization [108].On the other hand, for the nuclear modification to the transverse motion in nTMDPDFs, we follow the parameterization of Ref. [55] to have a nuclear modified Sudakov factor S b,A NP (b, Q 0 , √ ŝ).Specifically, we replace the g 1 parameter in Eq. (3.3), which accounts for the broadening effects of transverse momentum within the nucleus.Adopting the functional form obtained from the global extraction in [55], we take where g A 1 characterizes the transverse momentum width of partons inside the nucleus and is also proportional to the saturation scale in the small-x region [109].Thus the nuclear modified non-perturbative Sudakov factor is defined as

Numerical results
In this section, we present our numerical results for the pp and pA resummation formulas derived in the previous section.Specifically, we apply the theory formalisms in Eqs.(2.36) and (2.38) for pp and pA collisions, respectively and compare them with the existing experimental data.We discuss applications of this formalism to measuring nuclear modifications to collinear and transverse motions of partons in nTMDPDFs.We also provide predictions for the dijet production in the forward rapidity region in pPb collisions at the LHC, as well as in pAu collisions for the sPHENIX kinematics at the RHIC.A comprehensive investigation into the QCD resummation of azimuthal decorrelation in dijet production in pp collisions was carried out in [32] using the indirect method, as outlined in the introduction.The analysis successfully resummed the large logarithmic terms of azimuthal angle and jet radius at NLL and LL accuracy, respectively, while ignoring the contribution from NGLs.In our work, we present a resummation formula for azimuthal decorrelation in the direct method.This approach accounts for both the large logarithmic terms of the azimuthal angle and jet radius at NLL accuracy, including the contribution from NGLs.As a verification of the formula, we compare its theoretical predictions to measurements of dijet production in proton-proton collisions taken by the CMS collaboration at the LHC with √ s = 7 TeV, as presented in the left panel of Fig. 2. The QCD jets were reconstructed using the anti-k T algorithm [110] with a radius of R = 0.5 and the rapidities of each jet were limited to |y c,d | < 1.1.Additionally, to construct the denominator of the normalized ∆ϕ distribution, we use the LO expression for the cross section.The data, shown as black dots, covers five bins ranging from 80 GeV to 1 TeV for the jet transverse momentum p T .The theoretical results, displayed as lines of different colors, are found to agree well with the measurements in the back-to-back region across all p T bins.Besides, we also show the uncertainties from scale variations, which are given by the colored bands.
Here we vary the hard and jet scales by a factor of two around their default values as defined in Eq. (3.4), and the total uncertainty bands are obtained by the envelope of all the variations.Since the non-perturbative Sudakov factor in Eq. (3.3) is fitted at the canonical scale µ b * , we do not include uncertainties from its variations.It is noteworthy that the contribution of Glauber modes, which can potentially violate TMD factorization, is not considered in this analysis.Therefore, the magnitude of naive factorization breaking due to Glauber modes can be evaluated by comparing theoretical predictions with future high-precision experimental measurements.
On the right side of Fig. 2, we plot the azimuthal angle decorrelation in pPb collisions at √ s = 5.02 TeV from the CMS collaboration [61].The data is integrated within the region |y c,d | < 3 and the jets were reconstructed using an anti-k T algorithm with R = 0.3.In the theory calculation, we implement nTMDPDFs which encode the nuclear modification to the collinear (as in nPDFs f b/A (x, µ b * )) and transverse motion (as in the broadening parameter a N ) in Eqs.(2.37) and (3.6).The dashed blue theory curve is computed with the central fit of nPDFs in the EPPS16 parametrization and the broadening parameter a N in Eq. (3.5).The red band is the uncertainty from the nPDFs fit.Our calculations agree with the experimental data in the back-to-back region ∆ϕ ∼ π.We also observed in both plots of Fig. 2 that our theoretical prediction starts to deviate from the experimental data points away from the back-to-back region, i.e. when ∆ϕ moves away from π.This is expected since our formalism applies only to the resummation region.Such a discrepancy can be corrected by including the fixed order calculation for the dijet azimuthal angular decorrelation, see for instance [33].
In Fig. 3, we present a comparison between the NLL pQCD calculations of the dijet integrated angular decorrelation plotted as a function of the dijet pseudorapidity η = (y c + y d )/2 in pp collisions, respectively, and corresponding experimental measurement  taken by CMS [62].The data are categorized based on the transverse momentum (p T ) of the dijet system where the jet radius is R = 0.3.To enable an extensive comparison of the two datasets, the experimental measurements are superimposed onto the theoretical predictions, allowing us to evaluate the compatibility between the model and the experimental data.In the theoretical calculation, we integrate ∆ϕ from 2π/3 to π using (2.36) and to form the σ in the denominator of the integrated azimuthal angle decorrelation, we integrate over the pseudorapidity coverage of both jets following the experimental cuts.We observe that our theory calculations describes the experimental data quite well.In Fig. 4, we present our NLL calculation of the integrated angular decorrelation plotted as a function of the pseudorapidity in pPb collisions and the CMS experimental data in [62].To demonstrate the importance of nuclear modification to parton dynamics in the nucleus, we include a calculation where one only takes into account the isospin effect.In other words, going from pp to pA collisions, one only replaces the PDFs in the proton by the PDFs that include the isospin effect where Z is the atomic number of the nucleus while f i/p and f i/n denote the PDFs of the proton and neutron.We find that the calculations with the isospin effect alone undershoots the data rather significantly, especially in the mid-rapidity region 0 ≲ η ≲ 1 where an antishadowing effect is evident from the data [12,63].On the other hand, the central blue theory curve is computed with the central fit of nPDFs in the EPPS16 parametrization and the broadening parameter a N in Eq. (3.5).We further considered the uncertainty band associated with the collinear nPDFs as well as the broadening parameter a N .It is evident that our formalism with nuclear modification implemented in nTMDPDFs describe the CMS pPb collision data well though the size of the uncertainties from the broadening parameter a N is very small.This behavior is expected as the a N parameter acts to broaden the intrinsic width of the partons.At the large p T values of the CMS data, this broadening is small compared to the large transverse momentum that is generated from the resummation.Experimental data at smaller values of p T , which should be measurable at RHIC, will then depend more strongly on this parameter.However, the small dependence on a N indicates that both the integrated and unintegrated azimuthal angle decorrelation can be used to measure the collinear contribution to the nTMDPDFs.
To quantify the nuclear modification, we adopt the usual definition for the nuclear modification factor In Fig. 5, we present the nuclear modification factor R pA as a function of dijet rapidity η between our theory calculations and corresponding experimental data taken by the CMS collaboration at the LHC [62].In this plot, we have included a central curve as well as considered the uncertainty band associated with the nPDF and the broadening parameter  a N in nTMDPDFs.We have also included a prediction taking into account the isospin effect alone.The nuclear modification factor R pA with the isospin effect alone is almost unity as indicated by the dashed green curve.This is because the dijet production at this energy is mostly sensitive to the gluon distribution inside the nucleus and thus the isospin symmetry applied to u and d flavors does not play an important role here.On the other hand, we observe a strong consistency between the central curve of the NLL pQCD prediction with nTMDPDFs and the experimental data.However, we find that our calculations do not describe the strong suppression in the CMS data in the proton's forward region where η ≳ 2 and the probed parton momentum fraction x ∼ 10 −2 inside the nucleus.Since this modification in our nTMDPDFs formalism is mainly driven by the collinear nPDFs in the EPPS16 parametrization, as commented in [12,63], this remains an open question.As our formalism neglects all final-state interactions associated with Glauber interaction with the jets [111][112][113], we suspect that the cause of these discrepancies lies in these final-state effects [114].Addressing this discrepancy is vitally important for understanding the gluon distribution of the bound nucleons at this relatively small x region.
In the left and middle panels of Fig. 6 we present the results of our calculation for the azimuthal angular distribution in pp and pA collisions in forward rapidity regions at the ATLAS and ALICE kinematics at the LHC.In the right panel, we present the results of the decorrelation for the sPHENIX kinematics at the RHIC.In our study, we adopt the same kinematic cuts at the LHC as used in Ref. [115] and at the RHIC in Ref. [4], which are defined as follows: 1. 28 GeV < p T < 35 GeV and 2.7 < y * c,d < 4.0 for the FCal calorimeter of the ATLAS at the LHC , 2. p T > 10 GeV and 3.8 < y * c,d < 5.1 for the upgraded FoCal of the ALICE at the LHC , of nuclear modification effects.In future work, we see important applications of our formalism, e.g. in performing a simultaneous fit to both collinear and transverse momentum dependent contributions to the transverse momentum dependent distributions in nuclei.It would also be interesting to extend our results to other kinematic regions and incorporate the contributions from higher-order corrections, as well as to generalize our formalism to describe dijet production in the polarized scattering [116].

1 δϕFigure 1 .
Figure1.Definition of the azimuthal angular δϕ of dijet pair production in the x-y plane, where the transverse momentum of the leading jet j 1 is chosen to be aligned with the −y direction for convenience.

√ s = 7 5 |y | < 1 . 1 300Figure 2 .
Figure 2. Left: Comparison between theoretical calculations of the azimuthal decorrelation with the CMS data[60], where ∆ϕ is the difference in the azimuthal angle between two leading jets.The solid curves are the theoretical distributions, which are normalized by dividing LO cross section.The black dots are the CMS results, and the uncertainties of the data are smaller than the symbol size used in the plot.The colored bands indicate theoretical uncertainties from the variation of hard and jet scales.Right: A comparison of the dijet azimuthal angle decorrelation in pPb collisions from the CMS collaboration at the LHC[61].

55 < 3 ∆φFigure 3 .
Figure 3. Theoretical calculations for the dijet integrated angular decorrelation plotted as a function of the pseudorapidity η are compared with the CMS data[62] in proton-proton collisions for different kinematic cuts.The spectra were shifted by +0.465 to match the dijet pseudorapidity η range of the corresponding proton-lead collisions.

Figure 4 .
Figure 4. Theoretical calculations for dijet integrated angular decorrelation plotted as a function of the pseudorapidity η are compared with the CMS data [62] in proton-lead collisions for different kinematic cuts.

55 < 3 ∆φFigure 5 .
Figure 5. Theoretical calculations for the nuclear modification factor R pA plotted as a function of the pseudorapidity η are compared with the CMS data[62] for different jet transverse momentum cuts.

Figure 6 .
Figure 6.Top: The azimuthal angular distribution in pp (red curve) and pA (black curve) collisions for ATLAS (Left), ALICE (Middle), and sPHENIX (Right).In the lower panel, we plot the nuclear modification factor R pA .