Analysis of vector boson production within TMD factorization

TL;DR

This paper develops a new scheme for TMD scale setting, performs a global fit of unpolarized TMD PDFs using NNLO calculations, and evaluates various nonperturbative models with implications for high-energy collider data analysis.

Contribution

It introduces the $ zeta$-prescription for TMD scale fixation, enabling more accurate extraction of TMD PDFs from experimental data.

Findings

The $ zeta$-prescription reduces perturbative logarithm effects.

Most nonperturbative models are disfavored based on fit quality.

The arTeMiDe code facilitates TMD phenomenology.

Abstract

We present a comprehensive analysis and extraction of the unpolarized transverse momentum dependent (TMD) parton distribution functions, which are fundamental constituents of the TMD factorization theorem. We provide a general review of the theory of TMD distributions, and present a new scheme of scale fixation. This scheme, called the $\zeta$-prescription, allows to minimize the impact of perturbative logarithms in the large range of scales and does not generate undesired power corrections. Within $\zeta$-prescription we consistently include the perturbatively calculable parts up to next-to-next-to-leading order (NNLO), and perform the global fit of the Drell-Yan and Z-boson production, which include the data of E288, Tevatron and LHC experiments. The nonperturbative part of the TMDs are explored checking a large variety of models. We support the obtained results by a study of theoretical uncertainties, perturbative convergence, and dedicated study of the range of applicability of the TMD factorization theorem. The considered nonperturbative models present significant differences in the fitting behavior, which allow us to clearly disfavor most of them. The numerical evaluations are provided by the "arTeMiDe" code, which is introduced in this work and that can be used for current/future TMD phenomenology.

Tables

Table 1. Table 1: The perturbative orders studied in the fit. For each order we indicate the power of a s subscript 𝑎 𝑠 a_{s} of each piece that enters in the TMDs. Note, that the order of a s subscript 𝑎 𝑠 a_{s} and PDF set are related, since the values of a s subscript 𝑎 𝑠 a_{s} are taken from the PDF set.

Name	${\left|C_V\right|}^2$	$C_{f\leftarrow f^{\prime }}$	$\Gamma$	${\gamma }_V$	$𝒟$	PDF set	$a_s$(run)	${\zeta }_{\mu }$
NLL/LO	$a_s^0$	$a_s^0$	$a_s^2$	$a_s^1$	$a_s^2$	nlo	nlo	NLL
NLL/NLO	$a_s^1$	$a_s^1$	$a_s^2$	$a_s^1$	$a_s^2$	nlo	nlo	NLO
NNLL/NLO	$a_s^1$	$a_s^1$	$a_s^3$	$a_s^2$	$a_s^3$	nlo	nlo	NNLL
NNLL/NNLO	$a_s^2$	$a_s^2$	$a_s^3$	$a_s^2$	$a_s^3$	nnlo	nnlo	NNLO

Table 2. Table 2: The characteristics of the data measured at E288 experiment.

	E288 200	E288 300	E288 400
$\sqrt{s}$	19.4 GeV	23.8 GeV	27.4 GeV
process	p+Cu$\to \gamma \to {\mu }^{+}{\mu }^{-}$	p+Cu$\to \gamma \to {\mu }^{+}{\mu }^{-}$	p+Cu$\to \gamma \to {\mu }^{+}{\mu }^{-}$
$Q$ range	4-9 GeV	4-9 GeV	5-14 GeV
$\Delta Q$-bin	1 GeV	1 GeV	1 GeV
y	y=0.4	y=0.21	y=0.03
Observable	$E\frac{d^3\sigma }{d^{3}q}$	$E\frac{d^3\sigma }{d^{3}q}$	$E\frac{d^3\sigma }{d^{3}q}$
Ref.	Ito:1980ev	Ito:1980ev	Ito:1980ev

Table 3. Table 3: The characteristics of the data measured at CDF and D0 collaborations at run 1.

	CDF run I	D0 run I
$\sqrt{s}$	1.8 TeV	1.8 TeV
process	$p+\overline{p}\to Z\to e^{+}{e}^{-}$	$p+\overline{p}\to Z\to e^{+}{e}^{-}$
$M_{ll}$ range	66-116 GeV	75-105 GeV
y	y-integrated	y-integrated
Observable	$\frac{d\sigma }{d{q}_T}$	$\frac{d\sigma }{d{q}_T}$
Exp. ${\sigma }_{\text{tot}}$ [pb]	$248\pm 17$	$\sigma =221\pm 11$
${\sigma }_{\text{tot}}$[pb] Catani:2007vq ; Catani:2009sm ${}_{\text{NLO}}$ Catani:2007vq ; Catani:2009sm ${}_{\text{NNLO}}$:	$223.8\pm 0.05$ $237.63\pm 0.18$	$223.8\pm 0.05$ $237.63\pm 0.18$
Ref.	Affolder:1999jh	Abbott:1999wk ; Abbott:1999yd

Table 4. Table 4: The characteristics of the data measured at CDF and D0 collaborations at run 2.

	CDF run II	D0 run II
$\sqrt{s}$	1.96 TeV	1.96 GeV
process	$p+\overline{p}\to Z\to e^{+}{e}^{-}$	$p+\overline{p}\to Z\to e^{+}{e}^{-}$
$M_{ll}$ range	66-116 GeV	70-110 GeV
y	y-integrated	y-integrated
Observable	$\frac{d\sigma }{d{q}_T}$	$\frac{1}{\sigma }\frac{d\sigma }{d{q}_T}$
Exp. ${\sigma }_{\text{tot}}$ [pb]	$256\pm 2.91$	$\sigma =255$
${\sigma }_{\text{tot}}$[pb] Catani:2007vq ; Catani:2009sm ${}_{\text{NLO}}$ Catani:2007vq ; Catani:2009sm ${}_{\text{NNLO}}$:	$245.0\pm 0.06$ $259.77\pm 0.22$	$245.0\pm 0.06$ $259.77\pm 0.22$
Ref.	Aaltonen:2012fi	Abazov:2007ac

Table 5. Table 5: The characteristics of the Z-boson production data measured by ATLAS collaborations.

	ATLAS	ATLAS
$\sqrt{s}$	7 TeV	8 TeV
process	$pp\to Z\to ee+\mu \mu$	$pp\to Z\to \mu \mu$
$M_{ll}$ range	66 - 116 GeV	66 - 116 GeV
lepton cuts	$p_T>20$ GeV	$p_T>20$ GeV
$y$
Observable	$\frac{1}{\sigma }\frac{d\sigma }{d{q}_T}$	$\frac{1}{\sigma }\frac{d\sigma }{d{q}_T}$
Exp.${\sigma }_{\text{fid}}$[pb]	-	$537.1\pm 0.63\left(\pm 2.8\%\right)$
Theor.${\sigma }_{\text{fid}}$[pb]	Catani:2007vq ; Catani:2009sm ${}_{\text{NLO}}$: $448.56\pm 0.19$ Catani:2007vq ; Catani:2009sm ${}_{\text{NNLO}}$: $471.53\pm 0.94$	Catani:2007vq ; Catani:2009sm ${}_{\text{NLO}}$: $505.53\pm 0.21$ Catani:2007vq ; Catani:2009sm ${}_{\text{NNLO}}$: $531.39\pm 0.93$ Ridder:2016nkl : ${507.9}_{-0.7}^{+2.4}$
Ref.	Aad:2014xaa	Aad:2015auj

Table 6. Table 6: The characteristics of the data for the vector boson production off the Z-peak measured by ATLAS collaborations.

	ATLAS	ATLAS
$\sqrt{s}$	8 TeV	8 TeV
process	$pp\to Z/{\gamma }^{\ast }\to \mu \mu$	$pp\to Z/{\gamma }^{\ast }\to \mu \mu$
$M_{ll}$ range	46 - 66 GeV	116 - 150 GeV
lepton cuts	$p_T>20$ GeV	$p_T>20$ GeV
$y$
Observable	$\frac{1}{\sigma }\frac{d\sigma }{d{q}_T}$	$\frac{1}{\sigma }\frac{d\sigma }{d{q}_T}$
Exp.${\sigma }_{\text{fid}}$[pb]	$14.96\pm 2.62\left(\pm 2.8\%\right)$	$5.59\pm 1.52\left(\pm 2.8\%\right)$
Theor.${\sigma }_{\text{fid}}$[pb]	-	-
Ref.	Aad:2015auj	Aad:2015auj

Table 7. Table 7: The characteristics of the Z-boson production data measured by CMS collaborations.

	CMS	CMS
$\sqrt{s}$	7 TeV	8 TeV
process	$pp\to Z\to ee+\mu \mu$	$pp\to Z\to \mu \mu$
$M_{ll}$ range	60-120 GeV	60-120 GeV
lepton cuts	$p_T>20$ GeV	$p_T>15$ GeV
y
Observable	$\frac{1}{\sigma }\frac{d\sigma }{d{q}_T}$	$\frac{1}{\sigma }\frac{d\sigma }{d{q}_T}$
Norm. exp.	-	-
${\sigma }_{\text{fid}}$[pb] Catani:2007vq ; Catani:2009sm ${}_{\text{NLO}}$ Catani:2007vq ; Catani:2009sm ${}_{\text{NNLO}}$	$379.43\pm 0.16$ $398.27\pm 0.71$	$427.32\pm 0.53$ $448.04\pm 0.83$
Ref.	Chatrchyan:2011wt	Khachatryan:2016nbe

Table 8. Table 8: The characteristics of the Z-boson production data measured by LHCb collaborations.

	LHCb	LHCb	LHCb
$\sqrt{s}$	7 TeV	8 TeV	13 TeV
process	$pp\to Z\to \mu \mu$	$pp\to Z\to \mu \mu$	$pp\to Z\to \mu \mu$
$M_{ll}$ range	60-120 GeV	60-120 GeV	60-120 GeV
lepton cuts	$p_T>20$ GeV	$p_T>20$ GeV	$p_T>20$ GeV
y
Observable	$d\sigma \left(q_T\right)$	$d\sigma \left(q_T\right)$	$\frac{d\sigma }{d{q}_T}$
Norm. exp.	$\sigma =76.0\pm 3.1$ pb	$\sigma =95.0\pm 3.2$ pb	$\sigma =198.0\pm 13.3$ pb
${\sigma }_{\text{fid}}$[pb] Catani:2007vq ; Catani:2009sm ${}_{\text{NLO}}$ Catani:2007vq ; Catani:2009sm ${}_{\text{NNLO}}$	$69.85\pm 0.3$ $74.30\pm 0.21$	$88.98\pm 0.397$ $93.50\pm 0.3$	$185.0\pm 0.09$ $192.78\pm 0.82$
Ref.	Aaij:2015gna	Aaij:2015zlq	Aaij:2016mgv

Table 9. Table 9: The values of χ 2 / d . o . f formulae-sequence superscript 𝜒 2 𝑑 𝑜 𝑓 \chi^{2}/d.o.f for different single-parameter non-perturbative functions f N ​ P subscript 𝑓 𝑁 𝑃 f_{NP} , minimized on different data sets. The χ 2 / d . o . f formulae-sequence superscript 𝜒 2 𝑑 𝑜 𝑓 \chi^{2}/d.o.f values correspond to δ T = 0.2 subscript 𝛿 𝑇 0.2 \delta_{T}=0.2 and NNLL/NNLO.

data	$f_{NP}=e^{-\lambda b}$	$f_{NP}=e^{-\lambda b^2}$	$f_{NP}={\cosh }^{-1}\left(\lambda b\right)$
ATLAS	4.78	1.43	1.42
E288	2.70	5.68	3.64
E288+ATLAS	8.18	5.77	3.72

Table 10. Table 10: The number of points with q T < δ T ​ Q subscript 𝑞 𝑇 subscript 𝛿 𝑇 𝑄 q_{T}<\delta_{T}Q for each data set. In the majority of fits we use δ T = 0.2 subscript 𝛿 𝑇 0.2 \delta_{T}=0.2 , see explanation in the text.

${\delta }_T$	0.1	0.125	0.15	0.175	0.2	0.225	0.25	0.275	0.3
CDF+D0 run1	27	34	38	41	44	47	49	51	52
CDF+D0 run2	22	28	32	38	43	49	55	60	63
ATLAS Z-production (7+8 TeV)	10	12	14	16	18	20	21	23	24
ATLAS DY (8 TeV)	9	11	12	14	14	16	16	18	18
CMS (7+8 TeV)	8	10	10	12	14	16	16	18	18
LHCb (7+8+13 TeV)	18	21	24	27	30	30	33	33	33
High energy total	94	116	130	148	163	178	190	203	208
E288 200 GeV	16	20	24	29	35	36	41	44	47
E288 300 GeV	22	27	33	38	45	46	51	55	59
E288 400 GeV	33	40	49	57	66	69	76	82	85
Low energy total	71	87	106	124	146	151	168	181	191
Total	165	203	236	272	309	329	358	384	399

Table 11. Table 11: The normalization factors for the cross-section for each experiment. The dimensional-less numbers are ratios of partially integrated cross section over q T subscript 𝑞 𝑇 q_{T} ( 53 )(theory/experiment, i.e. 𝒩 − 1 superscript 𝒩 1 \mathcal{N}^{-1} ), for the data with the published value of total cross-section. For the data sets with unpublished values of total cross-section, the value of the total cross-section used for normalization is presented. The numbers are given for the model 1. The variation of the scales and models gives the change of numbers in the unrepresented digits. The numbers shown in bold are those which agree with the measured cross-section within the error bars.

order	ATLAS Z-boson 7TeV	ATLAS Z-boson 8TeV	ATLAS 46-66 8TeV	ATLAS 116-150 8TeV	CMS 7TeV	CMS 8TeV	LHCb 7TeV	LHCb 8TeV	LHCb 13TeV
NLL/NLO	438 pb	0.92	1.01	0.93	369 pb	407 pb	0.92	0.93	0.93
NNLL/NLO	438 pb	0.92	1.01	0.93	368 pb	407 pb	0.92	0.93	0.93
NNLL/NNLO	461 pb	0.97	1.08	0.98	387 pb	429 pb	0.97	0.99	0.98

Table 12. Table 12: The results of the χ 2 superscript 𝜒 2 \chi^{2} -minimization procedure with g K = 0 subscript 𝑔 𝐾 0 g_{K}=0 . The values of χ 2 superscript 𝜒 2 \chi^{2} are given including the theoretical error-band. The values of extracted parameters are given with statistical error-band (the first pair of numbers) and the theoretical error-band (the second pair of numbers). The visual presentation of this table is given in fig. 9 .

Order	$\frac{{\chi }^2}{d.o.f.}$	${\lambda }_1$	${\lambda }_2$
Model 1
NLL/NLO	$2.33{}_{-0.68}^{+2.76}$	${0.321}_{-0.007}^{+0.008}{}_{-0.100}^{+0.095}$	${0.271}_{-0.013}^{+0.014}{}_{-0.063}^{+0.155}$
NNLL/NLO	$1.76{}_{-0.48}^{+1.25}$	${0.289}_{-0.004}^{+0.004}{}_{-0.121}^{+0.007}$	${0.424}_{-0.045}^{+0.051}{}_{-0.139}^{+0.673}$
NNLL/NNLO	$1.34{}_{-0.20}^{+0.44}$	${0.271}_{-0.006}^{+0.007}{}_{-0.073}^{+0.076}$	${0.277}_{-0.012}^{+0.015}{}_{-0.042}^{+0.081}$
Model 2
NLL/NLO	$2.19{}_{-0.64}^{+2.34}$	${0.329}_{-0.008}^{+0.008}{}_{-0.101}^{+0.047}$	${0.289}_{-0.017}^{+0.019}{}_{-0.008}^{+0.276}$
NNLL/NLO	$1.65{}_{-0.39}^{+1.32}$	${0.236}_{-0.004}^{+0.005}{}_{-0.064}^{+0.070}$	${0.440}_{-0.044}^{+0.049}{}_{-0.126}^{+0.573}$
NNLL/NNLO	$1.36{}_{-0.18}^{+0.35}$	${0.284}_{-0.006}^{+0.007}{}_{-0.079}^{+0.074}$	${0.280}_{-0.017}^{+0.019}{}_{-0.034}^{+0.086}$

Table 13. Table 13: The results of the χ 2 superscript 𝜒 2 \chi^{2} -minimization procedure with non-zero g K subscript 𝑔 𝐾 g_{K} . The values of χ 2 superscript 𝜒 2 \chi^{2} are given with theoretical error-band. The values of extracted parameters are given with statistical error-band (the first pair of numbers) and the theoretical error-band (the second pair of numbers). The visual presentation of this table is given in fig. 9 .

Order	$\frac{{\chi }^2}{d.o.f.}$	${\lambda }_1$	${\lambda }_2$	$g_K\times {10}^{-2}$
Model 1
NLL/NLO	$1.17{}_{-0.07}^{+1.32}$	${0.189}_{-0.009}^{+0.009}{}_{-0.052}^{+0.114}$	${0.425}_{-0.045}^{+0.054}{}_{-0.250}^{+0.047}$	${2.31}_{-0.24}^{+0.25}{}_{-1.19}^{+1.44}$
NNLL/NLO	$1.21{}_{-0.02}^{+1.16}$	${0.175}_{-0.008}^{+0.008}{}_{-0.041}^{+0.089}$	${0.532}_{-0.067}^{+0.076}{}_{-0.203}^{+0.426}$	${1.27}_{-0.21}^{+0.22}{}_{-1.27}^{+1.19}$
NNLL/NNLO	$1.23{}_{-0.13}^{+0.30}$	${0.228}_{-0.013}^{+0.016}{}_{-0.060}^{+0.034}$	${0.306}_{-0.026}^{+0.031}{}_{-0.063}^{+0.265}$	${0.73}_{-0.23}^{+0.24}{}_{-0.73}^{+1.09}$
Model 2
NLL/NLO	$1.18{}_{-0.07}^{+1.31}$	${0.199}_{-0.010}^{+0.011}{}_{-0.062}^{+0.104}$	${0.443}_{-0.052}^{+0.061}{}_{-0.093}^{+0.503}$	${2.18}_{-0.25}^{+0.26}{}_{-1.06}^{+1.57}$
NNLL/NLO	$1.22{}_{-0.01}^{+1.16}$	${0.181}_{-0.009}^{+0.009}{}_{-0.045}^{+0.099}$	${0.562}_{-0.075}^{+0.092}{}_{-0.206}^{+0.468}$	${1.18}_{-0.21}^{+0.22}{}_{-1.18}^{+1.12}$
NNLL/NNLO	$1.29{}_{-0.18}^{+0.26}$	${0.244}_{-0.015}^{+0.016}{}_{-0.069}^{+0.035}$	${0.306}_{-0.029}^{+0.034}{}_{-0.050}^{+0.216}$	${0.59}_{-0.27}^{+0.24}{}_{-0.59}^{+1.01}$

Table 14. Table 14: Example of parameter extraction with the variation of c 1 , 2 , 3 , 4 subscript 𝑐 1 2 3 4 c_{1,2,3,4} constants, and evaluation of the theoretical error. Bold numbers in brackets represent the deviation of the parameter from its central value.

Variation	$\frac{{\chi }^2}{d.o.f.}$	${\lambda }_1$	${\lambda }_2$	$g_K\times {10}^{-2}$
Model 1 NNLL/NLO
$c_{1,2,3,4}=1$	$1.17$	$0.189$	$0.425$	$2.31$
$c_1=2$	$1.31\left(+0.14\right)$	$0.201\left(+0.012\right)$	$0.316\left(-0.109\right)$	$3.00\left(+0.69\right)$
$c_1=0.5$	$1.10\left(-0.07\right)$	$0.184\left(-0.005\right)$	$0.308\left(-0.117\right)$	$1.60\left(-0.71\right)$
$c_2=2$	$1.19\left(+0.02\right)$	$0.204\left(+0.015\right)$	$0.223\left(-0.202\right)$	$2.12\left(-0.19\right)$
$c_2=0.5$	$1.20\left(+0.03\right)$	$0.219\left(+0.030\right)$	$0.226\left(-0.199\right)$	$1.93\left(-0.38\right)$
$c_3=2$	$1.23\left(+0.06\right)$	$0.251\left(+0.062\right)$	$0.315\left(-0.110\right)$	$3.75\left(+1.44\right)$
$c_3=0.5$	$1.13\left(-0.04\right)$	$0.160\left(-0.029\right)$	$0.220\left(-0.205\right)$	$1.12\left(-1.19\right)$
$c_4=2$	$1.76\left(+0.59\right)$	$0.137\left(-0.052\right)$	$0.473\left(+0.046\right)$	$2.71\left(+0.40\right)$
$c_4=0.5$	$2.49\left(+1.32\right)$	$0.303\left(+0.114\right)$	$0.175\left(-0.250\right)$	$1.15\left(-1.16\right)$
Result	${1.17}_{-0.07}^{+1.32}$	${0.189}_{-0.052}^{+0.114}$	${0.425}_{-0.250}^{+0.047}$	${2.31}_{-1.19}^{+1.44}$
Model 1 N3LL/NNLO
$c_{1,2,3,4}=1$	$1.23$	$0.228$	$0.306$	$0.73$
$c_1=2$	$1.40\left(+0.17\right)$	$0.242\left(+0.014\right)$	$0.296\left(-0.010\right)$	$1.21\left(+0.48\right)$
$c_1=0.5$	$1.14\left(-0.09\right)$	$0.221\left(-0.007\right)$	$0.346\left(+0.020\right)$	$0.12\left(-0.61\right)$
$c_2=2$	$1.22\left(-0.01\right)$	$0.217\left(-0.011\right)$	$0.295\left(-0.011\right)$	$0.86\left(+0.13\right)$
$c_2=0.5$	$1.26\left(+0.03\right)$	$0.252\left(+0.024\right)$	$0.326\left(+0.020\right)$	$0.48\left(-0.25\right)$
$c_3=2$	$1.27\left(+0.04\right)$	$0.260\left(+0.032\right)$	$0.344\left(+0.038\right)$	$1.82\left(+1.09\right)$
$c_3=0.5$	$1.31\left(+0.08\right)$	$0.198\left(-0.030\right)$	$0.358\left(+0.052\right)$	$0.00\left(-0.73\right)$
$c_4=2$	$1.10\left(-0.13\right)$	$0.168\left(-0.060\right)$	$0.571\left(+0.265\right)$	$1.27\left(+0.54\right)$
$c_4=0.5$	$1.53\left(+0.30\right)$	$0.262\left(+0.034\right)$	$0.243\left(-0.063\right)$	$0.68\left(-0.05\right)$
Result	${1.23}_{-0.13}^{+0.30}$	${0.228}_{-0.060}^{+0.034}$	${0.306}_{-0.063}^{+0.265}$	${0.73}_{-0.73}^{+1.09}$

Table 15. Table 15: The values of χ 2 / p ​ o ​ i ​ n ​ t ​ s superscript 𝜒 2 𝑝 𝑜 𝑖 𝑛 𝑡 𝑠 \chi^{2}/points for individual data sets. The boxes indicate the values of partial χ 2 superscript 𝜒 2 \chi^{2} which are responsible for the increment of χ 2 / d . o . f . formulae-sequence superscript 𝜒 2 𝑑 𝑜 𝑓 \chi^{2}/d.o.f. from NLL/NLO to NNLL/NNLO.

Data set	point	Model 1			Model 2
		NLL/ NLO	NNLL/ NLO	NNLL/ NNLO	NLL/ NLO	NNLL/ NLO	NNLL/ NNLO
CDF run1	30	0.67	0.68	0.64	0.67	0.67	0.64
D0 run1	14	0.50	0.52	0.60	0.49	0.51	0.62
CDF run2	36	1.22	1.36	1.30	1.17	1.29	1.33
D0 run2	7	2.52	2.69	2.75	2.45	2.64	2.79
ATLAS (7TeV) Z-boson	9	1.54	1.55	2.01	1.60	1.59	2.27
ATLAS (8TeV) Z-boson	9	2.32	2.48	2.69	2.46	2.70	2.79
ATLAS (8TeV) 46-66 GeV	5	0.04	0.05	0.16	0.05	0.04	0.20
ATLAS (8TeV) 116-150 GeV	9	0.30	0.35	0.31	0.30	0.36	0.30
CMS (7 TeV)	7	1.38	1.39	1.36	1.38	1.38	1.36
CMS (8 TeV)	7	1.38	1.38	1.54	1.38	1.37	1.58
LHCb (7 TeV)	10	0.26	0.26	0.31	0.25	0.26	0.33
LHCb (8 TeV)	10	0.11	0.12	0.27	0.11	0.12	0.32
LHCb (13 TeV)	10	0.50	0.50	0.28	0.50	0.50	0.27
High energy data	163	0.95	1.00	0.94	0.94	1.00	1.04
E288(200) 4-5 GeV	5	3.86	4.28	3.86	4.25	4.59	4.30
E288(200) 5-6 GeV	6	3.00	3.03	1.92	3.05	3.07	1.92
E288(200) 6-7 GeV	7	1.68	1.68	0.84	1.66	1.67	0.79
E288(200) 7-8 GeV	8	1.10	1.10	0.93	1.13	1.11	1.00
E288(200) 8-9 GeV	9	1.83	1.84	0.78	1.89	1.87	1.87
E288(300) 4-5 GeV	5	1.93	2.20	4.09	2.24	2.44	4.90
E288(300) 5-6 GeV	6	1.15	1.18	1.15	1.19	1.21	1.21
E288(300) 6-7 GeV	7	0.84	0.83	0.66	0.85	0.83	0.69
E288(300) 7-8 GeV	8	1.18	1.17	0.90	1.16	1.17	0.86
E288(300) 8-9 GeV	9	1.13	1.14	1.13	1.11	1.36	1.10
E288(300) 11-12 GeV	12	1.08	1.08	1.00	1.11	1.10	1.04
E288(400) 5-6 GeV	6	2.11	2.04	1.12	1.94	1.92	1.01
E288(400) 6-7 GeV	7	2.59	2.68	2.55	2.59	2.64	2.55
E288(400) 7-8 GeV	8	0.83	0.97	2.02	0.99	1.07	2.44
E288(400) 8-9 GeV	9	1.36	1.31	1.37	1.37	1.32	1.54
E288(400) 11-12 GeV	12	1.08	1.06	1.25	1.05	1.05	1.17
E288(400) 12-13 GeV	12	0.88	0.88	1.10	0.87	0.88	1.14
E288(400) 13-14 GeV	12	0.39	0.38	0.72	0.39	0.39	0.71
Low energy data	146	1.38	1.41	1.35	1.50	1.48	1.49
Total	309	1.17	1.21	1.23	1.18	1.22	1.29

Table 16. Table 16: Exact values of the coefficients of the unpolarized TMD matching coefficient in the parameterization ( B.3 ), for the terms singular at x → 1 → 𝑥 1 x\to 1 .

	$A_1$	$A_2$	$A_3$	$A_4$	$A_5$	$A_6$
${\hat{C}}_{q\leftarrow q}$	${𝐋}_{\mu }^2\left(-14-\frac{64{\pi }^2}{27}+\frac{4N_f}{3}\right)$ $+{𝐋}_{\mu }\left(N_f\left(\frac{4}{3}+\frac{20{\pi }^2}{27}\right)-14-\frac{70{\pi }^2}{9}+\frac{16{\zeta }_3}{3}\right)$ $+N_f\left(\frac{352}{243}+\frac{10{\pi }^2}{27}+\frac{56{\zeta }_3}{27}\right)$ $-\frac{2416}{81}-\frac{67{\pi }^2}{9}+\frac{448{\zeta }_3}{9}+\frac{20{\pi }^4}{81}$	${𝐋}_{\mu }^2\left(-8+\frac{16N_f}{9}\right)$ $+{𝐋}_{\mu }\left(-\frac{1072}{9}+\frac{176{\pi }^2}{27}+\frac{160N_f}{27}\right)$ $+\frac{448N_f}{81}-\frac{3232}{27}+112{\zeta }_3$	$\frac{256}{9}{𝐋}_{\mu }^2$	$-\frac{256}{9}{𝐋}_{\mu }^2$ $-\frac{256}{9}{𝐋}_{\mu }+\frac{200}{9}$	$-\frac{64}{9}$	$0$
${\hat{C}}_{q\leftarrow g}$	$0$	$0$	$0$	$\frac{26}{3}{𝐋}_{\mu }^2-\frac{5}{3}$	$\frac{10}{3}{𝐋}_{\mu }$	$\frac{5}{9}$
${\hat{C}}_{g\leftarrow q}$	$0$	$0$	$0$	$\frac{208}{9}{𝐋}_{\mu }^2$ $+{𝐋}_{\mu }(\frac{32N_f}{9}-\frac{208}{9})$ $+\frac{32N_f}{27}-\frac{184}{9}$	$-\frac{80}{9}{𝐋}_{\mu }+\frac{8N_f}{9}-\frac{44}{9}$	$-\frac{40}{27}$
${\hat{C}}_{g\leftarrow g}$	$-12{\pi }^2{𝐋}_{\mu }^2$ $+{𝐋}_{\mu }\left(-96-108{\zeta }_3+\frac{32N_f}{3}\right)$ $+N_f\left(\frac{548}{27}+5{\zeta }_2-\frac{28{\zeta }_3}{3}\right)-\frac{56N_f^2}{81}$ $-112-\frac{201{\zeta }_2}{2}+154{\zeta }_3+\frac{225{\zeta }_4}{4}$	${𝐋}_{\mu }^2\left(66-4N_f\right)$ $+{𝐋}_{\mu }\left(-268+108{\zeta }_2+\frac{40N_f}{3}\right)$ $+\frac{112N_f}{9}-\frac{808}{3}+252{\zeta }_3$	$144{𝐋}_{\mu }^2$	$-144{𝐋}_{\mu }^2-144{𝐋}_{\mu }$ $+6-2N_f$	$-36$	$0$
${\hat{C}}_{q\leftarrow q^{\prime }}$	$0$	$0$	$0$	$0$	$0$	$0$
${\hat{C}}_{q\leftarrow \overline{q}}$	$0$	$0$	$0$	$0$	$0$	$0$

Table 17. Table 17: Exact values of the coefficients of the unpolarized TMD matching coefficient in the parameterization ( B.3 ), for the terms singular at x → 1 → 𝑥 1 x\to 1 .

	$B_1$	$B_2$	$B_3$	$B_4$	$B_5$
${\hat{C}}_{q\leftarrow q}$	$-\frac{32}{9}{𝐋}_{\mu }^2+{𝐋}_{\mu }\left(\frac{-152}{9}+\frac{16N_f}{9}\right)$ $+\frac{40N_f}{27}-8$	$-\frac{40{𝐋}_{\mu }}{9}-2+\frac{4N_f}{9}$	$-\frac{20}{27}$	0	0
${\hat{C}}_{q\leftarrow g}$	$\frac{14}{3}{𝐋}_{\mu }^2-\frac{22}{3}{𝐋}_{\mu }+\frac{34}{3}$	$\frac{14}{3}{𝐋}_{\mu }-\frac{7}{6}$	$\frac{7}{9}$	$4{𝐋}_{\mu }^2-\frac{51}{3}{𝐋}_{\mu }+\frac{172}{9}-\frac{4{\pi }^2}{3}$	$0$
${\hat{C}}_{g\leftarrow q}$	$-\frac{224}{9}{𝐋}_{\mu }^2+\frac{448}{9}{𝐋}_{\mu }-\frac{200}{9}$	$-\frac{224}{9}{𝐋}_{\mu }+\frac{112}{9}$	$-\frac{112}{27}$	$\left(\frac{104}{3}-\frac{64N_f}{9}\right){𝐋}_{\mu }^2$ $+\left(-16+16{\zeta }_2+\frac{320}{27}{N}_f\right){𝐋}_{\mu }$ $+\frac{896N_f}{81}-\frac{12640}{27}+\frac{352{\zeta }_2}{3}+192{\zeta }_3$	$-32{𝐋}_{\mu }^2$
${\hat{C}}_{g\leftarrow g}$	${𝐋}_{\mu }^2\left(-72+\frac{16N_f}{3}\right)$ $+{𝐋}_{\mu }\left(12+24N_f\right)$ $+\frac{74N_f}{3}-293$	${𝐋}_{\mu }\left(-72+\frac{16N_f}{3}\right)$ $+3+6N_f$	$-12+\frac{8N_f}{9}$	${𝐋}_{\mu }^2\left(-198-\frac{4N_f}{9}\right)$ $+{𝐋}_{\mu }\left(36{\zeta }_2+\frac{244N_f}{9}\right)$ $+\frac{226N_f}{9}-\frac{3160}{3}+264{\zeta }_2+432{\zeta }_3$	$-72{𝐋}_{\mu }^2$
${\hat{C}}_{q\leftarrow q^{\prime }}$	$\frac{8}{3}{𝐋}_{\mu }^2-\frac{8}{3}{𝐋}_{\mu }+\frac{8}{3}$	$\frac{8}{3}{𝐋}_{\mu }-\frac{2}{3}$	$\frac{4}{9}$	$\frac{16}{9}{𝐋}_{\mu }^2-\frac{208}{27}{𝐋}_{\mu }+\frac{688}{81}-\frac{32{\zeta }_2}{9}$	$0$
${\hat{C}}_{q\leftarrow \overline{q}}$	$\frac{16}{9}{𝐋}_{\mu }-\frac{4}{3}$	$\frac{8}{9}{𝐋}_{\mu }$	$\frac{4}{27}$	$0$	$0$

Equations

$$d\sigma =\frac{d^{4}q}{2s}\sum _{G,G^{\prime }=\gamma ,Z}{L}_{G{G}^{\prime }}^{\mu \nu }{W}_{\mu \nu }^{G{G}^{\prime }}{\Delta }_G(q){\Delta }_{G^{\prime }}(q),$$

$$W_{\mu \nu }^{G{G}^{\prime }}$$

$$L_{\mu \nu }^{G{G}^{\prime }}$$

$$W_{\mu \nu }^{G{G}^{\prime }}$$

$$x_1=\frac{\sqrt{Q^2+q_T^2}}{\sqrt{s}}{e}^y\simeq \frac{Q}{\sqrt{s}}{e}^y,x_2=\frac{\sqrt{Q^2+q_T^2}}{\sqrt{s}}{e}^{-y}\simeq \frac{Q}{\sqrt{s}}{e}^{-y}.$$

$$\frac{d\sigma }{d{Q}^{2}dyd(q_T^2)}$$

$$\frac{d\sigma }{d{Q}^{2}dyd(q_T^2)}=\frac{d{\sigma }^{\gamma \gamma }}{d{Q}^{2}dyd(q_T^2)}+\frac{d{\sigma }^{ZZ}}{d{Q}^{2}dyd(q_T^2)}+\frac{d{\sigma }^{\gamma Z}}{d{Q}^{2}dyd(q_T^2)},$$

$$z_{l{l}^{\prime }}^{\gamma \gamma }{z}_{f{f}^{\prime }}^{\gamma \gamma }$$

$$z_{l{l}^{\prime }}^{ZZ}{z}_{f{f}^{\prime }}^{ZZ}$$

$$z_{l{l}^{\prime }}^{Z\gamma }{z}_{f{f}^{\prime }}^{Z\gamma }$$

$$M_Z=91.2\text{GeV},{\Gamma }_Z=2.5\text{GeV},s_W^2=\mathrm{0.2313.}$$

$$F_{q\leftarrow h}(x,\mathbf{b};\zeta ,\mu )=\frac{Z_q(\zeta ,\mu ){\mathcal{R}}_q(\zeta ,\mu )}{2}$$

$$\times \sum _X\int \frac{d{\xi }^{-}}{2\pi }{e}^{-ix{p}^{+}{\xi }^{-}}⟨h|\left\{T{\left[{\bar{q}}_i{\tilde{W}}_n^T\right]}_a\left(\frac{\xi }{2}\right)|X⟩{\gamma }_{ij}^{+}⟨X|\bar{T}{\left[{\tilde{W}}_n^{T†}{q}_j\right]}_a\left(-\frac{\xi }{2}\right)\right\}|h⟩,$$

$${\mu }^2\frac{d}{d{\mu }^2}{F}_{f\leftarrow h}(x,\mathbf{b};\mu ,\zeta )$$

$$\zeta \frac{d}{d\zeta }{F}_{f\leftarrow h}(x,\mathbf{b};\mu ,\zeta )$$

$$\zeta \frac{d}{d\zeta }\frac{{\gamma }_F^f(\mu ,\zeta )}{2}={\mu }^2\frac{d}{d{\mu }^2}(-{\mathcal{D}}^f(\mu ,\mathbf{b}))=-\frac{{\Gamma }^f(\mu )}{2},$$

$${\gamma }_F^f={\Gamma }^f{\mathbf{l}}_{\zeta }-{\gamma }_V^f,{\mathcal{D}}^f=\sum _{n=1}^{\infty }{a}_s^n\sum _{k=0}^n{\mathbf{L}}_{\mu }^k{d}_f^{(n,k)}.$$

$${\mathbf{L}}_X=\ln \left(\frac{{\mathbf{b}}^{2}X}{4e^{-2{\gamma }_E}}\right),{\mathbf{l}}_X=\ln \frac{{\mu }^2}{X}.$$

$${\zeta }_1{\zeta }_2=(2p_1^{+}{p}_2^{-}{)}^2=(Q^2+q_T^2{)}^2\simeq Q^4.$$

$$F_{f\leftarrow h}(x,\mathbf{b};{\mu }_f,{\zeta }_f)=R^f[\mathbf{b};({\mu }_f,{\zeta }_f)\to ({\mu }_i,{\zeta }_i)]F_{f\leftarrow h}(x,\mathbf{b};{\mu }_i,{\zeta }_i),$$

$$R^f[\mathbf{b};({\mu }_f,{\zeta }_f)\to ({\mu }_i,{\zeta }_i)]=\exp \left[{\int }_P\left({\gamma }_F^f(\mu ,\zeta )\frac{d\mu }{\mu }-{\mathcal{D}}^f(\mu ,\mathbf{b})\frac{d\zeta }{\zeta }\right)\right].$$

$$:$$

$$R^f[\mathbf{b};({\mu }_f,{\zeta }_f)\stackrel{1}{\to }({\mu }_i,{\zeta }_i)]$$

$$R^f[\mathbf{b};({\mu }_f,{\zeta }_f)\stackrel{2}{\to }({\mu }_i,{\zeta }_i)]$$

$${\mathcal{D}}^f(\mu ,\mathbf{b})={\int }_{{\mu }_0}^{\mu }\frac{d{\mu }^{\prime }}{{\mu }^{\prime }}{\Gamma }^f+{\mathcal{D}}_{\text{pert}}^f({\mu }_0,\mathbf{b})+g_K{\mathbf{b}}^2,$$

$$R^f[\mathbf{b};({\mu }_f,{\zeta }_f)\to ({\mu }_i,{\zeta }_i);{\mu }_0]=$$

$$\exp \left[{\int }_{{\mu }_i}^{{\mu }_f}\frac{d\mu }{\mu }{\gamma }_F^f(\mu ,{\zeta }_f)-{\int }_{{\mu }_0}^{{\mu }_i}\frac{d\mu }{\mu }{\Gamma }^f(\mu )\ln \left(\frac{{\zeta }_f}{{\zeta }_i}\right)\right]{\left(\frac{{\zeta }_f}{{\zeta }_i}\right)}^{-{\mathcal{D}}_{\text{perp}}^f({\mu }_0,\mathbf{b})-g_K{\mathbf{b}}^2}.$$

$$F_{q\leftarrow h}(x,\mathbf{b};\mu ,\zeta )=\sum _{n=0}^{\infty }{\left(\frac{{\mathbf{b}}^2}{B^2}\right)}^n\sum _f\left(C_{q\leftarrow f}^{(n)}(\mathbf{b};\mu ,\zeta )\otimes f_{f\leftarrow h}^{(n)}(\mu )\right)(x),$$

$$b\ll B$$

$$b\sim B$$

$$C_{q\leftarrow q}^{\text{ren}}(x,\mathbf{b};\mu ,\zeta )$$

$$C_{q\leftarrow q}(x,{\mathbf{L}}_{\mu };\mu ,\zeta )=\delta (\bar{x})+$$

$$a_s(\mu )C_F\left[-2{\mathbf{L}}_{\mu }\left(\frac{2}{(1-x{)}_{+}}-1-x\right)+2\bar{x}+\delta (\bar{x})\left(-{\mathbf{L}}_{\mu }^2+2{\mathbf{L}}_{\mu }{\mathbf{l}}_{\zeta }-\frac{{\pi }^2}{6}\right)\right],$$

$${\mu }^2\frac{dF(x,\mathbf{b};\mu ,{\zeta }_{\mu })}{d{\mu }^2}=0.$$

Full text

1 Introduction

The transverse momentum dependent (TMD) distributions are universal functions that describe the interactions of partons in a hadron. The TMD distributions naturally appear within the TMD factorization theorem for the differential cross section of double-inclusive hard processes. A lot of efforts have been made to arrive to a comprehensive picture of TMD factorization (for the latest works see Collins:2011zzd ; GarciaEchevarria:2011rb ; Echevarria:2012js ; Echevarria:2014rua ; Gaunt:2014ska ; Becher:2010tm ; Chiu:2012ir ; Mantry:2010bi ). In this work we perform a detailed comparison of the experimental measurements with the theory expectations based on our studies of higher-order perturbative expansions and power corrections for unpolarized TMDPDFs made in refs. Echevarria:2015usa ; Echevarria:2015byo ; Echevarria:2016scs ; Scimemi:2016ffw .

Among many different spin (in)dependent TMD distributions, the unpolarized TMD parton distribution functions (TMDPDFs) play a central role. From the practical point of view, their precise knowledge is required to extract further TMD distributions and perform other precision measurements. The ideal process to study the unpolarized TMDPDFs is the unpolarized vector boson production. The data on the $q_{T}$-dependent cross-section for the Drell-Yan process are collected by many experiments, including the precise measurements done by Tevatron and LHC. The theoretical descriptions of Drell-Yan data were made by many groups within one or another form of TMD factorization (see e.g. Landry:1999an ; Qiu:2000hf ; Landry:2002ix ; Watt:2003vf ; Catani:2007vq ; Bozzi:2010xn ; Mantry:2010bi ; Su:2014wpa ; DAlesio:2014mrz ; Catani:2015vma ; Bacchetta:2017gcc ).

This work presents a number of differences with the respect to the previous literature. The collection of the improvements forms a completely new point of view in the TMD phenomenology. The main difference of the present work with respect to the more standard ones (here we consider as the most spread out, and de facto standard, analyses those based on the codes ResBosLandry:2002ix ; Guzey:2012jp and DYqT/DYRes Catani:2007vq ; Bozzi:2010xn ; Catani:2015vma ) are following: (i) We extract the parameters related to individual TMDPDFs, which are suitable for phenomenological description of other TMD-related processes. (ii) We consistently include the perturbative ingredients, such as coefficient functions and anomalous dimensions, at the next-to-next-to-leading order (NNLO), introducing and using the $\zeta$-prescription to solve the problem of perturbative convergence at large-$b$ (where $b$ is the transverse distance). (iii) The TMDPDF parameterization is based on and is consistent with the theory expectation on the TMD behavior with $b$. To our knowledge this is the first attempt to include in a fit both high and low energy data at NNLO precision. The extraction of TMDs takes into account (for the first time to our knowledge) also LHC data. All this represents for us a clear improvement with respect to the more classical analyses.

In a modern view, a TMD distribution is a cumbersome function of many factors, which mix up perturbative and non-perturbative information. In this context, the issue of the separation of perturbative and non-perturbative physics requires a fine analysis and it is open to different solutions. The $\zeta$-prescription proposed in this work, is an attempt to consider the perturbative input to a TMD distribution as it is, without artificial regulators. The $\zeta$-prescription is founded on the fact that the TMD factorization introduces two factorization scales, one for the collinear and one for the soft exchanges. These scales are usually fixed to the same point, while in the $\zeta$-prescription they are chosen to eliminate the problematic double-logarithmic contributions. In other words, the $\zeta$-prescription is based on the freedom to select the normalization and factorization scales, which is guaranteed by the structure of the perturbative theory. The $\zeta$-prescription is essentially different from other used schemes. In particular, it does not strictly solve the problem of the large logarithmic contributions at large-$b$. It only decreases the power of the logarithmic contributions. However, the $x$-dependence of the remaining logarithmic terms has a form which prevents the blow up of the perturbative series, which is not accidental, but the result of the charge conservation. In this way, the $\zeta$-prescription postpones the large logarithm problem to the very far domain of $b$-space, where other factors suppress a TMD distribution. The practical implementation of the $\zeta$-prescription shows that it is efficacious and it allows a very accurate and sound description of the data.

The description of the non-perturbative parts of TMD distributions is the most interesting task. It is highly non-trivial because the definition of the non-perturbative part is strongly affected by schemes and prescriptions used in the perturbative implementation. In this respect a full NNLO can be clarifying. As an example, we recall that the non-perturbative behavior of the TMDPDFs is often assumed to have a Gaussian shape (see e.g. discussions in Landry:2002ix ; Guzzi:2013aja ; Schweitzer:2010tt ; Bacchetta:2017gcc ). Although the Gaussian ansatz is widely used, it comes in the conflict with the usual picture of long-distance strong interaction fueled by light-meson exchanges. The typically expected behavior at long distances is exponential, which is confirmed also by model calculations Schweitzer:2012hh . However, the Gaussian shape is often introduced together with the $b^{*}$-prescription Collins:1981va . Notwithstanding many positive features, the $b^{*}$-prescription has a serious issue: it introduces undesired $b$-even power corrections. In turn, these power correction introduced by $b^{*}$ can easily simulate the Gaussian behavior (see also discussion in Collins:2014jpa ). Once the $b^{*}$-prescription is removed the Gaussian ansatz for the TMD shape is no more essential, according to what we find.

An additional remarkable point of the present study is the wide range of energies covered by the data that we have analyzed. The lowest energy measurements included in the fits have $(Q,\sqrt{s})=(4,19.4)\;\text{GeV}$ (E288 experiment Ito:1980ev ), while the most energetic have $(Q,\sqrt{s})=(116-150,8\cdot 10^{3})\;\text{GeV}$ (ATLAS collaboration Aad:2015auj ). Typically, the low- and high-energy data are considered separately. The main reason for a separate scan is the assumed physical picture of strong interactions, which describes different energies. The description of the high-energy data requires a precise perturbative input and it is expected to be less sensitive to the fine non-perturbative dynamics. The situation is the opposite for the low-energy measurements. Our experience shows that the inclusion of data of different energies is not only possible within TMD formalism, but it is also desired because it cuts away inappropriate models very sharply. We find also that the precision achieved by LHC is already sensitive enough to the non-perturbative structure of TMDs. We show that low and high energy data are sensitive to different regions of $b$-space, and consequently to different non-perturbative regimes of the TMDs: high energy data are better described by a Gaussian non-perturbative correction, while low energy data prefer an exponential type of non-perturbative models. The code (arTemiDe) that we have prepared allows to test all these hypotheses, and can be adaptded also to test different non-perturbative inputs for TMDs.

In order to extract the non-perturbative core of the TMDs, in the present study we choose a neutral tactic. We have scanned many possibilities such as a Gaussian/exponential behavior, with/without inclusion of power corrections, and so on. We have also studied the non-perturbative correction to the evolution kernel. During the examination of models we have prioritize the following criteria:

(i) Stability. The TMD factorization is valid at small-$q_{T}$ (the dilepton transverse momentum) up to a certain limit. Therefore, an acceptable model should produce a stable and good description within the allowed $q_{T}$-range. In other words, the value of $\chi^{2}$ should be valuable and values of parameters are stable independently on the number of included data points, as far as the points belong to the allowed range.

(ii) Convergence. The agreement with data should improve with the increase of the perturbative order. Given the current state of the art of the theory, we can define four successive perturbative orders, which is enough to test the perturbative convergence. Also, the value of the phenomenological non-perturbative constants that one extracts should converge to some central value.

(iii) $\chi^{2}$* minimization.* Naturally, among the models with similar behavior we select the model with the minimal $\chi^{2}$. We have found that it is difficult to find a model (with one or two parameters), which fulfills the demands (i) and (ii), and that at the same time provides a good $\chi^{2}$ value on the whole set of data points (although it is relatively easy to achieve this selecting a particular experiment). The models that we test consider a kind of minimal set of parameters which can be enlarged in future studies, refining the fitting hypotheses.

In the present fit, we have included the measurements of E288 at low-energies, Z-boson production at CDF, D0, ATLAS, CMS and LHCb, and Drell-Yan measurements from ATLAS. To our knowledge, this is the largest set of Drell-Yan data points ever simultaneously considered in a fit within the TMD formalism. We find also that the LHC data below the Z-boson peak and at small $q_{T}$ are very important for current/future TMD studies. In the article we present the most successful models that we have found, and discuss some popular models.

In order to numerically evaluate the theoretical expressions, we have produced the package arTeMiDe. arTeMiDe has a flexible module structure and can be used at any level of TMD theory description, from the evaluation of a single TMDPDF or evolution factor to an evaluation of differential cross-section. The arTeMiDe code is available at web and can be used to check our statements or test a possible future/alternative ansatz (for instance Qiu:2000hf ; Kang:2012am ). In arTeMiDe we have collected all recent achievements of TMD theory, including NNLO matching coefficient function, and N3LO TMD anomalous dimensions. In the current version, arTeMiDe evaluates only unpolarized TMDPDFs and related cross-sections, however, we plan to extend it further.

The body of the article is divided as in the following. In sec. 2 we review the theoretical construction of the Drell-Yan cross section and summarize the theoretical knowledge on unpolarized TMDPDFs. In this section, we also describe all the theoretical improvements which are original for this work. The main original point, namely $\zeta$-prescription is presented in sec. 2.4 and appendix B. The phenomenological studies are presented in sec. 3. This section includes also a dedicated discussion of the shape of the non-perturbative part of the TMD. The allowed range of validity of the TMD factorization is explored in sec. 3.4, the presentation of theoretical uncertainties is given in sec.3.5. The results of the final fit are presented in sec. 3.7. A final discussion and conclusions can be found in sec. 4.

2 Theoretical framework

We consider the Drell-Yan reaction $h_{1}+h_{2}\to G(\to ll^{\prime})+X$, where $G$ is the electroweak neutral gauge boson, $\gamma^{*}$ or $Z$. The incoming hadrons have momenta $p_{1}$ and $p_{2}$ with $(p_{1}+p_{2})^{2}=s$. The gauge boson decays to the lepton pair with momenta $k_{1}$ and $k_{2}$. The momentum of the gauge boson or equivalently the invariant mass of lepton pair is $Q^{2}=q^{2}=(k_{1}+k_{2})^{2}$. The differential cross-section for the Drell-Yan process can be written in the form Drell:1970wh ; Altarelli:1978id

[TABLE]

where $1/2s$ is the flux factor, $\Delta_{G}$ is the (Feynman) propagator for the gauge boson $G$. The hadron and lepton tensors are respectively

[TABLE]

where $J_{\mu}^{G}$ is the electroweak current.

The point of our interest is the $q_{T}$ dependence of the cross-section, where $q_{T}$ is the transverse component of the produced gauge boson in the center-of-mass frame. More precisely, we are interested in the regime $q_{T}\ll Q$, where the TMD factorization formalism can be applied. Within the TMD factorization, one obtains the following expression for the unpolarized hadron tensor (see e.g. Tangerman:1994eh )

[TABLE]

where $g_{T}$ is the transverse part of the metric tensor and the summation runs over the active quark flavors. The variable $\mu$ is the hard factorization scale. The variables $\zeta_{1,2}$ are the scale of soft-gluons factorization, and the subject of relation $\zeta_{1}\zeta_{2}\simeq Q^{4}$. In the following, we consider the symmetric point $\zeta_{1}=\zeta_{2}=\zeta=Q^{2}$. The variables $x_{1,2}$ are the longitudinal parts of parton momenta

[TABLE]

The factors $z_{ff^{\prime}}^{GG^{\prime}}$ are the electro-weak charges. The explicit form of factors $z$ is given in sec. 2.1. The factor $C_{V}$ is the matching coefficient of the QCD neutral current to the same current expressed in terms of collinear quark fields. The explicit expressions for $C_{V}$ can be found in Kramer:1986sg ; Matsuura:1988sm ; Idilbi:2006dg , and are also given in appendix A. The functions $F_{f\leftarrow h}$ are the unpolarized TMDPDFs for quark $f$ in the hadron $h$. They are universal non-perturbative functions and the main objects of our study. The details of their definition and their parametrization are given in sec. 2.3. Finally, the term $Y$ denotes the power corrections to the TMD factorization theorem (to be distinguished from the power corrections to the TMD operator product expansion). The $Y$-term is of the order $q_{T}/Q$ and composed of TMD distributions of the higher dynamical twist. In our study, we restrict ourself to the limit of low $q_{T}$ such that the $Y$-term can be dropped.

Evaluating the lepton tensor, and combining together all factors one obtains the cross-section for the unpolarized Drell-Yan process at the leading order of TMD factorization in the form Collins:1984kg ; Davies:1984sp ; Ellis:1997ii ; Becher:2010tm ; GarciaEchevarria:2011rb ; Collins:2011zzd

[TABLE]

where $y$ is the rapidity of the produced gauge boson. The factor $\mathcal{P}$ is a part of the lepton tensor and contains information on the fiducial cuts. It is discussed in details in sec. 2.6. In the rest of this section a more detailed description of the particular components is presented.

2.1 Expressions for cross-section for different produced bosons

In the case of neutral vector bosons production, the sum over $G$ and $G^{\prime}$ in eq. (2) has three terms

[TABLE]

which correspond to $\gamma^{*}$-production, $Z$-production and interference of $\gamma^{*}$-$Z$ production amplitudes. These three terms of the cross-sections differ from each other only due to the factors $z_{ff^{\prime}}^{GG^{\prime}}$ in eq. (2), which are

[TABLE]

where $M_{Z}$ and $\Gamma_{Z}$ are the mass and the width of the Z-boson, $s_{W}$ and $c_{W}$ are sine and cosine of the Weinberg angle. We use the following of values Olive:2016xmw

[TABLE]

In many studies (see e.g.DAlesio:2014mrz ; Su:2014wpa ; Landry:2002ix ; Bacchetta:2017gcc ; Bozzi:2008bb ) the contribution of $\gamma^{*}$ to the cross-section is neglected in the vicinity of the Z-peak, i.e. the zero-width approximation is used. Here, instead, we include the $\gamma^{*}$ and interference terms in the evaluation of the the cross-section. The inclusion of these terms is important for LHC (in particular ATLAS experiment), where the measurement precision often exceeds the theory precision.

2.2 TMD parton distributions: evolution

The quark unpolarized TMDPDFs are given by the matrix element Collins:2011zzd ; GarciaEchevarria:2011rb ; Echevarria:2016scs

[TABLE]

where $n$ is the light-cone vector along the large component of the hadron momentum, $\xi=\{0^{+},\xi^{-},\bm{b}\}$, $Z$ and $\mathcal{R}$ are the ultraviolet and rapidity divergence renormalization factors. The Wilson lines $W_{n}$ pointing along the direction $n$ to the infinity. For the detailed definition of all constituents in this expression we refer to Echevarria:2016scs .

The peculiar feature of the TMD operator is the presence of two types of divergences and, as a consequence, two renormalization factors $Z$ and $\mathcal{R}$. Firstly, we have ultraviolet divergences, which have their collinear counterpart in the coefficient function $C_{V}$. These divergences are the result of collinear factorization and give rise to the logarithms of the factorization scale $\mu$. Secondly, we have rapidity divergences, which arise in the factorization of the soft-gluon exchanges between partons. The singular soft-gluons exchanges can be collected into the soft factor, which in turn, can be written as a product of rapidity renormalization factors $\mathcal{R}$, see e.g. Echevarria:2015byo ; Echevarria:2016scs ; Vladimirov:2017ksc . This procedure introduces the rapidity factorization scale $\zeta$.

The dependence of TMDPDF on the factorization scales $\mu$ and $\zeta$ is given by the pair of evolution equations

[TABLE]

The TMD anomalous dimensions $\gamma$ and $\mathcal{D}$ are known up to order $a_{s}^{3}$ (see Moch:2005id for $\gamma_{V}$, and Vladimirov:2016dll ; Li:2016ctv ; Vladimirov:2017ksc for $\mathcal{D}$). They satisfy the consistency condition (Cauchy-Riemann condition), which guaranties the existence of the common solution for equations (11-12),

[TABLE]

where $\Gamma^{f}$ is the cusp anomalous dimension. This equation fixes the logarithmic part of the anomalous dimensions. So, the anomalous dimension $\gamma$ is linear in logarithm at all orders, while the rapidity anomalous dimension $\mathcal{D}$ has all powers of logarithms,

[TABLE]

Here and in the following, we use the following notation for logarithms

[TABLE]

The explicit expressions for the anomalous dimensions up to third-loop order can be found e.g. in the appendix of Echevarria:2016scs ; Vladimirov:2017ksc .

The initial values of the factorization scales are dictated by the kinematics of the considered process. In particular, the scales $\zeta_{1,2}$ are related to the momentum of hard partons as

[TABLE]

In the following, we use the symmetric normalization point, $\zeta_{1}=\zeta_{2}=\zeta=Q^{2}$. The $\mu$-dependence cancels between the parts of factorization formula, namely between hard coefficient function $|C_{V}|^{2}$ and the TMDPDFs. The natural choice of $\mu$ is such that logarithms appearing in $|C_{V}|^{2}$ are minimized, i.e. $\mu=Q$. Therefore, TMDPDFs enter in the cross-section in eq. (2) at the hard point $(\mu_{f},\zeta_{f})=(Q,Q^{2})$.

A typical construction of a model for a TMD distribution requires its evolution to a different set of scales. The evolution from $(\mu_{f},\zeta_{f})$ to $(\mu_{i},\zeta_{i})$ takes the form

[TABLE]

where

[TABLE]

Here, the $\int_{P}$ denotes the integration along the path $P$ in the $(\mu,\zeta)$-plane from the point $(\mu_{f},\zeta_{f})$ to the point $(\mu_{i},\zeta_{i})$. The integration can be done on an arbitrary path $P$, and the solution is independent on it, thanks to the Cauchy-Riemann condition eq. (13). At a finite perturbative order, the condition eq. (13) is violated by the next perturbative order. As a consequence the expression for the evolution factor $R$ is dependent on the path of integration. The two simplest choices of integration paths are the combinations of straight segments as

[TABLE]

These paths are depicted in fig. 1. The factor $R$ evaluated along these paths reads

[TABLE]

The numerical difference between these two expressions represents the value of the uncertainty at a given perturbative order.

The expressions for the evolution factor $R$ given in eqs. (2.2-19), contain the rapidity anomalous dimension $\mathcal{D}(\mu,\bm{b})$. The latter contains potentially large values of $\mathbf{L}_{\mu}$, which should be resummed with the help of eq. (13). Additionally, the rapidity anomalous dimension can acquire power corrections from the higher orders in the power expansion of the factorization theorem Becher:2013iya . These power corrections can be also observed in the renormalon structure described in Scimemi:2016ffw . The non-perturbative correction takes the form of a series of even powers of the transverse distance. Therefore, the practical expression for the rapidity anomalous $\mathcal{D}$ is

[TABLE]

where $g_{K}$ is an unknown parameter. Here, $\mathcal{D}_{\text{pert}}^{f}$ is the perturbative expression for $\mathcal{D}$. Correspondingly, the value $\mu_{0}$ should be chosen such that $\mathbf{L}_{\mu_{0}}$ is minimal in the perturbative region. Substituting this expression to the evolution factor, we obtain

[TABLE]

In this form, the evolution factor $R$ is independent on the path of evolution, as can be checked explicitly. The perturbative uncertainty which previously has been given by the variation of evolution path, now is represented by the dependence on the parameter $\mu_{0}$. Thus, using eq. (2.2) the uncertainties of the perturbative calculation can be measured by varying the scale $\mu_{0}$. In the following, we use the evolution factor as in eq. (2.2).

2.3 TMD parton distributions: $b$-space behavior

The TMDPDF is a genuine non-perturbative function, which is to be fitted by a certain ansatz, which covers the whole domain in $b$-space. Different intervals of $b$-space describe different regimes of strong interactions. The schematic picture of $b$-regions is shown in fig. 2. In order to construct an optimal and physically meaningful fitting ansatz, the behavior in every part of the $b$-space should be reproduced. In this section, we collect the main information on the $b$-dependence of TMDPDFs, as it is understood in the modern theoretical picture.

The starting point of our description of a TMD distribution is the small-$b$ operator product expansion (OPE), which results in the series

[TABLE]

where $f^{(n)}$ are PDFs of a $2(n+1)$-twist, $C^{(n)}$ are coefficient functions of OPE and the symbol $\otimes$ represents the convolution in momentum fractions of partons. The parameter $B$ is an unknown non-perturbative parameter which represents an intrinsic hadron scale.

Region 1: In the range $b\ll B$ the TMDPDF is dominated by the $n=0$ term of OPE eq. (22). The leading term is represented by the usual matching onto twist-2 PDFs and reads

[TABLE]

where $C$ is known up to two-loop order Echevarria:2016scs ; Gehrmann:2014yya .

There is a subregion $b\ll 1/Q$, which should be considered specially. While the TMD distribution is completely perturbative within this region, the contributions of this region to the cross-section strongly overlaps with the $Y$-term, eq. (2), which is formally $\mathcal{O}(1/(bQ))$. The behavior of TMD distributions within this tiny range together with the Y-term dictates the asymptotics of the cross-section at large $q_{T}$. As a consequence, it has a significant influence on the value of the total cross-section. In our current evaluation we restrict ourself to the range of small-$q_{T}$ (for a dedicated study of the applicability of this approximation in practice, see sec. 3.4). Therefore, we drop the $Y$-term and do not need any special treatment of $b\ll Q^{-1}$ region.

Region 2: In the range $b\lesssim B$ the OPE is still valid. However, one has to include the higher order terms in addition to the leading one. Very little is known about power suppressed terms of the small-$b$ OPE. Our recent study of the renormalon singularities Scimemi:2016ffw suggests several hints that can be used to model this region:

• (i)

The OPE contains only even powers of $b$. Moreover, the coefficient function of $n$’th order has a prefactor $x^{n}$. In other word, the natural scale of OPE is $x\bm{b}^{2}/B^{2}$ rather then just $\bm{b}^{2}/B^{2}$.

• (ii)

The higher order OPE contributions induced by renormalons, can be summed together to some effective non-perturbative function under the convolution integral.

Therefore, in this region the TMDPDF can be approximated by the form

[TABLE]

where the leading term of the power series in $b/B$ of $G$ is given by $C$. As the power $n$ grows, the sub-leading terms of OPE switch on, which is schematically presented by gray lines in fig. 2. The particular contributions at higher $n$ are not so important in the continuous TMD picture. However,

• (iii)

The $n=1$ contribution to OPE can be estimated by the leading renormalon contribution of order $\sim x\bm{b}^{2}$ Scimemi:2016ffw . It has the form

[TABLE]

where $\Lambda=\Lambda_{QCD}$ is the position of the Landau pole.

Region 3: At $b\gg B$ the small-$b$ OPE cannot be considered as a source of information, and the TMD is completely non-perturbative. Luckily, this region is suppressed by the evolution factor. As a consequence, the cross-section is not very sensitive to the fine structure of TMD distribution in this region, but the general behavior is important. We have tested several asymptotical forms of the TMDPDF, including Gaussian, exponential and power-like and found that the best agreement with the experimental data is achieved with exponential behavior. This observation is in agreement with the general physical intuition, that at high distances the strong forces are dominated by meson exchange, while the Gaussian and power-like asymptotics can not be produced in any simple way.

We should mention that the size of the parameter $B$, as well as, the order of convergence of the small-$b$ OPE, which influences the size of the intermediate region 2, are not known. Our estimations of these characteristic sizes are presented in sec. 4.

2.4 Definition of scaling parameters

The small-$b$ matching is the starting point for the construction of the majority of phenomenological ansatzes for TMD distributions. It can be considered as an additional collinear factorization, which is performed at some convenient set of scales $(\mu_{i},\zeta_{i})$. The difference of $(\mu_{i},\zeta_{i})$ from the initial (defined by process kinematic) scales of TMD distribution is compensated by the evolution factor in eq. (17). As usual, the all-order expression is independent on $(\mu_{i},\zeta_{i})$, but in practice, these scales are to be chosen such that the coefficient function $C_{f\leftarrow f^{\prime}}$ has good perturbative convergence. This procedure is alike the choice of hard-factorization scale, with one essential difference: the parameter $b$, which accompanies $\mu_{i}$ and $\zeta_{i}$ in the logarithms, has no fixed value. It varies from zero to infinity within the Fourier integral.

The choice of scales $(\mu_{i},\zeta_{i})$ is one of the central point of the TMD phenomenology. To make the discussion clearer, let us recall the expression for the coefficient function at NLO. It reads

[TABLE]

where the notation for the logarithms is defined in eq. (15). Ideally, the scales $\mu$ and $\zeta$ should be chosen such that no large perturbative contributions appear in the coefficient function. Clearly, it cannot be done at arbitrary $b$ due to the presence of $\mu$ in the coupling constant and in $\mathbf{L}_{\mu}$. However, such a strict choice is not required. The only requirement for scales is to keep the perturbative ansatz stable, i.e. to prevent its blowing up. There are several solutions of this problem. The most famous is the $b^{*}$-prescription Collins:1981va . Within the $b^{*}$-prescription the logarithms $\mathbf{L}_{\mu}$ are absent, an this fact allows the control of the perturbative series in the full region of $b$. However, the $b^{*}$-prescription introduces artificial power corrections to the small-$b$ OPE, which washes out any theoretical intuition. Another popular scheme Laenen:2000de ; Kulesza:2002rh is based on the re-expression of Hankel-integral as an integral in the complex $b$-plane. In this way, the logarithms $\mathbf{L}_{\mu}$ can be minimized by $\mu\sim b^{-1}$ and the Landau pole at large-$b$ is by-passed in the complex plane. The drawback of this scheme is the necessity of the analytical continuation into the complex $b$-plane, and the restriction to NNLO (since the analytical solution for running coupling at N3LO is unknown).

In this work we use another scheme which we call $\zeta$-prescription. It is a novel one (to our best knowledge), and it is described in the following.

The $\zeta$-prescription uses the fact that the TMD operator and hence its small-$b$ OPE depends on two scales $\mu$ and $\zeta$, which are entirely independent. This simple fact has been overlooked so far. Indeed, the first typical step is to fix $\zeta=C_{0}^{2}/\bm{b}^{2}$, or $\zeta=\mu^{2}$ Aybat:2011zv ; Collins:2011zzd ; Scimemi:2016ffw . It reduces the problem to a single variable problem, which looks simpler, but finally, it does not provide a simple solution for the appearance of large logarithms in the OPE.

The initial point of the $\zeta$-prescription is the observation that not all logarithms in the coefficient function are dangerous. So, the terms $\mathbf{L}^{2}_{\mu}$ and $\mathbf{L}_{\mu}\mathbf{l}_{\zeta}$ in eq. (26) are problematic, while the logarithm in first term is not. There are several reasons for it. First, the double logarithm contributions violate the normal perturbative counting and at large-$b$ grows faster than the single logarithms. Second, the first term of eq. (26) comes together with the DGLAP kernel, and thus, it preserves the area (say, the integral over $x$) of the TMDPDF, due to the conservation of the electromagnetic charge. We remind that logarithms accompanying the DGLAP kernel are related to PDF evolution, while the rest of logarithms are related to the TMD evolution. For this reason, the main problem of convergence is represented by the logarithms that are related to the TMD evolution. The logarithms related to the PDF evolution come with a particular $x$-dependent function. The probabilistic interpretation of PDF ensures their minimal contribution in the very large domain of $b$. Practically, this fact has been already demonstrated although not entirely realized in the fit DAlesio:2014mrz . In the realization of ref. DAlesio:2014mrz , the DGLAP logarithms were left unregulated and they do not influence the convergence of the fit.

The logarithms related to the TMD evolution can be eliminated completely by a particular choice of $\zeta=\zeta_{\mu}$. Along the curve $\zeta_{\mu}$, the TMD distributions are independent on $\mu$. In other words, the curve $\zeta_{\mu}$ is an equi-evolution curve in the plane $(\mu,\zeta)$. Such a curve satisfies the equation

[TABLE]

Using the definition of anomalous dimensions in eq. (11) we rewrite this equation as

[TABLE]

where $f(\mathbf{L}_{\mu})=\mathbf{l}_{\zeta_{\mu}}$. The perturbative solution is discussed and presented in the appendix B.1. The curve $\zeta_{\mu}$ is different for quark and for gluon TMDs, and it is expressed in terms of the TMD anomalous dimensions eq. (B.1). In our analysis, we need only the quark case. Up to NNLO it reads

[TABLE]

Note, that in eq. (2.4) we have set the boundary condition such that no terms singular at $\mathbf{L}_{\mu}\to 0$ appear in $\mathbf{l}_{\zeta}$ (see appendix B.1, for details). Also, in the current work we drop the power contributions to the rapidity anomalous dimension $\mathcal{D}$. The influence of these decisions should be investigated later. One can check that the leading term of $\zeta_{\mu}$ (i.e. $\mathbf{l}_{\zeta}=\mathbf{L}_{\mu}/2$) cancels leading powers of logarithms at all orders in perturbation theory (i.e. all terms $a_{s}^{n}\mathbf{L}_{\mu}^{2n}$). Then, including the next correction ($a_{s}\beta_{0}\mathbf{L}_{\mu}^{2}/12$) cancels subleading powers of logarithms at all orders of the perturbation theory (i.e. all terms $a_{s}^{n}\mathbf{L}_{\mu}^{2n-1}$) , and so on.

Substituting the leading term of the solution in eq. (2.4) to the quark small-$b$ coefficient function, we obtain

[TABLE]

This coefficient function is stable for any value of $\mathbf{L}_{\mu}$, which can be seen by considering its integral

[TABLE]

which is independent on $\mathbf{L}_{\mu}$.

The general expression for the coefficient of arbitrary flavour at NNLO has the form

[TABLE]

where $C^{(n,0)}$ is the finite part of the coefficient function at $n$’th perturbative order, and $P(x)=\sum a_{s}^{n}P^{(n)}$ is the DGLAP kernel. The detailed derivation of eq. (2.4) is presented in the appendix (B.2). Eq. (2.4) has the form of the usual coefficient function for an object without external evolution (e.g. coefficient function for DIS). In other words, it is straightforward to check that

[TABLE]

by direct differentiation of eq. (2.4). The integral of this function over $x$ is independent on $\mathbf{L}_{\mu}$ due to the charge conservation, and thus at least the area of TMDPDF is stable at large $b$.

A further positive point of the $\zeta$-prescription is that the scale $\mu$ remains unconstrained. Often, the scale $\mu$ is selected such that it behaves as $\sim 1/b$ at $b\to 0$. Such a choice minimizes the small-$b$ logarithms in small-$b$ OPE and in the evolution exponent. At large-$b$ the scale $\mu$ should be frozen to some fixed value (of the order of a few GeV’s), in order to avoid the Landau pole. We use the simplest function which satisfies these criteria

[TABLE]

where the low-energy fixed scale $2$ GeV is chosen due to the fact that it is the standard scale of PDF extractions.

Finally, we should also select the value for the parameter $\mu_{0}$ that enters expression for the evolution factor eq. (2.2). To keep our discussion simple, we set $\mu_{0}=\mu_{b}$.

2.5 Theoretical uncertainties and perturbative ordering

In the construction of the cross section, one finds several sources of perturbative uncertainties. The size of these uncertanties can be estimated by the variation of associated scales. We list here the ones that we have considered in the present work.

• •

Uncertainty associated with the perturbative matching of rapidity anomalous dimension : This uncertainty arises from the dependence (at the fixed perturbative order) on $\mu_{0}$, which should be compensated between the Sudakov factor and the boundary term $\mathcal{D}(\mu_{0})$ in the TMD evolution factor eq. (2.2). This uncertainty can be tested by changing $\mu_{0}\to c_{1}\mu_{0}$ and varying $c_{1}\in[0.5,2]$.

• •

Uncertainty associated with the hard factorization scale: This uncertainty arises from the dependence (at the fixed perturbative order) on the scale $\mu_{f}(\sim Q)$ which is to be compensated between the hard coefficient function $|C_{V}|^{2}$ and the TMD evolution factor. This uncertainty can be tested by changing $\mu_{f}\to c_{2}\mu_{f}$ and varying $c_{2}\in[0.5,2]$.

• •

Uncertainty associated with the TMD evolution factor: This uncertainty arises from the dependence (at the fixed perturbative order) on initial scale of TMD evolution $\mu_{i}$, which is to be compensated between the evolution integral and the $\mu$-dependence of $\zeta_{i}$ in eq. (2.2). This uncertainty can be tested by changing $\mu_{i}\to c_{3}\mu_{i}$ and varying $c_{3}\in[0.5,2]$.

• •

Uncertainty associated with the small-$b$ matching: This uncertainty arises from the dependence (at the fixed perturbative order) on the scale of the small-$b$ matching $\mu_{\text{OPE}}$ which is to be compensated between the small-$b$ coefficient function $C_{f\leftarrow f^{\prime}}$ and evolution of PDF. This uncertainty can be tested by changing $\mu_{\text{OPE}}\to c_{4}\mu_{\text{OPE}}$ and varying $c_{4}\in[0.5,2]$.

We remark that our definition of perturbative uncertainties $c_{1,2}$ is commonly used in the literature (as far as it can be compared among different schemes of calculation), see e.g. Catani:2015vma ; Bozzi:2008bb . The uncertainties $c_{3,4}$ are usually non distinguished and they are commonly varied simultaneously i.e. in the literature one finds discussions of errors for the case $c_{4}=c_{3}$. To our best knowledge, the distinction of the matching and evolution uncertainties is made here for the first time.

In this way, the general expression for the cross-section in eq. (2) with our choice of scales reads

[TABLE]

where the evolution factor $R$ is given in eq. (2.2) and the explicit expression for the $\zeta_{\mu}$ is given in eq. (2.4). The low-normalization point $\mu_{i}$ and the scale of small-b operator product expansion $\mu_{\text{OPE}}$ are fixed at the same point (34)

[TABLE]

The central value of the constants $c_{1,2,3,4}$ is $1$ and they are varied in order to estimate the theoretical uncertainties in the usual range $(0.5,2)$.

The perturbative orders of each constituent are to be combined consistently. Having at our disposal the NNLO expressions for coefficient function and N3LO expressions for anomalous dimensions, we can define four successive self-contained sets of ordering. This is reported in table 1. In our definition of orders we use the following logic: (i) The order of the $a_{s}$-running should be the same as the order of PDF set, since their extraction are correlated. (ii) The order of $\mathcal{D}$ should be the same as the order of $\Gamma$ since they enter the evolution kernel $R$ with the same counting of logarithms (i.e. $a_{s}^{n}\ln^{n+1}\mu$), and one-order higher then the order of $\gamma_{V}$, since it has counting $a_{s}^{n}\ln^{n}\mu$. (iii) The order of small-$b$ matching coefficient should be the same as the order of evolution of a PDF, because large logarithms of $b$ are to be compensated by the PDF evolution. (iv) The order of $\zeta_{\mu}$ should be such that no logarithms appear in the coefficient function, and the general logarithm counting coincides with the counting of the evolution factor. In table 1 the order of the $\zeta_{\mu}$ is defined as following: NLL is $\mathbf{l}_{\zeta}=\mathbf{L}_{\mu}/2$, NLO has in addition finite part at order $a_{s}^{0}$ (i.e. two first terms of eq. (2.4)), NNLL has in addition logarithmic part at order $a_{s}^{1}$ (i.e. the first line of eq. (2.4)), and NNLO is given by whole expression eq. (2.4).

To label the orders we use the nomenclature where the part with ’LO suffix designates the order of coefficient functions, and the part with ’LL suffix designates the order of the evolution factor in the $a_{s}\ln\mu\sim 1$ scheme. So, our highest order is NNLL/NNLO, which at the moment the highest available combination of the perturbative series. The order NLL/LO appears to be barely inconsistent, because it requires the LO PDF evolution to match the trivial coefficient function. Therefore, we exclude the NLL/LO from our phenomenological studies.

2.6 Implementation of lepton cuts

In modern experiments, the cross-section is often evaluated with the fiducial cuts on the dilepton momenta. That is, the lepton tensor in eq. (2) should be evaluated taking into account the experimental cut phase-space. At the leading order the lepton tensor takes the form

[TABLE]

where $\theta$-function restricts the lepton momenta to the allowed region.

In the limit $Q\to\infty$ and no restriction on the lepton pair phase space we obtain

[TABLE]

Substituting this expression to the cross-section we obtain the standard formula to the Drell-Yan cross-section within TMD factorization Collins:1984kg ; Davies:1984sp ; Ellis:1997ii ; Becher:2010tm ; GarciaEchevarria:2011rb ; Collins:2011zzd . In order to include the corrections due to a finite $Q$ and experimental cuts let us introduce a factor $\mathcal{P}$, i.e.

[TABLE]

which is consistent with the cross section expression presented in eq. (2). The function $\mathcal{P}$ in the absence of cuts reads

[TABLE]

In the presence of cuts the expression for $\mathcal{P}$ is involved. For example, at $q_{T}=0$ and $y=0$ it reads

[TABLE]

Generally, $\mathcal{P}$ cannot be evaluated analytically, but it is rather easy to evaluate numerically.

3 Comparison with experiment

3.1 Review of experimental data

In this section we present the experimental data that have been included in our fit. We have splitted the data into two large data sets with respect to a general energy scaling. They include the measurements from the following experiments:

• •

Low-energy data set

• –

E288: Drell-Yan process, at $4<Q<14$ GeV.

• •

High energy data set:

• –

CDF/D0: Z-boson production at $s=1.8,~{}1.96$ TeV.

• –

ATLAS/CMS/LHCb: Z-boson production at $s=7,8,13$ TeV.

• –

ATLAS: Vector boson production outside the Z-peak ($46<Q<66$ and $116<Q<150$ GeV) at $s=8$ TeV.

In the present study, we have not included the data of other experiments, such as E605, or R209. The main reason is that these data require some special discussion, due to possible internal inconsistencies (see e.g. DAlesio:2014mrz ), or due to low precision. We expect that the results of the present work are not sensibly affected by this omission.

In the following, we present each included measurement in more detail.

E288.

The E288 experiment Ito:1980ev presents a large number of low energy points. So, the total number of low-energy point is nearly equal the total number of points of high energy experiments. For convenience we have splitted this data set into three subsets with different center of mass energy $s$. The characteristics of the measurements are shown in table 2. Concerning these data we can comment the following:

• •

We exclude the data points in the range $9<Q<11$ GeV, because these data are dominated by the $\Upsilon$-resonance ($M_{\Upsilon}\simeq 9.5$ GeV). The description of $\Upsilon$-resonance production is beyond the scope of current TMD factorization approach.

• •

The E288 experiment is made on a copper target. To simulate the effects of copper nuclei we replace the proton PDFs by the following combinations

[TABLE]

where $Z=29$, $A=63$ and $N=A-Z=34$, are charge, atomic number and the number of neutrons in copper correspondingly.

• •

The absolute normalization of the E288 $p_{T}$-cross-section is unknown. Typically, one includes an additional normalization factor $N_{E288}$, as a parameter of the fit, see e.g. DAlesio:2014mrz ; Su:2014wpa ; Landry:1999an ; Landry:2002ix . There is no agreement on this factor values, it varies from $\sim 0.8$ DAlesio:2014mrz ; Landry:1999an ; Su:2014wpa to $\sim 1.2$ Landry:2002ix . In our analysis we fix $\mathcal{N}_{E288}=0.8$.

The theoretical uncertainties for low energy experiments are large, of the order $\pm 10\%$ at the best (see sec. 3.5). As a consequence, the value of the cross-section is very sensitive to the choice of the PDF set and the overall normalization factor. For example, we have checked that the E288 data can be fitted also with $N_{E288}=0.9$ with the same (or better) value of $\chi^{2}$ by an additional variation of $\mu_{b}$. However, we consider this as a bad practice and restrict ourself to $\mathcal{N}_{E288}=0.8$, as the most conventional solution.

• •

The data are splitted into different bins with different dilepton invariant mass. For each bin we evaluate the cross-section eq. (2) as

[TABLE]

where $Q_{\max,\min}$ are the boundary of the $Q$-bin.

CDF and D0.

The data on the Z-boson production measured by CDF and D0 collaborations at Tevatron Run 1 and Run 2 Affolder:1999jh ; Abbott:1999wk ; Abbott:1999yd ; Aaltonen:2012fi ; Abazov:2007ac have been used nearly in every fit of unpolarized TMDPDFs. They are summarized in tables 3-4. Concerning these data we can comment the following:

• •

There is a known tension between the values of total cross-section at run I of CDF and D0. Here we restrict ourself to the fit of the shape of the cross-section and normalize the theoretical points on the bin-by-bin integrals in the allowed range of $q_{T}$. I.e. we multiply the theoretical cross-section by the factor

[TABLE]

where $\Delta q_{T}$ is the size of $q_{T}$-bins. As we show in sec. 3.6 the obtained normalization factors are very close to one (at NNLO), and the values of partial cross-sections are in agreement with the experimental ones within error-bars. In the tables 3-4, we also present the values of the total cross-sections evaluated by DYNNLO code Catani:2007vq ; Catani:2009sm . In this calculation of the total-cross-section, we have used the same inputs as in the TMD fits, i.e. the PDF are taken from MMHT2014 set Harland-Lang:2014zoa .

• •

The experimental values for cross-section points are obtained by integrating over all values of $y$, integrating over measure range of $Q$ and averaging in $q_{T}$. Consequently, we have used the following expression for the cross-section

[TABLE]

where $y_{0}=\frac{1}{2}\ln(s/Q^{2})$, $q_{T,\text{low}}$ and $q_{T,\text{high}}$ are boundaries of $q_{T}$-bin, and $\Delta q_{T}$ is the size of the $q_{T}$-bin.

ATLAS.

The data by ATLAS collaboration in Aad:2014xaa ; Aad:2015auj cover a broad range of dilepton invariant masses for the Drell-Yan process with small experimental error-bands. So, this set provide the biggest constraints on TMD extraction coming from high energy data points. The characteristics of the measurements are resumed in tables 5-6 and here we comment the following:

• •

The data from the ATLAS detector at 8 TeV run are presented in several sets Aad:2015auj , which corresponds to different treatment of final-state photon radiation. We have considered the "dressed" set of the data.

• •

The values of cross-section have been calculated by the expression in eq. (48), where $y_{0}=2.4$, as it is presented in the tables 5-6.

• •

There is a known tension between the theoretical calculation of the integrated cross-section and the measured one, see e.g. Aad:2015auj ; Ridder:2016nkl . Moreover the theoretical cross-section for vector boson production is not known (at least, the DYNNLO package Catani:2007vq ; Catani:2009sm , which we have used for evaluation of total cross-sections, does not guarantee the accurate calculation outside the Z-peak). Therefore, we normalize the calculated cross-sections by a factor, as explained in more detail in the text around eq. (53). In sec. 3.6, we compare the obtained values of normalization to the total cross-section. We have found that the values of obtained normalization are practically independent on the non-perturbative input of the TMD model, and at NNLL/NNLO correctly reproduce (within the error-bars) the measured total cross-sectio.

• •

All data sets from LHC are presented within fiducial cross-sections. Therefore, we have implemented the cut leptonic tensor as it is discussed in sec. 2.6.

CMS and LHCb.

The CMS and LHCb collaborations provide data around the Z-boson peak in Chatrchyan:2011wt ; Khachatryan:2016nbe ; Aaij:2015gna ; Aaij:2015zlq ; Aaij:2016mgv , see tables 7,8. The treatment of these data is similar to the case of ATLAS data:

• •

The values of cross-section have been calculated by the expression in eq. (48), where the limits for $y$-integration $y_{0}$ are taken in accordance to the tables 7-8.

• •

Just as in the case of ATLAS data we have normalized the calculated cross-sections by the factor provided in eq. (47) discussed in sec.3.6. We have found a good agreement between the theoretical and experimental values for total cross-section for LHCb data.

• •

All data sets from LHC are fiducial cross-sections. Therefore, we have implemented the cut leptonic tensor as it is discussed in sec. 2.6.

Finally, we have considered only points which allow a consistent TMD treatment. I.e. the points with the value of $q_{T}<\delta_{T}Q$, where $\delta_{T}$ is sufficiently small. In the literature we have not found any special study on this limiting ratio. So, we present our study in sec. 3.4, and conclude that TMD factorization range is $q_{T}/Q<0.2$.

3.2 arTeMiDe

In order to evaluate the cross-sections we have prepared the program package arTeMiDe. The arTeMiDe package is a collection of FORTRAN modules that evaluates individual terms of the TMD factorization formalism, such as TMD evolution factors, TMDPDFs, and combines them into the differential cross-sections. arTeMiDe forms a flexible package for TMDPDF phenomenology based on the $\zeta$-prescription, as described in this article. It is publicly available at the web-page web .

arTeMiDe version 1.1 evaluates the quark and gluon unpolarized TMDPDFs (although in the discussed fit the gluon TMDPDFs are not necessary) for any given function $f_{NP}$, at any composition of perturbative orders from LO to NNLO, with or without renormalon-induced power corrections. For the current study, the input PDFs are taken from the MMHT2014 PDF set Harland-Lang:2014zoa .

The most time-consuming part of the numerical evaluation of the TMDPDFs, is the convolution integral in eq. (50), which is especially expensive at NNLL/NNLO. Within the arTeMiDe package the convolution integral is optimized using an approximate expression for NNLO coefficient functions. The approximate form of the NNLO coefficient function is (note, that NLO and renormalon coefficient functions can be presented in this form without approximation)

[TABLE]

where coefficients $A$, $B$ and $c$ are functions of $\mathbf{L}_{\mu}$. Such an approximate form is widely used in NNLO+ phenomenology of PDFs, see e.g. vanNeerven:1999ca . Here, the coefficients $A$ and $B$ represent the singular at $x\to 1$ and $x\to 0$ terms, and are evaluated exactly. The coefficients $c$ represent the smooth part of the coefficient function, which is reconstructed by appropriate values of $c_{i}$ with better then $1\%$ accuracy. The values of constants $A$, $B$ and $c$ are presented in the appendix B.3. In the convolution integral the main numerical contribution comes from the singular terms proportional to $A$ and $B$, which are exact. The relative difference between the convolution with exact coefficient function and approximate expression in eq. (3.2) is of order $10^{-6}$. This numerical error is compatible with the numerical error of integration procedure and far inside the theoretical error-bands.

The evaluation of the Hankel-type integral over $b$ is one of the main source of numerical errors. Typically, in order to obtain sufficient precision one should include a large number of points into the integral, which is very costly especially at NNLL/NNLO. arTeMiDe evaluates this integral with the Ogata quadratures Ogata_quadrature . The Ogata quadrature is a double exponential quadrature, whose nodes are the zeros of the Bessel function. It provides a fast and precise evaluation of Hankel-type integrals with the minimal number of integrand calls.

The fitting procedure has been performed by minimizing the $\chi^{2}$-function. The minimization of the $\chi^{2}$ distribution have been done using the MINUIT package from the CERN library James:1975dr ; James:1994vla . The estimation of the statistical uncertainties for non-perturbative parameters is made with the MINOS procedure, performing the variation of parameters in the range $\chi^{2}\pm\Delta\chi^{2}$, with $\Delta\chi^{2}$ corresponding to the 68% confidence level (i.e. $\Delta\chi^{2}\simeq\{1.03,2.32,3.55\}$ for 1,2,3 fitting parameters, correspondingly.) The sources of theoretical uncertainties have been pointed in sec 2.5, and parameterized by the constants $c_{1,2,3}$. The variation of these constants in the region $(0.5,2)$ produces the error-bands. The discussion on the individual contributions of theoretical uncertainties associated with different scales is given in sec. 3.5.

3.3 Models for non-perturbative part of TMDPDFs

The non-perturbative part of the TMDPDF in general needs some ansatz, the parameters of which are to be extracted from data. In our study we have tested different ansatzes of the following general form

[TABLE]

where $f_{f\leftarrow h}$ is PDF of the parton with the flavour $f$. The non-perturbative information of the TMDPDF, which is unreachable from the PDFs, is contained in $f_{NP}$. In order to match the perturbative regime, the function $f_{NP}$ should approach $1$ for $b\to 0$. Instead, the behavior of $f_{NP}$ for $b\to\infty$ is not so well established, which requests a test of different models. In the current study, we restrict ourself to flavor independent $f_{NP}$, i.e. $f_{NP}$ is common for TMDPDFs of different flavours. The flavour-dependence of TMDPDFs enters only through PDFs and coefficient functions, i.e. it is completely determined.

The large-$b$ behavior of TMD distributions is the key point of TMD parametrization and extraction. There is no common agreement on this behavior. Clearly, such an agreement cannot be achieved in general, since the $b$-shape of a TMD distribution is strongly dependent on the large-$b$ prescription. For example, the Gaussian behavior is typically observed in the models based on $b^{*}$-prescription. Moreover, the classical fits by ResBos package Landry:2002ix disfavor other non-perturbative behaviors, such an exponential one (for more recent discussion, see Guzzi:2013aja ). Also the Gaussian shape is used in DYRes code Catani:2015vma (together with $b^{*}$-prescription) and in DYqT code Bozzi:2010xn (together with the minimal prescription). Contrary, the fit made in ref. DAlesio:2014mrz , which does not employ the $b^{*}$-prescription, uses an exponential shape of $f_{NP}$ and also obtains an agreement with data. We point out that the use of LHC data for TMD extraction is made here for the first time (to our knowledge). Given the precision of LHC data, the consistency and/or goodness of all previous hypotheses have to be rediscussed.

In order to decide the best shape of $f_{NP}$ within $\zeta$-prescription, we have considered several subsets of the data. It appears important to include simultaneously both high-energy and low-energy data because they are sensitive to different parts of the $b$-space spectrum. We have found that the most optimal data subset is given by the E288 data and the ATLAS Z-boson production data, see tables 2 and 5. In this subset, the very small error-bands of ATLAS data are compensated by a large number of points in E288 data, and as a result, we have a certain equilibrium between low and high-energy inputs.

Using the E288/ATLAS subset we have performed multiple fits using several different functional forms of $f_{NP}$. Probably, the most informative preliminary test is the comparison of the pure Gaussian and exponential behavior for separate/joint low and high energy data points. In table 9 we demonstrate results of fit with some simple single-parameter models. According to this table, although the quality of the fit is still not optimal, the high-energy data clearly favor the Gaussian shape of $f_{NP}$, while the low-energy data favor the exponential behavior of $f_{NP}$. This difference is simply explained if we recall that at higher energies (and thus at generally higher $q_{T}$) the Fourier integral in eq. (2) is saturated by small values of $b$. At lower energies (and thus at generally smaller $q_{T}$) the Fourier integral in eq. (2) is affected by a wider interval of values of $b$. Therefore, the results presented in the table 9, suggest that $f_{NP}$ should be Gaussian at small-$b$ and exponential at large-$b$. This is in complete agreement with the theory expectations discussed in sec. (2.3). The expected $f_{NP}$ should be a function with a Taylor series expansion (around $b=0$) of even powers of $b$, with an exponential decay at $b\to\infty$. A simple representative of such functions is $\cosh^{-1}(\lambda b)$. The test of this $f_{NP}$ is given in the last columns of the table 9 which clearly shows that this function alone, although it works much better than a Gaussian or an exponential, is not able to describe both low and high energy data, and thus we need extra parameters.

The preliminary tests with simple one-parameter dependence for the $f_{NP}$ shape can be summerized by the following:

• (i)

The high and low energy data should be considered altogether, because they test different intervals of the $b$-space spectrum of $f_{NP}$.

• (ii)

The subset of data points E288 and ATLAS Z-boson, is very selective for the $f_{NP}$. A good fit of this subset guaranties the good fit for the whole set of data points. Nevertheless, in the following sections, we include all experiments, for consistency.

• (iii)

Both theoretically and phenomenologically, we argue that $f_{NP}$ should be a function of even powers of $b$ with an exponential asymptotic behavior at $b\to\infty$. Using a minimal set of two parameters (and the evolution parameter $g_{K}$) we find that one can easily fit the data with a $\chi^{2}/d.o.f\sim 1.2$-1.3. The addition of more parameters (say for the control of $b^{4}$ correction and/or flavor dependence) has the possibility to increase the quality of the fit. However, in this work, we do not consider extra parameters, since the current quality of the fit is already typical and reasonable for the modern TMD extraction (compare e.g. with Bacchetta:2017gcc ).

• (iv)

One needs at least two parameters (one to control $\sim b^{2}$ behavior at $b\to 0$ and another to control the asymptotics) to fit simultaneously low and high-energy data. However, the multiplication by polynomials (e.g. $f_{NP}\sim(1+\lambda b^{2})/\cosh(b)$) does not work well, which suggests that the asymptotic terms $\sim b^{2}e^{-b}$ are disfavored.

Based on this experience we have formulated some simple ansatzes for $f_{NP}$.

• •

Model 1: This ansatz uses the fact that the simplest even-$b$ function with exponent asymptotics is the hyperbolic cosine. The model reads

[TABLE]

where $\lambda_{1}[\text{GeV}]>0$ and $\lambda_{2}[\text{GeV}^{2}]>0$ are free parameters. The advantage of this model is its simplicity and independence of the Bjorken variable. The model 1 has a quadratic (Gaussian) behavior at small-$b$ $f_{NP}\sim e^{-\lambda_{2}b^{2}}$ and exponential behavior at large-$b$ $f_{NP}\sim e^{-\lambda_{1}b}$.

• •

Model 2: The model 2 reads

[TABLE]

where $\lambda_{1}[\text{GeV}]>0$ and $\lambda_{2}[\text{GeV}^{2}]>0$ are free parameters. In this model we attempt to incorporate the theoretical expectations on the $z$-dependence of $f_{NP}$. So, the model 2 has a $zb^{2}$-behavior at small-$b$ $f_{NP}\sim e^{-\lambda_{2}zb^{2}}$ and exponential behavior at large-$b$ $f_{NP}\sim e^{-\lambda_{1}b}$.

Both models have two parameters, which we include in the parameterization such that the parameter $\lambda_{1}[\text{GeV}]$ dictates the asymptotical behavior at large $b$. and the parameter $\lambda_{2}[\text{GeV}^{2}]$ gives the quadratic term. A priory, the parameter $\lambda_{1}$ should be of order of $m_{\pi}\sim 0.14\;\text{GeV}$, since it is the only natural scale of strong forces at large distances. The parameter $\lambda_{2}[\text{GeV}^{2}]$ roughly corresponds to the size of the leading power correction to small-$b$ OPE, see sec. 2.3. We can associate $\lambda_{2}$ with the scale $B$ as $\lambda_{2}\sim B^{-2}$. In ref. Scimemi:2016ffw we have estimated the size of this parameter in the large-$\beta_{0}$ approximation as $\lambda_{2}\sim 0.075\;\text{GeV}^{2}$.

Additionally, to the parameters $\lambda_{1,2}$ we have studied the parameter $g_{K}[\text{GeV}^{2}]>0$, which parametrizes the non-perturbative contribution to the rapidity evolution kernel $\mathcal{D}$ (see eq. (20)). The importance of this parameter is not clear from the literature. In ref. Scimemi:2016ffw we have estimated its size in the large-$\beta_{0}$ approximation as $0.01\pm 0.03\;\text{GeV}^{2}$, i.e. consistent with zero. Also, the fit of DAlesio:2014mrz shows a negligible influence of this parameter on the final results. Therefore, in the following we consider both possibilities $g_{K}=0$ and $g_{K}\neq 0$. In section 3.7, we demonstrate that the parameter $g_{K}$ is important at lower perturbative order, but its influence is negligible at NNLL/NNLO.

3.4 The domain of TMD factorization

The TMD factorization is restricted to the small-$q_{T}$ range. The size of the allowed $q_{T}$-region is a priory unknown. We have not found any phenomenological studies on this point but only some statement on the strong dependence of the fit on the $q_{T}$-window. A specific study on TeVatron Z-boson production data in ref. Becher:2011xn shows that the Y-term contribution is extremely marginal for $q_{T}<30$ GeV.

In order to make a qualitative study, we introduce the parameter $\delta_{T}$ and we consider all data points with $q_{T}<\delta_{T}Q$. The amount of data points which are allowed by such a restriction are shown in the table 10. In order to estimate the maximum value of $\delta_{T}$ we perform a series of fits with increasing values of $\delta_{T}$. Ideally, the $\chi^{2}/d.o.f.$ and the fitting parameters should be stable within and unstable outside of the allowed $\delta_{T}$ interval. In this way, considering the dependence on $\delta_{T}$ one should find an interval of $\delta_{T}$ for which the fit is not sensitive to the $Y$-term. This point indicates the region of TMD-factorization, and should not depend of the perturbative order.

We have performed such a test for high-energy data set with different one-parameter forms of $f_{NP}$. We have especially used the one parameter models to guarantee the absence of fine-tuning of the cross-section. For this reason we also exclude the E288 data, because it is impossible to describe high- and low-energy data with a single non-perturbative parameter. The result of the fits practically agrees for all tested models and orders. In fig. 3, we present some typical outcome of the fits.

In plots 3 one can see that all models reproduce the data very-well at very small $\delta_{T}$, which is expected since the TMD factorization is only valid at $q_{T}\ll Q$. Then the value of $\chi^{2}$ slightly grows but keeps less then one until $\delta_{T}=0.2$ and after this threshold it jumps to higher values. The next jump is seen at $\delta_{T}=0.25$. After $\delta_{T}=0.25$ the value of $\chi^{2}$ increases rapidly. We interpret this fact saying that at $\delta_{T}=0.2$ we become sensitive to $Y$-term, and at $\delta_{T}=0.25$ the $Y$-term starts to dominate the cross-section, i.e. we leave the domain of TMD factorization. We have found that the presented plots rather strongly depend on the set of pertubative scales. For some choice of these scales, one can obtain an ideally flat plateau of $\chi^{2}$ for $\delta_{T}\leqslant 0.2$. However, the values of the two important thresholds, namely, $\delta_{T}=0.2$ (where deviation form TMD factorization appears) and $\delta_{T}=0.25$ (where the TMD factorization is completely broken), are stable with perturbative scales.

As a result of these tests, in the following we use the data points with $q_{T}\lesssim 0.2\;Q$, or say $\delta_{T}=0.2$. The choice of $\delta_{T}$ that we make is consistent with Becher:2011xn . This range includes $163$ high-energy and $146$ low-energy data points (in total 309 data points). Comparing this number of points with the literature, we observe that, it is the largest set of points for Drell-Yan/Z-boson production used up to present in a simultaneous fit of TMDPDF (to our knowledge), which also has the largest considered range of energies from $(Q,\sqrt{s})=(4,19.4)\;\text{GeV}$ (from the E288 experiment) to $(Q,\sqrt{s})=(150,8000)\;\text{GeV}$ (from the ATLAS experiment).

3.5 Scale variations and theoretical uncertainties

The theoretical uncertainties of the perturbative inputs are tested by varying the perturbative scales around their central values, as it is discussed in sec. 2.4. The distribution of uncertainties through orders for a typical high energy experiment is shown in figs. 4 and 5, and for a typical low-energy experiments in fig. 6-7. The complete set of plots for every included experiment can be found in web .

The uncertainty associated with the TMD evolution factor is parameterized by the $c_{1}$-variation. This uncertainty drops down between NLL/NLO and NNLL/NLO orders, that is together with the increase of the perturbative order for $\mathcal{D}$ (see table 1). The size of the band is correlated with the energy of the process, that is, it is less significant for higher-energy experiments.

The uncertainty associated with the hard scale depends on the $c_{2}$-variation. This band is independent on $q_{T}$. This error is the main one at NLL/LO (which we do not present here), but becomes negligible at higher orders.

The uncertainty associated with the low-energy behavior of the evolution factor is parameterized by the $c_{3}$-variation. We have found that it significantly influences the shape of the cross-section and also it is rather large at small-$q_{T}$. As expected it is decreases going from NLL/NLO to NNLL/NNLO. At NNLL/NNLO it gives the main contribution to the uncertainty band for the cross-section.

The uncertainty associated with the small-$b$ matching of coefficients and PDFs is represented by the $c_{4}$-variation. It is the most interesting error because it checks the convergences of the $\zeta$-prescription. The corresponding error-band is larger at $q_{T}\to 0$, which corresponds to the contribution of large $\mathbf{L}_{\mu}$ (we remind that in $\zeta$-prescription, $\mathbf{L}_{\mu}$ grows unrestrictedly). The important observation is that the large uncertainty area significantly shrinks between NLL/NLO and NNLL/NNLO, although the NNLL/NNLO contains a higher power of $\mathbf{L}_{\mu}$. This shows a very good behavior of the $\zeta$-prescription. In total this error is dominant at NLL/NLO, but becomes smaller (although compatible) to the one coming from the $c_{3}$ variation at NNLL/NNLO.

The size of the theoretical error-band is significantly bigger at small-$Q$, as can be visually checked comparing figs. 4-5 to figs. 6-7. The uncertainties reduces when one increases the perturbative orders, both in high and low energy cases. However for the low energy case the error remains of order $\sim 20\%$ or higher even at NNLL/NNLO, which can be problematic for a precise description of these experiments. We additionally stress that at NLL/LO the uncertainties range from $50\%$ to $100\%$ and higher. This shows that this particular order has no prediction power, and should not be considered any serious for a well based extraction of TMDs. This is the main reason for excluding NLL/LO order from our analysis.

In order to provide a final definition of the theoretical error, we use all scale variations and we take the maximum deviation among them. We have found that our definition of uncertainties is close, as far as one can compare different theoretical expressions, to the common definition used e.g. in Catani:2015vma ; Bozzi:2008bb . In total, for the high-energy experiments we find that the theoretical uncertainty (at NNLO) is of the order $2-3$% at the peak. It grows to $\sim 5-6\%$ at maximum allowed $q_{T}$, and to $\sim 10\%$ at $q_{T}\to 0$. This value seems to be smaller (but comparable) to the typical values of uncertainties presented ResBos or DYRes. This is a definite positive point of the $\zeta$-prescription. Indeed, the main contribution (at high energies) to it comes from the $c_{3}$- and $c_{4}$-error-bands, which are controlled by $\zeta$-prescription. The $c_{4}$-band would be significantly larger in the presence of double-logarithms, which are absent due to the $\zeta$-prescription.

3.6 Normalization

As the TMD factorization approach describes the shape of the differential cross section only in a limited range of $q_{T}$, we need some extra input to normalize the curves. In order to compare with the data, we weight the differential cross-section by the total (or fiducial) cross-section. The values of the theory predictions for total cross-sections can be obtained from the studies of other groups. For example, one can use the DYNNLO code Catani:2007vq ; Catani:2009sm . Its predictions for the total cross-sections are presented in the tables 3-4. However, we found that such a strategy is unreliable, because even tiny disagreement in the normalization leads to huge effects in the $\chi^{2}$-minimization. This is especially important for LHC data sets, which have very small error-bands. Additionally, as we demonstrate later, the DYNNLO predictions are worse than that obtained using our normalization factors.

Therefore, to fit the high energy data set we introduce a normalization factor for each data set. This factor equals the partial integral over $q_{T}$ for experimental and theoretical cross-sections, and reads

[TABLE]

where $\Delta q_{T}$ is the size of the $q_{T}$ bin. In this way, we fit only the $q_{T}$-shape of cross-section, which is already very restricting, as we discussed in the previous section.

The values of $\mathcal{N}^{-1}$ resulting from the calculations are presented in table 11. It is clear that the deviation between the theory and experiment decreases with perturbative orders. For the majority of experiments (excluding the Z-boson production measured by ATLAS), we find a good agreement for the absolute value of the differential cross-section obtained from the data points and the TMD factorization. It is important that the values of $\mathcal{N}$ are very stable with respect to the change of non-perturbative model and to the scale variation. In particular, we do not present the error-band on the normalization values in the table 11, because they are smaller then the present precision.

The normalization of the data from E288 experiment is generally unknown. Most probably, the main source of discrepancy comes from the fiducial cuts made for E288 experiment, which cannot be restored nowadays. The small fiducial cuts do not seriously influence the $q_{T}$-shape of the differential cross-section, but can sizably decrease the total normalization. In our analysis, we change the common normalization of all E288 data points as

[TABLE]

This or close values have been used in different fits, see e.g. Landry:2002ix ; DAlesio:2014mrz . However, we do not seriously ground on it, e.g. we can switch to $0.85$ or $0.9$ without significant loss in $\chi^{2}$ (however, the value $1$ produces serious disagreement with our current input). One should take into account that the theoretical uncertainty at small$-Q$ is very large, see figures 6 and 7. It also implies that low-energy cross-sections are very sensitive to the choice of PDF set (in particular, our approximation of eq. (45) for nuclei PDF could be too crude). We have checked that the E288 data can be also fitted with $N_{E288}=1$ to the same values (or better) of $\chi^{2}$ by additional variation of $Q_{0}$ (similar to the fit made in Guzzi:2013aja ). Such an ambiguity represents a problem in the analysis of the low-energy data.

3.7 Results of the fits and TMD extraction

In this section, we present the results of the global fit for the complete data sets presented in sec. 3.1, which allows the extraction of the unpolarized TMDPDF. We have made two independent fits, with $g_{K}=0$ and with $g_{K}\neq 0$. The results of the $\chi^{2}$ minimization and the values of the extracted parameters are presented in tab. 12 and 13. The visual presentation is given in fig. 8 and 9.

We have estimated both statistical and theoretical errors on the fit parameters. The statistical errors are related to the uncertainty of the $\chi^{2}$-minimization and are induced by the experimental error-bands. The statistical errors have been estimated by the MINOS procedure of MINUIT package James:1994vla . The theoretical errors are related to the uncertainty of perturbation series. There is no common procedure for the estimation of the theoretical error. Therefore, we propose the method presented in the following.

The theoretical error is estimated by a set of independent fitting procedures for each variation of the scale constants $c_{1,2,3,4}\in[0.5,2]$, as discussed in sec. 2.5. In other words, we set, say, $c_{1}=2$ and perform the minimization of $\chi^{2}$. In this way, we obtain a new set of model parameters (and a new value of $\chi^{2}$). In total, we have 8 independent variations and hence have 8 values of parameters. The final theoretical error-band is given by the maximal positive and minimal deviations from the central value and the results are reported in tab. 14. A drawback of this procedure is the variation of a scale can lead to the serious increase in $\chi^{2}$. In other words changing the matching scales affects also the quality of the fit. In general, the size of the band for $\chi^{2}$ value represents the stability of the theoretical model, and they are also reported in tab. 14. One can see that the error for $\chi^{2}$ significantly drops with orders.

The values of the parameter $\lambda_{1}$, which parametrizes the asymptotics of TMDPDFs, extracted at different orders agree with each other within the error-band, that slightly reduces with the increase of the order. It has a natural size of the order of pion mass, $\lambda_{1}\sim 1.3m_{\pi}-2.3m_{\pi}$. The values of the parameter $\lambda_{2}$, which parameterizes the quadratic correction to the small-$b$ regime, are not so stable although the values at different orders are compatible within the errors. In particular, they have large error-bars at NNLL/NLO order. The behavior of $g_{K}$ is the most peculiar. It decreases with the increase of the perturbation order. Moreover, at NNLL (both /NLO and /NNLO) its error-band touches the zero. It can be interpreted as following: the parameter $g_{K}$ is very small (or even zero) but within the fit, it tends to compensate the missing higher perturbative orders of evolution exponent. We also observe that all extractions of $g_{K}$ agrees with the theoretical estimation $g_{K}=0.01\pm 0.03\;\text{GeV}^{2}$ made in Scimemi:2016ffw . One can see that both models produce very similar results both for $\chi^{2}$ and the parameters.

As expected the theoretical error is reduced with the increase of the perturbative order. In particular, the band on the value of $\chi^{2}$ is significantly smaller at NNLL/NNLO. The distribution of parameter values over perturbative orders presented in table 14 is typical. The variation of $c_{1}$ does not represent the main contribution to the error-band. It implies that the low-energy matching for the rapidity anomalous dimension is not so important (in comparison to other matchings), as typically expected.

The variation of $c_{2}$ is almost negligible. Here, however, we recall that $c_{2}$ influences only the common normalization factor, and thus the effect of its variation could be underestimated due to our fitting procedure. The variation of $c_{3}$ and $c_{4}$ produces the most part of the error-band and the strongest variation of $\chi^{2}$. At $g_{K}=0$ these variation are more-or-less equivalent. At $g_{K}\neq 0$ there is a clear pattern. In this case, the variation of $c_{3}$ gives the main error-band on $g_{K}$, while the variation of $c_{4}$ gives the main error-band on parameters $\lambda_{1,2}$. It is very natural since the variation of $c_{3}$ tests the low-energy normalization point of the evolution factor, and $c_{4}$ tests the uncertainties of perturbation determination of the TMDPDF.

In tab. 15 we present the distribution of values for $\chi^{2}$ among experiments. One can see that the most stringent constraints come from the Z-boson production data of ATLAS and D0 run2. This is due to the small experimental uncertainty of these data points. At the low-energy, the main tension is presented by the 4-5 GeV bins, while the rest are distributed more-or-less homogeneously. It probably indicates the influence of generic factorization violating terms. The plots of the theoretical curves (at NNLL/NNLO for model 1) and the data points for individual experiments are shown in figures 10-14. The plots for different models and at other orders can be found in web .

4 Conclusion

The unpolarized Drell-Yan process at small-$q_{T}$ offers the simplest application of the TMD factorization formalism, and as such it has been studied by many groups. In this work, we have revised the main points of the practical implementation of TMD factorization, and reveal some new aspects of the TMD phenomenology. Altogether it allows us to critically reanalyze the available Drell-Yan data and to extract consistently the unpolarized TMDPDFs, within some approximation. The primary aim of our analysis is to answer some general questions for the TMD approach such as: Up to which $q_{T}$ the TMD factorization works? What is the best asymptotical behavior of a TMD distributions? How convergent is TMD formalism at higher orders of perturbative expansion? The answers to these questions are naturally affected by the used prescriptions for the practical implementation of the TMD formalism. Even so, these important issues of TMD phenomenology are undiscussed in the literature or discussed very superficially. Implementing consistently the TMD factorization formalism, we are able to fit a large set of Drell-Yan data points which ranges from low ($Q=4\;\text{GeV}$) to high ($Q=116$-$150\;\text{GeV}$) dilepton invariant masses on a wide interval of center of mass energies and using a limited set of parameters (two for the non-perturbative part of TMDPDF and one for the non-perturbative part of the TMD evolution).

In this work, we have formulated and used the $\zeta$-prescription, which is one of the main new theoretical contributions of this article. The $\zeta$-prescription consists of a particular choice of the rapidity evolution scale $\zeta=\zeta_{\mu}$, which depends on $\mu$, $b$ and the parton flavor (quark or gluon). This choice corresponds to the equi-evolution line in the space of TMD scales, and thus a TMD distribution is $\mu$-independent along this line. As a consequence, all logarithms related to the TMD evolution, which are essentially double logarithms, are eliminated from the small-$b$ OPE. It significantly improves the perturbative convergence and the radius of convergences for the small-$b$ OPE. The value of $\zeta_{\mu}$ is dictated by the differential equation (28), which can be solved order-by-order in the perturbation theory. We stress that the $\zeta$-prescription does not strictly solve the problem of large-$b$ logarithms, which are still present in the matching coefficients. However, these logarithms are not related anymore to the TMD scales. Moreover, these logarithms are accompanied by the $x$-dependent coefficients which preserve the integral over $x$ in accordance with the probability interpretation of PDFs. Note, that $\zeta$-prescription is universal for all TMDs of the leading dynamical twist, due to the universality of TMD ultraviolet and rapidity renormalization factors. There are multiple possibilities to apply $\zeta$-prescription, see some discussion in appendix B.1. In this work, we have used the simplest one, which can certainly be improved. A further study of $\zeta$-prescription will be done elsewhere.

Within our implementation of TMD factorization and TMD distributions, which has a generic form, we have three independent perturbative series: one for hard matching, one for rapidity evolution, and one for small-$b$ matching. To defend the approach we provide the estimation of the perturbative uncertainty by variation of associated scales by factor 2, see sec. 3.5. We have considered three successive perturbative orders (from NLL/NLO to NNLL/NNLO, see table 1), and demonstrate that the theory uncertainties and the agreement with the data improve with the increase of the perturbative order. The agreement of the theory with the experiment resulting in our fit is a consequence of the $\zeta$-prescription to a large extent. The lowest possible combination of perturbative order, namely NLL/LO, produces very large theoretical error-bands and thus has been excluded from the present study.

Our analysis shows that data are very selective about the non-perturbative part of the TMDs and only well-behaved models can accommodate the fit. The best models for the non-perturbative part of TMDPDF that we have found are formulated in sec. 3.3. They have a common non-perturbative structure, namely

[TABLE]

where $f$ is the PDF, $C$ is the small-$b$ matching coefficient and $f_{NP}$ is a non-perturbative input. We have found that the best agreement with data is given when the function $f_{NP}$ behaves as

[TABLE]

We have considered two ansatzes which respect eq. (56), see sec. 3.3 and eq. (51)-(52). The models have different behavior in the intermediate $\bm{b}$ region, in particular, model 2 has $z$ dependence. Nonetheless, the models produce nearly identical values of $\chi^{2}$ and of parameters $\lambda_{1,2}$. It implies that the parameters $\lambda_{1,2}$ that largely determine the shape of TMDPDF have a precise physical meaning. The values of parameters are reasonable $\lambda_{1}\sim 1.5\;m_{\pi}$ and $\lambda_{2}\sim 0.5$ GeV2. We also study the influence of the parameter $g_{K}$, that parameterizes the non-perturbative part of TMD evolution. We have found that this parameter is significant at lower order (in our case NLL/NLO) and less significant at higher orders. Moreover, at NNLL/NNLO it is compatible with zero within the error-bars. We supplement our extraction with the estimation of the theoretical and statistical errorbars. The obtained theoretical uncertainty on the extracted parameters is notably larger in comparison to the experimental one (see fig. 8 and 9). It can indicate the purity of the models (say, the independence on the flavor, or insufficient parameterization), or the general weakness of the theory.

Another aspect that we point out, is the practical limitation of TMD factorization. To make the discussion quantitative we introduce the parameter $\delta_{T}$, which is the highest allowed ratio $q_{T}/Q$ accounted in the fit. Clearly, at very low $\delta_{T}$ the TMD formalism should perfectly work, e.g. provide small values of $\chi^{2}$-distribution. Our expectation is that within the domain of the TMD-factorization the value of $\chi^{2}/d.o.f.$ is largely constant and starts to grow outside of this domain. Indeed, for the best models, the observed picture agrees with the expectation. In this way, we have shown that TMD factorization as it is, i.e. in the absence of $Y$-term, is capable of describing the data with $q_{T}\lesssim 0.2\ Q$, i.e. $\delta_{T}=0.2$. With some risk, one can prolong it to $\delta_{T}=0.25$. After $\delta_{T}=0.25$ the TMD factorization loses any agreement with the experiment. This analysis is unique, or at least we do not know any analog in the literature.

The fit and the plots of the data has been done with the help of arTeMiDe, version 1.1, available at web . This is a code package for the numerical evaluation of TMD distributions and related cross-sections. It has a flexible structure and allows to consider an arbitrary combination of perturbative orders up to NNLO for coefficient functions and N3LO for evolution factors. In the current version, it evaluates only unpolarized TMDPDFs, but we expect to update it for polarized cases and TMD fragmentation functions, as well as, to include the $Y$-term, in the future.

Acknowledgements.

We thank Miguel G. Echevarria for the help on the initial stages of this work, and V. Bertone for useful comments on arTeMiDe code. I.S. also thanks Stefano Melis for discussions about some technical points. I.S. is supported by the Spanish MECD grant FPA2014-53375-C2-2-P and FPA2016-75654-C2-2-P and the group UPARCOS.

Appendix A Hard coefficient function

The hard coefficient function $|C_{V}(\mu,Q)|^{2}$ can be obtained from the expression for the quark form factor. At NNLO, it can found in Kramer:1986sg ; Matsuura:1988sm . Here, we present the combined expression in Drell-Yan kinematic

[TABLE]

Appendix B $\zeta$-prescription and expressions for coefficient functions

In this appendix, we elaborate the details of the $\zeta$-prescription and expression for the coefficient function. Throughout the paper, we use the following notation for logarithms

[TABLE]

For convenience we introduce the following notation for the perturbative coefficient of anomalous dimensions

[TABLE]

The LO terms are

[TABLE]

where $\beta_{0}=\frac{11}{3}C_{A}-\frac{2}{3}N_{f}$ is the leading order QCD $\beta$-function. The higher order terms can be found e.g. in ref. Echevarria:2016scs .

B.1 Derivation of $\zeta$-value

The $\zeta$-prescription is defined as a curve in $(\mu,\zeta)$-plane along which a TMD distribution has no evolution. In other words, at $\zeta=\zeta_{\mu}$

[TABLE]

This equation can be rewritten as

[TABLE]

Using the explicit expressions for the anomalous dimensions eq. (11), and introducing the intermediate function $\mathbf{l}_{\zeta_{\mu}}=f(\mathbf{L}_{\mu})$ we obtain

[TABLE]

Solving this equation order-by-order in perturbation theory, we obtain

[TABLE]

where $c_{1,2,3}^{f}$ are integration constants. To derive this expression we have used that the $\mathbf{L}_{\mu}$-dependence of rapidity anomalous dimension at NNLO has the form

[TABLE]

where the constants $d^{(2,0)}_{f}$ and $d^{(3,0)}_{f}$ can be found in Vladimirov:2016dll .

The integration constants $c_{1,2,3}$ that appears in eq. (B.1) are to be fixed by additional conditions, which would correspond to a selection of a particular curve among the family of equi-evolution curves in $(\mu,\zeta)$-plane. In our current analysis, we set all constants $c_{i}$ to zero for simplicity. It corresponds to the curve that passes though the point $(\mu,\zeta)=(0,0)$. We leave for the future a dedicated study of boundary condition and its influence on the phenomenology. Thus, the explicit NLO expression for $\mathbf{l}_{\zeta_{\mu}}$ for quark TMDPDF reads

[TABLE]

The explicit NLO expression for $\mathbf{l}_{\zeta_{\mu}}$ for gluon TMDPDF reads

[TABLE]

The expression for the $\zeta_{\mu}$ then reads

[TABLE]

B.2 Scale dependence and logarithmic part of coefficient function

The small-$b$ coefficient function satisfies the pair of equations

[TABLE]

where $P(x)$ is the DGLAP kernel of the PDF evolution, and $\otimes$ denotes the Mellin convolution in the variable $x$. Using these equations one finds the logarithmic part of the coefficient function. At NNLO the expression for the coefficient function reads

[TABLE]

where we omit the argument $(x)$ of DGLAP kernel $P(x)=\sum_{n}a_{s}^{n}P^{(n)}(x)$ and the finite part of the coefficient function $C^{(n,0)}(x)$, and $\delta_{ff^{\prime}}(\bar{x})=\delta_{ff^{\prime}}\delta(\bar{x})$.

Substituting the NLO expression for $\mathbf{l}_{\zeta}$, eq.  (B.1) into the coefficient function, eq. (B.2) we obtain at NNLO

[TABLE]

Note, that despite the fact that the solution for $\zeta$-prescription in eq. (B.1) has inverse powers of $\mathbf{L}_{\mu}$, the coefficient function has not. It is easy to check that this expression convoluted with PDF is renormalization-invariant,

[TABLE]

B.3 Expression for NNLO coefficient function in $\zeta$-prescription

The NNLO coefficient functions are cumbersome structures, which contain logarithms and polylogarithms of order 2 and 3 and their straight numerical evaluation is costly. To speed up the evaluation of convolutions within arTeMiDe, we have used an approximate expression for the coefficient function. A similar method for higher-order expressions has been suggested in ref. vanNeerven:1999ca and it is widely used in NNLO+ phenomenology of PDFs. We parameterize the NNLO coefficient function by 17 terms

[TABLE]

Here, the coefficients $A$ represent the singular at $x\to 1$ terms, and are evaluated exactly. The coefficients $B$ represent singular at $x\to 0$ terms, and also evaluated exactly. The coefficients $c$ represent interpolation functions between asymptotics. These coefficients are fit numerically. The relative precision of the approximation is $\sim 10^{-3}$. The convolution integral receives the main numerical contributions at singular points, while the rest are minor corrections. So, we find that the relative accuracy of the convolution is better then $10^{-6}$, which is far beyond any currently needed accuracy. The values of coefficients A, B, and c are given in the tables 16, 18 and 18.

[FIGURE:]

[FIGURE:]