The transverse momentum spectrum of weak gauge bosons at N$^3$LL+NNLO

TL;DR

This paper provides precise QCD predictions for the transverse momentum spectrum of electroweak gauge bosons at the LHC, combining NNLO calculations with N3LL resummation, reducing uncertainties and revealing stable ratios between related processes.

Contribution

It introduces a combined NNLO and N3LL approach for accurate pT spectrum predictions of gauge bosons at 13 TeV LHC, improving precision and understanding of process ratios.

Findings

Higher order corrections significantly alter the pT distribution shape.

Residual uncertainties are reduced to a few percent across the spectrum.

Ratios of charged- and neutral-current Drell-Yan pT distributions are highly stable.

Abstract

We present accurate QCD predictions for the transverse momentum pT spectrum of electroweak gauge bosons at the LHC for 13 TeV collisions, based on a consistent combination of a NNLO calculation at large pT and N3LL resummation in the small pT limit. The inclusion of higher order corrections leads to substantial changes in the shape of the differential distributions, and the residual perturbative uncertainties are reduced to the few percent level across the whole transverse momentum spectrum. We examine the ratio of pT distributions in charged- and neutral-current Drell-Yan production, and study different prescriptions for the estimate of perturbative uncertainties that rely on different degrees of correlation between these processes. We observe an excellent stability of the ratios with respect to the perturbative order, indicating a strong correlation between the corresponding QCD corrections.

Figures

Equations

$$\ln$$

$$=\sum _n\left\{\mathcal{O}\left({\alpha }_{\mathrm{s}}^n{L}^{n+1}\right)+\mathcal{O}\left({\alpha }_{\mathrm{s}}^n{L}^n\right)+\dots \right\}.$$

$${\sigma }_{X+\text{jet}}^{\text{NNLO}}=$$

$$+$$

$${\mathrm{\Sigma }}_{\mathrm{match}}^{{\mathrm{N}}^3\mathrm{LL}}(p_{\perp })=\frac{{\mathrm{\Sigma }}^{{\mathrm{N}}^3\mathrm{LL}}(p_{\perp })}{{\mathrm{\Sigma }}_{asym.}^{{\mathrm{N}}^3\mathrm{LL}}}{\left[{\mathrm{\Sigma }}_{asym.}^{{\mathrm{N}}^3\mathrm{LL}}\frac{{\mathrm{\Sigma }}^{{\mathrm{N}}^3\mathrm{LO}}(p_{\perp })}{{\mathrm{\Sigma }}_{exp.}^{{\mathrm{N}}^3\mathrm{LL}}(p_{\perp })}\right]}_{{\mathrm{N}}^3\mathrm{LO}},$$

$${\mathrm{\Sigma }}^{{\mathrm{N}}^3\mathrm{LO}}(p_{\perp })={\sigma }_{\mathrm{tot}}^{{\mathrm{N}}^3\mathrm{LO}}-{\int }_{p_{\perp }}^{\infty }\mathrm{d}{p}_{\perp }^{\prime }\frac{\mathrm{d}{\mathrm{\Sigma }}^{\mathrm{NNLO}}(p_{\perp }^{\prime })}{\mathrm{d}{p}_{\perp }^{\prime }},$$

$${\mathrm{\Sigma }}^{{\mathrm{N}}^3\mathrm{LL}}(p_{\perp })\underset{p_{\perp }\gg M}{\to }{\mathrm{\Sigma }}_{asym.}^{{\mathrm{N}}^3\mathrm{LL}}.$$

$${\mathrm{\Sigma }}_{\mathrm{match}}^{\mathrm{NNLL}}(p_{\perp })$$

$${\mathrm{\Sigma }}_{\mathrm{match}}^{\mathrm{NLL}}(p_{\perp })$$

$${\mathrm{\Sigma }}_{\mathrm{match}}^{{\mathrm{N}}^3\mathrm{LL}}(p_{\perp })$$

$${\mathrm{\Sigma }}_{\mathrm{match}}^{\mathrm{NNLL}}(p_{\perp })$$

$${\mathrm{\Sigma }}_{\mathrm{match}}^{\mathrm{NLL}}(p_{\perp })$$

$$M_Z=91.1876\mathrm{GeV},M_W=80.379\mathrm{GeV},$$

$${\mathrm{\Gamma }}_Z=2.4952\mathrm{GeV},{\mathrm{\Gamma }}_W=2.085\mathrm{GeV},$$

$$G_F=1.1663787\times 10^{-5}{\mathrm{GeV}}^{-2}.$$

$$|{\vec{p}}_{\perp }^{{\ell }^{\pm }}|>25\mathrm{GeV},|{\eta }^{{\ell }^{\pm }}|<2.5,$$

$$66\mathrm{GeV}<M_{\ell \ell }<116\mathrm{GeV},$$

$$|{\vec{p}}_{\perp }^{\ell }|>25\mathrm{GeV},|/{\vec{E}}_T|>25\mathrm{GeV},$$

$$|{\eta }^{\ell }|<2.5,m_T>50\mathrm{GeV},$$

$$m_T=\sqrt{{\left(|{\vec{p}}_{\perp }^{\ell }|+|/{\vec{E}}_T|\right)}^2-{\left(|{\vec{p}}_{\perp }^{\ell }+/{\vec{E}}_T|\right)}^2}.$$

$${\Delta }^{{\mathrm{N}}^3\mathrm{LL}}\equiv (\mathrm{d}{\Sigma }^{{\mathrm{N}}^3\mathrm{LL}+\mathrm{NNLO}}/\mathrm{d}{p}_{\perp }-\mathrm{d}{\Sigma }^{{\mathrm{N}}^3\mathrm{LL}}/\mathrm{d}{p}_{\perp })/{\sigma }_{\mathrm{tot}}^{\mathrm{NNLO}}$$

$$\frac{1}{2}\le \frac{x_{{\mu }_F}^{num.}}{x_{{\mu }_F}^{den.}}\le 2,$$

$$\frac{1}{2}\le \frac{x_{\mu }^{num.}}{x_{\mu }^{den.}}\le 2,$$

Full text

1 Introduction

The differential spectrum of electroweak gauge bosons, measured via their leptonic decays, is among the most prominent observables at the LHC.

Owing to the outstanding precision of their experimental measurement Chatrchyan:2011wt ; Aad:2012wfa ; Aad:2014xaa ; Aad:2015auj ; Aaij:2015gna ; Khachatryan:2015oaa ; Khachatryan:2015paa ; Aaij:2015zlq ; Aad:2016izn ; Khachatryan:2016nbe ; Aaij:2016mgv ; Sirunyan:2017igm ; Aaboud:2017svj ; Aaboud:2017ffb , such observables allow for a precise extraction of some of the Standard Model (SM) parameters — such as the $W$ boson mass Aaboud:2017svj , or parton densities Ball:2017nwa ; Boughezal:2017nla ; Bacchetta:2017gcc ; Bertone:2019nxa — as well as for the calibration of widely used event generators and analysis tools. For this reason, an accurate theoretical understanding of such observables is paramount to exploit the precise data and perform meticulous tests of the SM.

Inclusive and differential distributions for neutral and charged Drell-Yan (DY) production with lepton pair invariant mass $M$ are nowadays known with very high precision. The total cross section is known fully differentially in the Born phase space up to NNLO Hamberg:1990np ; vanNeerven:1991gh ; Anastasiou:2003yy ; Melnikov:2006di ; Melnikov:2006kv ; Catani:2010en ; Catani:2009sm ; Gavin:2010az ; Anastasiou:2003ds , while differential distributions in transverse momentum $p_{\perp}$ were recently computed up to NNLO both for $Z$- Ridder:2015dxa ; Ridder:2016nkl ; Gehrmann-DeRidder:2016jns ; Gauld:2017tww ; Boughezal:2015ded ; Boughezal:2016isb and $W$-boson Boughezal:2015dva ; Boughezal:2016dtm ; Gehrmann-DeRidder:2017mvr production. In the DY distributions, electroweak corrections become important especially at large transverse momenta, and they have been computed to NLO accuracy in Kuhn:2005az ; Kuhn:2007qc ; Denner:2009gj ; Denner:2011vu .

In kinematical regimes dominated by soft and collinear radiation, the fixed-order perturbative series for the differential $p_{\perp}$ distribution is affected by large logarithmic terms of the form $\alpha_{\mathrm{s}}^{n}L^{2n-1}/p_{\perp}$, with $L\equiv\ln(M/p_{\perp})$, which must be resummed to all orders for a reliable theoretical prediction. In such regimes, the perturbative (logarithmic) accuracy is defined in terms of the logarithm of the cumulative cross section ${\rm\Sigma}$ as

[TABLE]

One refers to the dominant terms $\alpha_{\mathrm{s}}^{n}L^{n+1}$ as leading logarithmic (LL), to terms $\alpha_{\mathrm{s}}^{n}L^{n}$ as next-to-leading logarithmic (NLL), to $\alpha_{\mathrm{s}}^{n}L^{n-1}$ as next-to-next-to-leading logarithmic (NNLL), and so on. The resummation of the $p_{\perp}$ spectrum of SM bosons has been studied in a multitude of theoretical formulations throughout the years Parisi:1979se ; Collins:1984kg ; Balazs:1997xd ; Landry:2002ix ; Becher:2010tm ; Bozzi:2010xn ; Becher:2011xn ; GarciaEchevarria:2011rb ; Monni:2016ktx ; Ebert:2016gcn ; Bizon:2017rah , and the current state of the art for phenomenological studies at the LHC reaches N3LL accuracy Bizon:2017rah ; Chen:2018pzu ; Bizon:2018foh ; Becher:2019bnm .

In this article, we reach a new milestone in the theoretical description of transverse momentum distributions in both neutral and charged DY production, aiming for per-cent level precision throughout the full kinematical range. This is achieved by matching the fixed-order NNLO QCD predictions with the N3LL resummation of large logarithmic corrections. We adopt the momentum-space formulation of refs. Monni:2016ktx ; Bizon:2017rah , in which the resummation is performed by generating the QCD radiation by means of a Monte Carlo (MC) algorithm. All the necessary ingredients for the N3LL $p_{\perp}$ resummation have been computed in refs. Catani:2012qa ; Gehrmann:2014yya ; Echevarria:2016scs ; Li:2016ctv ; Vladimirov:2016dll ; Moch:2018wjh ; Lee:2019zop . The combined framework enables fully differential N3LL+NNLO predictions for distributions that take proper account of the fiducial volume definitions used in the experimental measurements.

The article is organised as follows. In section 2 we briefly review the computation of the NNLO differential distributions in DY-pair production with the parton-level code NNLOjet, as well as the resummation for the $p_{\perp}$ distributions using the computer program RadISH. Section 3 describes our results for $13~{}\mathrm{TeV}$ LHC collisions. Finally, section 4 contains our conclusions.

2 Setup of the calculation

In this section we give a brief overview of the computational setup, and describe the ingredients of both the fixed order (Section 2.1) and the resummed (Section 2.2) calculations.

2.1 Fixed order

For the calculation of the DY process, we consider the off-shell production of either a pair of charged leptons (mediated by both a $Z$ boson and a virtual photon) or a charged lepton and a neutrino (mediated by $W^{\pm}$ bosons), in association with partonic jets. The jet requirement is replaced by a lower cut on the transverse momentum of the pair, that acts as an infrared regulator of the fixed-order calculation, hence preventing the radiation from being entirely unresolved.

The NNLO QCD predictions for neutral and charged DY production have been obtained in refs. Ridder:2015dxa ; Ridder:2016nkl ; Gehrmann-DeRidder:2016jns ; Gauld:2017tww ; Boughezal:2015ded ; Boughezal:2016isb ; Boughezal:2015dva ; Boughezal:2016dtm ; Gehrmann-DeRidder:2017mvr . Relative to the LO distribution, in which the leptonic system recoils against a single parton, the NNLO calculation receives contributions from configurations with two extra partons (RR: double-real corrections Campbell:2002tg ; Campbell:2003hd ; Hagiwara:1988pp ; Berends:1988yn ; Falck:1989uz ), with one extra parton and one extra loop (RV: real-virtual corrections Campbell:2002tg ; Campbell:2003hd ; Glover:1996eh ; Bern:1996ka ; Campbell:1997tv ; Bern:1997sc ) and with no extra partons but two extra loops (VV: double-virtual corrections Moch:2002hm ; Garland:2001tf ; Garland:2002ak ; Gehrmann:2011ab ). Each of the three contributions is separately infrared divergent either in an implicit manner from phase-space regions where the partonic radiation becomes unresolved (soft and/or collinear), or in a explicit manner from infrared poles in virtual loop corrections. Only the sum of the three contributions is finite.

We perform the calculation using the parton-level generator NNLOjet, which implements the antenna subtraction method GehrmannDeRidder:2005cm ; Daleo:2006xa ; Currie:2013vh to isolate infrared singularities and to enable their cancellation between different contributions prior to the numerical phase-space integration. The NNLO calculation can be structured as

[TABLE]

The antenna subtraction terms, $\mathrm{d}\sigma^{S,T,U}_{\text{NNLO}}$, are constructed from antenna functions GehrmannDeRidder:2005cm ; GehrmannDeRidder:2005aw ; Daleo:2009yj ; Boughezal:2010mc ; Gehrmann:2011wi ; GehrmannDeRidder:2012ja to cancel infrared singularities between the contributions of different parton multiplicities. The integrals are performed over the phase space ${\Phi_{X+1,2,3}}$ corresponding to the production of the colour singlet in association with one, two or three partons in the final state. The integration over the final-state phase space is fully differential such that any infrared-safe observable $\cal O$ can be studied through differential distributions as $\mathrm{d}\sigma^{\text{NNLO}}_{X+\text{jet}}/\mathrm{d}\cal O$.

The matching of the above NNLO prediction to a resummed calculation in the small $p_{\perp}$ limit is computationally very challenging. At small $p_{\perp}$, both the matrix elements and the subtraction terms grow rapidly in magnitude due to the presence of un-cancelled infrared singularities. This results in large numerical cancellations between them that ultimately challenge the stability of the final prediction. The finite remainder of such cancellations needs to be numerically stable in order to be consistently combined with a resummed calculation and extrapolated to the limit $p_{\perp}\to 0$. The stability of NNLOjet in this extreme regime has been tested thoroughly against the expansion of the N3LL resummations in refs. Chen:2018pzu ; Bizon:2018foh , where it is shown that the NNLO calculation can be reliably obtained down to very small $p_{\perp}$ values.

The residual infrared (logarithmic) divergences that persist in the $p_{\perp}\to 0$ limit are cancelled by combining the fixed-order computation with a resummed calculation, where the logarithms in the fixed-order prediction are subtracted and replaced by the sum of the corresponding enhanced terms to all orders in perturbation theory. This procedure is discussed in the following Section 2.2.

2.2 Resummation and matching

The resummation is performed in momentum space by means of the method formulated in refs. Monni:2016ktx ; Bizon:2017rah and implemented in the computer code RadISH. In this approach, the factorisation properties of the QCD matrix elements in the soft and collinear limits are exploited to devise a numerical procedure to generate the radiation at all perturbative orders. This allows for the resummation of the large logarithmic terms in a fashion that is similar in spirit to a Monte Carlo generator. A detailed technical description of the method can be found in Refs. Monni:2016ktx ; Bizon:2017rah , and the formulae up to N3LL accuracy are collected in Ref. Bizon:2018foh (Section 3 and Appendix B).

In order to have a reliable prediction across the whole $p_{\perp}$ spectrum, the fixed-order and resummed predictions must be consistently combined through a matching procedure. The matching is performed in such a way that it reduces to the resummed calculation at small $p_{\perp}$, while reproducing the fixed-order prediction at large transverse momentum. At a given perturbative order, one can adopt various schemes that differ from one another by terms beyond the considered order. In the present analysis we adopt the multiplicative scheme formulated in refs. Caola:2018zye ; Bizon:2018foh , in which the matching is performed at the level of the cumulative distribution (1) as follows:

[TABLE]

where ${\rm\Sigma}_{\rm exp.}^{\rm N^{3}LL}$ denotes the expansion of the resummation formula ${\rm\Sigma}^{\rm N^{3}LL}$ to ${\cal O}(\alpha_{s}^{3})$ (N3LO), and the whole squared bracket in Eq. (3) is expanded to N3LO. It should be recalled that ${\cal O}(\alpha_{s}^{3})$ corresponds to N3LO in the total (i.e. $p_{\perp}$-integrated) cross section and in any cumulative distribution, while being NNLO in the fixed-order $p_{\perp}$-differential cross section.

In the above equation, ${\rm\Sigma}^{\rm N^{3}LO}$ is the cumulative fixed-order distribution at N3LO, i.e.

[TABLE]

where $\sigma_{\rm tot}^{\rm N^{3}LO}$ is the total cross section for the charged or neutral DY processes at N3LO, and $\mathrm{d}{\rm\Sigma}^{\rm NNLO}/\mathrm{d}p_{\perp}^{\prime}$ denotes the corresponding NNLO $p_{\perp}$-differential distribution obtained with NNLOjet. Both of these quantities are accurate to ${\cal O}(\alpha_{s}^{3})$. Since the N3LO inclusive cross section for DY production is currently unknown, we approximate it with the NNLO cross section Hamberg:1990np ; vanNeerven:1991gh ; Anastasiou:2003yy ; Melnikov:2006di ; Melnikov:2006kv ; Catani:2010en ; Catani:2009sm ; Gavin:2010az ; Anastasiou:2003ds in the following. This approximation impacts only terms at N4LL order, and is thus beyond the accuracy considered in this study.

Finally, the quantity ${\rm\Sigma}^{\rm N^{3}LL}_{\rm asym.}$ is defined as the asymptotic ($p_{\perp}\gg M$) limit of the resummed cross section

[TABLE]

This prescription ensures that, in the $p_{\perp}\gg M$ limit, Eq. (3) reproduces by construction the fixed-order result, while in the $p_{\perp}\to 0$ limit it reduces to the resummed prediction. The detailed matching formulae can be found in Appendix A of ref. Bizon:2018foh .

In the next section, we will also report matched predictions at lower perturbative orders, NNLL+NLO and NLL+LO, that are obtained from the following matched cumulative distributions

[TABLE]

The above matching schemes guarantee that in the large-$p_{\perp}$ limit the matched cumulative cross sections reproduce, by construction, the following total cross sections

[TABLE]

We stress once more that the ${\rm\Sigma}_{\rm match}^{\rm N^{3}LL}$ reproduces the NNLO total cross section at large $p_{\perp}$ since the N3LO result for the inclusive DY process is currently unknown. The nominal accuracy of the predictions is unaffected by this choice.

The final normalised distributions that will be reported in Section 3 are obtained by differentiating Eqs. (3), (6) and (7), and dividing by the respective total cross sections of the right hand side of Eq. (8).

We recall that the resummed calculation contains a Landau singularity arising from configurations where the radiation is generated with transverse momentum scales $k_{\perp}\sim M\,\exp\left\{-1/(2\beta_{0}\alpha_{s})\right\}$ (with $\alpha_{s}=\alpha_{s}(M)$ and $\beta_{0}=(11\,C_{A}-2\,n_{f})/(12\pi)$). In the predictions presented in the following, we set the results to zero when the hardest radiation’s transverse momentum reaches the singularity. For the leptonic invariant masses studied here, this procedure acts on radiation emitted at very small transverse momentum that, due to the vectorial nature of the observable Parisi:1979se ; Bizon:2017rah , gives a very small contribution to the spectrum. We however stress that for a precise description of this kinematic regime, a thorough study of the impact of non-perturbative corrections is necessary.

3 Results at the LHC

In this section we report our numerical results for the neutral and charged DY transverse momentum distributions at N3LL+NNLO.

We consider $pp$ collisions at a centre-of-mass energy of $13~{}\mathrm{TeV}$, and we use the NNLO NNPDF3.1 parton distribution function set Ball:2017nwa with $\alpha_{\mathrm{s}}(M_{Z})=0.118$. The parton densities are evolved from a low scale $Q_{0}\sim 1~{}\mathrm{GeV}$ forwards with LHAPDF Buckley:2014ana , which correctly implements the heavy quark thresholds in the PDFs. All convolutions are handled with the Hoppet package Salam:2008qg . In the results reported below, we use the NNLO DGLAP evolution of the adopted PDF set for all perturbative orders shown in the figures. Although the NNLO corrections to the PDF evolution are formally of order N3LL, we include them also in the NLL and NNLL predictions in order to guarantee a consistent treatment of the quark thresholds in the parton densities. We note that this choice will lead to numerical differences in comparison to other NLL and NNLL results shown in the literature.

We adopt the $G_{\mu}$ scheme with the electro-weak parameters taken from the PDG Tanabashi:2018oca , that is

[TABLE]

Moreover, we set the CKM matrix equal to the identity matrix, and we have verified that this approximation is accurate at the few-permille level. For both neutral-current and charged-current DY we apply fiducial selection cuts that resemble the ones used by ATLAS in previous analyses Aad:2015auj .

The final state for the neutral DY production is defined by applying the following set of fiducial selection cuts on the leptonic pair:

[TABLE]

where $\vec{p}_{\perp}^{\ell^{\pm}}$ are the transverse momenta of the two leptons, $\eta^{\ell^{\pm}}$ are their pseudo-rapidities in the hadronic centre-of-mass frame, and $M_{\ell\ell}$ is the invariant mass of the di-lepton system. The central factorisation and renormalisation scales are chosen to be $\mu_{R}=\mu_{F}=\sqrt{M_{\ell\ell}^{2}+|\vec{p}_{\perp}^{Z}|^{2}}$ and the central resummation scale is set to $Q=M_{\ell\ell}/2$.

In the case of charged DY production, the fiducial volume is defined as

[TABLE]

where $\not{\vec{E}}_{T}$ is the missing transverse energy vector and

[TABLE]

The central factorisation and renormalisation scales are chosen to be $\mu_{R}=\mu_{F}=\sqrt{M_{\ell\nu}^{2}+|\vec{p}_{\perp}^{W}|^{2}}$ and the central resummation scale is again set to $Q=M_{\ell\nu}/2$.

In both processes, we assess the missing higher-order uncertainties by performing a variation of the renormalisation and factorisation scales by a factor of two around their respective central values whilst keeping $1/2\leq\mu_{R}/\mu_{F}\leq 2$. In addition, for central factorisation and renormalisation scales, we vary the resummation scale $Q$ by a factor of two in either direction. The final uncertainty is built as the envelope of the resulting nine-scale variation.

3.1 Fiducial distributions

We start by showing, in Figure 1, the comparison of the $Z$ and $W^{\pm}$ normalised distributions at NLL+LO (green), NNLL+NLO (blue), and N3LL+NNLO (red) in the fiducial volumes defined above. The lower inset of each panel of Figure 1 shows the ratio of all predictions to the previous state of the art (NNLL+NLO), with the same colour code as in the main panels. The difference between each prediction and the next order is of ${\cal O}(\alpha_{s})$, both in the large $p_{\perp}$ region and in the limit $p_{\perp}\to 0$ where $\alpha_{s}L\sim 1$.

In comparison to the NNLL+NLO result, we note that the N3LL+NNLO corrections lead to important distortions in the shape of the distributions, making the spectrum harder for $p_{\perp}\gtrsim 10~{}\mathrm{GeV}$, and softer below this scale. The perturbative errors are reduced by more than a factor of two across the whole $p_{\perp}$ range, and the leftover uncertainty is at the $5\%$ level. In general, we observe a good convergence of the perturbative description when the order is increased, although in some $p_{\perp}$ regions the N3LL+NNLO and the NNLL+NLO bands overlap only marginally. This feature can be understood by noticing that, as mentioned in Section 2.2, both predictions are normalised to the same NNLO total cross section. Since at large $p_{\perp}$ the NNLO corrections lead to an increase in the spectrum of about $10\%$, by unitarity of the matching procedure (that preserves the total cross section) this must be balanced by an analogous decrease in the spectrum in the region governed by resummation, as indeed observed in Figure 1. We stress, nevertheless, that the two orders are compatible within the quoted uncertainties.

In Figure 2, we show the comparison among the NNLO (green), the NNLL+NLO (blue), and N3LL+NNLO (red) predictions, where the bands are obtained as discussed above. Alongside these results, we also show the Monte Carlo predictions obtained using the Pythia8 generator Sjostrand:2014zea with the AZ tune Aad:2014xaa , that has been obtained from the $Z$-boson $p_{\perp}$ distribution at $7~{}\mathrm{TeV}$. At $7~{}\mathrm{TeV}$ and $8~{}\mathrm{TeV}$ the above tune is known to correctly describe the $Z$ transverse momentum spectrum within a few percent for $p_{\perp}\lesssim 50~{}\mathrm{GeV}$ Aad:2014xaa , and it has been employed in the extraction of the $W$-boson mass by the ATLAS collaboration Aaboud:2017svj . Although it is currently unknown how this tune performs at $13~{}\mathrm{TeV}$ in comparison to the data, we use the Pythia8 prediction for reference in the following plots. In particular, the lower inset of each panel of Figure 2 shows the ratio of all predictions to Pythia8. We observe a reasonable agreement between the N3LL+NNLO predictions and Pythia8 below $30~{}\mathrm{GeV}$, while it deteriorates for larger $p_{\perp}$ values. This feature is particularly visible in the case of $W^{\pm}$ production.

A comparison of the N3LL+NNLO band to the fixed-order one shows that the resummation starts making a significant difference for $p_{\perp}\lesssim 20~{}\mathrm{GeV}$, while above this scale the NNLO provides a reliable theoretical prediction. To further quantify the relative impact of the non-singular contributions in this region, we show in Fig. 3 the difference

[TABLE]

between the matched and the resummed predictions for the $Z$ and $W^{\pm}$ normalised distributions. In the lower panel of the plot we show the relative size of $\Delta^{\rm N^{3}LL}$ with respect to the matched N3LL+NNLO result, $\Delta^{\rm N^{3}LL}/(\Sigma^{\rm N^{3}LL+NNLO}/\mathrm{d}p_{\perp}/\sigma_{\rm tot}^{\rm NNLO})$. The non-singular contributions are somewhat larger for $W^{\pm}$; the relative size of $\Delta^{\rm N^{3}LL}$ with respect to the N3LL+NNLO result is smaller than $5\%$ ($10\%$) for $Z$ ($W^{\pm}$) for $p_{\perp}\lesssim 10$ GeV, and becomes larger than $10\%$ ($20\%$) for $p_{\perp}>20$ GeV.

3.2 Ratio of $Z/W$ and $W^{-}/W^{+}$ distributions

Another set of important quantities of interest are the ratios of the above distributions, which play a central role in recent extractions of the $W$-boson mass at the LHC Aaboud:2017svj . When taking ratios of perturbative quantities one has to decide how to combine the uncertainties in the numerator and denominator to obtain the final error.

One option is to try to identify the possible sources of correlation in the three processes considered here. From the point of view of the perturbative (massless) QCD description adopted in this study, one expects the structure of radiative corrections to such reactions to be nearly identical. This is certainly the case as far as resummation is concerned, since it is governed by the same anomalous dimensions and all-order structure in $W$ and $Z$ production. As a consequence, the resummation scale should be varied in a correlated manner in both predictions considered in the ratio. A similar argument can be made regarding the renormalisation scale $\mu_{R}$ and the factorisation scale $\mu_{F}$.

However, an important difference between $Z$, $W^{+}$, and $W^{-}$ production lies in the different combination of partonic channels probed by each process and, in particular, in the sensitivity to different heavy quark thresholds in the PDFs at small $p_{\perp}$. Therefore, it is not clear whether a fully correlated variation of the factorisation scale $\mu_{F}$ is physically justified. A more conservative uncertainty prescription is to vary the scales $\mu_{R}$ and $Q$ in numerator and denominator in a fully correlated way, while varying $\mu_{F}$ in an uncorrelated manner within the constraint Gehrmann-DeRidder:2017mvr

[TABLE]

where $x_{\mu_{F}}$ is the ratio of the factorisation scale to its central value. This corresponds to a total of $17$ scale combinations.

Finally, for comparison we also consider the uncorrelated variation of $\mu_{R}$ and $\mu_{F}$ in the ratio, while imposing

[TABLE]

where $x_{\mu}$ is the ratio of the scale $\mu$ to its central value, with $\mu\equiv\{\mu_{R},\mu_{F}\}$, together with a correlated variation of the resummation scale $Q$. This recipe amounts to taking the envelope of the predictions resulting from $33$ different combinations of scales in the ratio.

To examine the reliability of the above uncertainty schemes, in Figure 4 we analyse the convergence of the perturbative series for the ratios of distributions, by comparing the results at NLL+LO (green), NNLL+NLO (blue), and N3LL+NNLO (red). The three lower panels in each plot show the theory uncertainties obtained according to the three prescriptions outlined above, respectively, in comparison to the old state-of-the-art prediction at NNLL+NLO. In the case of the $Z/W^{+}$ ratio (shown in the upper plot of Figure 4), we observe that the different perturbative orders are very close to one another. The results are compatible even within the uncertainty bands obtained with the more aggressive error estimate, which in some bins is sensitive to minor statistical fluctuations due to the complexity of the NNLO calculation. This feature is strikingly evident in the case of the $W^{-}/W^{+}$ ratio (lower figure), where the excellent convergence of the series seems to indicate that either a fully correlated scale variation or the more conservative estimate of Eq. (14) is perfectly justified.

Figure 5 shows the comparison of the same two ratios ($Z/W^{+}$ and $W^{-}/W^{+}$) to the NNLO result (green), and to Pythia8. We observe that in both cases the N3LL+NNLO calculation leads to an important reduction of the theory uncertainty. In particular, even with the most conservative estimate of the theory error, our best prediction leads to errors of the order of $5\%$, with the exception of the first bin where the perturbative uncertainty is at the $10\%$ level. The kink around $p_{\perp}\sim 50-60~{}\mathrm{GeV}$ in the $Z/W^{+}$ ratio (upper plot in Figure 5) is due to the different fiducial selection cuts in the two processes. A change in the shape of the distributions around this scale is indeed visible in Figure 2, at slightly different $p_{\perp}$ values for $Z$ and $W^{+}$ production, respectively, that is reflected in the structure observed in Figure 5. We find a good agreement between our best predictions at N3LL+NNLO and the Pythia8 Monte Carlo in the small $p_{\perp}$ region of the ratios. However, the two predictions are not compatible within the quoted theory uncertainties if the scales are varied in a fully correlated manner. On the other hand, for $p_{\perp}\gtrsim 40~{}\mathrm{GeV}$, the Pythia8 result disagrees with the matched calculation. This behavior is not unexpected, since the nominal perturbative accuracy of Pythia8 is well below any of the matched calculations, and the AZ tune is optimised to describe the $Z$ spectrum in the region $p_{\perp}\leq 50~{}\mathrm{GeV}$ at $7~{}\mathrm{TeV}$.

4 Conclusions

In this work, we computed the transverse momentum distributions of electroweak gauge bosons at the LHC to N3LL+NNLO accuracy in QCD. This calculation opens up a new level of theoretical precision in the description of these observables. The new state-of-the-art prediction is obtained by combining the NNLO results from the NNLOjet program with the N3LL resummation performed with RadISH. Our phenomenological study adopts fiducial selection cuts similar to the setup adopted by ATLAS in previous studies. The numerical results we presented are made available in electronic format as additional material alongside this manuscript.

We find that, in comparison to the fixed-order prediction, the resummation effects become important for $p_{\perp}\lesssim 20~{}\mathrm{GeV}$. The effect of the N3LL+NNLO corrections with respect to the previous NNLL+NLO prediction is as large as $\sim 10\%$, and leads to significant shape distortions as well as to a substantial reduction in the perturbative uncertainty due to missing higher-order corrections. In particular, the distributions considered in this article are predicted with a residual uncertainty below the $5\%$ level across most of the $p_{\perp}$ spectrum. We also compared the results to the predictions obtained from the Pythia8 Monte Carlo with the AZ tune, that has been determined using the ATLAS experimental data for the $Z$ boson transverse momentum at $7~{}\mathrm{TeV}$ Aad:2014xaa .

Finally, we examined the ratios of the $Z$ to $W^{+}$, and $W^{-}$ to $W^{+}$ distributions, which play an important role in the $W$ mass extraction at the LHC. We consider different prescriptions for the estimate of perturbative uncertainties that rely on different degrees of correlation between the scales in the numerator and in the denominator. We find a remarkable convergence of the predictions for the ratios at different perturbative orders. This fact strongly indicates that the class of processes considered in this study feature very similar perturbative corrections suggesting that the perturbative sources of uncertainty are correlated to a large extent.

There are, however, additional sources of perturbative corrections to $W^{\pm}$ and $Z$ production that we ignored in our study. In particular, at the level of the residual theoretical errors obtained in our predictions, PDF theory uncertainties Harland-Lang:2018bxd ; AbdulKhalek:2019bux , QED corrections deFlorian:2018wcj ; Cieri:2018sfk , as well as a careful study of the impact of mass effects Aivazis:1993pi ; Thorne:1997ga ; Collins:1998rz ; Berge:2005rv ; Forte:2010ta ; Pietrulewicz:2017gxc ; CarloniCalame:2016ouw ; Krauss:2017wmx ; Figueroa:2018chn ; Forte:2018ovl ; Bagnaschi:2018dnh become necessary. The correlation pattern between the uncertainties due to such effects may well be different from what we have observed in this paper, and a dedicated study must be performed in order to reliably combine these effects with the N3LL+NNLO predictions presented here.

Acknowledgements

We are very grateful to Jan Kretzschmar for clarifying discussions and for kindly providing us with the Pythia8 setup used by ATLAS and with template predictions for the $W$ and $Z$ differential distributions. We would like to thank Paolo Torrielli for fruitful discussions on the topics presented here, and Stefano Camarda for useful correspondence on the impact of heavy flavours. We also thank Fabrizio Caola, Lucian Harland-Lang, Jonas Lindert, Gavin Salam, and Marek Schönherr for constructive comments on the uncertainty prescriptions. The work of PM is supported by the Marie Skłodowska Curie Individual Fellowship contract number 702610 Resummation4PS. LR is supported by the ERC Starting Grant REINVENT (714788) and acknowledges the CERN Theoretical Physics Department for hospitality and support during part of this work. This work has been supported in part by the Swiss National Science Foundation (SNF) under grant numbers 200021-172478 and 200020-175595, and by the Research Executive Agency (REA) of the European Union through the ERC Consolidator Grant HICCUP (614577) and the ERC Advanced Grant MC@NNLO (340983).