Mixed QCD (cid:10) QED corrections to exclusive Drell Yan production using the q T -subtraction method

: In this work we extend the q T -subtraction formalism, originally developed for QCD corrections, to the case of mixed QCD (cid:10) QED corrections, and apply it to the fully exclusive calculation of the O ( (cid:11) s (cid:11) ) contribution to the production of an o(cid:11)-shell Z boson in hadronic collisions. We present explicit results for the subtraction term and the hard factor, therefore providing all the ingredients needed for the application of the formalism up to O ( (cid:11) s (cid:11) ). To study the phenomenological impact we consider the decay of the o(cid:11)-shell Z boson into a pair of neutrinos, and present kinematical distributions for the (cid:12)nal-state leptons at LHC energies.


Introduction
The Large Hadron Collider (LHC) at CERN is reaching, for many observables, an impressive accuracy. The precision of the LHC measurements will be further enhanced in the next runs and even more in the High-Luminosity upgrade of the LHC (HL-LHC). Theory must be ready for this appointment, producing equally accurate instruments in order to interpret the high-precision data. In this framework, QCD corrections play a crucial role. However, higher order perturbative QCD corrections for many observables and benchmark processes are of the same order of QED/EW or mixed QCD-EW theoretical predictions. Next-to-next-to-leading order (NNLO) QCD theoretical predictions (the standard of theoretical precision) for observables measured at the LHC are quite accurate but, in many cases, they are not sufficient to match the current accuracy at the experimental level and the new precision that will be reached in the following years. This scenario motivates a new theoretical effort to go beyond NNLO QCD corrections by including the first QED/EW corrections, mixed QCD-EW contributions and even the next QCD perturbative order: the next-to-next-to-next-to-leading order (N 3 LO).
The Drell-Yan (DY) mechanism [1] (a benchmark process at modern colliders) constitutes a clear example of the statements of a precision observable. This process offers the possibility of studying fundamental electroweak (EW) parameters in a clean and accurate way. It also provides strong tests for QCD predictions and stringent information to determine parton distribution functions (PDFs) with high accuracy. The experimental precision for the DY mechanism at the LHC is at the percent level for the total cross section.Differential distributions and shapes are measured with an even greater accuracy (see for instance ref. [2]). The perturbative QCD corrections have been computed at nextto-leading order (NLO) in ref. [3], at NNLO for the inclusive cross section in refs. [4][5][6] and considering differential distributions in refs. [7][8][9][10][11][12]. In addition, threshold expansions have been also presented at N 3 LL accuracy in association with soft-virtual cross sections at N 3 LO in refs. [13,14]. Very recently, the N 3 LO QCD corrections have been obtained for the inclusive cross section for the production of a lepton pair via virtual photon exchange [15].

JHEP09(2020)155
However, computing several terms in the α s expansion is not enough to reach the ultimate accuracy goal, since the EW coupling α satisfies α ∼ α 2 s , and therefore NLO EW corrections, i.e. O(α), are expected to be of the same order as the NNLO QCD contributions for the total cross section. In the case of differential distributions, there are kinematical regions where the QED corrections could be as large (or even more) as the NLO QCD contributions.
The calculation of NLO EW corrections for the DY process has been addressed in refs. [16][17][18] and [19,20] for charged currents (CC) and neutral currents (NC), respectively. In order to improve our understanding of these EW effects, the calculation of their first order QCD corrections, i.e. the O(α s α), becomes necessary. These corrections represent the first term in the fixed order expansion that takes into account mixed effects from the strong and electroweak interactions. Different approaches have been followed in the literature in order to approximately combine the QCD and EW corrections [21][22][23][24][25][26][27][28][29], by either assuming the full factorisation or the additive combination of the strong and electroweak contributions. Monte Carlo event generators, which reach NLO (QCD+EW) accuracy [18,25,28] and contain all the leading logarithmically enhanced factors due to QCD and QED radiation, achieve a very good level of description of the shape of the differential distributions, but they clearly lack the NNLO QCD⊗QED accuracy for the total cross section or inclusive differential distributions as the rapidity or the invariant mass of the final state leptons.
A perturbative calculation of the Drell-Yan mechanism can be characterised by the following subsets: on one hand, purely factorisable terms that arise due to initial state (production, from the initial state partons) and final state (decay, from the final state leptons) emission and, on the other hand, non-factorisable terms originated by soft photon exchange between the production and the decay. Partial exclusive results have been presented for the resonance region, by relying on the pole approximation [30][31][32]. These non-factorisable O(α s α) terms have been shown [30][31][32] to have a negligible impact on the cross section, allowing to treat effectively Drell-Yan in the (resonant) limit of the decoupling between the production and decay processes, at least for the achieved experimental accuracy. Several steps towards the computation of the (inclusive) initial state QCD×EW corrections have been recently carried out in an analytical way [33][34][35]. The appearance of massive gauge bosons results in extra complications, so it seems natural to start by looking at the case of QED contributions instead.
The first computation of the mixed QCD⊗QED O(α s α) corrections to the inclusive on-shell production of a Z boson in hadronic collisions was achieved in ref. [36], by profiting from the available NNLO pure QCD corrections via the so-called abelianisation techniques [37,38]. Those contributions were shown to be of the order of the NNLO QCD corrections for LHC energies, which makes them relevant to reach an accurate theoretical description. Moreover, it would be highly desirable to evaluate their effect at a fully exclusive level.
A crucial ingredient in the calculation of fully differential distributions are the socalled subtraction methods. For the case of pure QCD corrections to the hadroproduction of colourless final states, the q T -subtraction method [39,40] has been extensively used in order to obtain NNLO-accurate predictions. In this work, we extend the q T -subtraction JHEP09(2020)155 formalism in order to apply it to the calculation of O(α s α) mixed corrections at a fully exclusive level. Our results are of value for transverse-momentum resummation at the corresponding logarithmic accuracy.
In particular, we will focus on the mixed QCD⊗QED corrections to the production of an off-shell Z boson decaying into a neutrino-antineutrino system. We consider the simplest case of uncharged particles in the decay of the off-shell Z boson as a way to directly address the relevance on initial state corrections to a number of exclusive observables. Note that a recent work [41] also considers the production of a Z boson, though in this case onshell, in a fully exclusive way, based on the abelianised version of the nested soft-collinear subtraction formalism [42].
This paper is organised as follows: in section 2 we present the relevant formulae for the extension of the q T -subtraction method to the QCD⊗QED case. In section 3 we present our numerical results and study the phenomenology of the corrections for different kinematical variables. Finally, in section 4 we present our conclusions.

Mixed order corrections with q T -subtraction formalism
We consider the inclusive hard-scattering reaction where the collision of the two hadrons h 1 and h 2 with momenta p 1 and p 2 produces the triggered generic final state F , without colour and electric charge, such as one or more neutral vector bosons (γ * , Z, ZZ, γγ, . . .), Higgs particles, and so forth. The observed final state F is accompanied by arbitrary and undetected radiation X, coming from the initial state only, which in this case, consists of either quarks, antiquarks, gluons or photons. The system F is composed by n final-state particles with momenta q 1 , q 2 , . . . , q n , and has total invariant mass M 2 = (q 1 + q 2 + · · · + q n ) 2 , transverse momentum q T and rapidity y. We use √ s to denote the centre-of-mass energy of the colliding hadrons, which are treated in the massless approximation (s = (p 1 + p 2 ) 2 = 2p 1 · p 2 ).
We start by considering the QCD⊗QED perturbative expansion of the (differential) cross section for the production of the final state F , by expanding in powers of the strong (α s ) and electromagnetic (α) couplings, where dσ  . Following a similar structure to the one valid in the pure QCD case [39,40], the basic formula for the q T -subtraction method in the case of mixed QCD⊗QED corrections can be expressed in the following way, 3)

JHEP09(2020)155
where dσ (1,1) F +jet corresponds to the F +jet production cross section at O(α s α). It is important to note that in this context 'jet' stands for either quarks, antiquarks, gluons or photons in the final state and all of them need to be considered in the initial state as well. The term inside the square bracket in eq. (2.3) is finite in the limit of vanishing transverse momentum of the F state, but the individual terms dσ (1,1) F +jet and dσ (1,1) F CT are separately divergent. In order to evaluate dσ (1,1) F +jet , we can make use of any NLO subtraction method (adapted, though, to the case of mixed QCD⊗QED corrections).
The subtraction counter-term dσ (1,1) F CT encodes the singular behaviour of the real scattering amplitudes in the small-q T region. The coefficient function H (1,1) F restores the correct normalisation to the total cross section and it has Born kinematics (e.g. it is proportional to δ(q T )). Both coefficient functions can be obtained, through the abelianisation procedures [37,38], from eqs. (63-70) in ref. [43]. We have checked (as a self-consistency check) that the same coefficient functions can be obtained from first principles, i.e. redefining eq. (6) in ref. [40] to take into account QED emissions and expanding it to a given fixed order.
We present in the following the explicit expression of all the required terms needed for the subtraction at O(α s α). These are constructed by convoluting the parton distributions with the corresponding partonic terms, which up to O(α s α) are given by In order to simplify the notation, we indicate by z the dependence on both partonic momentum fractions z 1 and z 2 . The explicit dependence on either z 1 and z 2 can be easily understood in terms on the dependence on the partonic label a and b, respectively. Also, it is implicit the dependence on the renormalisation (µ R ), factorisation (µ F ) and resummation (Q) scales. Note that, for the sake of generality, in the results contained in this section we keep the full dependence on the resummation scale [43]. This dependence is needed in the context of transverse-momentum resummation. The fixed-order cross-section is independent of this scale, and it is convenient to set Q = M to simplify the corresponding expressions.

JHEP09(2020)155
The contributions to the counter-term Σ F (i,j) cc←ab can be organized in the following way according to their power of logarithmic enhancement. The dependence on the transverse momentum is given by the known integrals [43] where b is the impact parameter, J 0 (x) is the 0th-order Bessel function and b 0 = 2e −γ E , with γ E representing the Euler number. Notice that we are using the "+1" prescription (see the argument of the logarithm inside eq. (2.9)), and therefore, the counter-terms vanish in the large-q T limit. More details about eq. (2.9) can be found in the appendix A of ref. [43].
The corresponding coefficients for the expansion of Σ cc←ab are more easily presented by considering their N -moments (Mellin) with respect to the variable z. At NLO in QCD and QED they are given by

JHEP09(2020)155
while for the mixed QCD⊗QED corrections at O(α s α) they are given by In the expressions above we have defined Q = ln M 2 /Q 2 and F = ln M 2 /µ 2 F , while γ

JHEP09(2020)155
and their explicit expression for quark-initiated case is given by

23)
Notice that we consider the electromagnetic coupling α as constant, in the sense that it is not running with any of the scales related to the process. For that reason the QED beta-function does not appear in the coefficients of eqs. (2.10)-(2.20). Eqs. (2.10)-(2.15) were derived for first time in ref. [44], where the transverse-momentum resummation for Z boson production combining QED and QCD was computed at NLO, considering the corresponding logarithmic accuracy: next-to-leading logarithmic (NLL). Although results in ref. [44] were only applied to the transverse momentum distribution of the Z boson, the formalism is fully differential and can be used to compute any infra-red safe exclusive distribution. It is worth noticing that for transverse-momentum resummation some novel mixed effects appear affecting the distribution already at leading logarithmic (LL) accuracy (see eqs. (7) and (11) of ref. [44]). Nevertheless, that contribution can only show up after performing the fixed order expansion up to O(α s α) (see eq. (3) in ref. [44]) if the electromagnetic coupling α is considered to be running, which is not the case in our current study. In this paper we confirm the results presented in ref. [44] and we extend the fixed order computation to the next perturbative order, the NNLO QCD⊗QED. Finally we present the collinear functions, again for c = q, and the hard-virtual coefficients, the latter specifically for the DY case as they are a process-dependent quantity. The separation between C and H coefficients is scheme dependent. Those presented here are obtained in the so-called hard scheme [40]. Up to NLO in QCD and QED, the hard-virtual coefficients take the form 24) and the collinear functions are given by The hard-virtual coefficient needed for the first order in the mixed QCD⊗QED expansion takes the following form,

JHEP09(2020)155
while the needed collinear functions are given by the following expressions: The results above provide all the ingredients needed for the application of the q Tsubtraction formalism to the calculation of mixed QCD⊗QED corrections. The same coefficients are required by the transverse-momentum resummation formalism, considering in this case the full dependence on the resummation scale Q.
In the following section we present our phenomenological results for the case of Z boson production.  α) contribution to the inclusive on-shell Z production cross section for the different partonic channels. The results obtained using q T -subtraction are compared to the inclusive predictions obtained in ref. [36]. Numerical uncertainties on the last digit are indicated in parenthesis for our predictions, while the uncertainties of the inclusive implementation are below the last digits shown. The category denoted by qq includes all combinations of quarks and anti-quarks.

Phenomenological results
In order to obtain quantitative results, our calculation is implemented in two independent parton-level Monte Carlo programs. One of them is based on MCFM [12] (including the NNLO QCD corrections), suitably modified to deal with mixed corrections and to apply the q T -subtraction formalism. The other is a private implementation, which relies on the FKS subtraction method [45] to deal with the NLO-type divergencies (adapted to the mixed QCD⊗QED case), and on analytic results for the relevant scattering amplitudes obtained from ref. [46], plus an explicit calculation of the tree-level all-quarks channels using FeynCalc 9.2.0 [47].
For our phenomenological analysis we consider n F = 5 massless quark flavours. We work in the G µ scheme for the EW couplings, using the input values G µ = 1.16639 × 10 −5 GeV −2 , M Z = 91.1876 GeV and M W = 80.385 GeV. The width of the Z boson is set to the value Γ Z = 2.4952 GeV. For the parton luminosities and strong coupling, we use the NNPDF3.1luxQED set with five flavours [48] through the LHAPDF interface [49], always at NNLO accuracy, regardless the order of the calculation. Both renormalisation and factorisation scales are set to the default value µ R = µ F = m 1 2 . For the cutoff parameter of the subtraction method, q T,cut , we choose the central value q T,cut = 0.2 GeV. We checked that our results are compatible within uncertainties when varying this parameter around its central value by a factor of 2.
As a first check of our implementation, we computed the inclusive cross section for the production of an on-shell Z boson, and compared to the predictions obtained from the analytic results presented in ref. [36]. The corresponding O(α s α) contributions to the cross section, split into quark-quark, quark-gluon, quark-photon and gluon-photon initiated channels, are shown in table 1. As can be seen from the table, we can reach sub-percent precision for these inclusive predictions, and we find full agreement with the analytic results from ref. [36]. As an additional validation, we have computed the NNLO QCD differential distributions using the public code Matrix [50], finding full agreement with our results.
For all of the differential distributions presented here, we consider the following set of cuts, p T, 1 > 25 GeV , p T, 2 > 20 GeV , |y| 1,2 < 2.5 , m 1 2 > 50 GeV, where 1 and 2 represent the final-state hard and soft leptons respectively, ordered according to their transverse momentum (p T, 1 > p T, 2 ). Since we consider only neutrinos in the final state, there is no need to recombine collinear leptons and photons. We start by presenting the transverse momentum distribution of the leptons in figure 1. The kinematical dependence of the mixed corrections is highly non trivial. This feature is also shared by the pure QCD and QED corrections, and it is expected due to the particular features that these two distributions present at fixed order in perturbation theory. At LO both leptons are back-to-back, and therefore the distributions are identical. The radiative corrections produce the change of shape that render the p T, 1 spectrum harder than the p T, 2 one, producing therefore sizeable distortions in the distribution. Furthermore, some regions of the phase space are almost not populated at LO, and therefore radiative corrections become more relevant. This is the case for the region of p T, 2 below the lower cut on p T, 1 , which is directly not allowed for Born kinematics, or the region above p T, 1,2 M Z /2, which does not receive contributions from the Z peak at LO.
From figure 1 we can observe that for p T, 1 < M Z /2 the mixed corrections are positive, representing an increase of about 0.5% with respect to the NLO prediction. The corrections then change sign, being of the order of −0.5% in the first bins after p T, 1 = M Z /2, which corresponds to the expected Sudakov shoulder near the kinematic boundaries mentioned in the previous paragraph [51]. The mixed corrections then result smaller at the tail of the distribution. With respect to the softer lepton, we can observe that the corrections become very large around and slightly above p T, 2 = M Z /2, a pattern shared by the NNLO QCD corrections. In this region, the effect of the mixed QCD⊗QED contribution can reach the dσ (2,0) dσ (1,1) × 50 pp → Z * → νν @ 13 TeV O(5%) with respect to the NLO QCD result. In addition, we can also observe a small (negative) peak in the corrections around p T, 2 = 25 GeV, which is related to the presence of a cut in p T, 1 , as mentioned before. Since the kinematic region around p T, 1 50 GeV (p T, 2 20 GeV) is affected by integrable singularities of soft origin, it is possible to mimic the effect of their corresponding all-order resummation by enlarging the bin size, recovering the reliability of the computation around the peak.
We continue by presenting the rapidity distributions of the leptons, again ordered according to their transverse momentum, in figure 2. In both cases, we can observe that the mixed corrections are extremely small, and show a very mild dependence on the corresponding kinematical variable. The reason for this particularly small value of the corrections is a very strong cancellation between the main partonic channels, that is the qq and qg initiated processes, over the whole rapidity range under consideration, a pattern that can also be observed for instance at the level of the total cross section. We note that this effect is even stronger with the set of cuts in eq. (3.1), compared to the fully inclusive case, with cancellations of about 90% between the different channels.
In figure 3 we present distributions for the lepton-pair system, specifically its transverse momentum and rapidity. The mixed corrections are negative below p T, 1 2 ∼ 15 GeV, and diverge in the p T, 1 2 → 0 limit. The sign of the mixed corrections in the low transverse momentum region is the same as the one of the NNLO QCD corrections, as one can infer from the sign of the logarithmic coefficient with highest power (see eq. (2.16) for the mixed corrections and eq. (66) of ref. [43] for NNLO QCD). Above p T, 1   region the NLO QED corrections are of the order of 0.5%. As it is well known, at low-q T , the large logarithmic corrections to the cross section have to be treated with transverse momentum resummation in order to recover the reliability of the prediction. This is true not only for the transverse momentum distribution but for any observable which presents a kinematical region directly related to q T = 0.
The mixed corrections for the lepton-pair rapidity present a kinematic dependence that is similar to the one of the NNLO QCD contribution. They are negative for small |y| 1 2 , and become positive for larger values of rapidity. The overall size of the mixed corrections is of course much smaller though, being of the order of 50 times smaller than the NNLO QCD corrections.
Finally, we present in figure 4 the φ * and cos θ * distributions, defined as [52] φ * = tan π − ∆Φ 2 sin θ * Since at LO the two leptons are back-to-back, the φ * distribution is trivial at that order, and contributions with φ * = 0 only start at NLO. As in the case of the transverse momentum, the small-φ * region is not well behaved at fixed order and it is necessary the use of transverse resummation in order to recover the reliability of the prediction in those kinematical regions. The pattern of corrections, not only for the mixed but also for JHEP09(2020)155 dσ (2,0) dσ (1,1) × 50 pp → Z * → νν @ 13 TeV the NNLO QCD and NLO QED contributions, is very similar to the one observed in the p T, 1 2 distribution, in particular with the mixed corrections being negative at small φ * and becoming positive for larger values, and about a factor of 2 smaller than the NLO QED corrections in the tail of the distribution.
In the case of cos θ * , the distribution is rather flat in the central region, and presents a strong suppression for cos θ * = ±1, which is only populated by events with very large and opposite rapidities of the corresponding leptons. This region is therefore particularly suppressed by the presence of the cuts on y 1,2 , which directly forbid the region above | cos θ * | ∼ 0.987. From the lower panel of the figure we can observe that the perturbative corrections are rather flat in the region where the bulk of the cross section is located, and therefore they follow a pattern similar to the one observed for the total cross section. In particular, the mixed QCD⊗QED corrections are extremely small, and become more relevant only close to the boundaries, where they reach the 0.6% level (note that the last bin of the distribution is larger and extends from | cos θ * | = 0.8 to 1).
Before going to the summary, it is interesting to compare the size of the mixed QCD⊗QED corrections computed here against the naive approximation in which QCD and QED corrections factorize. Specifically, defining for a given bin dσ (1,1) dσ (1,1) approx pp → Z * → νν @ 13 TeV Figure 5. Comparison between the mixed QCD⊗QED corrections (red) and the naive factorisation approximation (purple), for the transverse momentum of the hardest (left) and softer (center) lepton, and the rapidity of the pair (right).
The naive approximation shown in eq. (3.4) was implemented in several studies using fixed order numerical tools. Notice that experimental analyses, which use Monte Carlo event generators [18,25,28] rely on a combination of a QCD generator, convoluted at fully differential level, with QED parton showers, thus taking into account kinematic effects that the proposed naive factorisation misses.
In figure 5 we present the mixed QCD⊗QED corrections together with the approximation defined by eq. (3.4), for the transverse momentum of the two leptons and the rapidity of the pair. The results are normalized to the NLO QCD prediction, as in the lower panels of the previous figures. We can observe that, in all cases, the multiplicative approach is a rather poor approximation to the full results. This is in line with the observations made for the total cross section in ref. [36]. The discrepancies, however, can be strongly enhanced at the differential level. This can be seen for instance in the p T, 1 > M Z /2 region, where the exact O(α s α) corrections are at the per-mille level, while the factorisation approximation predicts ∼ 7% corrections. The reason for this big discrepancy is the presence of large K-factors at NLO (both in QCD and QED), associated to the fact that at LO this region is only populated by events that are away from the Z peak. In the case of the p T, 2 distribution, we can observe that the multiplicative approach has the wrong sign for p T, 2 < M Z /2 (note that the approximation is not well defined for p T, 2 < 25 GeV due to the cut in the hardest lepton), and fails to reproduce the correct size of the corrections around the peak located in p T, 2 ∼ M Z /2. Finally, for the rapidity of the lepton pair we can see that the factorisation approximation predicts a rather flat K-factor, failing to describe the kinematical dependence of the mixed corrections.

JHEP09(2020)155
The method can be applied to the fully exclusive calculation of the O(α s α) corrections for the production of a colourless and neutral final state (e.g. Z and Higgs bosons, photons, neutrinos). We have provided all the relevant formulas for its implementation at O(α s α). The coefficient functions and the hard virtual coefficients are also of value for transverse momentum resummation and our expressions contain the full dependence on the resummation scale Q.
We have applied the method to the production of an off-shell Z boson, and considered its decay into a pair of neutrinos. We presented differential distributions for the final-state leptons at the LHC, and found that the corrections can have a sizeable dependence on the kinematics, and not necessarily following the pattern of the NNLO QCD corrections for instance. The size of the corrections is typically very small and below 1%, though it can be enhanced in some particular phase space regions. We note that our predictions are in qualitative agreement with the corresponding results in ref. [41].
We have also compared the mixed QCD⊗QED contribution with the factorisation approximation based on the product of QCD and QED K-factors. We have found that this multiplicative approach is in general a bad approximation to the mixed corrections, and the disagreement can be quite extreme for some differential distributions.
As a final remark, it is interesting to point out that recent developments have allowed the application of the q T -subtraction method to the production of a heavy-quark pair at NNLO in QCD [53][54][55] (see also its related application to NLO EW corrections for massive lepton pair production in ref. [56]). Following similar abelianisation techniques to the ones used in the present paper, the method could be extended to also deal with the mixed QCD⊗QED corrections for the production of a massive charged (colourless) final state.