Pion-induced Drell-Yan processes within TMD factorization

We extract the pion transverse momentum dependent (TMD) parton distribu- tion by fitting the pion-induced Drell-Yan process within the framework of TMD factoriza- tion. The analysis is done at the next-to-next-to-leading order (NNLO) with proton TMD distribution and non-perturbative TMD evolution extracted earlier in the global fit. We observe the significant difference in the normalization of transverse momentum differential cross-section measured by E615 experiment and the theory prediction.


Introduction
Transverse momentum dependent (TMD) factorization theorem allows for a systematic study of partons transverse motions. Being equipped by the next-to-next-to-leading order (NNLO) evolution and matching, TMD factorization establishes an accurate framework for extractions of TMD distributions and production of trustful predictions for TMD crosssections. It has been recently demonstrated in ref. [1] by the global analysis of Drell-Yan process. In this work, I extend the analysis of ref. [1] by considering the pion-induced Drell-Yan process and extracting the pion unpolarized TMD parton distribution function (TMDPDF). Apart of the pure scientific interest this study is stimulated by the upcoming measurement of the pion-induced Drell-Yan process at COMPASS facility [2]. Formulated in [3,4], TMD factorization theorem has been proven at all orders of perturbation theory [5][6][7][8]. Within the modern construct, the TMD distributions are generic non-perturbative functions that obey the double-scale evolution [8,9] and match collinear distributions at small-b limit [6,7,[10][11][12][13][14][15]. The matching to the perturbative limit almost guaranties the agreement with the high-energy data and the collinear factorization. Simultaneously, it greatly constraints the value of TMD distributions in the numerically dominant part of cross-section formula. As a result, the TMD factorized cross-section has a great predictive power even at low-energies, where influence of non-perturbative corrections is larger. Let me note, that NNLO perturbative input is important to describe precise modern data [16].
The TMD factorized cross-section contains three non-perturbative functions. There are two TMD distributions and the non-perturbative evolution kernel. It is practically difficult to decorrelate these functions. In ref. [1] the large data-set has been considered with significant difference in energy (from 4 to 150 GeV), which allowed to reduce the correlation between non-perturbative evolution and TMDPDFs. In this work, the situation is simpler since the non-perturbative evolution and the proton TMDPDF are taken from [1]. Therefore, the extraction pion TMDPDF is direct, and can be considered as a part of global fit of Drell-Yan data.

JHEP10(2019)090
The numerical part of the work has been done by artemide package [17,18]. Artemide is the library of fotran modules related to different aspects of TMD factorization, from the small-b matching to the computation of the cross-section (including bin-integration and fiducial cuts, if required). The PDF sets are provided via the LHAPDF interface [19]. The artemide repository also includes sets for TMD distributions (and their evolution) together with distributions of replicas. The results of the current extraction are added to the repository, as Vpion19 set.
The pion-induced Drell-Yan process does not attract too much attention. For review of recent development see ref. [20]. Perhaps, the main reason is the low quality of the data. The last measurement has been done in the end of 80's at E615 at FermiLab [21]. In this work I have observed a systematic disagreement between E615-data and theory predictions in the normalization value. Currently, it is not possible to decide: is the disagreement a problem of theory or of the data. The similar problems have been observed recently in [22] in comparison of collinear factorization to low-energy Drell-Yan process. Hopefully, COMPASS results will resolve this issue.
The paper consists of two sections. The section 2, I briefly review the TMD factorization framework, with the emphasis on difference between this work and ref. [1], which consists in introduction of exact matching for zeta-line at large-b. The section 3 is devoted to the comparison of the theoretical prediction to the data and to the extraction of the pion TMDPDF. The significant part of section 3 is the discussion of the problem with the normalization for E615 measurement, and its possible origins. In the appendix A the derivation of the exact expression for the special null evolution line that was used for low-energy TMD evolution is presented.

Theoretical framework
The derivation of the cross-section for Drell-Yan process in the TMD factorization has been a subject of many studies, see e.g. refs. [5][6][7]16]. In this section I present only the main formulas used in this analysis. The theory framework coincides with refs. [1,16]. The particular points specific for discussed case are presented in details.
Cross-section within TMD factorization. The cross-section for (2.1) where q is the momentum of the photon, with the virtuality q 2 = Q 2 , and the transverse component q T . Variable x F is the Feynman x related to Bjorken x's and τ in a usual manner,

JHEP10(2019)090
The common factor for the cross-section is and the hard coefficient function is where sum over q runs though quarks and anti-quarks, L = ln(Q 2 /µ 2 ), C F = 4/3. The NNLO term of the hard coefficient function used in current evaluation can be found in [23].
The functions F f ←h (x, b; µ, ζ) in (2.1) are TMDPDF for parton f in the hadron h evaluated at the scale (µ, ζ).
Selection of scales and TMD evolution. The scales µ 2 , ζ 1 and ζ 2 are of order of Q 2 .
In ζ-prescription, TMDPDFs are defined at the line ζ = ζ(µ, b), which is a nullevolution line in the plane (µ, ζ). The optimal TMD distribution used in this work, belongs to the null-evolution line that passes through the saddle point of the evolution field. This boundary condition is very important for two reasons. First, there is only one saddle point in the TMD evolution field, and thus the special null-evolution line is unique. Second, the special null-evolution line is the only null-evolution line, which has finite ζ at all values of µ (with µ bigger than Λ QCD ). It follows from the definition of the saddle point, and guaranties the finiteness of perturbative series at each order. The optimal distribution is denoted as F f →h (x, b) (without scale arguments), what emphasizes its uniqueness and "naive" scale-invariance. The relation between the optimal TMD distribution and TMD distribution at the scale (µ, ζ) = (Q, Q 2 ) is where D is the rapidity anomalous dimension. The derivation of this simple expression and proof of its equivalence to the standard Sudakov exponent is given in [9]. The subscript NP on the rapidity anomalous dimension D NP and the special null-evolution line ζ NP stresses the presence of non-perturbative corrections in both objects.
Expression for TMDPDF. There are two places where non-perturbative physics enters TMD factorized cross-section. The first is the TMDPDFs F (x, b) that describes transverse motion of confined quarks in a hadron. The second is the rapidity anomalous dimension D NP (µ, b) that describes the long-range correlation of gluons in QCD vacuum. So, nonperturbative structure of these objects are related to different aspects of QCD dynamics JHEP10(2019)090 and are completely independent. At small values of b both F (x, b) and D NP (µ, b) could be calculated by means of the operator product expansion, see e.g. [11,12,14,15,24]. At large-b the values of these function should be calculated in non-perturbative models, as e.g. in refs. [25][26][27], or extracted from the data, as e.g. in refs. [1,16,[28][29][30].
The convenient ansatz that merges perturbative and non-perturbative part of TMD-PDF functions is where C is the perturbative coefficient function calculated at NNLO in [12,13], f 1 is unpolarized collinear PDF, and f NP is a non-perturbative modification function. The function f 1 must turn to 1 at b → 0. The selection of an ansatz for f NP is a delicate process, since it is the main source of biases, for a more detailed discussion see section 2 in ref. [1].
In the present work, the proton TMDPDF is taken from ref. [1], where it was extracted from the global fit of high-energy (Tevatron and LHC) and low-energy (FermiLab and PHENIX) Drell-Yan measurements. The analyzed measurements (E537 [31], E615 [21] and NA3 [39]) were made on a tungsten (E537, E615) and platinum (NA3) targets (Z = 74, A = 184 and Z = 78, A = 195). Therefore, the proton TMDPDF from [1] requires a modification to simulate the nuclear environment. It is done by the rotation of the iso-spin components only. For example, for u-quark the nuclear TMDPDF is and similar for d,ū andd distributions. The values of pion TMDPDF are fit to the data, as discussed in the following. The collinear pion PDF is taken JAM18pionPDFnlo-set [32]. The function f NP is taken similar to those used for proton in [1]. Taking into account the fact that typical values of x in the pion-induced Drell-Yan process are very high, the terms relevant for low-x values were dropped. The resulting function depends on three parameters and reads The parameters a 1,2,3 are to be fit to the data. Generally speaking, the non-perturbative function f NP depends on the flavor of parton. This dependence is ignored here, since the quality of analyzed data does not allow a flavor separation.
Expression for D NP . The non-perturbative expression for D NP has been extracted in [1] together with proton TMDPDF. It has the following form where D pert (µ, b) is the perturbative part of rapidity anomalous dimension calculated at NNLO in [5,33], and at N 3 LO in [34,35]. The function d NP (b) is the non-perturbative (2.10) The non-perturbative function d NP is The parameters B NP and c 0 are fit in [1], and reads B NP = 2.29 ± 0.43, c 0 = 0.022 ± 0.009. The only difference in the theory implementation between this work and ref. [1] is the expression for ζ NP . In ref. [1], ζ NP was modeled ζ NP (b) = ζ pert (b * (b)). It corresponds to the special null-evolution line derived for D NP = D pert (µ, b * (b)), ignoring the d NP -contribution. This choice is almost perfect for d NP D pert , but this model deviates from the exact ζ NP significantly for larger d NP (that happens at b > B NP ). The exact ζ NP is ζ NP determined by its differential equation (A.3). The deviation of ζ NP from its exact values at large-b could be seen as a part of non-perturbative model. However, it adds an undesired correlation between TMDPDFs and D NP at large-b. The only reason to use the model values for ζ NP in ref. [1] was the absence of a way to find exact ζ NP at large b, where the saddle point runs to µ < Λ QCD region. This problem has been solved recently by using the D NP as an independent variable that accumulates all non-perturbative information and b-dependence. In this case, expression for exact ζ NP can be found order-by-order in a s which is the only parameter. The details of calculation and explicit expression for ζ exact are given in appendix A.
One of the requirement of the small-b matching procedure (2.6) in the ζ-prescription is that the values of ζ-line should exactly match its pure perturbative expression at b → 0. Otherwise the exact cancellation of divergent ln(b 2 ) in the matching coefficient C(x, b, µ) in (2.6) does not take place [9]. In order to facilitate the cancellation, the following form for ζ NP has been used This form ζ NP exactly matches ζ pert at b B NP and smoothly turns to exact value.
Perturbative orders. Let me summarize the orders of perturbation theory are used in this work: up to a 2 s -terms inclusively) [23].
• The perturbative part of rapidity anomalous dimension D perp (µ, b) in (2.9) is taken at NNLO (i.e. up to a 2 s -terms inclusively) [33], in the resummed form [9,36]. Table 1. The synopsis on the data used in the work. N pt is the number of points in the data set after/before the application of TMD factorization cut. Typical statistical error is estimated from the first 3 points for each (Q, x F )-bin, and presented only for demonstration purposes. The data for NA3 is available only as a figure (figures 1 and 2 in ref. [39]).
• To evaluate expressions for last three points one needs cusp anomalous dimension and γ V anomalous dimension up to a 3 s -terms and a 2 s -terms, receptively. They could be found in [37,38].
Thus, the computation is done at complete NNLO perturbative accuracy.

Comparison to the data
Review of available data. There are three available measurements of transverse momentum cross-section for pion-induced Drell-Yan process. They were performed by NA3 [39], E537 [31] and E615 [21] experiments. The measurement by NA3 is presented only by a plot in ref. [39], and the exact values of data-points and their error-bars are not available. Therefore, only the visual comparison with NA3 is possible ( figure 6). The data tables for E537 [31] and E615 [21] can be can be found in [40].
Both experiments E537 and E615 have been performed in the same environment at different energies of the pion beam, P beam = 125 GeV for E537 and P beam = 252 GeV for E615, which corresponds to s = 235.4 GeV 2 and s = 473.6 GeV 2 , respectively. The data for both experiments are provided in two alternative binning: differential in x F , or differential in Q. In table 1, the summary of kinematics for each data set is shown. Let me mention that both of measurements are made at high values of x 1,2 . In particular, the lowest accessible value of x π is 0.26 (for E537) and 0.18 (for E615).

JHEP10(2019)090
Definition of χ 2 distribution. To estimate the theory-to-data agreement, I have used the χ 2 -test function contracted as usual where m i is the central value of i'th data-point, t i is the theory prediction for this data-point and i, j run through all points in the set. Both experiments provide an uncorrelated error for each point σ stat and a systematic error. The later is mainly generated by the luminosity uncertainty, and thus can be considered as a correlated error σ corr . The covariance matrix is be build according to general rules: As it is discussed below, E615 data-set has a problem with the general normalization (it could be also a theory problem). Due to it, the value of χ 2 calculated with covariance matrix (3.2) is extremely high, despite the errors of measurements are relatively large. It happens because the correlated part of χ 2 overweights the uncorrelated part by an order of magnitude. So, the fit procedure becomes impossible. To stabilize the values of χ 2 , I have split the data-set of E615 to subsets with the same values of (Q, x F ). Consequently, the correlated error has been adjusted to (Q, x F )-bin independently. In other words, the elements of covariance matrix that mixes different (Q, x F )-bins are set to zero. The possible sources of this problem are discussed below.
Selection of the data for the fit. The TMD factorization formula is derived in assumption that q T /Q is small. Practically, it is realized by considering the points with q T < δ · Q, where δ ≈ 0.25 as it has been derived in [16] from the global analysis of Drell-Yan measurements. For x F -differential measurements that have wide Q-bins, the center of Q-bin is used, what corresponds to q T 2.2 GeV.
There are 4 data-sets listed in table 1. Data-sets belonging to the same experiment could not be added to a single χ 2 , since it would imply a double counting of a measurement. Doubtless, it is preferable to consider x F -differential bins, since the Q-dependence is dictated by the evolution that is fixed from other data. Therefore, for the fit of nonperturbative parameters only the E615 differential in x F data-set has been used. Furthermore, the bins x F ∈ (0.8, 0.9) and x F ∈ (0.9, 1.0) have been excluded, because there x π ∼ 1 and thus, the threshold resummation must be applied. The resulting set has 80 points.
Dependence on collinear PDFs. According to (2.6), the values of TMDPDF depend on collinear PDF. The dependence is partially compensated by the non-perturbative parameters of TMDPDF, that are fit separately for each PDF set. Nonetheless, the values of TMDPDF based on different PDF sets could significantly vary. The choice of PDF set also affects the non-perturbative TMD evolution, although in a lesser amount. The original BSV19 extraction uses NNPDF3.1 set of collinear PDFs [41]. Additionally, the extraction of TMDPDFs and D NP based on different collinear PDFs were performed (the analysis of these results will be presented elsewhere), and they are available at [17,18]. In the present study, I have compared the predictions generated with proton TMD-PDFs (and D NP ) based on different collinear PDF, and found results alike. Particularly, χ 2 -minimization with proton TMDPDFs based on MMHT14 (nnlo) [42], NNPDF3.1 (nnlo) [41] and HERA20PDF (nnlo) [43] gives χ 2 /N p = 1.45, 1.70 and 1.44, correspondingly. Taking into account, that HERA20PDF set also shows better global χ 2 on the data-set from ref. [1], in the following the proton TMDPDF and non-perturbative part of TMD evolution is based on HERA20PDF are used. This set BSV19.HERA20PDF can be downloaded from artemide repository [17,18]. For pion collinear PDF JAM18pionPDF-set has been used [32].
Results of the fit. The minimization procedure for χ 2 -test yields the following values of non-perturbative parameters a 1 = 0.17 ± 0.11 ± 0.03, a 2 = 0.48 ± 0.34 ± 0.06, a 3 = 2.15 ± 3.25 ± 0.32. (3. 3) The first error-band is due to the uncertainty of data-points. It is estimated by the replica method, as in ref. [44], by minimization of χ 2 on 100 replicas of pseudodata. The second error is due to uncertainty in the proton TMDPDF and TMD evolution. It is estimated by the minimization of χ 2 on 100 of replicas of input distributions. Parameters a 1,2,3 are restricted to positive values. So, large error-bands in (3.3) are the result of very asymmetric distribution of parameters. Large error bands on parameters does not implies a significant point-by-point uncertainty for f NP , since all parameters are correlated. For example, at b ∼ 0.5 GeV −1 the uncertainty in f NP is 2-3%. However, this band is definitely biased by the ansatz (2.8). The plot for f NP is shown in figure 1 (left). The actual values of TMDPDF in b−space and k T -space (that is obtained by Fourier transformation) are shown in figure 1 (center, right). The pion TMDPDF obtained in this work together with distribution of 100 replicas is available in the artemide-repository [17,18] as Vpion19 TMDPDF set (for π − -meson).
The final values of χ 2 is χ 2 /N p = 1.44 (N p = 80). It can be compared with the result of fit in ref. [45] χ 2 /N p = 1.64, where almost the same data were used. The main contribution to the value of χ 2 comes from the systematic disagreement in the normalization between the data and the theory. In figure 2, 3, 4, 5 the comparison of the data to the theory prediction is shown together with the values of χ 2 /N p for a given subset of data-points. In  figure 6 the visual comparison of the theory to NA3 is shown. The plots for Q-differential bins made for the range of q T larger than it is allowed by the TMD factorization (the boundary q T 0.25Q is shown by the vertical dashed line). It is interesting to observe that the TMD factorization formula works unexpectedly well outside of this region.
Normalization issue. The main problem of presented analysis is the significant difference in the common value (normalization) between the theory prediction and E615 measurement. For a deeper understanding of this issue, it is instructive to perform the decomposition of χ 2 values as where χ 2 D (χ 2 λ ) represents the uncorrelated (correlated) part of χ 2 . Loosely speaking, the value of χ 2 D (χ 2 λ ) demonstrates the agreement in the shape (normalization) between the theory and the data. The decomposition (3.4) is done with the help of nuisance parameters [44,46]. As a by-product, this method allows determining the value of so-called   "systematics shifts" d i that are the deviation between the theory and the data due to the normalization only. The results of the nuisance-parameters-decomposition, as well as, average values of d i are presented in figure 2, 3, 4, 5 for each bin for E615 and common for E537. The decomposition of χ 2 for the selected data is The value χ 2 λ /N p = 0.77 is huge, accounting 16% systematic uncertainty. Indeed, figures 3 and 5 clearly demonstrates that the theory prediction is systematically below the data. For the first bins (the lowest x F and Q) the difference is practically factor 2. The comparison to E537 ( figure 2 and 4 does not show such a significant problem, but the quality of E537 measurement is much worse. The visual comparison to NA3 measurement (figure 6) also does not show any normalization problem. Neglecting the normalization part of the χ 2 the agreement between the data and the theory is almost perfect, which is also clear from comparison of dashed lines to data-points in figure 2, 3, 4, 5.  The analogous problem with the description of the transverse momentum spectrum for the Drell-Yan process has been recently discussed in ref. [22]. The authors of ref. [22] have observed that the data-points measured in the fixed-target experiments are significantly (2-3 times) above the theory expectations. The data analyzed in ref. [22] belong to the same kinematic domain as the data discussed here. The comparison has been done in the regime q T ∼ Q where the collinear factorization is well established. The authors have tested several ways to improve the theory predictions (threshold resummation, k T -smearing) but were not able to resolve the problem. In the TMD regime the same effect has been observed in [1] (for the same experiments that are considered in [22]), namely, about 40% deficit in the normalization that decreases with the increase of energy (see table 3 in [1]). Note, both analyses [22] and [1] have not a problem with the description of PHENIX data [47] that have a similar range of Q but measured in the collider regime. A similar problem was also observed in semi-inclusive deep-inelastic scattering (SIDIS) [48].
Previously, E615 measurement have been analyzed in the framework of TMD factorization in refs. [45] and [49]. In these articles, authors do not observe any problems with  the normalization. However, in both cases, functions used to fit the non-perturbatibatve parts include parameters that significantly influence the normalization. Therefore, it is possible that the normalization issue discussed here, was absorbed into model parameters in refs. [45,49]. The present situation could appear because of the problem with the theory. Let me list possible flaws of the current consideration • Nuclear effects. The nuclear effects are, for sure, presented in current measurements and goes beyond iso-spin modification (2.7). Generally, an extra factor R A i (x) for PDF should be added, see e.g. [50]. Typically, at x ∈ (0.1, 0.9) this factor provides ∼ 10% modification [50], which cannot compensate the gap between the theory and the data. Moreover, this effect should be much smaller for x-integrated bins (figure 5) due to the oscillation of R A i (x) between anti-shadowing and EMC regimes.  Figure 6. Comparison of the theory prediction (solid line) to NA3 measurement. The theory prediction is plot on the top of figures 1 and 2 by ref. [39]. The vertical dashed line shows the estimation of the boundary for TMD factorization approach.
• Effects of PDF. The collinear PDFs are poorly known at large-x, and values of PDF significantly differ between different sets. In particular, the difference between PDF values at large-x completely resolves the normalization issue (of the order of 5%) with LHCb Z-boson spectrum in [1]. I have checked that in the present kinematics the usage of different PDF sets could produce up to 20% difference at a point. Even so, it mainly affects the shape of the cross-section, whereas the normalization is affected only by 2-3%. Note, that the pion PDF were extracted mainly from the integrated over q T measurement by E615 [32], and in the present analysis TMDPDF accurately (at NNLO) matches collinear PDF.
• Threshold contributions. The large-x effects must be incorporated into the matching coefficient in (2.6). To my opinion, the ignorance of threshold effects leads to the disagreement in the shape of cross-section for bins with x F > 0.7 ( figure 2, 3). However, the effect of threshold resummation should be negligible at x ∼ 0.2 and Q ∼ 4-5 GeV, where the most significant deviation takes place. Also in ref. [22] a more accurate analysis has been performed, and it has been shown that the threshold resummation does not solve the problem • Power corrections. The TMD factorization theorem violates QED Ward identities and Lorentz invariance (it is typical for factorization theorems with several scales, see e.g., discussion in [51]). To restore it, one needs to account power corrections, which could be large. Nowadays, there are no systematic studies of power corrections to TMD factorization, and their size is unknown. Nonetheless, these corrections must vanish at q T /Q → 0, and so, their presence would be indicated in the deformation of the shape of cross-section, what is not observed.

JHEP10(2019)090
• Wrong shape for non-perturbative corrections. It could happen that the suggested ansatz for non-perturbative parts of TMD evolution (2.9) and TMDPDF (2.6) is essentially wrong, and confines the cross-section in improper domain. However, it looks very implausible because it agrees with known theory constraints, and nicely describe the proton-proton measurements [1].
• Resonance effects. The most problematic bins are the lower-Q bins. It could imply that the observed deficit in the normalization is produced by the interference of γ * with J/ψ, ψ resonances and their excitations that are located in the region Q ∼ 3-4 GeV. However, the post-resonance contamination typically looks exactly opposite, as an excess of the theory over the data.
In total, it is hard to imagine that any of these points (except resonance contamination) could change the value of cross-section normalization more than 5-10%. Unless the TMD factorization formula has a deep and systematic problem. Thus, I should conclude that probably the differential in q T data by E615 have an incorrect normalization. There are some details that further point to this possibility. First, there is a very good agreement in the shape of cross-sections. Second, the normalization issue is greater at smaller-Q and practically disappears at Q ∼ 9 GeV (the same with x Fdifferential bins since x F ∼ Q/ √ s). It could indicate the bad estimation of the background in the close-to-resonance region by E615 collaboration. Additionally, the traces of abnormal behavior in x F (for q T -spectrum) are already seen in the publication of E615 [21]. It was observed that q T -spectrum after subtraction of normalization has an extreme dependence on x F (see section V.B, and appendix A, in ref. [21]), which could not be explained within the perturbative QCD. Finally, the comparison to E537 and NA3 experiments has not a problem with normalization, although data-quality is significantly worse.

Conclusion
In the present work, the pion-induced Drell-Yan process has been studied, with the main aim to extract the values of pion unpolarized transverse momentum dependent parton distribution function (TMDPDF). The analysis is made in the TMD factorization framework with ζ-prescription [9] and compete next-to-next-to-leading (NNLO) perturbative input. To extract the values of pion TMDPDF, the measurements of E615 experiment have been used. I have used the differential in x F data for better sensitivity to x-dependence of TMD-PDF. The measurement of E615 differential in Q and measurements by E537 and NA3 were used for the cross-check of the fit. The resulting pion TMDPDFs are available as a part of artemide (model Vpion19) -the program package for TMD phenomenology [17,18].
During the fit procedure, I have faced the problem of systematic disagreement in the normalization between data and the theory. The measurements with low-Q and, correspondingly low-x F , are significantly higher (up to two times for Q ∼ 4-5 GeV) than the prediction. Simultaneously the shape of cross-sections is in an excellent agreement. The size of discrepancy in the normalization decreases with the increase of the Q. In the last part of section 3, I provide a discussion on possible sources of normalization disagreement and conclude that I do not see any possibility to obtain such a significant factor within the modern TMD factorization framework. There is a possibility that the observed normalization problem has an experimental origin. The comparison of the theory with E537 and NA3 has not such a problem, but the both experiments have much worse precision, and could not seriously compete with E615. A similar problem has been recently observed in the TMD spectrum of proton-nucleus Drell-Yan process in [22]. Within the nearest future, the COMPASS collaboration will repeat the analysis of the pion-induced Drell-Yan process in the similar kinematics regime. The announcement of this measurement is presented in [52]. In figure 7 (left) the comparison of the preliminary COMPASS data to the prediction made with Vpion19 is shown. Hopefully, the COMPASS measurement will resolve the problem with the normalization of E615 experiment.
A particularly engaging point to study pion TMDPDF is its comparison to proton TMDPDF since the confined motion of partons in mesons and baryons could be fundamentally different. However, any principal difference is not observed (at moderate x), see figure 7 (right). At high-x distributions looks different, but no conclusion can be done since high-x region is not well controlled both experimentally and theoretically. Definitely, the future measurements of TMD cross-section for pion-induced Drell-Yan process will shed light to this side of parton dynamics.

JHEP10(2019)090
expression for the special null-evolution line that exactly incorporates non-perturbative corrections.
A null-evolution line is defined as an equipotential line for the 2-dimensional field of anomalous dimensions E = (γ F (µ, ζ)/2, −D(µ, b)) in the plane (µ, ζ). The anomalous dimension γ F is the ultraviolet anomalous dimension of TMD operator. It has the following form where Γ cusp is the cusp-anomalous dimension, and γ V anomalous dimension of the vector form-factor. The rapidity anomalous dimension D(µ, b) is generally non-perturbative function, which can be computed perturbatively only at small-b, see e.g. [33,34] for NNLO and N 3 LO computations. It satisfies the renormalization group equation Due to this expression the field E is conservative. Parameterizing an equipotential line as (µ, ζ(µ, b)), one finds the following equation for ζ(µ, b) The special null-evolution line is the line that passes thorough the saddle point (µ 0 , ζ 0 ) of the evolution field. The saddle point is defined as Such boundary condition are very important for two reasons. First, there is only one saddle point in the evolution field, and thus, the special null-evolution line is unique. Second, the special null-evolution line is the only null-evolution line, which has finite ζ at all values of µ (bigger than Λ QCD ). It follows from the definition of the saddle point, and guaranties the finiteness of perturbative series order-by-order.
The field E, and consequently the equipotential line ζ(µ) and the position of the saddle point (µ 0 , ζ 0 ), depends on b, which is treated as a free parameter. It causes certain problems in the implementation of the ζ-prescription. The lesser problem is that additional numerical computations are required to determine the position of saddle-point and the values of the line for different non-perturbative models of D. The greater problem is that at larger b the value of µ 0 decreases and at some large value of b (typically b ∼ 3 GeV −1 ) µ 0 is smaller than Λ QCD . Due to this behavior, it is impossible to determine the special null-evolution line at large-b numerically. Note, that nonetheless the special null-evolution line is still uniquely defined by the continuation from smaller values of b. In ref. [1] the value of the special null-evolution line has been approximated by perturbative expression with b = f (b), which exactly matches true values at b → 0, and starts to significantly deviate from exact values at b ∼ 3-4 GeV −1 . This deviation has been considered as a part of non-perturbative model for evolution, which somewhat undermine universality of non-perturbative TMD evolution JHEP10(2019)090 kernel, and adds correlation between non-perturbative parts of TMD evolution kernel and TMDPDFs. Recently, I have found a simple solution for the problem of determination of the special null-evolution line, which is presented here.
The main breakthrough idea is to use the non-perturbative rapidity anomalous dimension as a generalized coordinate instead of the scale µ. It could not be done entirely, since scale µ also enters QCD coupling constant in anomalous dimensions Γ cusp and γ V . For values of µ large-enough the value of a s (µ) is small, and thus the solution could be evaluated order-by-order in a s (µ). Important, that the non-perturbative dependence is exactly accounted in such approach. The equation (A.3) can be rewritten where p = 2β 0 D/Γ 0 . Let me mention that NNLO term is exponentially grow at large-D (the N 3 LO term grows even faster as e 2p ). However, it is not a problem, since i) g enters the logarithm, ii) asymptotic regime takes place at very large values of b, iii) altogether such behavior only suppresses high-b tale of the evolution exponent. The expressions (A.7)-(A.9) provide a very accurate approximation, since a s is evaluated at µ = Q and typically a s = g 2 /(4π) 2 ∼ 10 −2 . The most important is that this expression is valid at all values of b, even then saddle point is below Λ QCD .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.