Three-dimensional imaging in nuclei

We perform the first simultaneous global QCD extraction of the transverse momentum dependent (TMD) parton distribution functions and the TMD fragmentation functions in nuclei. We have considered the world set of data from semi-inclusive electron-nucleus deep inelastic scattering and Drell-Yan di-lepton production. In total, this data set consists of 126 data points from HERMES, Fermilab, RHIC and LHC. Working at next-to-leading order and next-to-next-to-leading logarithmic accuracy, we achieve a $\chi^2/dof = 1.045$. In this analysis, we quantify the broadening of TMDs in nuclei comparing with those in free nucleons for the first time. We also make predictions for the ongoing JLab 12 GeV program and future EIC measurements.

Introduction. In recent years, quantum 3D imaging of the nucleon has become one of the hottest research topics in nuclear physics [1]. Such information is encoded in the transverse momentum dependent parton distribution functions (TMDPDFs) and significant progress has been made in extracting TMDPDFs for free nucleons from experimental data [2][3][4][5][6][7] as a central object in hadronic physics community. On the other hand, the corresponding quantum 3D imaging of a heavy nucleus is still at the primitive stage. Identifying the partonic structure of quarks and gluons in nuclei has remained as one of the most important challenges confronting the nuclear physics community since the pioneering EMC measurements in 1980s [8], and has been regarded as one of the major goals in future facilities of electron-ion colliders (EIC) [1,9,10]. Besides characterizing the non-trivial phenomena of nuclear modification of parton distribution inside bounded nucleons and the associated QCD dynamics, an accurate determination of such initial state nuclear effect is mandatory for providing precise benchmark information in searching for the signal of quarkgluon plasma created in heavy-ion collisions [11].
As demonstrated in both the generalized high-twist factorization formalism [32] and the dipole model [33,34], QCD multiple scattering in the nuclear medium is responsible for the difference between TMDPDFs in bound and free nucleons. This QCD multiple scattering leads to the so-called transverse momentum broadening effect, which manifests itself as the nuclear modification of the transverse momentum of the TMDPDFs in the framework of TMD factorization [35]. As such, while nTMD-PDFs represent the 3D partonic imaging of nuclei, they are also crucial for understanding the QCD dynamics of multiple scattering in nuclear medium. The accurate determination of nTMDPDFs is therefore one of the important objectives of the future EICs. Among the major goals of EICs, hadronization in medium is also of particular interest, which has been investigated experimentally such as in HERMES [36]. Such information is usually described by the nuclear modified fragmentation functions (nFFs) involved in the same collinear factorization as that in vacuum [37,38]. However, how the hadronization is influenced by medium in three-dimensional momentum space, i.e. nuclear modified transverse momentum dependent fragmentation functions (nTMDFFs), has never been explored.
The determination of nTMDPDFs and nTMDFFs (collectively called nTMDs) relies on the corresponding TMD QCD factorization theorem [35] for physical observables that involve two distinct scales, which are required to guarantee both the applicability of pQCD and the sensitivity to the parton's transverse motion. Two well-known observables are the transverse momentum distribution of semi-inclusive hadrons in leptonnucleus deep inelastic scattering (SIDIS) and of di-lepton in Drell-Yan (DY) processes in proton-nucleus (pA) collisions. These observables have been extensively mea-sured by HERMES [36], JLab [39,40], Fermilab [41,42], RHIC [43] and the LHC [44,45], and will be further measured by the future EIC [1,9,10] with unprecedented precision.
In this letter, we perform the first simultaneous QCD global analysis for nTMDPDFs and nTMDFFs using the world data from SIDIS and DY processes with nuclei. From our global analysis at next-to-leading order (NLO) and next-to-next-to-leading logarithmic (NNLL) accuracy, we demonstrate for the first time quantitatively the broadening of the transverse momentum of partons in bound nucleons.
TMD factorization formalism. To perform our global analysis, we select SIDIS and DY processes which have well-established TMD factorization formalism [35]. For ep SIDIS, e(l) + p(P ) → e(l ) + h(P h ) + X, the differential cross section at small hadron transverse momentum P h⊥ Q can be expressed as follows where, as in the standard TMD factorization, the result is written in the coordinate b-space that is conjugate to P h⊥ . We have dPS = dx dQ 2 dz d 2 P h⊥ with Q 2 = −(l − l) 2 , x and z the standard SIDIS kinematic variables, σ DIS 0 and H DIS are the Born cross section and the hard function for SIDIS. f q/p is the quark TMD-PDF inside a proton while D h/q denotes the TMDFF for q → h, with µ and ζ representing the renormalization and rapidity scales for TMDs. For the remainder of this paper, we always take µ = √ ζ 1 = √ ζ 2 = Q and replace the explicit dependence in the TMDs by the single scale Q. Within the so-called Collins-Soper-Sterman formalism [46], the evolved TMDs take the following form where C q←i andĈ i←q are the Wilson coefficient functions, ⊗ denotes the convolution, and f i/p (x, µ b * ) and D h/i (z, µ b * ) are the corresponding collinear PDFs and collinear FFs. Here µ b * = 2e −γ E /b * with γ E the Euler constant represents the natural scale for TMD evolution, while b * is the standard prescription. TMD evolution handles the evolution for both the longitudinal momentum fraction (x, z) and transverse component P h⊥ (or b in the coordinate space). The collinear functions in Eqs. (2) and (3) control the x (z) evolution via the usual DGLAP evolution equation. On the other hand, the perturbative (S pert ) and nonperturbative (S f,D NP ) Sudakov factors depend on b and Q, which control the corresponding perturbative (small b) and non-perturbative (large b) evolution on the parton's transverse momentum component and eventually resums logarithms in ln(Q 2 /P 2 h⊥ ) after the Fourier transform. While S pert is perturbatively calculable, the nonperturbative Sudakov factors have to be obtained by fitting experimental data and take the following form where g 2 (b) parameterizes the large-b behavior of socalled Collins-Soper evolution kernel and is both universal and independent of the species of external hadrons. We follow [47,48] to set g 2 (b) = g 2 ln(b/b * ). On the other hand, g q (g h ) represent the intrinsic transverse momentum of the TMDs at the initial scale Q 0 . In a simple Gaussian model, one typically has g q ∼ k 2 ⊥ /4 [49,50], likewise for g h . The parameters g 2 , g q , and g h in vacuum are all constrained in [47] with Q 0 = √ 2.4 GeV. For the DY process in pp collisions, p(P 1 ) + p(P 2 ) → γ * /Z(q) + X, the cross section in the TMD factorization region can be written as follows where dPS = dQ 2 dy d 2 q ⊥ with Q, y, q ⊥ the invariant mass, rapidity and transverse momentum of the vector boson, while c q (Q) denotes the quark coupling to the γ * /Z [4]. The term P takes into account the kinematic cuts on the transverse momentum, p ⊥ , and the rapidity, η, of the final state lepton pair [4,6,7]. In going from a proton to a nuclear target, we follow the same procedure [12] that is used for the nuclear collinear PDFs and FFs and make two assumptions. First, we assume that the TMD factorization takes exactly the same form as in Eq. (1), except that one replaces the TMDs by the nTMDs. Second, we assume that the perturbative physics for nTMDs and TMDs is the same. Using this assumption, the perturbative TMD evolution which is controlled by S pert would remain intact while the Wilson coefficient functions that enter into the OPE are also unchanged. Correspondingly, we would replace collinear functions in Eqs. (2) and (3) by their nuclear versions. In other words, these collinear functions would be modified at an initial scale Q = Q 0 and then evolved to the scale µ b * via the same DGLAP evolution. On the other hand, we modify the non-perturbative Sudakov S f,D NP to account for the effects from the nuclear medium. In principle, both g 2 (b) and the intrinsic components g q,h could be modified [51] due to the transverse momentum broadening through the parton multiple scattering in the nuclear target. We assume g 2 (b) to be the same as that for the proton target, and only replace g q,h by their corresponding nuclear version g A q,h . The g A q,h parameters would represent the parton's transverse momentum width inside a nucleus, which in general would depend on nuclear size (∝ A 1/3 ), momentum fraction x (or z) and the hard scale Q [52,53], denoted as g A q (x, Q) and g A h (z, Q). In the small-x or gluon saturation region, they would represent the typical size of saturation scale Q 2 s [34,54]. Global analysis. Considering the limited data available, in this paper we take the known parameters for collinear nPDFs f A i/p (x, Q 0 ) and nFFs D A h/i (z, Q 0 ), and perform the fit to extract g A q (x, Q) and g A h (z, Q). With more data in the future, one can simultaneously fit nuclear collinear functions and transverse modification encoded in g A q,h . Specifically, here we follow the EPPS16 [21] parameterization for collinear nPDFs with CT14nlo [55] for the proton PDFs, and we take LIKEn collinear nFFs in Ref. [56] for a nuclear target with the DSS14 parameterization [57] for the vacuum FFs. We perform our fit at NLO+NNLL accuracy. In the kinematic region probed by the current data, the x (or z) and Q-dependence is rather mild, which allows us to use two constant parameters a N and b N in the fit: where L = A 1/3 − 1. Such a modification is similar to the change of the saturation scale Q 2 s in the nucleus [33,34]. Thus, within our global analysis below, we have introduced the fit parameters a N and b N , which characterizes the nuclear broadening for the nTMDs.
For the data, we take SIDIS measurements from HER-MES collaboration and DY process from Fermilab, RHIC and the LHC. The HERMES collaboration [36], measured the hadron multiplicity ratio R A h = M A h /M D h , where the superscript A denotes the species of the nuclear target while D denotes a deuteron. On the other hand, , with the numerator given by the nuclear version of Eq. (1). The denominator is the inclusive DIS cross section, for which we use the APFEL library [58] at NLO with the collinear nPDFs. The CMS and ATLAS collaborations at the LHC directly measure the transverse momentum distribution for γ * /Z production, and we use the arTeMiDe library [4] to account for the phase space reduction in P. Finally, for Fermilab and RHIC, experimental measurements were performed for nuclear modification factor R AB = dσ A dPS / dσ B dPS , with A (B) the nuclear mass number of the heavy (lighter) nucleus.
In order to obtain the numerical values of the parameters a N and b N , we fit the experimental data using the Minuit package [59]. The normalization factors N of the LHC data are accounted for in the definition of the χ 2 according to the procedure of [21,57]. In Fig. 1, we plot the kinematic coverage of the world data and the kinematic coverage of future experimental data at JLab and the EIC. To select the HERMES data that is within the TMD region, we apply cuts P 2 h⊥ < 0.3 GeV 2 and z < 0.7. We note that in order to avoid correlations between the experimental data at HERMES, we choose to fit only one projection of the experimental data. Since the z dependent data provides the largest kinematic coverage, we chose to fit this experimental data. For the DY data, we enforce the standard kinematic cut q ⊥ /Q < 0.3. After performing these cuts, we are left with 126 points.
Results. The global analysis of these parameters re- sults in a χ 2 /d.o.f of 1.045 where the parameter values are given by a N = 0.0171 ± 0.003 GeV 2 and b N = 0.0144 ± 0.001 GeV 2 . The χ 2 distribution for the data sets used in our fit are provided in Table I. For the SIDIS data, we study only π production. The Baseline column represents the lighter nuclei used in the SIDIS multiplicity ratio and the DY nuclear modification factor. In our analysis, we have considered the uncertainty both from the result of our fit to a N and b N , and the uncertainty from the collinear distributions. In order to generate the fit uncertainties, we use the replica method in Refs. [60,61] with 200 replicas. However, at this point we note that we do not consider the collinear uncertainty when generating the replicas. In order to generate the uncertainty from the collinear nPDF and nFF, we use the prescription provided in Ref. [21]. The collinear and fit uncertainties are both displayed at 68%.
In Fig. 2, we plot the result of our fit against the experimental data. In the top row of this figure, we plot the comparison against: the multiplicity ratio measurement at HERMES [36] as a function of z (left two columns) and P h⊥ (third column from the left), and the DY q ⊥ distribution from the LHC (right column). We note that the P h⊥ dependent data in the third column is a prediction for those data points. Furthermore, for the LHC data [44,45], we have provided the N i for each of the data sets. In the left three columns of the second row, we plot the comparison against the R AB ratio for the E866 [42] and E772 [41] experiments. Finally, in the right column of this row, we plot the R AB at RHIC [43]. In each subplot, we have provided the uncertainty from our fit as a dark band, and the uncertainty from the collinear distributions as a light band.
In the top row of Fig. 3, we plot the ratio of the u-quark TMDPDF of a bound proton in a gold nucleus and that in a free proton as a function of x and k ⊥ . The curve along the plane for k ⊥ = 1 GeV demonstrates the shadowing, anti-shadowing, and the EMC effect which originate from the collinear distribution. The curves which lie in planes of constant x increase with increasing k ⊥ , which indicates the transverse momentum broadening effect, with a suppression at low k ⊥ and an enhancement at high k ⊥ , as expected from our model in Eq. (7). In the bottom row of this figure, we plot the ratio of the nTMDFF for u → π + in a Xe nucleus and that in vacuum as a function of z and p ⊥ . Analogous to the nTMDPDFs, we see that as p ⊥ grows, this ratio becomes larger, indicating that hadrons originating from fragmentation in the pres-ence of a nuclear medium will tend to have a broader distribution of transverse momentum relative to vacuum TMDFFs.
In Fig. 4, we plot our prediction for future JLab and EIC multiplicity ratio measurements as a function of P h⊥ for π + at z = 0.4. For the EIC, we plot the prediction at x = 0.05 and Q 2 = 4 GeV 2 (black) and Q 2 = 100 GeV 2 (red). For JLab, we plot the prediction at x = 0.4 and Q 2 = 2.5 GeV 2 (green). We expect the future measurements will provide a stringent constrain of nTMDs and test the QCD evolution shown in Fig. 4.

Summary.
We perform the first QCD global analysis of nuclear TMDs. For the processes with nuclei, assuming that TMD factorization and perturbative TMD evolution both take the same form as those in the vacuum, at the accuracy of NLO+NNLL we find that we can describe the global set of experimental data using a simple model which accounts for the nonperturbative TMD evolution. We demonstrate quantitatively that both the TMDPDFs and TMDFFs in the presence of the nuclear medium have a broader distribution of transverse momentum. We expect that the framework we have developed will have a large impact on the interpretation of future experimental data at JLab, RHIC, LHC, and the future EICs, allowing us to perform quantum 3D imaging of the nucleus.