Parton rescattering and gluon saturation in dijet production at EIC

Large angle gluon radiations induced by multiple parton scatterings contribute to dijet production in deeply inelastic scattering off a large nucleus at the Electron-Ion Collider. Within the generalized high-twist approach to multiple parton scattering, such contributions at the leading order in perturbative QCD and large Bjorken momentum fraction $x_B$ can be expressed as a convolution of the multiple parton scattering amplitudes and the transverse momentum dependent (TMD) two-parton correlation matrix elements. We study this medium-induced dijet spectrum and its azimuthal angle correlation under the approximation of small longitudinal momentum transfer in the secondary scattering and the factorization of two-parton correlation matrix elements as a product of quark and gluon TMD parton distribution function (PDF). Using a simple model for gluon saturation based on the parametrized gluon TMD PDF, we can calculate the $x_B$ and $Q^2$ dependence of the saturation scale and parton transport coefficient $\hat q$. Contributions to dijet cross section from double scattering are power-suppressed and only become sizable for mini-jets at small transverse momentum. We find that the total dijet correlation for these mini-jets, which also includes the contribution from single scattering, is sensitive to the transverse momentum broadening in the quark TMD PDF at large $x$ and saturation in the gluon TMD PDF at small $x$ inside the nucleus. The correlation from double scattering is also found to increase with the dijet rapidity gap and have a quadratic nuclear-size dependence because of the Landau-Pomeranchuk-Migdal (LPM) interference in gluon emission induced by multiple scattering. Experimental measurements of such unique features in the dijet correlation can shed light on the LPM interference in strong interaction and gluon saturation in large nuclei.


I. INTRODUCTION
Jet quenching in high-energy heavy-ion collisions caused by energy loss of energetic partons due to multiple parton scattering and induced gluon radiation in a hot quark-gluon plasma (QGP) has been the subject of intense theoretical and experimental studies [1][2][3][4][5][6][7] over the last several decades. Parton energy loss is dictated by the jet transport coefficientq, which is defined as the averaged transverse momentum transfer squared per unit length and related to the gluon density distribution in the medium [8][9][10][11]. Recent extraction of the temperature [12][13][14][15] and jet energy dependence [16][17][18][19] ofq through phenomenological studies of jet quenching can provide further insights into the properties of the QGP.
Multiple parton scattering and jet quenching can also occur in deeply inelastic scattering (DIS) off a large nucleus at the electron-ion collider (EIC) when the struck quark propagates inside the cold nuclear medium. Quark energy loss due to induced gluon radiation has been shown to cause the suppression of leading hadrons in semi-inclusive DIS (SIDIS) off nuclei [20][21][22][23][24][25][26]. The extracted jet transport coefficentq in cold nuclear matter from experimental data on the suppression of leading hadrons in SIDIS [27,28] and transverse momentum broadening [29] is about two orders of magnitude smaller * xnwang@lbl.gov † current address than that in the hot QGP in high-energy heavy-ion collisions. This is expected since the jet transport coefficent is directly related to the gluon number density of the medium which is much smaller in a confined cold nuclear matter than a deconfined hot QGP. It is, nevertheless, of great interest to explore the physics mechanism behind such a small value of jet transport coefficient in the cold nuclear matter. Several approaches have been developed to study parton energy loss in a medium with strong interaction over the last several decades.
These includes Baier-Dokshitzer-Mueller-Peigne-Schiff-Zakharov (BDMPS-Z) [8,9,30], Arnold-Moore-Yaffe (AMY) [31,32], Gyulassy-Levai-Vitev-Wiedemann (GLV-W) [33][34][35][36], high-twist (HT) [20,21,23,24,37] and more recent Soft Collinear Effective Theory gluon (SCET G ) [38][39][40] approaches with different approximations about the interaction between the propagating parton and the medium. Within the high-twist approach, one assumes the transverse momentum of the medium partons is small relative to that of the radiated gluon. Under such a collinear approximation an expansion in the medium parton transverse momentum is carried out and the collinear twist-four matrix elements can be factorized out [20,21,[41][42][43] in both DIS and Drell-Yan lepton pair production in p+A collisions. The collinear twist-four matrix elements can be related to the jet transport coefficient. Such an approach has been recently extended to a generalized high-twist approach [44][45][46] without collinear expansion in the transverse momentum of medium partons, by relaxing the approximation that medium parton transverse momentum is smaller than that of the radiated gluon. The final radiated gluon spectra can be expressed in terms of the transverse momentum dependent (TMD) jet transport coefficient which is just TMD gluon distribution density of the medium. In the soft radiation and small longitudinal momentum exchange limit, the final result becomes the same as that of the GLV approach under the first opacity approximation.
Within the generalized high-twist approach, one can also study multiple parton scattering in the opposite limit of the collinear exchanged gluon approximation. In this case, the initial transverse momentum from the medium gluon can be large and comparable to the transverse momentum of the final radiated gluon. Therefore, the final dijet spectra from multiple scattering will carry the information of the medium gluon TMD distribution. This will be the focus of this study. What's more, the hard splitting and inclusion of initial quark transverse momentum in the generalized High-Twist approach made the dijet a natural observable to entangle the medium gluon TMD distribution. The dijet with low transverse momentum and large angle capture the medium information better, since under our approximation final quark and gluon transverse momentum should be same order as initial quark and gluon transverse momentum. mong several approaches to induced gluon radiation due to multiple scattering in high-energy heavy-ion collisions, BDMPS-Z and GLV can be adapted for the application to DIS eA collisions. BDMPS-Z assumes multiple soft scattering while GLV considers the leading order in opacity expansion. Our approach is very similar to GLV and we can recover GLV results in soft and static scattering limit [46]. These two approaches are a better approximation in DIS since the number of multiple scatterings in DIS eA is very small. Our result should be valid at both small and large transverse momentum of the dijiet l T relative to the momentum imbalance k T . This is especially important because as we will show the nuclear enhanced nuclear modification of dijet cross section is measurable at l T a few GeV. At much larger l T , such nuclear modification is negligible for experimental observation.
Dijet production in proton-nucleus (p+A) collisions and electron-nucleus (e+A) DIS have been proposed to study the collinear nuclear parton distribution functions (PDF) [47,48], nuclear TMD PDF, multi-gluon correlations and gluon saturation in large nuclei [49][50][51][52][53][54][55]. In this paper, we will consider contributions to dijet production from multiple parton scattering in e+A DIS at the future EIC within the generalized high-twist approach. At large Bjorken x B , the momentum fraction of the struck quark from the virtual photon scattering in DIS, dijet production at the leading order (LO) is dominated by large angle gluon splitting from the struck quark after the photon-quark (single) scattering. Such large angle splitting can also be induced by a secondary scattering between the struck quark and a medium gluon from the nucleus target. Within the generalized high-twist approach, the spectrum of such medium-induced dijet production is related to the TMD PDF of soft gluons from the nuclear medium. The medium-induced dijet spectrum has a nuclear enhancement that is quadratic in the nuclear size R A and increases with the rapidity gap of the dijet, two unique features caused by the Landau-Pomeranchuk-Migdal (LPM) interference. We will study the sensitivity of this nuclear enhanced LO contribution to the medium gluon TMD PDF at small x, the scale of the gluon saturation and the transverse momentum broadening of the quark TMD PDF at large x B .
The remainder of this paper is organized as follows. In Section II, we calculate the dijet cross sections for both single and double scattering (medium-induced splitting) in e+A DIS at large x B . We will discuss the structure of the medium-induced splitting induced by different processes of multiple scattering with the LPM interference and how the final result depends on the two-parton correlation function which can be approximated as a product of quark and gluon TMD PDF. In Section III, we will introduce a simple model for implementing gluon saturation in the parameterized gluon TMD PDF inside the nucleus and its relation to the jet transport coefficientq. We will study the nuclear modifications to dijet spectra, from both single and double scattering through numerical calculations in Section IV. We will also explore the azimuthal angle, rapidity gap and nuclear size dependence of these nuclear modifications. A summary and an outlook are given in Section V.

II. DIJET PRODUCTION IN DIS
A. Single scattering We consider dijet production in DIS at large x B in this study. At LO, the dominant contribution in the single photon-quark scattering is from large angle gluon splitting from the struck quark as illustrated in Fig.1. The contribution from photon-gluon fusion γ * +g → q +q should be small because of the small gluon distribution at large x B . With the kinematics of dijet production we consider here, the initial state radiation will mainly contribute to the QCD evolution of the quark distribution in the nuclear target. Interference between initial and final state radiation is power-suppressed (k 2 T /Q 2 ) which we neglect in this study. Similar approximations are considered when we calculate contributions from multiple parton scattering.
The kinematics in the Breit frame (see Fig.1) by our convention are, where p is the momentum per nucleon of the nuclear target, q is the momentum of the virtual photon, l q is the momentum of the final quark, l is the momentum of the radiative gluon with the polarization vector (l) in an axial gauge A − = 0 and v ⊥ = l ⊥ + l q⊥ is the initial transverse momentum of the struck quark that determines the momentum imbalance of the final dijet in the single scattering. Under this convention and the collinear approximation Q l q⊥ , l ⊥ , v ⊥ , the hadronic tensor for gluon radiation from a single scattering in Fig.1 is, where, In the above equation, is the quark TMD PDF inside the nucleus target. This factorized form of the dijet cross section in DIS has been proven both within the traditional pQCD approach [56] and SCET approach [57,58] when the dijets are backto-back. Since we focus on nuclear modification of dijet cross section in the leading order in this study, we also assume that such factorization can be applied to the double hard scattering processes with large transverse momentum imbalance in the following section. The rigorous proof of the factorization and study of the associated soft functions arising from soft interactions is beyond the scope of this paper.
In the large x B region, we assume x B x L , x E . The quark TMD PDF can be approximated as q A (x B + x L + x E , l ⊥ + l q⊥ ) ≈ q A (x B , l ⊥ + l q⊥ ). Contracting the above hadronic tensor with the leptonic tensor, one can get the dijet cross section from gluon radiation associated with a single scattering in e+A DIS, where s = (p + q) 2 = 2p + q − is the center-of-mass energy squared of the photon-nucleon scattering, α em is the finestructure constant in quantum electrodynamics (QED) and e q is the electric charge of the quark. In the definition of the quark TMD PDF in Eq. (5), we have omitted the Wilson gauge link between two quark field operators due to multiple soft interaction between the quark and the nucleus target which guarantees the gauge invariance of the quark TMD PDF [59]. Such soft interactions in a nuclear target embedded in the quark TMD PDF leads to an effective transverse momentum broadening of the initial quark [11], where ρ A is the nuclear density distribution with the nor-malization, y − 0 is the light-cone coordinate in the Breit frame for the primary photon-quark scattering and b ⊥ is the impactparameter of the photon-nucleus interaction. The effective quark TMD PDF per nucleon inside a nucleus, with transverse momentum broadening from multiple soft collinear or eikonal interaction between the quark and the nucleus, is given at LO [11] by the convolution of the nucleon's quark TMD PDF in vacuum q 0 N (x, u ⊥ ) and a Gaussian broadening: The broadening is characterized by the quark transport coefficient [11], which is approximately proportional to the soft gluon dis- Here, α s is the strong coupling constant and φ N (x, k ⊥ ) is the gluon TMD PDF per nucleon inside the nucleus, We also take into account of the nuclear modification of the collinear or transverse-momentum-integrated PDF through a modification factor R q A (x, b ⊥ ) in Eq. (9). Such a modification factor has been parameterized through global fits to experimental data [60,61]. In the numerical studies in this paper we always consider x B = 0.2 and Q 2 = 200 GeV 2 . The modification factor R q A (x B , b ⊥ ) ≈ 1 in this kinematics.

B. Double scattering
In e+A DIS, as the struck quark from the photonquark interaction goes through multiple interactions inside the nucleus, a secondary hard scattering with another nucleon inside the nucleus can induce a gluon radiation leading to the medium-induced dijet production. We refer this secondary scattering as hard since the medium gluon carries small but finite momentum fraction as compared to the soft interactions that lead to the Wilson gauge link in the quark TMD PDF and the quark transverse momentum broadening. This process will interfere with the dijet production from the initial photon-quark scattering in which the final quark or radiated gluon scatters with another nucleon inside the nucleus, leading to the LPM interference. In Ref. [46], we have calculated the radiative gluon spectrum induced by double scattering in DIS including LPM interference. For this study, we extend the calculation to include the initial transverse momentum of the quark v ⊥ in the induced gluon spectra. In the large x B region, one can treat the longitudinal momentum transfer in the second scattering as small. When this momentum fraction comes from the medium gluon, the steep falling gluon TMD PDF effectively cuts off contributions when the momentum fraction is large. The hadronic tensor for gluon radiation induced by double scattering can be factorized in terms of the hard part of the photon-quark scattering H µν (0) , the transverse part of the induced gluon spectra per mean-free-path N g and the TMD quark-gluon correlation function T A qg : where the dijet momentum imbalance is l ⊥ + l q⊥ = v ⊥ + k ⊥ because of the momentum conservation, b ⊥ is the impact-parameter of the photon-nucleus collisions, y − 0 and y − 1 are the light-cone coordinates in the Breit frame of the primary photon-quark and the secondary quark-gluon scattering, respectively. The above result is the same as that in our last study [46] without the initial transverse momentum of the struck quark after the substitution The hadronic tensors for 23 different diagrams are listed in Appendix A. Take the central cut double scattering diagram in Fig. 2 as an example, its hadronic ten-sor has two terms, with two different TMD quark gluon correlation functions: The above result is obtained after integration over the light-cone momentum fractions of initial parton lines in the Feynman diagram x 0 , x 1 , x 2 and x 3 . One of these momentum fractions, x 3 , is fixed through momentum conservation x 1 + x 2 = x 0 + x 3 . The momentum fraction x 0 is fixed by the on-shell condition of the final quark which appears as a δ function in the hard part H0 µν . Two poles in the quark propagators effectively put two of the intermediate quark lines on shell through a contour integration, fixing both momentum fractions x 1 and x 2 , in terms of the momentum fractions defined as These typically are related to the momentum fractions that the initial partons need to carry in order to produce the final state gluon and quark with the given transverse momentum l T , k T and v T . The two contributions in Eq. (16) from Fig. 2 correspond to two different sets of poles in the contour integration: We refer readers to Refs. [21] and [46] for details about these contour integration and the different choices of the momentum fractions.
In the above result, we have assumed that the momentum fraction x B carried by the initial quark, usually referred to as the Bjorken variable, is much larger than other longitudinal momentum fractions carried by the medium gluon in different scattering amplitudes as shown in the Feynman diagrams in the Appendix A.
We consider the dominant case for a large nucleus target in which the primary photon-quark scattering happens at the light-cone coordinate (y − 0 , b ⊥ ) of one nucleon while the secondary quark-gluon scattering happens at (y − 1 , b ⊥ ) of another nucleon inside the nucleus in the Breit frame. We neglect the case that two scatterings happen inside the same nucleon which is not enhanced by the nucleus size. We further assume that the TMD quarkgluon correlation function in Eq. (13) can be factorized as a product of the initial quark TMD PDF of a nucleon and the soft gluon TMD PDF of another nucleon inside the nucleus [46,62], where q N (x, v ⊥ , b ⊥ ) is the effective quark TMD PDF per nucleon inside the nucleus as defined in Eq. (9) and φ N (x 2 , k ⊥ ) is the gluon TMD PDF per nucleon inside the nucleus as defined in Eq. (11). If one also considers nuclear effects on the medium gluon TMD distribution φ N (x 2 , k ⊥ ) such as gluon saturation, the gluon TMD dis-tribution also depends on the impact-parameter b ⊥ indirectly through the saturation scale as we will discuss later. This is implicitly assumed in the following calculations.
With the above factorized TMD quark and gluon correlation function, the medium-induced dijet cross section from double scattering in e+A DIS can be expressed as, The medium gluon in different processes of double scattering and induced radiation (see the Appendix A) carries different longitudinal momentum fractions as listed in Eq. (17). They are assumed to be small as compared to that of the primary struck quark x B . For we use x G to represent the small values of the longitudinal momentum fractions of the medium gluon in the medium-induced dijet cross section, Because of the interference among different radiation amplitudes, the final results for N qLPM g , N gLPM g , and N nonLPM g all vanish when k T → 0. This cancels the diverging behavior of the 1/k 2 T factor associated with the gluon TMD in the final dijet spectra in Eq. (19).
We have separated the transverse part of the induced gluon spectrum per mean-free-path or induced gluon spectrum rate N g into three different terms according to how the gluon is radiated in the symmetric central-cut diagrams. The symmetric central-cut diagrams for these three terms are illustrated in Figs Fig. 3(c) contains the final-state radiation from photon-quark scattering followed by gluon-gluon scattering and gluon radiation from the gluon propagator in the secondary quark-gluon scattering. In addition, one has to include all the interference among the three diagrams in Figs. 4(a), (b), and (c) and the interference between radiation amplitudes induced by single (vacuum radiation) and triple scattering (quark-photon scattering followed by two quark-gluon scatterings) in left/right-cut diagrams. Our separation of the three different contributions to the induced gluon spectrum rate and their physics interpretations are based on the central-cut diagrams.
The first term of the gluon spectrum rate N qLPM g in Eq. (20) comes from the final-state radiation of the photon-quark scattering in Fig. 4(b) and its interference with the amplitude of the final-state gluon radiation in Fig. 4(a). The LPM interference between gluon radiation induced by photon-quark scattering at y − 0 and quarkgluon scattering at y − 1 leads to the suppression factor which can be understood as the formation time of the radiated gluon from the struck quark in Fig. 4(b). The second term of the gluon spectrum rate N gLPM g in Eq. (21) comes from the gluon radiation off the gluon propagator in the quark-gluon scattering in Fig. 4(c) and its interference with the amplitudes of gluon radiation from the quark in Figs. 4(a) and (b). The formation time τ gf in the LPM suppression factor 1 − cos[(y − 1 − y − 0 )/τ gf ] depends on the transverse momentum k ⊥ of the medium gluon, which can be interpreted as the formation time of the intermediate gluon in Fig. 4(c). The third term of the induced gluon spectrum rate N nonLPM g in Eq. (22) comes from the final-state radiation in Fig. 4(a) minus the initial-state radiation of the quark-gluon scattering in Fig. 4(b), where the minus sign arises because of the space-like nature of the initial-state radiation. This term contains radiation amplitudes from the beginning and the end of multiple scatterings that do not participate in the LPM interference and is negligible when the number of scatterings is large. We keep this finite term since only two scatterings are considered in our study here.
In the soft radiation limit z → 1, both N nonLPM g and N qLPM g vanish, only N gLPM g remains and one recovers the GLV result [33,35] for induced gluon spectra in the leading opacity approximation. In our study of mediuminduced dijet production, we will keep all three contributions though N gLPM g is the most dominant contribution as we will show in the final numerical results.

III. JET TRANSPORT COEFFICIENT, TMD GLUON DISTRIBUTION AND SATURATION
Because of the momentum conservation, the transverse momentum imbalance of the final dijet is related to the transverse momentum of the initial quark and medium gluon, Therefore, the transverse momentum broadening in the nuclear quark TMD PDF will affect the LO pQCD result for the dijet spectrum due to single scattering in e+A DIS through the jet transport coefficientq F which is related to the gluon TMD PDF according to Eq. (10). The contribution from medium-induced dijet production due to double scattering, on the other hand, will depend on the jet transport coefficient through quark transverse momentum broadening as well as directly on the medium gluon TMD PDF φ N (x G , k ⊥ ) inside the nucleus. Therefore, the nuclear modification of the dijet spectrum in e+A DIS will be sensitive, both directly and indirectly, to the medium gluon TMD PDF inside the nucleus.

A. TMD PDF
In this study, we will use the TMDlib package [63,64] for the quark and gluon TMD PDF in nucleons and their scale evolution which are parameterized from global fits to experimental data. The gluon TMD distribution φ N (x, k ⊥ , µ 2 ) defined in our study is related to the TMD PDF xA(x, k ⊥ , µ 2 ) in TMDlib as, which are related to the collinear gluon PDF as This is an empirical formula for the relation between gluon TMD distribution and collinear distribution, since the connection between them is not very rigorous and fraught with theoretical subtleties [65][66][67]. Shown in Fig. 5 are the collinear gluon PDF from HERA2.0PDF parameterization "HERAPDF20-NLO-ALPHAS-118" in the LHAPDF package [68] (solid lines) and "PB-NLO-HERAI+II-2018-set1" [63,64] parameterization for gluon TMD PDF in the TMDlib package after integrating over the transverse momentum (dashed line) at different scales µ 2 = 100, 1000, 10000 GeV 2 . Numerically they agree with each other reasonably well. Note that the TMDlib package has a minimum starting evolution scale µ 2 0 for TMD PDF's. For the "PB-NLO-HERAI+II-2018-set1" that we use, µ 2 0 = 1.9 GeV 2 .

B. Gluon saturation
Parton distributions inside a nucleus will have nuclear modifications because of the multiple interaction in the initial state. These initial state interactions can lead to transverse momentum broadening [11], gluon saturation [69][70][71][72][73][74] and parton shadowing [75][76][77] as well as nuclear modification of the PDF at large x known as the EMC effect [78,79]. The nuclear modification of the quark TMD distribution in this study will be given by Eq. (9) with a transverse momentum broadening. The nuclear modification factor R q A (x B , b ⊥ ) [see Eq. (9)] for collinear quark PDF will be given by the EPS09 parameterization [61]. For gluon TMD PDF, we will consider the effect of gluon saturation.
To study the dijet spectrum in e+A DIS including the gluon splitting induced by double scattering, we propose a simple model for gluon saturation in the medium gluon TMD distribution, is the nucleon gluon TMD distribution in vacuum as given by the parameterizations in TMDlib and Q s is the saturation scale. This model is similar to the ansatz for unintegrated gluon distribution function with saturation in Refs. [80,81]. We will use the scale µ 2 = k 2 ⊥ in the above equation to include the scale dependence of the medium gluon TMD distribution that is involved in the secondary quark-gluon scattering. The typical small longitudinal momentum fraction carried by the medium gluon in the double scattering in e+A DIS can be written as, which is bounded by a lower limit x B Q 2 s /Q 2 when the transverse momentum becomes smaller than the saturation scale k ⊥ ≤ Q s .
Gluon saturation happens in scattering processes in a nucleus or nucleon when the gluon density at small x becomes high enough so that gluon fusion starts to overcome gluon splitting and the gluon density reaches a saturation limit. The saturation scale in this scenario can be related to the gluon density [74], where is the nuclear thickness function andq A is the gluon jet transport coefficient. The integration range for k ⊥ is bounded by the kinematic constraint ( We consider a hard-sphere model for the nuclear distribution in the rest frame, The nuclear thickness function is then, where R A = r 0 A 1/3 is the radius of the nucleus and r 0 = 1.12 fm. The solution for Q 2 s in Eq. (32) should scale approximately with the nuclear thickness function t A (b ⊥ ), which is proportional to the length of the quark propagation inside the nucleus. We can therefore approximately factor out the impact-parameter dependence of the saturation scale, Shown in Fig. 6 is the scaled saturation scale Q 2 s0 (x B , Q 2 ) as a function of Q 2 and x B from solving Eq. (32) with the nucleon gluon TMD PDF in vacuum given by the "PB-NLO-HERAI+II-2018-set1" in TMDlib. We can see in Fig. 6 there is a weak dependence on x B and scale Q 2 in the region of large x B ∼ 0.1 − 1 and moderate to large scale Q 2 ∼ 2 − 100 GeV 2 . For a large nucleus such as lead (Pb) (A = 208), the saturation scale in this region of kinematics is Q 2 s = Q 2 s0 A 1/3 = 0.25 − 1.29 GeV 2 in e+A DIS at zero impact-parameter b ⊥ = 0. We note that these values are much smaller than the saturation scale Q 2 s0 ∼ 1 GeV 2 [82] one expects at very small x B where the saturation scale in a large nucleus becomes large so that pQCD can be applied to calculate parton distributions with gluon saturation [72][73][74]. Within our simple model for gluon saturation in this study, the saturation scale Q 2 s0 can reach 1 GeV 2 at about x B ∼ 0.001 and Q 2 ∼ 100 GeV 2 .

C. Jet transport coefficient and p ⊥ broadening
Since we consider e+A DIS in the large x B region, the quark TMD PDF will not be affected by gluon saturation at small x. However, one should consider the effect of transverse momentum (p ⊥ ) broadening due to multiple soft interaction in addition to the EMC nuclear modification of quark PDF as we model in Eq. (9). This p ⊥ broadening is controlled by the quark transport coefficientq F .
The quark transport coefficient in cold nuclei has been extracted from phenomenological studies of p ⊥ broadening and suppression of the final-state hadrons in semiinclusive DIS (SIDIS). Analyses of the suppression of single inclusive hadrons in e+A SIDIS by Chang et al. [27] within the high-twist model of parton energy loss and medium modification of the fragmentation functions give the quark transport coefficient at the center of a cold nucleus in its rest frameq 0 F ≈ 0.02 GeV 2 /fm which was assumed to be independent of x B and Q 2 . Similar work by Li, Liu and Vitev [28] on the suppression of single inclusive hadrons in e+A SIDIS within the SCET G approach givesq 0 F ≈ 0.03 GeV 2 /fm. One can also extract the quark transport coefficientq F from the transverse momentum broadening of hadrons in e+A SIDIS. At LO and neglecting the radiative corrections [41][42][43][83][84][85][86][87], the average transverse momentum broadening of a quark can be related toq F according to Eq. (9) , in the hard-sphere model of the nuclear distribution, whereq 0 F is the quark transport coefficient at the center of a cold nucleus. The transverse momentum broadening of leading hadrons is ∆p 2 ⊥h eA = z 2 h ∆p 2 ⊥q eA and z h is the momentum fraction of hadrons in the quark fragmentation. A recent comprehensive analysis of the experimental data on the transverse momentum broadening of a variety of hadrons in e+A and Drell-Yan dilepton in p+A collisions [29] givesq 0 F ≈ 0.015 GeV 2 /fm in e+A DIS with a weak dependence on x B and the scale Q 2 in the range 0.05 < x B < 0.4 and 1 < Q 2 < 10 GeV 2 .
In general, one can define the TMD jet transport coefficient [10,46] for a parton with color representation R as,q where C R is the Casimir color factor, C F = (N 2 c −1)/2N c for a quark and C A = N c for a gluon. This definition of the gluon transport coefficientq A is the same as that for the gluon saturation scale Q 2 s in Eq. (32). After integratingq A over the parton propagation path, we get essentially the total transverse momentum broadening squared of a propagating gluon which is the same as the gluon saturation scale Q 2 s . Therefore, the gluon transport coefficientq 0 A at the center of a nucleus can be related to the scaled gluon saturation Q 2 s0 in Eq. (37), The quark transport coefficientq F is a factor C F /C A = 4/9 smaller than that of a gluonq A . Using this relation, one can also obtain the quark transport coefficient from the numerical solution to Eq. (32) for the gluon saturation scale as shown in Fig. 6. For x B = 0.1 − 0.4 and Q 2 = 2 − 6 GeV 2 , one getsq 0 F ≈ 0.013 − 0.023 GeV 2 /fm which is consistent with the valueq 0 F ≈ 0.015 GeV 2 /fm extracted from the transverse momentum broadening of hadrons in e+A SIDIS [29] andq 0 F ≈ 0.02 − 0.03 GeV 2 /fm from the suppression of single inclusive hadrons in e+A SIDIS [27,28]. Such momentum broadening can also be measured through azimuthal angle correlation between a single jet and the lepton in e+A DIS as proposed recently by Liu et al. [88]. For self-consistency, we will use the simple model for both the gluon TMD PDF in Eq. (30) with saturation andq 0 F as obtained from the gluon saturation scale in Eq. (40) for the quark transverse momentum broadening in Eq. (9) in our following calculation of dijet spectrum in e+A DIS. Using the transverse momentum broadening and gluon saturation for TMD PDF inside nuclei as modeled in the above section, we will evaluate the nuclear modification of the dijet spectrum at LO in e+A DIS. We will focus on the region of relative large x B ≥ 0.2, which is within the kinematic coverage in experiments at the proposed EIC [89,90] at the Brookhaven National Laboratory as shown in Fig. 7. We assume the proposed EIC will have FIG. 7: The Q 2 and x coverage of EIC with the electron beam energy E e = 10 GeV and the ion beam energy per nucleon E N = 100 GeV [89,90].
the electron beam energy E e = 10 GeV, the highest ion beam energy per nucleon E N = 100 GeV and the centerof-mass energy is √ s eN = 63.2 GeV. We will assume a typical set of kinematics x B = 0.2 and Q 2 = 200 GeV 2 for all the numerical calculations of the dijet spectrum in the following, unless specifically stated otherwise. We will work in the Breit frame of the initial quark in which the rapidity of the radiated gluon y l and the final quark y lq are, respectively, according to Eq. (1). To have two well separated jets in the dijet production, we require (y lq − y l ) 2 + ∆φ 2 > ∆R 2 , where ∆φ is the azimuthal angle difference between the two jets. Such a requirement will constrain the range of the transverse momenta l ⊥ and l q⊥ , momentum fraction z and azimuthal angle difference ∆φ.
In the calculation of the hadronic tensor, we have made the collinear approximation which requires Q 2 ( l ⊥ + l q⊥ ) 2 (Q 2 ≥ 4( l ⊥ + l q⊥ ) 2 ). This provides an additional constraint on the dijet azimuthal angle ∆φ and transverse momenta l ⊥ and l q⊥ where our calculations are applicable. Among the longitudinal momentum fractions of the medium gluon, the largest fraction ). This will provide an upper bound for the jet transverse momentum l 2 ⊥ ≤ (1 − z)zQ 2 /2 . In the transverse part of the induced gluon spectrum per mean-free-path N g in Eqs. (21)- (22), there are three collinear divergences in the denominators of the radiation amplitudes: The first and the third divergence are canceled by the LPM interference factors in N qLPM g and N gLPM g , respectively. The second divergence is regulated by the angular separation in the kinematics of the dijet. However, the first divergence still remains in N nonLPM g and N gLPM g . This divergence arises when the radiated [ Fig.4 (b)] or intermediate gluon [ Fig.4 (c)] with the momentum l = (0, (1−z)q − , (1−z) v ⊥ ) becomes collinear to the quark as the emitter with the momentum (0, q − , v ⊥ ). This divergence is normally absorbed into the renormalized TMD quark-gluon correlation function for The factorization scale µ f will serve to regularize the collinear divergence at l ⊥ − (1 − z) v ⊥ = 0 in the dijet spectrum.
We also include the running strong coupling constant α s (µ) for both gluon radiation and the secondary scattering in the dijet cross section and in the calculation of the saturation scale Q 2 s or the quark jet transport coefficientq F . We will set the scale µ = l ⊥ in the running coupling constant α s (µ) associated with the gluon radiation, µ = max(k ⊥ , Q s ) in the secondary scattering and the factorization scale µ f = 1 GeV. The final result is found not sensitive to µ f for µ f ≤ 1 GeV.

B. Dijet spectrum
The dijet spectrum in e+A DIS from single and double scattering in Sec.II [Eqs. (6) and (19)] can be expressed as where dσ 0 ep /dQ 2 ≡ 2πα 2 em /Q 4 , and dN dijet is the dijet spectrum per target nucleon that we will discuss in the remainder of this section. In e+p DIS, only single scattering without p ⊥ broadening contributes to the dijet production at LO. Fig. 8 are contributions to the dijet spectrum from single (dotted), double scattering (dashed) and their sum (solid) in e+Pb as well as e + p DIS (dotdashed) as a function of the azimuthal angle difference ∆φ (or the transverse momentum imbalance | l ⊥ + l q⊥ | on the top legend) between the two jets that have the same transverse momentum l ⊥ = l q⊥ . The dijet spectra from single scattering in both e+Pb and e + p DIS peak at the back-to-back direction (∆φ = π). The change in the single scattering from e+Pb to e + p is caused by the quark transverse momentum broadening due to multiple soft scattering inside the nucleus target.

Plotted in
For large values of the jet transverse momenta l ⊥ and l q⊥ , contributions from the double hard scattering are power-suppressed relative to the single scattering. Therefore they are much smaller than the contribution from single scattering. One has to go to smaller values of the jet transverse momentum in order to see the effect of double hard scattering. These two mini-jets don't peak at the back-to-back direction because of the secondary hard scattering. The double scattering contribution can become negative in some region because of the destructive interference between the initial and final-state radiation in N nonLPM g . However, the total dijet cross section (single + double scattering) (solid lines) is still positive because it is dominated by the contribution from single scattering.

C. Nuclear modification
To quantify the nuclear modification to the dijet spectrum in e+A DIS, we evaluate numerically the modification factor which is defined as the ratio of the dijet differential cross sections in e+A and e + p DIS. The modification factor has contributions from both single where dP ≡ dx B dQ 2 dzd 2 l ⊥ d 2 l q⊥ . We assume that we do not distinguish quark from gluon jet in experiments. The above dijet cross sectionsσ are the sum of the cross sections with the kinematics of quark and gluon exchanged, We will examine the azimuthal angle difference ∆φ, the momentum fraction z or rapidity gap |y l − y lq | and the nuclear size R A dependence of the nuclear modification factor in this section. Since the calculation of I S eA is straightforward, which contains just the effect of transverse momentum broadening of the initial quark on the process of single scattering, we will focus on the behavior of I D eA and its dependence on the azimuthal angle, transverse momentum imbalance, rapidity and nuclear size.
One can see that in the dijet spectrum from double scattering, the contribution containing N qLM g is negligible since it is suppressed by both the color factor and the LPM interference. The term containing N gLPM g is the most dominant in which the gluon is emitted from the gluon propagator. The contribution from N nonLPM g is small and finite, but its relative importance increases with large momentum fraction z or small rapidity gap |y l − y lq |. The magnitude of the dijet cross section from double scattering, however, decreases with the increase of l ⊥ when the transverse momenta of the initial quark and gluon become negligible and the dijet cross section from double scattering is power-suppressed relative to the LO cross section of single scattering. We therefore have to limit ourselves to mini-jets if we want to observe the contributions from double scattering in the dijet spectrum.    We have included both the transverse momentum broadening in the initial quark TMD PDF and the gluon saturation in the medium gluon TMD PDF in the calculation of the dijet cross section from double scattering according to Eq. (19). Inside a Pb nucleus, the saturation scale in our simple model is Q s (0 ⊥ ) = 1.2 GeV at zero impact-parameter. The corresponding total quark transverse momentum broadening squared is The combined effect of the initial quark transverse momentum broadening and gluon saturation inside a nucleus leads to the peak structure in the azimuthal angle distribution of the nuclear modification factor from double scattering at | l ⊥ + l q⊥ | ≈ 3 GeV/c as seen in Fig. 9.
The transverse momentum broadening in quark TMD PDF inside a nucleus can also lead to nuclear modifica-tion of the dijet spectrum from single hard scattering in e+A DIS. Since single scattering dominates in the total dijet cross section, such a nuclear modification purely due to quark transverse momentum broadening also dominates the nuclear modification of the total dijet spectrum. To illustrate the relative importance of nuclear modification in single and double scattering, we show in Fig. 10(a) the nuclear modification factor of the dijet angular distribution from single (dashed line) and single + double (solid line) scattering with l ⊥ = 2 GeV/c and z = 0.04. The quark TMD PDF [91] in a nucleon that we use has a Gaussian form at small transverse momentum and transits to a power-law form at large transverse momentum. The transverse momentum broadening in a nucleus target according to Eq. (9) will therefore suppress the nuclear modification factor for the dijet spec-trum from single scattering at v ⊥ = | l ⊥ + l q⊥ | ≈ 0, but enhance the modification factor at intermediate value of | l ⊥ + l q⊥ | reaching a peak before falling back asymptotically to 1 at large | l ⊥ + l q⊥ |. Such a nuclear modification of the dijet angular distribution in Fig. 10 is very similar to the typical Cronin effect of transverse momentum broadening in hadron spectra in p+A collisions [92]. For Q s (0 ⊥ ) = 1.2 GeV and ∆p 2 ⊥q = 0.48 GeV 2 in a Pb nucleus within our simple model for gluon saturation, the nuclear modification factor from single scattering peaks at | l ⊥ + l q⊥ | ≈ 1.3 GeV/c as seen in Fig. 10(a) (dashed line). From Fig. 9, we know that contributions from double scattering peak at | l ⊥ + l q⊥ | ≈ 3 GeV/c in these ranges of kinematics. Contributions from double scattering, therefore, further enhance the total dijet spectrum at large | l ⊥ + l q⊥ | (small azimuthal angle ∆φ) as we see in Fig. 10(a) (solid line).
Since the peak of the total nuclear modification (enhancement) is caused mainly by the quark transverse momentum broadening which in turn is determined by the gluon saturation scale in our model, it will shift to a larger value of | l ⊥ + l q⊥ | (smaller angle ∆φ) when the saturation scale is increased to Q s (0 ⊥ ) = 2.0 GeV as shown in Fig. 10(b) (dashed line), which is the solution to the self-consistent equation in Eq. (32) when the medium gluon density is artificially increased by a factor of 5. Correspondingly, the contribution to the dijet spectrum from double scattering is also bigger with its peak moving to a higher value of | l ⊥ + l q⊥ | due to the increased quark p ⊥ broadening and gluon saturation scale Q s as shown in Fig. 11. This would further enhance the total nuclear modification factor at large | l ⊥ + l q⊥ | as we see in Fig. 10(b) (solid line).

Rapidity gap dependence
In the Breit frame, the rapidity gap of the dijet is related to the longitudinal momentum fraction z of the radiated gluon according to Eq. (41), Shown in Fig. 12 are (a) different contributions to the dijet nuclear modification factor from double scattering and (b) the nuclear modification factor for the dijet cross section with (solid) and without (dashed) double scattering in e+Pb DIS as a function of the rapidity gap |y lq −y l | or momentum fraction z. We have chosen l ⊥ = l q⊥ = 2 GeV/c and ∆φ = 1 in these calculations where double scattering has the maximum contribution. The rapidity gap |y lq − y l | or momentum fraction z dependence of (a) the nuclear modification factor from double scattering I D eA (l ⊥ , l q⊥ , ∆φ, z) (b) the nuclear modification factor with (solid) I S+D eA (l ⊥ , l q⊥ , ∆φ, z) and without (dashed) double scattering I S eA (l ⊥ , l q⊥ , ∆φ, z) in e+Pb DIS.
Since the nuclear modification factor of the dijet spectrum due to single scattering is mainly caused by the quark transverse momentum broadening, it will be approximately independent of the longitudinal momentum fraction z or the rapidity gap |y lq −y l | as seen in Fig. 12 The contributions to the dijet spectrum from double scattering, however, have some unique z or rapidity gap |y lq − y l | dependence. Among the three contributions from double scattering in Eqs. (12)- (21), one can neglect N qLPM g [dot-dashed line in Fig. 12(a)] since it is suppressed by both a color factor and the LPM interference. For given transverse momentum l ⊥ (l q⊥ ) and azimuthal angle ∆φ, the small but finite contribution N nonLPM g [dashed line in Fig.12 (a)] vanishes as z → 0. The most dominant contribution is from N gLPM g [dotted line in Fig. 12(a)] which contains the LPM interference factor, with y − 01 = y − 1 − y − 0 . As z → 0 or at increasing rapidity gap, the formation time for the medium-induced gluon splitting becomes small such that the destructive LPM interference disappears, leading to an increased contribution due to incoherent dijet production induced by double scattering. This happens as the formation time τ gf becomes small such that 2R A /τ gf > ∼ 2π, or For the kinematics we use in Fig. 12, Q 2 = 200 GeV 2 , x B = 0.2, l ⊥ = 2 GeV/c and R A = 6.6 fm in a Pb nucleus, this leads to z incoh < ∼ 0.04 when medium-induced dijet production becomes incoherent. This is consistent with what we observe in Fig. 12. One can therefore consider the increase of the nuclear modification factor for dijet cross section with the rapidity gap, as illustrated by the solid line in Fig. 12(b), as a unique feature caused by the LPM interference in medium-induced dijet production.

Nuclei size dependence
In the single scattering process, the dijet cross section depends on the nuclear quark TMD PDF which is proportional to the atomic number A of the nucleus and the effective nucleon quark TMD PDF according to Eqs. (7) and (9). The corresponding nuclear modification will be determined by the quark transverse momentum broadening. At large transverse momentum v ⊥ = l ⊥ + l q⊥ , if the dijet spectrum has an effective power-law form 1/v 2n ⊥ , and the p ⊥ broadening is given by Eq. (38), the nuclear modification factor from single scattering has an enhancement, that is linear in the nuclear size R A . The numerical result on the nuclear size dependence of I S eA in Fig. 13(b) ( dashed line) indeed shows such an approximate linear dependence.
In addition to the nuclear quark TMD PDF, the dijet spectrum induced by double scattering is also proportional to a path integral of the differential induced gluon splitting rate (per mean-free-path) over the total length of the quark propagation, according to Eq. (19). The nuclear modification factor for dijet spectrum due to double scattering should then have a linear nuclear size R A dependence, from contribution N nonLPM g that does not have LPM interference. This length dependence comes from the integration over the position of the nucleon y 1 involved in the secondary scattering.
For contributions N qLPM g and N gLPM g that have the LPM interference, the integration over the position of the nucleon involved in the secondary scattering has to be weighted with the LPM interference factor. In the case when the photon-quark scattering occurs at the center of the nucleus, this leads to a nuclear size dependence, of the nuclear modification factor for the mediuminduced dijet spectrum, which is approximately quadratic when R A /τ f < 1. Shown in Fig. 13(a) are the nuclear modification factors of the dijet spectrum due to double scattering for l ⊥ = l q⊥ = 2 GeV/c, z = 0.04 at the azimuthal angle ∆φ = 1.45 where the medium-induced dijet spectrum is the largest (see Fig. 9). Again, the contribution from N gLMP g dominates. Including contributions from both single and double scattering, the nuclear size dependence of the modification factor for the total dijet spectrum in e+Pb DIS, shown as the solid line in Fig. 13(b), has a quadratic component due to double scattering on top of a linear dependence due to the transverse momentum broadening in the single scattering. Such a quadratic component in the nuclear size dependence of the nuclear modification factor is another unique feature due to the LPM interference in the dijet production induced by multiple scattering in the cold nuclear medium.

V. CONCLUSION AND DISCUSSIONS
In this study, we calculate the dijet spectrum at LO in pQCD within the framework of a generalized hightwist approach to multiple parton scattering in e+A DIS at EIC. We have specifically considered contributions to the dijet spectrum from both single and double scattering and examined in detail the dijet angular correlation which will be influenced by the transverse momentum broadening in the quark TMD PDF and saturation in the gluon TMD PDF in a large nucleus. We have employed a simple model in which we can determine the gluon saturation scale Q 2 s and the quark transport coefficientq F or transverse momentum broadening squared per unit length.
In the single hard scattering, the quark transverse momentum broadening dominates the nuclear modification of the dijet angular correlation which resembles the Cronin effect in the nuclear modification of the hadron transverse momentum spectra in p+A collisions. Corrections to the dijet cross section from double hard scattering is relatively small reaching to about 10% at jet transverse momentum l ⊥ ≈ 2 GeV/c and at the azimuthal angle where contributions from double scattering reach a peak value due to quark transverse momentum broadening and the saturation of medium gluon TMD PDF. Therefore the nuclear modification of the di-minijet angular correlation is sensitive to the gluon saturation scale in cold nuclei.
We have also examined the dependence of the nuclear modification of the dijet angular correlation on the dijet rapidity gap and on the nuclear size. We found that the dijet correlation increases with the dijet rapidity gap. The nuclear size dependence has a quadratic component on top of a linear dependence. Both of these two features are unique consequences of the LPM interference in the gluon splitting induced by double scattering.
We should note that our calculations are in LO of pQCD. Though we have used TMD PDFs that include a power-law tail at high transverse momentum due to QCD evolution, these are not exactly higher order corrections and resummation of soft gluon radiations which will lead to higher order corrections to the dijet angular correlation in both e + p and e+A DIS. Their nuclear modification and effects on the nuclear modification of the dijet spectra need more careful and quantitative investigations. Going into higher order corrections, one also needs to consider jet radius and jet algorithm dependence.
In principle, gluon saturation occurs when its coherence length 1/x G p + becomes larger than the nuclear size 2R A m N /p + or x G < x A ≡ 1/2m N R A . In this saturation limit, the effective gluon distribution per nucleon in Eq. (30) should be bounded from below by φ 0 N (x A , k ⊥ , µ 2 )/A 1/3 . To estimate the effect of such gluon saturation due to large coherence length, we can restrict the region of integration over the phase space in the dijet spectra to require x G > x A , effectively setting the saturated gluon distribution per nucleon to be zero when the coherence length is larger than the nuclear size. As shown in Appendix B, the modified dijet spectra are slightly smaller than that without the restriction x G > x A . Since the saturated gluon distribution below x G < x A should be non-zero, what is estimated in Appendix B is just a lower bound. For a more realistic estimate of this effect in the future, one should consider a gluon saturation model that takes into account of the coherence length of small-x gluons.
As we have shown by our numerical calculations, effects of double scattering and the LPM interference on the nuclear modification of dijets are only significant and measurable for minijets at EIC. Identification and reconstruction of these minijets is, however, rather challenging if not impossible. It is more straightforward to measure the correlation of dihadrons with moderately high transverse momentum. We expect the nuclear modification of dijet correlation due to multiple scatterings and induced gluon splittings should also be applicable to dihdadron correlation. This can be calculated by convoluting the dijet cross section with TMD fragmentation functions. This will be our next step in a follow-up study. ACKNOWLEDGEMENT We thank Z. B. Kang  In this appendix we list the hadronic tensors for all diagrams (Fig. 14 to Fig. 22) of dijet production induced by double scattering. According to Eqs. (12) and (18), these hadronic tensors can be expressed in the following form, In the following we list W according to the labeling of the corresponding cut diagrams. We also suppress the impactparameter b ⊥ dependence of the effective quark and gluon TMD PDF inside W. The definitions of momentum fractions x L , x S , ... are given in Eq. (17).
14: Central Cut 11 and Central Cut 22        From the azimuthal angle dependence in Fig. 23, the rapidity gap in Fig. 24 and nuclei size dependence in Fig. 25, we see that the restriction x G > x A reduces the dijet cross section from double parton scattering slightly. However, the rapidity gap dependence (Fig. 24) and nuclear size dependence (Fig. 25) show the same pattern when such coherence is not taken into account. The results here only illustrate the upper bound on the effect of the coherence length of small-x gluons since the saturated gluon distribution per nucleon should be non-zero, lying between the φ 0 N (x A , k ⊥ , µ 2 ) and the lower bound φ 0 N (x A , k ⊥ , µ 2 )/A 1/3 . | l ⊥ + l q⊥ | (GeV/c) The rapidity gap |y lq − y l | dependence of the nuclear modification factor from double scattering I D eA (l ⊥ , l q⊥ , ∆φ, z) (a) without (b) and with the constraint x G > x A , the nuclear modification factor I S+D eA (l ⊥ , l q⊥ , ∆φ, z) with (solid) and I S eA (l ⊥ , l q⊥ , ∆φ, z) without (dashed) double scattering (c) without (d) and with constraint x G > x A in e+Pb DIS. The nuclear size R A dependence of the nuclear modification factor from double scattering I D eA (l ⊥ , l q⊥ , ∆φ, z) (a) with (b) and without (b) constraint x G > x A , the nuclear modification factor I S+D eA (l ⊥ , l q⊥ , ∆φ, z) with (solid) and without (dashed) double scattering I S eA (l ⊥ , l q⊥ , ∆φ, z) without (c) and with (d) the constraint x G > x A in e+Pb DIS.