# π - N Drell-Yan Process in TMD Factorization

Authors: Xiaoyu Wang; Zhun Lu
Journal: Advances in High Energy Physics
 (2019)

---

## Abstract

This article presents the review of the current understanding on the pion-nucleon Drell-Yan process from the point of view of the TMD factorization. Using the evolution formalism for the unpolarized and polarized TMD distributions developed recently, we provide the theoretical expression of the relevant physical observables, namely, the unpolarized cross section, the Sivers asymmetry, and thecos⁡2ϕ asymmetry contributed by the double Boer-Mulders effects. The corresponding phenomenology, particularly at the kinematical configuration of the COMPASS πN Drell-Yan facility, is displayed numerically.


---

## Body

## 1. Introduction

After the first observation of theμ+μ- lepton pairs produced in pN collisions [1], the process was interpreted that a quark and an antiquark from each initial hadron annihilate into a virtual photon, which in turn decays into a lepton pair [2]. This explanation makes the process an ideal tool 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 and therefore cannot serve as a target in deep inelastic scattering processes. Several pion-induced experiments have been carried out, such as the NA10 experiment at CERN [3–6], the E615 [7], E444 [8], and E537 [9] experiments at Fermilab three decades ago. These experimental measurements have provided plenty of data, which have been used to considerably constrain the distribution function of the pion meson. Recently, a new pion-induced Drell-Yan program with polarized target was also proposed [10] at the COMPASS of CERN, and the first data using a high-intensity π beam of 190 GeV colliding on a NH3 target has already come out [11].Bulk of the events in the Drell-Yan reaction are from the region where the transverse momentum of the dileptonq⊥ is much smaller than the mass Q of the virtual vector boson; thus the intrinsic transverse momenta of initial partons become relevant. It is also the most interesting regime where a lot of intriguing physics arises. Moreover, in the small q⊥ region (q⊥~ΛQCD), the fixed-order calculations of the cross sections in the collinear picture fail, leading to large double logarithms of the type αln2(q⊥2/Q2). It is necessary to resum such logarithmic contributions to all orders in the strong coupling αs to obtain a reliable result. The standard approach for such resummation is the Collins-Soper-Sterman (CSS) formalism [12], originated from previous work on the Drell-Yan process and the e+e- annihilation three decades ago. In recent years the CSS formalism has been successfully applied to develop a factorization theorem [13–15] in which the gauge-invariant [16–19] transverse momentum dependent (TMD) parton distribution functions or fragmentation functions (collectively called TMDs) [20, 21] play a central role. From the point of view of TMD factorization [12, 13, 15, 22], physical observables can be written as convolutions of a factor related to hard scattering and well-defined TMDs. 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 in the perturbative region, 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 [12, 15, 23, 24]. The TMD factorization has been widely applied to various high energy processes, such as the semi-inclusive deep inelastic scattering (SIDIS) [14, 15, 22, 23, 25, 26], e+e- annihilation [15, 27, 28], Drell-Yan [15, 29], and W/Z production in hadron collision [12, 15, 30]. The TMD factorization can be also extended to the moderate q⊥ region where an equivalence [31, 32] between the TMD factorization and the twist-3 collinear factorization is found.One of the most important observables in the polarized Drell-Yan process is the Sivers asymmetry. It is contributed by the so-called Sivers function [33], a time-reversal-odd (T-odd) distribution describing the asymmetric distribution of unpolarized quarks inside a transversely polarized nucleon through the correlation between the quark transverse momentum and the nucleon transverse spin. Remarkably, QCD predicts that the sign of the Sivers function changes in SIDIS with respect to the Drell-Yan process [16, 34, 35]. The verification of this sign change [36–41] is one of the most fundamental tests of our understanding of the QCD dynamics and the factorization schemes, and it is also the main pursue of the existing and future Drell-Yan facilities [10, 11, 42–45]. The advantage of the πN Drell-Yan measurement at COMPASS is that almost the same setup [11, 46] is used in SIDIS and Drell-Yan processes, which may reduce the uncertainty in the extraction of the Sivers function. In particular, the COMPASS Collaboration measured for the first time the transverse-spin-dependent azimuthal asymmetries [11] in the π-N Drell-Yan process.Another important observable in the Drell-Yan process is thecos⁡2ϕ angular asymmetry, where ϕ corresponds to the azimuthal angle of the dilepton. The fixed-target measurements from the NA10 and E615 collaborations showed that the unpolarized cross section possesses large cos⁡2ϕ asymmetry, which violates the Lam-Tung relation [47]. Similar violation has also been observed in the pp colliders at Tevatron [48] and LHC [49]. It has been explained from the viewpoints of higher-twist effect [50–53], the noncoplanarity effect [30, 54], and the QCD radiative effects at higher order [55, 56]. Another promising origin [57] for the violation of the Lam-Tung relation at low transverse momentum is the convolution of the two Boer-Mulders functions [58] from each hadron. The Boer-Mulders function is also a TMD distribution. As the chiral-odd partner of the Sivers function, it describes the transverse-polarization asymmetry of quarks inside an unpolarized hadron [57, 58], thereby allowing the probe of the transverse spin physics from unpolarized reaction.This article aims at a review on the current status of our understanding on the Drell-Yan dilepton production at low transverse momentum, especially from theπN collision, based on the recent development of the TMD factorization. We will mainly focus on the phenomenology of the Sivers asymmetry as well as the cos⁡2ϕ asymmetry from the double Boer-Mulders effect. In order to quantitatively understand various spin/azimuthal asymmetries in the πN Drell-Yan process, a particularly important step is to know in high accuracy the spin-averaged differential cross section of the same process with azimuthal angles integrated out, since it always appears in the denominator of the asymmetries’ definition. Thus, the spin-averaged cross section will be also discussed in great details.The remained content of the article is organised as follows. In Section2, we will review the TMD evolution formalism of the TMDs, mostly following the approach established in [15]. Particularly, we will discuss in detail the extraction of the nonperturbative Sudakov form factor for the unpolarized TMD distribution of the proton/pion as well as that for the Sivers function. In Section 3, putting the evolved result of the TMD distributions into the TMD factorization formulae, we will present the theoretical expression of the physical observables, such as the unpolarized differential cross section, the Sivers asymmetry, and the cos⁡2ϕ asymmetry contributed by the double Boer-Mulders effect. In Section 4, we present the numerical evolution results of the unpolarized TMD distributions and the Boer-Mulders function of the pion meson, as well as that of the Sivers function of the proton. In Section 5, we display the phenomenology of the physical observables (unpolarized differential cross section, the Sivers asymmetry, and the cos⁡2ϕ asymmetry) in the πN Drell-Yan with TMD factorization at the kinematical configuration of the COMPASS experiments. We summarize the paper in Section 6.


## 2. The TMD Evolution of the Distribution Functions

In this section, we present a review on the TMD evolution of the distribution functions. Particularly, we provide the evolution formalism for the unpolarized distribution functionf1, transversity h1, Sivers function f1⊥, and the Boer-Mulders function h1⊥ of the proton, as well as f1 and h1⊥ of the pion meson, within the Collins-11 TMD factorization scheme [15].In general, it is more convenient to solve the evolution equations for the TMD distributions in the coordinate space (b space) other than that in the transverse momentum k⊥ space, with b conjugate to k⊥ via Fourier transformation [12, 15]. The TMD distributions F~(x,b;μ,ζF) in b space have two kinds of energy dependence, namely, μ is the renormalization scale related to the corresponding collinear PDFs, and ζF is the energy scale serving as a cutoff to regularize the light-cone singularity in the operator definition of the TMD distributions. Here, F is a shorthand for any TMD distribution function and the tilde denotes that the distribution is the one in b space. If we perform the inverse Fourier transformation on F~(x,b;μ,ζF), we recover the distribution function in the transverse momentum space Fq/H(x,k⊥;μ,ζF), which contains the information about the probability of finding a quark with specific polarization, collinear momentum fraction x, and transverse momentum k⊥ in a specifically polarized hadron H.


### 2.1. TMD Evolution Equations

The energy evolution for theζF dependence of the TMD distributions is encoded in the Collins-Soper (CS) [12, 15, 63] equation:(1)∂ln⁡F~x,b;μ,ζF∂ζF=K~b;μ,while the μ dependence is driven by the renormalization group equation as(2)dK~dln⁡μ=-γKαsμ,(3)dln⁡F~x,b;μ,ζFdln⁡μ=γFαsμ;ζF2μ2,with αs being the strong coupling at the energy scale μ, K~ being the CS evolution kernel, and γK, γF being the anomalous dimensions. The solutions of these evolution equations were studied in detail in [15, 63, 64]. Here, we will only discuss the final result. The overall structure of the solution for F~(x,b;μ,ζF) is similar to that for the Sudakov form factor. More specifically, the energy evolution of TMD distributions from an initial energy μ to another energy Q is encoded in the Sudakov-like form factor S by the exponential form exp(-S)(4)F~x,b,Q=F×e-S×F~x,b,μ,where F is the factor related to the hard scattering. Hereafter, we will set μ=ζF=Q and express F~(x,b;μ=Q,ζF=Q2) as F~(x,b;Q).As theb-dependence of the TMDs can provide very useful information regarding the transverse momentum dependence of the hadronic 3D structure through Fourier transformation, it is of fundamental importance to study the TMDs in b space. In the small b region, the b dependence is perturbatively calculable, while in the large b region, the dependence turns to be nonperturbative and may be obtained from the experimental 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 bmax serving as the boundary between the two regions. The prescription also allows for a smooth transition from perturbative to nonperturbative regions and avoids the Landau pole singularity in αs(μb). A b-dependent function b∗ is defined to have the property b∗≈b at low values of b and b∗≈bmax at large b values. In this paper, we adopt the original CSS prescription [12]:(5)b∗=b1+b2/bmax2,bmax<1ΛQCD.The typical value of bmax is chosen around 1GeV-1 to guarantee that b∗ is always in the perturbative region. Besides the CSS prescription, there were several different prescriptions in literature. In [65, 66] a function bmin(b) decreasing with increasing 1/Q was also introduced to match the TMD factorization with the fixed-order collinear calculations in the very small b region.In the smallb region 1/Q≪b≪1/ΛQCD, the TMD distributions at fixed energy μ can be expressed as the convolution of the perturbatively calculable coefficients and the corresponding collinear PDFs or the multiparton correlation functions [22, 67](6)F~q/Hx,b;μ=∑iCq←i⊗Fi/Hx,μ.Here, ⊗ stands for the convolution in the momentum fraction x(7)Cq←i⊗f1i/Hx,μ≡∫x1dξξCq←ixξ,b;μf1i/Hξ,μand fi/H(x,μ) is the corresponding collinear counterpart of flavor i in hadron H at the energy scale μ. The latter one could be a dynamic scale related to b∗ by μb=c0/b∗, with c0=2e-γE and the Euler Constant γE≈0.577 [22]. The perturbative hard coefficients Cq←i, independent of the initial hadron type, have been calculated for the parton-target case [23, 68] as the series of (αs/π) and the results have been presented in [67] (see also Appendix A of [23]).


### 2.2. Sudakov Form Factors for the Proton and the Pion

The Sudakov-like form factorS in (4) can be separated into the perturbatively calculable part SP and the nonperturbative part SNP(8)S=SP+SNP.According to the studies in [26, 39, 69–71], the perturbative part of the Sudakov form factor SP has the same result among different kinds of distribution functions, i.e., SP is spin-independent. It has the general form(9)SPQ,b∗=∫μb2Q2dμ¯2μ¯2Aαsμ¯ln⁡Q2μ¯2+Bαsμ¯.The coefficients A and B in(9) can be expanded as the series of αs/π:(10)A=∑n=1∞Anαsπn,(11)B=∑n=1∞Bnαsπn.Here, we list A(n) to A(2) and B(n) to B(1) up to the accuracy of next-to-leading-logarithmic (NLL) order [12, 23, 26, 69, 72, 73]:(12)A1=CF(13)A2=CF2CA6718-π26-109TRnf(14)B1=-32CF.For the nonperturbative form factor SNP, it can not be analytically calculated by the perturbative method, which means it has to be parameterized to obtain the evolution information in the nonperturbative region.The general form ofSNP(Q;b) was suggested as [12](15)SNPQ;b=g2bln⁡QQ0+g1b.The nonperturbative functions g1(b) and g2(b) are functions of the impact parameter b and depend on the choice of bmax. To be more specific, g2(b) contains the information on the large b behavior of the evolution kernel K~. Also, according to the power counting analysis in [74], g2(b) shall follow the power behavior as b2 at small-b region, which can be an essential constraint for the parameterization of g2(b). The well-known Brock-Landry-Nadolsky-Yuan (BLNY) fit parameterizes g2(b) as g2b2 with g2 a free parameter [72]. We note that g2(b) is universal for different types of TMDs and does not depend on the particular process, which is an important prediction of QCD factorization theorems involving TMDs [15, 23, 39, 75]. The nonperturbative function g1(b) contains information on the intrinsic nonperturbative transverse motion of bound partons, namely, it should depend on the type of hadron and the quark flavor as well as x for TMD distributions. As for the TMD fragmentation functions, it may depend on zh, the type of the produced hadron, and the quark flavor. In other words, g1(b) depends on the specific TMDs.There are several extractions forSNP in literature, we review some often-used forms below.The original BLNY fit parameterizedSNP as [72](16)g1+g2ln⁡Q2Q0+g1g3ln100x1x2b2,where x1 and x2 are the longitudinal momentum fractions of the incoming hadrons carried by the initial state quark and antiquark. The BLNY parameterization proved to be very reliable to describe Drell-Yan data and W±,Z boson production [72]. However, when the parameterization is extrapolated to the typical SIDIS kinematics in HERMES and COMPASS, the transverse momentum distribution of hadron can not be described by the BLNY-type fit [76, 77].Inspired by [72, 78], a widely used parameterization of SNP for TMD distributions or fragmentation functions was proposed [39, 67, 72, 78–80](17)SNPpdf/ff=b2g1pdf/ff+g22ln⁡QQ0,where the factor 1/2 in front of g2 comes from the fact that only one hadron is involved for the parameterization of SNPpdf/ff, while the parameter in [78] is for pp collisions. The parameter g1pdf/ff in (17) depends on the type of TMDs, which can be regarded as the width of intrinsic transverse momentum for the relevant TMDs at the initial energy scale Q0 [23, 73, 81]. Assuming a Gaussian form, one can obtain(18)g1pdf=k⊥2Q04,g1ff=pT2Q04z2,where 〈k⊥2〉Q0 and 〈pT2〉Q0 represent the relevant averaged intrinsic transverse momenta squared for TMD distributions and TMD fragmentation functions at the initial scale Q0, respectively.Since the original BLNY fit fails to simultaneously describe Drell-Yan process and SIDIS process, in [77] the authors proposed a new form for SNP which releases the tension between the BLNY fit to the Drell-Yan (such as W, Z and low energy Drell-Yan pair productions) data and the fit to the SIDIS data from HERMES/COMPASS in the CSS resummation formalism. In addition, the x-dependence in (16) was separated with a power law behavior assumption: (x0/x)λ, where x0 and λ are the fixed parameters as x0=0.01 and λ=0.2. The two different behaviors (logarithmic in (16) and power law) will differ in the intermediate x regime. Reference [76] showed that a direct integration of the evolution kernel from low Q to high Q led to the form of ln(Q) term as ln(b/b∗)ln(Q) and could describe the SIDIS and Drell-Yan data with Q values ranging from a few GeV to 10 GeV. Thus, the g2(b) term was modified to the form of ln(b/b∗) and the functional form of SNP extracted in [77] turned to the form(19)g1b2+g2ln⁡bb∗ln⁡QQ0+g3b2x0x1λ+x0x2λ.At small b region (b is much smaller than bmax), the parameterization of the g2(b) term g2ln(b/b∗) can be approximated as b2/(2bmax2), which satisfied the constraint of the b2 behavior for g2(b). However, at large b region, the logarithmic behavior will lead to different predictions on the Q2 dependence, since the Gaussian-type parameterization suggests that it is strongly suppressed [82]. This form has been suggested in an early research by Collins and Soper [83], but has not yet been adopted in any phenomenological study until the study in [77]. The comparison between the original BLNY parameterization and this form with the experimental data of Drell-Yan type process has shown that the new form of SNP can fit with the data as equally well as the original BLNY parameterization.In [66], the g2(b) function was parameterized as g2b2, following the BLNY convention. Furthermore, in the function g1(b), the Gaussian width also depends on x. The authors simultaneously fit the experimental data of SIDIS process from HERMES and COMPASS Collaborations, the Drell-Yan events at low energy, and the Z boson production with totally 8059 data points. The extraction can describe the data well in the regions where TMD factorization is supposed to hold.To study the pion-nucleon Drell-Yan data, it is also necessary to know the nonperturbative Sudakov form factor for the pion meson. In [59], we extended the functional form for the proton TMDs [77] to the case of the pion TMDs:(20)SNPf1q/π=g1q/πb2+g2q/πln⁡bb∗ln⁡QQ0,with g1q/π and g2q/π the free parameters. Adopting the functional form of SNP in (20), for the first time, we performed the extraction [59] of the nonperturbative Sudakov form factor for the unpolarized TMD PDF of pion meson using the experimental data in the π-p Drell-Yan process collected by the E615 Collaboration at Fermilab [7, 84]. The data fitting was performed by the package MINUIT [85, 86], through a least-squares fit:(21)χ2α=∑i=1M∑j=1Nitheoq⊥ij,α-dataij2errij2.The total number of data in our fit is N=∑i8Ni=96. Since the TMD formalism is valid in the region q⊥≪Q, we did a simple data selection by removing the data in the region q⊥>3GeV. We performed the fit by minimizing the chi-square in (21), and we obtained the following values for the two parameters: (22)g1q/π=0.082±0.022,g2q/π=0.394±0.103,with χ2/d.o.f=1.64.Figure1 plots the q⊥-dependent differential cross section (solid line) calculated from the fitted values for g1q/π and g2q/π in (22) at the kinematics of E615 at different xF bins. The full squares with error bars denote the E615 data for comparison. As Figure 1 demonstrates, a good fit is obtained in the region xF<0.8.Figure 1
The fitted cross section (solid line) of pion-nucleon Drell-Yan as functions ofq⊥, compared with the E615 data (full square), for different xF bins in the range 0<xF<0.8. The error bars shown here include the statistical error and the 16% systematic error. Figure from [59].From the fitted result, we find that the value of the parameterg1q/π is smaller than the parameter g1q/p extracted in [77] which used the same parameterized form. For the parameter g2q/π we find that its value is very close to that of the parameter g2q/p for the proton [77] (here g2q/p=g2/2=0.42). This may confirm that g2 should be universal, e.g., g2 is independent on the hadron type. Similar to the case of the proton, for the pion meson g2π is several times larger than g1π. We note that a form of SNPf1,q/π motivated by the NJL model was given in [87].


### 2.3. Solutions for Different TMDs

After solving the evolution equations and incorporating the Sudakov form factor, the scale-dependent TMD distribution functionF~ of the proton and the pion in b space can be rewritten as(23)F~q/px,b;Q=e-1/2SPQ,b∗-SNPFq/pQ,bFαsQ∑iCq←i⊗Fi/px,μb,(24)F~q/πx,b;Q=e-1/2SPQ,b∗-SNPFq/πQ,bFαsQ∑iCq←i⊗Fi/πx,μb.Here, Fi/H(x,μb) is the corresponding collinear distributions at the initial energy scale μb. To be more specific, for the unpolarized distribution function f1,q/H and transversity distribution function h1,q/H, the collinear distributions Fi/H(x,μb) are the integrated distribution functions f1,q/H(x,μb) and h1,q/H(x,μb). As for the Boer-Mulders function and Sivers function, the collinear distributions are the corresponding multiparton correlation functions. Thus, the unpolarized distribution function of the proton and pion in b space can be written as(25)f~1,q/px,b;Q=e-1/2SpertQ,b∗-SNPf1,q/pQ,bFαsQ∑iCq←i⊗f1,i/px,μb(26)f~1,q/πx,b;Q=e-1/2SpertQ,b∗-SNPf1,q/πQ,bFαsQ∑iCq←i⊗f1,i/πx,μb.If we perform a Fourier transformation on the f~1,q/H(x,b;Q), we can obtain the distribution function in k⊥ space as(27)f1,q/px,k⊥;Q=∫0∞dbb2πJ0k⊥bf~1,q/px,b;Q,(28)f1,q/πx,k⊥;Q=∫0∞dbb2πJ0k⊥bf~1,q/πx,b;Q.where J0 is the Bessel function of the first kind, and k⊥=k⊥.Similarly, the evolution formalism of the proton transversity distribution inb space and k⊥-space can be obtained as [75](29)h~1,q/px,b;Q=e-1/2SPQ,b∗-SNPh1Q,bHαsQ∑iδCq←i⊗h1,i/px,μb,(30)h1,q/px,k⊥;Q=∫0∞dbb2πJ0k⊥bh~1,q/px,b;Q,where H is the hard factor, and δCq←i is the coefficient convoluted with the transversity. The TMD evolution formalism in (30) has been applied in [75] to extract the transversity distribution from the SIDIS data.The Sivers function and Boer-Mulders function, which are T-odd, can be expressed as follows inb-space [39](31)f~1T,q/H⊥αDYx,b;μ,ζF=∫d2k⊥e-ik→⊥·b→k⊥αMpf1T,q/H⊥DYx,k⊥;μ,h~1,q/H⊥αDYx,b;μ,ζF=∫d2k⊥e-ik→⊥·b→k⊥αMph1,q/H⊥DYx,k⊥;μ.Here, the superscript “DY" represents the distributions in the Drell-Yan process. Since QCD predicts that the sign of the distributions changes in the SIDIS process and Drell-Yan process, for the distributions in SIDIS process, there has to be an extra minus sign regard to f1T,q/H⊥(DY) and h1,q/H⊥(DY).Similar to what has been done to the unpolarized distribution function and transversity distribution function, in the lowb region, the Sivers function f~1T,q/H⊥α(DY) can also be expressed as the convolution of perturbatively calculable hard coefficients and the corresponding collinear correlation functions as [69, 88](32)f~1T,q/H⊥αDYx,b;μ=-ibα2∑iΔCq←iT⊗fi/p3x′,x′′;μ.Here, fi/p3(x′,x′′) denotes the twist-three quark-gluon-quark or trigluon correlation functions, among which the transverse spin-dependent Qiu-Sterman matrix element Tq,Fx′,x′′ [89–91] is the most relevant one. Assuming that the Qiu-Sterman function Tq,F(x,x) is the main contribution, the Sivers function in b-space becomes(33)f~1T,q/H⊥αDYx,b;Q=-ibα2FSivαsQ∑iΔCq←iT⊗Ti/H,Fx,x;μbe-SNPsiv-1/2SP,where FSiv is the factor related to the hard scattering. The Boer-Mulders function in b-space follows the similar result for the Sivers function as:(34)h~1,q/Hα⊥DYx,b;Q=-ibα2HBMαsQ∑iCq←iBM⊗Ti/H,Fσx,x;μbe-SNPBM-1/2SP,Here, Cq←iBM stands for the flavor-dependent hard coefficients convoluted with Ti/H,F(σ), HBM the hard scattering factor and Ti/H,F(σ)(x,x;μb) denotes the chiral-odd twist-3 collinear correlation function. After performing the Fourier transformation back to the transverse momentum space, one can get the Sivers function and the Boer-Mulders function as(35)k⊥MHf1T,q/H⊥x,k⊥;Q=∫0∞dbb22πJ1k⊥bFSivαsQ∑iΔCq←iT⊗f1T,i/H⊥1x,μbe-SNPsiv-1/2SP,(36)k⊥MHh1,q/H⊥x,k⊥;Q=∫0∞dbb22πJ1k⊥bHBMαsQ∑iCq←iBM⊗h1,i/H⊥1x;μbe-SNPBM-1/2SP,and Tq,F(x,x;μb) and Ti/H,F(σ)(x,x;μb) are related to Sivers function and Boer-Mulders function as [69, 88](37)Tq/H,Fx,x=∫d2k⊥k⊥2MHf1T,q/H⊥DYx,k⊥=2MHf1T,q/H⊥1DYx,(38)Tq/H,Fσx,x=∫d2k⊥k⊥2MHh1,q/H⊥DYx,k⊥=2MHh1,q/H⊥1DYx.The TMD evolution formalism in (35) has been applied to extract [39, 70, 81, 92, 93] the Sivers function. The similar formalism in (36) could be used to improve the previous extractions of the proton Boer-Mulders function [62, 94–96] and future extraction of the pion Boer-Mulders function.


## 2.1. TMD Evolution Equations

The energy evolution for theζF dependence of the TMD distributions is encoded in the Collins-Soper (CS) [12, 15, 63] equation:(1)∂ln⁡F~x,b;μ,ζF∂ζF=K~b;μ,while the μ dependence is driven by the renormalization group equation as(2)dK~dln⁡μ=-γKαsμ,(3)dln⁡F~x,b;μ,ζFdln⁡μ=γFαsμ;ζF2μ2,with αs being the strong coupling at the energy scale μ, K~ being the CS evolution kernel, and γK, γF being the anomalous dimensions. The solutions of these evolution equations were studied in detail in [15, 63, 64]. Here, we will only discuss the final result. The overall structure of the solution for F~(x,b;μ,ζF) is similar to that for the Sudakov form factor. More specifically, the energy evolution of TMD distributions from an initial energy μ to another energy Q is encoded in the Sudakov-like form factor S by the exponential form exp(-S)(4)F~x,b,Q=F×e-S×F~x,b,μ,where F is the factor related to the hard scattering. Hereafter, we will set μ=ζF=Q and express F~(x,b;μ=Q,ζF=Q2) as F~(x,b;Q).As theb-dependence of the TMDs can provide very useful information regarding the transverse momentum dependence of the hadronic 3D structure through Fourier transformation, it is of fundamental importance to study the TMDs in b space. In the small b region, the b dependence is perturbatively calculable, while in the large b region, the dependence turns to be nonperturbative and may be obtained from the experimental 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 bmax serving as the boundary between the two regions. The prescription also allows for a smooth transition from perturbative to nonperturbative regions and avoids the Landau pole singularity in αs(μb). A b-dependent function b∗ is defined to have the property b∗≈b at low values of b and b∗≈bmax at large b values. In this paper, we adopt the original CSS prescription [12]:(5)b∗=b1+b2/bmax2,bmax<1ΛQCD.The typical value of bmax is chosen around 1GeV-1 to guarantee that b∗ is always in the perturbative region. Besides the CSS prescription, there were several different prescriptions in literature. In [65, 66] a function bmin(b) decreasing with increasing 1/Q was also introduced to match the TMD factorization with the fixed-order collinear calculations in the very small b region.In the smallb region 1/Q≪b≪1/ΛQCD, the TMD distributions at fixed energy μ can be expressed as the convolution of the perturbatively calculable coefficients and the corresponding collinear PDFs or the multiparton correlation functions [22, 67](6)F~q/Hx,b;μ=∑iCq←i⊗Fi/Hx,μ.Here, ⊗ stands for the convolution in the momentum fraction x(7)Cq←i⊗f1i/Hx,μ≡∫x1dξξCq←ixξ,b;μf1i/Hξ,μand fi/H(x,μ) is the corresponding collinear counterpart of flavor i in hadron H at the energy scale μ. The latter one could be a dynamic scale related to b∗ by μb=c0/b∗, with c0=2e-γE and the Euler Constant γE≈0.577 [22]. The perturbative hard coefficients Cq←i, independent of the initial hadron type, have been calculated for the parton-target case [23, 68] as the series of (αs/π) and the results have been presented in [67] (see also Appendix A of [23]).


## 2.2. Sudakov Form Factors for the Proton and the Pion

The Sudakov-like form factorS in (4) can be separated into the perturbatively calculable part SP and the nonperturbative part SNP(8)S=SP+SNP.According to the studies in [26, 39, 69–71], the perturbative part of the Sudakov form factor SP has the same result among different kinds of distribution functions, i.e., SP is spin-independent. It has the general form(9)SPQ,b∗=∫μb2Q2dμ¯2μ¯2Aαsμ¯ln⁡Q2μ¯2+Bαsμ¯.The coefficients A and B in(9) can be expanded as the series of αs/π:(10)A=∑n=1∞Anαsπn,(11)B=∑n=1∞Bnαsπn.Here, we list A(n) to A(2) and B(n) to B(1) up to the accuracy of next-to-leading-logarithmic (NLL) order [12, 23, 26, 69, 72, 73]:(12)A1=CF(13)A2=CF2CA6718-π26-109TRnf(14)B1=-32CF.For the nonperturbative form factor SNP, it can not be analytically calculated by the perturbative method, which means it has to be parameterized to obtain the evolution information in the nonperturbative region.The general form ofSNP(Q;b) was suggested as [12](15)SNPQ;b=g2bln⁡QQ0+g1b.The nonperturbative functions g1(b) and g2(b) are functions of the impact parameter b and depend on the choice of bmax. To be more specific, g2(b) contains the information on the large b behavior of the evolution kernel K~. Also, according to the power counting analysis in [74], g2(b) shall follow the power behavior as b2 at small-b region, which can be an essential constraint for the parameterization of g2(b). The well-known Brock-Landry-Nadolsky-Yuan (BLNY) fit parameterizes g2(b) as g2b2 with g2 a free parameter [72]. We note that g2(b) is universal for different types of TMDs and does not depend on the particular process, which is an important prediction of QCD factorization theorems involving TMDs [15, 23, 39, 75]. The nonperturbative function g1(b) contains information on the intrinsic nonperturbative transverse motion of bound partons, namely, it should depend on the type of hadron and the quark flavor as well as x for TMD distributions. As for the TMD fragmentation functions, it may depend on zh, the type of the produced hadron, and the quark flavor. In other words, g1(b) depends on the specific TMDs.There are several extractions forSNP in literature, we review some often-used forms below.The original BLNY fit parameterizedSNP as [72](16)g1+g2ln⁡Q2Q0+g1g3ln100x1x2b2,where x1 and x2 are the longitudinal momentum fractions of the incoming hadrons carried by the initial state quark and antiquark. The BLNY parameterization proved to be very reliable to describe Drell-Yan data and W±,Z boson production [72]. However, when the parameterization is extrapolated to the typical SIDIS kinematics in HERMES and COMPASS, the transverse momentum distribution of hadron can not be described by the BLNY-type fit [76, 77].Inspired by [72, 78], a widely used parameterization of SNP for TMD distributions or fragmentation functions was proposed [39, 67, 72, 78–80](17)SNPpdf/ff=b2g1pdf/ff+g22ln⁡QQ0,where the factor 1/2 in front of g2 comes from the fact that only one hadron is involved for the parameterization of SNPpdf/ff, while the parameter in [78] is for pp collisions. The parameter g1pdf/ff in (17) depends on the type of TMDs, which can be regarded as the width of intrinsic transverse momentum for the relevant TMDs at the initial energy scale Q0 [23, 73, 81]. Assuming a Gaussian form, one can obtain(18)g1pdf=k⊥2Q04,g1ff=pT2Q04z2,where 〈k⊥2〉Q0 and 〈pT2〉Q0 represent the relevant averaged intrinsic transverse momenta squared for TMD distributions and TMD fragmentation functions at the initial scale Q0, respectively.Since the original BLNY fit fails to simultaneously describe Drell-Yan process and SIDIS process, in [77] the authors proposed a new form for SNP which releases the tension between the BLNY fit to the Drell-Yan (such as W, Z and low energy Drell-Yan pair productions) data and the fit to the SIDIS data from HERMES/COMPASS in the CSS resummation formalism. In addition, the x-dependence in (16) was separated with a power law behavior assumption: (x0/x)λ, where x0 and λ are the fixed parameters as x0=0.01 and λ=0.2. The two different behaviors (logarithmic in (16) and power law) will differ in the intermediate x regime. Reference [76] showed that a direct integration of the evolution kernel from low Q to high Q led to the form of ln(Q) term as ln(b/b∗)ln(Q) and could describe the SIDIS and Drell-Yan data with Q values ranging from a few GeV to 10 GeV. Thus, the g2(b) term was modified to the form of ln(b/b∗) and the functional form of SNP extracted in [77] turned to the form(19)g1b2+g2ln⁡bb∗ln⁡QQ0+g3b2x0x1λ+x0x2λ.At small b region (b is much smaller than bmax), the parameterization of the g2(b) term g2ln(b/b∗) can be approximated as b2/(2bmax2), which satisfied the constraint of the b2 behavior for g2(b). However, at large b region, the logarithmic behavior will lead to different predictions on the Q2 dependence, since the Gaussian-type parameterization suggests that it is strongly suppressed [82]. This form has been suggested in an early research by Collins and Soper [83], but has not yet been adopted in any phenomenological study until the study in [77]. The comparison between the original BLNY parameterization and this form with the experimental data of Drell-Yan type process has shown that the new form of SNP can fit with the data as equally well as the original BLNY parameterization.In [66], the g2(b) function was parameterized as g2b2, following the BLNY convention. Furthermore, in the function g1(b), the Gaussian width also depends on x. The authors simultaneously fit the experimental data of SIDIS process from HERMES and COMPASS Collaborations, the Drell-Yan events at low energy, and the Z boson production with totally 8059 data points. The extraction can describe the data well in the regions where TMD factorization is supposed to hold.To study the pion-nucleon Drell-Yan data, it is also necessary to know the nonperturbative Sudakov form factor for the pion meson. In [59], we extended the functional form for the proton TMDs [77] to the case of the pion TMDs:(20)SNPf1q/π=g1q/πb2+g2q/πln⁡bb∗ln⁡QQ0,with g1q/π and g2q/π the free parameters. Adopting the functional form of SNP in (20), for the first time, we performed the extraction [59] of the nonperturbative Sudakov form factor for the unpolarized TMD PDF of pion meson using the experimental data in the π-p Drell-Yan process collected by the E615 Collaboration at Fermilab [7, 84]. The data fitting was performed by the package MINUIT [85, 86], through a least-squares fit:(21)χ2α=∑i=1M∑j=1Nitheoq⊥ij,α-dataij2errij2.The total number of data in our fit is N=∑i8Ni=96. Since the TMD formalism is valid in the region q⊥≪Q, we did a simple data selection by removing the data in the region q⊥>3GeV. We performed the fit by minimizing the chi-square in (21), and we obtained the following values for the two parameters: (22)g1q/π=0.082±0.022,g2q/π=0.394±0.103,with χ2/d.o.f=1.64.Figure1 plots the q⊥-dependent differential cross section (solid line) calculated from the fitted values for g1q/π and g2q/π in (22) at the kinematics of E615 at different xF bins. The full squares with error bars denote the E615 data for comparison. As Figure 1 demonstrates, a good fit is obtained in the region xF<0.8.Figure 1
The fitted cross section (solid line) of pion-nucleon Drell-Yan as functions ofq⊥, compared with the E615 data (full square), for different xF bins in the range 0<xF<0.8. The error bars shown here include the statistical error and the 16% systematic error. Figure from [59].From the fitted result, we find that the value of the parameterg1q/π is smaller than the parameter g1q/p extracted in [77] which used the same parameterized form. For the parameter g2q/π we find that its value is very close to that of the parameter g2q/p for the proton [77] (here g2q/p=g2/2=0.42). This may confirm that g2 should be universal, e.g., g2 is independent on the hadron type. Similar to the case of the proton, for the pion meson g2π is several times larger than g1π. We note that a form of SNPf1,q/π motivated by the NJL model was given in [87].


## 2.3. Solutions for Different TMDs

After solving the evolution equations and incorporating the Sudakov form factor, the scale-dependent TMD distribution functionF~ of the proton and the pion in b space can be rewritten as(23)F~q/px,b;Q=e-1/2SPQ,b∗-SNPFq/pQ,bFαsQ∑iCq←i⊗Fi/px,μb,(24)F~q/πx,b;Q=e-1/2SPQ,b∗-SNPFq/πQ,bFαsQ∑iCq←i⊗Fi/πx,μb.Here, Fi/H(x,μb) is the corresponding collinear distributions at the initial energy scale μb. To be more specific, for the unpolarized distribution function f1,q/H and transversity distribution function h1,q/H, the collinear distributions Fi/H(x,μb) are the integrated distribution functions f1,q/H(x,μb) and h1,q/H(x,μb). As for the Boer-Mulders function and Sivers function, the collinear distributions are the corresponding multiparton correlation functions. Thus, the unpolarized distribution function of the proton and pion in b space can be written as(25)f~1,q/px,b;Q=e-1/2SpertQ,b∗-SNPf1,q/pQ,bFαsQ∑iCq←i⊗f1,i/px,μb(26)f~1,q/πx,b;Q=e-1/2SpertQ,b∗-SNPf1,q/πQ,bFαsQ∑iCq←i⊗f1,i/πx,μb.If we perform a Fourier transformation on the f~1,q/H(x,b;Q), we can obtain the distribution function in k⊥ space as(27)f1,q/px,k⊥;Q=∫0∞dbb2πJ0k⊥bf~1,q/px,b;Q,(28)f1,q/πx,k⊥;Q=∫0∞dbb2πJ0k⊥bf~1,q/πx,b;Q.where J0 is the Bessel function of the first kind, and k⊥=k⊥.Similarly, the evolution formalism of the proton transversity distribution inb space and k⊥-space can be obtained as [75](29)h~1,q/px,b;Q=e-1/2SPQ,b∗-SNPh1Q,bHαsQ∑iδCq←i⊗h1,i/px,μb,(30)h1,q/px,k⊥;Q=∫0∞dbb2πJ0k⊥bh~1,q/px,b;Q,where H is the hard factor, and δCq←i is the coefficient convoluted with the transversity. The TMD evolution formalism in (30) has been applied in [75] to extract the transversity distribution from the SIDIS data.The Sivers function and Boer-Mulders function, which are T-odd, can be expressed as follows inb-space [39](31)f~1T,q/H⊥αDYx,b;μ,ζF=∫d2k⊥e-ik→⊥·b→k⊥αMpf1T,q/H⊥DYx,k⊥;μ,h~1,q/H⊥αDYx,b;μ,ζF=∫d2k⊥e-ik→⊥·b→k⊥αMph1,q/H⊥DYx,k⊥;μ.Here, the superscript “DY" represents the distributions in the Drell-Yan process. Since QCD predicts that the sign of the distributions changes in the SIDIS process and Drell-Yan process, for the distributions in SIDIS process, there has to be an extra minus sign regard to f1T,q/H⊥(DY) and h1,q/H⊥(DY).Similar to what has been done to the unpolarized distribution function and transversity distribution function, in the lowb region, the Sivers function f~1T,q/H⊥α(DY) can also be expressed as the convolution of perturbatively calculable hard coefficients and the corresponding collinear correlation functions as [69, 88](32)f~1T,q/H⊥αDYx,b;μ=-ibα2∑iΔCq←iT⊗fi/p3x′,x′′;μ.Here, fi/p3(x′,x′′) denotes the twist-three quark-gluon-quark or trigluon correlation functions, among which the transverse spin-dependent Qiu-Sterman matrix element Tq,Fx′,x′′ [89–91] is the most relevant one. Assuming that the Qiu-Sterman function Tq,F(x,x) is the main contribution, the Sivers function in b-space becomes(33)f~1T,q/H⊥αDYx,b;Q=-ibα2FSivαsQ∑iΔCq←iT⊗Ti/H,Fx,x;μbe-SNPsiv-1/2SP,where FSiv is the factor related to the hard scattering. The Boer-Mulders function in b-space follows the similar result for the Sivers function as:(34)h~1,q/Hα⊥DYx,b;Q=-ibα2HBMαsQ∑iCq←iBM⊗Ti/H,Fσx,x;μbe-SNPBM-1/2SP,Here, Cq←iBM stands for the flavor-dependent hard coefficients convoluted with Ti/H,F(σ), HBM the hard scattering factor and Ti/H,F(σ)(x,x;μb) denotes the chiral-odd twist-3 collinear correlation function. After performing the Fourier transformation back to the transverse momentum space, one can get the Sivers function and the Boer-Mulders function as(35)k⊥MHf1T,q/H⊥x,k⊥;Q=∫0∞dbb22πJ1k⊥bFSivαsQ∑iΔCq←iT⊗f1T,i/H⊥1x,μbe-SNPsiv-1/2SP,(36)k⊥MHh1,q/H⊥x,k⊥;Q=∫0∞dbb22πJ1k⊥bHBMαsQ∑iCq←iBM⊗h1,i/H⊥1x;μbe-SNPBM-1/2SP,and Tq,F(x,x;μb) and Ti/H,F(σ)(x,x;μb) are related to Sivers function and Boer-Mulders function as [69, 88](37)Tq/H,Fx,x=∫d2k⊥k⊥2MHf1T,q/H⊥DYx,k⊥=2MHf1T,q/H⊥1DYx,(38)Tq/H,Fσx,x=∫d2k⊥k⊥2MHh1,q/H⊥DYx,k⊥=2MHh1,q/H⊥1DYx.The TMD evolution formalism in (35) has been applied to extract [39, 70, 81, 92, 93] the Sivers function. The similar formalism in (36) could be used to improve the previous extractions of the proton Boer-Mulders function [62, 94–96] and future extraction of the pion Boer-Mulders function.


## 3. Physical Observables inπN Drell-Yan Process within TMD Factorization

In this section we will set up the necessary framework for physical observables inπ-N Drell-Yan process within TMD factorization by considering the evolution effects of the TMD distributions, following the procedure developed in [15].In Drell-Yan process(39)HAPπ+HBPN→γ∗q+X→l+l+l-l′+X,Pπ/N and q denote the momenta of the incoming hadron π/N and the virtual photon, respectively; q is a time-like vector, namely, Q2=q2>0, which is the invariant mass square of the final-state lepton pair. One can define the following useful kinematical variables to express the cross section:(40)s=Pπ+PN2,xπ/N=Q22Pπ/N·q,xF=2qLs=xπ-xN,τ=Q2s=xπxN,y=12ln⁡q+q-=12ln⁡xπxN,where s is the center-of-mass energy squared; xπ/N is the light-front momentum fraction carried by the annihilating quark/antiquark in the incoming hadron π/N; qL is the longitudinal momentum of the virtual photon in the c.m. frame of the incident hadrons; xF 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π/N is expressed as the function of xF, τ and of y, τ(41)xπ/N=±xF+xF2+4τ2,xπ/N=τe±y.


### 3.1. Differential Cross Section for Unpolarized Drell-Yan Process

The differential cross section formulated in TMD factorization is usually expressed in theb-space to guarantee conservation of the transverse momenta of the emitted soft gluons. Later on it can be transformed back to the transverse momentum space to represent the experimental observables. We will introduce the physical observables in the following part of this section.The general differential cross section for the unpolarized Drell-Yan process can be written as [12](42)d4σUUdQ2dyd2q⊥=σ0∫d2b2π2eiq→⊥·b→W~UUQ;b+YUUQ,q⊥where σ0=4παem2/3NCsQ2 is the cross section at tree level with αem the fine-structure constant, W~(Q;b) is the structure function in the b-space which contains all-order resummation results and dominates in the low q⊥ region (q⊥≪Q); and the Y term provides necessary correction at q⊥~Q. 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 (42).In general, TMD factorization [15] aims at separating well-defined TMD distributions such that they can be used in different processes through a universal way and expressing the scheme/process dependence in the corresponding hard factors. Thus, W~(Q;b) can be expressed as [97](43)W~UUQ;b=HUUQ;μ∑q,q¯eq2f~q/πsubxπ,b;μ,ζFf~q/psubxp,b;μ,ζF,where f~q/Hsub is the subtracted distribution function in the b space and HUU(Q;μ) is the factor associated with hard scattering. The superscript “sub" represents the distribution function with the soft factor subtracted. The subtraction guarantees the absence of light-cone singularities in the TMDs and the self-energy divergencies of the soft factors [15, 22]. However, the way to subtract the soft factor in the distribution function and the hard factor HUU(Q;μ) depends on the scheme to regulate the light-cone singularity in the TMD definition [12, 14, 15, 22, 98–103], leading to the scheme dependence in the TMD factorization. In literature, several different schemes are used [97]: the CSS scheme [12, 22], the Collins-11 (JCC) scheme [15], the Ji-Ma-Yuan (JMY) scheme [13, 14], and the lattice scheme [103]. Although different schemes are adopted, the final results of the structure functions W~(Q;b) as well as the differential cross section should not depend on a specific scheme. In the following we will apply the JCC and JMY schemes to display the scheme-independence of the unpolarized differential cross section.The hardHUU(Q;μ) have different forms in the JCC and JMY schemes:(44)HJCCQ;μ=1+αsμ2πCF3ln⁡Q2μ2-ln2Q2μ2+π2-8,(45)HJMYQ;μ,ρ=1+αsμ2πCF1+ln⁡ρ2ln⁡Q2μ2-ln⁡ρ2+ln2ρ+2π2-4.Like ζF, here ρ is another variable to regulate the light-cone singularity of TMD distributions. The scheme dependence of the distribution function is manifested in the hard factor F(αs(Q)), which has the following forms in different schemes:(46)F~JCCαsQ=1+Oαs2,(47)F~JMYαsQ,ρ=1+αs2πCFln⁡ρ-12ln2ρ-π22-2,The C coefficients in (25) and (26) do not depend on the types of initial hadrons and are calculated for the parton-target case [23, 68] with the results presented in [67] (see also Appendix A of [23])(48)Cq←q′x,b;μ,ζF=δqq′δ1-x+αsπCF21-x,(49)Cq←gx,b;μ,ζF=αsπTRx1-x,where CF=(NC2-1)/(2NC), TR=1/2.One can absorb the scheme-dependent hard factorsHUU(Q;μ) and F of the TMD distributions into the C-functions using(50)Cj←i′=Cj←i×F×HUUQ;μ=Q.The results for the splitting to quark are(51)Cq←q′′x,b;μb=δqq′δ1-x+αsπCF21-x+CF4π2-8δ1-x,(52)Cq←g′x,b;μb=αsπTRx1-x.The new C-coefficients turn out to be scheme independent (independent on ρ) [104] but process dependent [105, 106].With the newC-coefficients in hand, one can obtain the structure functions W~UU(Q;b) in b-space as(53)W~UUQ;b=e-SpertQ2,b-SNPf1q/πQ2,b-SNPf1q/pQ2,b×∑q,q¯eq2Cq←i′⊗fi/π-x1,μbCq¯←j′⊗fj/px2,μb.After performing the Fourier transformation, we can get the differential cross section as(54)d4σdQ2dyd2q⊥=σ0∫0∞dbb2πJ0q⊥b×W~UUQ;b,where J0 is the zeroth order Bessel function of the first kind.


### 3.2. The Sivers Asymmetry

In the Drell-Yan process with aπ beam colliding on the transversely polarized nucleon target, an important physical observable is the Sivers asymmetry, as it can test the sign change of the Sivers function between SIDIS and Drell-Yan processes, a fundamental prediction in QCD. The future precise measurement of the Sivers asymmetry in πN Drell-Yan in a wide kinematical region can be also used to extract the Sivers function. The Sivers asymmetry is usually defined as [39](55)AUT=d4Δσ/dQ2dyd2q⊥d4σ/dQ2dyd2q⊥,where d4σ/dQ2dyd2q⊥ and d4Δσ/dQ2dyd2q⊥ are the spin-averaged (unpolarized) and spin-dependent differential cross section, respectively. The latter one has the general form in the TMD factorization [39, 69, 88](56)d4ΔσdQ2dyd2q⊥=σ0ϵ⊥αβS⊥α∫d2b2π2eiq→⊥·b→W~UTβQ;b+YUTβQ,q⊥.Similar to (42), W~UT(Q,b) denotes the spin-dependent structure function in the b-space and dominates at q⊥≪Q, and YUTβ provides correction for the single-polarized process at q⊥~Q. The antisymmetric tensor ϵ⊥αβ is defined as ϵαβμνPπμPpν/Pπ·Pp, and S⊥ is the transverse-polarization vector of the proton target.The structure functionW~UT(Q,b) can be written in terms of the unpolarized distribution function of pion and Sivers function of proton as(57)W~UTαQ;b=HUTQ;μ∑q,q¯eq2f~1q¯/πxπ,b;μ,ζFf~1Tq/p⊥αDYxp,b;μ,ζF,with f~1Tq/p⊥α(DY)(xp,b;μ,ζF) given in (33). Similar to the unpolarized case, the scheme-dependent hard factors can be absorbed into the C-coefficients, leading to [69, 88](58)ΔCq←q′Tx,b;μb=δqq′δ1-x+αsπ-14Nc1-x+CF4π2-8δ1-x.The spin-dependent differential cross section in (56) thus has the form(59)d4ΔσdQ2dyd2q⊥=σ04π∫0∞dbb2J1q⊥b∑q,i,jeq2ΔCq←iT⊗Ti,Fxp,xp;μbCq¯←j⊗f1,j/πxπ,μbe-SNPSiv+SNPf1q/π+SP.Combing (55), (54), (59), one can get the Sivers asymmetry in the Drell-Yan process with a π beam colliding on a transversely polarized proton target.


### 3.3. Thecos⁡2ϕ Asymmetry in the Unpolarized Drell-Yan from Double Boer-Mulders Effect

The angular differential cross section for unpolarized Drell-Yan process has the following general form(60)1σdσdΩ=34π1λ+31+λcos2θ+μsin⁡2θcos⁡ϕ+ν2sin2θcos⁡2ϕ,where θ is the polar angle and ϕ is the azimuthal angle of the hadron plane with respect to the dilepton plane in the Collins-Soper (CS) frame [107]. The coefficients λ, μ, ν in (60) describe the sizes of different angular dependencies. Particularly, ν stands for the asymmetry of the cos⁡2ϕ azimuthal angular distribution of the dilepton.The coefficientsλ, μ, ν have been measured in the process π-N→μ+μ-X by the NA10 Collaboration [5, 6] and the E615 Collaboration [7] for a π- beam with energies of 140, 194, 286 GeV, and 252 GeV, with N denoting a nucleon in the deuterium or tungsten target. The experimental data showed a large value of ν, near 30% in the region QT~3 GeV. This demonstrates a clear violation of the Lam-Tung relation [47]. In the last decade the angular coefficients were also measured in the pN Drell-Yan processes in both the fixed-target mode [108, 109] and collider mode [48, 49]. The origin of large cos⁡2ϕ asymmetry–or the violation of the Lam-Tung relation– observed in Drell-Yan process has been studied extensively in literature [30, 50–57, 110–114]. Here we will only consider the contribution from the coupling of two Boer-Mulders functions from each hadron, denoted by νBM. It might be measured through the combination 2νBM≈2ν+λ-1, in which the perturbative contribution is largely subtracted.Thecos⁡2ϕ asymmetry coefficient νBM contributed by the Boer-Mulders function can be written as(61)νBM=2∑qF2h^·k⊥h^·p⊥-k⊥·p⊥h1,q/π⊥h1,q¯/p⊥/MπMp∑qFf1,q/πf1,q¯/p,where the notation(62)Fωff¯=eq2∫d2k⊥d2p⊥δ2k⊥+p⊥-q⊥ωfxπ,k⊥2f¯xp,p⊥2has been adopted to express the convolution of transverse momenta. q⊥, k⊥ and p⊥ are the transverse momenta of the lepton pair, quark, and antiquark in the initial hadrons, respectively. h^ is a unit vector defined as h^=q⊥/|q⊥|=q⊥/q⊥. One can perform the Fourier transformation from q⊥ space to b space on the delta function in the notation of (62) to obtain the denominator in (61) as(63)Ff1,q/πf1,q¯/p=∑qeq2∫d2b2π2∫d2k⊥d2p⊥eiqT-k⊥-p⊥·bf1,q/πxπ,k⊥2f1,q¯/pxp,p⊥2=∑qeq2∫d2b2π2eiq⊥·bf~1,q/πxπ,b;Qf~1,q¯/pxp,b;Qwhere the unpolarized distribution function in b space is given in (25) and (26). Similar to the treatment of the denominator, using the expression of the Boer-Mulders function in (34) the numerator can be obtained as(64)F2h^·k⊥h^·p⊥-k⊥·p⊥h1,q/π⊥h1,q¯/p⊥MπMp=∑qeq2∫d2b2π2∫d2k⊥d2p⊥eiqT-k⊥-p⊥·b2h^·k⊥h^·p⊥-k⊥·p⊥h1,q/π⊥h1,q¯/p⊥MπMp=∑qeq2∫d2b2π2eiq⊥·b2h^αh^β-gαβ⊥h~1,q/πα⊥xπ,b;Qh~1,q¯/pβ⊥xp,b;Q=∑qeq2∫0∞dbb38πJ2q⊥bTq/π,Fσxπ,xπ;μbTq¯/p,Fσxp,xp;μbe-SNPf1,q/p+SNPf1,q/π+SP.Different from the previous two cases, the hard coefficients CiBM and HBM for the Boer-Mulders function have not been calculated up to next-to-leading order (NLO), and still remain in leading order (LO) as Cq←i=δqiδ(1-x) and H=1.


## 3.1. Differential Cross Section for Unpolarized Drell-Yan Process

The differential cross section formulated in TMD factorization is usually expressed in theb-space to guarantee conservation of the transverse momenta of the emitted soft gluons. Later on it can be transformed back to the transverse momentum space to represent the experimental observables. We will introduce the physical observables in the following part of this section.The general differential cross section for the unpolarized Drell-Yan process can be written as [12](42)d4σUUdQ2dyd2q⊥=σ0∫d2b2π2eiq→⊥·b→W~UUQ;b+YUUQ,q⊥where σ0=4παem2/3NCsQ2 is the cross section at tree level with αem the fine-structure constant, W~(Q;b) is the structure function in the b-space which contains all-order resummation results and dominates in the low q⊥ region (q⊥≪Q); and the Y term provides necessary correction at q⊥~Q. 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 (42).In general, TMD factorization [15] aims at separating well-defined TMD distributions such that they can be used in different processes through a universal way and expressing the scheme/process dependence in the corresponding hard factors. Thus, W~(Q;b) can be expressed as [97](43)W~UUQ;b=HUUQ;μ∑q,q¯eq2f~q/πsubxπ,b;μ,ζFf~q/psubxp,b;μ,ζF,where f~q/Hsub is the subtracted distribution function in the b space and HUU(Q;μ) is the factor associated with hard scattering. The superscript “sub" represents the distribution function with the soft factor subtracted. The subtraction guarantees the absence of light-cone singularities in the TMDs and the self-energy divergencies of the soft factors [15, 22]. However, the way to subtract the soft factor in the distribution function and the hard factor HUU(Q;μ) depends on the scheme to regulate the light-cone singularity in the TMD definition [12, 14, 15, 22, 98–103], leading to the scheme dependence in the TMD factorization. In literature, several different schemes are used [97]: the CSS scheme [12, 22], the Collins-11 (JCC) scheme [15], the Ji-Ma-Yuan (JMY) scheme [13, 14], and the lattice scheme [103]. Although different schemes are adopted, the final results of the structure functions W~(Q;b) as well as the differential cross section should not depend on a specific scheme. In the following we will apply the JCC and JMY schemes to display the scheme-independence of the unpolarized differential cross section.The hardHUU(Q;μ) have different forms in the JCC and JMY schemes:(44)HJCCQ;μ=1+αsμ2πCF3ln⁡Q2μ2-ln2Q2μ2+π2-8,(45)HJMYQ;μ,ρ=1+αsμ2πCF1+ln⁡ρ2ln⁡Q2μ2-ln⁡ρ2+ln2ρ+2π2-4.Like ζF, here ρ is another variable to regulate the light-cone singularity of TMD distributions. The scheme dependence of the distribution function is manifested in the hard factor F(αs(Q)), which has the following forms in different schemes:(46)F~JCCαsQ=1+Oαs2,(47)F~JMYαsQ,ρ=1+αs2πCFln⁡ρ-12ln2ρ-π22-2,The C coefficients in (25) and (26) do not depend on the types of initial hadrons and are calculated for the parton-target case [23, 68] with the results presented in [67] (see also Appendix A of [23])(48)Cq←q′x,b;μ,ζF=δqq′δ1-x+αsπCF21-x,(49)Cq←gx,b;μ,ζF=αsπTRx1-x,where CF=(NC2-1)/(2NC), TR=1/2.One can absorb the scheme-dependent hard factorsHUU(Q;μ) and F of the TMD distributions into the C-functions using(50)Cj←i′=Cj←i×F×HUUQ;μ=Q.The results for the splitting to quark are(51)Cq←q′′x,b;μb=δqq′δ1-x+αsπCF21-x+CF4π2-8δ1-x,(52)Cq←g′x,b;μb=αsπTRx1-x.The new C-coefficients turn out to be scheme independent (independent on ρ) [104] but process dependent [105, 106].With the newC-coefficients in hand, one can obtain the structure functions W~UU(Q;b) in b-space as(53)W~UUQ;b=e-SpertQ2,b-SNPf1q/πQ2,b-SNPf1q/pQ2,b×∑q,q¯eq2Cq←i′⊗fi/π-x1,μbCq¯←j′⊗fj/px2,μb.After performing the Fourier transformation, we can get the differential cross section as(54)d4σdQ2dyd2q⊥=σ0∫0∞dbb2πJ0q⊥b×W~UUQ;b,where J0 is the zeroth order Bessel function of the first kind.


## 3.2. The Sivers Asymmetry

In the Drell-Yan process with aπ beam colliding on the transversely polarized nucleon target, an important physical observable is the Sivers asymmetry, as it can test the sign change of the Sivers function between SIDIS and Drell-Yan processes, a fundamental prediction in QCD. The future precise measurement of the Sivers asymmetry in πN Drell-Yan in a wide kinematical region can be also used to extract the Sivers function. The Sivers asymmetry is usually defined as [39](55)AUT=d4Δσ/dQ2dyd2q⊥d4σ/dQ2dyd2q⊥,where d4σ/dQ2dyd2q⊥ and d4Δσ/dQ2dyd2q⊥ are the spin-averaged (unpolarized) and spin-dependent differential cross section, respectively. The latter one has the general form in the TMD factorization [39, 69, 88](56)d4ΔσdQ2dyd2q⊥=σ0ϵ⊥αβS⊥α∫d2b2π2eiq→⊥·b→W~UTβQ;b+YUTβQ,q⊥.Similar to (42), W~UT(Q,b) denotes the spin-dependent structure function in the b-space and dominates at q⊥≪Q, and YUTβ provides correction for the single-polarized process at q⊥~Q. The antisymmetric tensor ϵ⊥αβ is defined as ϵαβμνPπμPpν/Pπ·Pp, and S⊥ is the transverse-polarization vector of the proton target.The structure functionW~UT(Q,b) can be written in terms of the unpolarized distribution function of pion and Sivers function of proton as(57)W~UTαQ;b=HUTQ;μ∑q,q¯eq2f~1q¯/πxπ,b;μ,ζFf~1Tq/p⊥αDYxp,b;μ,ζF,with f~1Tq/p⊥α(DY)(xp,b;μ,ζF) given in (33). Similar to the unpolarized case, the scheme-dependent hard factors can be absorbed into the C-coefficients, leading to [69, 88](58)ΔCq←q′Tx,b;μb=δqq′δ1-x+αsπ-14Nc1-x+CF4π2-8δ1-x.The spin-dependent differential cross section in (56) thus has the form(59)d4ΔσdQ2dyd2q⊥=σ04π∫0∞dbb2J1q⊥b∑q,i,jeq2ΔCq←iT⊗Ti,Fxp,xp;μbCq¯←j⊗f1,j/πxπ,μbe-SNPSiv+SNPf1q/π+SP.Combing (55), (54), (59), one can get the Sivers asymmetry in the Drell-Yan process with a π beam colliding on a transversely polarized proton target.


## 3.3. Thecos⁡2ϕ Asymmetry in the Unpolarized Drell-Yan from Double Boer-Mulders Effect

The angular differential cross section for unpolarized Drell-Yan process has the following general form(60)1σdσdΩ=34π1λ+31+λcos2θ+μsin⁡2θcos⁡ϕ+ν2sin2θcos⁡2ϕ,where θ is the polar angle and ϕ is the azimuthal angle of the hadron plane with respect to the dilepton plane in the Collins-Soper (CS) frame [107]. The coefficients λ, μ, ν in (60) describe the sizes of different angular dependencies. Particularly, ν stands for the asymmetry of the cos⁡2ϕ azimuthal angular distribution of the dilepton.The coefficientsλ, μ, ν have been measured in the process π-N→μ+μ-X by the NA10 Collaboration [5, 6] and the E615 Collaboration [7] for a π- beam with energies of 140, 194, 286 GeV, and 252 GeV, with N denoting a nucleon in the deuterium or tungsten target. The experimental data showed a large value of ν, near 30% in the region QT~3 GeV. This demonstrates a clear violation of the Lam-Tung relation [47]. In the last decade the angular coefficients were also measured in the pN Drell-Yan processes in both the fixed-target mode [108, 109] and collider mode [48, 49]. The origin of large cos⁡2ϕ asymmetry–or the violation of the Lam-Tung relation– observed in Drell-Yan process has been studied extensively in literature [30, 50–57, 110–114]. Here we will only consider the contribution from the coupling of two Boer-Mulders functions from each hadron, denoted by νBM. It might be measured through the combination 2νBM≈2ν+λ-1, in which the perturbative contribution is largely subtracted.Thecos⁡2ϕ asymmetry coefficient νBM contributed by the Boer-Mulders function can be written as(61)νBM=2∑qF2h^·k⊥h^·p⊥-k⊥·p⊥h1,q/π⊥h1,q¯/p⊥/MπMp∑qFf1,q/πf1,q¯/p,where the notation(62)Fωff¯=eq2∫d2k⊥d2p⊥δ2k⊥+p⊥-q⊥ωfxπ,k⊥2f¯xp,p⊥2has been adopted to express the convolution of transverse momenta. q⊥, k⊥ and p⊥ are the transverse momenta of the lepton pair, quark, and antiquark in the initial hadrons, respectively. h^ is a unit vector defined as h^=q⊥/|q⊥|=q⊥/q⊥. One can perform the Fourier transformation from q⊥ space to b space on the delta function in the notation of (62) to obtain the denominator in (61) as(63)Ff1,q/πf1,q¯/p=∑qeq2∫d2b2π2∫d2k⊥d2p⊥eiqT-k⊥-p⊥·bf1,q/πxπ,k⊥2f1,q¯/pxp,p⊥2=∑qeq2∫d2b2π2eiq⊥·bf~1,q/πxπ,b;Qf~1,q¯/pxp,b;Qwhere the unpolarized distribution function in b space is given in (25) and (26). Similar to the treatment of the denominator, using the expression of the Boer-Mulders function in (34) the numerator can be obtained as(64)F2h^·k⊥h^·p⊥-k⊥·p⊥h1,q/π⊥h1,q¯/p⊥MπMp=∑qeq2∫d2b2π2∫d2k⊥d2p⊥eiqT-k⊥-p⊥·b2h^·k⊥h^·p⊥-k⊥·p⊥h1,q/π⊥h1,q¯/p⊥MπMp=∑qeq2∫d2b2π2eiq⊥·b2h^αh^β-gαβ⊥h~1,q/πα⊥xπ,b;Qh~1,q¯/pβ⊥xp,b;Q=∑qeq2∫0∞dbb38πJ2q⊥bTq/π,Fσxπ,xπ;μbTq¯/p,Fσxp,xp;μbe-SNPf1,q/p+SNPf1,q/π+SP.Different from the previous two cases, the hard coefficients CiBM and HBM for the Boer-Mulders function have not been calculated up to next-to-leading order (NLO), and still remain in leading order (LO) as Cq←i=δqiδ(1-x) and H=1.


## 4. Numerical Estimate for the TMD Distributions

Based on the TMD evolution formalism for the distributions set up in Section2, we will show the numerical results for the TMD distributions. Particularly attention will be paid on those of the pion meson, as those of the proton have been studied numerically in [23, 71, 77].


### 4.1. The Unpolarized TMD Distribution of the Pion Meson

In [59], the authors applied (26) and the extracted parameters g1π and g2π to quantitatively study the scale dependence of the unpolarized TMD distributions of the pion meson with the JCC Scheme. For the collinear unpolarized distribution function of the pion meson, the NLO SMRS parameterization [115] was chosen. The results are plotted in Figure 2, with the left and right panels showing the subtracted distribution in b space and k⊥ space, for fixed xπ=0.1, at three different energy scales: Q2=2.4GeV2 (dotted line), 10GeV2 (solid line), 1000GeV2 (dashed line). From the b-dependent plots, one can see that at the highest energy scale Q2=1000GeV2, the peak of the curve is in the low b region where b<bmax, since in this case the perturbative part of the Sudakov form factor dominates. However, at lower energy scales, e.g., Q2=10GeV2 and Q2=2.4GeV2, the peak of the b-dependent distribution function moves towards the higher b region, indicating that the nonperturbative part of the TMD evolution becomes important at lower energy scales. For the distribution in k⊥ space, at higher 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⊥. It is interesting to point out that the shapes of the pion TMD distribution at different scales are similar to those of the proton, namely, Fig. 8 in [77].Figure 2
Subtracted unpolarized TMD distribution of the pion meson for valence quarks inb-space (left panel) and k⊥-space (right panel), at energies: Q2=2.4GeV2 (dotted lines), Q2=10GeV2 (solid lines), and Q2=1000GeV2 (dashed lines). From [59].


### 4.2. The Sivers Function of the Proton

The scale dependence of the T-odd distributions, such as the Sivers function and the Boer-Mulders function, is more involved than that of the T-even distributions. This is because their collinear counterparts are the twist-3 multiparton correlation functions [39, 60, 61, 69, 88], for which the exact evolution equations are far more complicated than those for the unpolarized distribution function. In numerical calculation, some approximations on the evolution kernels are usually adopted.In [39], the Qiu-Sterman function was assumed to be proportional to f1, namely, it follows the same evolution kernel as that for f1. A different choice was adopted in [60], where the homogenous terms of the exact evolution kernel for the Qiu-Sterman function [88, 116–124] were included to deal with the scale dependence of Qiu-Sterman function:(65)PqqQS≈Pqqf1-Nc21+z21-z-Ncδ1-z,with Pqqf1 the evolution kernel of the unpolarized PDF(66)Pqqf1=431+z21-z++32δ1-z.To solve the QCD evolution numerically, we resort to the QCD evolution package HOPPET [125] and we custom the code to include the splitting function in (65). For a comparison, in Figure 3 we plot the TMD evolution of the Sivers function for proton in b space and the k⊥ space using the above-mentioned two approaches [60]. In this estimate, the next leading order C-coefficients ΔCq←iT was adopted from [69, 88] and the nonperturbative Sudakov form factor for the Sivers function of proton was adopted as the form in (17). The Sivers functions are presented at three different energy scales: Q2=2.4GeV2, Q2=10GeV2 and Q2=100GeV2. Similar to the result for f1, one can conclude from the curves that the perturbative Sudakov form factor dominated in the low b region at higher energy scales and the nonperturbative part of the TMD evolution became more important at lower energy scales. However, the k⊥ tendency of the Sivers function in the two approaches is different, which indicates that the scale dependence of the Qiu-Sterman function may play a role in the TMD evolution.Figure 3
Subtracted Sivers function for the up quarks in Drell-Yan inb-space (upper panels) and k⊥-space (lower panels), at energies: Q2=2.4GeV2 (solid lines), Q2=10GeV2 (dashed lines), and Q2=100GeV2 (dotted lines). The left panel shows the result from the same evolution kernel as that for f1 in (66); the right panel shows the result from an approximate evolution kernel for the Qiu-Sterman function in (65), respectively. Figure from [60].


### 4.3. The Boer-Mulders Function of the Pion Meson

The evolution of the Boer-Mulders function for the valence quark insideπ meson has been calculated from (34) and (36) in [61], in which the collinear twist-3 correlation function Tq,F(σ) at the initial energy scale was obtained by adopting a model result of the Boer-Mulders function of the pion meson calculated from the light-cone wave functions [126]. For the scale evolution of Tq,F(σ), the exact evolution effect has been studied in [116]. For our purpose, we only consider the homogenous term in the evolution kernel(67)PqqTq,Fσx≈ΔTPqqx-NCδ1-x,with ΔTPqq(x)=CF2z/(1-z)++3/2δ(1-x) being the evolution kernel for the transversity distribution function h1(x). We customize the original code of QCDNUM [127] to include the approximate kernel in (67). For the nonperturbative part of the Sudakov form factor associated with Boer-Mulders function, the information still remains unknown. The assumption that SNP for Boer-Mulders function is same as that for f1 can be a practical way to access the information of TMD evolution for Boer-Mulders function.We plot theb-dependent and kT-dependent Boer-Mulders function at x=0.1 in the left and right panels of Figure 4, respectively. In calculating h~1,q/π⊥(x,b;Q) in Figure 4, we have rewritten the Boer-Mulders function in b space as (68)h~1,q/π⊥x,b;Q=ibαπh~1,q/πα⊥x,b;Q.The three curves in each panel correspond to three different energy scales: Q2=0.25GeV2 (solid lines), Q2=10GeV2 (dashed lines), Q2=1000GeV2 (dotted lines). From the curves, we find that the TMD evolution effect of the Boer-Mulders function is significant and should be considered in phenomenological analysis. The result also indicates that the perturbative Sudakov form factor dominates in the low b region at higher energy scales and the nonperturbative part of the TMD evolution becomes more important at lower energy scales.Figure 4
The Boer-Mulders function foru-quark in b space (left panel) and kT space (right panel) considering three different energy scales: Q2=2.4GeV2 (solid lines), Q2=10GeV2 (dashed lines), and Q2=1000GeV2 (dotted lines). Figures from [61].In conclusion, we find that the tendency of the distributions is similar: the distribution is dominated by perturbative region inb space at large Q2, while at lower Q2 the distribution shifts to the large b region, indicating that the nonperturbative effects of TMD evolution become important. For the distributions in k⊥ space, as the value of Q2 increases, the distributions become wider with a perturbative tail, while at low values of Q2, the distributions resemble Gaussian-type parameterization. However, the widths of the transverse momentum differ among different distributions.


## 4.1. The Unpolarized TMD Distribution of the Pion Meson

In [59], the authors applied (26) and the extracted parameters g1π and g2π to quantitatively study the scale dependence of the unpolarized TMD distributions of the pion meson with the JCC Scheme. For the collinear unpolarized distribution function of the pion meson, the NLO SMRS parameterization [115] was chosen. The results are plotted in Figure 2, with the left and right panels showing the subtracted distribution in b space and k⊥ space, for fixed xπ=0.1, at three different energy scales: Q2=2.4GeV2 (dotted line), 10GeV2 (solid line), 1000GeV2 (dashed line). From the b-dependent plots, one can see that at the highest energy scale Q2=1000GeV2, the peak of the curve is in the low b region where b<bmax, since in this case the perturbative part of the Sudakov form factor dominates. However, at lower energy scales, e.g., Q2=10GeV2 and Q2=2.4GeV2, the peak of the b-dependent distribution function moves towards the higher b region, indicating that the nonperturbative part of the TMD evolution becomes important at lower energy scales. For the distribution in k⊥ space, at higher 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⊥. It is interesting to point out that the shapes of the pion TMD distribution at different scales are similar to those of the proton, namely, Fig. 8 in [77].Figure 2
Subtracted unpolarized TMD distribution of the pion meson for valence quarks inb-space (left panel) and k⊥-space (right panel), at energies: Q2=2.4GeV2 (dotted lines), Q2=10GeV2 (solid lines), and Q2=1000GeV2 (dashed lines). From [59].


## 4.2. The Sivers Function of the Proton

The scale dependence of the T-odd distributions, such as the Sivers function and the Boer-Mulders function, is more involved than that of the T-even distributions. This is because their collinear counterparts are the twist-3 multiparton correlation functions [39, 60, 61, 69, 88], for which the exact evolution equations are far more complicated than those for the unpolarized distribution function. In numerical calculation, some approximations on the evolution kernels are usually adopted.In [39], the Qiu-Sterman function was assumed to be proportional to f1, namely, it follows the same evolution kernel as that for f1. A different choice was adopted in [60], where the homogenous terms of the exact evolution kernel for the Qiu-Sterman function [88, 116–124] were included to deal with the scale dependence of Qiu-Sterman function:(65)PqqQS≈Pqqf1-Nc21+z21-z-Ncδ1-z,with Pqqf1 the evolution kernel of the unpolarized PDF(66)Pqqf1=431+z21-z++32δ1-z.To solve the QCD evolution numerically, we resort to the QCD evolution package HOPPET [125] and we custom the code to include the splitting function in (65). For a comparison, in Figure 3 we plot the TMD evolution of the Sivers function for proton in b space and the k⊥ space using the above-mentioned two approaches [60]. In this estimate, the next leading order C-coefficients ΔCq←iT was adopted from [69, 88] and the nonperturbative Sudakov form factor for the Sivers function of proton was adopted as the form in (17). The Sivers functions are presented at three different energy scales: Q2=2.4GeV2, Q2=10GeV2 and Q2=100GeV2. Similar to the result for f1, one can conclude from the curves that the perturbative Sudakov form factor dominated in the low b region at higher energy scales and the nonperturbative part of the TMD evolution became more important at lower energy scales. However, the k⊥ tendency of the Sivers function in the two approaches is different, which indicates that the scale dependence of the Qiu-Sterman function may play a role in the TMD evolution.Figure 3
Subtracted Sivers function for the up quarks in Drell-Yan inb-space (upper panels) and k⊥-space (lower panels), at energies: Q2=2.4GeV2 (solid lines), Q2=10GeV2 (dashed lines), and Q2=100GeV2 (dotted lines). The left panel shows the result from the same evolution kernel as that for f1 in (66); the right panel shows the result from an approximate evolution kernel for the Qiu-Sterman function in (65), respectively. Figure from [60].


## 4.3. The Boer-Mulders Function of the Pion Meson

The evolution of the Boer-Mulders function for the valence quark insideπ meson has been calculated from (34) and (36) in [61], in which the collinear twist-3 correlation function Tq,F(σ) at the initial energy scale was obtained by adopting a model result of the Boer-Mulders function of the pion meson calculated from the light-cone wave functions [126]. For the scale evolution of Tq,F(σ), the exact evolution effect has been studied in [116]. For our purpose, we only consider the homogenous term in the evolution kernel(67)PqqTq,Fσx≈ΔTPqqx-NCδ1-x,with ΔTPqq(x)=CF2z/(1-z)++3/2δ(1-x) being the evolution kernel for the transversity distribution function h1(x). We customize the original code of QCDNUM [127] to include the approximate kernel in (67). For the nonperturbative part of the Sudakov form factor associated with Boer-Mulders function, the information still remains unknown. The assumption that SNP for Boer-Mulders function is same as that for f1 can be a practical way to access the information of TMD evolution for Boer-Mulders function.We plot theb-dependent and kT-dependent Boer-Mulders function at x=0.1 in the left and right panels of Figure 4, respectively. In calculating h~1,q/π⊥(x,b;Q) in Figure 4, we have rewritten the Boer-Mulders function in b space as (68)h~1,q/π⊥x,b;Q=ibαπh~1,q/πα⊥x,b;Q.The three curves in each panel correspond to three different energy scales: Q2=0.25GeV2 (solid lines), Q2=10GeV2 (dashed lines), Q2=1000GeV2 (dotted lines). From the curves, we find that the TMD evolution effect of the Boer-Mulders function is significant and should be considered in phenomenological analysis. The result also indicates that the perturbative Sudakov form factor dominates in the low b region at higher energy scales and the nonperturbative part of the TMD evolution becomes more important at lower energy scales.Figure 4
The Boer-Mulders function foru-quark in b space (left panel) and kT space (right panel) considering three different energy scales: Q2=2.4GeV2 (solid lines), Q2=10GeV2 (dashed lines), and Q2=1000GeV2 (dotted lines). Figures from [61].In conclusion, we find that the tendency of the distributions is similar: the distribution is dominated by perturbative region inb space at large Q2, while at lower Q2 the distribution shifts to the large b region, indicating that the nonperturbative effects of TMD evolution become important. For the distributions in k⊥ space, as the value of Q2 increases, the distributions become wider with a perturbative tail, while at low values of Q2, the distributions resemble Gaussian-type parameterization. However, the widths of the transverse momentum differ among different distributions.


## 5. Numerical Estimate for the Physical Observables inπ-N Drell-Yan Process

Based on the general TMD factorization framework provided in Section3, we present several physical observables in π-N Drell-Yan process in this section.QCD predicts that the T-odd PDFs present generalized universality, i.e., the sign of the Sivers function measured in Drell-Yan process should be opposite to its sign measured in SIDIS [16, 34, 35] process. The verification of this sign change [36–41] 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, 42–45]. The COMPASS Collaboration has reported the first measurement of the Sivers asymmetry in the pion-induced Drell-Yan process, in which a π- beam was scattered off the transversely polarized NH3 target [11]. The polarized Drell-Yan data from COMPASS together with the previous measurement of the Sivers effect in the W- and Z-boson production from p↑p collision at RHIC [45] will provide the first evidence of the sign change of the Sivers function. As COMPASS experiment has almost the same setup [11, 46] for SIDIS and Drell-Yan process, it will provide the unique chance to explore the sign change since the uncertainties in the extraction of the Sivers function from the two kinds of measurements can be reduced.


### 5.1. The Normalized Cross Section for Unpolarizedπ-N Drell-Yan Process

The very first step to understand the Sivers asymmetry in theπ-N Drell-Yan process is to quantitatively estimate the differential cross section in the same process for unpolarized nucleon target with high accuracy, since it always appears in the denominator of the asymmetry definition. The differential cross section for unpolarized Drell-Yan process has been given in (54). Applying the extracted nonperturbative Sudakov form factor for pion meson in [59] and the extracted SNP for nucleon in [77], we estimated the normalized transverse momentum spectrum of the dimuon production in the pion-nucleon Drell-Yan process at COMPASS for different q⊥ bins with an interval of 0.2 GeV. The result is plotted in Figure 5. From the curves, one can find the theoretical estimate on the q⊥ distribution of the dimuon agreed with the COMPASS data fairly well in the small q⊥ region where the TMD factorization is supposed to hold. The comparison somehow confirms the validity of extraction of the nonperturbative Sudakov form factor for the unpolarized distribution f1π of pion meson, within the TMD factorization. This may indicate that the framework can also be extended to the study of the azimuthal asymmetries in the πN Drell-Yan process, such as the Sivers asymmetry and Boer-Mulders asymmetry. We should point out that at larger q⊥, the numerical estimate in [59] cannot describe the data, indicating that the perturbative correction from the YUU term may play an important role in the region q⊥~Q. Further study on the Y term is needed to provide a complete picture of the q⊥ distribution of lepton pairs from πN Drell-Yan in the whole q⊥ range.Figure 5
The transverse spectrum of lepton pair production in the unpolarized pion-nucleon Drell-Yan process, with anNH3 target at COMPASS. The dashed line represents the theoretical calculation in [59]. The solid line shows the experimental measurement at COMPASS [11]. Figure from [59].


### 5.2. The Sivers Asymmetry

In [39], the authors adopted the Gaussian form of the nonperturbative Sudakov form factor SNP in (17) and the leading order C coefficients to perform a global fit on the Sivers function from the experimental data at HERMES [128], COMPASS [129, 130], and Jefferson Lab (JLab) [131]. With the extracted Sivers function from SIDIS process at hand, they made predictions for the Sivers asymmetry in Drell-Yan lepton pair and W production at future planned Drell-Yan facilities at COMPASS [10], Fermilab [42, 43], and RHIC [44, 132], which can be compared to the future experimental measurements to test the sign change of the Sivers functions between SIDIS and Drell-Yan processes. The predictions were presented in Fig. 12 and 13 of [39].The TMD evolution effect of the Sivers asymmetry in SIDIS andpp Drell-Yan at low transverse momentum has also been studied in [88], in which a framework was built to match SIDIS and Drell-Yan and cover the TMD physics with Q2 from several GeV2 to 104GeV2 (for W/Z boson production). It has shown that the evolution equations derived by a direct integral of the CSS evolution kernel from low to high Q can describe well the transverse momentum distribution of the unpolarized cross sections in the Q2 range from 2 to 100 GeV2. With this approach, the transverse moment of the quark Sivers functions can be constrained from the combined analysis of the HERMES and COMPASS data on the Sivers asymmetries in SIDIS. Based on this result, [88] provided the predictions for the Sivers asymmetries in pp Drell-Yan, as well as in π-p Drell-Yan. The latter one has been measured by the COMPASS Collaboration, and the comparison showed that the theoretical result is consistent with data (Fig. 6 in [11]) within the error bar.With the numerical results of the TMD distributions in (57), the Sivers asymmetry AUTSiv as function of xp, xπ, xF, and q⊥ in π-p↑→μ+μ-+X in the kinematics of COMPASS Collaboration was calculated in [60], as shown in Figure 6. The magnitude of the asymmetry is around 0.05÷0.10, which is consistent with the COMPASS measurement (full squares in Figure 6) [11] within the uncertainties of the asymmetry. The different approaches dealing with the energy dependence of Qiu-Sterman function lead to different shapes of the asymmetry. Furthermore, the asymmetry from the approximate evolution kernel has a fall at larger q⊥, which is more compatible to the shape of q⊥-dependent asymmetry of measured by the COMPASS Collaboration. The study may indicate that, besides the TMD evolution effect, the scale dependence of the Qiu-Sterman function will also play a role in the interpretation of the experimental data.Figure 6
The Sivers asymmetry within TMD factorization for aπ- scattering off transversely polarized proton Drell-Yan process as functions of xp (upper left), xπ (upper right), xF (lower left), and q⊥ (lower right), compared with the COMPASS data [11]. Figure from [60].


### 5.3. Thecos⁡2ϕ Azimuthal Asymmetry

Using (64), the cos⁡2ϕ azimuthal asymmetry contributed by the double Boer-Mulders effect in the πN Drell-Yan process was analyzed in [61], in which the TMD evolution of the Boer-Mulders function was included. In this calculation, the Boer-Mulders function of the proton was chosen from the parameterization in [62] at the initial energy Q02=1GeV2. As mentioned in Section 4.3, the Boer-Mulders function of the pion was adopted from the model calculation in [126]. Here we plot the estimated asymmetry νBM as function of xp,xπ,xF and q⊥ in the kinematical region of COMPASS in Figure 7. The bands correspond to the uncertainty of the parameterization of the Boer-Mulders function of the proton [62]. We find from the plots that, in the TMD formalism, the cos⁡2ϕ azimuthal asymmetry in the unpolarized π-p Drell-Yan process contributed by the Boer-Mulders functions is around several percent. Although the uncertainty from the proton Boer-Mulders functions is rather large, the asymmetry is firmly positive in the entire kinematical region. The asymmetries as the functions of xp,xπ,xF show slight dependence on the variables, while the q⊥ dependent asymmetry shows increasing tendency along with the increasing q⊥ in the small q⊥ range where the TMD formalism is valid. The result in Figure 7 indicates that precise measurements on the Boer-Mulders asymmetry νBM as functions of xp,xπ,xF, and q⊥ can provide an opportunity to access the Boer-Mulders function of the pion meson. Furthermore, the work may also shed light on the proton Boer-Mulders function since the previous extractions on it were mostly performed without TMD evolution.Figure 7
Thecos⁡2ϕ azimuthal asymmetries νBM as the functions of xp (upper left), xπ (upper right), xF (lower left), and q⊥ (lower right) for the unpolarized πp Drell-Yan process considering the TMD evolution in the kinematical region of COMPASS. The shadow areas correspond to the uncertainty of the parameters in the parameterization of the Boer-Mulders function for proton in [62]. Figures from [61].


## 5.1. The Normalized Cross Section for Unpolarizedπ-N Drell-Yan Process

The very first step to understand the Sivers asymmetry in theπ-N Drell-Yan process is to quantitatively estimate the differential cross section in the same process for unpolarized nucleon target with high accuracy, since it always appears in the denominator of the asymmetry definition. The differential cross section for unpolarized Drell-Yan process has been given in (54). Applying the extracted nonperturbative Sudakov form factor for pion meson in [59] and the extracted SNP for nucleon in [77], we estimated the normalized transverse momentum spectrum of the dimuon production in the pion-nucleon Drell-Yan process at COMPASS for different q⊥ bins with an interval of 0.2 GeV. The result is plotted in Figure 5. From the curves, one can find the theoretical estimate on the q⊥ distribution of the dimuon agreed with the COMPASS data fairly well in the small q⊥ region where the TMD factorization is supposed to hold. The comparison somehow confirms the validity of extraction of the nonperturbative Sudakov form factor for the unpolarized distribution f1π of pion meson, within the TMD factorization. This may indicate that the framework can also be extended to the study of the azimuthal asymmetries in the πN Drell-Yan process, such as the Sivers asymmetry and Boer-Mulders asymmetry. We should point out that at larger q⊥, the numerical estimate in [59] cannot describe the data, indicating that the perturbative correction from the YUU term may play an important role in the region q⊥~Q. Further study on the Y term is needed to provide a complete picture of the q⊥ distribution of lepton pairs from πN Drell-Yan in the whole q⊥ range.Figure 5
The transverse spectrum of lepton pair production in the unpolarized pion-nucleon Drell-Yan process, with anNH3 target at COMPASS. The dashed line represents the theoretical calculation in [59]. The solid line shows the experimental measurement at COMPASS [11]. Figure from [59].


## 5.2. The Sivers Asymmetry

In [39], the authors adopted the Gaussian form of the nonperturbative Sudakov form factor SNP in (17) and the leading order C coefficients to perform a global fit on the Sivers function from the experimental data at HERMES [128], COMPASS [129, 130], and Jefferson Lab (JLab) [131]. With the extracted Sivers function from SIDIS process at hand, they made predictions for the Sivers asymmetry in Drell-Yan lepton pair and W production at future planned Drell-Yan facilities at COMPASS [10], Fermilab [42, 43], and RHIC [44, 132], which can be compared to the future experimental measurements to test the sign change of the Sivers functions between SIDIS and Drell-Yan processes. The predictions were presented in Fig. 12 and 13 of [39].The TMD evolution effect of the Sivers asymmetry in SIDIS andpp Drell-Yan at low transverse momentum has also been studied in [88], in which a framework was built to match SIDIS and Drell-Yan and cover the TMD physics with Q2 from several GeV2 to 104GeV2 (for W/Z boson production). It has shown that the evolution equations derived by a direct integral of the CSS evolution kernel from low to high Q can describe well the transverse momentum distribution of the unpolarized cross sections in the Q2 range from 2 to 100 GeV2. With this approach, the transverse moment of the quark Sivers functions can be constrained from the combined analysis of the HERMES and COMPASS data on the Sivers asymmetries in SIDIS. Based on this result, [88] provided the predictions for the Sivers asymmetries in pp Drell-Yan, as well as in π-p Drell-Yan. The latter one has been measured by the COMPASS Collaboration, and the comparison showed that the theoretical result is consistent with data (Fig. 6 in [11]) within the error bar.With the numerical results of the TMD distributions in (57), the Sivers asymmetry AUTSiv as function of xp, xπ, xF, and q⊥ in π-p↑→μ+μ-+X in the kinematics of COMPASS Collaboration was calculated in [60], as shown in Figure 6. The magnitude of the asymmetry is around 0.05÷0.10, which is consistent with the COMPASS measurement (full squares in Figure 6) [11] within the uncertainties of the asymmetry. The different approaches dealing with the energy dependence of Qiu-Sterman function lead to different shapes of the asymmetry. Furthermore, the asymmetry from the approximate evolution kernel has a fall at larger q⊥, which is more compatible to the shape of q⊥-dependent asymmetry of measured by the COMPASS Collaboration. The study may indicate that, besides the TMD evolution effect, the scale dependence of the Qiu-Sterman function will also play a role in the interpretation of the experimental data.Figure 6
The Sivers asymmetry within TMD factorization for aπ- scattering off transversely polarized proton Drell-Yan process as functions of xp (upper left), xπ (upper right), xF (lower left), and q⊥ (lower right), compared with the COMPASS data [11]. Figure from [60].


## 5.3. Thecos⁡2ϕ Azimuthal Asymmetry

Using (64), the cos⁡2ϕ azimuthal asymmetry contributed by the double Boer-Mulders effect in the πN Drell-Yan process was analyzed in [61], in which the TMD evolution of the Boer-Mulders function was included. In this calculation, the Boer-Mulders function of the proton was chosen from the parameterization in [62] at the initial energy Q02=1GeV2. As mentioned in Section 4.3, the Boer-Mulders function of the pion was adopted from the model calculation in [126]. Here we plot the estimated asymmetry νBM as function of xp,xπ,xF and q⊥ in the kinematical region of COMPASS in Figure 7. The bands correspond to the uncertainty of the parameterization of the Boer-Mulders function of the proton [62]. We find from the plots that, in the TMD formalism, the cos⁡2ϕ azimuthal asymmetry in the unpolarized π-p Drell-Yan process contributed by the Boer-Mulders functions is around several percent. Although the uncertainty from the proton Boer-Mulders functions is rather large, the asymmetry is firmly positive in the entire kinematical region. The asymmetries as the functions of xp,xπ,xF show slight dependence on the variables, while the q⊥ dependent asymmetry shows increasing tendency along with the increasing q⊥ in the small q⊥ range where the TMD formalism is valid. The result in Figure 7 indicates that precise measurements on the Boer-Mulders asymmetry νBM as functions of xp,xπ,xF, and q⊥ can provide an opportunity to access the Boer-Mulders function of the pion meson. Furthermore, the work may also shed light on the proton Boer-Mulders function since the previous extractions on it were mostly performed without TMD evolution.Figure 7
Thecos⁡2ϕ azimuthal asymmetries νBM as the functions of xp (upper left), xπ (upper right), xF (lower left), and q⊥ (lower right) for the unpolarized πp Drell-Yan process considering the TMD evolution in the kinematical region of COMPASS. The shadow areas correspond to the uncertainty of the parameters in the parameterization of the Boer-Mulders function for proton in [62]. Figures from [61].


## 6. Summary and Prospects

It has been a broad consensus that the study on the TMD observables will provide information on the partons’ intrinsic transverse motions inside a hadron. In the previous sections we have tried to substantiate this statement mainly focusing on the unpolarized and single-polarizedπp Drell-Yan process within the TMD factorization. In particular, we reviewed the extraction of the nonperturbative function from the Drell-Yan and SIDIS data in the evolution formalism of the TMD distributions. We also discussed the further applications of the TMD factorization in the phenomenology of unpolarized cross section, the Sivers asymmetry, and the cos⁡2ϕ azimuthal asymmetry in the πp Drell-Yan process. In summary, we have the following understanding on the πN Drell-Yan from the viewpoint of the TMD factorization:(i)
The extraction of nonperturbative Sudakov form factor from theπN Drell-Yan may shed light on the evolution (scale dependence) of the pion TMD distribution. The prediction on transverse momentum distribution of the dilepton in the small q⊥ region is compatible with the COMPASS measurement and may serve as a first step to study the spin/azimuthal asymmetry in the πN Drell-Yan process at COMPASS.(ii)
The precise measurement on the single-spin asymmetry in the kinematical region of COMPASS can provide great opportunity to access the Sivers function. Besides the TMD evolution effect, the choice of the scale dependence of the Qiu-Sterman function can affect the shape of the asymmetry and should be considered in the future extraction of the Sivers function.(iii)
Sizablecos⁡2ϕ asymmetry contributed by the convolution of the Boer-Mulders functions of the pion meson and the proton can still be observed at COMPASS after the TMD evolution effect is considered. Future data with higher accuracy may provide further constraint on the Boer-Mulders function of the pion meson as well as that of the proton.Although a lot of progress on the theoretical framework of the TMD factorization and TMD evolution has been made, the improvement is still necessary both from the perturbative and nonperturbative aspects. In the future, the study ofSNP based on more precise experimental data is needed, such as including the flavor dependence and hadron dependence on the functional form for SNP. From the viewpoint of the perturbative region, higher-order calculation of the hard factors and coefficients will improve the accuracy of the theoretical framework. Moreover, most of the numerical calculations are based on the approximation that the Y-term correction is negligible in the small transverse momentum region; the inclusion of this term in the future estimate could be done to test the magnitude of the term. In addition, the TMD factorization is suitable to describe the small transverse momentum physics, while the collinear factorization is suitable for the large transverse momentum or the integrated transverse momentum. The matching between the two factorization schemes to study the unpolarized and polarized process over the whole transverse momentum region may be also necessary [65, 133].



---

*Source: 6734293-2019-01-22.xml*