#### 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 the asymmetry contributed by the double Boer-Mulders effects. The corresponding phenomenology, particularly at the kinematical configuration of the COMPASS Drell-Yan facility, is displayed numerically.

#### 1. Introduction

After the first observation of the lepton pairs produced in 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 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 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 dilepton is much smaller than the mass 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 region (), the fixed-order calculations of the cross sections in the collinear picture fail, leading to large double logarithms of the type . It is necessary to resum such logarithmic contributions to all orders in the strong coupling 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 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], 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 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 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 Drell-Yan process.

Another important observable in the Drell-Yan process is the 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 asymmetry, which violates the Lam-Tung relation [47]. Similar violation has also been observed in the 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 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 asymmetry from the double Boer-Mulders effect. In order to quantitatively understand various spin/azimuthal asymmetries in the 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 Section 2, 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 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 asymmetry) in the 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 function , transversity , Sivers function , and the Boer-Mulders function of the proton, as well as and 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 ( space) other than that in the transverse momentum space, with conjugate to via Fourier transformation [12, 15]. The TMD distributions in space have two kinds of energy dependence, namely, is the renormalization scale related to the corresponding collinear PDFs, and is the energy scale serving as a cutoff to regularize the light-cone singularity in the operator definition of the TMD distributions. Here, is a shorthand for any TMD distribution function and the tilde denotes that the distribution is the one in space. If we perform the inverse Fourier transformation on , we recover the distribution function in the transverse momentum space , which contains the information about the probability of finding a quark with specific polarization, collinear momentum fraction , and transverse momentum in a specifically polarized hadron .

##### 2.1. TMD Evolution Equations

The energy evolution for the dependence of the TMD distributions is encoded in the Collins-Soper (CS) [12, 15, 63] equation:while the dependence is driven by the renormalization group equation aswith being the strong coupling at the energy scale , being the CS evolution kernel, and , 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 is similar to that for the Sudakov form factor. More specifically, the energy evolution of TMD distributions from an initial energy to another energy is encoded in the Sudakov-like form factor by the exponential form where is the factor related to the hard scattering. Hereafter, we will set and express as .

As the -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 space. In the small region, the dependence is perturbatively calculable, while in the large region, the dependence turns to be nonperturbative and may be obtained from the experimental data. To combine the perturbative information at small with the nonperturbative part at large , a matching procedure must be introduced with a parameter 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 . A -dependent function is defined to have the property at low values of and at large values. In this paper, we adopt the original CSS prescription [12]:The typical value of is chosen around to guarantee that is always in the perturbative region. Besides the CSS prescription, there were several different prescriptions in literature. In [65, 66] a function decreasing with increasing was also introduced to match the TMD factorization with the fixed-order collinear calculations in the very small region.

In the small region , 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]Here, stands for the convolution in the momentum fraction and is the corresponding collinear counterpart of flavor in hadron at the energy scale . The latter one could be a dynamic scale related to by , with and the Euler Constant [22]. The perturbative hard coefficients , independent of the initial hadron type, have been calculated for the parton-target case [23, 68] as the series of 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 factor in (4) can be separated into the perturbatively calculable part and the nonperturbative part According to the studies in [26, 39, 69–71], the perturbative part of the Sudakov form factor has the same result among different kinds of distribution functions, i.e., is spin-independent. It has the general formThe coefficients and in(9) can be expanded as the series of :Here, we list to and to up to the accuracy of next-to-leading-logarithmic (NLL) order [12, 23, 26, 69, 72, 73]:For the nonperturbative form factor , 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 of was suggested as [12]The nonperturbative functions and are functions of the impact parameter and depend on the choice of . To be more specific, contains the information on the large behavior of the evolution kernel . Also, according to the power counting analysis in [74], shall follow the power behavior as at small- region, which can be an essential constraint for the parameterization of . The well-known Brock-Landry-Nadolsky-Yuan (BLNY) fit parameterizes as with a free parameter [72]. We note that 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 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 for TMD distributions. As for the TMD fragmentation functions, it may depend on , the type of the produced hadron, and the quark flavor. In other words, depends on the specific TMDs.

There are several extractions for in literature, we review some often-used forms below.

The original BLNY fit parameterized as [72]where and 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 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 for TMD distributions or fragmentation functions was proposed [39, 67, 72, 78–80]where the factor in front of comes from the fact that only one hadron is involved for the parameterization of , while the parameter in [78] is for collisions. The parameter 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 [23, 73, 81]. Assuming a Gaussian form, one can obtainwhere and represent the relevant averaged intrinsic transverse momenta squared for TMD distributions and TMD fragmentation functions at the initial scale , 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 which releases the tension between the BLNY fit to the Drell-Yan (such as , 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 -dependence in (16) was separated with a power law behavior assumption: , where and are the fixed parameters as and . The two different behaviors (logarithmic in (16) and power law) will differ in the intermediate regime. Reference [76] showed that a direct integration of the evolution kernel from low to high led to the form of term as and could describe the SIDIS and Drell-Yan data with values ranging from a few GeV to 10 GeV. Thus, the term was modified to the form of and the functional form of extracted in [77] turned to the formAt small region ( is much smaller than ), the parameterization of the term can be approximated as , which satisfied the constraint of the behavior for . However, at large region, the logarithmic behavior will lead to different predictions on the 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 can fit with the data as equally well as the original BLNY parameterization.

In [66], the function was parameterized as , following the BLNY convention. Furthermore, in the function , the Gaussian width also depends on . The authors simultaneously fit the experimental data of SIDIS process from HERMES and COMPASS Collaborations, the Drell-Yan events at low energy, and the 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:with and the free parameters. Adopting the functional form of 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 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:The total number of data in our fit is . Since the TMD formalism is valid in the region , we did a simple data selection by removing the data in the region . We performed the fit by minimizing the chi-square in (21), and we obtained the following values for the two parameters: with .

Figure 1 plots the -dependent differential cross section (solid line) calculated from the fitted values for and in (22) at the kinematics of E615 at different 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 .

From the fitted result, we find that the value of the parameter is smaller than the parameter extracted in [77] which used the same parameterized form. For the parameter we find that its value is very close to that of the parameter for the proton [77] (here ). This may confirm that should be universal, e.g., is independent on the hadron type. Similar to the case of the proton, for the pion meson is several times larger than . We note that a form of 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 function of the proton and the pion in space can be rewritten asHere, is the corresponding collinear distributions at the initial energy scale . To be more specific, for the unpolarized distribution function and transversity distribution function , the collinear distributions are the integrated distribution functions and . 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 space can be written asIf we perform a Fourier transformation on the , we can obtain the distribution function in space aswhere is the Bessel function of the first kind, and .

Similarly, the evolution formalism of the proton transversity distribution in space and -space can be obtained as [75]where is the hard factor, and 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 in -space [39]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 and .

Similar to what has been done to the unpolarized distribution function and transversity distribution function, in the low region, the Sivers function can also be expressed as the convolution of perturbatively calculable hard coefficients and the corresponding collinear correlation functions as [69, 88]Here, denotes the twist-three quark-gluon-quark or trigluon correlation functions, among which the transverse spin-dependent Qiu-Sterman matrix element [89–91] is the most relevant one. Assuming that the Qiu-Sterman function is the main contribution, the Sivers function in -space becomeswhere is the factor related to the hard scattering. The Boer-Mulders function in -space follows the similar result for the Sivers function as:Here, stands for the flavor-dependent hard coefficients convoluted with , the hard scattering factor and 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 asand and are related to Sivers function and Boer-Mulders function as [69, 88]

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 Drell-Yan Process within TMD Factorization

In this section we will set up the necessary framework for physical observables in - 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 and denote the momenta of the incoming hadron and the virtual photon, respectively; is a time-like vector, namely, , 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:where is the center-of-mass energy squared; is the light-front momentum fraction carried by the annihilating quark/antiquark in the incoming hadron ; is the longitudinal momentum of the virtual photon in the c.m. frame of the incident hadrons; is the Feynman variable, which corresponds to the longitudinal momentum fraction carried by the lepton pair; and is the rapidity of the lepton pair. Thus, is expressed as the function of , and of ,

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

The differential cross section formulated in TMD factorization is usually expressed in the -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]where is the cross section at tree level with the fine-structure constant, is the structure function in the -space which contains all-order resummation results and dominates in the low region (); and the term provides necessary correction at . In this work we will neglect the -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, can be expressed as [97]where is the subtracted distribution function in the space and 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 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 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 hard have different forms in the JCC and JMY schemes:Like , 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 , which has the following forms in different schemes:The 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])where , .

One can absorb the scheme-dependent hard factors and of the TMD distributions into the -functions usingThe results for the splitting to quark areThe new -coefficients turn out to be scheme independent (independent on ) [104] but process dependent [105, 106].

With the new -coefficients in hand, one can obtain the structure functions in -space asAfter performing the Fourier transformation, we can get the differential cross section aswhere 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 Drell-Yan in a wide kinematical region can be also used to extract the Sivers function. The Sivers asymmetry is usually defined as [39]where and 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]Similar to (42), denotes the spin-dependent structure function in the -space and dominates at , and provides correction for the single-polarized process at . The antisymmetric tensor is defined as , and is the transverse-polarization vector of the proton target.

The structure function can be written in terms of the unpolarized distribution function of pion and Sivers function of proton aswith given in (33). Similar to the unpolarized case, the scheme-dependent hard factors can be absorbed into the -coefficients, leading to [69, 88]

The spin-dependent differential cross section in (56) thus has the formCombing (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. The 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 formwhere 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 azimuthal angular distribution of the dilepton.

The coefficients , , have been measured in the process 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 denoting a nucleon in the deuterium or tungsten target. The experimental data showed a large value of , near 30% in the region GeV. This demonstrates a clear violation of the Lam-Tung relation [47]. In the last decade the angular coefficients were also measured in the Drell-Yan processes in both the fixed-target mode [108, 109] and collider mode [48, 49]. The origin of large 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 . It might be measured through the combination , in which the perturbative contribution is largely subtracted.

The asymmetry coefficient contributed by the Boer-Mulders function can be written aswhere the notationhas been adopted to express the convolution of transverse momenta. , and are the transverse momenta of the lepton pair, quark, and antiquark in the initial hadrons, respectively. is a unit vector defined as . One can perform the Fourier transformation from space to space on the delta function in the notation of (62) to obtain the denominator in (61) aswhere the unpolarized distribution function in 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 asDifferent from the previous two cases, the hard coefficients and for the Boer-Mulders function have not been calculated up to next-to-leading order (NLO), and still remain in leading order (LO) as and .

#### 4. Numerical Estimate for the TMD Distributions

Based on the TMD evolution formalism for the distributions set up in Section 2, 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 and 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 space and space, for fixed , at three different energy scales: (dotted line), (solid line), (dashed line). From the -dependent plots, one can see that at the highest energy scale , the peak of the curve is in the low region where , since in this case the perturbative part of the Sudakov form factor dominates. However, at lower energy scales, e.g., and , the peak of the -dependent distribution function moves towards the higher region, indicating that the nonperturbative part of the TMD evolution becomes important at lower energy scales. For the distribution in space, at higher energy scale the distribution has a tail falling off slowly at large , while at lower energy scales the distribution function falls off rapidly with increasing . 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].

##### 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 , namely, it follows the same evolution kernel as that for . 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:with the evolution kernel of the unpolarized PDFTo 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 space and the space using the above-mentioned two approaches [60]. In this estimate, the next leading order -coefficients 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: , and . Similar to the result for , one can conclude from the curves that the perturbative Sudakov form factor dominated in the low region at higher energy scales and the nonperturbative part of the TMD evolution became more important at lower energy scales. However, the 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.

##### 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 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 , the exact evolution effect has been studied in [116]. For our purpose, we only consider the homogenous term in the evolution kernelwith being the evolution kernel for the transversity distribution function . 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 for Boer-Mulders function is same as that for can be a practical way to access the information of TMD evolution for Boer-Mulders function.

We plot the -dependent and -dependent Boer-Mulders function at in the left and right panels of Figure 4, respectively. In calculating in Figure 4, we have rewritten the Boer-Mulders function in space as The three curves in each panel correspond to three different energy scales: (solid lines), (dashed lines), (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 region at higher energy scales and the nonperturbative part of the TMD evolution becomes more important at lower energy scales.

In conclusion, we find that the tendency of the distributions is similar: the distribution is dominated by perturbative region in space at large , while at lower the distribution shifts to the large region, indicating that the nonperturbative effects of TMD evolution become important. For the distributions in space, as the value of increases, the distributions become wider with a perturbative tail, while at low values of , 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 - Drell-Yan Process

Based on the general TMD factorization framework provided in Section 3, we present several physical observables in - 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 target [11]. The polarized Drell-Yan data from COMPASS together with the previous measurement of the Sivers effect in the - and -boson production from 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 - Drell-Yan Process

The very first step to understand the Sivers asymmetry in the - 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 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 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 distribution of the dimuon agreed with the COMPASS data fairly well in the small 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 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 Drell-Yan process, such as the Sivers asymmetry and Boer-Mulders asymmetry. We should point out that at larger , the numerical estimate in [59] cannot describe the data, indicating that the perturbative correction from the term may play an important role in the region . Further study on the term is needed to provide a complete picture of the distribution of lepton pairs from Drell-Yan in the whole range.

##### 5.2. The Sivers Asymmetry

In [39], the authors adopted the Gaussian form of the nonperturbative Sudakov form factor in (17) and the leading order 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 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 and 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 from several to (for boson production). It has shown that the evolution equations derived by a direct integral of the CSS evolution kernel from low to high can describe well the transverse momentum distribution of the unpolarized cross sections in the range from 2 to 100 . 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 Drell-Yan, as well as in 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 as function of , , , and in in the kinematics of COMPASS Collaboration was calculated in [60], as shown in Figure 6. The magnitude of the asymmetry is around , 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 , which is more compatible to the shape of -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.

##### 5.3. The Azimuthal Asymmetry

Using (64), the azimuthal asymmetry contributed by the double Boer-Mulders effect in the 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 . 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 as function of and 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 azimuthal asymmetry in the unpolarized 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 show slight dependence on the variables, while the dependent asymmetry shows increasing tendency along with the increasing in the small range where the TMD formalism is valid. The result in Figure 7 indicates that precise measurements on the Boer-Mulders asymmetry as functions of , and 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.

#### 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 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 azimuthal asymmetry in the Drell-Yan process. In summary, we have the following understanding on the Drell-Yan from the viewpoint of the TMD factorization:(i)The extraction of nonperturbative Sudakov form factor from the 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 region is compatible with the COMPASS measurement and may serve as a first step to study the spin/azimuthal asymmetry in the 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)Sizable 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 of based on more precise experimental data is needed, such as including the flavor dependence and hadron dependence on the functional form for . 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 -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].

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work is partially supported by the NSFC (China) Grant 11575043 and by the Fundamental Research Funds for the Central Universities of China. X. Wang is supported by the NSFC (China) Grant 11847217 and the China Postdoctoral Science Foundation under Grant no. 2018M640680.