The impact of the LHC Z-boson transverse momentum data on PDF determinations

TL;DR

This paper assesses how recent precise LHC measurements of Z-boson transverse momentum distributions, combined with NNLO QCD predictions, improve the accuracy of parton distribution functions, especially reducing uncertainties in Higgs production cross sections.

Contribution

It demonstrates the significant impact of incorporating latest $p_T^Z$ data into NNLO PDF fits, reducing uncertainties in key parton distributions and Higgs production predictions.

Findings

Uncertainty in gluon-fusion Higgs cross section reduced by ~30%.

Inclusion of $p_T^Z$ data improves PDF precision.

Central values of PDFs remain stable after data inclusion.

Abstract

The LHC has recently released precise measurements of the transverse momentum distribution of the Z-boson that provide a unique constraint on the structure of the proton. Theoretical developments now allow the prediction of these observables through next-to-next-to-leading order (NNLO) in perturbative QCD. In this work we study the impact of incorporating these latest advances into a determination of parton distribution functions (PDFs) through NNLO including the recent ATLAS and CMS 7 TeV and 8 TeV $p_T^Z$ data. We investigate the consistency of these measurements in a global fit to the available data and quantify the impact of including the $p_T^Z$ distributions on the PDFs. The inclusion of these new data sets significantly reduces the uncertainties on select parton distributions and the corresponding parton-parton luminosities. In particular, we find that the $p_T^Z$ data ultimately leads to a reduction of the PDF uncertainty on the gluon-fusion and vector-boson fusion Higgs production cross sections by about 30%, while keeping the central values nearly unchanged.

Tables

Table 1. Table 1: χ 2 superscript 𝜒 2 \chi^{2} per degree of freedom for the normalized ATLAS 7 TeV p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} on-peak distributions in the separate rapidity bins before their inclusion in the fit. As input PDFs we use the NNPDF3.0 set with α S ​ ( M Z ) = 0118 subscript 𝛼 𝑆 subscript 𝑀 𝑍 0118 \alpha_{S}(M_{Z})=0118 . The computation is done at NLO, NNLO and NNLO QCD + NLO EW, with (N)NLO PDF set for (N)NLO computations. Results for CT14, MMHT2014 and ABMP16 are also shown.

Bin	Order	$N_{\mathrm{dat}}$	${\chi }_{d.o.f.}^2$ (NN30)	${\chi }_{d.o.f.}^2$(CT14)	${\chi }_{d.o.f.}^2$ (MMHT14)	${\chi }_{d.o.f.}^2$(ABMP16)
	NLO	14	10	21	9.2	n.a.
	NNLO	14	2.2	3.8	4.3	11
	NNLO+EW	14	1.3	2.3	2.6	9.1
	NLO	14	13	18	12	n.a.
	NNLO	14	5.6	8.2	9.3	15.
	NNLO+EW	14	3.9	6.0	6.8	12.
	NLO	14	7.0	7.1	6.0	n.a.
	NNLO	14	7.0	8.2	8.7	11.
	NNLO+EW	14	5.9	7.1	7.5	9.5
All bins	NLO	42	9.9	15	9.1	n.a.
	NNLO	42	4.9	6.7	7.4	13.
	NNLO+EW	42	3.7	5.2	5.6	12.

Table 2. Table 2: Same as Table  1 for the ATLAS 8 TeV p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} distributions in the low and high invariant-mass bins before their inclusion in the fit.

Bin	Order	$N_{\mathrm{dat}}$	${\chi }_{d.o.f.}^2$ (NN30)	${\chi }_{d.o.f.}^2$ (CT14)	${\chi }_{d.o.f.}^2$ (MMHT14)	${\chi }_{d.o.f.}^2$ (ABMP16)
GeV	NLO	8	2.5	2.2	1.9	n.a.
	NNLO	8	1.7	1.6	1.7	1.1
	NNLO+EW	8	1.2	1.1	1.2	0.9
GeV	NLO	8	2.3	2.6	2.3	n.a
	NNLO	8	1.3	1.2	1.2	2.1
	NNLO+EW	8	1.1	1.1	1.0	2.1
GeV	NLO	8	1.4	1.3	1.0	n.a
	NNLO	8	1.1	1.0	1.0	0.7
	NNLO+EW	8	0.8	0.8	0.8	0.6
GeV	NLO	10	1.9	1.9	1.5	n.a
	NNLO	10	0.8	0.8	0.8	1.0
	NNLO+EW	10	0.4	0.4	0.4	0.5
GeV	NLO	10	2.3	2.1	1.6	n.a
	NNLO	10	1.2	1.0	1.0	1.3
	NNLO+EW	10	1.2	1.1	1.0	1.5
All bins	NLO	44	1.3	1.2	1.1	n.a
	NNLO	44	1.0	0.9	0.9	1.1
	NNLO+EW	44	1.0	1.0	1.0	1.4

Table 3. Table 3: Same as Table  1 for the ATLAS 8 TeV p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} on-peak distributions in the separate rapidity bins before their inclusion in the fit.

Bin	Order	$N_{\mathrm{dat}}$	${\chi }_{d.o.f.}^2$ (NN30)	${\chi }_{d.o.f.}^2$ (CT14)	${\chi }_{d.o.f.}^2$ (MMHT14)	${\chi }_{d.o.f.}^2$ (ABMP16)
	NLO	10	4.0	3.2	2.4	n.a.
	NNLO	10	2.7	2.7	2.6	2.7
	NNLO+EW	10	3.4	3.2	3.1	5.4
	NLO	10	5.6	4.6	3.8	n.a.
	NNLO	10	5.4	5.2	5.3	3.3
	NNLO+EW	10	4.0	3.9	3.7	3.8
	NLO	10	5.8	3.8	3.0	n.a.
	NNLO	10	4.7	4.0	4.3	2.1
	NNLO+EW	10	2.3	2.0	1.9	1.7
	NLO	10	4.5	3.2	2.5	n.a.
	NNLO	10	5.1	4.0	4.6	3.0
	NNLO+EW	10	3.3	2.6	2.7	2.5
	NLO	10	4.4	3.2	2.4	n.a.
	NNLO	10	5.4	4.3	5.0	3.7
	NNLO+EW	10	3.9	3.2	3.4	3.0
	NLO	10	4.1	3.2	2.4	n.a.
	NNLO	10	3.4	3.1	3.3	3.2
	NNLO+EW	10	2.6	2.3	2.4	2.5
All bins	NLO	60	3.4	2.0	1.9	n.a.
	NNLO	60	4.5	4.0	4.4	2.6
	NNLO+EW	60	2.9	2.5	2.7	2.4

Table 4. Table 4: Same as Table  1 for the CMS 8 TeV p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} distributions in the separate rapidity bins before their inclusion in the fit.

Bin	Order	$N_{\mathrm{dat}}$	${\chi }_{d.o.f.}^2$ (NN30)	${\chi }_{d.o.f.}^2$ (CT14)	${\chi }_{d.o.f.}^2$ (MMHT14)
GeV	NLO	9	3.1	2.6	2.2
	NNLO	9	2.2	2.4	2.3
	NNLO+EW	9	1.4	1.4	1.3
GeV	NLO	9	2.4	1.9	1.5
	NNLO	9	1.4	1.4	1.3
	NNLO+EW	9	1.9	1.9	1.7
GeV	NLO	9	2.3	2.0	1.6
	NNLO	9	1.6	1.4	1.4
	NNLO+EW	9	1.2	1.2	1.0
GeV	NLO	9	1.6	2.0	1.9
	NNLO	9	1.4	1.4	1.3
	NNLO+EW	9	2.8	3.1	2.8
All bins	NLO	36	3.3	3.2	3.3
	NNLO	36	3.5	3.5	3.7
	NNLO+EW	36	3.9	4.0	4.0

Table 5. Table 5: Overview of fits run with HERA-only as a baseline. For each fit, we indicate which measurements from ATLAS and CMS has been included, whether an uncorrelated uncertainty has been added to the χ 2 superscript 𝜒 2 \chi^{2} (in brackets unless it is set to 0).

	(a)	(b)	(c)	(d)	(e)	(f)	(g)	(h)	(i)	(j)
HERA	y	y	y	y	y	y	y	y	y	y
ATLAS8TEV	n	y(1%)	n	y(1%)	y(0.5%)	n	y(0.5%)	y	n	y
CMS8TEV	n	n	y(1%)	y(1%)	n	y(0.5%)	y(0.5%)	n	y	y

Table 6. Table 6: Fully correlated χ d . o . f . 2 subscript superscript 𝜒 2 formulae-sequence d o f \chi^{2}_{{\rm d.o.f.}} for the fits described in Table  5 . The numbers in brackets correspond to the χ 2 superscript 𝜒 2 \chi^{2} for experiments which are not fitted. In particular the ATLAS 7 TeV data are not fitted in any of these fits. The total χ 2 superscript 𝜒 2 \chi^{2} is computed over all baseline HERA data the included p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} distributions.

fit id	extra $\Delta$	${\chi }_{\mathrm{ATLAS7tev}}^2$	${\chi }_{\mathrm{ATLAS8tev},m}^2$	${\chi }_{\mathrm{ATLAS8tev},y}^2$	${\chi }_{\mathrm{CMS8tev}}^2$	${\chi }_{\mathrm{tot}}^2$
(a)	1%	(21.8)	(1.00)	(1.56)	(1.55)	1.168
(b)	1%	(19.6)	0.91	0.70	(1.61)	1.146
(c)	1%	(16.2)	(1.04)	(1.56)	1.21	1.176
(d)	1%	(18.0)	0.90	0.77	1.42	1.156
(a)	0.5%	(27.6)	(1.10)	(2.83)	(2.46)	1.168
(e)	0.5%	(23.0)	0.99	1.05	(3.01)	1.168
(f)	0.5%	(20.5)	(1.13)	(3.15)	1.91	1.198
(g)	0.5%	(21.4)	0.99	1.29	2.44	1.207
(a)	no	(30.6)	(1.15)	(4.65)	(3.46)	1.168
(h)	no	(25.5)	1.02	1.66	(4.79)	1.193
(i)	no	(19.5)	(1.28)	(5.44)	2.51	1.225
(j)	no	(24.5)	1.03	2.09	3.59	1.251

Table 7. Table 7: Overview of fits run with HERA-only as a baseline including the normalized ATLAS 7 TeV along with the other data sets. For each fit, we indicate which measurements from ATLAS and CMS has been included, whether an uncorrelated uncertainty has been added to the χ 2 superscript 𝜒 2 \chi^{2} (in brackets unless it is set to 0).

	(a)	(k)	(l)	(m)	(n)	(o)	(p)
HERA	y	y	y	y	y	y	y
ATLAS7TEV	n	y(1%)	y(1%)	y(0.5%)	y(0.5%)	y	y
ATLAS8TEV	n	n	y(1%)	n	y(0.5%)	n	y
CMS8TEV	n	n	y(1%)	n	y(0.5%)	n	y

Table 8. Table 8: Fully correlated χ d . o . f . 2 subscript superscript 𝜒 2 formulae-sequence d o f \chi^{2}_{{\rm d.o.f.}} for the fits described in Table  7 . The numbers in brackets correspond to the χ 2 superscript 𝜒 2 \chi^{2} for experiments which are not fitted. The total χ 2 superscript 𝜒 2 \chi^{2} is computed over all baseline HERA data the included p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} distributions.

fit id	extra $\Delta$	${\chi }_{\mathrm{ATLAS7tev}}^2$	${\chi }_{\mathrm{ATLAS8tev},m}^2$	${\chi }_{\mathrm{ATLAS8tev},y}^2$	${\chi }_{\mathrm{CMS8tev}}^2$	${\chi }_{\mathrm{tot}}^2$
(a)	1%	(21.8)	(1.00)	(1.56)	(1.55)	1.168
(k)	1%	1.39	(1.39)	(2.04)	(1.41)	1.176
(l)	1%	1.64	1.05	1.17	1.27	1.171
(a)	0.5%	(27.6)	(1.10)	(2.83)	(2.46)	1.168
(m)	0.5%	1.58	(1.54)	(3.36)	(2.11)	1.186
(n)	0.5%	2.13	1.18	1.98	2.21	1.253
(a)	no	(30.6)	(1.15)	(4.65)	(3.46)	1.168
(o)	no	1.74	(1.69)	(4.79)	(3.06)	1.185
(p)	no	2.35	1.24	2.81	3.19	1.301

Table 9. Table 9: Fully correlated χ 2 superscript 𝜒 2 \chi^{2} for the experiments in the HERA + p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} 8 TeV fit.

fit id	extra $\Delta$	Theory	${\chi }_{\mathrm{ATLAS7tev}}^2$	${\chi }_{\mathrm{ATLAS8tev},m}^2$	${\chi }_{\mathrm{ATLAS8tev},y}^2$	${\chi }_{\mathrm{CMS8tev}}^2$	${\chi }_{\mathrm{tot}}^2$
(e)	1%	NNLO	(18)	0.90	0.77	1.42	1.156
(q)	1%	NNLO+EW	(16)	1.00	0.87	1.72	1.182

Table 10. Table 10: Fully-correlated χ 2 superscript 𝜒 2 \chi^{2} per degree of freedom when the p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} data is added to the global fits. The numbers in brackets correspond to the χ 2 superscript 𝜒 2 \chi^{2} for experiments which are not fitted. The total χ 2 superscript 𝜒 2 \chi^{2} is computed over all data in the baseline fit and the included p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} distributions. We have labeled our slightly-modified NNPDF3.0 global baseline as NN30red in the table below.

fit	${\chi }_{\mathrm{ATLAS7tev}}^2$	${\chi }_{\mathrm{ATLAS8tev},\mathrm{mdist}}^2$	${\chi }_{\mathrm{ATLAS8tev},\mathrm{ydist}}^2$	${\chi }_{\mathrm{CMS8tev}}^2$	${\chi }_{\mathrm{tot}}^2$
NN30red	(6.93)	(0.98)	(1.06)	(1.41)	1.17677
NN30red + 8 TeV	(7.87)	0.96	0.88	1.32	1.17690

Table 11. Table 11: Predictions for the Higgs cross sections in 13 TeV p ​ p 𝑝 𝑝 pp collisions before and after inclusion of the p T Z superscript subscript 𝑝 𝑇 𝑍 p_{T}^{Z} data in the global fits. The indicated errors are the PDF errors computed according to the NNPDF prescription.

	Before $p_T^Z$ data	After $p_T^Z$ data
${\sigma }_{gg\to H}$ [pb]	$48.22\pm 0.89$ (1.8%)	$48.61\pm 0.61$ (1.3%)
${\sigma }_{\mathrm{VBF}}$ [pb]	$3.92\pm 0.06$ (1.5%)	$3.96\pm 0.04$ (1.0%)

Equations

$$\begin{array}{rl}{M}_Z& =91.1876\text{GeV},{\Gamma }_Z=2.4925\text{GeV},\\ G_{\mu }& =1.11639\times 10^{-5}{\text{GeV}}^2,M_W=80.398\text{GeV}.\end{array}$$

$${\mu }_R={\mu }_F=\sqrt{(p_T^Z{)}^2+M_{ll}^2}.$$

$$\frac{d{\sigma }^{(2,1)}}{d{p}_T^Z}=\frac{d{\sigma }^{(2,0)}}{d{p}_T^Z}\times \frac{d{\sigma }^{(0,1)}/d{p}_T^Z}{d{\sigma }^{(0,0)}/d{p}_T^Z}.$$

$$p_T^Z>30\mathrm{GeV}$$

$$|y_Z|<1.6(\mathrm{CMS}\mathrm{only}).$$

$$p_T^Z<150\mathrm{GeV}(\mathrm{ATLAS}8\mathrm{TeV},\mathrm{peak}\mathrm{region})$$

$$p_T^Z<170\mathrm{GeV}(\mathrm{CMS}8\mathrm{TeV}).$$

$$p_T^Z>30\mathrm{GeV}$$

$$p_T^Z<500\mathrm{GeV},$$

Full text

††arxiv: 1705.00343

1 Introduction

The production of a $Z$-boson that subsequently decays into a pair of leptons is a benchmark Standard Model (SM) process at the Large Hadron Collider (LHC). Thanks to its large production rate and clean experimental signature, it can be be measured very accurately by the LHC experiments. It can also be calculated to high accuracy within the Standard Model, with the first prediction to next-to-next-to-leading order (NNLO) in the strong coupling constant appearing more than two decades ago Hamberg:1990np , and predictions for differential cross sections appearing over one decade ago Anastasiou:2003ds ; Melnikov:2006kv ; Catani:2009sm ; Gavin:2010az . This combination of precise experimental data and highly-developed theory allows this process to be used to determine quantities of fundamental importance to our understanding of high-energy phenomena, such as parton distribution functions (PDFs).

Among the many distributions in $Z$-boson production that have been measured, the transverse momentum ($p_{T}$) distribution stands out as an especially interesting one. First of all, the $Z$-boson $p_{T}$ spectrum is sensitive to the gluon and the light-quark PDFs in the not-so-well constrained intermediate Bjorken-$x$ region, which makes it a promising observable for constraining these distributions. The fact that the Higgs production cross section at the LHC is also sensitive to the same PDF combinations in the same region of Bjorken-$x$, makes the measurement of this process of direct importance to the search for beyond-the-SM phenomena in the Higgs sector. Second, the transverse momentum spectrum of the $Z$-boson is sensitive to both soft QCD radiation (at small $p_{T}$) and to large electroweak (EW) Sudakov logarithms (at large $p_{T}$). Given that PDF fits typically rely on fixed-order perturbative QCD, it is interesting to test how well fixed-order QCD predictions can describe this data. This has direct impact on which range of data can be included into PDF fits.

The potential for $p_{T}^{Z}$ measurements to provide valuable constraints on PDF determinations has been considered previously, both on general grounds Malik:2013kba ; Rolph:2015hoa , and when considering a recent measurement performed by the CMS collaboration Chatterjee:2016sxt . Both of these studies, which are based on NLO QCD, show the potential of these measurements. At the same time, they also stress the importance of including the full NNLO QCD corrections to the $Z$-boson transverse momentum distribution in order to fully exploit the constraining power of the data.

In present global PDF determinations, the gluon distribution at medium and large $x$ is primarily constrained by the inclusive-jet $p_{T}$ spectrum measurements. The full NNLO prediction for this observable has been recently calculated in the leading-color approximation Currie:2016bfm , but results have not yet been made available for all jet data sets included in PDF fits. This deficiency motivates the study of other cross sections known to NNLO for this purpose, such as the $Z$-boson $p_{T}$ spectrum, or top-pair production. For the latter, studies have appeared that explored in great detail the possibility of making use of the total cross section Beneke:2012wb ; Czakon:2013tha and more recently of the differential distribution Czakon:2016olj measurements. In particular, it was shown that differential distributions from top-pair production provide significant constraints on the large-$x$ gluon that are comparable to those obtained from inclusive jet production data.

The importance of including NNLO corrections is especially clear in the case of the $Z$-boson transverse momentum distribution given the recent experimental progress in measuring this observable. The data sets from the 7 and 8 TeV LHC runs from both ATLAS and CMS feature percent-level experimental errors, clearly requiring predictions beyond NLO in order to achieve a comparable theoretical precision.

It is our intent in this manuscript to investigate the inclusion of the $p_{T}^{Z}$ data from the LHC into a global PDF fit. We perform this study in a framework based on the NNPDF3.0 global analysis Ball:2014uwa . The data sets we consider in our work are the 7 TeV measurement of the $Z$-boson $p_{T}$ by the ATLAS collaboration Aad:2014xaa , and the 8 TeV measurements from both ATLAS and CMS Aad:2015auj ; Khachatryan:2015oaa . These data sets include doubly-differential distributions in both the rapidity and invariant mass of the lepton pair coming from the $Z$-boson decay. Our theoretical predictions are based on the NNLO QCD calculation of Ref. Boughezal:2015ded . We also study the impact of including approximate NLO electroweak corrections, as described later in the text. The major findings of our study are summarized below.

• •

The inclusion of the NNLO QCD corrections generally improves the agreement of theory with the experimental data. This conclusion is consistent with previous observations Ridder:2016nkl ; Gehrmann-DeRidder:2016jns . The simultaneous inclusion of the NLO electroweak contributions together with NNLO QCD, done here for the first time, further improves the data/theory agreement at high $p_{T}$.

• •

The experimental errors, particularly in the higher-luminosity 8 TeV measurements from ATLAS and CMS, have dropped to the percent level. With the data becoming so precise, a very careful accounting of both experimental and theoretical errors is needed. We observe difficulties in fitting the data without the introduction of an additional uncorrelated error in the fit. This can come from a combination of Monte Carlo integration errors on the theory calculation, residual theoretical uncertainties in the prediction, or from underestimated experimental errors. We expect this issue to become increasingly prevalent in future PDF fits as data becomes more precise.

• •

We observe difficulties when attempting to simultaneously fit the 7 TeV and 8 TeV LHC data. The ATLAS 7 TeV data is provided only in terms of normalized distributions, while the 8 TeV measurements are also provided as absolute, unnormalized distributions. The normalization to the fiducial cross section performed for the ATLAS 7 TeV data introduces correlations between the low-$p_{T}^{Z}$ bins and the $p_{T}^{Z}>30$ GeV region to which we must restrict our fit due to the appearance of large logarithms in the low-$p_{T}^{Z}$ region that require resummation. The covariance matrix provided for the whole data set then turns out to be incorrect when used for fitting a subset of the data. This prevents us from consistently including the ATLAS 7 TeV data in the fit. To validate this hypothesis, in Sec. 5.3 we perform a fit including the normalized ATLAS 8 TeV data rather than the unnormalized ones but, in analogy to what is done for the 7 TeV data, using the covariance matrix provided for the whole data set, and explore the differences in the fit results. It would be interesting to revisit this issue if the unnormalized data for the 7 TeV measurement were released or if the experimental covariance matrix for the $p_{T}^{Z}>30$ GeV region was available. Attempting to include resummed predictions for the low-$p_{T}^{Z}$ region is also possible, although this would introduce additional theoretical uncertainties.

• •

When adding the 8 TeV LHC $Z$-boson $p_{T}$ data to the global NNPDF3.0-like fit, we observe a significant decrease of the gluon PDF uncertainty in the Bjorken-$x$ region $10^{-3}$ to $10^{-1}$ as well as a reduction of the uncertainty for light quarks. This leads to a reduction of the PDF uncertainty on the gluon-fusion and Vector Boson Fusion (VBF) Higgs boson cross section of 30%, while the central value prediction for both processes increases by roughly 1%.

Our manuscript is organized as follows. In Section 2 we describe the experimental measurements of $p_{T}^{Z}$ that we include in our fit. We also present the baseline fits that do not include these data that we use to assess their impact. In Section 3 we discuss the details of the theoretical calculation and settings that we use in the fit. We give a comparison of theory with the $p_{T}^{Z}$ data in Section 4. We discuss the agreement observed upon using NLO QCD, NNLO QCD or a combined NNLO QCD + NLO EW prediction, and also consider several different global PDF sets. Our fit to the $p_{T}^{Z}$ data and several baseline fits is described in Section 5. We briefly discuss the phenomenological impact of the new fits on the Higgs cross section in Section 6. Finally, we conclude in Section 7.

2 Descriptions of the experimental data and fit settings

In this Section we first discuss the features of the available experimental measurements. We then describe the methodology and settings of our fit to the parton distribution functions including these data.

2.1 $p_{T}^{Z}$ measurements from the LHC

In this work we consider the most recent differential cross section measurements of the $Z$-boson transverse momentum spectrum from ATLAS Aad:2014xaa ; Aad:2015auj and CMS Khachatryan:2015oaa , both with$\sqrt{s}=7$ TeV and $\sqrt{s}=8$ TeV .

The ATLAS measurement of the $Z$-boson transverse momentum spectrum at the centre-of-mass energy of $\sqrt{s}$=7 TeV Aad:2014xaa is performed in the $Z\to e^{+}e^{-}$ and $Z\to\mu^{+}\mu^{-}$ channels, using data based on an integrated luminosity of 4.7 ${\rm fb}^{-1}$. The results from each channel are combined for transverse momenta up to 800 GeV. The measurement is provided both inclusive in the $Z$-boson rapidity up to 2.4, and separated into three rapidity bins: $0.0<|y_{Z}|<1.0$, $1.0<|y_{Z}|<2.0$ and $2.0<|y_{Z}|<2.4$. In order to maximize the constraints on PDFs, we include the data in the three exclusive rapidity bins in our analysis. In the experimental paper only the normalized distributions are provided. The measurement is very accurate, with statistical and systematical uncertainties below 1% in all $p_{T}^{Z}$ bins up to 150 GeV and for central rapidities ($|y_{Z}|<2.0$), and about 3% for the largest rapidity bin.

In the ATLAS measurement at $\sqrt{s}=8$ TeV Aad:2015auj , the transverse momentum distribution is based on the full 8 TeV data set, with 20.3 fb*-1* of integrated luminosity. Measurements are performed in the electron-pair and muon-pair channels and then combined. Compared to the 7 TeV measurement Aad:2014xaa , this measurement has higher statistics and an improved control of experimental systematics. Measurements are performed in six invariant mass bins: four bins at low invariant mass below the $Z$-peak, one on-peak invariant mass bin, and one bin at high invariant mass above the $Z$-peak, reaching up to $M_{ll}=$ 150 GeV. Results for the off-peak bins are provided in one inclusive rapidity bin ($0.0<|y_{Z}|<2.4$), while the $Z$-peak measurement results are given both inclusive over the whole rapidity range $0.0<|y_{Z}|<2.4$ and separated in six rapidity bins $0.0<|y_{Z}|<0.4$, $0.4<|y_{Z}|<0.8$, $0.8<|y_{Z}|<1.2$, $1.2<|y_{Z}|<1.6$, $1.6<|y_{Z}|<2.0$ and $2.0<|y_{Z}|<2.4$. Again, in order to maximize the constraints on PDF, we include the on-peak exclusive rapidity bins in our analysis.

The measurement by the CMS collaboration at the center-of-mass energy $\sqrt{s}=8$ TeV Khachatryan:2015oaa is performed differentially in five rapidity bins: $0.0<|y_{Z}|<0.4$, $0.4<|y_{Z}|<0.8$, $0.8<|y_{Z}|<1.2$, $1.2<|y_{Z}|<1.6$ and $1.6<|y_{Z}|<2.0$. The analysis uses the data sample of $pp$ collisions collected with the CMS detector at the LHC in 2012, which corresponds to an integrated luminosity of 19.7 fb*-1*. The $Z$-boson is identified via its decay to a pair of muons. We only include the measurements exclusive in the muon rapidities up to $|y_{Z}|=1.6$, given that the data in the highest rapidity bin display a significant incompatibility with respect to the corresponding ATLAS measurement. We leave this issue to further investigation by the experimental collaborations.

2.2 Settings for the PDF analysis

The PDF fits presented in this work are based on the NNPDF3.0 global analysis Ball:2014uwa framework. As in the NNPDF3.0 fit, both PDF evolution and DIS structure functions are evaluated in the fit using the public APFEL library Bertone:2013vaa ; Carrazza:2014gfa ; Bertone:2016lga , with heavy-quark structure functions computed in the FONLL-C general-mass variable-flavor-number scheme Forte:2010ta with pole masses and with up to $n_{f}$=5 active flavors. The DGLAP evolution equations are solved up to NNLO using a truncated solution, and the input parametrization scale is taken to be $Q_{0}=1$ GeV. The strong coupling $\alpha_{s}$ is set to $\alpha_{s}(M_{Z})=0.118$, in accordance with the PDG average Agashe:2014kda . The charm and bottom PDFs are generated perturbatively from light quarks and gluons and the value of the heavy-quark masses are set to $m_{c}=1.51$ GeV and $m_{b}=4.92$ GeV, corresponding to the values recommended by the Higgs Cross Section Working Group Dittmaier:2011ti . Note that these values are different from the ones used in NNPDF3.0, which were instead set to the PDG value of the $\overline{\text{MS}}$ masses. These values will be used in the forthcoming NNPDF3.1 release nnpdf31 . The dependence of the fit on the values of the heavy quark masses is moderate, and in particular is negligible for the observables under consideration.

In the analysis performed in this work, we consider two baseline data sets. One consists of all available HERA deep inelastic scattering (DIS) data. An important difference with respect to the NNPDF3.0 HERA-only baseline is that the HERA inclusive structure function data, which in NNPDF3.0 were separated into the HERA-II measurements from H1 and ZEUS Aaron:2012qi ; Collaboration:2010ry ; Abramowicz:2012bx , have been replaced by the HERA legacy combination Abramowicz:2015mha that has become available recently. This data is supplemented by the combined measurements of the charm production cross section $\sigma_{cc}^{\rm red}$ Abramowicz:1900rp , and the H1 and ZEUS measurement of the bottom structure function $F_{2}^{b}(x,Q^{2})$ Aaron:2009af ; Abramowicz:2014zub . The other baseline, a global one, contains all data mentioned in the paragraph above along with the other data analyzed in the NNPDF3.0 global fit: fixed-target neutral-current DIS structure functions from NMC Arneodo:1996kd ; Arneodo:1996qe , BCDMS bcdms1 ; bcdms2 , and SLAC Whitlow:1991uw ; charged-current structure functions from CHORUS inclusive neutrino DIS Onengut:2005kv and from NuTeV dimuon production data Goncharov:2001qe ; MasonPhD ; fixed-target E605 Moreno:1990sf and E866 Webb:2003ps ; Webb:2003bj ; Towell:2001nh DY production data; Tevatron collider data including the CDF Aaltonen:2010zza and D0 Abazov:2007jy $Z$ rapidity distributions; and LHC collider data including ATLAS Aad:2011dm ; Aad:2013iua ; Aad:2011fp , CMS Chatrchyan:2012xt ; Chatrchyan:2013mza ; Chatrchyan:2013uja ; CMSDY and LHCb Aaij:2012vn ; Aaij:2012mda vector boson production measurements, adding up to a total of $N_{\rm dat}=3530$ data points. A further difference from the global baseline (on top of the use of the HERA combined measurements) is that in order to ensure a consistent treatment of NNLO corrections, we exclude jet production measurements Abulencia:2007ez ; Aad:2011fc ; Chatrchyan:2012bja ; Aad:2013lpa from the global baseline data set. Only the leading color approximation has been made available at NNLO for this process Currie:2016bfm and $K$-factors are not yet available for all data sets included in global PDF determinations.

3 Description of the theoretical calculation

For our study we have calculated the $Z$-boson transverse momentum distribution through next-to-next-to-leading order in perturbative QCD. This computation uses a recent result for the related process of $Z$-boson in association with a jet Boughezal:2015ded ; Boughezal:2016isb based on the $N$-jettiness subtraction scheme for NNLO calculations Boughezal:2015dva ; Boughezal:2015eha ; Gaunt:2015pea . As the $Z$-boson obtains its transverse momentum through recoil against jets, these two processes are identical in perturbation theory as long as the cuts on the final-state jets are relaxed sufficiently so that the entire hadronic phase space is integrated over for the $Z$-boson $p_{T}$ values under consideration. Since at most three jets can recoil against the $Z$-boson at NNLO, we take the lower cut on the leading-jet $p_{T}$ to be less than $1/3$ times the lowest $Z$-boson $p_{T}$ included in our study. We have confirmed that our predictions are not sensitive to the exact choice of this jet cut. We furthermore remove completely any constraints on the pseudorapidities of final-state jets. We note that the low transverse momentum region of $Z$-boson production requires the resummation of large logarithmic corrections of the form $(\alpha_{s}\text{ln}^{2}(M_{Z}/p_{T}^{Z}))^{n}$ to all orders in perturbation theory for a proper theoretical description. This resummation is not present in our fixed-order calculation. We consequently restrict our attention to the region $p_{T}^{Z}>30$ GeV when comparing our predictions to the experimental data. In Sect. 5.3 we study the effect of raising the cut on $p_{T}^{Z}$ to 50 GeV and observe that results are stable upon the choice of the $p_{T}^{Z}$ cut.

We compare the theoretical predictions against both the unnormalized $p_{T}$ spectra provided by the 8 TeV ATLAS and CMS measurements, and also to the distributions normalized to the fiducial $Z$-boson production cross section provided by the 7 TeV ATLAS measurement. For the normalized distributions we compute the fiducial $Z$-boson production cross section using the $N$-jettiness subtraction scheme as implemented in MCVM v8.0 Boughezal:2016wmq . We cross-check this result against FEWZ Melnikov:2006kv ; Gavin:2010az . For the normalized distributions we do not expand the ratio in the strong coupling constant; i.e., we compute both the numerator and denominator through relative ${\cal O}(\alpha_{s}^{2})$.

We make the following choices for the electroweak input parameters in our calculation:

[TABLE]

We use the $G_{\mu}$ electroweak renormalization scheme. All other couplings are therefore derived using the parameters above, including the electromagnetic couplings and the weak mixing angle. We choose the following dynamical scale choices for both the renormalization and factorization scales:

[TABLE]

Here, $M_{ll}$ denotes the invariant mass of the final-state lepton pair. We note that our calculation includes both the $Z$-boson production and decay to lepton pairs, the contribution from virtual photons, as well as all interferences. The residual theoretical uncertainty on the prediction as estimated by independently varying $\mu_{R}$ and $\mu_{F}$ around this central value is at the few-percent level.

As we will see later it is also important when describing the high-$p_{T}$ data to include the effect of electroweak perturbative corrections. The exact NLO electroweak corrections to the $Z$-boson transverse momentum spectrum, including the leptonic decay of the $Z$ boson, are known in the literature Denner:2011vu ; Hollik:2015pja . However, they are not publicly available in the form of a numerical code. To account for their effect in our calculation we instead utilize the approximate expressions presented in Refs. Kuhn:2004em ; Kuhn:2005az . These include all one-loop weak corrections up to terms power-suppressed by the ratio $M_{Z}^{2}/((p_{T}^{Z})^{2}+M_{Z}^{2})$, and additionally the leading two-loop electroweak Sudakov logarithms. These expressions are strictly valid only after inclusive integration over the final-state lepton phase space; we apply them also to the cross sections with fiducial cuts on the leptons. For the $Z$-boson peak region in 8 TeV collisions we have checked that these approximations reproduce the numerical magnitude of the exact electroweak corrections to 2% or better in the high-$p_{T}^{Z}$ range where the EW effects become relevant. Since the electroweak corrections themselves do not exceed 10% for the entire region studied, this furnishes an approximation to the distributions we study that is good to the few-per-mille level or better, which is sufficient for our purposes 111The exact electroweak corrections were used in the ATLAS analysis of 8 TeV data; we thank A. Denner, S. Dittmaier, and A. Mueck for providing us with these results to cross-check our approximations. When we study normalized distributions, the NLO electroweak corrections to the fiducial $Z$-boson cross section are obtained from FEWZ Li:2012wna . To combine the electroweak and QCD corrections we assume that the two effects factorize, leading to a multiplicative combination. Denoting the differential cross sections at the $m$-th order in the strong coupling constant relative to the LO result and the $n$-th order in the QED coupling constant relative to the LO result as $d\sigma^{(m,n)}$, we assume that

[TABLE]

This factorization of the electroweak and QCD corrections is supported by a calculation of the dominant mixed ${\cal O}(\alpha\alpha_{s})$ corrections in the resonance region Dittmaier:2015rxo .

The experimental errors in the $Z$-peak region have reached an unprecedented level for a high-energy collider experiment, approaching the per-mille level over two orders of magnitude in transverse momentum. Numerous effects that were previously not relevant may now come into play, and it is worthwhile to briefly discuss the theoretical issues that arise when attempting to reach this precision. While we can not currently address these issues, they should be kept in mind when considering these data sets.

• •

The uncalculated N3LO perturbative QCD corrections may be needed to further improve the agreement between theory and experimental data. As we will see in a later section the theoretical predictions are generally below the experimental measurements. The inclusion of the NNLO corrections greatly improves the agreement between theory and experiment, but one may expect a further increase from the N3LO corrections.

• •

The electroweak corrections become important for $p_{T}^{Z}\sim 100$ GeV, reaching the percent level at this point and continuing to grow as $p_{T}^{Z}$ is increased. While we assume that the electroweak and QCD corrections factorize, this assumption should be addressed, particularly in the high-$p_{T}^{Z}$ region. Non-factorizing ${\cal O}(\alpha\alpha_{s})$ effects could conceivably affect the cross section at the percent level.

• •

Finally, at this level of precision non-perturbative QCD effects that shift the $p_{T}^{Z}$ distribution must be considered. 222We thank G. Salam for raising this issue, and for discussions.. Since the $Z$-boson transverse momentum distribution is generated by recoil against a final-state jet, there may be linear non-perturbative power correction of the form $\Lambda_{\text{QCD}}/p_{T}^{Z}$. Simple Monte Carlo estimates indicate that this could reach the half-per-cent level Gavintalk .

4 Comparison of theory with LHC data

In this Section we compare the theoretical predictions for the $p_{T}^{Z}$ spectrum to the experimental measurements described in Sect. 2. We assess the impact of NNLO QCD and NLO electroweak corrections and quantify the agreement between data and theory by computing the fully-correlated $\chi^{2}$ for each of the experiments that we include in our analysis using as input the most recent public releases of four PDF determinations: ABMP16 Alekhin:2017kpj CT14 Dulat:2015mca , MMHT2014 Harland-Lang:2014zoa and NNPDF3.0 Ball:2014uwa .

In Fig. 1 we compare the NLO and NNLO predictions to the experimental measurements performed by the 7 TeV ATLAS measurements, described in Ref. Aad:2014xaa , after imposing the additional cut of $p_{T}^{Z}>30$ GeV discussed earlier. We also include the NLO EW corrections as described in Section 3. All three rapidity bins measured by ATLAS are shown.

We observe that the NNLO corrections significantly increase the NLO predictions, bringing them closer to the measured values of the distribution. The NNLO corrections are approximately constant as a function of $p_{T}^{Z}$. The EW corrections become significant only for the last three $p_{T}^{Z}$ bins. The quantitative agreement with the theory is summarized in table 1, in which the fully-correlated $\chi^{2}$ is provided, for each bin separately and for the three bins together.

For MMHT2014, CT14 and NNPDF3.0 the agreement is improved for central rapidities after the inclusion of NNLO QCD corrections, with a further improvement observed upon including NLO electroweak corrections. For ABMP16 only the NNLO fit is available, so in this case we can only test that the agreement is improved upon adding electroweak corrections. In the highest rapidity bin this improvement is only observed for NNPDF3.0. The CT14 $\chi^{2}_{\rm d.o.f.}$ remains unchanged after including NNLO QCD+NLO electroweak, while the result for MMHT2014 becomes slightly worse. For all PDF sets the $\chi^{2}_{\rm d.o.f.}$ is much larger than one, indicating a poor agreement between theory and data (before the fit) even after including higher-order corrections.

In Figs. 2 and 3 a similar comparison is performed for the off $Z$-peak bins of the 8 TeV ATLAS measurement Aad:2015auj . The NNLO QCD corrections again provide a positive shift of the NLO result that is approximately independent of $p_{T}^{Z}$, with NLO electroweak corrections causing a approximatively constant upwards (downwards) shift for the bins below (above) the $Z$-peak. While the NNLO predictions are in better agreement with the data than the NLO ones, the data are again higher than the theoretical predictions. The quantitative comparison of the NNPDF3.0, MMHT2014, CT14 and ABMP16 PDF sets using the $\chi^{2}_{\rm d.o.f.}$ defined previously is shown in Table 2. In all cases an improvement is seen upon inclusion of the NNLO QCD corrections, while the incorporation of the NLO electroweak corrections as well further improves the agreement in all individual bins below the $Z$ peak.

We next consider the 8 TeV ATLAS data on the $Z$-peak divided into rapidity bins. The comparisons of NLO, NNLO and NNLO+EW theory with data are shown in Figs. 4 and 5, while the $\chi^{2}_{\rm d.o.f.}$ results for NNPDF3.0, MMHT2014, CT14 and ABMP16 are shown in Table 3. The general trends observed in this comparison are similar to those seen in the ATLAS 7 TeV comparison and in the comparison of the invariant mass binned data: the NNLO corrections increase the NLO predictions by an amount almost independent of $p_{T}^{Z}$, bringing theory closer to data. The quantitative comparison of $\chi^{2}_{\rm d.o.f.}$ in Table 3 reveals that NNLO improves upon the NLO description in four of the six rapidity bins for NNPDF3.0, while NNLO+EW improves upon NLO for all six bins. For CT14 NNLO+EW improves upon NLO for five of the six bins, while for MMHT the improvement is only observed for two bins. One reason that the inclusion of the NNLO corrections does not improve the theory/data agreement as significantly as in the other data sets is because the experimental error in this case is very small, and is dominated by the correlated systematic error. Even if NNLO reduces the normalization difference between theory and experiment, remaining shape differences between the predictions and data prevent a large improvement in $\chi^{2}_{\rm d.o.f.}$ from being obtained. This issue will arise again when we attempt to add this data set to the PDF fit.

Finally, in Figs. 6 and 7, we show the comparison of the various theoretical predictions with the CMS 8 TeV data divided into rapidity bins Khachatryan:2015oaa . The $\chi^{2}_{\rm d.o.f.}$ is shown in Table 4. As discussed when describing the data in Sect. 2, we focus on the region $|y_{Z}|<1.6$. Including NNLO corrections improves the agreement between theory and data in all four rapidity bins, while adding NLO EW corrections further improves the comparison in all but the highest rapidity bins. We note that the CMS relative errors are larger than those found by ATLAS, and the issues seen in the $\chi^{2}_{\rm d.o.f.}$ comparison are not as pronounced as for the ATLAS 8 TeV data set. Interestingly, even though each individual rapidity bin is improved upon including NNLO, the $\chi^{2}_{\rm d.o.f.}$ combining all bins is slightly worsened at NNLO, again showing the impact of the correlated uncertainties when attempting to describe these very precise data sets. Fitting the data modifies the PDF shape, thus significantly improving the data description.

5 Inclusion of the $p_{T}^{Z}$ distribution in PDF fits

In this Section we first look at the correlation between the measured distributions and the various PDF combinations, which provides a first intuition for what parton distributions and at what value of $x$ we should expect to observe the largest impact when including these data in the fit. We then add each data set separately to a DIS HERA-only fit to determine basic compatibility of different data sets and to assess the impact of including EW corrections. Finally, we perform a fit adding $p_{T}^{Z}$ data to a global data set to estimate the impact of including these data in a realistic PDF determination.

5.1 Correlations between PDFs and $p_{T}^{Z}$ measurements

To determine the specific PDFs and regions in $x$ for which the $Z$-boson transverse momentum distribution measurements from ATLAS and CMS provide the most stringent constraints we study the correlation coefficient as a function of $x$ ($\rho(x)$), between PDFs at a given scale $Q$ and each bin of the measurements included in the present analysis. In Figure 8 we plot the correlations, computed using the SMPDF code Carrazza:2016htc , of the gluon, up-quark and down-quark distributions with the lowest invariant mass bin of the ATLAS 8 TeV measurement, and with the on-peak 8 TeV measurement of ATLAS, for the lowest rapidity bin. Each line corresponds to one $p_{T}^{Z}$ bin. These are representative examples, the pattern of correlations found for the other measurements is similar. We observe a strong correlation between the gluon distribution in the region $x\approx 10^{-3}-10^{-2}$ with the $p_{T}^{Z}$ measurements, with the correlation coefficient reaching nearly 90%. Slightly weaker correlations of approximately 60% are found for the up-quark and down-quark distributions. These plots make it clear that these data sets have a strong potential to improve our knowledge of PDFs in the $10^{-3}-10^{-2}$ region. The largest $p_{T}^{Z}$ bins are correlated to the $10^{-2}-10^{-1}$ region, thus an increase in the experimental statistics in that region would provide a stronger constraint also in the large-$x$ region.

5.2 Impact of the $p_{T}^{Z}$ data on a DIS HERA-only fit

We begin by assessing the quality of a fit to the HERA DIS data upon inclusion of the available $p_{T}^{Z}$ data at 8 TeV. The inclusion of the normalized ATLAS 7 TeV data is problematic and we discuss it separately in Sect. 5.3. We perform several fits that add the individual ATLAS and CMS data sets to HERA separately and in various combinations. As discussed in previous Sections we impose the following cuts on the $p_{T}^{Z}$ data:

[TABLE]

These constraints leave us with 60 data points for the ATLAS 8 TeV doubly-differential distributions in rapidity and $p_{T}$ on the $Z$-peak, 44 data points for the ATLAS 8 TeV doubly-differential distributions in the dilepton invariant mass and $p_{T}$, and 36 data points for the CMS 8 TeV doubly-differential distributions in rapidity and $p_{T}$ on the $Z$-peak. Additionally, we consider fits using pure NNLO QCD theory and fits with NNLO QCD and NLO EW corrections combined. In the pure NNLO fits we remove the $p_{T}^{Z}$ bins for which the EW corrections are larger than the sum quadrature of the statistical and uncorrelated systematic uncertainty of that data point to avoid fitting EW effects. This imposes the additional constraints

[TABLE]

These cuts reduce the number of data points to 48 for the ATLAS 8 TeV doubly-differential distributions in rapidity and $p_{T}$ on the $Z$-peak, 44 data points for the ATLAS 8 TeV doubly-differential distributions in invariant mass and $p_{T}$, and 28 data points for the CMS 8 TeV doubly-differential distributions in rapidity and $p_{T}$ in the $Z$-peak region.

Since we have considered numerous combinations of the available data and several different settings, we begin by summarizing the fits in Table 5. These are labelled (a)-(j). Our baseline fit with only HERA data is labelled (a). Fits (b) and (c) add individually the ATLAS 8 TeV data and CMS 8 TeV data sets. Fit (d) adds all 8 TeV data sets. A new feature we find necessary in our analysis is the inclusion of an additional uncorrelated uncertainty. This uncertainty can be due to a combination of Monte-Carlo integration uncertainties on the computationally expensive NNLO theoretical calculations, to residual theoretical uncertainties not accounted for in our fit, and to possible underestimated experimental errors. The need and approximate size of this contribution to the uncertainty can be inferred from an analysis based on modelling the NLO and NNLO theoretical predictions and their fluctuations along the lines of the one described for inclusive jet production in Carrazza:2017bjw . The addition of this new effect is needed to obtain a good $\chi^{2}$ in our fit, as shown later in this section. To study the stability of our fit with respect to this uncertainty we consider the values 0%, 0.5%, and 1%. Fits (b)-(d) use a 1% uncorrelated uncertainty, while fits (e)-(g) use 0.5%. This uncertainty is removed in fits (h)-(j). We will see later that the fitted PDFs are insensitive to the value of this parameter.

The results of fits (a)-(j) are summarized in Table 6. For each fit the $\chi^{2}$ per degree of freedom ($\chi^{2}_{{\rm d.o.f.}}$) of the experiments included in the fit, and of the prediction for the observables not included in the fit (in brackets), are displayed. The additional uncorrelated uncertainty added to the fit is denoted by $\Delta$. We have repeated the baseline HERA-only fit (a) at the beginning of each Table section for ease of comparison. A few things are apparent from the table.

• •

The addition of $\Delta$ improves the description of the ATLAS 8 TeV on-peak and CMS 8 TeV data. The $\chi^{2}_{{\rm d.o.f.}}$ decreases from 1.66 to 0.77 for the ATLAS 8 TeV set and from 2.51 to 1.21 for the CMS 8 TeV set as $\Delta$ is changed from 0% to 1% in the baseline fit. This effect is less noticeable for the invariant-mass binned ATLAS data due to the slightly larger errors for this set.

• •

Comparing fit (b) (where only the ATLAS 8 TeV data is fit along with HERA) to fit (c) (where only the CMS 8 TeV data is fit together with HERA) shows that the ATLAS 8 TeV data is slightly more consistent with HERA than CMS. The $\chi^{2}_{{\rm d.o.f.}}$ is below one for the ATLAS sets in fit (b) after including them in the fit, while it is at 1.21 in (c) when CMS is combined with HERA.

• •

Fit (d) shows that it is possible to obtain a reasonably good fit of ATLAS 8 TeV data, CMS 8 TeV data, and HERA with the inclusion of a $\Delta=1\%$ additional uncorrelated uncertainty. Reducing this uncertainty to 0.5% in fit (g) leads to a noticeably worse description of the CMS data. Both the CMS and on-peak ATLAS 8 TeV data sets get a worse $\chi^{2}_{{\rm d.o.f.}}$ if $\Delta$ is removed completely, as in fit (j).

• •

It is clear from the table that the ATLAS 7 TeV measurement is inconsistent with the other data sets. We discuss this further in Section 5.3.

We now study the implications of these fits for the PDF sets. All plots have been done by using the on-line interface of APFEL Bertone:2013vaa . We consider the gluon and the singlet-quark combination. To avoid too large a proliferation of plots we focus on the $\Delta=1\%$ and $\Delta=0\%$ cases. In Fig. 9 we display the impact of the inclusion of these data on the gluon and singlet-quark PDFs by adding them with an additional uncertainty $\Delta=1\%$. As can be seen from the upper left panel of of Fig. 9, including either the ATLAS 8 TeV and CMS 8 TeV data sets leads to a gluon consistent with the HERA result but with a slightly smaller uncertainty. The upper right panel shows that HERA+8 TeV gives a gluon similar to HERA-only but with a significantly smaller uncertainty for $x>10^{-3}$.

The situation for the singlet-quark distribution is similar. However the ATLAS and CMS data seem to pull in slightly different directions, the former preferring a harder singlet in the $x=10^{-1}$ region, as it can be observed in the lower-left panel. The lower-right panel shows that the ATLAS data have a stronger pull in the fit and that the simultaneous inclusion of the ATLAS and CMS data at 8 TeV leads to a significantly reduced uncertainty.

The effects of the $\Delta=1\%$ fits on the down-quark and up-quark distributions is similar to the effect on the singlet and thus is not shown here: the PDF errors when HERA and the 8 TeV data sets are simultaneously fit decreases significantly for both the up and down distributions.

In Fig. 10 we show the results for the PDFs assuming no additional uncertainty, $\Delta=0\%$. The observed patterns of PDF shifts when 8 TeV data sets are included is very similar to those seen for $\Delta=1\%$, with only small differences in the estimated PDF errors in certain $x$ regions.

5.3 Normalized versus unnormalized distributions

In this Section we focus on the inclusion of the normalized ATLAS 7 TeV data and give details on the tension we observe with the 8 TeV data. We consider a NNLO fit, applying the following cuts

[TABLE]

where the latter is motivated by the fact that in the last $p_{T}^{Z}$ bin the EW corrections are larger than the sum in quadrature of the statistical and uncorrelated systematic uncertainties of the data. We are left then with 39 data points for the ATLAS 7 TeV distribution.

We summarize the fits in Table 7. These are labelled (k)-(p). The baseline is the same as the one presented in the previous section. Fits (k), (m) and (o) add individually the ATLAS 7 TeV data by adding an uncorrelated uncertainty of 1%, 0.5% and none respectively. Fits (l), (n) and (p) add them along with the unnormalized ATLAS and CMS data at 8 TeV with an extra uncorrelated uncertainty of 1%, 0.5% and none respectively.

The results of fits (k)-(p) are summarized in Table 8. For each fit the $\chi^{2}$ per degree of freedom ($\chi^{2}_{{\rm d.o.f.}}$) of the experiments included in the fit, and of the prediction for the observables not included in the fit (in brackets), are displayed. The additional uncorrelated uncertainty added to the fit is denoted by $\Delta$. Again, we have repeated the baseline HERA-only fit (a) at the beginning of each Table section for ease of comparison.

A few things are apparent from the table.

• •

The ATLAS 7 TeV data is inconsistent with the HERA-only fit, with a $\chi^{2}_{{\rm d.o.f.}}$ over 20 regardless of the $\Delta$ chosen. A primary reason for this is that the ATLAS 7 TeV data is normalized to the fiducial cross section in each rapidity bin, while the 8 TeV data sets are unnormalized. The normalization performed for the ATLAS 7 TeV data introduces correlations between the low-$p_{T}^{Z}$ bins and the $p_{T}^{Z}>30$ GeV region to which we must restrict our fit due to the theoretical considerations discussed earlier. Due to this cut on the data the covariance matrix provided by the experiments for the whole data set cannot be used to consistently include the 7 TeV data in the fit. It would be interesting to revisit this issue if the unnormalized data became available.

• •

Studying fits (l), (n) and (p) shows that it is hard to simultaneously fit the ATLAS 7 TeV data with the 8 TeV data sets. In table 6 we observed that fitting the 8 TeV data leads to a $\chi^{2}_{{\rm d.o.f.}}$ of 18 for the ATLAS 7 TeV data in fit (d). In table 8 we see that the $\chi^{2}_{{\rm d.o.f.}}$ of the 8 TeV data deteriorates when we attempt to include the 7 TeV too.

We now study the implications of these fits for the PDF sets. We consider the gluon, up-quark and down-quark distributions and focus on the $\Delta=1\%$ case only, as we have seen that PDFs remain basically unchanged upon a reduction of $\Delta$. In Fig. 11 we display the impact of the inclusion of these data on the gluon, up and down quark PDFs by adding them with an additional uncertainty $\Delta=1\%$. An important feature of these plots is the difference between the impact of the ATLAS 7 TeV data on the gluon, compared to the impact of the 8 TeV data sets. As can be seen from the upper left panel of Fig. 11, including either the ATLAS 8 TeV and CMS 8 TeV data sets leads to a gluon consistent with the HERA result but with a slightly smaller uncertainty. Adding the ATLAS 7 TeV data leads to an increased gluon distribution for $x>5\cdot 10^{-3}$. The upper right panel shows that HERA+8 TeV gives a gluon similar to HERA-only but with a significantly smaller uncertainty for $x>10^{-2}$. Attempting to fit both 7 TeV and ATLAS 8 TeV data leads to an increased uncertainty, which is barely visible. The tension present between the ATLAS 7 TeV data, and the combined HERA+8 TeV data observed for the gluon PDF is also observed for the up and down distributions. The middle right panel shows that the error on the up-quark PDF is greatly increased for $x\approx 10^{-3}$ when we attempt to simultaneously fit all data. The reason for this can be seen from the left middle panel. The ATLAS 7 TeV data prefers a peak in the up-quark distribution at this value. In contrast, the upper right panel shows a decrease in the PDF error when HERA and the 8 TeV data sets are simultaneously fit. A similar pattern is observed for the down-quark distribution, as is shown in the lower two panels of Fig. 11.

In order to confirm that the origin of the anomalous behaviour of PDFs upon the inclusion of the 7 TeV data is due to the fact that they are normalized, we perform an additional fit including the precise ATLAS 8 TeV normalized data in the $Z$-peak region in a HERA-only fit. As far as the quality of the fit is concerned, we observe that these data are harder to fit than the 7 TeV ones, as the $\chi^{2}_{\rm d.o.f.}$ ranges from 9 (for a fit with $\Delta=0$%) to 2.1 (for a fit with $\Delta=1$%). As far as PDFs are concerned, in Fig. 12 it is apparent that, while the inclusion of the on-peak ATLAS 8 TeV unnormalized data reduces the uncertainty of the gluon and up-quark distributions, the inclusion of the on-peak ATLAS 8 TeV normalized data inflates their uncertainties, thus pointing to their inconsistency with respect to the baseline.

5.4 Impact of NLO EW corrections

Another interesting aspect that we can investigate is the impact of electroweak corrections on the obtained PDFs. To probe this we perform fits to the HERA and 8 TeV data sets, with NNLO QCD corrections and both with and without EW corrections. We recall that in the pure NNLO QCD fit we remove bins where the EW corrections are larger than the combined uncorrelated uncertainty, as explained previously. We first display the gluon, singlet, down-quark and up-quark distributions with and without EW corrections in Fig. 13. The EW corrections have a small but noticeable effect on the PDFs, lowering both the gluon and singlet distributions in the intermediate-$x$ regions. The $\chi^{2}_{\rm d.o.f.}$ is shown in Table 9. The quality of the fit deteriorates slightly upon including EW corrections. This results primarily not because EW corrections worsen the agreement between theory and data, but because with EW corrections included we are able to include additional high-$p_{T}^{Z}$ bins in the fit that were excluded in the pure NNLO QCD fit, and these bins are slightly more discrepant than the lower-$p_{T}^{Z}$ ones. The agreement with the 7 TeV data is marginally improved upon including EW corrections, although it is still inconsistent with the HERA+8 TeV combined fit.

5.5 Impact of the $p_{T}^{Z}$ data on a global fit

Having investigated the impact of the LHC $p_{T}^{Z}$ data in a fit consisting of only HERA data, which allowed us to consider several aspects of this new data in detail, we turn to their inclusion in a global fit of the available measurements. We follow the NNPDF3.0 analysis with the modifications explained in Section 2.2. We set the additional uncorrelated error to $\Delta=1\%$, and, having established that we cannot consistently include the normalized 7 TeV data in a PDF fit, we only add the unnormalized 8 TeV data to the global baseline. The results for the $\chi^{2}$ per degree of freedom of each fit is shown in Table 10. The $\chi^{2}_{\rm d.o.f}$ of the fitted $p_{T}^{Z}$ distributions reveals a mild tension between the CMS and ATLAS data sets, with $\chi^{2}_{\rm d.o.f}$ of the CMS set reaching 1.32, while the ATLAS 8 TeV sets give a $\chi^{2}_{\rm d.o.f}$ below one. We notice that when including the 8 TeV data the $\chi^{2}_{\rm d.o.f}$ of the (not-fitted) ATLAS 7 TeV data deteriorates.

In Fig. 14 we display the agreement of the NNLO predictions and the data before and after the fit. We observe that the agreement improves and uncertainties shrink.

In Fig. 15 and 16 we show the impact of the precise 8 TeV $p_{T}^{Z}$ data on the various PDFs determined from the global fit of the available data. The observed shifts of the PDFs are similar to those seen in the HERA-only fit. The reduction of the uncertainty is milder but still significant. The new PDFs obtained after including the 8 TeV $p_{T}^{Z}$ data are consistent with those found in the baseline.

It is interesting to compare our results with those presented in Czakon:2016olj , in which a similar baseline was used and the impact of including top-pair production differential distributions in PDF fits was studied in detail for the first time. The gluon is pulled in the same direction by both data sets, thus displaying a perfect compatibility between these two complementary measurements. The inclusion of the $p_{T}^{Z}$ data decreases the uncertainties on the gluon PDF more than the top-pair data in the intermediate-$x$ region between $10^{-3}$ and $10^{-2}$. The impact of the top-pair data is much stronger for $x>10^{-2}$. This result follows the correlation patterns presented in Section 5.1 for $p_{T}^{Z}$ and in Czakon:2016olj for top-quark differential distributions, from which it is clear that the latter are strongly correlated with the gluon in the large-$x$ region, while the former are mostly correlated with the gluon (and slightly less with the light-quark distributions) in the intermediate-$x$ region. Given that these two observables provide such strong and complementary constraints, we expect that their impact in a joint fit will be stronger than the impact of the jet data, which were traditionally thought to be the best probe of the gluon in the intermediate and large-$x$ regions.

To conclude, we explore the stability of our results upon increasing the $p_{T}^{Z}$ cut from 30 GeV to 50 GeV. Both the gluon and singlet central values are very stable, with uncertainties that are larger when a larger $p_{T}^{Z}$ cut is used. We note that the number of $p_{T}^{Z}$ data points in the fit decreases from 48 to 40 for the ATLAS 8 TeV on-peak data, from 44 to 36 for the ATLAS 8 TeV off-peak data and from 28 to 24 for the CMS 8 TeV on-peak data. Thus an increase in the PDF uncertainty when the cut is raised is expected. Everything else is consistent with expectations.

6 Phenomenological implications

Having derived a new global fit of PDFs with the 8 TeV $p_{T}^{Z}$ data included, it is interesting to investigate the impact of these new measurements on quantities of phenomenological interest.

Parton luminosities directly show the impact of the inclusion of a given data set on the computation of processes. A comparison of the 13 TeV parton-parton luminosities before the $p_{T}^{Z}$ data, and after including the unnormalized 8 TeV data, is presented in Fig. 18. The uncertainties significantly decrease in all three luminosities, while their central values remain nearly the same as before.

Furthermore, we present below the 13 TeV predictions for both the gluon-fusion Higgs production cross section and the VBF Higgs production cross section before and after the inclusion of the $p_{T}^{Z}$ data in our global baseline fit. For the gluon-fusion production cross section we set $m_{H}=125$ GeV and $\mu_{R}=\mu_{F}=m_{H}/2$ and use the code ggHiggs v3.5 Bonvini:2016frm to compute the result through N3LO in QCD perturbation theory Anastasiou:2016cez . The result below includes no charm or bottom quarks running in the loop, and no quark mass effects beyond leading order. The impact on the Higgs production cross section uncertainties is significant. The error on the gluon-fusion production cross section is reduced by 30%, following the corresponding improvement in the gluon-gluon-luminosity observed in Fig. 18. The central value is increased by only 1%, indicating consistency with the cross section obtained using the previous global fit. For Higgs production in Vector Boson Fusion we compute the total cross section to N3LO in QCD using the proVBFH-inclusive code provbfh based on the computation presented in Cacciari:2015jma ; Dreyer:2016oyx .

7 Conclusions

In this manuscript we have included for the first time the precision $p_{T}^{Z}$ measurements from the LHC into a global fit of parton distribution functions to next-to-next-to-leading order in QCD. This result is made possible by the recent theoretical predictions of this process to the necessary order. We have performed a detailed study of the impact of various perturbative corrections, including higher-order QCD and electroweak corrections, on the agreement between theory and data. To asses in detail the impact of these new data we have tested the effect of adding them to several baseline fits, including a DIS HERA-only PDF determination and a global fit with settings closely following those of NNPDF3.0.

The major findings of our study are summarized below. In their current form the normalized ATLAS 7 TeV data cannot be fit simultaneously with the 8 TeV $p_{T}^{Z}$ data. It also cannot be fit together with HERA data, nor in a global fit. The normalization performed on the 7 TeV data ties together the low and high $p_{T}^{Z}$ regions. When we perform the fit on the high$-p_{T}^{Z}$ region needed for a stable fixed-order QCD prediction, thus on a region in $p_{T}^{Z}$ which is different from the one used to normalise the data, the correlations between the bins are lost. The inclusion of this data requires either the experimental covariance matrix for the $p_{T}^{Z}>30$ GeV range only, the unnormalized data, or the inclusion of low-$p_{T}^{Z}$ resummation in the theoretical prediction. This last option would introduce an additional theoretical uncertainty into the fit.

The extreme precision of the 8 TeV $p_{T}^{Z}$ data binned in rapidity, with uncertainties at the few-per-mille level for the majority of bins, necessitates the introduction of an additional uncorrelated uncertainty for a fit with a low $\chi^{2}$ per degree of freedom. This additional parameter is meant to cover the residual theoretical uncertainty and the Monte-Carlo integration uncertainty on the theoretical prediction, as well as possible under-reported experimental errors. While the introduction of this extra uncertainty improves the $\chi^{2}$ per degree of freedom of the fit, we have varied the chosen value of this parameter to check that it has little impact on the actual PDFs obtained from the fit.

Including the 8 TeV $p_{T}^{Z}$ data into a global fit based on the NNPDF3.0 settings results in a significant reduction of the 13 TeV gluon-gluon, quark-gluon and quark-antiquark luminosity errors. To quantify this we have computed the gluon-fusion Higgs production using our NNPDF3.0 baseline, before and after including the $p_{T}^{Z}$ data in the fit. We find that the PDF uncertainty on the Higgs cross section decreases by 30%, while the central value of the prediction increases by 1%, within the previously-estimated uncertainty. We caution that this quantitative estimate of uncertainty reduction holds upon including only the $p_{T}^{Z}$ data into the NNPDF3.0 baseline fit. If additional data sets are included as well, these numbers will change. However, given the power of the $p_{T}^{Z}$ data found in our study, we expect that future global fits using this data will observe similar results.

Acknowledgements

We thank Jan Kretzschmar, Maarten Boonekamp and Stefano Camarda for insightful discussions on the ATLAS experimental data. We thank Gavin Salam for discussions on possible non-perturbative corrections to the $p_{T}^{Z}$ distribution, and Alexander Mueck for sending results for the NLO electroweak corrections. We thank Alex Mitov for useful discussions on the comparison with the inclusion of top data. We acknowledge extensive discussions and cross checks of results and PDF fits with the members of the NNPDF Collaboration, in particular Zahari Dim, Luca Rottoli, Stefano Carrazza, Nathan Hartland and Juan Rojo. We thank the Kavli Institute for Theoretical Physics in Santa Barbara for hosting the authors during the completion of this manuscript. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915 to the Kavli Institute of Theoretical Physics in Santa Barbara. R. B. is supported by the DOE contract DE-AC02-06CH11357. F. P. is supported by the DOE grants DE-FG02- 91ER40684 and DE-AC02-06CH11357. M. U. is supported by a Royal Society Dorothy Hodgkin Research Fellowship and partially supported by the STFC grant ST/L000385/1. A. G. is supported by the European Union s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 659128 - NEXTGENPDF. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.