Transverse Energy-Energy Correlators in the Color-Glass Condensate at the Electron-Ion Collider

We investigate the transverse energy-energy correlators (TEEC) in the small-$x$ regime at the upcoming Electron-Ion Collider (EIC). Focusing on the back-to-back production of electron-hadron pairs in both $ep$ and $eA$ collisions, we establish a factorization theorem given in terms of the hard function, quark distributions, soft functions, and TEEC jet functions, where the gluon saturation effect is incorporated. Numerical results for TEEC in both $ep$ and $eA$ collisions are presented, together with the nuclear modification factor $R_A$. Our analysis reveals that TEEC observables in deep inelastic scattering provide a valuable approach for probing gluon saturation phenomena. Our findings underscore the significance of measuring TEEC at the EIC, emphasizing its efficacy in advancing our understanding of gluon saturation and nuclear modifications in high-energy collisions.

Numerous endeavors are dedicated toward the investigation of event-shape observables within the context of DIS.In this context, our focus is directed towards the Transverse Energy-Energy Correlation (TEEC) event shape observable in DIS.TEEC, as introduced in [22], originates as an extension of the Energy-Energy Correlation (EEC) [23,24] that was introduced for e + e − collisions to characterize global event shapes.In the environment of hadronic colliders, the event shape observable can be extended by considering the transverse energy of the hadrons [25,26].In the realm of DIS, the generalization of TEEC occurs through the application of the transverse energy correlation between the lepton and hadrons in the final state in the lab frame of lepton-proton collisions, which was initially conducted in Ref. [27].As demonstrated in Ref. [27], with the angle ϕ defined as the azimuthal angle difference between the produced electron and hadron transverse momentum, resummed predictions in the limit of back-to-back ϕ → π configurations can be obtained with high accuracy, allowing for reliable calculations of the distribution of ϕ across the entire range of [0, π].EEC and TEEC present a notable advantage in that the contribution from soft radiation is effectively suppressed due to its low energy.Consequently, the impact of hadronization effects is anticipated to be comparatively small when contrasted with other event-shape observables.Another advantage of the TEEC lies in the fact that the collision kinematics can be accurately reconstructed in the lab frame as pointed out in Ref. [28], and thus the TEEC can serve as great probes for the transverse-momentum dependent structures of the proton [27,29].In DIS, TEEC also offers a precise approach for determining the strong coupling, like the analyses in Refs.[30][31][32], and facilitates the exploration of nuclear dynamics as discussed in Refs.[33,34].
On the other hand, it has long been realized that the extracted parton distribution functions (PDFs) from experimental data, particularly the gluon distribution, exhibit a rapid increase as the partonic longitudinal momentum fraction, x, diminishes.The evolution of the gluon density at high energies, under the condition of fixed momentum transfer Q 2 , is encapsulated by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution equa-tion [35][36][37].The BFKL equation, a linear evolution equation, describes the evolution of the gluon distribution in terms of x.Its solution manifests a sharp increase as x decreases.Nonetheless, the gluon density is constrained from escalating indefinitely at high energies.In experimental observations, compelling evidence has emerged, especially at diminutive x values, indicating the presence of a distinct QCD regime known as the saturation regime.This regime eludes comprehensive explication through conventional linear QCD evolution frameworks [38][39][40][41].
Searching for the gluon saturation phenomenon [38,[41][42][43][44][45][46][47] is one of the scientific goals of the future EIC.The saturation physics refers to a phenomenon where the gluon density becomes so dominating that the interactions among gluons become significant, leading to a saturation of parton densities at small values of the partonic longitudinal momentum fraction x.Namely, this saturation occurs at high energy and small x, characterized by a saturation scale, denoted as Q s .Traditional linear QCD evolution equations, such as the BFKL equation, no longer accurately describe the dynamics in this regime [38,42].One then needs the non-linear extension of the BFKL equation, the Balitsky-Kovchegov (BK) equation [48,49].This non-linear dynamic phenomenon can be characterized better when a nuclear target is involved, wherein the interaction extends across a longitudinal distance approximately equal to or greater than the size of the nucleus.Under these conditions, the individual nucleons positioned at the same impact parameter become indistinguishable.Gluons originating from distinct nucleons have the potential to magnify the overall transverse gluon density by a factor of A 1/3 with A being the mass number of the target.Therefore, a substantial alteration in the TEEC is expected when the target hadron is substituted from a proton to a heavy nucleus like gold.Consequently, this novel observable, when explored at the forthcoming EIC, has the potential to provide further compelling evidence for parton saturation.
The rest of the paper is structured as follows.Section II provides the theoretical formalism for TEEC in DIS.We explain each component in the factorization, including the quark distribution in the small-x region and a detailed construction of the TEEC jet function.Section III presents our phenomenological study to demonstrate the potential of TEEC observables for probing gluon saturation and nuclear modification effects using ep/eA collisions.Finally, we conclude our work in Section IV.

II. THEORETICAL FORMALISM
In this section, following the theoretical formalism of TEEC in Deep Inelastic Scattering [27], we study the transverse energy-energy correlation between the lepton and hadrons in the final state: where the scattered electron and final hadron are produced in a back-to-back configuration in the transverse plane.The TEEC is illustrated in Fig. 1 and defined as: where the sum runs over all the hadrons in the final state, and we define the variable τ as: Here ϕ is the azimuthal angle between the final-state lepton e and hadron h as shown in the right panel of Fig. 1.We have also defined the angle which is a small angle under the back-to-back limit, ϕ → π.Correspondingly, we have τ ≪ 1.As we have mentioned in the Introduction, we analyze the event in the center-of-mass frame of the lepton and proton collisions, with the proton (or the nucleus) moving in the +z direction while the incoming lepton moving in the −z direction, as shown in Fig. 1.
The TMD factorization theorem for the TEEC observable in the back-to-back region (i.e.τ ≪ 1) is given by [27,50]: Even though the TEEC is the cross section weighted by the hadron momentum fraction as in Eq. ( 2), we abuse the notation a bit by still denoting it as dσ.Here y e and p e T are the rapidity and transverse momentum of the produced lepton in the laboratory frame with respect to the beam direction, and we take the outgoing lepton to lie along the y-axis.On the other hand, f (u) q x, b, µ, ζ/ν 2 is the "unsubtracted" TMD quark distribution, where b is the x-component of the b vector in the standard quark TMD distribution as probed e.g. in semi-inclusive DIS [51,52].In other words, we have b ≡ (b x , b y ) = (b, 0) and thus the integration limits are given by b ∈ (−∞, ∞) in the first line of Eq. ( 4).It is important to realize that the cross section is differential in variable τ (i.e.azimuthal angle ϕ), which is related to the x component of the transverse momentum of the final observed hadron, where z is the momentum fraction of the quark carried by the hadron fragmenting from it.Consequently, we have a one-dimensional Fourier transform, i.e.only the x component of the conjugated coordinate variable b is relevant.This has been derived clearly in [27,50,53].S nn h (b, µ, ν) is the soft function representing the contribution from soft gluon radiation, and H(Q, µ) is the hard function.At the same time, J q (b, µ, ζ ′ /ν 2 ) is the "unsubtracted" TEEC jet function, which has a close relation with the TMD fragmentation functions as given below.On the second line of Eq. ( 4), taking the advantage that the functions f (u) q , S nn h , and J (u) q are all even function of b as they depend on b 2 , we further simplify the integration to be in the region b ∈ (0, ∞).
Finally, the well-known prefactor σ 0 is the leadingorder (LO) partonic cross section for lepton-quark scattering where α em is the fine structure constant, s is the centerof-mass energy squared of the incoming lepton and the proton beam, Q 2 represents the photon virtuality.In the back-to-back lepton-hadron production region, the partonic Mandelstam variables ŝ, t and û are connected to the Bjorken-x and other kinematic variables: For convenience, we also list here the Bjorken x and inelasticity y written in terms of other kinematical variables of interest: where we have used the momentum conservation relation ŝ + t + û = 0.
In the following subsections, we will identify all the components in the factorization theorem as given in Eq. (4).

A. Quark distribution
In this subsection, we provide a short overview of TMD quark distribution and discuss its expansion in terms of gluon dipole distribution in the small-x limit.
For the "unsubtracted" TMD quark distribution f (u) q x, b, µ, ζ/ν 2 , we have the Collins-Soper scale ζ [51,54,55] and a rapidity scale ν [56].The rapidity divergence in f (u) q can be canceled by subtracting a square root of the standard soft function S nn (b, µ, ν) whose result at the next-to-leading order (NLO) is given by where µ b is defined as It is worth noting that in this work we have applied the 4 − 2ϵ space-time dimensions and the rapidity regulator η [56].As a consequence, we further defined the "subtracted" parton distribution f q (x, b, µ, ζ) without a rapidity divergence [55]: TMD evolution for the "subtracted" TMD quark distribution is governed by two equations, the Collins-Soper evolution associated with the Collins-Soper scale ζ [51,55] and the renormalization group equation related to the scale µ.They are given by where K(b, µ) denotes the Collins-Soper evolution kernel [51,55,57,58] and γ q µ α s (µ), ζ/µ 2 is given by: where Γ q cusp and γ q µ are the cusp and non-cusp anomalous dimensions.They can be perturbatively expanded as: Solving the renormalization group equations on ζ and µ and taking into account the non-perturbative contribution at the large b ≫ 1/Λ QCD region, we obtain the TMD quark distribution as where we evolve the TMD quark distribution Throughout this paper, we will work at the nextto-leading logarithmic (NLL) level, where we have On the other hand, S NP (b, Q 0 , ζ) is a non-perturbative Sudakov factor for the TMD quark distribution, see e.g.
In this work, in order to explore the gluon saturation, following Refs.[68,69], we expand this TMD quark distribution at the initial scale where S ⊥ is the averaged transverse area of the target hadron and S x (r) represents the dipole scattering matrix with the dipole transverse size r.We consider two different models for S x (r).The first is the Golec-Biernat-Wüsthoff (GBW) model [70,71] which can be written as: where the saturation scale Q s reads: The free parameters in this model are chosen as λ = 0.29, x 0 = 3 × 10 −4 and S ⊥ = 1/2 × 23 mb for proton targets following Ref.[70].The other model we consider is based on the McLerran-Venugopalan (MV) initial condition [45,46,72] which is then evolved with a running-coupling BK (rcBK) equation to smaller values in x.Specifically, we use the MV e initial condition [73]: and the rcBK equation: For the kernel K(r, r ′ ), we use the Balitsky prescription [74]: with the coordinate-space running coupling We shall call this the rcBK model.The values of the parameters for proton targets are taken from Ref. [73], with the transverse size being S ⊥ = σ 0 /2 in terms of the parameters presented there.We also note that these parameter values are very close to the more recent ones in Ref. [75] determined using Bayesian inference.
Finally, one has the following expression for TMD quark distributions in the CGC formalism, with f q (x, b, µ b * , µ 2 b * ) at the small-x region provided in Eq. ( 24) and the perturbative Sukadov factor given in Eq. ( 20).In comparison with the standard TMD quark distribution in Eq. ( 19), we ignore the non-perturbative Sudakov factor S NP (b, Q 0 , ζ).This is because in principle the small-x formula for the TMD quark distribution in Eq. ( 24) has already contained the non-perturbative contribution in the large-b region [69].

B. Hard and Soft functions
The hard function H(Q, µ) with the renormalized expression at the one-loop is given by [76][77][78]: The natural scale for the hard function is given by µ ∼ Q.
On the other hand, the soft function S nn h (b, µ, ν) in DIS for TEEC at the NLO is given by: where n • n h = 1 − tanh(y) with y the rapidity of the final-state hadron.This soft function can be related to the soft function for EEC in e + e − , namely the standard soft function S nn (b, µ, ν) in Eq. ( 12) by [27,53]: C. TEEC jet function and factorization In Eq. ( 4), the function denoted by J (u) q b, µ, ζ ′ /ν 2 is the unsubtracted TEEC jet function [79] which is related to the "unsubtracted" transverse-momentum-dependent fragmentation functions (TMD FFs) via: where the D (u) 1,h/q (z, b, µ, ζ ′ /ν 2 ) are the TMD FFs in the b-space.To simplify the notation, here we introduce the "subtracted" TMD FFs as: Using the results for the soft functions S nn (b, µ, ν) and S nn h (b, µ, ν) given in Eqs. ( 12) and ( 33), we find that the Collins-Soper scale for the "subtracted" TMD FFs D 1,h/q will be given by ζ Note that in the rapidity regulator we adopt [56], for TMD PDFs, the Collins-Soper scale is √ ζ/2 = xP + A , and for TMD FFs, one has √ ζ ′ /2 = P − h /z.Thus: Namely, we find that ζ ζ = Q 4 , and thus one can choose ζ = ζ = Q 2 as a natural scale choice for the TMDs involved in the factorization formalism.Subsequently, the corresponding "subtracted" TEEC jet function J q b, µ, ζ can be further written as: The TMD FFs D 1,h/q z, b, µ, ζ with QCD evolution is given by where the matching coefficients C i←q can be found in Refs.[59,[65][66][67].The corresponding non-perturbative Sudakov factor is given by: with g 2 = 0.84 and g D 1 = 0.042 GeV 2 [59,60].Plugging Eq. ( 39) into (38), one thus obtains a general form for the TEEC jet function.If it were not for the z-dependence in the non-perturbative Sudakov term exp −g D 1 b 2 /z 2 in Eq. ( 39), one could decouple the zintegral in Eq. ( 38) with the y-integral in Eq. ( 39): Here in the second line, we change the integration variable u = z/y, and in the third line, we apply the momentum sum rule, This result is consistent with [79].Unfortunately, the explicit z-dependence in the non-perturbative Sudakov factor S NP z, b, Q 0 , ζ makes the TEEC jet function more complicated.
To proceed, we choose the coefficient function at the leading order C i←q (z, b) = δ iq δ(1 − z) in Eq. (39), and thus the TEEC jet function J q in Eq. ( 38) can be written as: Next, to prepare for the phenomenological study, we proceed by specifying a model for the TEEC jet function.
Following [80], we perform a fit to obtain a simple form for the z-integrated expression.Specifically, we define: and use the DSS parameterization [81] for all the hadrons h = h + , h − , h 0 in the fit.We find that the following functional form works very well: The fitted parameters are given by g e 0 = 0.226 GeV 1/2 , g e 1 = 0.463 GeV, and g e 2 = 0.033 GeV 2 .This fitted result is slightly different from what was obtained in Ref. [80].Therefore, one has the TEEC jet function given by Eventually, one can write the factorization theorem in Eq. ( 4) in terms of "subtracted" quark distributions and TEEC jet functions as: In the phenomenological section below, we choose the As indicated in the introduction, when changing from ep to eA collisions, one takes nuclear modification effects into consideration and substitutes the saturation scale Q s in Eq. ( 25) by the nuclear saturation scale Q 2 s,A ∼ A 1/3 Q 2 s .More details about the numerical values of relevant parameters will be discussed in Section III.

III. PHENOMENOLOGY
In this section, we make numerical predictions for the TEEC at the future EIC for both ep and eA collisions.
With the factorization of the TEEC jet function given in Section II C, we are now ready to perform numerical predictions for the TEEC at the future EIC.We choose the highest center-of-mass energy √ s = 140 GeV for electron-proton collisions.We work in the frame where the proton is moving along the +z direction, and the electron moves along the −z direction.In order to probe the small-x region, we need to choose a proper lepton rapidity and transverse momentum.As an example, we choose y e = −2 and p e T = 2, 4, 6 GeV.This corresponds to the probed x values between 2 × 10 −3 and 8.5 × 10 −3 .In Fig. 2, we plot the TEEC as a function of τ for these three different p e T values.The solid curves are from the rcBK parameterizations while the dashed ones are based on the GBW model.The red curves are for p e T = 2 GeV, the blue ones are for p e T = 4 GeV, and the green ones are for p e T = 6 GeV.We find that the numerical results based on the rcBK and the GBW model for the TEEC observables can differ by a factor of two, especially at the small p e T = 2 GeV, indicating the TEEC at the EIC can be a good observable constraining the dipole gluon distribution.
To study the nuclear modification in e + A collisions in comparison with the e + p scatterings, we define the nuclear modification factor R A as follows: where A is the atomic mass of the nuclear target.Below, we choose the gold nucleus with A = 197.To go from the proton to the nucleus beam, we change the proton saturation scale to the nuclear saturation scale Q s,A (x) or Q 2 s,0,A for the GBW and the rcBK models, respectively: where Q s (x) and Q s,0 are the proton saturation scales for the GBW and the rcBK models, respectively.The parameter c is chosen in the range 0.5 < c < 1.0 [69,82].Correspondingly, we also change the active nuclear transverse area S ⊥ → S ⊥,A = 1/c × A 2/3 S ⊥ .Having the same scaling constant c for both the saturation scale and the transverse area can be motivated by the smooth nucleus approach [83], where in the dilute region we can integrate over the impact parameter and write future, we plan to implement the impact parameter dependence [83][84][85] directly inside the saturation formalism and thus provide more accurate predictions.
In Fig. 3, we plot the nuclear modification factor R A as a function of τ for p e T = 2 GeV (top panel), 4 GeV (middle panel), and 6 GeV (bottom panel).We choose √ s = 140 GeV and lepton rapidity y e = −2.The bands correspond to the uncertainty in the parameter 0.5 < c < 1.0.The red bands are for the GBW model, and the blue bands are for the rcBK calculations.It shows that nuclear modifications on the order of 15%−20% can be expected in the small τ region, for both the rcBK and the GBW model.On the other hand, the nuclear modification factor starts to approach 1 as the τ value increases.Such behavior is a manifestation of the cos(2b √ τ p e T ) modulation in Eq. ( 46).In the large τ region, the integration is dominated by the small-b region where the dipole size is small and thus the saturation effect is less important and one expects R A → 1.On the other hand, in the small τ region, one would receive more contribution from the larger dipole size (large b region) and correspondingly stronger nuclear modification.This indicates that the TEEC is a good observable for gluon saturation.

IV. CONCLUSIONS
In this paper, we explore the transverse energy-energy correlators in the small-x regime for the future EIC.For the production of electron-hadron pairs in the back-toback region in the transverse plane where the azimuthal angle difference ϕ → π between the final-state lepton and the hadron, we provide a factorization theorem that incorporates the gluon saturation effects.We present numerical results for TEEC in both e+p and e+A collisions, alongside evaluations of the nuclear modification factor R A .We find that the TEEC observables in e + p collisions are significantly influenced by different models of the dipole gluon distribution, emphasizing the potential of TEEC at the EIC as a robust observable for constraining the dipole gluon distribution in the small-x region.We introduce the variable τ = (1 + cos ϕ)/2, and our results indicate that the nuclear modification factor R A for TEEC exhibits a suppression in the range of 15% − 20% in the small τ region.Conversely, as τ increases, R A tends toward unity.This trend aligns with expectations, as larger τ values correspond to smaller dipole sizes being probed by TEEC, resulting in reduced nuclear modifications.
The demonstrated potential of measuring TEEC at the EIC underscores its importance in improving our understanding of gluon saturation and nuclear modifications.As the EIC becomes operational, we anticipate that the insights gained from TEEC measurements will play a pivotal role in refining our understanding of the fundamental aspects of strong interaction physics.

FIG. 1 .
FIG. 1. Illustration of TEEC for DIS in the lab frame (left panel).The incoming proton momentum PA and electron momentum ℓ define the z-axis.We align the transverse momentum of the outgoing electron p e T with the +y-direction to define the xy-plane (right panel).

FIG. 2 .
FIG.2.The TEEC plotted as a function of τ for e+p collisions at the future EIC.We choose the center-of-mass energy √ s = 140 GeV and the lepton rapidity ye = −2.The solid curves are from rcBK parameterizations while the dashed ones are based on the GBW model.The red, blue and green curves correspond to p e T = 2, 4 and 6 GeV, respectively.
FIG.3.Nuclear modification factor RA from Eq. (47) is plotted as a function of for p e T = 2 GeV (top panel), 4 GeV (middle panel), and 6 GeV (bottom panel).We choose √ s = 140 GeV and lepton rapidity ye = −2.The red bands are for the GBW model, and the blue bands are for the rcBK calculations.