Probing impact-parameter dependent nuclear parton densities from double-parton scatterings in heavy-ion collisions

We propose a new method to determine the spatially or impact-parameter dependent nuclear parton distribution functions (nPDFs) using the double-parton scattering (DPS) processes in highenergy heavy-ion (proton-nucleus and nucleus-nucleus) collisions. We derive a simple generic DPS formula in nuclear collisions by accommodating both the nuclear collision geometry and the spatiallydependent nuclear modification effect, under the assumption that the impact-parameter dependence of nPDFs is only related to the nuclear thickness function. While the geometric effect is widely adopted, the impact of the spatially-dependent nuclear modification on DPS cross sections is overlooked so far, which can however be significant when the initial nuclear modification is large. In turn, the DPS cross sections in heavy-ion collisions can provide useful information on the spatial dependence of nPDFs. They can be in general obtained in minimum-bias nuclear collisions, featuring the virtue of independence of the Glauber modeling.


I. INTRODUCTION
Multiple particle production at high-energy hadron colliders, such as at the LHC, is dominated by simultaneous multiple interactions between partons from the intial hadrons. Such multiple-parton interactions (MPIs) are indispensable in scrutinizing many event activities and hadron multiplicities at colliders. One particularly interesting case is that when more-than-one reactions in a collision are lying at hard scales, the perturbative QCD (pQCD) approach based on the factorization theorem [1] or conjecture applies. The studies of the so-called hard multiple-parton scattering processes can deepen our understanding of QCD and the possible multiple-parton correlations in a nucleon (see e.g. Refs. [2][3][4][5][6]). They provide new means to access the information of the nonperturbative structure of hadrons, which is complementary to the one obtained from nucleon one-body distributions. Due to the power counting of the cross sections, the leading multiple-parton scattering mechanism is the double-parton scattering (DPS), where only two partonic scattering subprocesses happen at the same time.
In the LHC era, we have witnessed the rapid theoretical developments , the vast phenomenological applications  and the impressive experimental measurements [69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84][85] of DPS and even triple-parton scattering (TPS) in the last decade, concentrating on protonproton (pp) collisions. Moreover, following the poineering work [86] by Strikman and Treleani, it is suggested that DPS cross sections in proton-nucleus (pA) and nucleusnucleus (AB) collisions will be largely enhanced thanks to the collision geometry. The DPS pA cross sections scale by three times of the number of nucleons A in a nucleus [87][88][89][90], while the single-parton scattering (SPS) cross sections in pA only scale by A (modulo other nuclear matter effects). The geometrical enhancement is more pronounced in nucleus-nucleus collisions. The DPS cross sections in AA collisions scale as A 3.3 /5 [91,92], while those for SPS are scaling as A 2 . However, one should bear in mind that such quantitive estimates are based on the assumption that the nuclear matter effects (thermal or non-thermal) are independent of the collision geometry, which is in fact not always justified. One counter example is the nuclear modification of the intial parton flux encoded in the nuclear parton distribution functions (nPDFs). Although the additional geometric effect from the nuclear modifications for DPS is anyway insignificant if the sizes of such modifications are small, the nuclear parton densities can deviate significantly from their free-nucleon counterparts at a scale of a few GeV (see e.g. Refs. [93][94][95]). Therefore, the existing DPS formula in heavy-ion collisions should be revised in order to incorporate the extra geometric/spatial effect from nPDFs 1 .
In addition, the understanding of the spatial (or impact-parameter) dependence of nPDFs is also essential to interpret the nuclear observables measured in different centrality classes, where one usually has to use the (optical or Monte Carlo) Glauber model to link the impact parameter and the centrality. Given the modest amount of available nuclear hard-reaction-process data available in the global fits, almost all considered nPDFs nowadays are only spatially averaged. They can only be directly used in minimum-bias nuclear collisions. The only one exception is that the authors of Ref. [105] determined the two impact-parameter dependent nPDFs 2 from the A dependencies of the two spatially averaged nPDFs. New means of extracting the spatial dependencies of the nPDFs from experimental data are therefore 1 To the best of our knowledge, none of the existing DPS phenomenological applications in heavy-ion collisions [89][90][91][92][96][97][98][99][100][101][102][103][104] considers such an effect. 2 There are also a few attempts to obtain the spatial form of the nPDFs based on phenomenological models (see e.g. Refs. [106,107]).
desired not only because the spatially dependent nPDFs in Ref. [105] are more-or-less obsolete but also due to the fact that a second independent validation is always valuable. The primary goal of the paper is to derive the generic expression for the DPS cross section in heavy-ion collisions by accommodating both the nuclear collision geometry and the spatially-dependent nuclear modification effect. Besides, we also suggest that the measurements of the DPS processes in minimum-bias collisions are able to constrain the impact-parameter dependent nPDFs. The remainder of the context is organized as follows. After introducing the nucleon density and the thickness function in section II, we derive a new generic formula for the DPS cross section in nucleus-nuclues collisions in section III. We explore the possibility of using the DPS cross sections in pA to determine the spatial-dependence of the nPDFs in section IV. A short summary is presented in section V. Finally, the appendix A discusses the transverse parton profile and the overlap function.

II. THE NUCLEON DENSITY AND THE THICKNESS FUNCTION
The nucleon number density in a nucleus is usually parameterized by a Woods-Saxon nucleon density function where r ≡ | − → r |, ρ 0,A corresponds to the nucleon density in the center of the nucleus, R A is the radius of the nucleus A, a A is the skin thickness, and w A characterizes deviations from a spherical shape. The concrete values of these parameters can be found in e.g. Ref. [108]. Such a function should work well for nuclei with A ≥ 4. An alternative simpler density function is the so-called hardsphere model, i.e.
where θ(x) is the Heaviside function. The normalization for ρ A is For convenience, we also define the nucleon probability A , which is normalized to unity. In the case of a single nucleon, we can write with δ n () is a n-dimensional Dirac delta function. For other A < 4 nuclei, one should utilize other nucleon density profiles. For instance, a profile for the deuterium was suggested in Ref. [105]. Let us consider a heavy-ion collision of A (the target) and B (the projectile), which is schematically shown in Fig. 1. Their transverse displacement is a 2-dimensional vector b . One considers two flux tubes located at the transverse displacement s and s − b with respect to the centers of the target A and the projectile B, respectively. In our notation, a three-dimensional vector − → x can be decomposed into a two-dimensional transverse part x and the longitudial part x z , i.e. − → x = ( x , x z ). The two-dimensional nucleon number density or the thickness function is where the nucleon probability density per unit transverse area isT The integrations of T A ( b ) andT A ( b ) over b in the whole two-dimensional area result into the nucleon number A and unity. One can define the thickness function T AB ( b ) and thickness probability functionT AB ( b ) for AB collisions as follows which are normalized to AB and unity after integrating over the transverse plane b . One can interpret the thickness probability functionT AB ( b ) as the effective overlap area for which a nucleon in A can meet with a nucleon in B, which is a pure geometrical factor.

III. DPS IN NUCLEUS-NUCLEUS COLLISIONS
The DPS cross section for a generic reaction AB → where the relevant geometry is shown in Fig. 1.σ f1 ik andσ f2 jl are the two partonic cross sections for ik → f 1 and jl → f 2 with the intial partons i, j, k, l being either (anti)quarks or gluons. δ f1f2 is the Kronecker delta function to take into account the symmetry of the final states for the generalized double parton distribution (GDPD) of the nucleus A 3 . They are whereΓ ij N/A andΓ i N/A are the isospin-averaged GDPD and the isospin-averaged generalized single parton distribution (GSPD) for a nucleon in A. The isospin average is written as for a nucleus A with Z protons and A − Z neutrons. We call Γ ij N A and Γ i N A the generalized double parton distribution in a nucleon and the generalized single parton distribution in a nucleon, respectively. We have used the bounded nucleon (bounded proton, bounded neutron) in a nucleus A as N A (p A , n A ), while a free-nucleon, a freeproton and a free-neutron are denoted as N, p and n respectively. The first term in the right-hand side of Eq. (9) represents that the two partons i and j belong to the same nucleon in the nuclues A with the impact parameter of the nucleon s 1 . The second term means that the two partons are from two distinct nucleons, where the prefactor A−1 A takes into account the difference between the number of nucleon pairs and the number of different nucleon pairs. Such a factor is essential to guarantee the correct normalization of Γ ij A , which was first noticed in Refs. [87,89]. If A is a nucleon (A = 1), the second term is zero because of the prefactor. Finally, we also have the decomposition for the GDPD Γ kl B of the nucleus B akin to Eq.(9) for the nucleus A.
If we use the factorized ansatz for the remaining nucleon GDPD Γ ij N A and the nucleon GSPDs where we have used t N A ( u ) as the transverse parton profile in the bounded nucleon N A and g i N A (x, s ) the impact-parameter s dependent nPDF for the parton i. A similar factorized ansatz is widely used in double parton scattering processes in (free) nucleon-nucleon collisions. It assumes the vanishing parton-parton correlations in DPS and yields the well-known "pocket formula" for the cross section in nucleon-nucleon N 1 N 2 collisions with the effective cross section as and the overlap function Usually, one assumes that the transverse parton profile t N ( u ) is independent of the type of the free-nucleon N , which is either a proton or a neutron. Then, we are left with one single effective cross section parameter σ eff,N1N2 = σ eff,pp , ∀ N i ∈ {p, n}. In addition, it is reasonable to assume that the transverse parton profile is not affected by the surrounding nucleons in a nucleus, i.e. t N A ( u ) = t N ( u ). In the following, we will use such two simplifications and retain only one t p and one F pp as the unique transverse parton profile and the overlap function. The further discussions about these two functions can be found in the appendix A. Motivated by the shadowing at small x and the Gribov-Glauber modeling [106] of the nPDFs, we can assume the nuclear matter effects encoded in nPDFs are only depending on the thickness function T A . We can introduce the general expression as where g i N (x) is the free-nucleon N PDF for the parton i and g i N A (x) is the spatially averaged nucleon N PDF for the parton i in A. G() can be an arbitary function with the normalization condition The simple form G T A ( s ) is the most frequently used one in the literature [109][110][111][112][113][114] and also in the HIJING event generator [115]. However, such a simple form conflicts with the A dependence of the nPDF global fit [105]. A study based on the A dependencies of the nPDFs reveals that a polynomial function G() with terms up to T A ( s ) can reproduce the nPDF A-dependence over the entire x range. In the following, for simplicity, we will use the abbreviations . Therefore, Eq.(15) can be reformulated as After applying the ansatz Eq.(11) and the relation Eq.(17) into Eq. (8), we arrive at the main result of the paper (2) nB,m3m4 where we have used Besides, the symbolsT are defined aŝ Because of the normalization relations, we havê The interpretation of the four terms in the brackets of Eq.(18) is straightforward. They represent three different DPS contributions from nucleus-nucleus interactions. The first term is from the two pairs of the colliding partons belong to the same pair of incident nucleons. The second and the third terms are originating from that two partons from a nucleon in a nucleus interact with the two partons from two different nucleons in anthoer nucleus. The last term is the contribution of the two pairs of partons belonging to two different nucleons from both nuclei. A few special situations worth being explored. When we take the identity of the impact-parameter dependent nPDF g i N A (x, s ) and the spatially-averaged nPDF g i N A (x) via G A,1 ( s ) = 1 and G A,2 ( s ) = 0, we can recover the well-known DPS formula in AB collisions (see e.g. Eqs. (1,2) in Ref. [102]), which however does not take into account the spatially-dependent initial nuclear modifications. Moreover, if we set g i N A (x) = g i N (x) (i.e. zero nuclear modification), we have g i The final expression is independent of G() as it must be. Finally, if we take the nucleus B as a proton, which amounts to setting B = 1, G B,1 ( s ) = 1, G B,2 ( s ) = 0, Eq.(18) is reduced to This gives rise a DPS formula in pA (or Ap) collisions. If A 1, we can impose a good approximation F pp ( v ) ≈ δ 2 ( v ). This can be understood because a nucleon in a heavy nucleus looks like a point in space.
With such a simplification, the transverse parton profile t p will only enter into σ eff,pp in Eq. (18) and Eq. (22). The goodness for the above approximation will be validated in the appendix A with a few concrete modelings of F pp ( v ).

IV. IMPACT-PARAMETER DEPENDENT NUCLEAR PDF FROM DPS
From Eq.(18) and Eq. (22), we know that the DPS cross sections in nuclear collisions depend on the function G() characterizing the impact-parameter dependence of nPDFs as introduced in Eq. (15). In turn, we can view DPS as a probe to determine the spatial dependence G(). The task of DPS cross section extraction from experimental data is however far non-trivial due to the presence of the contamination from the SPS contribution. An ideal case is to look for a process, in which the SPS contribution is suppressed. A few such examples are same-sign open charm [79], J/ψ+charm [79], Υ+charm [80] and J/ψ + Υ [41,74] production. In order to avoid the complications from the final-state nuclear effects, we only take the pA collisions as an example here. The above mentioned processes are dominated by gluon-gluon initial state at the LHC energies, which is blind with the isospin effect. For these gluon-induced processes, Eq. (22) can be further simplified. The nuclear modification factor is expressed as Aσ pp→f . One reasonable approximation we can take is that the nucleon number density follows the hard-sphere form Eq. (2). Then, the thickness function is The analytical expres- In such a case, the nuclear modification factor becomes It is easy to check that when a = 0, we are left with the first term proportional to R f1 pA R f2 pA . Let us take lead (Pb) beam with A = 208, R A = 6.624 fm and σ eff,pp = 15 mb as a special example. Such a beam is available at the LHC. Different numbers of the power a in G(x) ∝ x a predict quite different values of the nuclear modification factor R DPS pA→f1f2 , as reported in Fig. 2. The curves corresponding to five different values of R f pA = R f1 pA = R f2 pA are displayed. R DPS pA→f1f2 dramatically increases when a > 1.5, 2.0, 3.0, 1.0 for R f pA = 0.4, 0.6, 0.8, 1.2. As anticipated, the curve of R f pA = 1.0 (no nuclear modification) is independent of G(x) (or a).
As realistic examples, R f pA from the single-f inclusive processes with f being either the open charm or J/ψ mesons at the LHC pPb collisions were precisely measured to be close to 0.6 in the forward rapidity region (see e.g. Fig.1 in Ref. [93]). The a = 0, 1, 2, 3 predict R DPS pA→f1f2 = 1.27, 1.22, 1.07, 7.26 in the same kinematic regime. These numbers can be refined by using the Woods-Saxon density (cf. Eq.(1)) and with a concrete parton overlap function F pp ( v ) (cf. Eq.(14)). The numerical differences with respect to what we have shown should be minor though. From this example, we have clearly show that the nuclear modification factors of J/ψ plus open charm and same-sign charm production in proton-lead collisions will provide the precious inputs for determining the impact-parameter dependent nPDFs. Such measurements are independent of the centrality-based measurements, where the latter ones are crucially depending on the Glauber modeling (see e.g. Refs. [116,117] In this paper, for the first time, we have considered both the nuclear collision geometry and the impactparameter-dependent nuclear modification in the nPDFs for DPS processes in heavy-ion collions. A simple generic equation Eq. (18) has been derived for evaluating the DPS cross sections in nucleus-nucleus collisions, while its pA counterpart is given in Eq. (22). Both of the above effects are important in scrutinizing the DPS heavy-ion data. The latter is particular relevant when the nuclear modification encoded in the nPDFs is significant (e.g. the open/hidden charm and beauty production [93] at the LHC). In turn, we can also extract the spatial dependence of the nPDFs by measuring DPS cross sections in minimum-bias nuclear collisions. We take the gluoninduced charm and beauty production processes as an example. σ eff,pp can be determined from their pp data (e.g. Ref. [79]), and R f1 pA and R f2 pA are measured in their single inclusive processes. The measurements of the DPS nuclear modification factor R DPS pA→f1f2 can be readily used to pin down the spatial function G() entering into the impact-parameter-dependent nPDFs. Such an approach has the virtue of independence of the Glauber modeling. Several empirical functional forms of the transverse parton profile t p ( u ) in a nucleon were suggested in the literature [118][119][120][121]. They are collected in Tab. I. The "Dipole" profile is equivalent to the "Exponential" profile as long as we take r −1 0 = m g . Both of them are proportional to the modified Bessel function K 1 (). The analytical expressions of the mean three-dimensional radius squared − → r 2 and the mean two-dimensional radius squared u 2 can be found in Tab. II. Due to the spatial symmetry, we always have − → r 2 = 3 2 u 2 . We have also evaluated the analytic functions of the overlap function F pp ( v ) and of σ eff,pp for all these profiles in Tab. III.
As an illustration, in the following, we take the − → r 2 = (0.875 fm) 2 for all profiles, where 0.875 fm is the proton charge radius. Note that such values are not necessary agree with other tunings. For instance, Ref.
[119] took m 2 g = 1.1 GeV 2 from the analysis of the exclusive J/ψ photo-(or electro-) production. Such a value results into the value of − → r 2 1.5 times smaller than 0.875 fm. For the "Double Gaussian" profile, we adopt the values of β = 0.5, r0,1 r0,2 = 5 as suggested in Ref. [118]. In such a circumstance, we can predict the numerical values of σ eff,pp shown in the second column of Tab. IV. The setup results into pretty large values of σ eff,pp , ranging from 35 mb with "Double Gaussian" to 70 mb with "Top Hat". Alternatively, we can also fix the value of σ eff,pp to extract the parameters. The predicted − → r 2 are displayed in the third column of Tab. IV by using σ eff,pp = 15 mb. The − → r 2 values are generally 1.5 to 2.0 times smaller than 0.875 fm.
In Fig. 3, we have shown the comparisons between for the lead A = 208. Both the Woods-Saxon and hard-sphere ρ A have been used with the parameters R A = 6.624 fm, a = 0.549 fm and w = 0. We have tried the two transverse parton profiles "Hard Sphere" (HS) and "Dipole" (D). The approxima-      The predictions of σ eff,pp after imposing the mean three-dimensional radius squared 0.875 2 fm 2 (second column), and the values of the square root of the mean three-dimensional radius squared by fixing σ eff,pp = 15 mb (third column).