1 Introduction

The non perturbative transverse structure of hadrons has attracted recently much attention and the issue of extracting transverse momentum dependendent parton distributions (TMDs) from data taken in different processes in present and forthcoming high-luminosity facilities represents an important goal of nowadays hadronic Physics. In particular, Drell–Yan (DY) pair production [1], discussed in this paper, and semi-inclusive deep inelastic scattering are the main processes under investigation [2].

The cross section for DY pair production, differential in the transverse momentum of the pair, \(q_T\), is a particularly suitable observable for this kind of studies. In particular at small \(q_T\), where the TMD formalism is formulated, fixed order calculation of this process show large logarithmic corrections due to an incomplete cancellation of soft and collinear singularities between real and virtual contributions and need to be resummed to all orders to recover the predictivity of the theory [3,4,5].

The description of the \(q_T\) DY spectrum in pp collisions has reached a high degree of sophistication [6]. On one side, theoretical improvements have increased the perturbative accuracy of the predictions [7,8,9,10,11]. On the other side, global fits of DY production at different energies have given access to the non perturbative proton transverse structure [12, 13]. Both aspects have received increasing attention due to the formalisation of new and old concepts in the TMD language [14,15,16,17]. While there are differences between the language used in the modern and the older TMD approaches, physical results should not depend on it. A detailed comparison of the formalisms can be found in Refs. [18, 19].

At high energy colliders, this improved knowledge aims to an increasingly better description of electroweak bosons production, with the Higgs \(q_T\) spectrum being the highlighted case. Measurements of \(q_T\) spectrum of the DY process, at lower centre of mass energies, are instead more sensitive to the hadronic non perturbative transverse structure.

DY pair production in pion–nucleus scattering is a unique probe of pion parton distribution functions (PDFs) and, as such represents a source of information on the pion parton structure. In particular for the \(q_T\) spectrum this was realized long time ago by the authors of Ref. [20]. More recently, phenomenological analyses have appeared [21]. A fit to the \(q_T\) spectrum of DY pairs produced in pion–nucleus collisions has been recently presented in Ref. [22].

Pion TMDs, which could be extracted in principle in a next generation of pion–nucleus DY experiments [23], have received recently considerable theoretical interest [21, 24,25,26, 28,29,30,31]. In this paper we study the DY unpolarized pair production in pion–nucleus scattering, to next-to-leading logarithmic (NLL) perturbative accuracy, up to a transverse momentum of the produced lepton pair of 2 GeV. As non perturbative inputs, we use, for the bound nucleons, a longitudinal structure which takes into account nuclear effects and, for the transverse structure, a well established parameterization obtained through a phenomenological fit to proton–proton DY data (called, from now on, KN05 prescription) [12]. For the pion, we use TMDs obtained in a recent calculation [30], within a Nambu–Jona Lasinio (NJL) framework [32], with Pauli–Villars regularization. The corresponding RGE scale of the model is determined in a novel way by comparing the DY unpolarized cross section, integrated over \(q_T\), described by evolved pion PDFs evaluated in the NJL model to the data.

The aim of the present paper is to study the performances of the NJL model, widely used to describe the non-perturbative meson structure, against DY differential cross section data for the first time. We also analyze to what extent this process can be used to obtain information on the pion transverse structure in momentum space, as it happens for the proton in the corresponding process.

The paper is structured as follows. In the next section, we present the set-up of the calculation and introduce the ingredients used to describe the proton and pion structure. In the third section, we discuss the results of the calculation of DY cross sections in the kinematics of presently available data for pion-tungsten scattering. Eventually, we draw our conclusions in the last section.

2 Setting-up the calculation

2.1 Drell–Yan cross section

In the following we will be interested in the process of the type

$$\begin{aligned} h_1(p_1) \; h_2(p_2)\rightarrow \gamma ^*(q)+X, \end{aligned}$$
(1)

in which a virtual photon is produced with large invariant mass \(Q^2\) and transverse momentum \(q_T\) in the collisions of two hadrons at a centre-of-mass energy \(s = (p_1 + p_2)^2\), with \(p_{1,2}\) the four momentum of hadrons \(h_{1,2}\), respectively. When \(q_T^2\) becomes small compared to \(Q^2\), large logarithmic corrections of the form of \(\alpha _s^n \log ^m(Q^2/q_T^2)\) with \(0 \le m \le 2n-1\) appear in fixed order results, being n the order of the perturbative calculation. These large logarithmic corrections can be resummed to all orders by using the Collins-Soper-Sterman (CSS) formalism [6]. In this limit, of interest for the present analysis and neglecting finite corrections in the \(q_T\sim Q\) region, the cross-section can be written as

$$\begin{aligned} \frac{d\sigma }{dq_T^2 d\tau dy}= & {} \sum _{a,b} \sigma _{q\bar{q}}^{(LO)} \int _0^\infty db \frac{b}{2} J_0(b \, q_T)\,S_q(Q,b)\, S_{NP}^{h_1 h_2}(b) \nonumber \\&\times \left[ (f_{a/h_1}\otimes C_{qa})\left( x_1,\frac{b_0^2}{b^2}\right) (f_{b/h_2}\otimes C_{\bar{q}b})\left( x_2,\frac{b_0^2}{b^2}\right) \right. \nonumber \\&\left. \phantom {\left[ (f_{a/h_1}\otimes C_{qa})\left( x_1,\frac{b_0^2}{b^2}\right) (f_{b/h_2}\otimes C_{\bar{q}b})\left( x_2,\frac{b_0^2}{b^2}\right) \right. }+ q \leftrightarrow {\bar{q}} \right] . \end{aligned}$$
(2)

where \(b_0=2 e^{-\gamma _e}\), the symbol \(\otimes \) stands for convolution and \(\sigma _{q\bar{q}}^{(LO)}\) is the leading-order total partonic cross section for producing a lepton pair, \(\sigma (q {\bar{q}} \rightarrow l^+ l^-)\), and it is given by

$$\begin{aligned} \sigma _{q\bar{q}}^{(LO)} = {4 \pi \alpha _{em}^2 \over 9 Q^2} e_q^2. \end{aligned}$$
(3)

In Eq. (2), the a, b indices run on quark and gluons, \(J_0(b \, q_T)\) is the Bessel function of first kind and \(f_{i/h}\) corresponds to the distribution of a parton i in a hadron h. The cross section in Eq. (2) is differential in \(\tau =Q^2/s\) and y, the rapidity of the DY pair. Momentum fractions appearing in parton distribution functions can be expressed in terms of these variables as

$$\begin{aligned} x_{1(2)} = \sqrt{\tau } e^{\pm y}, \quad y=\frac{1}{2} \ln \frac{x_1}{x_2}. \end{aligned}$$
(4)

Cross sections differential in \(x_F=x_1-x_2= 2 q_\parallel / \sqrt{s}\), the longitudinal momentum of the pair in the hadronic centre of mass system, can be obtained from those differential in rapidity y by a suitable transformation. By defining \(A=\sqrt{x_F^2 + 4 \tau }\) one gets

$$\begin{aligned} x_1= {x_F \over 2} + {A \over 2}, \quad x_2= -{x_F \over 2} + {A \over 2}, \quad dy = dx_F / A . \end{aligned}$$
(5)

Momentum conservation further imposes that \(|x_F|<1-\tau \). The large logarithmic corrections are conventiently exponentiated in b-space in the Sudakov perturbative form factor

$$\begin{aligned}&S_q(Q,b)\nonumber \\&\quad =\exp \left\{ -\int _{b_0^2/b^2}^{Q^2} \frac{dq^2}{q^2} \left[ A(\alpha _s(q^2)) \;\ln \frac{Q^2}{q^2} + B(\alpha _s(q^2)) \right] \right\} . \end{aligned}$$
(6)

The functions \(C_{ab}\) in Eq. (2) and A, B in Eq. (6) have perturbative expansions in \(\alpha _s\),

$$\begin{aligned} A(\alpha _s)= & {} \sum _{n=1}^\infty \left( \frac{\alpha _s}{2\pi } \right) ^n A^{(n)}, \nonumber \\ B(\alpha _s)= & {} \sum _{n=1}^\infty \left( \frac{\alpha _s}{2\pi } \right) ^n B^{(n)}, \end{aligned}$$
(7)
$$\begin{aligned} C_{ab}(\alpha _s,z)= & {} \delta _{ab} \,\delta (1-z) + \sum _{n=1}^\infty \left( \frac{\alpha _s}{2\pi } \right) ^n C_{ab}^{(n)}(z). \end{aligned}$$
(8)

At present, the perturbative Sudakov form factor can be evaluated at next-to-next-to-leading logarithmic (NNLL) accuracy [11]. In the \(q\bar{q}\) annihilation channel pertinent to Drell–Yan production, the evaluation of the Sudakov form factor at next-to-leading logarithmic (NLL) accuracy, the one reached in the present analysis, involves the coefficients

$$\begin{aligned} A^{(1)}=2 C_F \quad B^{(1)}=-3 C_F, \end{aligned}$$
(9)

which are the coefficient of the singular \((1-z)^{-1}\) and \(\delta (1-z)\) terms of the one-loop splitting function \(P_{qq}^{(0)}(z)\) and

$$\begin{aligned} A^{(2)}= K A^{(1)}, \quad K=C_A \left( \frac{67}{18}- \frac{\pi ^2}{6} \right) - n_f T_R \frac{10}{9}, \end{aligned}$$
(10)

which is the coefficient of the singular term of the two-loop splitting function \(P_{qq}^{(1)}(z)\) in the \(z\rightarrow 1\) limit [34]. The general expression for \(C_{ab}^{(1)}\) are given by [11, 35]

$$\begin{aligned} C^{(1)}_{qa} (z)= & {} C^{(1)}_{\bar{q}b} (z)\nonumber \\= & {} \delta _{qa} C_F\, (1-z) + \delta _{qa}\, \delta (1-z) C_F \left( -4 + \frac{\pi ^2}{2} \right) \,,\nonumber \\ C^{(1)}_{qg} (z)= & {} C^{(1)}_{\bar{q}g} (z)= 2 T_R\, z (1-z). \end{aligned}$$
(11)

Color factors in the previous equations are given by \(C_A=3\), \(C_F=4/3\), \(T_R=1/2\) with \(n_f\) being the number of active flavours. Together with the use of NLO pdfs, this guarantees the evaluation of the cross section at small \(q_T\) at NLL accuracy. The last ingredient in Eq. (2) is the non perturbative form factor, \(S_{NP}^{h_1 h_2}(b)\), which encodes the transverse structure of both the colliding hadrons. The latter is either fixed by comparison with data or parametrized with the help of hadronic models, as we shall do in this paper.

2.2 Proton structure

Predictions for the transverse momentum spectrum of DY pairs produced in pion–proton collisions do rely on the knowledge of the proton NP form factor. The latter is extracted from the transverse momentum spectrum of DY pairs produced in proton–proton (pp) and proton–nucleus (pA) collisions. Quite recent analyses [15, 16] have appeared which address such an extraction. Since our aim here is to establish the possibility of studying the pion transverse non perturbative structure in pion–nucleus DY experiments, we here intend to minimize the uncertainity coming from the proton structure part of the calculation. We use the well known and widely accepted results of Konychev and Nadolsky (KN05) [12] obtained within the CSS formalism [6] where \(S_{NP}^{pp}(b)\) is extracted from global fit to Z-boson and low mass DY data, updating the results presented in Ref. [13]. The latter is parametrised as

$$\begin{aligned}&S_{NP}^{pp}(b)\\&=\exp \{-[a_1 + a_2 \ln (M/(3.2 \,\text{ GeV })) + a_3 \ln (100 x_1 x_2)] b^2\}.\nonumber \end{aligned}$$
(12)

The \(a_i\) parameters appearing in Eq. (12) are determined by a minimisation procedure against data and are given by [12]

$$\begin{aligned} a_1=0.201 \pm 0.011,\quad&\;\; a_2=0.184 \pm 0.018, \nonumber \\ a_3=&-0.026 \pm 0.007. \end{aligned}$$
(13)

The fit is fully specified once a prescription for the treatment of the non perturbative, large-b, region both in the Sudakov form factor, Eq. (6), and the parton distributions is given. The authors of Ref. [12] adopt the so-called \(b_\star \)-prescription, substituting b with

$$\begin{aligned} b_\star (b,b_{max})=\frac{b}{\sqrt{1+\left( \frac{b}{b_{max}}\right) ^2}}, \end{aligned}$$
(14)

and setting \(b_{max}=1.5 \, \text{ GeV }^{-1}\) in the perturbative form factor. In principle, the same setting should be used in PDFs, which are evaluated at the factorisation scale \(\mu _F=b_0/b_*\). However this choice for \(b_{max}\) may imply a call to a specific PDFs parameterization below their lowest available scale, \(Q_{in}\). Since in Ref. [12] cross sections are evaluated with the NLO CTEQ6M PDFs [36], whose lowest Q accessible is \(Q_{in}=1.3\) GeV, the \(b_\star \)-prescription entering PDFs calls is used with \(b_{max}=b_0/Q_{in} \simeq 0.86 \, \text{ GeV }^{-1}\) which always guarantees \(\mu _F>Q_{in}\). It is important to remark that the non perturbative form factor is determined not only by fitting the parameters of the chosen functional form, but also by the specific regularisation prescription and its associated parameters adopted to deal with the infrared region. In general all these ingredients have been found to be highly correlated.

Fig. 1
figure 1

Theoretical predictions obtained with the KN05 model [12] compared to DY transverse momentum spectra in pA collisions [37] in bins of the invariant mass of the pair, M, expressed in GeV, for different incident beam energies and DY pairs rapidities. Solid lines indicate predictions in the phase space region included in the KN05 fit whereas dashed ones indicate predictions in an extrapolation regime. The error band corresponds to the \(a_i's\) error propagation

Full size image

In order to present a benchmark of our code and to gauge how theory performs in extrapolation regions, we compare predictions from KN05 to the pA data of Ref. [37]. An additional \(\pm 25 \%\) normalisation error is assigned to the data [37]. In the original KN05 analysis, only the data at \(p_{lab}=400\) GeV, \(q_T<1.4\) GeV, \(5<M/\text{ GeV }<9\) were included in the fit. In such a restricted region indeed the theory (solid lines) performs well offering a good benchmark of our code, as shown in the first row of Fig. 1. Since the \(\pi W\) data to be analyzed in the following are at \(p_{lab}\)=252 GeV, it is important to check how well the theory performs in extrapolation regions at lower \(\sqrt{s}\) and higher DY rapidity. Therefore we present in the second and third rows of Fig. 1 the KN05 benchmark (dashed lines) versus data [37] at \(p_{lab}\)= 200 and 300 GeV, which were not included in the KN05 fit. By using Eq. (4) and Eq. (5) and assuming the invariant mass values indicated on the plots, the rapidity coverage of these data can be converted to the range \(0<x_F<0.3\). In both cases we find good agreement between data and theory up to \(q_T \sim 2\) GeV giving us confidence that the KN05 model can be successfully used in this (\(x_F,q_T\)) range at the \(\sqrt{s}\) of interest in this analysis.

2.3 Pion structure

A calculation of pion TMDs in a NJL framework, with Pauli–Villars regularisation, has been recently presented in Ref. [30]. Model calculations of meson partonic structure within this approach have a long story of successful predictions [38,39,40,41,42,43]. Collinear parton distributions obtained within a model have to be associated to a low momentum scale \(Q_0^2\) and, in order to be used to predict measured quantities, have to be evolved to higher momentum scales according to perturbative QCD (pQCD).

In Ref. [30] the unpolarized NJL parton TMD has been obtained. Among its good properties, we stress that, upon integration over the intrinsic quark transverse momentum \(\varvec{k_T}\), the pion PDF q(x) is properly recovered with correct normalisation and the momentum sum rule is exactly satisfied. This is due to the fact that NJL is a field theoretical scheme and the correct support of the PDF, \(0\le x\le 1\), is not imposed but arises naturally. In particular, the momentum sum rule reads \(\int dx~x~q\left( x\right) =0.5\), i.e. the fraction of momentum carried by each quark is one half of the total momentum, since at the scale of the model only valence quarks are present. The dependence on \(\varvec{k_T}\) of the TMD obtained in Ref. [30] is very important for the present study. It is worth stressing that, in this approach, the \(\varvec{k_T}\) dependence is automatically generated by the NJL dynamics and it is not imposed by using any educated guess. This is an important feature of the results of Ref. [30], not found in other approaches [21, 27]. In this paper we will use the pion TMD obtained in Ref. [30] in the chiral limit, the latter allowing for a factorisation of the x and \(k_T\) dependence at the low but undetermined scale \(Q_0^2\) associated to the model:

$$\begin{aligned} f^{q/\pi }(x_{\pi },\varvec{k_T},Q_0^2)= q(x_{\pi },Q_0^2) T(\varvec{k_T})~, \end{aligned}$$
(15)

where one has (in \(\pi ^-\), of interest here):

$$\begin{aligned} q(x_{\pi },Q_0^2)=d_v(x_{\pi },Q_0^2)=\bar{u}(x_{\pi },Q_0^2)=1. \end{aligned}$$
(16)

The function T is given by

$$\begin{aligned} T(\varvec{k_T})= \frac{3}{4 \pi ^3} {\left( m \over f_\pi \right) ^2}~ \sum _{i=0,2} \frac{c_i}{k_T^2+m_i^2}\,, \end{aligned}$$
(17)

which, due to a proper combination of the \(c_i\) [30], behaves as \(k_T^{-6}\) for asymptotic values of \(k_T=|\varvec{k_T}|\) and satisfies the normalisation

$$\begin{aligned} \int d^2 \varvec{k_T} \, T(\varvec{k_T})= 1. \end{aligned}$$
(18)

Since the distribution Eq. (17) depends only upon \(k_T^2\), its Fourier transform can be cast in the form

$$\begin{aligned} S_{NP}^{\pi }(b)= & {} \frac{3}{2 \pi ^2} {\left( m \over f_\pi \right) ^2}~ \sum _{i=0,2} \int \, dk_T \, k_T \, J_0(b k_T) \frac{c_i}{k_T^2+m_i^2} \nonumber \\= & {} \frac{3}{2 \pi ^2} {\left( m \over f_\pi \right) ^2}~ \sum _{i=0,2} c_i K_0 ( m_i \, b)\,, \end{aligned}$$
(19)

where \(K_0\) is the modified Bessel function of the second kind. The parameters used in Eq. (19) are given in Ref. [30] and read

$$\begin{aligned} m_0^2= & {} m^2=(0.238\;\text{ GeV })^2, \\ m_1^2= & {} m^2+ \Lambda ^2, \quad m_2^2=m^2+2\Lambda ^2 ,\\ \Lambda= & {} 0.860 \;\text{ GeV },\\ c_0= & {} 1, \quad c_1=-2, \quad c_2=1, \\ f_\pi= & {} 0.0924 \; \text{ GeV }. \end{aligned}$$

We notice that, beyond the chiral limit, the factorised expression (15) is violated only slightly due to a non trivial x-dependence (see, e.g., Ref. [46]). We have therefore used the expression valid in the chiral limit to avoid further irrelevant complications in an evaluation which is already rather involved.

Fig. 2
figure 2

Drell–Yan pairs production in \(\pi ^- W\) collisions. Next-to-leading order cross sections obtained by using evolved NJL pion PDFs for three values of \(Q_0^2\) are compared to data of Ref. [44]

Full size image

As noted above, the NJL pion model corresponds to a low hadronic scale \(Q_0^2\). Such a low scale has been determined previously by directly comparing the second moment of the pion PDF evaluated in NJL model with the results from the analysis of Ref. [48]. The procedure gives a value of \(Q_0^2=0.18\) GeV\(^2\) at NLO [45, 46].Footnote 1 In the present paper we use a different strategy: we consider \(Q_0^2\) a free parameter of the NJL model which is then fixed with a minimisation procedure, outlined in the following, of the theoretical \(\pi ^- W\) DY cross sections, differential in \(\sqrt{\tau }\) and \(x_F\), against the corresponding experimental ones [44]. Theoretical cross sections are calculated according to

$$\begin{aligned}&\frac{d^2 \sigma }{dQ^2 dx_F}=\frac{4\pi \alpha _{em}^2}{9 Q^2 s}\nonumber \\&\sum _{ij} e_i^2 \int _{x_1}^{1} dt_1 \int _{x_2}^{1} dt_2 \frac{d^2 \hat{\sigma }^{ij}}{dQ^2 dx_F} f_{i/\pi }(t_1,Q^2) f_{j/p}(t_2,Q^2)\,,\nonumber \\ \end{aligned}$$
(20)

where the partonic cross sections \(d\hat{\sigma }^{ij}\) are calculated at NLO accuracy by using the results of Ref. [48]. An additional correction takes into account the correct number of gluon polarisations in the \(\overline{\text{ MS }}\) in dimensional regularisation [49]. The NJL pion PDFs are evolved to NLO accuracy in the Variable Flavor Number Scheme, with the initial condition given in Eq. (16), with the help of the QCDNUM [50] evolution code. The QCD parameters are those of the NLO CTEQ6M parameterisation [36]. In particular we set the NLO running coupling to \(\alpha _s^{(n_f=5)}(M_Z)=0.118\) at the Z-boson mass, \(M_Z\). Since the data we are comparing to are obtained on a tungsten target, we take into account nuclear effects by using nuclear PDFs of Ref. [51]. We have carried out a \(\chi ^2\) study to establish the hadronic scale of the model that describes the best the data at NLO in pQCD. Two cases have been considered: an evaluation of the \(\chi ^2\) for the full range of \(x_F\) and another one with a cut \(x_F<0.4\), since the NJL model is expected to better reproduce the pion valence distributions, expected to populate the range of large and positive \(x_F\) . The scales thus determined are

$$\begin{aligned} Q_{0,\, \text{ no } \text{ cut }}^2= & {} 0.212_{-0.012}^{+0.011} \; \text{ GeV }^2, \quad \nonumber \\ Q_{0,\, \text{ cut }}^2= & {} 0.209_{-0.009}^{+0.008} \; \text{ GeV }^2, \end{aligned}$$
(21)

and correspond to a chisquare value of \(\chi ^2/\)d.o.f.\(=2.1\) and 1.9, respectively. The quoted errors correspond to a variation of one unit in \(\chi ^2\), i.e. 1-\(\sigma \). Those results are compatible with each other. We will therefore refer to \(Q_0^2=0.21\) GeV\(^2\), as the scale associated to the pion NJL model. The other two curves in Fig. 2, corresponding to \(Q_0^2=0.19\) GeV\(^2\) and \(Q_0^2=0.25\) GeV\(^2\) respectively, are added, in order to show the sensitivity to this particular choice of infrared \(Q_0^2\). It is worth noticing that the results show an acceptable agreement, both in shape and in normalisation. More in detail, a tendency of the theory to undershoot the data is identified in the range of small \(x_F\) (\(-0.2<x_F<0.2\)). This deficiency is not unexpected since, in the mentioned kinematic region, the dominant contribution to the cross sections involves sea quarks and gluons which are absent at \(Q_0^2\) and are radiatively generated by QCD evolution. This is a typical drawback of models which contain only valence contributions at the hadronic scale. At this point we would like to mention that the theoretical description of the \(x_F\)-spectra at large \(x_F\) and the determination of pion parton distributions can be further improved employing resummation techniques presented in Refs. [49, 52, 53]. It is worth noticing that, as shown in those papers, threshold NLL resummation of the Wilson coefficients leads to larger cross sections at large x with respect to NLO ones. This, in turn, implies softer pion PDFs at large x. In the present context, this fact would imply a scale \(Q_0^2\) for the NJL model lower than the one already determined by using NLO Wilson coefficients in Eq. (20).

Fig. 3
figure 3

Transverse profile in b-space for the NJL pion, Eq. (19), compared, in the top panel, to the profile of the WLS pion [22] and, in the bottom panel, to that of the KN05 proton, \(\sqrt{S_{NP}^{pp}(b)}\), both evaluated for different values of the scale M

Full size image

3 Predictions for \(\pi W\) collisions data

Predictions for the \(\pi W\) Drell–Yan cross sections are obtained once appropriate modifications are implemented in Eq. (2). Evolved NJL pion parton distributions replace proton PDFs for hadron 1. Moreover the non-perturbative form factor \(S_{NP}^{h_1 h_2}(b)\) depends on the particle species initiating the reaction. Therefore in \(\pi W\) collisions the latter is written as follows:

$$\begin{aligned} S_{NP}^{\pi W}(b)=S_{NP}^{\pi }(b) \, \sqrt{S_{NP}^{pp}(b)} \,, \end{aligned}$$
(22)

where \(S_{NP}^{\pi }(b)\) is given in Eq. (19) and the square root on \(S_{NP}^{pp}(b)\), given in Eq. (12), takes into account that now only one proton is involved in the process. It is instructive to directly compare the proton and pion non perturbative transverse distributions used in the calculation. It is important to remark that the NJL pion transverse distribution in Eq. (19) differs from the corresponding proton factor in Eq. (12) in that it does not contain any explicit dependence neither on hard scale M nor on parton fractional momenta. Such a comparison is meaningful at the typical scale for which the transverse form factors and the longitudinal momentum part factorize. For the pion case this happens at the scale \(Q_0^2\) determined in the previous section. For the proton TMD such a scale is ambiguously defined and, according to KN05 analysis, ranges between \(Q_{in}^2\) and \((b_0/b_{max}^{KN05})^2\). Therefore we choose \(M=Q_{in}=1.3\) GeV in Eq. (12) and fix the product \(x_1 x_2=M^2/s\), see Eq. (4), exploiting the \(\pi ^- W\) kinematics with s calculated according to a beam energy of \(p_{lab}=252\) GeV. This comparison is presented in Fig. 3, where our result for the pion non-perturbative form factor, \(S_{NP}^{\pi }(b)\), is also compared to the parametrisation of the non-perturbative pion form factor of Ref. [22] (called hereafter WLS). The approach of Ref.  [22] is rather different from ours, both in the spirit and in the physical ingredients used. As a matter of fact, in that paper the proton non-perturbative form factor, \(S_{NP}^{pp}(b)\), has a structure similar to that of our Eq. (12) and for the pion the same form has been assumed, with the corresponding parameters obtained from a fit of the same cross section data used in the present paper. As a result, the pion non-perturbative form factor depends on both the hard scale and the parton momenta. A fit is then performed up to \(q_T\simeq 3\) GeV. We reiterate that our goal here is not to fit but rather to assume a well known structure for the proton non-perturbative form factor and to test the pure NJL predictions for the pion against the data. The purpose of the comparison with the WLS parametrisation is therefore mainly illustrative and quantitative conclusions can be hardly reached. All the distributions presented in Fig. 3 reduce to unity in the \(b \rightarrow 0\) limit, since they are all normalised to unity in transverse momentum space. In the top panel of Fig. 3 we compare the NJL transverse distribution to the pion parametrisation of Ref. [22] (called hereafter WLS) obtained from a fit of the same cross section data used in the present paper. One may notice that, for this model, the width of the distribution is smaller with respect to the NJL one, implying a larger average transverse momentum. In the bottom panel of Fig. 3 one may notice that the NJL pion transverse distribution develops a larger tail with respect to the gaussian drop of the proton distributions. Moreover the b-space width of the KN05 proton with \(M=1.3\) GeV is larger with respect to the pion one. When transformed back in \(k_T\) space, this implies that the intrinsic transverse momentum in the pion is larger than the one in the proton, in agreement with the general expectations, since the pion is a much smaller system with respect to the proton. It is worth mentioning that both the KN05 and WLS non perturbative form factors have an explicit, althought slightly different, dependence upon the hard scale M, in both cases set equal to the invariant mass of the dilepton pair. Therefore we plot in each panels, as a representative case, the curves corresponding to both form factors evaluated with the scale set to \(M=4\) GeV. Comparing the latter curves to the ones with \(M \sim 1\) GeV, we conclude that the M-dependence generates a sizable non perturbative evolution of the form factor which is more pronounced for KN05 proton model than for the WLS pion model.

We now turn to the discussion of the perturbative part of the Sudakov form factor, Eq. (6). The latter, at variance with its non perturbative counter part, does not depend upon the type of initial state hadrons involved in the scattering process. In principle, the same regularisation procedure should be used both in the Sudakov and in the PDFs. This optimum indeed faces some technical problem, for example the call to PDFs to values outside the boundary of the grid in which they are defined and the different scales at which the transverse distributions are assumed to factorise on the proton and pion side, respectively. In order to accomodate all these different settings, we find useful to split the perturbative form factor in Eq. (6) in a form which allows to use distinct \(b_{max}\) on the proton and pion side:

$$\begin{aligned}&S_q(Q,b) \nonumber \\&\quad \equiv S_q(Q,b_*,b_{max}^p,b_{max}^\pi )\nonumber \\&\quad =\exp \left\{ -\int _{\frac{b_0^2}{b^2_*(b_{max}^p)}}^{Q^2} \frac{dq^2}{2\,q^2} \left[ A(\alpha _s(q^2)) \;\ln \frac{Q^2}{q^2} + B(\alpha _s(q^2)) \right] \right\} \, \nonumber \\&\qquad \times \exp \left\{ -\int _{\frac{b_0^2}{b^2_*(b_{max}^\pi )}}^{Q^2} \frac{dq^2}{2\,q^2} \left[ A(\alpha _s(q^2)) \;\ln \frac{Q^2}{q^2} + B(\alpha _s(q^2)) \right] \right\} .\nonumber \\ \end{aligned}$$
(23)

For the proton parameters, we stick to KN05 settings since the \(a_i\)’s in Eq. (12) optimized for \(b_{max}^p=1.5\) Ge\(\text{ V }^{-1}\). On the pion side there is some freedom in adjusting \(b_{max}^{NJL}\). However the pion TMD shows a \(x-k_T\) factorised structure only at \(Q_0^2\), whose numerical value has been determined in the previous section. Therefore we can expect the \(b_*\)-prescription to involve \(b_{max}^\pi \) values of the order \(b_0/Q_0\sim 2.44\) Ge\(\text{ V }^{-1}\), which will be our default value to be used both in the Sudakov and in NJL pion parton distributions regularisation. We now turn to the comparison to lepton pair \(q_T\)-spectra collected in tables D92–D97 of Refs. [44, 54], measured in \(\pi W\) collisions. Such data actually refer to less differential cross sections with respect to the one appearing in Eq. (2). In this case differential cross sections are integrated over additional variables according to values specified in experimental analyses. We start presenting our results showing, in Fig. 4, cross sections differential in \(q_T\) integrated in \(0<x_F<1\) in various bins of the invariant mass of the pair, M. The comparison is performed up to a \(q_T \sim 2\) GeV, where we have checked that the KN05 gives an adequate description of pp data. All three different predictions, to be discussed in the following, capture the normalisation of the data and share a tendency to slightly overestimate the data at very small \(q_T\) and to underestimate them at larger \(q_T\). This effect progressively disappears increasing the mass of the lepton pair. Comparing the two curves corresponding to \(b_{max}^{NJL}= 2.44\) Ge\(\text{ V }^{-1}\) and \(b_{max}^{NJL}= 1.5\) Ge\(\text{ V }^{-1}\), one may notice a substantial stability upon variation of the regulators on the pion side. On the same plot, in order to investigate the sensitivity to the pion transverse distribution, we additionally show the predictions obtained by substituting the pion transverse factor, Eq. (19), with \(\sqrt{S_{NP}^{pp}(b)}\). As shown in Fig. 3, the non perturbative transverse distributions for the proton and pion differ at low scales. The corresponding curve, indicated with pp on the plot, is barely distinguishable from the other two. Such a comparison supports the hypothesis that the effect of the perturbative evolution, driven by Eq. (23), is to wash away differences in the non perturbative structure found at the hadronic scale. This result implies a reduced sensitivity to non perturbative structure. A quantitative analysis of the results yields that the \(\chi ^2/d.o.f.\) of the shown distributions slowly decreases from values of order 4 to values of order 1 with increasing M from 4 GeV to 13 GeV.Footnote 2 These numbers are remarkably good if one considers that we are presenting pure model predictions without a fitting of parameters a posteriori. Besides, we observe that the agreement with data of the NJL distributions is slightly better than those obtained with a proton-like non-perturbative form factor for the pion, in any M bin, and, more importantly, that the difference in \(\chi ^2/d.o.f.\) reaches 30 % for the lowest values of M. The region of low (but still perturbative) M is therefore selected as the most promising to access non-perturbative details of the pion transverse structure.

Fig. 4
figure 4

Predictions compared to cross sections in various invariant mass bins of the pair integrated in \(0<x_F<1\). Data from Refs. [44, 54]

Full size image
Fig. 5
figure 5

Predictions compared to differential cross sections in various \(x_F\) bins integrated in the mass range \(4<M<8.55\) GeV. Data from Refs. [44, 54]

Full size image

We proceed our discussion presenting in Fig. 5 the comparison between theory predictions against the same data, now integrated in the mass range \(4<M<8.55\) GeV in a number of \(x_F\) bins. We remind the reader that we have verified that the KN05 model gives a satisfactory description of pp data up to \(x_F \sim 0.3\). Up to this \(x_F\) value, as shown in the first row of Fig. 5, the description \(\pi ^-W \) data is fair, as already observed in Fig. 4. Beyond that range, however, the width of the theoretical curves decreases more rapidly than observed in the data, with data substantially undershooted beyond \(q_T \sim 1\) GeV. This effect is more pronounced as \(x_F\) increases. In this region of relatively large pion fractional momenta it would be tempting to invoke, in order to describe the data, an x-dependent non perturbative structure. Such an interesting hypothesis, however, cannot be tested unless fixed order contributions at finite \(q_T\) are included in the calculation. Moreover, given our working assumptions, the failure to agree with data at large values of \(q_T\) was somehow expected. Whether it is due to the breakdown of the \(x-k_T\) factorization, the inclusion of a more complex NP Sudakov form factor or the matching with the so called Y-term, further studies shall be pursued to answer that question.

On the theoretical side, we would like to mention that, in this range of quite large pion parton fractional momenta, the theoretical description of the \(q_T\)-spectrum can be further improved employing joint resummation techniques described in Refs. [55, 56]. In order to better appreciate how the width of theoretical predictions evolves with \(x_F\) (and therefore with \(x_{\pi }\)) and the invariant mass of the lepton pair, we show in Fig. 6 the average transverse momentum of the pair, \(\langle q_T^2 \rangle \), calculated as

$$\begin{aligned} \langle q_T^2 \rangle =\frac{\displaystyle \int _{x_F^{ \text{ min }}}^{x_F^{ \text{ max }}} dx_F \int _{\tau _{ \text{ min }}}^{\tau _{ \text{ max }}} d\tau \int _{0}^{q_T^{2, \text{ max }}} dq_T^2 \; q_T^2 \; \frac{d^3 \sigma }{dx_F d\tau d q_T^2}}{\displaystyle \int _{x_F^{ \text{ min }}}^{x_F^{ \text{ max }}} dx_F \int _{\tau _{ \text{ min }}}^{\tau _{ \text{ max }}} d\tau \int _{0}^{q_T^{2, \text{ max }}} dq_T^2 \; \frac{d^3 \sigma }{dx_F d\tau d q_T^2}} . \end{aligned}$$
(24)

Integration limits are provided by experimental conditions. For data, indicated by black lines in Fig. 6, the phenomenological parametrisation presented in Ref. [44] is used. For both theory and data, the \(\langle q_T^2 \rangle \) is calculated with a maximum value of \(q_T^{max}=2\) GeV. Theory predictions tend to undershoot the data but, overall, a good shape agreement is found. By comparing lines with and without TMD evolution (for the latter the perturbative Sudakov \(S_q\) is removed from the evaluation of Eq. (2)) one can appreciate its large impact on the amount of generated \(\langle q_T \rangle \). On the same plot, in order to investigate the sensitivity to the pion transverse structure, we additionally show the predictions obtained by substituting the pion transverse factor, Eq. (19), with \(\sqrt{S_{NP}^{pp}(b)}\). As already seen in Fig. 4, differences are minimal, implying a reduced sensitivity to details of the non perturbative transverse factor. Therefore if one aims to better appreciate the strictly non perturbative form factor, one has to confine in corners where TMD evolution is minimised, but still in a perturbative range. These phase space regions can be identified by extrapolation from the right plot as the one at the lowest, but still perturbative, values of the invariant masses of the pair, as already noticed above while discussing Fig. 4.

Fig. 6
figure 6

Top panel: lepton pair average transverse momentum, \(\langle q_T^2 \rangle \), as a function of \(x_F\) integrated in the mass range \(4.0<M/\)GeV\( <8.55\). Bottom panel: \(\langle q_T^2 \rangle \) as a function of M integrated in the range \(0<x_F<1\). Averaged values are obtained integrating both predictions and the phenomenological parametrisation of the data up to \(q_T^{max}=2\) GeV

Full size image

4 Conclusions

A thorough analysis of DY pair production in pion–nucleus scattering has been presented. The main goal of our work has been the test of model predictions, obtained within the Nambu–Jona–Lasinio model for the transverse pion structure. In particular we have focused on the study of differential transverse momentum spectra of DY pairs produced in pA collisions calculated in the CSS framework at NLL accuracy borrowing from the literature the longitudinal and transverse proton structure. The pion is treated in the Nambu–Jona–Lasinio model. No further assumption has been made: even the momentum scale associated to the model is obtained via a minimization procedure of NLO theory to DY experimental longitudinal spectra. The latter turns out to be a low one, in line with that normally used, which could be predicted within the spirit of the model without fitting “a posteriori”. The agreement found between our pion–nucleus theoretical cross sections and experimental data is rather successful, confirming the predictive power of the NJL model, for both the longitudinal pion parton distributions and its transverse structure. We notice that the theory tends to systematically undershoot the data on the higher end of the considered \(q_T\) interval. All interpretations of this effect, however, are not conclusive without the inclusion of the finite, fixed order, contributions which populate the \(q_T \sim Q\) region and are neglected in our calculation.

The possibility to distinguish between different non perturbative transverse momentum distributions in DY data appears instead more questionable. In this complicated scenario, a possible strategy would be the measurement of DY pion–nucleus \(q_T\)-spectra, in bins of \(x_F\), at low values of the mass of the pair, as the present study suggests to look into this kinematical window to emphasize the non-perturbative content of the pion. Further analyses of the pion non-perturbative form factor, as a function of the hard scale, should be pursued so we could progress on that point. In the very same window, new data could allow a deeper investigation of the dependence of the non perturbative form factor upon the hard scale of the process.