The $q_T$ and $\Delta\phi$ spectra in W and Z production at the LHC at N$^3$LL'+N$^2$LO

The production of weak gauge bosons, $W^{\pm}$ and $Z$, are at the core of the LHC precision measurement program. Their transverse momentum spectra as well as their pairwise ratios are key theoretical inputs to many high-precision analyses, ranging from the $W$ mass measurement to the determination of parton distribution functions. Owing to the different properties of the $W$ and $Z$ boson and the different accessible fiducial regions for their measurement, a simple one-dimensional correlation is insufficient to capture the differing vector and axial-vector dynamics of the produced lepton pair. We propose to correlate them in two observables, the transverse momentum $q_T$ of the lepton pair and its azimuthal separation $\Delta\phi$. Both quantities are purely transverse and therefore accessible in all three processes, either directly or by utilising the missing transverse momentum of the event. We calculate all the single-differential $q_T$ and $\Delta\phi$ as well as the double-differential $(q_T,\Delta\phi)$ spectra for all three processes at N$^3$LL'+N$^2$LO accuracy, resumming small transverse momentum logarithms in the soft-collinear effective theory approach and including all singlet and non-singlet contributions. Using the double-differential cross sections we build the pairwise ratios $\mathcal{R}_{W^+/Z},~\mathcal{R}_{W^-/Z},~\mathcal{R}_{W^+/W^-}$ and determine their uncertainties assuming fully correlated, partially correlated, and uncorrelated uncertainties in the respective numerators and denominators. In the preferred partially correlated case we find uncertainties of less than 1% in most phase space regions and up to 3% in the lowest $q_T$ region.


Introduction
One of the most important observables at the LHC is the differential spectrum of the electroweak gauge bosons in their leptonic decay channels [1][2][3][4][5][6][7][8][9][10][11]. The extraordinary precision reached by the ATLAS and CMS collaborations in their measurements enables the precision extraction of the parameters of the Standard Model (SM), such as the W boson mass [12] and parton densities [13][14][15][16]. In order to exploit this precision data to its fullest, however, theory calculations of equal precision are indispensable. Of particular interest here are angular observables of the final state leptons, such as φ * in Z production [17,18], as the angular resolution of charged objects is much more precise than their energy resolution, which is needed for p Ttype observables. The purely transverse azimuthal decorrelation ∆φ of the lepton pair carries the same experimental advantages as φ * , while being less favoured theoretically as it is not weighted by the scattering angle and, thus, less sensitive to the vector boson transverse momentum q T . Unfortunately, the measurement of W production always involves the determination of the event's missing momentum as a proxy for the inaccessible neutrino momentum. Furthermore, only the transverse part of the missing momentum can be determined at hadron colliders due to the composite nature of the incident protons and the incomplete detector geometry, and thus all W observables have to be constructed solely in the transverse plane. As such, observables such as φ * can not be used, in contrast to its simpler version ∆φ, the azimuthal decorrelation of the lepton and the missing transverse momentum. For its lepton ingredient it has the same experimental advantages as φ * , depending only on the lepton direction but not its momentum. Conversely, however, the missing transverse momentum's transverse direction resolution, being determined by the sum of all other measurable particles' momentum vectors, is not significantly improved as compared to its magnitude [19][20][21][22]. Still, in combination, a better resolution for ∆φ should be achievable as compared to the transverse momentum q T of the reconstructed W boson. Similarly, due to the purely transverse nature of measurable missing momentum vector as well as the rapidity limitations of the physical lepton detectors (electromagnetic calorimeter and muon chambers), the fiducial regions for lepton-neutrino and lepton-pair final states differ by definition. Thus, any correlation of the two production processes that aims for the precise extrapolation from one to the other, as is paramount in the W mass measurements for example, should take into account detector-acceptance-induced difference of the internal dynamics of both systems. Hence, also multidifferential precision predictions of the W + , W − and Z production cross sections and their ratios are needed. On the theory side, the Drell-Yan processes have drawn extensive attention for decades. The total cross section at LO [23] was one of the earliest processes calculated in the Standard Model, and its NLO QCD corrections have been derived soon after [24][25][26][27][28]. The NNLO QCD corrections [29][30][31][32][33][34][35][36] and, very recently, the third-order QCD results [37,38] are known as well. Similarly, the transverse momentum spectrum of the gauge boson at finite q T is known also up to O(α 3 s ) [39][40][41][42][43][44][45][46][47][48][49][50]. In addition to the QCD corrections higher-order electro-weak (EW) effects [51][52][53][54][55][56][57][58][59][60][61][62] and QCD-EW mixed corrections emerging at two-loop order [63][64][65][66][67][68][69][70][71][72][73] are of importance at this level of accuracy as well. In addition, the small transverse momentum region of the gauge boson is of particular interest as it contains the bulk of the production cross section. In light of its sensitivity to the soft and collinear radiations, fixedorder calculations are dominated by the powers of large logarithms of the form ln λ T , with λ T ≡ q T /m and m being the mass of the produced (off-shell) gauge boson, spoiling the convergence of the perturbative expansion. It is thus imperative to resum these logarithms to all perturbative orders. Based on the infrared-collinear (IRC) properties of QCD, the exploration of the small-q T exponentiation has been a topic of investigation since the formulation of the theory [74][75][76][77][78][79][80][81][82][83], and the first all-order proof was achieved by Collins, Soper and Sterman [84]. After that, a formalism of recombining the occurring ingredients has been proposed in Refs. [85][86][87], such as to arrive at a process-independent Sudakov form factor. In the recent decades, many efforts have been devoted to this theme and alternative schemes have been proposed and implemented, such as the distributional space resummation [88], the direct momentumspace resummation [89][90][91][92] and a number of variants within the soft-collinear effective theory (SCET) [93][94][95]. As one of the more popular factorisation techniques, SCET enables a formal but flexible way to explore the factorisation properties in the small-q T domain [96][97][98][99][100][101][102]. With the progresses made in the framework development, the resummation accuracy has increased steadily. Throughout this work, we take the following counting rule for the resummation results, where the exponent indicates the anomalous dimension level required by the resummation accuracy, the desired perturbative level for the fixed-order functions is presented in the curly brackets. We have taken α s ln λ T ∼ O(1) here. In the previous investigations, the N 2 LL calculations were carried out in Refs. [96,[103][104][105][106][107][108][109], whilst the authors in Refs. [91,92,[110][111][112] have accomplished N 3 LL very recently. Thanks to the developments in fixed-order perturbative calculations, the cusp anomalous dimension has achieved four-loop accuracy [113] and the fixed-order ingredients (including the non-singlet quark form factor [114], soft [102,115] and beam functions [116][117][118]) at N 3 LO are available now. The non-cusp anomalous dimensions can be extracted from the latter sectors. Consequently, we present in this work N 3 LL accurate calculations where the fixed-order functions are improved by one power of α s order with respect to the unprimed accuracy. 1 During our calculation, not only are the ingredients mentioned above assembled, the singlet contributions in the neutral Drell-Yan process are also addressed for completeness. In spite of its expected smallness, the singlet contribution in fact acts as the essential constituent starting from N 2 LL or N 3 LL. To include it, the low energy effective field theory (LEEFT) [121][122][123][124][125][126] resulting from integrating out the top quark has been employed in this work. With the help of the availability of the complete four-momenta of the final state leptons and neutrinos, we compute the q T spectrum in an experimentally accessible fiducial region. Additionally, this work will also project the resummation of small-q T logarithms on the azimuthal decorrelation ∆φ, and compute single-and doubledifferential cross sections and their ratios. The paper is organised as follows: In Sec. 2 we detail the ingredients of our calculation, emphasising on the specifics of the resummation. Sec. 3 then presents the results for off-shell Z, W + and W − production and the respective W ± /Z ratios, double differential in (q T , ∆φ). Sec. 4 summarises our findings. Finally, the appendices collect the details on the specific size and impact of the singlet contributions and leptonic power corrects. They also detail the process-specific hard functions.

Details of the calculation
In this section we detail the construction of our resummed results using the SCET formalism for both the q T and ∆φ observables. These expressions, fully differential in the lepton momenta, are then matched to the respective fixed-order calculation.

Factorisation and fixed-order functions
From the QCD factorization theorem [127], the differential cross section for the Drell-Yan (DY) process can be expressed as where Y L and M L stand for the rapidity and the invariant mass of the final state lepton pair, respectively. Ω L represents the solid angle of one of the final leptons in the lepton-pair rest system.σ ij denotes the partonic cross section which depends on the renormalisation scale µ R as well as the factorisation scale µ F . ff ij is the effective parton luminosity function. It is defined as where f i/N is the parton distribution function (PDF) for the parton i out of the nucleon N ± traveling in the ±z direction. In addition, eq. (2.1) also involves the parameter τ ≡ŝ/s, where s andŝ are the square of the hadronic and partonic colliding energies, respectively. Its minimal value is a function of the invariant mass and transverse momentum to be produced, Therein, M T = M 2 L + q 2 T . In the small q T regime particularly concerned, the partonic cross sectionσ ij can be factorised further. In the context of SCET, one can in principle utilise the decoupling transformation [94] to expressσ ij as a convolution of hard, collinear and soft sectors. However, neither the collinear nor soft sector at this stage is well-defined due to the appearance of the rapidity singularity [128]. Various different regulators have been proposed in the recent years, such as the analytic regulator [96,98], the ∆ regulator [99,100] and exponential regulators [101,102]. In light of its particular performance in the fixedorder calculations, the framework with the exponential regulator will be employed in this work. In this approach, the factorisation formula of eq. (2.1) can be rewritten as [101,102], Therein the light-cone decomposition has been carried out upon the impact parameter b, i.e., Here n ± are two reference vectors satisfying n 2 ± = 0 and n + · n − = 2. For later reference, we also introduce the light-cone coordinate b The integrand of eq. (2.4) comprises three kinds of ingredients. S( b T , µ R , ν) is the soft function encoding all the soft quantum fluctuations surrounding the beam. As a result of the ultraviolet (UV) and rapidity divergences, it possesses explicit dependences on the virtuality scale µ R and the rapidity scale ν. The definition of S reads [101] is the soft renormalisation constant, Y ± stands for the incoming Wilson line along the n ± direction, and τ here denotes the rapidity regulator and b 0 = 2e γ E . Currently, the soft function S( b T , µ R , ν) is known at N 3 LO accuracy [102]. The beam functions B ± contain the collimated contributions along the ±z direction. Introducing the momentum fractions η ± = M L exp(±Y L )/ √ s, the field-operator definition for B + is [101,129] where Z −1 B+ is the collinear renormalisation constant, Z −1 0 is the zero-bin subtractor to remove the softcollinear overlapping contribution, χ +,i is the gauge-invariant building block for the collinear quark q i [130,131] in +z direction. The expression for B − or the anti-quark case can be obtained through changing the light-cone components or field operators in eq. (2.7), respectively. Working in the hierarchy q T Λ QCD , the beam function can be further factorised into the following form [101,129], where the factorisation scale µ F emerges in the right-hand side in the usual way. Considering the µ F dependence of the PDFs will cancel against that in I ij order by order, we suppress µ F in the B ± arguments. Currently, the hard-collinear function I ij is known at the N 3 LO accuracy [116,129,132]. In addition to the soft and beam functions, the factorisation formula in eq. (2.4) also involves the hard sector H V ij (V = γ/Z, W ± ). H V ij can be calculated from the UV renormalised partonic amplitudes, where the sum runs over all the colours and polarisations of the initial partons and includes the appropriate averaging factors. As the result of renormalisation, |M| 2 is free of UV divergences but presents manifest IRC singularities. From the strategy of asymptotic expansion [133,134], these IRC behaviours should be exactly removed by the product of Z B+ , Z B− and Z S and thus H V ij is left as a finite quantity. In contrast to B i ± and S, which are universal within the three processes under consideration in this paper, the hard function depends on the process and encodes the specifics of the hard partonic interaction. In Sec. 2.2 their structure will be detailed.

Hard function: non-singlet and singlet contributions
The Feynman diagrams contributing to the three hard processes under consideration in this paper, off-shell Z or W ± production, can be grouped into two categories according to whether or not the incident quark lines are connected to the EW vertex or not: the singlet and non-singlet contributions (see Fig. 1). Here, the non-singlet contribution of Fig. 1a collects all configurations that connect the external quark lines to the EW vertex while all other configurations are part of the singlet contribution. The latter can be subdivided according to whether they couple to the vector or axial-vector part of the EW vertex, as depicted in Fig. 1b and Fig. 1c, respectively.
Z production. Here, the quarks coupling to the EW vertex are always of the same flavour. Hence, both singlet and non-singlet amplitudes contribute. In the following, their construction and embedding in the resummation framework is detailed. Their numerical impact is fully examined in App. A. In essence, while (a) Non-singlet amplitude.
(b) Singlet amplitude induced by a vector current.
(c) Singlet amplitude induced by an axial-vector current. Figure 1: Representative Feynman diagrams which contribute to the hard functions. The curly and straight lines denote the gluon and quark propagators, respectively. The crossed circle represents the electroweak current operator.
at N 2 LO the impact of the singlet contributions is numerically small, it reaches the same size at N 3 LO as the standard non-singlet contribution and it is essential to include it. All needed amplitudes can in principle be calculated in the full SM and evaluated loop by loop. However, the presence of multiple mass scales complicates the loop integrations considerably, rendering this method somewhat involved. Alternatively, in this work we will employ LEEFT resulting from the SM by integrating out the top quark. For the vector current, due to its conservative nature, the matching from SM onto LEEFT amounts to the re-definitions of the strong coupling α s as well as the field operators [121,122]. 2 As a results, one can collect the effective vector currents as It is apparent that the expressions for SM light quark vector currents are formally retained in V µ γ/Z . For later convenience the currents induced by γ * qq and Z * qq vertices are listed separately. In absence of the EW corrections, one can always distinguish them. The coupling factors g qi γ,V collect the EW coupling constants, (2.11) Here we utilise e as the electromagnetic coupling and s w and c w as the sine and cosine of the weak-mixing angle θ w . T 3 q and Q q are the third component of the weak isospin and charge for the quark q, respectively. The axial-vector current on the other hand, which exists only for Z boson exchange, requires a more delicate matching procedure. As the renormalisability of SM relies on the anomaly cancellation within each quark generation, the removal of the top quark by brute force will break the renormalisation group invariance (RGI) explicitly and thus it necessitates an additional renormalisation constant in the LEEFT for its restoration. To this end, the following operator basis is proposed for the axial-vector effective current [123][124][125][126], where the coupling factor g A again collects the EW coupling constant, i.e.
It is independent of the quark flavour. ∆ ns i stands for the non-singlet operators of the ith quark generation, (2.14) For the first two quark generations, the operators ∆ ns 1,2 are formally identical to those in SM, whilst for the third one, ∆ ns 3 is artificially constructed to facilitate the anomalous cancellation. Here N F represents the number of massless quarks, taking the value N F = 5 throughout this work. O s is the singlet operator, which is defined as In addition, A µ Z also comprises a novel structure C t O s beyond what is present in the SM. In presence of the axial-anomaly, the renormalised singlet O s operator exhibits the explicit dependence on the renormalisation group transformations (RGT). However, this RGT-dependence will be eliminated order by order by the Wilson coefficient C t such that the RGI is still maintained in LEEFT. To calculate C t , one needs to compare the SM amplitudes induced bytγ µ γ 5 t −bγ µ γ 5 b with those from O s in the limit M L m t . Here m t denotes the pole mass of the top quark. In this work we follow the conventions of [125], Here the tree-level result (−1/N F ) balances the O s /N F term in ∆ ns 3 . The singlet contribution starts from the two-loop level and the O(α 2 s ) results can be either straightforwardly read from the axial-anomaly form factor in Ref. [135], or extracted from [124,136]. The third order expressions can be found in Ref. [123]. Particular attention needs to be noted that in Refs. [123,124] the O s operator is renormalised in a different prescription from those in Refs. [125,135]. Their conversion is essential and we present the corresponding details in App. C. Now that we have all ingredients at hand we can finally calculate the hard function. From eq. (2.9), the hard sector can be addressed from the square of IRC-subtracted on-shell amplitudes. In absence of EW corrections, the amplitudes involved can be naturally decomposed into two parts: the hadronic and the leptonic currents. The hadronic current can be expressed using the effective vector and axial-vector currents, V µ Z/γ and A µ Z , respectively. More explicitly, we have where g Σ γ/V is introduced to collect the EW coupling, namely, (2.18) V µ ij and A µ ij are born level amplitudes induced by the vector and axial-vector vertices, respectively. Their expressions read, As the massless QCD interactions conserve chirality, one can always factor them out after the renormalisation and IRC-pole subtraction. In addition, we have also utilised a set of hard coefficients in eq. (2.17) encoding the loop corrections. C ns contains all the non-singlet contributions (see Fig. 1a) and can be extracted from the γ * qq form factor. In the recent years, the γ * qq form factor has been calculated up to three-loop level [114,137,138]. Further, C V s stems from the singlet contribution to the vector current (see Fig. 1b). As Furry's theorem forbids contributions at two-loop order, the lowest order result enters at O(α 3 s ), see [114,138]. Similarly, C A s represents the QCD corrections induced by O s . Its N 2 LO expression can be extracted from Refs. [114,135,137,138], while the logarithmic dependences at third-loop order accuracy can be obtained from the anomalous dimensions. The specific expressions for all three coefficient functions C ns , C V s and C A s are listed in App. C.
Besides the hadronic contributions, the hard function also comprises the leptonic ones, namely the leptonic currents, including the vector boson propagators, (2.20) Therein, µ Z stands for the complex mass of Z boson, which will be introduced properly in Sec. 3.1. In absence of EW corrections, the couplings of the leptonic sector g γ,V,A can be inferred from the leading order vertices, Putting all pieces together, the hard function can be expressed as W ± production. In this process, as the quarks participating in the EW vertex are always of different flavours, only non-singlet diagrams will contribute. In analogy to eqs. (2.17)-(2.22), we write the hard functions as

(2.24)
Here V ji signals the Cabibbo-Kobayashi-Maskawa (CKM) matrix, and µ W represents the complex mass of the W boson.

Resummation
With the help of eq. (2.4), the scales relevant to the small q T regime can be successfully separated. However, in presence of the scale hierarchy, such as q T M L , any fixed-order expansion of eq. (2.4) will suffer from the large logarithmic terms reducing the perturbativity. To address this issue, we employ the solutions of RGEs as well as the rapidity group equations (RaGEs) to carry out the scale evolution. In this way, the fixed-order functions contribute only at their intrinsic scales and the large logarithmic terms arising from the scale hierarchy can collectively be resummed within the exponential functions of the evolution. The procedure is detailed in the following. We begin with the RGEs and RaGEs for each sector. The RGE for C t reads where γ t is the anomalous dimension and its N 3 LO result can be found in Ref. [125]. The evolution equations for C ns , C A s and C V s can be expressed as where Γ cusp stands for the cusp anomalous dimension. Its result is known at the three-loop accuracy [139], whilst the four-loop evaluation has been accomplished very recently [113]. γ h represents the hard anomalous dimensions and can be extracted from [114]. The remaining equations for the soft and beam functions, S and B, read [101] ∂ ∂ ln µ 2 ln S( b T , µ, ν) = Γ cusp ln and (2.28) Therein, γ s is the soft anomalous dimension, on which the accuracy has already arrived at N 3 LO [115]. γ b is the non-cusp anomalous dimension for the beam function, which can be derived from the consistency relationship γ b = (γ s −2γ h )/2. In addition to the ingredients of the RGEs, the rapidity anomalous dimension, γ r , is involved within the RaGEs. Its universality (or the rapidity regulator independence) has been discussed in Refs. [101,140], and the N 3 LO results are detailed in Refs. [102,140]. Table 1: Needed accuracy of the fixed-order inputs to achieve a given logarithmic accuracy of the resummation, in accordance with eq. (1.1).
With all ingredients collected, we can now carry out the resummation. In this work, the resummed distributions are calculated at NLL , N 2 LL , and N 3 LL level. For each level of the logarithmic accuracy, the perturbative order for each ingredient can be found out in eq. (1.1) and has been also summarised in Tab. 1. Note that as the accuracy of the fixed-order functions is counted with respect to Born cross section, only the O(α 3 s ) terms for C V s are needed in this paper. The resummed spectrum is then obtained by substituting the solutions of the RGEs as well as RaGEs into eq. (2.4), giving where the kernels U V,R , respectively, effect the virtuality and rapidity evolutions. Their explicit expressions read . (2.30) Finally, H V,res ij is the part of the hard function participating in the resummation. In W ± production, as M L is the only physical scale in the hard sector, the natural choice of scale µ h can be easily identified (assuming it is related to M L ). H V,res ij then coincides directly with H W ij from eq. (2.23), giving Conversely, the hard function in Z production, contains a second intrinsic scale in addition to M L , m t , originating in the singlet contributions. In order to resum the logarithm associated with their scale hierarchy, ln[m 2 t /M 2 L ], we substitute the Wilson coefficient C t in eq. (2.11) with its resummed version, and then evaluate all other Wilson coefficients C ns , C V s and C A s , at the scale µ h . As a result, the part of the hard function participating in the resummation can be expressed as . (2.33)

Power corrections
To extrapolate eq. (2.29) to the entire phase space, it is crucial to properly account for power corrections stemming from the truncation of the asymptotic series in q T /M L (recalling that eq. (2.4) only contains the leading power contributions). They can impact both the leptonic currents L µ W,Z,γ (in eq. (2.20) and eq. (2.24)) as well as the hadronic ones. To improve the lepton currents, we substitute the L µ W,Z,γ by their pre-expanded results, i.e.
where V = W, Z, γ. Λ µ ,ν accounts for the Lorentz transformation from the leptonic centre of mass reference frame to the rest frame of the initial proton pair or lab frame. Noting that in eq. (2.20) and eq. (2.24), as a result of the q T /M L expansion in the hard sector, L µ V is defined in the lepton centre of mass frame. Thus, the replacement in eq. (2.34) restores the leptonic currents to their form in the exact SM calculation. Their impact is assessed in App. B. It is worth noting that even though the substitution in eq. (2.34) manifestly breaks energy-momentum conservation when coupling the hadronic and the corrected leptonic tensor, it induces no gauge violations in practice. This is because the longitudinal component of the γ, Z, or W ± propagators can always be eliminated by the attached leptonic currents in the massless limit. Alternative methods with the analogous effects can also be found in Refs. [108,119,[141][142][143]. Very recently, a systematic classification has been carried out in Ref. [112] as to the leptonic power corrections. On the other hand, it is also necessary to handle the power corrections to the hadronic currents. In principle, one can incorporate them power by power with the aid of the subleading SCET Lagrangians and Hamiltonians [95,130,144]. Nevertheless, not only can more structures beyond eq. (2.4) be introduced by the subleading vertices, the inhomogeneity in asymptotic series, which arises from the rapidity regulation, may further complicate the investigations [145][146][147]. So for simplicity, this work only considers the leading power factorisation and the according resummation in the hadronic sector, treating the power corrections through the matching to fixed-order.

Matching to fixed-order QCD and observable calculation
Finally, we match the resummation to the exact QCD fixed-order calculation. We therefore introduce the following matching procedure, where σ exp+lpc is the fixed-order expansion of σ res+lpc and the "lpc" addition in the subscript signals the inclusion of the leptonic power corrections. σ f.o. denotes the exact fixed-order perturbative result at the conventional renormalisation and factorisation scales µ R and µ F . Of course, both σ exp+lpc and σ f.o. must be expanded to the same order. dΦ stands for the differential phase space element d 2 q T dY L dM 2 L dΩ L . The transition function f is introduced to assure that the resummation is only active in the asymptotic region, more explicitly, The numerical choice of the matching scale µ Q will be discussed in Sec. 3. The base of the exponent a can be thought of as a free parameter governing the shape of the transition function. We choose to link it to µ Q and define it as a = exp(µ Q /q ref T ), q ref T is taken to be 1 GeV. It is important to note that although f is smaller than unity for all physical q T , the precise choice of µ Q will ensure that it differs from unity by amounts much smaller than one permille in a wide region around the Sudakov peak where the resummation effects dominate. Likewise, q T = µ Q does not mark the endpoint of the resummation region, but where the resummation effects are faded out to half their intrinsic size, f (1) = 1 2 . The functional form then induces a symmetric reduction of the influence of the resummation in the region [µ Q − ∆q T , µ Q + ∆q T ], wherein ∆q T = q ref T ln((1 − f )/f ) and f is the fraction of reduction. Thus, with the above value for q ref T the resummation effects are reduced from 99% to 1% of their actual size in the interval µ Q ∓ 4.6 GeV. The functional form of the base of the exponent a ensures that the width of the transition region, ∆q T , is directly proportional to q ref T , with q ref T taking the role of an easily interpretable shape parameter. When q T enters the hadronic regime Λ QCD ∼ q T M L , the non-perturbative contributions become non-negligible. Hence, a modification factor parametrising the non-perturbative dynamics of this regime, S mod (b 2 T ), encoding all the hadronic contributions within the low q T region [148,149], should be introduced. In the recent years, a number of parametrisations and input extractions of S mod (b 2 T ) have been discussed in the literatures [15,16,[150][151][152][153][154]. However, as the main concern of this work is with the perturbative domain, we ignore S mod (b 2 T ) here and leave its incorporation to the future investigations. Taking the integral of eq. (2.35) together with the customised theta and delta functions, one is able to implement the fiducial cuts and carry out the observable calculations, more explicitly, where G is a function of the phase space encoding the fiducial cuts. O stands for a generic observable and F O is its corresponding functional definition. This work focusses on the q T and ∆φ spectra, giving Here p 1(2),T denotes the transverse momentum of the final (anti-)lepton in the lab frame. In particular, both variables are purely transverse and are thus well-defined for both W and Z production in a realistic detector environment where only the transverse part of the neutrino's momentum is observable through as missing transverse momentum. Further, in calculating the phase-space integral of eq. (2.37), we retain the full kinematic dependences in G(Φ) and F O (Φ) and do not carry out any multipole-expansions.

Numerical Results
In this section we are discussing numerical results obtained with the methods outlined in the previous section for a representative fiducial region.

Setup and fiducial region
Throughout this paper, all distributions for W ± production are calculated in the G µ scheme while those for Z production are calculated in the α(m Z ) scheme. We use the following input parameters [155] The width of the top quark as well as the mass and the width of the Higgs boson have no phenomenological impact and are neglected. All other particles are considered massless. We work in the complex-mass scheme [156], with the complex masses and mixing angles defined through The electromagnetic coupling constant is defined as in the G µ and α(m Z ) schemes, respectively. Other input schemes, in particular for Z boson production, may be chosen in order to improve the description of the correlations of the lepton pair system and the hadronic recoil, e.g. by using the weak mixing angle as an input parameter [157]. Furthermore, we use the NNPDF31 nnlo as 0118 [13] parton distribution functions, interfaced through LHAPDF [158], with the corresponding value of the strong coupling α s (m Z ) = 0.118. In accordance with this choice we neglect all photon induced contributions. The CKM matrix is chosen as diagonal, i.e. V ij = δ ij . To carry out the resummation entering our matched results of eqs. (2.35) and (2.37), we employ the CUBA library [159,160] performing the respective integrations and make use of LHAPDF evolving the factorisation scale. To tackle the harmonic poly-logarithms participating in the I ij function (see eq. (2.8)), the package HPOLY [161] is used. The LO hard functions in eqs. (2.23) and (2.22) are computed with FeynCalc [162][163][164] and FeynArts [165,166]. All remaining parts of the hard function, in particular the non-singlet and singlet contributions C ns , C V s , and C A s as well as the top-quark contribution C t , are implemented analytically. The fixed-order contribution to eq. (2.35) and, hence, the matched result is computed using SHERPA [167,168] in combination with OPENLOOPS [169,170]. In this framework, renormalised virtual amplitudes are provided by OPENLOOPS, which uses COLLIER tensor reduction library [171] as well as CUTTOOLS [172] together with the ONELOOP library [173]. All remaining tasks, i.e. the bookkeeping of partonic subprocesses, phase-space integration, and the subtraction of all infrared singularities, are provided by SHERPA using the matrix element generator AMEGIC [174][175][176]. In estimating the theoretical uncertainties, we emphasise two aspects. The first one is the uncertainty arising from higher order corrections, which can be estimated by examining the sensitivity of the results to the scale variations. As presented in Sec. 2.5 the matched result depends on a set of auxiliary scales, originating both in the resummation and the fixed-order calculation, {µ, ν} ≡ {µ t , µ h , µ b+ , µ b− , µ s , ν b+ , ν b− , µ R , µ F }. During the numerical implementation, we take µ b = µ b± as well as ν b = ν b± in accordance with identical initial state particles, and set µ R = µ F throughout to simplify the matching procedure of the fixed-order full QCD calculation to the results obtained in the soft-collinear effective theory. We set their default values, well away from the non-perturbative regime, to Please note that while M L = m in Z production, it is equal to invariant chanrged-lepton-neutrino mass in W production. For the sake of reducing the logarithmic contribution in the H, B, S functions, the default intrinsic scales are taken as It is worth noting that with the choice of µ def can not be carried out straightforwardly due to presence of the Landau singularity at small momentum transfers, or large impact parameters, of the strong coupling constant α s . To cope with this issue, the prescription in Ref. [177] is adopted in this work to suppress the higher-b T influences. For estimating the uncertainties from the choices in eqs. (3.3) and (3.4), we vary the scales µ t,h,b,s , µ R,F and ν b,s independently to twice and half their default values, and then combine the deviations from the results using their above defined default values in the quadrature. The thus constructed uncertainty estimate is denoted δ scale hereafter. This leaves the uncertainty originating in the matching parameter µ Q . As our default we take µ Q = 16 GeV. The effectiveness of this choice will be illustrated and discussed in Sec. 3.2. For the estimation of its error, we alter the value in the interval [15,17] GeV and denote the uncertainty from this by δ mat . Together with δ scale , the total uncertainty can be evaluated as δ tot = δ 2 scale + δ 2 mat . We compute our results in a fiducial phase space that is loosely modelled after the respective ATLAS and CMS measurement regions. It differs slightly depending on the process and the physics objects. Leptons are required to have a transverse momentum p T of larger than 20 GeV and an absolute pseudo-rapidity of less than 2.4. In the dilepton channel we require an opposite-sign same-flavour pair with an invariant mass m of more than 80 GeV and less than 100 GeV to effectively restrict it to the Z pole. In both lepton-plusmissing-transverse-momentum channels we require one such lepton of the respective charge and a missing transverse momentum p T of more than 20 GeV. For simplicity, the missing transverse momentum is equated with transverse momentum of the neutrino. In addition, the missing transverse momentum and the charged lepton have to have a transverse mass m T = 2 p T p T (1 − cos ∆φ) of more than 40 GeV, where ∆φ is the opening angle between the two in the transverse plane. These definitions are summarised in Tab. 2. As we work at strict leading order in the electroweak sector, questions about the precise lepton definition do not arise. While the resummation part of the matched result is calculated for the specific observable under consideration, the fixed-order part is generated as conventional collider events and analysed using RIVET [178,179] Fig. 2 now details the distribution at N 2 LO, that is including terms up to O(α 2 s ), in the transverse momentum of the leptonic system q T and the opening angle between the two leptons in the transverse plane ∆φ. A regularising cut of q T > 1 GeV is applied. As can be seen, the different sets of fiducial cuts necessitated by the differing measurable physics objects in each channel, induce differing coverages of the q T -∆φ-plane. Although the fiducial regions overlap well at very small q T and ∆φ, they illustrate that ratios between Z and W channels for the purpose of reweighting one channel to the other should not be taken single-differentially in q T but need to take the constrained internal dynamics of the lepton system into account. Finally, Fig. 3 shows the qualitative behaviour of the single-differential behaviour of the matched cross section in q T and ∆φ. In both distributions the typical resummed shape can be observed, exhibiting a Sudakov    peak close to the singular point (q T , ∆φ) = (0, 180 • ). It must be noted though, that the expression of the Sudakov peak in the ∆φ spectrum strongly depends on the precise definition of the fiducial region. With these observations at hand, we define three different regions of interest for both spectra: (I) a region between the singular point and the Sudakov peak, (II) a region containing the Sudakov peak, and (III) a region beyond the Sudakov peak that includes the part of the spectrum where the resummation becomes unimportant.

Validation
To begin with the examination of our result we first present a comparison of the resummed result, expanded to next-to-leading order (NLO s ) and next-to-next-to-leading order (N 2 LO s ), i.e. truncating at O(α s ) and O(α 2 s ), respectively, to the exact full QCD result at the same order, NLO and N 2 LO. It needs to be noted that, contrary to the complete resummed result, the (expanded) fixed-order results diverge as q T → 0. We thus, impose a cut of q T > 1 GeV for the validation in this section, both for the q T and ∆φ spectra discussed in the following. The N 3 LO s expansion of the resummed result will be compared with the corresponding full QCD computation in a future study. Fig. 4 shows this comparison for the q T spectra in the three slices in ∆φ introduced earlier. In all three slices the agreement between the exact QCD result and the expanded resummed result derived from the softcollinear effective theory is excellent for all three processes, Z, W + , and W − production, under investigation. It is interesting to note that within a few percent the approximation in the soft-collinear effective theory holds to sizeable distances from the singular point at (q T , ∆φ) = (0 GeV, 180 • ). In particular, deviations do not exceed 2% for q T up to 10 GeV for all ∆φ ranges, neither at NLO nor N 2 LO. These observations hold independent of the factorisation and renormalisation scale, as is indicated by the coincidence of the shown        scale variation bands. The situation is slightly different for the ∆φ distributions in the three q T slices depicted in Fig. 5. Here, the fixed-order expansion of the resummation very well coincides with the exact QCD calculation only in the first two slices for q T smaller than 6 GeV. In the last q T slice, well away from the singular point, the agreement is worse, ranging from around 5% to 10%, with the best agreement unsurprisingly found on the Sudakov peak. Small effects due to limited statistics can be observed near the phase space boundaries, but do not impact our findings.  Figure 5: Comparison of the fixed-order ∆φ spectra in all three processes. N (2) LO denotes the fixed-order full QCD perturbative result, while N (2) LO s is the fixed-order expansion of the SCET-based resummation.

Resummation improved results
We now turn to examine our full resummation improved results, calculated according to eq. (2.37). Contrary to the fixed-order evaluation of the previous section, the behaviour of the cross section as we approach the singular point at (q T , ∆φ) = (0, 180 • ) is regular and we can evaluate the complete (q T , ∆φ) plane. To arrive at our resummation improved results of this section, however, the governing eq. (2.35) still contains two terms that are separately diverging as the singular point is approached, originating in the behaviour of  Figure 6: Single-differential cross section in q T for all three processes. We compare the full QCD NLO and N 2 LO distributions to the resummation improved results, NLL +NLO, N 2 LL +N 2 LO, and N 3 LL +N 2 LO.
the full QCD fixed-order calculation and the fixed-order expansion of the resummation in the soft-collinear effective theory. To regulate their behaviour, following the results of Sec. 3.2, we set them to be exactly equal for q T < 1 GeV, leaving only the resummed result similar to the treatment in Refs. [111,180]. Having said that, the most important question to answer, is the choice of matching scale µ Q . For this, we follow two arguments to guide us to our choice of matching scale. Firstly, we want to restrict the resummation to the asymptotic regime where its intrinsic approximations are valid. Recalling that the resummation in eq. (2.29) is derived from the factorisation of eq. (2.4) where only the leading contributions in an expansion in q T /M L are kept, eq. (2.29) is valid only in the small q T regime and thus should be disabled beyond it. To this end, the transition function f (x) of eq. (2.36) is employed to provide a smooth transition from the resummed spectra to the fixed-order contribution in the matched result of eq. (2.35) within the interval µ Q ∓ ∆q T with ∆q T ≈ 4.6 GeV with our choice of scales. We note that the value of µ Q is related to the effective range of the q T /M L expansion. Following the spirit of [180] we determine that range by comparing in Fig. 4 the N 2 LO result with the corresponding expansion of the resummation N 2 LO s . To be precise, we require the fixed-order expansion of the SCET approximation to deviate from the exact result by no more than 20%. For all three processes and ∆φ slices we extract similar values, leading to a common choice of matching scale of µ Q = 16 GeV. Secondly, we want the additional corrections introduced by the resummation with respect to the fixed-order calculation to be small or negligible at the matching scale. To evaluate this requirement, Fig. 6 is of particular interest. Here we observe that at values around the chosen matching scale µ Q = 16 GeV the resummation improved results coincide with the pure fixed-order one to better than 3%. At this point, the reader is reminded that, although we are not comparing the pure resummation with the full QCD result but instead a result where the resummation at the scale µ Q is already subjected to the suppression function f (q T ) of eq. (2.36), the suppression function has the value f (µ Q ) = 0.5. Therefore, we still find that the resummation and the fixed-order result still agree to better than 5% and resummation effects are no longer important. It is interesting to note that this observation holds for both N 2 LL +N 2 LO and N 3 LL +N 2 LO. The situation is slightly different for NLL +NLO. However, since this result is mainly included to illustrate the progression of the increased accuracy of our calculation we choose the same value for µ Q . Analytically this can be understood in the following. In the above argument we are essentially evaluating the relative size of the contribution the resummation is supplying beyond the accuracy of the fixed-order calculation. These terms are of O(α 2 s L 4 + α 2 s L 3 + . . .) for the NLL +NLO matched result, while they are of O(α 3 s L 6 + α s L 5 + . . .) for the N 2 LL +N 2 LO and N 3 LL +N 2 LO calculations, L = log(q T /M L ). Now while q T is small, these contributions are of the same order. Choosing a µ Q sufficiently removed from the singular point, such that the ratio q T /M L is of O(1), L follows a different power counting. Thus, the additional terms induced by the resummation with respect to the fixed-order calculation are indeed of higher-order in N 2 LL +N 2 LO and N 3 LL +N 2 LO than in NLL +NLO. With this choice of resummation scale the single-differential distributions in the leptonic transverse opening angle ∆φ similarly receive substantial resummation effects, as shown in Fig. (7). Since the suppression function f acts in another variable, no clear transition from one regime to the other can be observed at any order. We observe, however, that while at our highest order, N 3 LL +N 2 LO, all three processes behave  Figure 7: Single-differential cross section in ∆φ for all three processes. We compare the full QCD NLO and N 2 LO distributions to the resummation improved results, NLL +NLO, N 2 LL +N 2 LO, and N 3 LL +N 2 LO.
extremely similarly and receive very similar resummation induced corrections to the N 2 LO result, this is markedly different at NLL +NLO accuracy. Here, the W production channels receive somewhat larger corrections in the region between ∆φ = 150 • and ∆φ = 165 • than in Z production.
With this global picture examined, we now turn to the aim of our study, the double-differential distributions. We again begin with the q T spectra in our three chosen slices of ∆φ in Fig. 8. While only the first slice for ∆φ > 175.5 • contains the singular point, all slices are close enough to the singularity that a Sudakov peak is formed. The description of this resummation region depends strongly on the order of logarithmic corrections included. While the central values do not change significantly order-by-order, validating our choice for the central scales involved in the resummation, the estimated uncertainty steadily decreases when going to higher orders. To be definite, the uncertainty in our lowest accuracy calculation (NLL ) ranges from −20% to +15% around the central value around the Sudakov peak, independent of the ∆φ slice, and involves a sizeable shape uncertainty as well. This shape uncertainty in particular differs in the region below the Sudakov peaks depending on ∆φ. Both uncertainties are reduced greatly when higher-order logarithmic corrections are included. At N 2 LL they amount to −4% to +7% while at N 3 LL they are reduced to −2% to +4%. Moving towards larger q T , the matching uncertainty tends to become comparable to the perturbative uncertainties in the resummation, dominating in particular NLL +NLO calculation around µ Q . As discussed above, the matching scale was chosen based on arguments for the highest precision calculation in this study and significant contributions of the resummation beyond the NLO fixed-order accuracy were found. Hence, this finding is not surprising. On the contrary, the matching uncertainty is substantially reduced in both the N 2 LL +N 2 LO and N 3 LL +N 2 LO calculations, not exceeding ±1% in the latter. At even higher transverse momenta the fixed-order calculation dominates the spectrum and its usual behaviour and uncertainties are recovered. Finally, we examine the resummation improved results for the ∆φ spectra in the three chosen q T regions.
Recall that the first region with transverse momenta smaller than 2 GeV resides entirely below the Sudakov peak in the transverse momentum spectrum, while the second one contains the peak, and the third region with q T > 6 GeV resides entirely beyond the Sudakov peak. Consequently, very different behaviour can be observed despite the per-bin-cross sections being of similar orders of magnitude. The results of our computation for the ∆φ spectra are shown in Fig. 9. As all three processes exhibit a very similar behaviour we continue to discuss them simultaneously. In the lowest q T region, which probes the region below the Sudakov peak in the q T spectrum, near back-toback topologies are favoured unsurprisingly and resummation effects dominate the calculation throughout. Consequently, as observed before, the central values of our three predictions of increasing logarithmic accuracies agree very well. Their uncertainties have nearly no ∆φ dependence and are steadily decreasing with the increasing accuracy of the resummation, reaching −1% to +2% in the N 3 LL case.
In the intermediate q T slice, focussing on the region around the Sudakov peak in the q T spectrum, a peaked structure is developing. Still, the distribution is dominated by resummation effects, leading again very well agreeing central values, their uncertainties being dictated by the order of the logarithms included in the exponent. We, thus, find a variation of ±1% in our most accurate calculation. This time, however, there is  Figure 8: Double-differential cross section in q T and three slices of ∆φ for all three processes. We present the resummation improved results at NLL +NLO, N 2 LL +N 2 LO, and N 3 LL +N 2 LO accuracy.
a slight shape to these scale uncertainties, predominantly in the back-to-back region. The third q T slice, located entirely above the q T Sudakov peak, now sees a stronger impact of the full QCD fixed-order calculations. Hence we find the central values no longer agree, partially reflecting the difference in cross section between the NLO and N 2 LO predictions for the spectrum. The increasing importance of the full QCD fixed-order part for smaller ∆φ likewise explains the shape corrections at higher order. Still, the uncertainties are much reduced in our highest precision calculation at N 3 LL +N 2 LO accuracy, being smaller than ±1% throughout. Finally, please note that the smallest ∆φ bin in the lower two q T slices carries almost no cross section and, thus, suffers from larger statistical uncertainties.

The W ± /Z and W + /W − correlations
With the double-differential resummation-improved cross sections of the previous section at hand we finally turn towards cross section ratios. They are useful to obtain high-precision W + and W − production data, by measuring the fully and precisely reconstructible cross sections in Z production and applying the following theory predictions, see e.g. Refs. [12]. Similarly, the W + to W − ratio is of interest for PDF extractions, see e.g. [13]. A key question, however, is how to determine the uncertainties of such a ratio.  Figure 9: Double-differential cross section in ∆φ and three slices of q T for all three processes. We present the resummation improved results at NLL +NLO, N 2 LL +N 2 LO, and N 3 LL +N 2 LO accuracy.
is the fundamental lack of statistical interpretation of the theoretical uncertainties on (multi-)differential cross sections as presented so far. Thus, in particular various assumptions about their (non-)correlation in the numerator and denominator of the R W ± /Z have to be made. In the following, we present results for the following three correlation assumptions: • Uncorrelated. The scales of the W ± and Z processes are assumed to be completely uncorrelated. All scales in the numerator and denominator are varied independently. This corresponds to the assumption that both processes have no common structure in the form of their higher-order corrections or input functions, such as the PDFs. As this is known not to be the case, the uncertainties on the ratio obtained this way are likely to be severely overestimated.
• Fully correlated. All scales of the W ± and Z processes are assumed to be completely correlated. They are thus varied by a common factor in the numerator and denominator simultaneously. This corresponds to the assumption that both processes have exactly the same structure in the form of their higher-order corrections and input functions, such as the PDFs. As this is known not to be the case, the uncertainties on the ratio obtained this way are likely to be severely underestimated.
• Partially correlated. A careful analysis of the internal structure of the higher-order corrections to W ± and Z production allows to carefully assess which corrections have identical (or at least very similar)   Figure 11: Double-differential cross section ratios in ∆φ and three slices of q T for all three processes. We present the resummation improved results at NLL +NLO, N 2 LL +N 2 LO, and N 3 LL +N 2 LO accuracy.  R W−/Z Figure 12: Double-differential cross section ratios R W + /Z and R W − /Z in q T and the two ∆φ slices closest to the singular region. We present the resummed results at N 3 LL accuracy for the default parameters and fiducial regions used in all physical predictions (black), for the default parameters and the Z production phase space G Z adapted to the W ± production phase space G W ± (green), for the Z boson mass set equal to the W boson mass and default fiducial regions (blue), and for both adaptations of the Z production parameters (red), see text for details.
structures and which differ. This allows, to first approximation, to select a subset of scales to fully correlate, uncorrelating the rest. Following a detailed analysis of the derivations of Sec. 2, we find that the singlet contributions to Z production are numerically small, see App. A. In consequence, both the hard and soft functions for W ± and Z decays show the same dependence on the respective scales. Differences, however, occur for the beam functions, arising in the different composition of initial states in all three processes. We thus choose to fully correlate the variation of all scales except for µ b and ν b , which we fully uncorrelate.
Figs. 10 and 11 now show the ratios R W + /Z , R W − /Z , and R W + /W − for all three definitions of their uncertainties detailed above. In the following, we will discuss the different features of both their central values and uncertainties separately.
Central values. The central values of the cross section ratios are of course unaffected by the precise definition of the uncertainty band. Instead, they are largely determined by the slightly different location of the Sudakov peak induced by the mass difference, the different x-dependence of the contributing parton distributions, and the slightly different fiducial phase spaces. In addition, they only exhibit a small dependence on the perturbative order at which they are calculated. In fact, with respect to the perturbative stability of these ratios, we observe only minor corrections on the level of up to 2% in R W + /Z and R W − /Z , and much smaller in R W + /W − , when increasing the intrinsic accuracy of the resummation-improved calculation from NLL +NLO to N 3 LL +N 2 LO. Further, the cross section ratios depend only weakly on the transverse momentum of the reconstructed vector boson for q T > 5 GeV. Below that value, in the vicinity and left flank of the Sudakov peak, all three ratios exhibit a marked increase. As alluded to earlier, this increase is induced by the differing precise locations of the Sudakov peak in each process. Fig. 12 investigates this phenomenon more closely, by a) setting complex mass of the Z boson to that of the W boson in the propagator only, keeping all other parameters at the default values, b) replacing the fiducial phase space definition for the Z production channel, G Z , by that of the W ± channel, G W ± with the ∓ taking the role of the neutrino, and c) applying both modifications a) and b). R W−/Z Figure 13: Double-differential cross section ratios R W + /Z and R W − /Z in ∆φ and the two q T slices closest to the singular region. We present the resummed results at N 3 LL accuracy for the default parameters and fiducial regions used in all physical predictions (black), for the default parameters and the Z production phase space G Z adapted to the W ± production phase space G W ± (green), for the Z boson mass set equal to the W boson mass and default fiducial regions (blue), and for both adaptations of the Z production parameters (red), see text for details.
This leaves the differences due to the participating parton fluxes and the different spin-structures in the underlying EW couplings. We find that in the q T spectra the majority of the effect is induced by the differing fiducial regions, with smaller additional corrections stemming from the different W and Z boson masses. With both effects accounted for, the ratios are nearly q T and ∆φ independent, with the remaining small deviations attributed to the PDFs and the different spin-structures of the underlying EW coupling. It needs to be noted, though, that the difference in fiducial phase spaces between W ± and Z measurements in the W -boson mass measurement [12] was larger than the one used here. Hence, the effect can be estimated to have been larger in that phase space as well.
The ∆φ spectra show a stronger variation of the central value of the cross section ratio, increasing to up to a factor of three above their value far away from the Sudakov peak in both q T and ∆φ, in particular in the first two regions. This increase, however, appears on the far side of the peak away from the back-to-back region, in contrast to the increase observed in the q T spectra. Nonetheless, its origin can be traced to the same factors as for the q T spectra in Fig. 13, the different definitions of the fiducial phase space in W and Z production and the different W and Z boson masses. To be specific, we observe that when both effects are accounted for the ratio is nearly independent of both q T and ∆φ. This also means that the remaining PDF dependence is small. At this point it is important to note that although both the fiducial phase spaces in W + and W − appear to be the same, they are not. The reason is that, in terms of spin-correlation the anti-neutrino produced in the decay of the W − takes the role of the charged lepton in the decay of the W + , but not in the observable definition. This is compounded by the fact that, out of the three processes under consideration here, W + and W − show the largest divergence of the contributing partonic fluxes, imparting differing rapidity distributions on the produced boson and, thus, slightly different effects of the fiducial cuts. These factors add up to explain the remaining small, but non-negligible q T and ∆φ dependence of R W + /W − . Finally, the ratios R W + /Z , R W − /Z , and R W + /W − are nearly flat in the third region containing events with q T > 6 GeV, its only structure being induced mostly by the difference in the fiducial phase space in W and Z boson production.
Uncertainties. The uncertainty of the cross section ratios follows the pattern laid out in their definition: while the fully-correlated case leads to vanishingly small uncertainties, smaller than 1% in most regions (in fact, the largest surviving uncertainty is related to the matching procedure), the fully-uncorrelated case lies on the opposite end of the spectrum with uncertainties of ±4% for our best calculation at N 3 LL +N 2 LO accuracy both in the fixed-order region and the resummation region. The partially-correlated ansatz so far yields the, in our judgement, most reliable result, ranging from ±1% in the fixed-order region where the respective scales are correlated and ±3% in the resummation region. In particular, it is interesting to note that the beam function uncertainties are the driving force of the resummation uncertainties overall, reinforcing the difference in contributing parton fluxes as a driving factor for the details of the ratio overall. Similarly, we observe that, apart from the fully-correlated uncertainty estimate, the uncertainty for R W + /W − largely follows the pattern of R W + /Z and R W − /Z both qualitatively and quantitatively.

Conclusions
In this paper we have computed the single-differential q T and ∆φ as well as the double-differential (q T , ∆φ) spectra for inclusive Z, W + , and W − production in the experimentally accessible fiducial phase space up to N 3 LL +N 2 LO accuracy resumming small transverse momentum logarithms. Besides the essential inclusion of the third-order soft and beam functions, the resummation features the incorporation of leptonic power corrections and the singlet contribution into the hard sector. The leptonic power corrections have been found to extend the region of validity of the approximate SCET result. The singlet contributions, on the other hand, characterised by topologies where the external quarks do not directly couple to the electroweak gauge boson, enter the Z boson production process at second and third order in α s , and have been found to yield corrections of similar size as the non-singlet third order ones, and are thus non-negligible at N 3 LL . In our numerical evaluation we first confronted the approximate results derived from the SCET with the exact ones and excellent agreement has been observed in the asymptotic regime. We then computed the resummation-improved single-and multidifferential distributions at NLL +NLO, N 2 LL +N 2 LO and N 3 LL +N 2 LO and found excellent perturbative convergence in the asymptotic regime, i.e. the higher order predictions and their estimated uncertainties are fully contained in the lower order uncertainty band. Further, the respective uncertainties themselves are systematically reduced to the level around 4% or less at N 3 LL +N 2 LO. In addition, we computed the ratios R W + /Z , R W − /Z , and R W + /W − of these calculations and estimated their uncertainty assuming no correlation, full correlation, and, as our best prediction, a partial correlation of the scale variations in the numerator and denominator making up the uncertainty. For the partial correlation case in particular, a careful assessment of the internal structure of our calculation allowed to identify identical (or very similar) components and structures that differed between the three different processes. Consequently, the scales used to estimate the uncertainties originating in similar components were correlated while scales used to estimate the uncertainties originating in differing structures were varied independently. The main driver for the differences were identified to be related to the incident beams and the partonic fluxes initiating the respective processes. In summary, we determine the ratios with relative uncertainties of less than 1%, rising to 3-4% in the resummation region at N 3 LL +N 2 LO accuracy. The shape of the ratios, although largely perturbatively stable, is not constant but shows a strong impact of the fact that the Sudakov peak is located at slightly different positions in all three processes. This observation is the consequence of three main factors: 1) the difference in the W and Z boson masses, 2) the different partonic fluxes contributing to the three processes and the different (x, Q)-dependence of the u and d valence quarks in particular, and 3) the difference in the fiducial regions for + − and ± + p T final states.
The latter is, in fact, the dominating factor in the ∆φ and low-q T -dependence of the R W + /Z and R W − /Z ratios.
With the presented calculation at hand, precise predictions for both absolute fiducial multidifferential cross sections and their ratios that are vital for the LHC's precision measurements programme can be made. Nonetheless, important contributions from higher-order corrections originating in the electroweak sector of the Standard Model are not included in our calculation yet, and we leave their investigation to a future publication. Similarly, we have for now not included non-perturbative corrections that are expected to give non-negligible contributions in particular at small transverse momenta or large azimuthal separations. Their  reliable modelling is intricate, including both flavour and x-dependent contributions, and goes beyond the scope of this paper. It will also be addressed in a future publication.

Acknowledgements
WJ is grateful to Xuan Chen and Markus Ebert for helpful discussions on the matching procedure. WJ also thanks Konstantin G. Chetyrkin

A Impact of the singlet contributions
In this appendix we present a quantitative discussion of the singlet contributions to the hard function. To this end, we re-express the hadronic current H µ,ij of (2.17) in the following form,     Z,V for the uū and dd partonic channels. As the singlet contribution to the vector current is forbidden at two loop level by C-parity, we plot only the third order therein. The numerical results computed in three different ways are displayed: 1) the O(α 3 s ) correction with both the non-singlet coefficient C ns and the singlet coefficient C V s included; 2) the same O(α 3 s ) correction including only the non-singlet contribution C ns ; and 3) the O(α 3 s ) singlet contribution C V s . It is seen that for the entire µ range, the both the real and the imaginary part of the full O(α 3 s ) correction (in blue lines) almost coincide with the pure non-singlet ones (in red lines), and the singlet contributions (in green lines) are negligible. This indicates that at least up to O(α 3 s ), the role played by the singlet terms in the vector hadronic current is of little phenomenological impact. Moreover, it can also be observed that the full and non-singlet results of uū initial state are curved in the opposite direction to the dd ones. This is caused by the different charges and weak isospins of the u and d quarks. As shown in eq. (2.17), given the negligible singlet vector terms, the vector current hard coefficients between different initial states can be related as follows Here the flipping signs account for the different curvatures as observed, and the proportions derived in the last step give the relationship of the magnitudes of the uū and dd initiated results. In Fig. 15, the magnitudes of the real and imaginary parts of axial-vector current hard coefficients are depicted. As the singlet contribution here starts at O(α 2 s ), we show both h Z,A and h Z,A . In complete analogy to the vector current coefficients, we here also graph the three types of outputs therein: 1) the full result calculated as eq. (2.17); 2) the non-singlet one including C ns only; and 3) the singlet one which is the difference between the previous two cases. Contrary to Fig. 14, the axial-vector singlet terms give an unignorable contribution here for all values of µ, but in particular for values of µ close to 0.5 or 1.5 where either the real or imaginary part of the non-singlet contribution vanishes. In the region µ ∼ m Z , which is the default hard scale taken in the resummation in this paper (see eq. (3.4)), the singlet terms can account for around 20% of the real contributions of the full results in O(α 2 s ), and nearly 40% in O(α 3 s ). In the imaginary parts, it is also noted that more than 10% of [h (2),uū Z,A ] and [h (2),dd Z,A ] is made up of the singlet terms for the majority of the µ range, and although at the third order accuracy the singlet terms experience zeros in the vicinity of µ = m Z , its proportions can recover to ∼ 10% when µ/m Z exceeds 2. Furthermore, one interesting phenomenon in Fig. 15 is that the singlet contributions remain the same in the uū and dd initial states, but the non-singlet lines are curved in the different directions. The reason can be found in eq. (2.17). Therein the non-singlet contribution C ns is directly proportional to the weak isospin of the initial parton, T qi , whilst the singlet terms C t C A s + C ns /N F are universal for all initial states. Therefore, after excluding the singlet terms, the h Z,A s in Fig. 16. It is seen that for the majority of the µ range, the singlet contribution takes the comparable magnitude to the third-order nonsinglet result. Especially in the µ ∼ m Z region, the singlet terms approach [h  Table. 1, the resummation at N 2 LL requires the O(α 2 s ) hard function, whilst the N 3 LL accuracy needs those up to the third order. So if only a precision of the order of 10% of the O(α 2 s ) correction is needed, the singlet terms could be neglected at the N 2 LL accuracy in a numerically approximate sense. Nevertheless, as the complete singlet corrections are of the same magnitude as, if not larger than, the third-order non-singlet corrections, their inclusion is mandatory for a meaningful and robust N 3 LL resummation.

B Impact of leptonic power corrections
In eq. (2.34), the Lorenz-transformation matrix Λ ( q T ) has been incorporated into the hard functions for including the leptonic power corrections. In this part, we will investigate its numerical influences on the double-differential observable d 2 σ/(dq T d∆φ). To this end, we define the ratio of the approximate results to the exact ones as follows,  Figure 17: Numeric impacts of the leptonic power corrections on the q T spectra.
or without Λ ( q T ) contributions. While the result experienced eq. (2.34) is still named N (m) LO s , same as those in Figs. (4)(5), those excluding Λ ( q T ) effects are labelled as "N (m) LO s w/o. lpc". Here the superscript "m" specifies the expansion order in α s . Throughout this section, the results in former case will be depicted in the solid lines, and to distinguish, the dashed ones illustrate the later case. In Fig. 17, we exhibit the numerical results for κ(q T , ∆φ) after integrating out ∆φ over the following four intervals: In the small q T regime, the κ distributions in the four slices all approach the unity and the differences due to the leptonic power corrections are insensible. This phenomenon is in agreement with the observations in Sec. 3.2 and also validates the leading power factorization in eq. (2.4). However, with the increase in the q T , the power corrections start to manifest themselves. For instance, in the intermediate region q T ∼ 20 GeV, one can find that the κ spectra obviously deviate from the unity grid and furthermore the discrepancies arising from the Λ ( q T ) incorporation emerge. One interesting phenomenon is that differing from the first three slices, the κ spectra in the last slice are fairly sensitive to the leptonic power corrections. To be specific, around the point q T = 25 GeV, the Λ ( q T ) incorporations in the first three slices improve the dashed lines by nearly 1% towards the "f.o." ones, whilst it escalates to 5% for the last slice. To interpret this, note the topological configuration of final leptons in the fourth slice considerably differs from the others. For the first three slices, as required by the small value of (π − ∆φ), the final leptons are almost in the back to back configuration. So for the q T ≥ 25 GeV region, the direction of q T tends to be pp → Z/γ * → µ + µ − (13 TeV LHC)  Figure 18: Numeric impacts of the leptonic power corrections on the ∆φ spectra.
aligned withˆ p ,T in the transverse plane, so as to avoid the energetic recoil enhancing the (π−∆φ) value. Herê p ,T represents the (anti-)lepton momentum in the rest frame of the lepton pair. Given the limit q T ˆ p ,T , the Λ ( q T ) matrix acts only on the time-like and the longitudinal components of L V , whereas the later case is only able to polarize perpendicularly toˆ p ,T in the massless limit. 3 After the contraction, the Λ ( q T ) matrix is effectively reduced to the metric tensor and therefore gives rise to mild influences in the first three slices. The situation in the last slice is however different. A variety of (π − ∆φ) therein can encourage the stronger recoils against lepton momenta, so that q T is capable of developing sizable perpendicular component with respect toˆ p ,T , which in turn permits Λ ( q T ) non-trivially coupled with L V . With particular attention is that except for the boundary bins which suffer from statistical issues, the κ values are all enhanced towards the unity grid after incorporating the leptonic power correction matrix Λ ( q T ). This indicates that Λ ( q T ) indeed compensates for the power expansions and therefore demonstrates the effectiveness of our strategy. Additionally, Fig. 18 illustrates κ(q T , ∆φ)s after integrating out q T over the four intervals as follows: 1) 1 GeV< q T < 2 GeV; 2) 2 GeV< q T < 6 GeV; 3) 6 GeV< q T < 20 GeV; 4) 20 GeV< q T < 40 GeV. 3 Generically, one can decompose the transformation matrix as Λ µ ,ν ( q T ) ≡ Λ µ z Λ µ x Λ µ y . Here Λ µ z represents the boosttransformation along the colliding beam direction. Λ µ x(y) stands for those in the transverse plane. Here we only consider the transformation in the transverse plane, since the boost along the beam direction has vanishing impacts after being contracted against the leading power hadronic current.
For the first two slices, since the q T value is constrained to be fairly small and thus the power corrections are significantly suppressed therein, it is challenging to observe any impacts from Λ ( q T ). However, with the value of q T growing in the next two slices, it emerges that the appreciable discrepancies between the results of N (m) LO s and those without Λ ( q T ) insertions. Instructively, it is observed that the Λ ( q T ) influences are gradually corrupted when ∆φ moves towards 180 • . For example, Λ ( q T ) can account for 3% ∼ 5% contribution in the ∆φ ∼ 150 • region, whilst the proportion decreases to less than 1% in the vicinity of ∆φ = 180 • . This phenomenon can be interpreted by the arguments above: as to the non-vanishing q T , the greater ∆φ value, the less opportunities that Λ ( q T ) is able to couple with the leptonic current. In this way, the Λ ( q T ) contributions in the right end of the graphs are relatively weaker than the left one, which therefore produces the inverse relationships between the Λ ( q T ) influences and ∆φ as illustrated in the ratio plots of the last two slices.

C Fixed-order functions
Axial Wilson coefficient C t As illustrated in (2.12), the axial-vector effective current comprises a novel structure C t O s to restore the RGI in LEEFT and encode the contributions induced by the top loops. In general, C t can be determined by matching the SM amplitudes induced bytγ µ γ 5 t −bγ µ γ 5 b onto those from O s in the limit M L m t . However, owing to differences in the O s renormalisation, the expressions for C t differ. Two prescriptions exist in the literatures: 1) Larin's [125] and 2) Chetyrkin's [123,124,126]. In the latter case, the results for C t up to four-loop accuracy have been computed [123,124,126,[181][182][183][184]. In the former scheme, the axial-anomaly form factor induced by a top quark loop has been calculated in [135] at the two-loop level. The C t expressions in different schemes can be related as where the superscripts represent the schemes of O s renormalisation. This paper employs α s which is renormalized with N F = 5 active quarks throughout. In the Chetyrkin's prescription, the finite renormalisation constant Z f, Chetyrkin 5 is actually equal to the non-singlet case Z f,Larin ns in Larin's scheme. The Z f,Larin ns expression up to O(α 3 s ) can be found in [123,185]. As to the Z f,Larin 5 expression, the first two order investigation was carried out some time ago in [125], while the third order result is given in a very recent publication [186]. In this work, we use Larin's scheme. The corresponding expression for C t reads, where L t = ln(µ 2 /m 2 t ). Here the tree-level result (−1/N F ) balances the O s /N F term in ∆ ns 3 . The singlet contribution starts from the two-loop level and the O(α 2 s ) results can be either straightforwardly read from the axial-anomaly form factor in Ref. [135], or extracted from [124,136] with the aid of eq. (C.1). We can confirm that both methods result in the same O(α 2 s ) expression. The third order expression is obtained from [123] after multiplying the converter in eq. (C.1). We have checked that the C t expression here indeed satisfies the RGE in (2.25) up to O(α 3 s ), as expected by the Larin's renormalisation scheme [125].
Non-singlet and singlet functions C ns , C V s , and C A s As shown in eq. (2.17), the hadronic currents H γ , H Z,V and H Z,A involve a set of coefficients C ns , C V s and C A s encoding the hard contributions in the loop integrals. In practice, they can be extracted from the UVrenormalized and IRC-subtracted quark form factors. To cope with the possible ambiguities arising from the γ 5 manipulation, the form factors with the Larin's prescription will be adopted throughout, in accordance with the choice in C t . In the following paragraphs, their expressions will be presented.
First, we specify the non-singlet function C ns . As illustrated in eq. (2.17), C ns participates in both the vector and axial-vector hadronic sectors. Due to the appearance of γ 5 , one may expect that the non-singlet vector quark form factor would differ from those in the axial vector case. However, since (at least) in Larin's prescription the anticommutativity of γ 5 is effectively restored for the massless QCD, the axial-vector form factor coincides with the vector one after the renormalisation [125]. Therefore, we utilize C ns to represent both cases here. Without any loss of generosity, the perturbative expansion for C ns can be defined as According to the calculations on γ * qq amplitudes, the first three coefficients can be given as [114,137,138] C (0) ns =1, where L H = ln[(−M 2 L − i )/µ 2 ]. Next, the expression for the vector singlet contribution C V s starts at the third-loop level and can be given as [114,138]  Note that as there are no -poles confronted by C V s at O(α 3 s ), the expression in eq. (C.5) contains no logarithmic terms. At last, the axial vector function C A s induced by the O s operator will be determined. Similarly, we also introduce the perturbative expansion for C A s here, where C A,(i) s encodes the O s contribution in each order. Topologically, O s can conduct both the singlet and non-singlet Feynman diagrams (see Fig. 1). Considering that the singlet part will start to work at O(α 2 s ), C A s should be identical to C ns up to NLO. Hence we have (C.7) From N 2 LO, C A s however comprises both the singlet and non-singlet contributions. To discuss, it is convenient to subtract the C ns from C A s , and define the pure singlet contribution as C A ps ≡ (C A s − C ns )/N F (see Fig. 1c). For now, the N 2 LO calculation on C A ps has been carried out in Ref. [135] based on the prescription in [125], while the C ns expressions have been specified above. After the combination, we get In comparison to C ns , it is observed that the terms proportional to L 2 H and those with the higher power remain the same, while the participation of singlet contribution has modified the others. This phenomenon is in agreement with the expectation of eq. (2.26), where the cusp anomalous dimension receives no changes but the non-cusp one endures an extra term, γ t . We also have checked that up to O(α 2 s ), C A s indeed satisfies the corresponding RGEs in eq. (2.26). In order to obtain the third order contribution, we resort to the perturbative solution of the RGE in eq. (2.26), which gives, s .

(C.9)
It is also seen that the coefficients in front of L 4 H , L 5 H as well as L 6 H all stay still with respect to C ns , while γ t make differences in the others. Also note that there is one constant term c