QCD$\oplus$QED NNLO corrections to Drell Yan production

We compute the QCD$\times$QED (${\cal{O}}(\alpha_s \alpha)$) mixed and QED$^2$ (${\cal{O}}(\alpha^2)$) corrections to the production of an on-shell $Z$ boson in hadronic collisions. We obtain them by profiting from the calculation of the pure QCD terms after taking the corresponding abelian limits. Therefore, we extend the available knowledge up to complete next-to-next-to leading order precision in QCD$\oplus$QED. We present explicit results for the perturbative coefficients and perform the phenomenological analysis at different collider energies with particular emphasis on the mixed corrections. We study the contribution from the different channels and discuss the scale dependence stabilisation effect. We consider a factorisation approximation for the mixed order terms and show that it fails to reproduce the exact result. We find that the contributions are small, typically at the few per mille level, but that under some kinematical conditions they can compete with the pure QCD NNLO ones.


Introduction
In recent years the development of high precision experiments in particle physics demanded a theoretical upgrade to match the accuracy achieved at the LHC. The state-of-the-art in fixed order computations for processes with up to two hard partons in the final state is reaching nextto-next-to-leading-order (NNLO), i.e. O(α 2 s ). Since α 2 s ∼ α, it becomes necessary to include also the corresponding NLO ElectroWeak (EW) corrections, that for many observables exceed the few percent level (e.g. [1,2]) and become quantitatively important for an accurate description. Both precise measurements and calculations are essential to test different aspects of the Standard Model (SM) and to discern between them and possible new physics evidences due to Beyond the Standard Model (BSM) effects.
In this sense, inclusive massive lepton pair production (Drell-Yan process) has worked as an important testground of perturbative quantum chromodynamics (QCD). On one hand, it has offered a sensitive way to study parton distribution functions (PDFs) [3][4][5]. From weak bosons production, charge asymmetry measurements and invariant mass dependences have helped to extract precise information on both the quarks valence structure functions and the separation of quarks flavours, as well. On the other hand, knowing the behaviour of charged-current (CC) and neutral-current (NC) processes has allowed to perform high-precision measurements of fundamental ElectroWeak parameters, like Z and W widths and masses and the EW mixing angle. With W and Z bosons in final or intermediate states of hadronic processes, some clean production channels can be characterised with great level of accuracy by studying leptonic decay modes. This has also served to accurately calibrate detector components, which was important for further measurements [6].
In addition, the Drell-Yan process is not only relevant to test SM predictions, but also to evaluate alternative BSM theories, where W and Z bosons usually appear as final or intermediate states in the decay of particles predicted in new physics models, like new gauge interactions, supersymmetry or heavy resonances [7]. For all these reasons, the study of Drell-Yan processes, and in general, theoretical and experimental exploration of multiple EW gauge boson production, becomes a priority in the development of current high-energy physics.
In this sense, and considering that the improvement in statistics over the last years has made higher-order corrections experimentally noticeable, having access to QCD (and QED) corrections to these processes has become of great importance to put the previous predictions on a firmer ground. For instance, it is important to achieve a better comprehension of the so-called Sudakov regime, where large logarithmic EW factors play an important role, and O(αα s ), that represent the first QCD corrections to these large EW factors, or even O(α 2 ) corrections, might be sizable.
Furthermore, recent work has been performed to include QED effects in the evolution of parton distributions, by providing explicit expressions for splitting kernels up to O(αα s ) [8] and O(α 2 ) [9] and by the determination of precise photon distributions in the proton within the LUXqed approach [10,11] . The availability of these QED corrected parton distributions is essential to match the theoretical calculations at the partonic level.
From the point of view of partonic cross-sections, the O(αα s ) and O(α) corrections represent the first EW and mixed order contributions to Drell-Yan pair production in the general expansion where pure EW dσ (0,j) and QCD dσ (i,0) corrections, as well as mixed order contributions, which combine effects of the two interactions, arise.
On the other hand, concerning the EW contributions, exclusive computations for NLO-EW corrections to CC-DY are available in refs. [24][25][26] and for NC-DY, in refs. [27,28]. Finally, progress towards the computation of NNLO-EW has been accomplished in recent years too [29][30][31][32]. Due to the lack of the full calculation of the NNLO mixed-order terms O(αα s ), different approaches have been followed to approximately combine the QCD and QED/EW corrections [33][34][35][36][37], by either assuming the full factorisation or the additive combination of the strong and electroweak contributions. Particularly, recent partial exclusive results have been presented for the resonance region, by using the pole approximation [38][39][40].
The contributions for a general (i.e. including the decay of the gauge boson) perturbative calculation of Drell-Yan can be roughly characterised into 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. The non-factorisable O(αα s ) terms have been shown [38][39][40] 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. The results presented in [40] also rely on the assumption that the missing initial-initial state factorisable O(αα s ) contributions are very small.
The computation of the so far unknown mixed QCD×QED O(αα s ) corrections to the inclusive on-shell production of a Z boson in hadronic collisions is exactly the main goal of this paper † . Those contributions are by themselves a gauge-invariant set of the complete Drell-Yan cross section calculation at O(αα s ), even for CC-DY. Furthermore, counting with analytical expressions for the total cross section can be useful to establish a subtraction method to compute differential distributions for different observables at O(αα s ) by extending, for example, the q T − subtraction method [41] originally developed for pure QCD corrections.
In principle a full computation of QCD×QED O(αα s ) terms involves, as in any NNLO calculation, the evaluation of double-virtual, single-virtual plus one parton emission and double parton emission contributions, where parton in general refers to quarks, antiquarks, gluons, and photons.
Most of the needed double real contributions where recently presented in [42], including the case of W boson production which is not discussed in this paper ‡ ,while the master integrals for the two-loop calculation were obtained in [43].
Instead of following the path of a dedicated calculation for each term, in this work we profit from the available computation of NNLO pure QCD corrections O(α 2 s ) presented in [13] and, by pointing out the abelian component, we extract the corresponding QCD×QED O(αα s ) contributions and asses the phenomenological impact for the inclusive cross section at different hadronic energies.
Furthermore, by following the same procedure we also present the QED 2 O(α 2 ) corrections, completing, therefore, the set of NNLO contributions in QCD⊕QED (i.e. all terms that correspond to i + j = 2 in Eq.(1)). This paper is organised as follows: In Section 2 we present the method used to compute the QCD×QED O(αα s ) and QED 2 O(α 2 ) contributions from the pure QCD corrections O(α 2 s ). In Section 3 we present the results and study the detailed phenomenology of the corrections at different collider energies. Finally, in Section 4 we present our conclusions. The explicit NNLO coefficients are presented in the Appendices.

Abelianisation procedure
In order to achieve the O(αα s ) contributions we profit from the availability of previous NNLO QCD calculations [13]. In summary, we analyse the contributing diagrams for each interaction and take the corresponding abelian limit from the existing NNLO QCD calculation [8,9]. To schematise the procedure used to accomplish the mixed order contribution, we describe the algorithm of gluon-photon interchange taking as an example the most relevant qq channel.
First of all we analyse all the topologies inherent to each partonic subprocess in QCD and determine the corresponding colour factors. Then we replace a gluon by a photon in each diagram and profit from the similarity in QED and QCD fermion coupling factors to pin down the abelian parts of both calculations. Furthermore, given that the photon has no colour, non-abelian terms (i.e. all the terms arising from NNLO QCD diagrams that contain three or four gluon vertices) do not contribute to QED corrections and therefore can be thrown out. Thus, we recalculate the colour factors for the mixed topologies attained identifying changes and substitutions to be made in previous QCD results in order to obtain mixed order QCD×QED correction terms.
For the explicit case of the qq channel, in order to show how the method works and given that the colour factors in the partonic cross section only depend on the QCD structure, we concentrate on diagrams with double real emission which appear only to tree level (the same considerations can be applied to the two-loop and one-loop plus real emission contributions). As can be observed in Fig.1 there are four kinds of diagrams to consider, which we will label with the supra-index (k, l) according to the total number of QCD (k) and QED (l) vertices for each topology. It is also important to note that the order of external momenta does affect the colour structure of the diagram. In this sense, we will refer as (a) to the diagram showed in Fig.1, while (a ′ ) represents the corresponding diagram obtained after crossing the final state parton lines. There we can recognise ‡ Only the contribution from the interference of QCD and QED qq → qq diagrams is missing in [42]. For example, terms corresponding to the calculation of (a) (2,0) 2 in NNLO QCD result proportional to 1 in the denominator arises due to the average over the colour factor of the incoming quarks and the symmetry factor 1/2 is due to the appearance of two identical gluons in the final state. For the case of (a) (2,0) × (a ′ * ) (2,0) both abelian and non-abelian contributions appear resulting in a factor 1 ) and, when considering terms from (b) (2,0) × (a * ) (2,0) , they result proportional to 1 , a purely non-abelian contribution.
Once colour factors are characterised for each term, we choose a gluon in the diagram, replace it by a photon and recalculate the colour structure, thus obtaining modified diagrams with the corresponding new factors for QCD×QED corrections. These are shown in Fig.2. Naturally, all the diagrams of type (b) (i.e. topologies containing at least one 3-gluon-vertex), which always contribute to second order for the NNLO QCD calculation in this process, vanish when considering the abelian limit.
Taking this into account, we find that the modified factors for (a) (1,1) 2 and (a) (1,1) × (a ′ * ) (1,1) are both given by where we have included the charge of the quark for the QED coupling, while the non-abelian one obviously vanishes. Here we may notice that all the colour factors proportional to C A , which corresponds to non-abelian part of the calculation, could be thrown out when considering the abelian limit, while the ones proportional to C 2 F are to be replaced by 2e 2 q C F , thus obtaining QCD×QED factors in each case. It is worth noticing by performing the same analysis for the topology shown in Fig.1c, i.e. the production of a qq pair, (1,1) Figure 2: Diagrams that result after applying the abelianisation procedure to the real NNLO QCD corrections in Fig.1.
that also the colour factor T R vanishes when the similar contribution is analysed in the QCD×QED case, since the result for (c) ( . Therefore, since terms proportional to both C A and T R are vanishing, the same occurs for terms proportional to β QCD 0 in the original pure QCD calculation, consistent with the fact that no renormalisation is needed at this order either for the QED or QCD couplings § . Same wise, only a few contributions survive in the products of the type (c) (2,0) × (d * ) (0,2) and (d) (2,0) × (d ′ * ) (0,2) , i.e. the interference of amplitudes with one photon and with one gluon exchange.
This strategy can be extended for all the topologies in qq. In Table 2 we show the different colour factors (after factorising an overall factor of 1/2N C ) for diagrams contributing to σ (2,0) , and the resulting ones after the abelianisation procedure corresponding to σ (1,1) . The replacements in the colour structures needed to go from the NNLO QCD coefficients to the QCD×QED ones can be directly read from the entries in Table 2.
As an important feature, this method shows to be versatile in order to obtain NNLO QED corrections to Drell-Yan as well (i.e. the calculation of σ (0,2) ), if a deeper abelian limit is considered in this case. Here, by turning two gluons into photons from the topologies of NNLO QCD calculation one can recover correction terms up to second order in α, thus completing the set of QCD⊕QED NNLO corrections to Drell-Yan, in the sense of Eq.(1). The corresponding colour factors (including electric charges of both quarks and leptons that might appear in the final state) are also shown in Table 2 for the qq channel.
The same occurs for other channels, after treating carefully the initial flux factor, which depends on the colour properties of initial state particles. For instance, both qγ and qg contributions § As stated above, we consider the Born coupling between the quarks and the Z in the sense of an effective coupling. Table 1: Colour factors corresponding to qq channel for each contribution to NNLO QCD⊕QED corrections to Drell-Yan, up to an overall 1 2N C factor. Focusing on α 2 factors, the third column includes sums over sets of quark (Q) and lepton (L) final state charges, while e i and e j refer to different quark flavour charges in the scattering.

Colour factors in qq
to σ (1,1) can be obtained from the qg calculation for NNLO QCD corrections, by choosing the initial or final state gluon, respectively, to perform the abelianisation and following the procedure detailed above. Particularly, in the case of γg channel, we have performed the explicit calculation of the fixed order corrections, finding perfect agreement with the result obtained by applying the abelianisation procedure.

Results and Phenomenology
In general the cross section can be written as where σ Z is the point-like LO cross section, √ S is the hadronic centre-of-mass energy, Q the invariant mass of the produced Z, τ = Q 2 S and W Z (τ, Q 2 ) is the hadronic structure function. The point-like cross section that appears in Eq.(2) is defined as where N C = 3 is the number of quark colours, θ W is the weak mixing angle (with sin 2 θ W = 0.23), M Z = 91.187 GeV and Γ Z are the mass and width of the Z, and Γ Z→X is the partial width due to the decay of the Z to X (e.g. for leptonic decay, X = ℓl). The narrow-width approximation used along this paper consists on making the following replacement ensuring the decoupling of the production and decay mechanisms. The hadronic structure function appearing in (2) can be written as a sum of contributions of different orders where the dependence on the factorisation µ F and renormalisation µ R scales is implicit.
The analytic expressions for the inclusive cross section of Drell-Yan Z-production at QCD⊕QED NNLO are presented in the Appendices. In this section we study the phenomenology of the total inclusive cross section, i.e. in all the decay channels of the Z, within the narrow-width approximation. To this end, a specific code was written which makes use of the LHAPDF [44] package to interpolate sets of parton distribution functions.  to NNLO (QCD) accuracy [3][4][5]45] with the corresponding QED corrections from LUXqed [10,11]. In Fig.3 we plot the K-factors for different orders as a way to quantify the size of the QED and QCD corrections to Drell-Yan at different centre-of-mass energies.
Here the K-factor is defined as the ratio of the cross-section computed at a given order over the previous one, i.e.
As can be observed, the NNLO QCD corrections are of the same (∼ 5 per mille level) order, but typically with the opposite sign, as the NLO QED corrections, as expected from the simple counting α 2 s ∼ α. The mixed QCD×QED turn out to be positive and below the per mille level over the whole range of energies spanned in the plot. Interestingly, due to the particular dependence of the NNLO QCD corrections with the energy, with a sign change around √ S ∼ 18 TeV, for the LHC at √ S ∼ 14 TeV the mixed QCD×QED corrections are only a factor of ∼ 3.5 smaller than the pure NNLO QCD contributions. Furthermore, for lower centre-of-mass energies √ S ∼ 2 TeV the mixed terms almost reach the per mille level and are just a factor of 5 smaller than the NLO QED ones, showing that the elementary counting of couplings can fail under certain kinematical conditions. The pure NNLO QED terms, also plotted in Fig.3, are negative but the corrections always remain at the O(10 −5 ) level.
Even though for this particular observable the mixed QCD×QED contributions are small, it is interesting to study how well they can be approximated by the factorisation assumption on QED plus QCD corrections, where it is assumed that σ (0,0) σ (0,0) , compared to the exact case κ mixed = αα s σ (1,1) σ (0,0) . For that purpose, in Fig.4 we plot the following quantity which is the ratio between the exact and the approximated factorised contribution. As it can be observed, the factorisation approach fails to reproduce the correct behaviour of the mixed contribution typically by a factor of two or more. Of course, given the size of the corrections, the effect of the factorised treatment of these contributions is small at the level of the cross section, as shown in the inset plot of Fig.4, where we show the ratio between the cross section computed exactly and within the factorisation approach, but the situation might not hold for other observables or even for more exclusive distributions in Drell-Yan. Figure 5: Contribution to the mixed QCD×QED K-factor from the different channels. Here the label q accounts for both quarks and antiquarks and qq represents the sum of qq and qq.
In Fig.5 we show the contribution to the mixed QCD×QED K-factor from the different channels. It is noticeable that the photon initiated contributions are rather small, mostly due to the size of the photon pdf in the proton, as can be observed by comparing qγ and qg contributions, which share the same partonic coefficient apart from the colour factor. It is also clear that the different signs of qq (fully dominated by the born level qq channel and exceeding the per mille level) and qg contributions conspire to reduce the effect of the mixed QCD×QED corrections to the Drell-Yan cross section. Again, in more exclusive distributions this partial cancellation might be spoiled by some kinematical cuts, resulting in an increase of the mixed order corrections. Finally, we discuss the effect of the higher order contributions in the stabilisation of the perturbative expansion in terms of the scale dependence for √ S = 13 TeV (very similar behaviours are observed for other values of √ S). In Fig.6 we show the LO (σ (0,0) ), NLO (σ (0,0) +α σ (0,1) +α s σ (1,0) ) and NNLO (σ (0,0) + α σ (0,1) + α s σ (1,0) + αα s σ (1,1) + α 2 σ (0,2) + α 2 s σ (2,0) ) cross sections for different values of the factorisation and renormalisation scales µ R = µ F = µ, normalised by the corresponding value at the central scale µ = M Z . From the slope of the different curves, it is clearly visible the reduction in the scale dependence when including higher order corrections, mostly due to the dominant QCD effects but also thanks to the inclusion of the QED and mixed contributions.

Conclusions
In this article, mixed QCD×QED as well as pure QED 2 NNLO corrections to the total Drell-Yan Zproduction cross section were presented for the first time. This was achieved via an abelianisation procedure that profits from the available pure QCD NNLO result and proved to be a versatile technique.
We performed the phenomenological analysis finding that the mixed corrections are of the order of per mille at the LHC, but only a factor of ∼ 3.5 smaller than the pure QCD NNLO due to a sign change that occurs in the latter at √ S ∼ 14 TeV. Pure QED NNLO terms are shown to be negative corrections of the order of 10 −5 . The full QCD⊕QED NNLO corrections are found to further stabilise the scale dependence of the final result.
The exact K-factor at order O(αα s ) was compared with the naive factorisation approximation, which consists of the mixed order term of the product of QCD and QED NLO K-factors. It was shown that the latter fails to reproduce the exact result by a factor of two or more. Although in this case the difference is not significant, due to the smallness of the overall contribution, this result hints that the factorisation approximation may not work for other processes nor for more exclusive measurements of the one presented herein.

Appendix A: Definitions and Notation
In this appendix we present for the first time, the coefficients w (1,1) and w (0,2) in the sense of Eq. (5). They include the NNLO mixed QED×QCD and QED corrections, respectively.
where Q (Q) and L are the sets of (anti)quarks and leptons considered (e.g. Q = {u, d, c, s, b}, L = {e, µ, τ }), e k is the charge of the particle k in units of the electron charge, µ F and µ R are the factorisation and renormalisation scales, c i = v 2 i + a 2 i with v i and a i defined as the vector and axial couplings of particle i: and replicated through families, the QED and QCD beta functions are given by n F is the number of quark flavours considered (n F = #Q), and the various ∆(x) corrections functions are defined in Appendix B.
The other coefficients needed for the full NNLO calculation, w (i,0) for i ≤ 2 and w (0,1) , were presented in [13] and [27,28] respectively. Here they are rewritten in this notation for the sake completeness and as a reference for the full set of NNLO corrections to the Drell-Yan Z production.

Appendix B: Correction Terms
In order to present the expressions for the different corrections, we define the following distributions which appear in the soft terms regularising the divergence of soft emission (x ≈ 1) and defined as usual by We also define an auxiliary variable to write the dependence on the factorisation scale, where µ F is the factorisation scale and Q the invariant mass of the produced Z.
The corrections needed for the NLO result are For the second NNLO, several correction terms are introduced. We denote with a C A superscript the corrections coming from the non-abelian part of the contributions (only relevant for the QCD NNLO contribution), with n F the ones that involve a sum over fermion families (relevant for the QCD and QED NNLO result), and with C F the rest of the abelian contributions.
The non-abelian contributions on the quark-antiquark channel are Here we amend the result for ∆ (2) qq,A 2 given in Eq.(B.11) of ref [14] by adding the corresponding missing x factor to term 103x above.
The singlet contributions for the (anti)quark-(anti)quark channel include terms arising from identical initial quarks, and non-identical ones. These are given by the following expressions.