Pion-induced Drell-Yan processes within TMD factorization

TL;DR

This paper extracts the pion's TMD parton distribution using NNLO analysis of Drell-Yan data within TMD factorization, highlighting normalization discrepancies with experimental measurements.

Contribution

It provides a next-to-next-to-leading order extraction of the pion TMD distribution within the TMD factorization framework, incorporating previous global fit results.

Findings

Significant normalization difference observed in E615 data

Successful NNLO extraction of pion TMD distribution

Insights into non-perturbative TMD evolution

Abstract

We extract the pion transverse momentum dependent (TMD) parton distribution by fitting the pion-induced Drell-Yan process within the framework of TMD factorization. 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.

Figures

Tables

Table 1. Table 1: The synopsis on the data used in the work. N p ​ t subscript 𝑁 𝑝 𝑡 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 ) 𝑄 subscript 𝑥 𝐹 (Q,x_{F}) -bin, and presented only for demonstration purposes. The data for NA3 is available only as a figure (figs. 1 and 2 in ref. Badier:1982zb ).

Experiment	$\sqrt{s}$[GeV]	$Q$[GeV]	$x_F$	$N_{pt}$	corr.err.	Typical stat.err.
E537 ($Q$-diff.)	15.3	in 10 bins		60/146	8%	$\sim 20\%$
E537 ($x_F$-diff.)	15.3		in 11 bins	110/165	8%	$\sim 20\%$
E615 ($Q$-diff.)	21.8	in 10(8) bins		51/155	16%	$\sim 5\%$
E615 ($x_F$-diff.)	21.8		in 10 bins	90/159	16%	$\sim 5\%$
NA3	16.8 , 19.4 22.9		$y>0$(?)	–	15%	–

Equations

$$\frac{d\sigma }{d{Q}^{2}d{x}_{F}d{q}_T^2}={\sigma }_0\sum _{f_1,f_2}{H}_{f_1{f}_2}(Q,\mu ){\int }_0^{\infty }\frac{bdb}{2}{J}_0(b{q}_T)F_{f_1\leftarrow h_1}(x_1,b;\mu ,{\zeta }_1)F_{f_2\leftarrow h_2}(x_2,b;\mu ,{\zeta }_2),$$

$$x_{1,2}=\frac{\pm x_F+\sqrt{x_F^2+4\tau }}{2},\tau =x_1{x}_2=\frac{Q^2+{\mathbf{q}}_T^2}{s}.$$

$${\sigma }_0=\frac{4\pi {\alpha }_{\text{em}}^2(Q)}{9Q^{2}s\sqrt{x_F^2+4\tau }},$$

$$H_{f_1,f_2}(Q,\mu )=\sum _q{\delta }_{q{f}_1}{\delta }_{\bar{q}{f}_2}{e}_q^2\left[1+2a_s(\mu )C_F\left(-L^2+3L-8+\frac{7{\pi }^2}{6}\right)+O(a_s^2)\right],$$

$$F_{f\leftarrow h}(x,b;Q,Q^2)={\left(\frac{Q^2}{{\zeta }_{\text{NP}}(Q,b)}\right)}^{-{\mathcal{D}}_{\text{NP}}(b,Q)}{F}_{f\leftarrow h}(x,b),$$

$$F_{f\leftarrow h}(x,b)=\sum _{f^{\prime }}{\int }_x^1\frac{dy}{y}{C}_{f\leftarrow f^{\prime }}(y,b,\mu )f_{1,f^{\prime }\leftarrow h}\left(\frac{x}{y},\mu \right)f_{\text{NP}}(x,b),$$

$$F_{u\leftarrow A}(x,b)=\frac{Z}{A}{F}_{u\leftarrow p}(x)+\frac{A-Z}{A}{F}_{u\leftarrow n}(x)=\frac{Z}{A}{F}_{u\leftarrow p}(x)+\frac{A-Z}{A}{F}_{d\leftarrow p}(x),$$

$$f_{\text{NP}}^{\pi }(x,b)=\exp \left(-\frac{(a_1+(1-x{)}^2{a}_2)b^2}{\sqrt{1+a_3{b}^2}}\right).$$

$${\mathcal{D}}_{\text{NP}}(\mu ,b)={\mathcal{D}}_{\text{pert}}(\mu ,b^{\ast }(b))+d_{\text{NP}}(b),$$

$$\displaystyle b^{*}(b)=b\Big{/}\sqrt{1+\frac{b^{2}}{B_{\text{NP}}^{2}}}.$$

$$d_{\text{NP}}(b)=c_{0}b{b}^{\ast }(b).$$

$${\zeta }_{\text{NP}}(\mu ,b)={\zeta }_{\text{pert}}(\mu ,b)e^{-{\mathbf{b}}^2/B_{\text{NP}}^2}+{\zeta }_{\text{exact}}(\mu ,b)\left(1-e^{-{\mathbf{b}}^2/B_{\text{NP}}^2}\right).$$

$${\chi }^2=\sum _{ij}(m_i-t_i)V_{ij}^{-1}(m_j-t_j),$$

$$V_{ij}={\sigma }_{\text{stat},i}^2{\delta }_{ij}+{\sigma }_{\text{corr},i}{\sigma }_{\text{corr},j}.$$

$$a_1=0.17\pm 0.11\pm 0.03,a_2=0.48\pm 0.34\pm 0.06,a_3=2.15\pm 3.25\pm \mathrm{0.32.}$$

$${\chi }^2={\chi }_D^2+{\chi }_{\lambda }^2,$$

$${\chi }^2/N_p=0.67+0.77=\mathrm{1.44.}$$

$${\gamma }_F(\mu ,\zeta )={\Gamma }_{\text{cusp}}(\mu )\ln \left(\frac{{\mu }^2}{\zeta }\right)-{\gamma }_V(\mu ),$$

$${\mu }^2\frac{d\mathcal{D}(\mu ,b)}{d{\mu }^2}=\frac{{\Gamma }_{\text{cusp}}(\mu )}{2}.$$

$${\Gamma }_{\text{cusp}}(\mu )\ln \left(\frac{{\mu }^2}{\zeta (\mu ,b)}\right)-{\gamma }_V(\mu )=2\mathcal{D}(\mu ,b)\frac{d\ln \zeta (\mu ,b)}{d\ln {\mu }^2}.$$

$$\mathcal{D}({\mu }_0,b)=0,{\gamma }_F({\mu }_0,{\zeta }_0)=0.$$

$$2\mathcal{D}\left(1+\beta (a_s)\frac{\partial g(a_s,\mathcal{D})}{\partial a_s}-\frac{{\Gamma }_{\text{cusp}}(a_s)}{2}\frac{\partial g(a_s,\mathcal{D})}{\partial \mathcal{D}}\right)-{\Gamma }_{\text{cusp}}(a_s)g(a_s,\mathcal{D})+{\gamma }_V(a_s)=0,$$

$$g(a_s,\mathcal{D})=\frac{1}{a_s}\sum _{n=0}^{\infty }{a}_s^n{g}_n(\mathcal{D}),$$

$$\beta (a_s)=\sum _{n=0}^{\infty }{a}_s^{n+2}{\beta }_n,{\Gamma }_{\text{cusp}}(a_s)=\sum _{n=0}^{\infty }{a}_s^{n+1}{\Gamma }_n,{\gamma }_V(a_s)=\sum _{n=1}^{\infty }{a}_s^n{\gamma }_n.$$

$$g_0$$

$$g_1$$

$$g_2$$

Full text

1 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 cross-sections. It has been recently demonstrated in ref.Bertone:2019nxa by the global analysis of Drell-Yan process. In this work, I extend the analysis of ref.Bertone:2019nxa 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 Gautheron:2010wva .

Formulated in Collins:1981va ; Collins:1984kg , TMD factorization theorem has been proven at all orders of perturbation theory Becher:2010tm ; Collins:2011zzd ; GarciaEchevarria:2011rb ; Vladimirov:2017ksc . Within the modern construct, the TMD distributions are generic non-perturbative functions that obey the double-scale evolution Vladimirov:2017ksc ; Scimemi:2018xaf and match collinear distributions at small-b limit Collins:2011zzd ; GarciaEchevarria:2011rb ; Becher:2011xn ; Aybat:2011zv ; Echevarria:2016scs ; Gehrmann:2014yya ; Scimemi:2019mlf ; Scimemi:2019gge . 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 higher. Let me note, that NNLO perturbative input is important to describe precise modern data Scimemi:2017etj .

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.Bertone:2019nxa 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 Bertone:2019nxa . Therefore, the extraction pion TMDPDF is direct, and can be considered as a part of global fit of Drell-Yan data.

The numerical part of the work has been done by artemide package web . 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 Buckley:2014ana . 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.Wang:2018syo . 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 Conway:1989fs . 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 Bacchetta:2019tcu 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 sec.2, I briefly review the TMD factorization framework, with the emphasis on difference between this work and ref.Bertone:2019nxa , 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 sec. 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.

2 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.Becher:2010tm ; Collins:2011zzd ; GarciaEchevarria:2011rb ; Scimemi:2017etj . In this section I present only the main formulas used in this analysis. The theory framework coincides with refs.Scimemi:2017etj ; Bertone:2019nxa . The particular points specific for discussed case are presented in details.

Cross-section within TMD factorization.

The cross-section for $h_{1}+h_{2}\to\gamma^{*}(\to ll^{\prime})+X$ is

[TABLE]

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 $\tau$ in a usual manner,

[TABLE]

The common factor for the cross-section is

[TABLE]

and the hard coefficient function is

[TABLE]

where sum over $q$ runs though quarks and anti-quarks, $L=\ln(Q^{2}/\mu^{2})$, $C_{F}=4/3$. The NNLO term of the hard coefficient function used in current evaluation can be found in Gehrmann:2010ue . The functions $F_{f\leftarrow h}(x,b;\mu,\zeta)$ in (1) are TMDPDF for parton $f$ in the hadron $h$ evaluated at the scale $(\mu,\zeta)$.

Selection of scales and TMD evolution.

The scales $\mu^{2}$, $\zeta_{1}$ and $\zeta_{2}$ are of order of $Q^{2}$. To be specific, we fix $\mu=Q$, so $L=0$ in the expression for hard coefficient function (4), and $\zeta_{1}=\zeta_{2}=Q^{2}$, so $\zeta_{1}\zeta_{2}=Q^{4}$ as it is defined within TMD factorization Collins:2011zzd ; Vladimirov:2017ksc ; GarciaEchevarria:2011rb . From the hard scale $(\mu,\zeta)=(Q,Q^{2})$ TMDPDFs are evolved to the defining scale with the help of the TMD evolution, Aybat:2011zv ; Collins:2011zzd ; Scimemi:2018xaf . The defining scale for TMDPDF is selected in accordance to $\zeta$-prescription Scimemi:2017etj ; Scimemi:2018xaf .

In $\zeta$-prescription, TMDPDFs are defined at the line $\zeta=\zeta(\mu,b)$, which is a null-evolution line in the plane $(\mu,\zeta)$. 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 $\zeta$ at all values of $\mu$ (with $\mu$ bigger than $\Lambda_{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\to 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 $(\mu,\zeta)=(Q,Q^{2})$ is

[TABLE]

where $\mathcal{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 Scimemi:2018xaf . The subscript NP on the rapidity anomalous dimension $\mathcal{D}_{\text{NP}}$ and the special null-evolution line $\zeta_{\text{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 $\mathcal{D}_{\text{NP}}(\mu,b)$ that describes the long-range correlation of gluons in QCD vacuum. So, non-perturbative structure of these objects are related to different aspects of QCD dynamics and are completely independent. At small values of $b$ both $F(x,b)$ and $\mathcal{D}_{\text{NP}}(\mu,b)$ could be calculated by means of the operator product expansion, see e.g. Aybat:2011zv ; Bacchetta:2013pqa ; Echevarria:2016scs ; Scimemi:2019mlf ; Scimemi:2019gge . At large-$b$ the values of these function should be calculated in non-perturbative models, as e.g. in refs.Schweitzer:2012hh ; Lorce:2014hxa ; Noguera:2015iia , or extracted from the data, as e.g. in refs.Su:2014wpa ; DAlesio:2014mrz ; Bacchetta:2017gcc ; Scimemi:2017etj ; Bertone:2019nxa .

The convenient ansatz that merges perturbative and non-perturbative part of TMDPDF functions is

[TABLE]

where $C$ is the perturbative coefficient function calculated at NNLO in Gehrmann:2014yya ; Echevarria:2016scs , $f_{1}$ is unpolarized collinear PDF, and $f_{\text{NP}}$ is a non-perturbative modification function. The function $f_{1}$ must turn to $1$ at $b\to 0$. The selection of an ansatz for $f_{\text{NP}}$ is a delicate process, since it is the main source of biases, for a more detailed discussion see sec.2 in ref.Bertone:2019nxa .

In the present work, the proton TMDPDF is taken from ref. Bertone:2019nxa , 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 Anassontzis:1987hk , E615 Conway:1989fs and NA3 Badier:1982zb ) 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 Bertone:2019nxa 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

[TABLE]

and similar for $d$, $\bar{u}$ and $\bar{d}$ distributions.

The values of pion TMDPDF are fit to the data, as discussed in the following. The collinear pion PDF is taken JAM18pionPDFnlo-set Barry:2018ort . The function $f_{\text{NP}}$ is taken similar to those used for proton in Bertone:2019nxa . Taking into account the 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

[TABLE]

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 $\mathcal{D}_{\text{NP}}$.

The non-perturbative expression for $\mathcal{D}_{\text{NP}}$ has been extracted in Bertone:2019nxa together with proton TMDPDF. It has the following form

[TABLE]

where $\mathcal{D}_{\text{pert}}(\mu,b)$ is the perturbative part of rapidity anomalous dimension calculated at NNLO in Becher:2010tm ; Echevarria:2015byo , and at N3LO in Vladimirov:2016dll ; Li:2016ctv . The function $d_{\text{NP}}(b)$ is the non-perturbative correction. In (9) the resummed version of $\mathcal{D}_{\text{pert}}(\mu,b)$ Echevarria:2012pw ; Scimemi:2018xaf is used. The resummed expression for $\mathcal{D}_{\text{pert}}$ contains Landau pole at large values of $b$. To avoid it, the parameter $b$ is replaced by $b^{*}(b)$ in (9),

[TABLE]

The non-perturbative function $d_{\text{NP}}$ is

[TABLE]

The parameters $B_{\text{NP}}$ and $c_{0}$ are fit in Bertone:2019nxa , and reads $B_{\text{NP}}=2.29\pm 0.43$, $c_{0}=0.022\pm 0.009$.

The only difference in the theory implementation between this work and ref.Bertone:2019nxa is the expression for $\zeta_{\text{NP}}$. In ref.Bertone:2019nxa , $\zeta_{\text{NP}}$ was modeled $\zeta_{\text{NP}}(b)=\zeta_{\text{pert}}(b^{*}(b))$. It corresponds to the special null-evolution line derived for $\mathcal{D}_{\text{NP}}=\mathcal{D}_{\text{pert}}(\mu,b^{*}(b))$, ignoring the $d_{\text{NP}}$-contribution. This choice is almost perfect for $d_{NP}\ll\mathcal{D}_{\text{pert}}$, but this model deviates from the exact $\zeta_{\text{NP}}$ significantly for larger $d_{\text{NP}}$ (that happens at $b>B_{\text{NP}}$). The exact $\zeta_{\text{NP}}$ is $\zeta_{\text{NP}}$ determined by its differential equation (20). The deviation of $\zeta_{\text{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 $\mathcal{D}_{\text{NP}}$ at large-$b$. The only reason to use the model values for $\zeta_{\text{NP}}$ in ref.Bertone:2019nxa was the absence of a way to find exact $\zeta_{\text{NP}}$ at large $b$, where the saddle point runs to $\mu<\Lambda_{QCD}$ region. This problem has been solved recently by using the $\mathcal{D}_{\text{NP}}$ as an independent variable that accumulates all non-perturbative information and $b$-dependence. In this case, expression for exact $\zeta_{\text{NP}}$ can be found order-by-order in $a_{s}$ which is the only parameter. The details of calculation and explicit expression for $\zeta_{\text{exact}}$ are given in appendix A.

One of the requirement of the small-b matching procedure (6) in the $\zeta$-prescription is that values of $\zeta$-line should exactly match its pure perturbative expression at $b\to 0$. Otherwise the exact cancellation of divergent $\ln(b^{2})$ in the matching coefficient $C(x,b,\mu)$ in (6) does not take place Scimemi:2018xaf . In order to facilitate the cancellation, the following form for $\zeta_{\text{NP}}$ has been used

[TABLE]

This form $\zeta_{\text{NP}}$ exactly matches $\zeta_{\text{pert}}$ at $b\ll B_{\text{NP}}$ and smoothly turns to exact value.

Perturbative orders.

Let me summarize the orders of perturbation theory are used in this work:

• •

Hard coefficient function $H_{ff^{\prime}}(\mu,Q)$ in (1) is taken at NNLO (i.e. up to $a_{s}^{2}$-terms inclusively) Gehrmann:2010ue .

• •

Matching coefficient for unpolarized TMDPDF $C_{f\to f^{\prime}}(x,b)$ in (6) is taken at NNLO (i.e. up to $a_{s}^{2}$-terms inclusively) Echevarria:2016scs , in $\zeta$-prescription Scimemi:2017etj .

• •

The perturbative part of rapidity anomalous dimension $\mathcal{D}_{\text{perp}}(\mu,b)$ in (9) is taken at NNLO (i.e. up to $a_{s}^{2}$-terms inclusively) Echevarria:2015byo , in the resummed form Echevarria:2012pw ; Scimemi:2018xaf .

• •

The $\zeta_{\text{pert}}(\mu,b)$ in (12) is taken at NNLO (i.e. up to $a_{s}^{2}$-terms inclusively) Scimemi:2018xaf .

• •

The $\zeta_{\text{exact}}(\mu,b)$ in (12) is taken at NNLO (i.e. up to $a_{s}^{1}$-terms inclusively), see eqns.(23-26).

• •

To evaluate expressions for last three points one needs cusp anomalous dimension and $\gamma_{V}$ anomalous dimension up to $a_{s}^{3}$-terms and $a_{s}^{2}$-terms, receptively. They could be found in Moch:2004pa ; Moch:2005tm .

Thus, the computation is done at complete NNLO perturbative accuracy.

3 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 Badier:1982zb , E537 Anassontzis:1987hk and E615 Conway:1989fs experiments. The measurement by NA3 is presented only by a plot in ref.Badier:1982zb , and the exact values of data-points and their error-bars are not available. Therefore, only the visual comparison with NA3 is possible (fig.6). The data tables for E537 Anassontzis:1987hk and E615 Conway:1989fs can be can be found in Stirling:1993gc .

Both experiments E537 and E615 have been performed in the same environment at different energies of the pion beam, $P_{\text{beam}}=125$GeV for E537 and $P_{\text{beam}}=252$ GeV for E615, which corresponds to $s=235.4$GeV2 and $s=473.6$GeV2, 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_{\pi}$ is $0.26$ (for E537) and $0.18$ (for E615).

Definition of $\chi^{2}$ distribution.

To estimate the theory-to-data agreement, I have used the $\chi^{2}$-test function contracted as usual

[TABLE]

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 $\sigma_{\text{stat}}$ and a systematic error. The later is mainly generated by the luminosity uncertainty, and thus can be considered as a correlated error $\sigma_{\text{corr}}$. The covariance matrix is be build according to general rules:

[TABLE]

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 $\chi^{2}$ calculated with covariance matrix (14) is extremely high, despite the errors of measurements are relatively large. It happens because the correlated part of $\chi^{2}$ overweights the uncorrelated part by an order of magnitude. So, the fit procedure becomes impossible. To stabilize the values of $\chi^{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}<\delta\cdot Q$, where $\delta\approx 0.25$ as it has been derived in Scimemi:2017etj 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}\lesssim 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 $\chi^{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 non-perturbative parameters only the E615 differential in $x_{F}$ data-set has been used. Furthermore, the bins $x_{F}\in(0.8,0.9)$ and $x_{F}\in(0.9,1.0)$ have been excluded, because there $x_{\pi}\sim 1$ and thus, the threshold resummation must be applied. The resulting set has 80 points.

Dependence on collinear PDFs.

According to (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 Ball:2017nwa . Additionally, the extraction of TMDPDFs and $\mathcal{D}_{\text{NP}}$ based on different collinear PDFs were performed (the analysis of these results will be presented elsewhere TOBE ), and they are available at web .

In the present study, I have compared the predictions generated with proton TMDPDFs (and $\mathcal{D}_{\text{NP}}$) based on different collinear PDF, and found results alike. Particularly, $\chi^{2}$-minimization with proton TMDPDFs based on MMHT14(nnlo) Harland-Lang:2014zoa , NNPDF3.1 (nnlo) Ball:2017nwa and HERA20PDF (nnlo) Abramowicz:2015mha gives $\chi^{2}/N_{p}=1.45$, $1.70$ and $1.44$, correspondingly. Taking into account, that HERA20PDF set also shows better global $\chi^{2}$ on the data-set from ref.Bertone:2019nxa , 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 web . For pion collinear PDF JAM18pionPDF-set has been used Barry:2018ort .

Results of the fit.

The minimization procedure for $\chi^{2}$-test yields the following values of non-perturbative parameters

[TABLE]

The first error-band is due to the uncertainty of data-points. It is estimated by the replica method, as in ref.Ball:2008by , by minimization of $\chi^{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 $\chi^{2}$ on 100 of replicas of input distributions.

Parameters $a_{1,2,3}$ are restricted to positive values. So, large error-bands in (15) 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_{\text{NP}}$, since all parameters are correlated. For example, at $b\sim 0.5$GeV*-1* the uncertainty in $f_{\text{NP}}$ is $2-3\%$. However, this band is definitely biased by the ansatz (8). The plot for $f_{\text{NP}}$ is shown in fig.1(left). The actual values of TMDPDF in $b-$space and $k_{T}$-space (that is obtained by Fourier transformation) are shown in fig.1(center,right). The pion TMDPDF obtained in this work together with distribution of 100 replicas is available in the artemide-repository web as Vpion19 TMDPDF set (for $\pi^{-}$-meson).

The final values of $\chi^{2}$ is $\chi^{2}/N_{p}=1.44$ $(N_{p}=80)$. It can be compared with the result of fit in ref.Wang:2017zym $\chi^{2}/N_{p}=1.64$, where almost the same data were used. The main contribution to the value of $\chi^{2}$ comes from the systematic disagreement in the normalization between the data and the theory. In fig.2,3,4,5 the comparison of the data to the theory prediction is shown together with the values of $\chi^{2}/N_{p}$ for a given subset of data-points. In fig.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}\simeq 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 $\chi^{2}$ values as

[TABLE]

where $\chi_{D}^{2}$($\chi_{\lambda}^{2}$) represents the uncorrelated(correlated) part of $\chi^{2}$. Loosely speaking, the value of $\chi_{D}^{2}$($\chi_{\lambda}^{2}$) demonstrates the agreement in the shape (normalization) between the theory and the data. The decomposition (16) is done with the help of nuisance parameters Ball:2008by ; Ball:2012wy . 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 fig.2,3,4,5 for each bin for E615 and common for E537.

The decomposition of $\chi^{2}$ for the selected data is

[TABLE]

The value $\chi_{\lambda}^{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 (fig.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 (fig.6) also does not show any normalization problem. Neglecting the normalization part of the $\chi^{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 fig.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.Bacchetta:2019tcu . The authors of ref.Bacchetta:2019tcu 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.Bacchetta:2019tcu belong to the same kinematic domain as the data discussed here. The comparison has been done in the regime $q_{T}\sim 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 Bertone:2019nxa (for the same experiments that are considered in Bacchetta:2019tcu ), namely, about 40% deficit in the normalization that decreases with the increase of energy (see table 3 in Bertone:2019nxa ). Note, both analyses Bacchetta:2019tcu and Bertone:2019nxa have not a problem with the description of PHENIX data Aidala:2018ajl 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) Gonzalez-Hernandez:2018ipj .

Previously, E615 measurement have been analyzed in the framework of TMD factorization in refs.Wang:2017zym and Ceccopieri:2018nop . 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.Wang:2017zym ; Ceccopieri:2018nop .

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 (7). Generally, an extra factor $R_{i}^{A}(x)$ for PDF should be added, see e.g.Eskola:2016oht . Typically, at $x\in(0.1,0.9)$ this factor provides $\sim 10\%$ modification Eskola:2016oht , which cannot compensate the gap between the theory and the data. Moreover, this effect should be much smaller for $x$-integrated bins (fig.5) due to the oscillation of $R_{i}^{A}(x)$ between anti-shadowing and EMC regimes.

• •

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 Bertone:2019nxa ; TOBE . 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 Barry:2018ort , 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 (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$ (fig.2,3). However, the effect of threshold resummation should be negligible at $x\sim 0.2$ and $Q\sim 4-5$GeV, where the most significant deviation takes place. Also in ref.Bacchetta:2019tcu 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 Braun:2011dg ). 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\to 0$, and so, their presence would be indicated in the deformation of the shape of cross-section, what is not observed.

• •

Wrong shape for non-perturbative corrections. It could happen that the suggested ansatz for non-perturbative parts of TMD evolution (9) and TMDPDF (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 Bertone:2019nxa .

• •

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 $\gamma^{*}$ with $J/\psi$, $\psi^{\prime}$ resonances and their excitations that are located in the region $Q\sim 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\sim 9$GeV (the same with $x_{F}$-differential bins since $x_{F}\sim Q/\sqrt{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 Conway:1989fs . It was observed that $q_{T}$-spectrum after subtraction of normalization has an extreme dependence on $x_{F}$ (see sec.V.B, and appendix A, in ref.Conway:1989fs ), 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.

4 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 $\zeta$-prescription Scimemi:2018xaf 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 TMDPDF. 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 web .

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\sim 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 sec.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 Bacchetta:2019tcu .

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 Aghasyan:2017jop . In fig.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 fig.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.

Acknowledgements.

I thank Wen-Chen Chang for the correspondence that initiated this work, and for critical remarks and suggestions.

Appendix A Special null-evolution line at large $b$

The concept of the special null-evolution line plays the central role in $\zeta$-prescription. The $\zeta$-prescription, the double evolution and properties of TMD evolution have been elaborated in ref.Scimemi:2018xaf , where I refer for further details. In this appendix, I derive the (perturbative) 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 $\mathbf{E}=(\gamma_{F}(\mu,\zeta)/2,-\mathcal{D}(\mu,b))$ in the plane $(\mu,\zeta)$. The anomalous dimension $\gamma_{F}$ is the ultraviolet anomalous dimension of TMD operator. It has the following form

[TABLE]

where $\Gamma_{\text{cusp}}$ is the cusp-anomalous dimension, and $\gamma_{V}$ anomalous dimension of the vector form-factor. The rapidity anomalous dimension $\mathcal{D}(\mu,b)$ is generally non-perturbative function, which can be computed perturbatively only at small-b, see e.g.Vladimirov:2016dll ; Echevarria:2015byo for NNLO and N3LO computations. It satisfies the renormalization group equation

[TABLE]

Due to this expression the field $\mathbf{E}$ is conservative. Parameterizing an equipotential line as $(\mu,\zeta(\mu,b))$, one finds the following equation for $\zeta(\mu,b)$

[TABLE]

The special null-evolution line is the line that passes thorough the saddle point $(\mu_{0},\zeta_{0})$ of the evolution field. The saddle point is defined as

[TABLE]

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 $\zeta$ at all values of $\mu$ (bigger than $\Lambda_{QCD}$). It follows from the definition of the saddle point, and guaranties the finiteness of perturbative series order-by-order.

The field $\mathbf{E}$, and consequently the equipotential line $\zeta(\mu)$ and the position of the saddle point $(\mu_{0},\zeta_{0})$, depends on $b$, which is treated as a free parameter. It causes certain problems in the implementation of the $\zeta$-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 $\mathcal{D}$. The greater problem is that at larger $b$ the value of $\mu_{0}$ decreases and at some large value of $b$ (typically $b\sim 3$GeV*-1*) $\mu_{0}$ is smaller than $\Lambda_{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.Bertone:2019nxa 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\to 0$, and starts to significantly deviate from exact values at $b\sim 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 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 $\mu$. It could not be done entirely, since scale $\mu$ also enters QCD coupling constant in anomalous dimensions $\Gamma_{\text{cusp}}$ and $\gamma_{V}$. For values of $\mu$ large-enough the value of $a_{s}(\mu)$ is small, and thus the solution could be evaluated order-by-order in $a_{s}(\mu)$. Important, that the non-perturbative dependence is exactly accounted in such approach. The equation (20) can be rewritten

[TABLE]

where $g(\mu,b)=\ln(\mu^{2}/\zeta(\mu,b))$, and $\beta$ is QCD beta-function. In this terms the boundary condition turns into finiteness of function $g$ at $\mathcal{D}=0$. The equation (22) can be easily solved order-by-order in $a_{s}$. Denoting

[TABLE]

[TABLE]

I obtain

[TABLE]

where $p=2\beta_{0}\mathcal{D}/\Gamma_{0}$. Let me mention that NNLO term is exponentially grow at large-$\mathcal{D}$ (the N3LO 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 (24-26) provide a very accurate approximation, since $a_{s}$ is evaluated at $\mu=Q$ and typically $a_{s}=g^{2}/(4\pi)^{2}\sim 10^{-2}$. The most important is that this expression is valid at all values of $b$, even then saddle point is below $\Lambda_{QCD}$.