$\pi-N$ Drell-Yan process in TMD factorization

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


I. INTRODUCTION
After the first observation of the µ + µ − lepton pairs produced in p N collisions [1], the process was interpreted as a quark and an antiquark from each initial hadron annihilating into a virtual photon, which in turn decays into a lepton pair [2]. This explanation makes the process an ideal tool in order to explore the internal structure of both the beam and target hadrons. Since then, a wide range of studies on this (Drell-Yan) process have been carried out. In particular, the πN Drell-Yan process has the unique capability to pin down the partonic structure of the pion, which is an unstable particle 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][4][5][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 constrain considerably 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 NH 3 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 q ⊥ is much smaller than the mass Q of the virtual vector boson, thus 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 α ln 2 (q 2 ⊥ /Q 2 ). 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][14][15] in which the gauge-invariant [16][17][18][19] transverse momentum dependent (TMD) parton distribution functions or fragmentation functions (collectively called as 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 scatering (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][37][38][39][40][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][43][44][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 the cos 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 possess 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][51][52], the non-coplanarity effect [30,53] and the QCD radiative effects at higher order [54,55]. Another promising origin [56] for the violation of the Lam-Tung relation at low transverse momentum is the convolution of the two Boer-Mulders functions [57] 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 [56,57], 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 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 differential spin-averaged 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 Sec. II, we will review the TMD evolution formalism for several distributions, mostly following the approach established in Ref. [15]. Particularly, we will discuss in details 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. III, 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. IV, 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. V, 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 VI.

II. 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 f 1 , transversity h 1 , Sivers function f ⊥ 1 and the Boer-Mulders function h ⊥ 1 of the proton, as well as f 1 and h ⊥ 1 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 distributionsF (x, b; µ, ζ F ) in the 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 the b space. If we perform the reverse Fourier Transformation onF (x, b; µ, ζ F ), we recover the distribution function in the transverse momentum space F q/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.

A. TMD evolution equations
The energy evolution for the ζ F dependence of the TMD distributions is encoded in the Collins-Soper (CS) [12,15,58] equation: while the µ dependence is driven by the renormalization group equation as with α s the strong coupling at the energy scale µ,K the CS evolution kernel, and γ K , γ F the anomalous dimensions. The solutions of these evolution equations were studied in details in Ref. [15,58,59]. Here, we will only discuss the final result. The overall structure of the solution forF (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) where F is the factor related to the hard scattering. Hereafter, we will set µ = √ ζ F = Q and expressF (x, b; µ = Q, ζ F = Q 2 ) asF (x, b; Q).
As the b-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 the b space. In the small b region, the b dependence is perturbatively calculable, while in the large b region, the dependence turns to 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 b max 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 * ≈ b max at large b values. In this paper, we adopt the original CSS prescription [12]: The typical value of b max is chosen around 1 GeV −1 to guarantee that b * is always in the perturbative region. Besides the CSS prescription, there were several different prescriptions in literature. In Refs. [60,61] a function b min (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 small b 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,62] Here, ⊗ stands for the convolution in the momentum fraction x and F i/H (ξ, µ) 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 = c 0 /b * , with c 0 = 2e −γE and the Euler Constant γ E ≈ 0.577 [22]. The perturbative hard coefficients C q←i , independent of the initial hadron type, have been calculated for the parton-target case [23,63] as the series of (α s /π) and the results have been presented in Ref. [62] (see also Appendix A of Ref. [23]).

B. Sudakov form factors for the proton and the pion
The Sudakov-like form factor S in Eq. (4) can be separated into the perturbatively calculable part S P and the nonperturbative part S NP S = S P + S NP .
According to the studies in Refs. [26,39,[64][65][66], the perturbative part of the Sudakov form factor S P has the same result among different kinds of distribution functions, i.e., S P is spin-independent. It has the general form The coefficients A and B in Eq.(9) can be expanded as the series of α s /π.
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,64,67,68]: For the nonperturbative form factor S NP , 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 S NP (Q; b) was suggested as [12] The nonperturbative functions g 1 (b) and g 2 (b) are functions of the impact parameter b and depend on the choice of b max . To be more specific, g 2 (b) contains the information on the large b behavior of the evolution kernelK. Also, according to the power counting analysis in Ref. [69], g 2 (b) shall follow the power behavior as b 2 at small-b region, which can be an essential constraint for the parametrization of g 2 (b). The well-known Brock-Landry-Nadolsky-Yuan (BLNY) fit parameterizes g 2 (b) as g 2 b 2 with g 2 a free parameter [67]. We note that g 2 (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,70]. The nonperturbative function g 1 (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 z h , the type of the produced hadron, and the quark flavor. In other words, g 1 (b) depends on the specific TMDs. There are several extractions for S NP in literature, we review some often-used forms below. The original BLNY fit parameterized S NP as [67] (g 1 + g 2 ln(Q/2Q 0 ) + g 1 g 3 ln(100 where x 1 and x 2 are the longitudinal momentum fractions of the incoming hadrons carried by the initial state quark and antiquark. BLNY parameterization proved to be very reliable to describe Drell-Yan data and W ± , Z boson production [67]. However, when the parametrization 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 [71,72]. Inspired by Refs. [67,73], a widely used parametrization of S NP for TMD distributions or fragmentation functions was proposed [39,62,67,[73][74][75] S pdf/ff NP where the factor 1/2 in front of g 2 comes from the fact that only one hadron is involved for the parametrization of S pdf/ff NP , while the parameter in Ref. [73] is for pp collisions. The parameter g pdf/ff 1 in Eq. (17) depends on the type of TMDs, which can be regarded as the intrinsic transverse momentum width for the relevant TMDs at the initial energy scale Q 0 [23,68,76]. Assuming a Gaussian form, one can obtain where k 2 ⊥ Q0 and p 2 T Q0 represent the relevant averaged intrinsic transverse momenta squared for TMD distributions and TMD fragmentation functions at the initial scale Q 0 , respectively.
Since the original BLNY fit fails to simultaneously describe Drell-Yan process and SIDIS process, in Ref. [72] the authors proposed a new form for S NP which releases the tension between the BLNY fit to the Drell-Yan (such as W , Z and low energy Drell-Yan pair productions) data with the SIDIS data from HERMES/COMPASS in the CSS resummation formalism. In addition, the x-dependence in Eq. (16) was separated with a power law behavior assumption: (x 0 /x) λ , where x 0 and λ are the fixed parameters as x 0 = 0.01 and λ = 0.2. The two different behaviors (logarithmic in Eq. (16) and power law) will differ in the intermediate x regime. Ref. [71] 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 g 2 (b) term was modified to the form of ln(b/b * ) and the functional form of S NP extracted in Ref. [72] turned to the form At small b region (b is much smaller than b max ), the parametrization of the g 2 (b) term g 2 ln(b/b * ) can be approximated as b 2 /(2b 2 max ), which satisfied the constraint of the b 2 behavior for g 2 (b). However, at large b region, the logarithmic behavior will lead to different predictions on the Q 2 dependence, since the Gaussian-type parametrization suggests that it is strongly suppressed [77]. This form has been suggested in an early research by Collins and Soper [78], but has not yet been adopted in any phenomenological study until the study in Ref. [72]. The comparison between the original BLNY parametrization and this form with the experimental data of Drell-Yan type process has shown that the new form of S NP can fit with the data as equally well as the original BLNY parametrization.
In Ref. [61], the g 2 (b) function was parameterized as g 2 b 2 , following the BLNY convention. Furthermore, in the function g 1 (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 Z boson production with totally 8059 data points in the fitting procedure. 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 Ref. [79], we extended the functional form for the proton TMDs [72] to the case of the pion TMDs: with g q/π 1 and g q/π 2 the free parameters. Adopting the functional form of S NP in Eq. (20), for the first time, we performed the extraction [79] 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,80]. The data fitting was performed by the package MINUIT [81,82], through a least-squares fit: The total number of data in our fit is N = Since the TMD formalism is valid in the region q ⊥ ≪ Q, we do a simple data selection by removing the data in the region q ⊥ > 3 GeV. We performed the fit by minimizing the chisquare in Eq. (21), and we obtained the following values for the two parameters: with χ 2 /d.o.f = 1.64. Fig. 1 plots the q ⊥ -dependent differential cross section (solid line) calculated from the fitted values for g q/π 1 and g q/π 2 in Eq. (22) at the kinematics of E615 at different x F bins. The full squares with error bars denote the E615 data for comparison. As Fig. 1 demonstrates, a good fit is obtained at the region x F < 0.8.  From the fitted result, we find that the value of the parameter g q/π 1 is smaller than the parameter g q/p 1 extracted in Ref. [72] which used the same parameterized form. For the parameter g q/π 2 we find that its value is very close to that of the parameter g q/p 2 for the proton [72] (here g q/p 2 = g 2 /2 = 0.42). This may confirm that g 2 should be universal, e.g., g 2 is independent on the hadron type. Similar to the case of the proton, for the pion meson g π 2 is several times larger than g π 1 . We note that a form of S f 1,q/π NP motivated by the NJL model was given in Ref. [83].

C. 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 Here, is the corresponding collinear distributions at the initial energy scale µ b . To be more specific, for the unpolarized distribution function f 1,q/H and transversity distribution function h 1,q/H , the collinear distributions As for the Boer-Mulders function and Sivers function, the collinear distributions are the corresponding multi-parton correlation functions. Thus, the unpolarized distribution function of the proton and pion in the b space can be written as If we perform a Fourier transformation on thef 1,q/H (x, b; Q), we can obtain the distribution function in the transverse momentum space as where J 0 is the Bessel function of the first kind, and k ⊥ = |k ⊥ |. Similarly, the evolution formalism of the proton transversity distribution in the b-space and the k ⊥ -space can be obtained as [70] where H is the hard factor, and δC q←i is the coefficient convoluted with the transversity. The TMD evolution formalism in Eq. (30) has been applied in Ref. [70] 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 the b-space [39] f ⊥α(DY) 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 f 1T,q/H can also be expressed as the convolution of perturbatively calculable hard coefficients and the corresponding collinear correlation functions as [64,84] f ⊥α(DY) Here, f i/p (x ′ , x ′′ ) denotes the twist-three quark-gluon-quark or trigluon correlation functions, among which the transverse spin-dependent Qiu-Sterman matrix element T q,F (x ′ , x ′′ ) [85][86][87] is the most relevant one. Assuming that the Qiu-Sterman function T q,F (x, x) is the main contribution, the Sivers function in the b-space becomes where F Siv is the factor related to the hard scattering. The Boer-Mulders function in the b-space follows the similar result for the Sivers function as: Here, C BM q←i stands for the flavor-dependent hard coefficients convoluted with T and T q,F (x, x; µ b ) and T The TMD evolution formalism in Eq. (35) has been applied to extract [39,65,76,88,89] the Sivers function. The similar formalism in Eq. (36) could be used to improved the previous extractions of the proton Boer-Mulders function [90][91][92][93] and future extraction of the pion Boer-Mulders function.

III. 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 Ref. [15].
In Drell-Yan process P π/N and q denote the momenta of the incoming hadron π/N and the virtual photon, respectively; q is a time-like vector, namely, Q 2 = q 2 > 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: 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 ; q L is the longitudinal momentum of the virtual photon in the c.m. frame of the incident hadrons; x F is the Feynman x variable, which corresponds to the longitudinal momentum fraction carried by the lepton pair; and y is the rapidity of the lepton pair. Thus, x π/N is expressed as the function of x F , τ and of y, τ A. Differential cross section for unpolarized Drell-Yan process The differential cross section formulated in TMD factorization is usually expressed in the b-space to guarantee conservation of the transverse momenta of the emitted soft gluons. Later on it can be transformed back to 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 σ 0 = 4πα 2 em 3NC sQ 2 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 Eq. (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 [94] wheref sub q/H is the subtracted distribution function in b-space and H UU (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 H UU (Q; µ) depend on the scheme to regulate the light-cone singularity in the TMD definition [12,14,15,22,[95][96][97][98][99][100], leading to the scheme dependence in the TMD factorization. In literature, several different schemes are used [94]: the CSS scheme [12,22], the Collins-11 (JCC) scheme [15], the Ji-Ma-Yuan (JMY) scheme [13,14], and the lattice scheme [100]. 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 hard H UU (Q; µ) have different forms in the JCC and JMY schemes: 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:F The C coefficients in Eqs. (25) and (26) do not depend on the types of initial hadrons and are calculated for the parton-target case [23,63] with the results presented in Ref. [62] (see also Appendix A of Ref. [23]) where C F = (N 2 C − 1)/(2N C ), T R = 1/2. One can absorb the scheme-dependent hard factors H UU (Q; µ) and F of the TMD distributions into the C-functions using The results for the quark splitting are The new C-coefficients turn out to be scheme independent (independent on ρ) [101] but process dependent [102,103].
With the new C-coefficients in hands, one can obtain the structure functions W UU (Q; b) in b-space as After performing the Fourier transformation, we can get the differential cross section as where J 0 is the zeroth order Bessel function of the first kind.

B. The Sivers asymmetry
In the Drell-Yan process with π 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] where d 4 σ dQ 2 dyd 2 q ⊥ and d 4 ∆σ dQ 2 dyd 2 q ⊥ are the spin-averaged (unpolarized) and spin-dependent differential cross section, respectively. The latter one has the general form in the TMD factorization [39,64,84] Similar to Eq. (42), W UT (Q, b) denotes the spin-dependent structure function in b-space and dominates at q ⊥ ≪ Q, and Y β UT provides correction for the single polarized process at q ⊥ ∼ Q. The antisymmetric tensor ǫ αβ ⊥ is defined as ǫ αβµν P µ π P ν p /P π · P p , and S ⊥ is the transverse polarization vector of the proton target.
The structure function W UT (Q, b) can be written in terms of the unpolarized distribution function of pion and Sivers function of proton as withf ⊥α(DY) 1T q/p (x p , b; µ, ζ F ) given in Eq. (33). Similar to the unpolarized case, the scheme-dependent hard factors can be absorbed into the C-coefficients, leading to [64,84] The spin-dependent differential cross section in Eq. (56) thus has the form Combing Eqs. (55), (54), (59), one can get the Sivers asymmetry in the Drell-Yan process with a π beam colliding on a transversely polarized proton target.

C. The cos 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 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 [104]. The coefficients λ, µ, ν in Eq. (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 Q T ∼ 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 p N Drell-Yan processes in both the fixed target mode [105,106] 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][51][52][53][54][55][56][107][108][109][110][111]. 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.
The cos 2φ asymmetry coefficient ν BM contributed by the Boer-Mulders function can be written as where the notation has been adopted to express the convolution of transverse momenta. q T , k T and p T are the transverse momenta of the lepton pair, quark and antiquark in the initial hadrons, respectively.ĥ is a unit vector defined asĥ = qT |qT | = qT qT . One can perform the Fourier transformation from the q T space to b space on the delta function in the notation of Eq. (62) to obtain the denominator in Eq. (61) as where the unpolarized distribution function in the b space is given in Eqs. (25) and (26). Similar to the treatment of the denominator, using the expression of the Boer-Mulders function in Eq. (34) the numerator can be obtained as Different from the previous two cases, the hard coefficients C BM i and H BM for the Boer-Mulders function have not been calculated up to next-to-leading order (NLO), which still remained in leading order (LO) as C q←i = δ qi δ(1 − x) and H = 1.

IV. NUMERICAL ESTIMATE FOR THE TMD DISTRIBUTIONS
Based on the TMD evolution formalism for the distributions set up in Sec. II, 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 has the been studied numerically in Refs. [23,66,72].

A. The unpolarized TMD distribution of the pion meson
In Ref. [79], the authors applied Eq. (26) and the extracted parameters g π 1 and g π 2 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 parametrization [112] was chosen. The results are plotted in Fig. 2, with the left and right panels showing the subtracted distribution in the b−space and k ⊥ -space, for fixed x π = 0.1, at three different energy scales: Q 2 = 2.4 GeV 2 (dotted line), 10 GeV 2 (solid line), 1000 GeV 2 (dashed line). From the b-dependent plots, one can see that at the highest energy scale Q 2 = 1000 GeV 2 , the peak of the curve is in the low b region where b < b max , since in this case the perturbative part of the Sudakov form factor dominates. However, at lower energy scales, e.g., Q 2 = 10 GeV 2 and Q 2 = 2.4 GeV 2 , the peak of the b−dependent distribution function moves towards the higher b region, indicating that the nonperturbative part of the TMD evolution becomes important at lower energy scales. For the distribution in transverse momentum k ⊥ space, at the high energy scale the distribution has a tail falling off slowly at large k ⊥ , while at lower energy scales the distribution function falls off rapidly with increasing k ⊥ . 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 Ref. [72].

B. 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,64,84,113,114], 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 Ref. [39], the Qiu-Sterman function was assumed to be proportional to f 1 , namely, it follows the same evolution kernel as that for f 1 . A different choice was adopted in Ref. [114], where the homogenous terms of the exact evolution kernel for the Qiu-Sterman function [84,[115][116][117][118][119][120][121][122][123] was included to deal with the scale dependence of Qiu-Sterman function: with P f1 qq the evolution kernel of the unpolarized PDF P f1 qq = 4 3 To solve the QCD evolution numerically, we resort to the QCD evolution package HOPPET [124] and we custom the code to include the splitting function in Eq. (65). For comparison we plot the TMD evolution of the Sivers function for proton in the b space and the k ⊥ space using the above mentioned two approaches [114]. In this estimate, the next leading order C-coefficients ∆C T q←i was adopted from Ref. [64,84] and the nonperturbative Sudakov form factor for the Sivers function of proton was adopted as the form in Eq. (17). The Sivers functions are presented at three different energy scales: Q 2 = 2.4GeV 2 , Q 2 = 10GeV 2 and Q 2 = 100GeV 2 . Similar to the result for f 1 , 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.

C. 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 Eqs. (34) and (36) in Ref. [113], in which the collinear twist-3 correlation function T (σ) q,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 [125]. For the scale evolution of T (σ) q,F , the exact evolution effect has been studied in Ref. [115]. For our purpose, we only consider the homogenous term in the evolution kernel with ∆ T P qq (x) = C F 2z (1−z)+ + 3 2 δ(1 − x) being the evolution kernel for the transversity distribution function h 1 (x). We customize the original code of QCDNUM [126] to include the approximate kernel in Eq. (67). Different to the unpolarized distributions f 1 , transversity distribution h 1 and Sivers function, the C-coefficients for Boer-Mulders function have not been calculated up to the next leading order. Thus, the C-coefficients in the numerical procedure stayed at leading order being delta functions. As for the nonperturbative part of the Sudakov form factor associated with Boer-Mulders function, the information still remains unknown. The assumption that S NP for Boer-Mulders function is same as that for f 1 can be a practical way to access the information of TMD evolution for Boer-Mulders function.
We plot the b dependent and k T -dependent Boer-Mulders function at x = 0.1 in the left and right panels of Fig. 4, respectively. In calculatingh ⊥ 1,q/π (x, b; Q) in Fig. 4, we have rewritten the Boer-Mulders function in b space as The three curves in each panel correspond to three different energy scales: Q 2 = 0.25GeV 2 (solid lines), Q 2 = 10GeV 2 (dashed lines), Q 2 = 1000GeV 2 (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. As a conclusion, we find that the tendency of the distributions is similar: the distribution is dominated by perturbative region in the b space at large Q 2 , while at lower Q 2 the distribution shifts to the large b region, indicating that the nonperturbative effects of TMD evolution become important. For the distributions in the k ⊥ space, as the value of Q 2 increasing, the distributions become wider with a perturbative tail, while at low values of Q 2 , the distributions resemble gaussian-type parametrization. However, the widths of the transverse momentum differ among different distributions.

V. NUMERICAL ESTIMATE FOR THE PHYSICAL OBSERVABLES IN π-N DRELL-YAN PROCESS
Based on the general TMD factorization framework provided in Sec. III, 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][37][38][39][40][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][43][44][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 NH 3 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 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.
A. 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 Eq. (54). Applying the extracted nonperturbative Sudakov form factor for pion meson in Ref. [79] and the extracted S NP for nucleon in Ref. [72], 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 Fig. 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 supposed to hold. The comparison somehow confirms the validity of extraction of the nonperturbative Sudakov form factor for the unpolarized distribution f 1π 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 ⊥ values, the numerical estimate in Ref. [79] cannot describe the data, indicating that the perturbative correction from the Y UU term may play an important role in the region q ⊥ ∼ Q. The further study on Y term is needed to provide a complete picture The transverse spectrum of lepton pair production in the unpolarized pion-nucleon Drell-Yan process, with an NH3 target at COMPASS. The dashed line represents the theoretical calculation in Ref. [79]. The solid line shows the experimental measurement at COMPASS [11]. Figure from Ref. [79].
of the q ⊥ distribution of lepton pairs from πN Drell-Yan in the whole q ⊥ range.

B. The Sivers asymmetry
In Ref. [39], the authors adopted the Gaussian form of the nonperturbative Sudakov form factor S NP in Eq. (17) and the leading order C coefficients to perform a global fit on the Sivers asymmetry in SIDIS with the experimental data from HERMES [127], COMPASS [128,129], and Jefferson Lab (JLab) [130]. 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,131], 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 Ref. [39].
The TMD evolution effect of the Sivers asymmetry in SIDIS and pp Drell-Yan at low transverse momentum has also been studied in Ref. [84], in which a framework was build to match SIDIS and Drell-Yan and cover the TMD physics with Q 2 from several GeV 2 to 10 4 GeV 2 (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 Q 2 range from 2 to 100 GeV 2 . 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, Ref. [84] 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 Ref. [11]) within the error bar.
With the numerical results of the TMD distributions in Eq. (57), the Sivers asymmetry A Siv UT as function of x p , x π , x F and q ⊥ in π − p ↑ → µ + µ − + X in the kinematics of COMPASS Collaboration was calculated in Ref. [114], as shown in Fig. 6. The magnitude of the asymmetry is around 0.05 ÷ 0.10, which is consistent with the COMPASS measurement (full squares in Fig. 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.  [11]. Figure from Ref. [114].
C. The cos 2φ azimuthal asymmetry Using Eq. (64), the cos 2φ azimuthal asymmetry contributed by the double Boer-Mulders effect in the πN Drell-Yan process was analyzed in Ref. [113], 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 parametrization in Ref. [91] at the initial energy Q 2 0 = 1GeV 2 . As mentioned in Sec. IV C, the Boer-Mulders function of the pion is adopted from the model calculation in Ref. [125]. Here we plot the estimated asymmetry ν BM as function of x p , x π , x F and q T in the kinematical region of COMPASS in Fig. 7. The bands correspond to the uncertainty of the parametrization of the Boer-Mulders function of the proton [91]. 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 x p , x π , x F show slight dependence on the variables, while the q T dependent asymmetry shows increasing tendency along with the increasing q ⊥ in the small q ⊥ range where the TMD formalism is valid. The result in Figs. 7 indicates that, precise measurements on the Boer-Mulders asymmetry ν BM as functions of x p , x π , x F and q T 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.  [91]. Figures from Ref. [113].

VI. 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: • 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.
• 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 th Qiu-Sterman function can affect the shape of the asymmetry and should be considered in the future exaction of the Sivers function.
• Sizable cos 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 of S NP based on more precise experimental data is needed, such as including the flavor dependence and hadron dependence on the functional form for S NP . 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 numerical estimate could be done to test how negligible the term is. In addition, the TMD factorization is suitable to describe the small transverse momentum physics, while the collinear factorization 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 [60,132].