Transverse momentum spectrum of dilepton pair in the unpolarized $\pi^-N$ Drell-Yan process within TMD factorization

We study the transverse momentum spectrum of dilepton produced in the unpolarized $\pi^- N$ Drell-Yan process, using transverse momentum dependent factorization up to next-to-logarithmic order of QCD. We extract the nonperturbative Sudakov form factor for the pion in the evolution formalism of the unpolarized TMD distribution function, by fitting the experimental data collected by the E615 Collaboration at Fermilab. With the extracted Sudakov factor, we calculate the normalized differential cross section with respect to transverse momentum of the dimuon and compare it with the recent measurement by the COMPASS Collaboration.


Introduction
After the first observation of the µ + µ − lepton pairs produced in p N collisions [1], the process was interpreted as a quark and an antiquark from each initial hadron annihilating into a virtual photon, which in turn decays into a lepton pair [2]. This explanation makes the process an ideal tool in order to explore the internal structure of both the beam and target hadrons. Since then, a wide range of studies on this (Drell-Yan) process have been carried out. In particular, the πN Drell-Yan process has the unique capability to pin down the partonic structure of the pion, which is an unstable particle, it therefore cannot serve as a target in deep inelastic scattering processes. Several pion induced experiments have been carried out, such as NA10 at CERN [3][4][5][6], E615 [7], E444 [8] and E537 [9] at Fermilab. These measurements have provided plenty of data, which have been used to constrain considerably the pion distribution function. A new pion-induced Drell-Yan program was also proposed to be executed at the COMPASS facility [10]. In particular, the COMPASS collaboration measured for the first time the transverse-spin-dependent azimuthal asymmetries in the π − N Drell-Yan process, using a high-intensity π beam of 190 GeV colliding on a transversely polarized NH 3 proton target [11].
One of the main goals of the COMPASS Drell-Yan program is to measure the Sivers asymmetry or the Sivers function with high precision. The Sivers function is a transverse momentum dependent (TMD) parton distribution function (PDF) describing the asymmetric density of unpolarized quarks inside a transversely polarized nucleon. Remarkably, QCD predicts that the sign of the Sivers function changes in SIDIS with respect to the Drell-Yan process [12][13][14]. The verification of this sign change [15][16][17][18][19][20] is one of the most fundamental tests of our understanding of the QCD dynamics and the factorization scheme, and it is also the main pursue of the existing and future Drell-Yan facilities [10,11,[21][22][23][24]. The advantage of the π N Drell-Yan measurement at COMPASS is that almost the same setup [11,25] is used in SIDIS and Drell-Yan, which may reduced the uncertainty in the extraction of the Sivers function.
In order to quantitatively understand the Sivers asymmetry in the πN Drell-Yan process at COMPASS, a first step is to know with high accuracy the differential cross-section of the same process for an unpolarized nucleon, since the spin-averaged cross-section is what appears in the denominator of the asymmetry definition. This motivates us to study the unpolarized πN Drell-Yan process in a comprehensive way. Since most of the data accumulated by the COMPASS collaboration are at low transverse momentum (of the dilepton q ⊥ ), a natural choice for the analysis framework is TMD factorization, which is valid in the region where q ⊥ is much smaller than the hard scale Q. The TMD factorization has been widely applied to various high energy processes, such as the SIDIS [26][27][28][29][30][31][32], e + e − annihilation [27,33,34], Drell-Yan [27,35,36] and W/Z production in hadron collision [27,37,38]. From the point of view of TMD factorization [26,27,37,39], physical observables can be written as convolutions of a factor related to hard scattering and well-defined TMD distribution functions or fragmentation functions. We will calculate the cross-section of the process up to the next-to-leading logarithm (NLL) and also take into account the evolution effects of the TMD distributions.
In contrast to what is done in the collinear PDFs case, the evolution of TMDs is usually performed in b space [27,37], which is conjugate to the transverse momentum k ⊥ . Moreover, the TMD PDFs depend on two energy scales. One is related to the corresponding collinear PDFs and the other is related to the definition of the TMD PDFs. After solving the evolution equations, the TMDs at fixed energy scale can be expressed as a convolution of their collinear counterparts and perturbatively calculable coefficients, and the evolution from one energy scale to another energy scale is included in the exponential factor of the so-called Sudakov-like form factors [27,37,40]. The Sudakov form factor can be separated into a perturbatively calculable part and a nonperturbative part, the later one can be fitted to experimental data. In this work, we perform the first extraction of the nonperturbative Sudakov form factor for the unpolarized TMD PDF f 1 of pion meson in the π − P Drell-Yan process, which sheds light on the partonic structure of the pi meson.
The rest of the paper is organized as follows. In Sec. 2 we set up the necessary theoretical framework of the TMD factorization and evolution effects in the unpolarized π − N Drell-Yan process. We also propose a prescription for the nonperturbative Sudakov form factor in the case of the spin-independent pion distributions. In Sec. 3 we calculate numerically the transverse momentum dependent differential cross section, based on the framework given in Sec. 2, and we fit the nonperturbative Sudakov form factor using the q ⊥ -dependent experimental data obtained by the E615 collaboration at Fermilab [7]. In Sec. 4 we present the prediction for the q ⊥ spectrum of the dilepton pair produced in the π − N collision at the COMPASS kinematics, applying the extracted Sudakov form factor obtained in Sec. 3.
In Sec. 5 we summarize the results of the paper and give some conclusions.

Framework of unpolarized Drell-Yan process in TMD factorization
In this section we will set up the necessary framework of TMD factorization, following the procedure presented in Ref. [27], in order to obtain the differential cross section for unpolarized Drell-Yan process considering TMD evolution effects.

Kinematics and general formula for Drell-Yan Process
The process under study is the unpolarized π − P Drell-Yan process where P π , P p and q denote the momenta of the π − meson, the proton and the virtual photon, respectively. In contrast to the SIDIS process, in the Drell-Yan process q is a time-like vector, namely, Q 2 = q 2 > 0, which is the invariant mass square of the lepton pair. In order to express the differential cross section of the process, we define the following kinematical variables where s is the center-of-mass energy squared; x π and x p are the light-front momentum fraction carried by the annihilating quarks in both the π − and the proton, respectively; q L is the longitudinal momentum of the virtual photon in the c.m. frame of the incident hadrons; x F is the Feynman x variable, which corresponds to the longitudinal momentum fraction carried by the lepton pair; and y is the rapidity of the lepton pair. Thus, x π and x p can be expressed as functions of x F , τ and of y, τ Since it is much more convenient to solve the TMD evolution equations in coordinate space (b-space) than in the transverse momentum space q ⊥ , where b is conjugate to q ⊥ via Fourier Transformation, the differential cross section is usually expressed as a b−dependent function, formulated in TMD factorization. Later on it can be transformed back to transverse momentum space to represent the experimental observables. Thus, the differential cross section for the unpolarized π − -proton Drell-Yan process described in (2.1) has the form [37] where σ 0 = 4πα 2 em 3N C sQ 2 is the cross section at tree level. Hereafter, we use the tilde to denote b−space terms. The structure function W (Q; b) contains all-order resummation results and dominates at low q ⊥ ≪ Q value, while the term Y U U provides the necessary correction at moderate q ⊥ ∼ Q values. In this work we will neglect the Y -term, which means that we will only consider the first term on the r.h.s of Eq. (2.4).
In general, TMD factorization [27] aims at separating well defined TMD distributions, such that they can be used in different processes, and including the scheme/process dependence in the corresponding hard factors. Thus, W (Q; b) can be expressed as [41] wheref sub q/H is the subtracted distribution function in b−space with the superscript "sub" representing the fact that the soft factor is subtracted, H U U (Q; µ) is the factor associated with hard scattering, µ is renormalization scale in the case of collinear PDFs, and ζ F denotes an energy scale related to the cutoff of the TMD distributions. We note that the distribution functionf sub 1 q/H is universal, since it can be used in different processes. However, the way to subtract the soft factor in the distribution function depends on the scheme to regulate the light-cone singularity in the TMD definition [26], whereas the hard factor H U U (Q; µ) is also scheme-dependent. In the literature, two different schemes are usually used: the Collins-11 scheme [27] and the Ji-Ma-Yuan scheme [28,39]. If we perform a Fourier Transformation onf sub 1 q/H (x, b; µ, ζ F ), we obtain the distribution function in the transverse momentum space which contains the information about the probability of finding a quark with collinear momentum fraction x and transverse momentum k ⊥ .

TMD evolution and Sudakov form factor
The TMD evolution for the ζ F dependence of TMD PDFs is encoded in a Collins-Soper (CS) [27] equation through while the µ dependence is derived from the renormalization group equation as whereK is the CS evolution kernel, and γ K and γ F the anomalous dimensions. The solutions of these evolution equations are studied in detail in Ref. [27]. Here, we will only discuss the final result. The overall solution structure forf 1 (x, b; µ, ζ F ) is the same as that for the Sudakov form factor. Namely, the energy evolution of TMDs from an initial energy µ to another energy Q is encoded in the Sudakov-like form factor S through the exponential Here, F is the hard factor, which depends on the scheme that we choose. The b-dependence can determine the transverse momentum dependence of experimental observable through Fourier Transformation, which makes understanding the b-dependence become quite important. In the small b region, the dependence is perturbatively calculable, while at large b the dependence turns nonperturbative and should be obtained from experiment data. To combine the perturbative information at small b with the nonperturbative part at large b, a matching procedure must be introduced with a parameter b max serving as boundary between the perturbative and nonperturbative regions. A b-dependent function b * is defined, which has the property b * ≈ b at low values of b and b * ≈ b max at large b values. A typical value of b max is chosen around 1 GeV −1 , such that b * is always at the perturbative region. There are several different b * prescriptions in the literature [42,43]. In this work we adopt the original prescription introduced in Ref. [27], In the small b region 1/Q ≪ b ≪ 1/Λ, the defined TMDs at fixed energy µ can be expressed as a convolution of perturbatively calculable hard coefficients and the corresponding collinear PDFs [26,44] where ⊗ stands for the convolution in the momentum fraction x is the corresponding collinear PDF of the i flavor in hadron H at the energy scale µ, which can be a dynamic scale related to b * by µ b = c 0 /b * , with c 0 = 2e −γ E and γ E ≈ 0.577, the Euler Constant [26]. In addition, the sum Σi runs over all parton flavors. Independently on the type of initial hadrons, the perturbative hard coefficients have been calculated for the parton-target case [29,45] and the results have been presented in Ref. [44] (see also Appendix A of Ref. [29]) as , and α s is the strong coupling constant at the energy scale µ b . The b * prescription prevents α s (µ b ) from hitting the so-called Landau pole at large b regime. We note that the C-coefficients do not depend on the TMD scheme and on the initial hadrons. Namely, the C-coefficients in Eq. (2.12) are universal for both different schemes and initial hadrons.
The Sudakov-like form factor in Eq. (2.9) can be separated into a perturbatively calculable part and a nonperturbative part (2.14) The perturbative part of S has the form The coefficients A and B in Eq.(2.15) can be expanded as a α s /π series: In this work, we take A (n) to A (2) and B (n) to B (1) up to the accuracy of next-to-leadinglogarithmic (NLL) order [29,32,37,[46][47][48]: For the nonperturbative form factor S NP in pp Drell-Yan process, a general parameterization associated with the TMD PDF has been proposed in Ref. [49] and it has the form In Ref. [49] the parameters g 1 , g 2 , g 3 are fitted from the nucleon-nucleon Drell-Yan process data at the initial scale Q 2 0 = 2.4 GeV 2 with b max = 1.5 GeV −1 , x 0 = 0.01 and λ = 0.2. With the extracted parameters, the fitted results are g 1 = 0.212, g 2 = 0.84, g 3 = 0.
Since the nonperturbative form factor S NP for quarks and antiquarks satisfies the following relation [41] 22) and the quark and antiquark contributions to S NP are the same: the S NP associated with the TMD distribution function of one of the initial protons can be expressed as With all the ingredients presented above, we can rewrite the b−space unpolarized distribution function f 1 (e.g., for an up quark inside the proton), which depends on x, b, Q, asf where g(x, µ b ) stands for the distribution function of the gluon and µ b turns to If we perform a Fourier Transformation onf sub 1q/p (x, b; Q), we can obtain the distribution function in transverse momentum space where J 0 is the Bessel function of the first kind. The nonperturbative Sudakov form factor for the pion TMD distributions has never been obtained. Here we assume that it has the same structure as for the proton TMD distributions, with parameters g π 1 and g π 2 S f q/π 1 Similarly, we can express the evolved result for the pion TMD distributions as and also in transverse momentum space as

Expression for the differential cross section
We note that the hard factors F(α s (Q)), which appears in Eq.(2.25) and Eq.(2.29), and H U U (Q; µ) in Eq. (2.5), are scheme dependent. In the two different schemes, the Collins-11 scheme [27] and the Ji-Ma-Yuan scheme [28,39], the hard factors F(α s (Q)) and H U U (Q; µ) have different forms: where ρ is another parameter to regulate the light-cone singularity of TMD distributions. However, if one absorbs both of the hard factors H and F into the definition of the Ccoefficients, the C-coefficients turn out to be ρ-independent (scheme-independent) and can be expresses as [50] C Substituting the above new C-coefficients into Eq. (2.25) and Eq. (2.29), we can obtain the evolved TMD distributions used in the calculation of the differential cross section for the unpolarized Drell-Yan process. Finally, the structure function in b-space can be written as After performing the Fourier Transformation, we can get the differential cross section, written in Eq. (2.4), as

Extracting the nonperturbative Sudakov form factor for the pion TMD distribution
Using the framework set up above, in this section we fit the theoretical cross-section of π − N Drell-Yan at the kinematics of E615 measurement, to the corresponding experimental data [7], in order to extract the parameters of the Sudakov form factor for the pion TMD distributions.

The differential cross section at E615 experiment
The E615 experiments at FermiLab produced a π − beam with P π = 252 GeV, which collided on a tungsten target (74 protons and 110 neutrons). The differential cross-section was given as 2-fold data depending on x F and q T [51], while the formula in Eq. (2.38) denotes the 4-fold differential cross section as a function of Q 2 , y, q T . Thus, we apply the kinematics at E615 to perform the integration over Q 2 from 4.05 2 GeV 2 to 8.55 2 GeV 2 , and we transform the differential with respect to x F to arrive at We adopt the NLO set of the CT10 parametrization [52] (central PDF sets) for the unpolarized distribution function f 1 (x) of the proton. For the unpolarized distribution function of the pion meson, we use the NLO SMRS parametrization [53]. In our calculation we ignore any nuclear effects from the tungsten target. The values of the strong coupling α s are consistently obtained at 2-loop order as with fixed n f = 5 and Λ QCD = 0.225 GeV. We note that the running coupling in Eq. (3.3) satisfies α s (M 2 Z ) = 0.118. The experimental data that we use in the fit are the q ⊥ -dependent differential cross sections of the unpolarized π − N Drell-Yan from the E615 experiment, as listed in Table. 36 of Ref. [51]. The data span the region 0 < x F < 1 with ten x F bins in a interval of 0.1. However, we find that a good fit can be obtained if we exclude the data in the bins 0.8 < x F < 0.9 and 0.9 < x F < 1. Therefore, in the fit we disregard the data in the region x F > 0.8. To estimate the result in each x F bin, we calculate the q ⊥ -dependent cross section at the mean value of the x F bin.

The fitting procedure
The data fitting is performed by the package MINUIT [54,55], through a least-squares fit in which the chisquare is defined as In the definition above, i stands for the i−th x F bin out of the total of M bins, each one with N i data points, α is the vector of the free parameters: g π 1 and g π 2 , q ⊥ij denotes the j−th value of q ⊥ measured in the i−th x F bin, theo(q ⊥ij , α) is the theoretical estimate at q ⊥ij and α, and exp ij is the data measured at q ⊥ij with error ∆err ij . The total number of data in our fit is N = 8 i N i = 96. We note that, apart from the statistical error in Table. 36 of Ref. [51], there is also an overall systematic error of 16%. Thus, in our fit the error is the total error which has the form ∆err tot = err 2 stat + err 2 sys . Since the TMD formalism is valid in the region q ⊥ ≪ Q, we do a simple data selection by removing the data in the region q ⊥ > 3 GeV. We perform the fit by minimizing the chisquare in Eq. (3.4), and we obtain the following values for the two parameters: g π 1 = 0.082 ± 0.022, g π 2 = 0.394 ± 0.103. with χ 2 /d.o.f = 1.64. From the fitted result, we find that the values of the parameters g 1 and g 2 for the pion meson are smaller than those for the proton. Similar to the case of the proton, for the pion meson g π 2 is several times larger than g π 1 . In Fig. 1, we plot the q ⊥ -dependent differential cross section at different x F bins, using the fitted values for g π 1 and g π 2 in Eq. (3.5). The experimental data (full square) at E615 with error bars are also depicted for comparison. As Fig. 1 demonstrates, a good fit is obtained at the region x F < 0.8. We also check the theoretical result in the region x F > 0.8. It turns out that, using the two parameters of the Sudakov form factor for the pion meson, the theoretical calculation underestimates the experimental data in that region, particularly when x F > 0.9. The reason or this is that at large x F the higher twist effects dominates, and therefore the leading-twist TMD factorization breaks down.

Results of TMDs for the pion at different scales
Using Eq. (2.29) and the extracted parameters g π 1 and g π 2 , we quantitatively study the scale dependence of the pion TMD distributions. In the left panel of Fig. 2 we plot the subtracted valence distribution function in b−spacef sub 1 u/π + (x, b; Q) (multiplied by b), for fixed x π = 0.1, at three different energy scales: Q 2 = 2.4 GeV 2 (dotted line), 10 GeV 2 (solid line), 1000 GeV 2 (dashed line). In the right panel of Fig. 2 we plot similar curves, but for the distribution f 1 u/π + in transverse momentum space. In this calculation, we have applied the ρ-independent C−coefficients presented in Eqs. (2.35) and (2.36). From the b-dependent plots, we see that at the highest energy scale Q 2 = 1000 GeV 2 , the peak of the curve is in the low b region where b < b max , since in this case the perturbative part of the Sudakov form factor dominates. However, at relatively lower energy scales, e.g., Q 2 = 10 GeV 2 and Q 2 = 2.4 GeV 2 , the peak of the b−dependent distribution function moves towards the higher b region, indicating that the nonperturbative part of the TMD evolution becomes more important at this lower energy scales. For the distribution in transverse momentum k ⊥ space, at the high energy scale the distribution has a tail falling off slowly at large k ⊥ , while at lower energy scales the distribution function falls off rapidly with increasing k T . It is interesting to point out that the shapes of the pion TMD distribution at different scales are similar to those of the proton [56].  4 Prediction of the q ⊥ spectrum of the lepton pair at COMPASS As a cross check, we calculate the q ⊥ spectrum of the lepton pair produced in the unpolarized pion-proton Drell-Yan process at the kinematics of COMPASS, using the TMD framework and applying the nonperturbative Sudakov form factor extracted in the previous section. The COMPASS Collaboration recently reported the measurements on the transverse-spindependent azimuthal asymmetries with three different modulations in the π − N Drell-Yan process [11], in which a π − beam with P π = 190 GeV collides on a NH 3 target (10 protons and 7 neutrons). These asymmetries provide a great opportunity to access the proton TMD distributions, such as the Sivers function, the transversity and prezelozity distributions, as well as the pion Boer-Mulders function. Besides the measurement on the spin asymmetries, the COMPASS collaboration also presented the distribution of dimuon pairs as a function of the transverse momentum q ⊥ of the dimuon, up to q ⊥ = 5.4 GeV, which corresponds to the result in the case of an unpolarized target, since the azimuthal angles are integrated out. The study of the unpolarized cross-section of πN Drell-Yan is equally important, since it appears in the denominator of various asymmetries. In particular, the study on the transverse momentum spectrum may also shed light on the unpolarized TMD distribution of the pion.
We use the following kinematics at COMPASS to calculate the q ⊥ distribution of the dimuon in piN Drell-Yan process 0.05 < x N < 0.4, 0.05 < x π < 0.9, 4.3 GeV < Q < 8.5 GeV.  Figure 3. The transverse spectrum of lepton pair production in the unpolarized pion-nucleon Drell-Yan process, with an NH 3 target at COMPASS. The dashed line is our theoretical calculation using the extracted Sudakov form factor for the pion TMD PDF. The solid line shows the experimental measurement at COMPASS.
The differential cross section at COMPASS can be calculated via Eq. (3.2) by replacing the integration limit of Q. As shown in Ref. [11], the events of the dimuon production were measured for different q ⊥ bins with an interval of 0.20 GeV. It can be also written as N = σ × L, where σ is the total cross section and L is the integrated luminosity of incident hadrons. To cancel the uncertainty of the integrated luminosity, we resort to the normalized differential cross section defined as [57] 1 σ Here, N i refers to the measured events at the i−th q ⊥ bin, and ∆q ⊥ represents the width of q ⊥ bin. For the cross section dσ dq ⊥ as the function of q ⊥ , we perform the integration over the entire x F region at COMPASS and estimate the value for different q ⊥ bins. The total cross section σ is calculated by integrating over q ⊥ and x F . To make the TMD factorization valid, we choose the upper limit of q ⊥ in our calculation to be 2.6 GeV.
In Fig. 3, we plot the normalized transverse momentum spectrum of the dimuon pair produced in the pion-nucleon Drell-Yan process at COMPASS, at different q ⊥ bins, with an interval of 0.2 GeV. The dashed line depicts our theoretical estimate which implements the extracted nonperturbative Sudakov form factor for the pion TMD distribution. The solid line corresponds to the recent experimental measurement by the COMPASS collaboration [11]. Comparing the two curves, we find that our theoretical estimate on the q ⊥ distribution of the dimuon agrees with the COMPASS data fairly well in the small q ⊥ region, in both the shape and the size. In the region q ⊥ < 2 GeV, the difference between calculation and data is less then 10%. This validates our extraction of the nonperturbative Sudakov form factor for the pion distribution f 1π , within the TMD factorization. At relatively larger q ⊥ values, our calculation cannot describe the data, indicating that the perturbative correction from the Y U U term plays an important role in the region q ⊥ ∼ Q. Further study is needed in order to provide a complete picture of the q ⊥ distribution of dilepton pairs from πN Drell-Yan in the whole q ⊥ range.

Conclusion
In this work, we applied the framework of the TMD factorization to study the transverse momentum distribution of the dilepton in the unpolarized π − N Drell-Yan process. We took into account the TMD evolution of the distributions and considered the nonperturbative Sudakov form factor for the unpolarized TMD distribution function of the pion meson. We parameterized the form factor in a form analogous to that of the proton TMD distribution, and we performed a fit using the experimental data for the differential cross section as a function of q ⊥ at E615. Adopting the extracted values of the parameters g π 1 and g π 2 , we presented the evolved results of the unpolarized pion TMD distributions f 1π at different energy scales. We found a similarity between the TMD evolution of the pion distribution and that of the the proton distribution. Finally, we calculated the transverse momentum spectrum of dimuon pairs produced in unpolarized pion-nucleon Drell-Yan processes at the COMPASS kinematics, and compared it with the COMPASS data to verify our extraction. As a result, the extracted nonperturbative Sudakov form factor is compatible with the COMPASS measurement at small q ⊥ region with q ⊥ ≪ Q, indicating that our approach can be used as a first step to study the Drell-Yan process at COMPASS. Our study may also provide a better understanding on the pion TMD distribution as well as its role in Drell-Yan process. Furthermore, the framework applied in this work can also be extended to the study of the azimuthal asymmetries in the πN Drell-Yan process, which we reserve for a future study.