Extraction of partonic transverse momentum distributions from semi-inclusive deep-inelastic scattering, Drell-Yan and Z-boson production

TL;DR

This paper performs a global fit to extract unpolarized partonic transverse momentum distributions (TMDs) from diverse experimental data, utilizing TMD evolution at next-to-leading logarithmic accuracy, focusing on low transverse momentum regions.

Contribution

It introduces a novel simultaneous fitting method of TMDs across different processes using TMD evolution, marking a first attempt at a comprehensive global analysis.

Findings

Successful extraction of TMDs from multiple data sets

Implementation of TMD evolution at low scales (~1 GeV^2)

Focus on low transverse momentum without high-momentum matching

Abstract

We present an extraction of unpolarized partonic transverse momentum distributions (TMDs) from a simultaneous fit of available data measured in semi-inclusive deep-inelastic scattering, Drell-Yan and Z-boson production. To connect data at different scales, we use TMD evolution at next-to-leading logarithmic accuracy. The analysis is restricted to the low-transverse-momentum region, with no matching to fixed-order calculations at high transverse momentum. We introduce specific choices to deal with TMD evolution at low scales, of the order of 1 GeV$^2$. This could be considered as a first attempt at a global fit of TMDs.

Figures

Tables

Table 1. Table 1: SIDIS proton-target data ( Hermes  experiment).

	Hermes	Hermes	Hermes	Hermes
	$p\to {\pi }^{+}$	$p\to {\pi }^{-}$	$p\to K^{+}$	$p\to K^{-}$
Reference	Airapetian et al. (2013)
Cuts	$Q^2>1.4{\text{ GeV}}^2$
	GeV
Points	190	190	189	187
Max. $Q^2$	$9.2{\text{ GeV}}^2$
$x$ range

Table 2. Table 2: SIDIS deuteron-target data ( Hermes  and Compass  experiments).

	Hermes	Hermes	Hermes	Hermes	Compass	Compass
	$D\to {\pi }^{+}$	$D\to {\pi }^{-}$	$D\to K^{+}$	$D\to K^{-}$	$D\to h^{+}$	$D\to h^{-}$
Reference	Airapetian et al. (2013)				Adolph et al. (2013)
Cuts	$Q^2>1.4{\text{ GeV}}^2$
	GeV
Points	190	190	189	189	3125	3127
Max. $Q^2$	$9.2{\text{ GeV}}^2$				$10{\text{ GeV}}^2$
$x$ range
Notes					Observable: $m_{\text{norm}}(x,z,{𝑷}_{hT}^2,Q^2)$, Eq. (41)

Table 3. Table 3: Low energy Drell–Yan data collected by the E288 and E605 experiments at Tevatron, with different center-of-mass energies.

	E288 200	E288 300	E288 400	E605
Reference	Ito et al. (1981)	Ito et al. (1981)	Ito et al. (1981)	Moreno et al. (1991)
Cuts	GeV
Points	45	45	78	35
$\sqrt{s}$	19.4 GeV	23.8 GeV	27.4 GeV	38.8 GeV
$Q$ range	4-9 GeV	4-9 GeV	5-9, 11-14 GeV	7-9, 10.5-11.5 GeV
Kin. var.	$\eta$=0.40	$\eta$=0.21	$\eta$=0.03	$x_F=0.1$

Table 4. Table 4: Z 𝑍 Z boson production data collected by the CDF and D0 experiments at Tevatron, with different center-of-mass energies.

	CDF Run I	D0 Run I	CDF Run II	D0 Run II
Reference	Affolder et al. (2000)	Abbott et al. (2000)	Aaltonen et al. (2012)	Abazov et al. (2008)
Cuts	GeV
Points	31	14	37	8
$\sqrt{s}$	1.8 TeV	1.8 TeV	1.96 TeV	1.96 TeV
Normalization	1.114	0.992	1.049	1.048

Table 5. Table 5: Total number of points analyzed, number of free parameters and χ 2 superscript 𝜒 2 \chi^{2} values.

Points	Parameters	${\chi }^2$	${\chi }^2/$d.o.f.
8059	11	$12629\pm 363$	$1.55\pm 0.05$

Table 6. Table 6: Number of points analyzed and χ 2 superscript 𝜒 2 \chi^{2} values for SIDIS off a proton target.

	Hermes	Hermes	Hermes	Hermes
	$p\to {\pi }^{+}$	$p\to {\pi }^{-}$	$p\to K^{+}$	$p\to K^{-}$
Points	190	190	189	187
${\chi }^2/$points	$4.83\pm 0.42$	$2.47\pm 0.28$	$0.91\pm 0.14$	$0.82\pm 0.17$

Table 7. Table 7: Number of points analyzed and χ 2 superscript 𝜒 2 \chi^{2} values for SIDIS off a deuteron target.

	Hermes	Hermes	Hermes	Hermes	Compass	Compass
	$D\to {\pi }^{+}$	$D\to {\pi }^{-}$	$D\to K^{+}$	$D\to K^{-}$	$D\to h^{+}$	$D\to h^{-}$
Points	190	190	189	189	3125	3127
${\chi }^2/$points	$3.46\pm 0.32$	$2.00\pm 0.17$	$1.31\pm 0.26$	$2.54\pm 0.57$	$1.11\pm 0.03$	$1.61\pm 0.04$

Table 8. Table 8: Number of points analyzed and χ 2 superscript 𝜒 2 \chi^{2} values for fixed-target Drell–Yan experiments at low energy. The labels in square brackets were introduced in Sec.  III.2 .

	E288 [200]	E288 [300]	E288 [400]	E605
Points	45	45	78	35
${\chi }^2/$points	$0.99\pm 0.09$	$0.84\pm 0.10$	$0.32\pm 0.01$	$1.12\pm 0.08$

Table 9. Table 9: Number of points analyzed and χ 2 superscript 𝜒 2 \chi^{2} values for Z 𝑍 Z boson production at Tevatron.

	CDF Run I	D0 Run I	CDF Run II	D0 Run II
Points	31	14	37	8
${\chi }^2/$points	$1.36\pm 0.00$	$1.11\pm 0.02$	$2.00\pm 0.02$	$1.73\pm 0.01$

Table 10. Table 10: Values of parameters common to TMD PDFs and TMD FFs.

	${\xi }_{\max }$ [GeV-1]	${\xi }_{\min }$ [GeV-1]	$g_2$ [GeV2]
	(fixed)	(fixed)
All replicas	$2e^{-{\gamma }_E}$	$2e^{-{\gamma }_E}/Q$	$0.13\pm 0.01$
Replica 105	$2e^{-{\gamma }_E}$	$2e^{-{\gamma }_E}/Q$	$0.128$

Table 11. Table 11: 68% confidence intervals of best-fit values for parametrizations of TMDs at Q = 1 𝑄 1 Q=1 GeV.

TMD PDFs	$g_1$	$\alpha$	$\sigma$		$\lambda$
	[GeV2]				[GeV-2]
All replicas	$0.28\pm 0.06$	$2.95\pm 0.05$	$0.17\pm 0.02$		$0.86\pm 0.78$
Replica 105	$0.285$	$2.98$	$0.173$		$0.39$
TMD FFs	$g_3$	$\beta$	$\delta$	$\gamma$	${\lambda }_F$	$g_4$
	[GeV2]				[GeV-2]	[GeV2]
All replicas	$0.21\pm 0.02$	$1.65\pm 0.49$	$2.28\pm 0.46$	$0.14\pm 0.07$	$5.50\pm 1.23$	$0.13\pm 0.01$
Replica 105	$0.212$	$2.10$	$2.52$	$0.094$	$5.29$	$0.135$

Table 12. Table 12: χ 2 superscript 𝜒 2 \chi^{2} /d.o.f. for Hermes  data with and without normalization to the value of the first bin in P h ​ T subscript 𝑃 ℎ 𝑇 P_{hT} .

	$p\to {\pi }^{+}$	$p\to {\pi }^{-}$	$p\to K^{+}$	$p\to K^{-}$	$D\to {\pi }^{+}$	$D\to {\pi }^{-}$	$D\to K^{+}$	$D\to K^{-}$
Original	5.18	2.67	0.75	0.78	3.63	2.31	1.12	2.27
Normalized	1.94	1.13	0.57	0.29	1.59	0.80	0.47	0.97

Equations

$$\ell (l)+N(P)\to \ell (l^{\prime })+h(P_h)+X.$$

$$x$$

$$m_N^h(x,z,|{\mathbf{P}}_{hT}|,Q^2)=\frac{d{\sigma }_N^h/(dxdzd|{\mathbf{P}}_{hT}|d{Q}^2)}{d{\sigma }_{\text{DIS}}/(dxd{Q}^2)},$$

$$m_N^h(x,z,|{\mathbf{P}}_{hT}|,Q^2)=\frac{2\pi |{\mathbf{P}}_{hT}|F_{UU,T}(x,z,{\mathbf{P}}_{hT}^2,Q^2)+2\pi \epsilon |{\mathbf{P}}_{hT}|F_{UU,L}(x,z,{\mathbf{P}}_{hT}^2,Q^2)}{F_T(x,Q^2)+\epsilon F_L(x,Q^2)},$$

$$\epsilon$$

$$F_{UU,T}(x,z,{\mathbf{P}}_{hT}^2,Q^2)$$

$$\displaystyle\times x\int d^{2}\bm{k}_{\perp}\,d^{2}\bm{P}_{\perp}\,f_{1}^{a}\big{(}x,\bm{k}_{\perp}^{2};Q^{2}\big{)}\,D_{1}^{a\mbox{\raisebox{-4.5pt}[0.0pt][0.0pt]{$\vec{\phantom{v}}$}}h}\big{(}z,\bm{P}_{\perp}^{2};Q^{2}\big{)}\,\delta^{(2)}\big{(}z{\bm{k}}_{\perp}-{\bm{P}}_{hT}+{\bm{P}}_{\perp}\big{)}$$

$$\displaystyle+Y_{UU,T}\big{(}Q^{2},\bm{P}_{hT}^{2}\big{)}+\mathcal{O}\big{(}M^{2}/Q^{2}\big{)}\,.$$

$$F_{UU,T}(x,z,{\mathbf{P}}_{hT}^2,Q^2)$$

$$\displaystyle\tilde{f}_{1}^{a}\big{(}x,\xi_{T}^{2};Q^{2}\big{)}$$

$$\displaystyle\tilde{D}_{1}^{a\mbox{\raisebox{-4.5pt}[0.0pt][0.0pt]{$\vec{\phantom{v}}$}}h}\big{(}z,\xi_{T}^{2};Q^{2}\big{)}$$

$$A(P_A)+B(P_B)\to [{\gamma }^{\ast }/Z+X\to ]{\ell }^{+}(l)+{\ell }^{-}(l^{\prime })+X.$$

$$\eta=\frac{1}{2}\log\bigg{(}\frac{q^{0}+q_{z}}{q^{0}-q_{z}}\bigg{)}\ .$$

$$\frac{d\sigma }{d{Q}^{2}d{q}_T^{2}d\eta }$$

$${\sigma }_0^{\gamma }$$

$$x_A$$

$$\eta$$

$$F_{UU}^1(x_A,x_B,{\mathbf{q}}_T^2,Q^2)$$

$$\displaystyle\times x_{A}x_{B}\int d^{2}\bm{k}_{\perp A}\,d^{2}\bm{k}_{\perp B}\,f_{1}^{a}\big{(}x_{A},\bm{k}_{\perp A}^{2};Q^{2}\big{)}\,f_{1}^{\bar{a}}\big{(}x_{B},\bm{k}_{\perp B}^{2};Q^{2}\big{)}\,\delta^{(2)}\big{(}{\bm{k}}_{\perp A}-{\bm{q}}_{T}+{\bm{k}}_{\perp B}\big{)}$$

$$\displaystyle+Y_{UU}^{1}\big{(}Q^{2},\bm{q}_{T}^{2}\big{)}+\mathcal{O}\big{(}M^{2}/Q^{2}\big{)}\,.$$

$${\mathcal{H}}_{UU,\gamma }^{1a}(Q^2)$$

$$V_a$$

$$F_{UU}^1(x_A,x_B,{\mathbf{q}}_T^2,Q^2)$$

$${\tilde{f}}_1^a(x,{\xi }_T^2;Q^2)$$

$${\tilde{D}}_1^{a\to h}(z,{\xi }_T^2;Q^2)$$

$${\mu }_b$$

$${\bar{\xi }}_{\ast }$$

$${\xi }_{\max }$$

$$\displaystyle\bigl{(}C_{a/i}\otimes f_{1}^{i}\bigr{)}(x,\bar{\xi}_{\ast},\mu_{b}^{2})$$

$$\displaystyle\bigl{(}\hat{C}_{a/i}\otimes D_{1}^{i\to h}\bigr{)}(z,\bar{\xi}_{\ast},\mu_{b}^{2})$$

$$\displaystyle C_{a/i}\Big{(}\frac{x}{u},\bar{\xi}_{\ast},\alpha_{S}\big{(}\mu_{b}^{2}\big{)}\Big{)}$$

$${\tilde{f}}_1^a(x,{\xi }_T^2;Q^2)$$

$${\tilde{D}}_1^{a\to h}(z,{\xi }_T^2;Q^2)$$

$$S(\mu_{b}^{2},Q^{2})=-\int_{\mu_{b}^{2}}^{Q^{2}}{d\mu^{2}\over\mu^{2}}\bigg{[}A\Big{(}\alpha_{S}(\mu^{2})\Big{)}\ln\bigg{(}{Q^{2}\over\mu^{2}}\bigg{)}+B\Big{(}\alpha_{S}(\mu^{2})\Big{)}\bigg{]}\ ,$$

Full text

††preprint: JLAB THY 17-2437

1. I Introduction

2. II Formalism

3. II.1 Semi-inclusive DIS

4. II.2 Drell–Yan and Z production

5. II.3 TMDs and their evolution

6. III Data analysis

7. III.1 Semi-inclusive DIS data

8. III.1.1 Hermes data

9. III.1.2 Compass data

10. III.2 Low-energy Drell–Yan data

11. III.3 Z-boson production data

12. III.4 The replica method

13. IV Results

14. IV.1 Agreement between data and theory

15. IV.2 Transverse momentum dependence at 1 GeV

16. IV.3 Stability of our results

17. V Conclusions

I Introduction

Parton distribution functions describe the internal structure of the nucleon in terms of its elementary constituents (quarks and gluons). They cannot be easily computed from first principles, because they require the ability to carry out Quantum Chromodynamics (QCD) calculations in its nonperturbative regime. Many experimental observables in hard scattering experiments involving hadrons are related to parton distribution functions (PDFs) and fragmentation functions (FFs), in a way that is specified by factorization theorems (see, e.g., Refs. Collins et al. (1989); Collins (2013)). These theorems also elucidate the universality properties of PDFs and FFs (i.e., the fact that they are the same in different processes) and their evolution equations (i.e., how they get modified by the change in the hard scale of the process). Availability of measurements of different processes in different experiments makes it possible to test factorization theorems and extract PDFs and FFs through so-called global fits. On the other side, the knowledge of PDFs and FFs allows us to make predictions for other hard hadronic processes. These general statements apply equally well to standard collinear PDFs and FFs and to transverse-momentum-dependent parton distribution functions (TMD PDFs) and fragmentation functions (TMD FFs). Collinear PDFs describe the distribution of partons integrated over all components of partonic momentum except the one collinear to the parent hadron; hence, collinear PDFs are functions of the parton longitudinal momentum fraction $x$. TMD PDFs (or TMDs for short) include also the dependence on the transverse momentum $\bm{k}_{\perp}$. They can be interpreted as three-dimensional generalizations of collinear PDFs. Similar arguments apply to collinear FFs and TMD FFs Angeles-Martinez et al. (2015).

There are several differences between collinear and TMD distributions. From the formal point of view, factorization theorems for the two types of functions are different, implying also different universality properties and evolution equations Rogers (2016). From the experimental point of view, observables related to TMDs require the measurement of some transverse momentum component much smaller than the hard scale of the process Bacchetta (2016); Radici (2016). For instance, Deep-Inelastic Scattering (DIS) is characterized by a hard scale represented by the 4-momentum squared of the virtual photon ($-Q^{2}$). In inclusive DIS this is the only scale of the process, and only collinear PDFs and FFs can be accessed. In semi-inclusive DIS (SIDIS) also the transverse momentum of the outgoing hadron ($P_{hT}$) can be measured Mulders and Tangerman (1996); Bacchetta et al. (2007). If $P_{hT}^{2}\ll Q^{2}$, TMD factorization can be applied and the process is sensitive to TMDs Collins (2013).

If polarization is taken into account, several TMDs can be introduced Mulders and Tangerman (1996); Boer and Mulders (1998); Bacchetta and Mulders (2000); Mulders and Rodrigues (2001); Boer et al. (2016). Attempts to extract some of them have already been presented in the past Bacchetta and Radici (2011); Anselmino et al. (2012); Echevarria et al. (2014); Anselmino et al. (2016); Lu and Schmidt (2010); Barone et al. (2015); Lefky and Prokudin (2015); Anselmino et al. (2013); Kang et al. (2016). In this work, we focus on the simplest ones, i.e., the unpolarized TMD PDF $f_{1}^{q}(x,k_{\perp}^{2})$ and the unpolarized TMD FF $D_{1}^{q\to h}(z,P_{\perp}^{2})$, where $z$ is the fractional energy carried by the detected hadron $h$, $k_{\perp}$ is the transverse momentum of the parton with respect to the parent hadron, and $P_{\perp}$ is the transverse momentum of the produced hadron with respect to the parent parton. Despite their simplicity, the phenomenology of these unpolarized TMDs presents several challenges Signori (2016a): the choice of a functional form for the nonperturbative components of TMDs, the inclusion of a possible dependence on partonic flavor Signori et al. (2013), the implementation of TMD evolution Bacchetta et al. (2015); Rogers (2016), the matching to fixed-order calculations in collinear factorization Collins et al. (2016).

We take into consideration three kinds of processes: SIDIS, Drell–Yan processes (DY) and the production of $Z$ bosons. To date, they represent all possible processes where experimental information is available for unpolarized TMD extractions. The only important process currently missing is electron-positron annihilation, which is particularly important for the determination of TMD FFs Bacchetta et al. (2015). This work can therefore be considered as the first attempt at a global fit of TMDs.

The paper is organized as follows. In Sec. II, the general formalism for TMDs in SIDIS, DY processes, and Z production is briefly outlined, including a description of the assumptions and approximations in the phenomenological implementation of TMD evolution equations. In Sec. III, the criteria for selecting the data analyzed in the fit are summarized and commented. In Sec. IV, the results of our global fit are presented and discussed. In Sec. V, we summarize the results and present an outlook for future improvements.

II Formalism

II.1 Semi-inclusive DIS

In one-particle SIDIS, a lepton $\ell$ with momentum $l$ scatters off a hadron target $N$ with mass $M$ and momentum $P$. In the final state, the scattered lepton momentum $l^{\prime}$ is measured together with one hadron $h$ with mass $M_{h}$ and momentum $P_{h}$. The corresponding reaction formula is

[TABLE]

The space-like momentum transfer is $q=l-l^{\prime}$, with $Q^{2}=-q^{2}$. We introduce the usual invariants

[TABLE]

The available data refer to SIDIS hadron multiplicities, namely to the differential number of hadrons produced per corresponding inclusive DIS event. In terms of cross sections, we define the multiplicities as

[TABLE]

where $d\sigma_{N}^{h}$ is the differential cross section for the SIDIS process and $d\sigma_{\text{DIS}}$ is the corresponding inclusive one, and where $\bm{P}_{hT}$ is the component of $\bm{P}_{h}$ transverse to $\bm{q}$ (we follow here the notation suggested in Ref. Boer et al. (2011a)). In the single-photon-exchange approximation, the multiplicities can be written as ratios of structure functions (see Ref. Bacchetta et al. (2007) for details):

[TABLE]

where

[TABLE]

In the numerator of Eq. (4) the structure function $F_{XY,Z}$ corresponds to a lepton with polarization $X$ scattering on a target with polarization $Y$ by exchanging a virtual photon in a polarization state $Z$. In the denominator, only the photon polarization is explicitly written ($T$, $L$), as usually done in the literature.

The semi-inclusive cross section can be expressed in a factorized form in terms of TMDs only in the kinematic limits $M^{2}\ll Q^{2}$ and $\bm{P}_{hT}^{2}\ll Q^{2}$. In these limits, the structure function $F_{UU,L}$ of Eq. (4) can be neglected Bacchetta et al. (2008a). The structure function $F_{L}$ in the denominator contains contributions involving powers of the strong coupling constant $\alpha_{S}$ at an order that goes beyond the level reached in this analysis; hence, it will be consistently neglected (for measurements and estimates of the $F_{L}$ structure function see, e.g., Refs. Chekanov et al. (2009); Andreev et al. (2014) and references therein).

To express the structure functions in terms of TMD PDFs and FFs, we rely on the factorized formula for SIDIS Collins and Soper (1981); Collins et al. (1985); Ji and Yuan (2002); Ji et al. (2005); Collins (2013); Aybat and Rogers (2011); Echevarria et al. (2012, 2013); Collins and Rogers (2013) (see Fig. 1 for a graphical representation of the involved transverse momenta):

[TABLE]

Here, $\mathcal{H}_{UU,T}$ is the hard scattering part; $f_{1}^{a}(x,\bm{k}_{\perp}^{2};Q^{2})$ is the TMD PDF of unpolarized partons with flavor $a$ in an unpolarized proton, carrying longitudinal momentum fraction $x$ and transverse momentum $\bm{k}_{\perp}$. The D_{1}^{a\mbox{\raisebox{-4.5pt}[0.0pt][0.0pt]{\vec{\phantom{v}}}}h}(z,\bm{P}_{\perp}^{2};Q^{2}) is the TMD FF describing the fragmentation of an unpolarized parton with flavor $a$ into an unpolarized hadron $h$ carrying longitudinal momentum fraction $z$ and transverse momentum $\bm{P}_{\perp}$ (see Fig. 1). TMDs generally depend on two energy scales Collins (2013), which enter via the renormalization of ultraviolet and rapidity divergencies. In this work we choose them to be equal and set them to $Q^{2}$. The term $Y_{UU,T}$ is introduced to ensure a matching to the perturbative fixed-order calculations at higher transverse momenta.

In our analysis, we neglect any correction of the order of $M^{2}/Q^{2}$ or higher to Eq. (6). At large $Q^{2}$ this is well justified. However, fixed-target DIS experiments typically collect a large amount of data at relatively low $Q^{2}$ values, where these assumptions should be all tested in future studies. The reliability of the theoretical description of SIDIS at low $Q^{2}$ has been recently discussed in Refs. Boglione et al. (2017); Moffat et al. (2017).

Eq. (6) can be expanded in powers of $\alpha_{S}$. In the present analysis, we will consider only the terms at order $\alpha_{S}^{0}$. In this case $\mathcal{H}^{a}_{UU,T}(Q^{2})\approx e_{a}^{2}$ and $Y_{UU,T}\approx 0$. However, perturbative corrections include large logarithms L\equiv\log\big{(}z^{2}Q^{2}/P_{hT}^{2}\big{)}, so that $\alpha_{S}L\approx 1$. In the present analysis, we will take into account all leading and Next-to-Leading Logarithms (NLL).111We remark that formulas at NNLL are available in the literature Echevarria et al. (2016).

In these approximations ($\alpha_{S}^{0}$ and NLL), only the first term in Eq. (6) is relevant (often in the literature this has been called $W$ term). We expect this term to provide a good description of the structure function only in the region where $P_{hT}^{2}\ll Q^{2}$. It can happen that $Y_{UU,T}$, defined in the standard way (see, e.g., Ref. Collins et al. (1985)), gives large contributions also in this region, but it is admissible to redefine it in order to avoid this problem Collins et al. (2016). We leave a detailed treatment of the matching to the high $P_{hT}^{2}\approx Q^{2}$ region to future investigations.

To the purpose of applying TMD evolution equations, we need to calculate the Fourier transform of the part of Eq. (6) involving TMDs. The structure function thus reduces to

[TABLE]

where we introduced the Fourier transforms of the TMD PDF and FF according to

[TABLE]

II.2 Drell–Yan and Z production

In a Drell–Yan process, two hadrons $A$ and $B$ with momenta $P_{A}$ and $P_{B}$ collide at a center-of-mass energy squared $s=(P_{A}+P_{B})^{2}$ and produce a virtual photon or a $Z$ boson plus hadrons. The boson decays into a lepton-antilepton pair. The reaction formula is

[TABLE]

The invariant mass of the virtual photon is $Q^{2}=q^{2}$ with $q=l+l^{\prime}$. We introduce the rapidity of the virtual photon/Z boson

[TABLE]

where the $z$ direction is defined along the momentum of hadron A (see Fig. 2).

The cross section can be written in terms of structure functions Boer and Vogelsang (2006); Arnold et al. (2009). For our purposes, we need the unpolarized cross section integrated over $d\Omega$ and over the azimuthal angle of the virtual photon,

[TABLE]

The elementary cross sections are

[TABLE]

where $\theta_{W}$ is Weinberg’s angle, $M_{Z}$ is the mass of the $Z$ boson, and $B_{R}(Z\rightarrow\ell^{+}\ell^{-})$ is the branching ratio for the $Z$ boson decay in two leptons. We adopted the narrow-width approximation, i.e., we neglect contributions for $Q^{2}\neq M_{Z}^{2}$. We used the values $\sin^{2}\theta_{W}=0.2313$, $M_{Z}=91.18$ GeV, and $B_{R}(Z\rightarrow\ell^{+}\ell^{-})=3.366$ Patrignani et al. (2016). Similarly to the SIDIS case, in the kinematic limit $q_{T}^{2}\ll Q^{2}$ the structure function $F_{UU}^{2}$ can be neglected (for measurement and estimates of this structure function see, e.g., Ref. Lambertsen and Vogelsang (2016) and references therein).

The longitudinal momentum fractions of the annihilating quarks can be written in terms of rapidity in the following way

[TABLE]

Some experiments use the variable $x_{F}$, which is connected to the other variables by the following relations

[TABLE]

The structure function $F_{UU}^{1}$ can be written as (see Fig. 2 for a graphical representation of the involved transverse momenta)

[TABLE]

As in the SIDIS case, in our analysis we neglect the $Y_{UU}$ term and we consider the hard coefficients only up to leading order in the couplings, i.e.,

[TABLE]

where222We remind the reader that the value of weak isospin $I_{3}$ is equal to $+1/2$ for $u$, $c$, $t$ and $-1/2$ for $d$, $s$, $b$.

[TABLE]

The structure function can be conveniently expressed as a Fourier transform of the right-hand side of Eq. (16) as

[TABLE]

II.3 TMDs and their evolution

Evolution equations quantitatively describe the connection between different values for the energy scales. In the following we will set their initial values to $\mu_{b}^{2}$ and their final values as $Q^{2}$, so that only $Q^{2}$ and $\mu_{b}^{2}$ need to be specified in a TMD distribution. Following the formalism of Refs. Collins (2013); Aybat and Rogers (2011), the unpolarized TMD distribution and fragmentation functions in configuration space for a parton with flavor $a$ at a certain scale $Q^{2}$ can be written as

[TABLE]

We choose the scale $\mu_{b}$ to be

[TABLE]

where $\gamma_{E}$ is the Euler constant and

[TABLE]

This variable replaces the simple dependence upon $\xi_{T}$ in the perturbative parts of the TMD definitions of Eqs. (20), (21). In fact, at large $\xi_{T}$ these parts are no longer reliable. Therefore, the $\bar{\xi}_{\ast}$ is chosen to saturate on the maximum value $\xi_{\rm max}$, as suggested by the CSS formalism Collins (2013); Aybat and Rogers (2011). On the other hand, at small $\xi_{T}$ the TMD formalism is not valid and should be matched to the fixed-order collinear calculations. The way the matching is implemented is not unique. In any case, the TMD contribution can be arbitrarily modified at small $\xi_{T}$. In our approach, we choose to saturate $\bar{\xi}_{\ast}$ at the minimum value $\xi_{\rm min}\propto 1/Q$. With the appropriate choices, for $\xi_{T}=0$ the Sudakov exponent vanishes, as it should Parisi and Petronzio (1979); Altarelli et al. (1984). Our choice partially corresponds to modifying the resummed logarithms as in Ref. Bozzi et al. (2011) and to other similar modifications proposed in the literature Boer and den Dunnen (2014); Collins et al. (2016). One advantage of these kind of prescriptions is that by integrating over the impact parameter $\xi_{T}$, the collinear expression for the cross section in terms of collinear PDFs is recovered, at least at leading order Collins et al. (2016). We remind the reader that there are different schemes available to deal with the high-$\xi_{T}$ region, such as the the so-called “complex-$\xi$ prescription” Laenen et al. (2000) or an extrapolation of the perturbative small-$\xi_{T}$ calculation to the large $\xi_{T}$ region based on dynamical power corrections Qiu and Zhang (2001).

The values of $\xi_{\rm max}$ and $\xi_{\rm min}$ could be regarded as arbitrary scales separating perturbative from nonperturbative regimes. We choose to fix them to the values

[TABLE]

The motivations are the following:

• •

with the above choices, the scale $\mu_{b}$ is constrained between 1 GeV and $Q$, so that the collinear PDFs are never computed at a scale lower than 1 GeV and the lower limit of the integrals contained in the definition of the perturbative Sudakov factor (see Eq. (30)) can never become larger than the upper limit;

• •

at $Q=Q_{0}=1$ GeV, $\xi_{\rm max}=\xi_{\rm min}$ and there are no evolution effects; the TMD is simply given by the corresponding collinear function multiplied by a nonperturbative contribution depending on $k_{\perp}$ (plus possible corrections of order $\alpha_{S}$ from the Wilson coefficients).

At NLL accuracy, for our choice of scales $K(\bar{\xi}_{\ast},\mu_{b})=0$. Similarly, the $C$ and $\hat{C}$ are perturbatively calculable Wilson coefficients for the TMD distribution and fragmentation functions, respectively. They are convoluted with the corresponding collinear functions according to

[TABLE]

In the present analysis, we consider only the leading-order term in the $\alpha_{S}$ expansion for $C$ and $\hat{C}$ , i.e.,

[TABLE]

As a consequence of the choices we made, the expression for the evolved TMD functions reduces to

[TABLE]

The Sudakov exponent $S$ can be written as

[TABLE]

where the functions $A$ and $B$ have a perturbative expansions of the form

[TABLE]

To NLL accuracy, we need the following terms Davies and Stirling (1984); Collins et al. (1985)

[TABLE]

We use the approximate analytic expression for $\alpha_{S}$ at NLO with the $\Lambda_{\text{QCD}}=$ 340 MeV, 296 MeV, 214 MeV for three, four, five flavors, respectively, corresponding to a value of $\alpha_{S}(M_{Z})=0.117$. We fix the flavor thresholds at $m_{c}=1.5$ GeV and $m_{b}=4.7$ GeV. The integration of the Sudakov exponent in Eq. (30) can be done analytically (for the complete expressions see, e.g., Refs. Frixione et al. (1999); Bozzi et al. (2006); Echevarria et al. (2013)).

Following Refs. Nadolsky et al. (2000); Landry et al. (2003); Konychev and Nadolsky (2006), for the nonperturbative Sudakov factor we make the traditional choice

[TABLE]

with $g_{2}$ a free parameter. Recently, several alternative forms have been proposed Aidala et al. (2014); Collins and Rogers (2015). Also, recent theoretical studies aimed at calculating this term using nonperturbative methods Scimemi and Vladimirov (2017). All these choices should be tested in future studies. In Ref. D’Alesio et al. (2014a), a good agreement with data was achieved even without this term, but this is not possible when including data at low $Q^{2}$.

In this analysis, for the collinear PDFs $f_{1}^{a}$ we adopt the GJR08FFnloE set Gluck et al. (2008) through the LHAPDF library Buckley et al. (2015), and for the collinear fragmentation functions the DSS14 NLO set for pions de Florian et al. (2015) and the DSS07 NLO set for kaons de Florian et al. (2007).333After the completion of our analysis, a new set of kaon fragmentation function was presented in Ref. de Florian et al. (2017). We will comment on the use of other PDF sets in Sec. IV.3.

We parametrize the intrinsic nonperturbative parts of the TMDs in the following ways

[TABLE]

After performing the anti-Fourier transform, the $f_{1\rm NP}$ and $D_{1\rm NP}$ in momentum space correspond to

[TABLE]

The TMD PDF at the starting scale is therefore a normalized sum of a Gaussian with variance $g_{1}$ and the same Gaussian weighted by a factor $\lambda\bm{k}_{\perp}^{2}$. The TMD FF at the starting scale is a normalized sum of a Gaussian with variance $g_{3}$ and a second Gaussian with variance $g_{4}$ weighted by a factor $\lambda_{F}\bm{P}_{\perp}^{2}/z^{2}$. The choice of this particular functional forms is motivated by model calculations: the weighted Gaussian in the TMD PDF could arise from the presence of components of the quark wave function with angular momentum $L=1$ Bacchetta et al. (2008b); Pasquini et al. (2008); Avakian et al. (2010); Bacchetta et al. (2010); Burkardt and Pasquini (2016). Similar features occur in models of fragmentation functions Bacchetta et al. (2002); Bacchetta et al. (2008b); Matevosyan et al. (2012).

The Gaussian width of the TMD distributions may depend on the parton flavor $a$ Signori et al. (2013); Matevosyan et al. (2012); Schweitzer et al. (2013). In the present analysis, however, we assume they are flavor independent. The justification for this choice is that most of the data we are considering are not sufficiently sensitive to flavor differences, leading to unclear results. We will devote attention to this issue in further studies.

Finally, we assume that the Gaussian width of the TMD depends on the fractional longitudinal momentum $x$ according to

[TABLE]

where $\alpha,\,\sigma,$ and $N_{1}\equiv g_{1}(\hat{x})$ with $\hat{x}=0.1$, are free parameters. Similarly, for fragmentation functions we have

[TABLE]

where $\beta,\,\gamma,\,\delta,$ and $N_{3,4}\equiv g_{3,4}(\hat{z})$ with $\hat{z}=0.5$ are free parameters.

The average transverse momentum squared for the distributions in Eq. (36) and (37) can be computed analytically:

[TABLE]

III Data analysis

The main goals of our work are to extract information about intrinsic transverse momenta, to study the evolution of TMD parton distributions and fragmentation functions over a large enough range of energy, and to test their universality among different processes. To achieve this we included measurements taken from SIDIS, Drell–Yan and $Z$ boson production from different experimental collaborations at different energy scales. In this chapter we describe the data sets considered for each process and the applied kinematic cuts.

Tab. 1 refers to the data sets for SIDIS off proton target (Hermes experiment) and presents their kinematic ranges. The same holds for Tab. 2, Tab. 3, Tab. 4 for SIDIS off deuteron (Hermes and Compass experiments), Drell–Yan events at low energy and $Z$ boson production respectively. If not specified otherwise, the theoretical formulas are computed at the average values of the kinematic variables in each bin.

III.1 Semi-inclusive DIS data

The SIDIS data are taken from Hermes Airapetian et al. (2013) and Compass Adolph et al. (2013) experiments. Both data sets have already been analyzed in previous works, e.g., Refs. Signori et al. (2013); Anselmino et al. (2014), however they have never been fitted together, including also the contributions deriving from TMD evolution.

The application of the TMD formalism to SIDIS depends on the capability of identifying the current fragmentation region. This task has been recently discussed in Ref. Boglione et al. (2017), where the authors point out a possible overlap among different fragmentation regions when the hard scale $Q$ is sufficiently low. In this paper we do not tackle this problem and we leave it to future studies. As described in Tabs. 1 and 2, we identify the current fragmentation region operating a cut on $z$ only, namely $0.2<z<0.74$.

Another requirement for the applicability of TMD factorization is the presence of two separate scales in the process. In SIDIS, those are the $Q^{2}$ and $P_{hT}^{2}$, which should satisfy the condition $P_{hT}^{2}\ll Q^{2}$, or more precisely $P_{hT}^{2}/z^{2}\ll Q^{2}$. We implement this condition by imposing $P_{hT}<\min[0.2\ Q,0.7\ Qz]+0.5$ GeV. With this choice, $P_{hT}^{2}$ is always smaller than $Q^{2}/3$, but in a few bins (at low $Q^{2}$ and $z$) $P_{hT}^{2}/z^{2}$ may become larger than $Q^{2}$. The applicability of TMD factorization in this case could be questioned. However, as we will explain further in Sec. IV.3, we can obtain a fit that can describe a wide region of $P_{hT}$ and can also perform very well in a restricted region, where TMD factorization certainly holds.

All these choices are summarized in Tabs. 1 and 2.

III.1.1 Hermes data

Hermes hadron multiplicities are measured in a fixed target experiment, colliding a $27.6$ GeV lepton beam on a hydrogen ($p$) or deuterium ($D$) gas target, for a total of 2688 points. These are grouped in bins of $(x,z,Q^{2},P_{hT})$ with the average values of $(x,Q^{2})$ ranging from about $(0.04,1.25\text{ GeV}^{2})$ to $(0.4,9.2\text{ GeV}^{2})$. The collinear energy fraction $z$ in Eq. (2) ranges in $0.1\leq z\leq 0.9$. The transverse momentum of the detected hadron satisfies $0.1\text{ GeV}\leq|P_{hT}|\leq 1.3\text{ GeV}$. The peculiarity of Hermes SIDIS experiment lies in the ability of its detector to distinguish between pions and kaons in the final state, in addition to determining their momenta and charges. We consider eight different combinations of target ($p,\,D$) and detected charged hadron ($\pi^{\pm},\,K^{\pm}$ ). The Hermes collaboration published two distinct sets, characterized by the inclusion or subtraction of the vector meson contribution. In our work we considered only the data set where this contribution has been subtracted.

III.1.2 Compass data

The Compass collaboration extracted multiplicities for charge-separated but unidentified hadrons produced in SIDIS off a deuteron (${}^{6}\text{LiD}$) target Adolph et al. (2013). The number of data points is an order of magnitude higher compared to the Hermes experiment. The data are organized in multidimensional bins of $(x,z,Q^{2},P_{hT})$, they cover a range in $(x,Q^{2})$ from about $(0.005,1.11\text{ GeV}^{2})$ to $(0.09,7.57\text{ GeV}^{2})$ and the interval $0.2\leq z\leq 0.8$. The multiplicities published by Compass are affected by normalization errors (see the erratum to Ref. Adolph et al. (2013)). In order to avoid this issue, we divide the data in each bin in ($x,z,Q^{2}$) by the data point with the lowest $P_{hT}^{2}$ in the bin. As a result, we define the normalized multiplicity as

[TABLE]

where the multiplicity $m_{N}^{h}$ is defined in Eq. (3). When fitting normalized multiplicities, the first data point of each bin is considered as a fixed constraint and excluded from the degrees of freedom.

III.2 Low-energy Drell–Yan data

We analyze Drell–Yan events collected by fixed-target experiments at low-energy. These data sets have been considered also in previous works, e.g., in Ref. Landry et al. (2001, 2003); Konychev and Nadolsky (2006); D’Alesio et al. (2014b). We used data sets from the E288 experiment Ito et al. (1981), which measured the invariant dimuon cross section $Ed^{3}\sigma/dq^{3}$ for the production of $\mu^{+}\mu^{-}$ pairs from the collision of a proton beam with a fixed target, either composed of Cu or Pt. The measurements were performed using proton incident energies of $200$, $300$ and $400$ GeV, producing three different data sets. Their respective center of mass energies are $\sqrt{s}=19.4,23.8,27.4$ GeV. We also included the set of measurements $Ed^{3}\sigma/dq^{3}$ from E605 Moreno et al. (1991), extracted from the collision of a proton beam with an energy of $800$ GeV ($\sqrt{s}=38.8$ GeV) on a copper fixed target .

The explored $Q$ values are higher compared to the SIDIS case, as can be seen in Tab. 3. E288 provides data at fixed rapidity, whereas E605 provides data at fixed $x_{F}=0.1$. We can apply TMD factorization if $q_{T}^{2}\ll Q^{2}$, where $q_{T}$ is the transverse momentum of the intermediate electroweak boson, reconstructed from the kinematics of the final state leptons. We choose $q_{T}<0.2\ Q+0.5$ GeV. As suggested in Ref. Ito et al. (1981), we consider the target nuclei as an incoherent ensemble composed $40\%$ by protons and $60\%$ by neutrons.

As we already observed, results from E288 and E605 experiments are reported as $\frac{Ed^{3}\sigma}{d^{3}q}$; this variable is related to the differential cross section of Eq. (12) in the following way:

[TABLE]

where $\phi$ is the polar angle of $q_{T}$ and the third term is the average over $\phi$. Therefore, the invariant dimuon cross section can be obtained from Eq. (12) integrating over $Q^{2}$ and adding a factor $1/\pi$ to the result

[TABLE]

Numerically we checked that integrating in $Q^{2}$ only the prefactor $\sigma_{q}^{\gamma}$ (see Eq. (13)) introduces only a negligible error in the theoretical estimates. We also assume that $\alpha_{\rm em}$ does not change within the experimental bin. Therefore, for Drell–Yan we obtain

[TABLE]

where $Q_{i,f}$ are the lower and upper values in the experimental bin.

III.3 Z-boson production data

In order to reach higher $Q$ and $q_{T}$ values, we also consider $Z$ boson production in collider experiments at Tevatron. We analyze data from CDF and D0, collected during Tevatron Run I Affolder et al. (2000); Abbott et al. (2000) at $\sqrt{s}=1.8\text{ TeV}$ and Run II Aaltonen et al. (2012); Abazov et al. (2008) at $\sqrt{s}=1.96\text{ TeV}$. CDF and D0 collaborations studied the differential cross section for the production of an $e^{+}e^{-}$ pair from $p\bar{p}$ collision through an intermediate $Z$ vector boson, namely $p\bar{p}\rightarrow Z\rightarrow e^{+}e^{-}+X$.

The invariant mass distribution peaks at the $Z$-pole, $Q\approx M_{Z}$, while the transverse momentum of the exchanged $Z$ ranges in $0<q_{T}<20\text{ GeV}$. We use the same kinematic cut applied to Drell–Yan events: $q_{T}<0.2\ Q+0.5$ GeV $=18.7$ GeV, since $Q$ is fixed to $M_{Z}$.

The observable measured in CDF and D0 is

[TABLE]

apart from the case of D0 Run II, for which the published data refer to $1/\sigma\times d\sigma/dq_{T}$. In order to work with the same observable, we multiply the D0 Run II data by the total cross section of the process $\sigma_{exp}=255.8\pm 16\text{ pb}$ Abulencia et al. (2007). In this case, we add in quadrature the uncertainties of the total cross section and of the published data.

We normalize our functional form with the factors listed in Tab. 4. These are the same normalization factors used in Ref. D’Alesio et al. (2014b), computed by comparing the experimental total cross section with the theoretical results based on the code of Ref. Catani et al. (2009). These factors are not precisely consistent with our formulas. In fact, as we will discuss in Sec. IV.3 a 5% increase in these factors would improve the agreement with data, without affecting the TMD parameters.

III.4 The replica method

Our fit is based on the replica method. In this section we describe it and we give a definition of the $\chi^{2}$ function minimized by the fit procedure. The fit and the error analysis are carried out using a similar Monte Carlo approach as in Refs. Bacchetta et al. (2013); Signori et al. (2013); Radici et al. (2015) and taking inspiration from the work of the Neural-Network PDF (NNPDF) collaboration (see, e.g., Refs. Forte et al. (2002); Ball et al. (2009, 2010)). The approach consists in creating $\mathcal{M}$ replicas of the data points. In each replica (denoted by the index $r$), each data point $i$ is shifted by a Gaussian noise with the same variance as the measurement. Each replica, therefore, represents a possible outcome of an independent experimental measurement, which we denote by $m_{N,r}^{h}(x,z,\bm{P}_{hT}^{2},Q^{2})$. The number of replicas is chosen so that the mean and standard deviation of the set of replicas accurately reproduces the original data points. In this case 200 replicas are sufficient for the purpose. The error for each replica is taken to be equal to the error on the original data points. This is consistent with the fact that the variance of the $\mathcal{M}$ replicas should reproduce the variance of the original data points.

A minimization procedure is applied to each replica separately, by minimizing the following error function:

[TABLE]

The sum runs over the $i$ experimental points, including all species of targets $N$ and final-state hadrons $h$. In each $z$ bin for each replica the values of the collinear fragmentation functions D_{1}^{a\mbox{\raisebox{-4.5pt}[0.0pt][0.0pt]{\vec{\phantom{v}}}}h} are independently modified with a Gaussian noise with standard deviation equal to the theoretical error \Delta D_{1}^{a\mbox{\raisebox{-4.5pt}[0.0pt][0.0pt]{\vec{\phantom{v}}}}h}. In this work we rely on different parametrizations for D_{1}^{a\mbox{\raisebox{-4.5pt}[0.0pt][0.0pt]{\vec{\phantom{v}}}}h}: for pions we use the DSEHS analysis de Florian et al. (2015) at NLO in $\alpha_{S}$; for kaons we use the DSS parametrization de Florian et al. (2007) at LO in $\alpha_{S}$. The uncertainties \Delta D_{1}^{a\mbox{\raisebox{-4.5pt}[0.0pt][0.0pt]{\vec{\phantom{v}}}}h} are estimated from the plots in Ref. Epele et al. (2012); they represents the only source of uncertainty in $\Delta m_{N,{\rm theo}}^{h}$. Statistical and systematic experimental uncertainties $\Delta m_{N,\mbox{\tiny stat}}^{h}$ and $\Delta m_{N,\mbox{\tiny sys}}^{h}$ are taken from the experimental collaborations. We do not take into account the covariance among different kinematic bins.

We minimize the error function in Eq. (46) with Minuit James and Roos (1975). In each replica we randomize the starting point of the minimization, to better sample the space of fit parameters. The final outcome is a set of $\mathcal{M}$ different vectors of best-fit parameters, $\{p_{0r}\},\;r=1,\ldots\mathcal{M}$, with which we can calculate any observable, its mean, and its standard deviation. The distribution of these values needs not to be necessarily Gaussian. In fact, in this case the $1\sigma$ confidence interval is different from the 68% interval. The latter can simply be computed for each experimental point by rejecting the largest and the lowest 16% of the $\mathcal{M}$ values.

Although the minimization is performed on the function defined in Eq. (46), the agreement of the $\mathcal{M}$ replicas with the original data is expressed in terms of a $\chi^{2}$ function defined as in Eq. (46) but with the replacement $m_{N,r}^{h}\to m_{N}^{h}$, i.e., with respect to the original data set. If the model is able to give a good description of the data, the distribution of the $\mathcal{M}$ values of $\chi^{2}$/d.o.f. should be peaked around one.

IV Results

Our work aims at simultaneously fitting for the first time data sets related to different experiments. In the past, only fits related either to SIDIS or hadronic collisions have been presented. Here we mention a selection of recent existing analyses.

In Ref. Signori et al. (2013), the authors fitted Hermes multiplicities only (taking into account a total of 1538 points) without taking into account QCD evolution. In that work, a flavor decomposition in transverse momentum of the unpolarized TMDs and an analysis of the kinematic dependence of the intrinsic average square transverse momenta were presented. In Ref. Anselmino et al. (2014) the authors fitted Hermes and Compass multiplicities separately (576 and 6284 points respectively), without TMD evolution and introducing an ad-hoc normalization for Compass data. A fit of SIDIS data including TMD evolution was performed on measurements by the H1 collaboration of the so-called transverse energy flow Nadolsky et al. (2000); Aid et al. (1995).

Looking at data from hadronic collisions, Konychev and Nadolsky Konychev and Nadolsky (2006) fitted data of low-energy Drell–Yan events and $Z$-boson production at Tevatron, taking into account TMD evolution at NLL accuracy (this is the most recent of a series of important papers on the subject Ladinsky and Yuan (1994); Landry et al. (2001, 2003)). They fitted in total 98 points. Contrary to our approach, Konychev and Nadolsky studied the quality of the fit as a function of $\xi_{\text{max}}$. They found that the best value for $\xi_{\text{max}}$ is $1.5$ GeV*-1* (to be compared to our choice $\xi_{\text{max}}\approx 1.123$ GeV*-1*, see Sec. II.3). Comparisons of best-fit values in the nonperturbative Sudakov form factors are delicate, since the functional form is different from ours. In 2014 D’Alesio, Echevarria, Melis, Scimemi performed a fit D’Alesio et al. (2014b) of Drell–Yan data and $Z$-boson production data at Tevatron, focusing in particular on the role of the nonperturbative contribution to the kernel of TMD evolution. This is the fit with the highest accuracy in TMD evolution performed up to date (NNLL in the Sudakov exponent and ${\cal O}(\alpha_{S})$ in the Wilson coefficients). In the same year Echevarria, Idilbi, Kang and Vitev Echevarria et al. (2014) presented a parametrization of the unpolarized TMD that described qualitatively well some bins of Hermes and Compass data, together with Drell–Yan and Z-production data. A similar result was presented by Sun, Isaacson, Yuan and Yuan Sun et al. (2014).

In the following, we detail the results of a fit to the data sets described in Sec. III with a flavor-independent configuration for the transverse momentum dependence of unpolarized TMDs. In Tab. 5 we present the total $\chi^{2}$. The number of degrees of freedom (d.o.f.) is given by the number of data points analyzed reduced by the number of free parameters in the error function. The overall quality of the fit is good, with a global $\chi^{2}$/d.o.f. $=1.55\pm 0.05$. Uncertainties are computed as the $68\%$ confidence level (C.L.) from the replica methodology.

IV.1 Agreement between data and theory

The partition of the global $\chi^{2}$ among SIDIS off a proton, SIDIS off a deuteron, Drell–Yan and $Z$ production events is given in Tab. 6, 7, 8, 9 respectively.

Semi-inclusive DIS

For SIDIS at Hermes off a proton, most of the contribution to the $\chi^{2}$ comes from events with a $\pi^{+}$ in the final state. In Ref. Signori et al. (2013) the high $\chi^{2}$ was attributed to the poor agreement between experiment and theory at the level of the collinear multiplicities. In this work we use a newer parametrization of the collinear FFs (DSEHS de Florian et al. (2015)), based on a fit which includes Hermes collinear pion multiplicities. In spite of this improvement, the contribution to $\chi^{2}$ from Hermes data is higher then in Ref. Signori et al. (2013), because the present fit includes data from other experiments (Hermes represents less than 20% of the whole data set). The bins with the worst agreement are at low $Q^{2}$. As we will discuss in Sec. IV.3, we think that the main reason for the large $\chi^{2}$ at Hermes is a normalization difference. This may also be due to the fact that we are computing our theoretical estimates at the average values of the kinematic variables, instead of integrating the multiplicities in each bin. Kaon multiplicities have in general a lower $\chi^{2}$, due to the bigger statistical errors and the large uncertainties for the kaon FFs.

For pion production off a deuteron at Hermes the $\chi^{2}$ is lower with respect to the production off a proton, but still compatible within uncertainties. For kaon production off a deuteron the $\chi^{2}$ is higher with respect to the scattering off a proton. The difference is especially large for $K^{-}$.

SIDIS at Compass involves scattering off deuteron only, $D\to h^{\pm}$, and we identify $h\equiv\pi$. The quality of the agreement between theory and Compass data is better than in the case of pion production at Hermes. This depends on at least two factors: first, our fit is essentially driven by the Compass data, which represent about 75% of the whole data set; second, the observable that we fit in this case is the normalized multiplicity, defined in Eq. (41). This automatically eliminates most of the discrepancy between theory and data due to normalization.

Fig. 3 presents the agreement between the theoretical formula in (3) and the Hermes multiplicities for production of pions off a proton and a deuteron. Different $\langle x\rangle$, $\langle z\rangle$ and $\langle Q^{2}\rangle$ bins are displayed as a function of the transverse momentum of the detected hadron $P_{hT}$. The grey bands are an envelope of the $200$ replica of best-fit curves. For every point in $P_{hT}$ we apply a $68\%$ C.L. selection criterion. Points marked with different symbols and colors correspond to different $\langle z\rangle$ values. There is a strong correlation between $\langle x\rangle$ and $\langle Q^{2}\rangle$ that does not allow us to explore the $x$ and $Q^{2}$ dependence of the TMDs separately. Studying the contributions to the $\chi^{2}$/points as a function of the kinematics, we notice that the $\chi^{2}(Q^{2})$ tends to improve as we move to higher $Q^{2}$ values, where the kinematic approximations of factorization are more reliable. Moreover, usually the $\chi^{2}(z)$ increases at lower $z$ values.

Fig. 4 has same contents and notation as in Fig. 3 but for kaons in the final state. In this case, the trend of the agreement as a function of $Q^{2}$ is not as clear as for the case of pions: good agreement is found also at low $Q^{2}$.

In Fig. 5 we present Compass normalized multiplicities (see Eq. (41)) for production of $\pi^{-}$ off a deuteron for different $\langle x\rangle$, $\langle z\rangle$, and $\langle Q^{2}\rangle$ bins as a function of the transverse momentum of the detected hadron $P_{hT}$. The open marker around the first $P_{hT}$ point in each panel indicates that the first value is fixed and not fitted. The correlation between $x$ and $Q^{2}$ is less strong than at Hermes and this allows us to study different $\langle x\rangle$ bins at fixed $\langle Q^{2}\rangle$. For the highest $Q^{2}$ bins, the agreement is good for all $\langle x\rangle$, $\langle z\rangle$ and $P_{hT}^{2}$. In bins at lower $Q^{2}$, the descriptions gets worse, especially at low and high $z$. For fixed $\langle Q^{2}\rangle$ and high $\langle z\rangle$, a good agreement is recovered moving to higher $\langle x\rangle$ bins (see, e.g., the third line from the top in Fig. 5).

Fig. 6 has same content and notation as in Fig. 5, but for $h^{+}\equiv\pi^{+}$. The same comments on the agreement between theory and the data apply.

Drell–Yan and $Z$ production

The low energy Drell–Yan data collected by the E288 and E605 experiments at Fermilab have large error bands (see Fig. 7). This is why the $\chi^{2}$ values in Tab. 8 are rather low compared to the other data sets.

The agreement is also good for $Z$ boson production, see Tab. 9. The statistics from Run-II is higher, which generates smaller experimental uncertainties and higher $\chi^{2}$, especially for the CDF experiment.

Fig. 7 displays the cross section for DY events differential with respect to the transverse momentum $q_{T}$ of the virtual photon, its invariant mass $Q^{2}$ and rapidity $y$. As for the case of SIDIS, the grey bands are the $68\%$ C.L. envelope of the 200 replicas of the fit function. The four panels represents different values for the rapidity $y$ or $x_{F}$ (see Eq. (15)). In each panel, we have plots for different $Q^{2}$ values. The lower is $Q$, the less points in $q_{T}$ we fit (see also Sec. III.2). The hard scale lies in the region $4.5<\langle Q\rangle<13.5$ GeV. This region is of particular importance, since these “moderate” $Q$ values should be high enough to safely apply factorization and, at the same time, low enough in order for the nonperturbative effects to not be shaded by transverse momentum resummation.

In Fig. 8 we compare the cross section differential with respect to the transverse momentum $q_{T}$ of the virtual $Z$ (namely Eq. (12) integrated over $\eta$) with data from CDF and D0 at Tevatron Run I and II. Due to the higher $Q=M_{Z}$, the range explored in $q_{T}$ is much larger compared to all the other observables considered. The tails of the distributions deviate from a Gaussian behavior, as it is also evident in the bins at higher $Q^{2}$ in Fig. 7. The band from the replica methodology in this case is much narrower, due to the reduced sensitivity to the intrinsic transverse momenta at $Q=M_{Z}$ and to the limited range of best-fit values for the parameter $g_{2}$, which controls soft-gluon emission. As an effect of TMD evolution, the peak shifts from $\sim 1$ GeV for Drell–Yan events in Fig. 7 to $\sim 5$ GeV in Fig. 8. The position of the peak is affected both by the perturbative and the nonperturbative part of the Sudakov exponent (see Sec. II.3 and Signori (2016a)). Most of the contributions to the $\chi^{2}$ comes from normalization effects and not from the shape in $q_{T}$ (see Sec. IV.3).

IV.2 Transverse momentum dependence at 1 GeV

The variables $\xi_{\rm min}$ and $\xi_{\rm max}$ delimit the range in $\xi_{T}$ where transverse momentum resummation is computed perturbatively. The $g_{2}$ parameter enters the nonperturbative Sudakov exponent and quantifies the amount of transverse momentum due to soft gluon radiation that is not included in the perturbative part of the Sudakov form factor. As already explained in Sec. II.3, in this work we fix the value for $\xi_{\rm min}$ and $\xi_{\rm max}$ in such a way that at $Q=1$ GeV the unpolarized TMDs coincide with their nonperturbative input. We leave $g_{2}$ as a fit parameter.

Tab. 10 summarizes the chosen values of $\xi_{\rm min}$, $\xi_{\rm max}$ and the best-fit value for $g_{2}$. The latter is given as an average with $68\%$ C.L. uncertainty computed over the set of 200 replicas. We also quote the results obtained from replica 105, since its parameters are very close to the mean values of all replicas. We obtain a value $g_{2}=0.13\pm 0.01$, smaller than the value ($g_{2}=0.184\pm 0.018$) obtained in Ref. Konychev and Nadolsky (2006), where however no SIDIS data was taken into consideration, and smaller than the value ($g_{2}=0.16$) chosen in Ref. Echevarria et al. (2014). We stress however that our prescriptions involving both $\xi_{\text{min}}$ and $\xi_{\text{max}}$ are different from previous works.

Tab. 11 collects the best-fit values of parameters in the nonperturbative part of the TMDs at $Q=1$ GeV (see Eqs. (34) and (35)); as for $g_{2}$, we give the average value over the full set of replicas and the standard deviation based on a $68\%$ C.L. (see Sec. III.4), and we also quote the value of replica 105.

In Fig. 9 we compare different extractions of partonic transverse momenta. The horizontal axis shows the value of the average transverse momentum squared for the incoming parton, \big{\langle}\bm{k}_{\perp}^{2}\big{\rangle}(x=0.1) (see Eq. (40)). The vertical axis shows the value of \big{\langle}\bm{P}_{\perp}^{2}\big{\rangle}(z=0.5), the average transverse momentum squared acquired during the fragmentation process (see Eq. (40)). The white square (label 1) indicates the average values of the two quantities obtained in the present analysis at $Q^{2}=1$ GeV2. Each black dot around the white square is an outcome of one replica. The red region around the white square contains the $68\%$ of the replicas that are closest to the average value. The same applies to the white circle and the orange region around it (label 2), related to the flavor-independent version of the analysis in Ref. Signori et al. (2013), obtained by fitting only Hermes SIDIS data at an average $\langle Q^{2}\rangle=2.4$ GeV2 and neglecting QCD evolution. A strong anticorrelation between the transverse momenta is evident in this older analysis. In our new analysis, the inclusion of Drell–Yan and $Z$ production data adds physical information about TMD PDFs, free from the influence of TMD FFs. This reduces significantly the correlation between \big{\langle}\bm{k}_{\perp}^{2}\big{\rangle}(x=0.1) and \big{\langle}\bm{P}_{\perp}^{2}\big{\rangle}(z=0.5). The $68\%$ confidence region is smaller than in the older analysis. The average values of \big{\langle}\bm{k}_{\perp}^{2}\big{\rangle}(x=0.1) are similar and compatible within error bands. The values of \big{\langle}\bm{P}_{\perp}^{2}\big{\rangle}(z=0.5) in the present analysis turn out to be larger than in the older analysis, an effect that is due mainly to Compass data. It must be kept in mind that the two analyses lead also to differences in the $x$ and $z$ dependence of the transverse momentum squared. This dependence is shown in Fig. 10 (a) for \big{\langle}\bm{k}_{\perp}^{2}\big{\rangle}(x) and Fig. 10 (b) for \big{\langle}\bm{P}_{\perp}^{2}\big{\rangle}(z). The bands are computed as the $68\%$ C.L. envelope of the full sets of curves from the 200 replicas. Comparison with other extractions are presented and the legend is detailed in the caption of Fig. 9.

IV.3 Stability of our results

In this subsection we discuss the effect of modifying some of the choices we made in our default fit. Instead of repeating the fitting procedure with different choices, we limit ourselves to checking how the $\chi^{2}$ of a single replica is affected by the modifications.

As starting point we choose replica 105, which, as discussed above, is one of the most representative among the whole replica set. The global $\chi^{2}/$d.o.f. of replica 105 is 1.51. We keep all parameters fixed, without performing any new minimization, and we compute the $\chi^{2}/$d.o.f. after the modifications described in the following.

First of all, we analyze Hermes data with the same strategy as Compass, i.e., we normalize Hermes data to the value of the first bin in $P_{hT}$. In this case, the global $\chi^{2}/$d.o.f. reduces sharply to 1.27. The partial $\chi^{2}$ for the different SIDIS processes measured at Hermes are shown in Table 12. This confirms that normalization effects are the main contribution to the $\chi^{2}$ of SIDIS data and have minor effects on TMD-related parameters. In fact, even if we perform a new fit with this modification, the $\chi^{2}$ does not improve significantly and parameters do not change much.

We consider the effect of changing the normalization of the $Z$-boson data: if we increase the normalization factors quoted in the last row of Tab. 4 by 5%, the $\chi^{2}$ quoted in the last row of Tab. 9 drops to 0.66, 0.52, 0.65, 0.68. This effect is also already visible by eye in Fig. 8: the theoretical curves are systematically below the experimental data points, but the shape is reproduced very well.

We consider the sensitivity of our results to the parameterizations adopted for the collinear quark PDFs. The $\chi^{2}/\text{d.o.f.}$ varies from its original value 1.51, obtained with the NLO GJR 2008 parametrization Gluck et al. (2008), to 1.84 using NLO MSTW 2008 Martin et al. (2009), and 1.85 using NLO CJ12 Owens et al. (2013). In both cases, the agreement with Hermes and $Z$ boson data is not affected significanlty, the agreement with Compass data becomes slightly worse, and the agreement with DY data becomes clearly worse.

An extremely important point is the choice of kinematic cuts. Our default choices are listed in Tabs. 1–4. We consider also more stringent kinematic cuts on SIDIS data: $Q^{2}>1.5$ GeV2 and $0.25<z<0.6$ instead of $Q^{2}>1.4$ GeV2 and $0.2<z<0.7$, leaving the other ones unchanged. The number of bins with these cuts reduces from 8059 to 5679 and the $\chi^{2}/\text{d.o.f.}$ decreases to the value 1.23. In addition, if we replace the constraint $P_{hT}<\text{Min}[0.2\,Q,0.7\,Qz]+0.5$ GeV with $P_{hT}<\text{Min}[0.2\,Q,0.5\,Qz]+0.3$ GeV, the number of bins reduces to 3380 and the $\chi^{2}/$d.o.f. decreases further to 0.96. By adopting the even stricter cut $P_{hT}<0.2\,Qz$, the number of bins drops to only 477, with a $\chi^{2}$/d.o.f. =1.02. We can conclude that our fit, obtained by fitting data in an extended kinematic region, where TMD factorization may be questioned, works extremely well also in a narrower region, where TMD factorization is expected to be under control.

V Conclusions

In this work we demonstrated for the first time that it is possible to perform a simultaneous fit of unpolarized TMD PDFs and FFs to data of SIDIS, Drell–Yan and $Z$ boson production at small transverse momentum collected by different experiments. This constitutes the first attempt towards a global fit of $f_{1}^{a}(x,k_{\perp}^{2})$ and $D_{1}^{a\to h}(z,P_{\perp}^{2})$ in the context of TMD factorization and with the implementation of TMD evolution at NLL accuracy.

We extracted unpolarized TMDs using 8059 data points with 11 free parameters using a replica methodology. We selected data with $Q^{2}>1.4$ GeV2 and $0.2<z<0.7$. We restricted our fit to the small transverse momentum region, selecting the maximum value of transverse momentum on the basis of phenomenological considerations (see Sec. III). With these choices, we included regions where TMD factorization could be questioned, but we checked that our results describe very well the regions where TMD factorization is supposed to hold. The average $\chi^{2}$/d.o.f. is $1.55\pm 0.05$ and can be improved up to 1.02 restricting the kinematic cuts, without changing the parameters (see Sec. IV.3). Most of the discrepancies between experimental data and theory comes from the normalization and not from the transverse momentum shape.

Our fit is performed assuming that the intrinsic transverse momentum dependence of TMD PDFs and FFs can be parametrized by a normalized linear combination of a Gaussian and a weighted Gaussian. We considered that the widths of the Gaussians depend on the longitudinal momenta. We neglected a possible flavor dependence. For the nonperturbative component of TMD evolution, we adopted the choice most often used in the literature (see Sec. II.3).

We plan to release grids of the parametrizations studied in this work via TMDlib Hautmann et al. (2014) to facilitate phenomenological studies for present and future experiments.

In future studies, different functional forms for all the nonperturbative ingredients should be explored, including also a possible flavor dependence of the intrinsic transverse momenta. A more precise analysis from the perturbative point of view is also needed, which should in principle make it possible to relax the tension in the normalization and to describe data at higher transverse momenta. Moreover, the description at low transverse momentum should be properly matched to the collinear fixed-order calculations at high transverse momentum.

Together with an improved theoretical framework, in order to better understand the formalism more experimental data is needed. It would be particularly useful to extend the coverage in $x$, $z$, rapidity, and $Q^{2}$. The $12$ GeV physics program at Jefferson Lab Dudek et al. (2012) will be very important to constrain TMD distributions at large $x$. Additional data from SIDIS (at Compass, at a future Electron-Ion Collider), Drell–Yan (at Compass, at Fermilab), $Z/W$ production (at LHC, RHIC, and at A Fixed-Target Experiment at the LHC Brodsky et al. (2013)) will be very important. Measurements related to unpolarized TMD FFs at $e^{+}e^{-}$ colliders (at Belle-II, BES-III, at a future International Linear Collider) will be invaluable, since they are presently missing.

Our work focused on quark TMDs. We remark that at present almost nothing is known experimentally about gluon TMDs Mulders and Rodrigues (2001); Echevarria et al. (2015), because they typically require higher-energy scattering processes and they are harder to isolate as compared to quark distributions. Several promising measurements have been proposed in order to extract both the unpolarized and linearly polarized gluon TMDs inside an unpolarized proton. The cleanest possibility would be to look at dijet and heavy quark pair production in electron-proton collisions at a future EIC Boer et al. (2011b); Pisano et al. (2013). Other proposals include isolated photon-pair production at RHIC Qiu et al. (2011) and quarkonium production at the LHC Boer and Pisano (2012); den Dunnen et al. (2014); Signori (2016b); Lansberg et al. (2017).

Testing the formalism of TMD factorization and understanding the structure of unpolarized TMDs is only the first crucial step in the exploration of the 3D proton structure in momentum space and this work opens the way to global determinations of TMDs. Building on this, we can proceed to deepen our understanding of hadron structure via polarized structure function and asymmetries (see, e.g., Refs. Aschenauer et al. (2016); Boglione and Prokudin (2016) and references therein) and, at the same time, to test the impact of hadron structure in precision measurements at high-energies, such as at the LHC. A detailed mapping of hadron structure is essential to interpret data from hadronic collisions, which are among the most powerful tools to look for footprints of new physics.

Acknowledgements.

Discussions with Giuseppe Bozzi are gratefully acknowledged. This work is supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 647981, 3DSPIN). AS acknowledges support from U.S. Department of Energy contract DE-AC05-06OR23177, under which Jefferson Science Associates, LLC, manages and operates Jefferson Lab. The work of AS has been funded partly also by the program of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO).