Mixed NNLO QCD x electroweak corrections of O(N_f \alpha_s \alpha) to single-W/Z production at the LHC

First results on the radiative corrections of order O(N_f\alpha_s\alpha) are presented for the off-shell production of W or Z bosons at the LHC, where N_f is the number of fermion flavours. These corrections comprise all diagrams at O(\alpha_s\alpha) with closed fermion loops, form a gauge-invariant part of the next-to-next-to-leading-order corrections of mixed QCD x electroweak type, and are the ones that concern the issue of mass renormalization of the W and Z resonances. The occurring irreducible two-loop diagrams, which involve only self-energy insertions, are calculated with current standard techniques, and explicit analytical results on the electroweak gauge-boson self-energies at O(\alpha_s\alpha) are given. Moreover, the generalization of the complex-mass scheme for a gauge-invariant treatment of the W/Z resonances is described for the order O(\alpha_s\alpha). While the corrections, which are implemented in the Monte Carlo program RADY, are negligible for observables that are dominated by resonant W/Z bosons, they affect invariant-mass distributions at the level of up to 2% for invariant masses of>500 GeV and are, thus, phenomenologically relevant. The impact on transverse-momentum distributions is similar, taking into account that leading-order predictions to those distributions underestimate the spectrum.


Introduction
The production of charged leptons via an electroweak (EW) gauge boson in hadronic collisions, known as Drell-Yan-like W/Z production, is among the most important processes at the LHC [ 1,2,3,4] owing to its clean experimental signature and high cross section.Both luminosity monitoring and detector calibration are possible using Drell-Yan-like processes, the former by using the total cross section and the latter by performing measurements of the mass and width of the Z boson.On the theoretical side, the Drell-Yan (DY) production of lepton pairs is among the best understood processes, and in combination with the distinct experimental signature it is possible to use them to constrain parton distribution functions (PDFs) [ 5] via the W charge asymmetry and the Z rapidity distribution.Furthermore, DY production can be used to measure EW precision observables such as the W-boson mass [ 6] or the effective weak mixing angle sin 2 θ lept eff [ 7].
A natural next step is the calculation of mixed QCD×EW NNLO O(α s α) corrections which are assumed to be the largest unknown fixed-order part.Given the complexity of the full calculation, several approximations were applied to get a handle on these corrections.The so-called pole approximation (PA) [ 43,44] (see also [ 45] and references therein for the general concept) is based on a systematic expansion of the cross section about the W/Z resonance, allowing for a split of the O(α s α) corrections into well-defined, gauge-invariant parts and a classification of these parts according to their impact on the production and decay subprocesses.To be precise, in the PA the corrections are split into factorizable and non-factorizable contributions, where the former incorporate radiative corrections to the production or decay mode and the latter non-factorizable corrections originate from contributions including soft photon exchange between production and decay.In Refs.[ 43,44] these subsets were calculated (and implemented in the program Rady, which is the basis of the NLO corrections discussed in Refs.[ 11,18,19]) except for the "initial-initial" factorizable contributions, which contain double-real and two-loop corrections involving only the initial state and are expected to be small.In contrast to the narrow-width approximation (NWA), which treats the intermediate W/Z bosons as stable, the PA describes off-shell effects of the W/Z bosons in the vicinity of the resonance.Using the NWA, in Ref. [ 46] the QCD×QED corrections to the total DY-like Z-production cross section were obtained by an abelianisation procedure of the known QCD NNLO results.Inclusive results for the mixed QCD-EW corrections to on-shell Z production were calculated in [ 47,48] and fully differential results in Refs.[ 49,50,51].The two-loop formfactor for Z-boson production in quark-antiquark annihilation was calculated in Ref. [ 52].
Since physics beyond the SM might also show up in the tails of invariant-mass or transverse-momentum distributions outside the resonance regions, it is important to provide information about the size of O(α s α) corrections beyond the PA or NWA.To this end, first technical steps have been made.In Ref. [ 53] results for the two-loop integrals needed for DY-like W/Z-boson production were given in terms of iterated integrals, and recently it has been shown that it is indeed possible to write the needed integrals in terms of multiple polylogarithms [ 54,55].A first step towards the full O(α s α) corrections to off-shell DY processes is the calculation of the gauge-invariant O(N f α s α) two-loop corrections to single W/Z-boson production which are enhanced by the number of fermion flavours N f in the Standard Model (SM) and result from diagrams including closed fermion loops and additional gluon exchange or radiation.The necessary genuine two-loop O(α s α) corrections to the vector-boson self-energies were first calculated in Refs.[ 56,57,58,59,60,61] a long time ago.
In this paper, we present first results of an evaluation of the O(N f α s α) corrections to DY-like W/Z-boson production including a reevaluation of the occurring two-loop self-energies by reducing the two-loop integrals with current standard methods [ 62,63] to a set of master integrals suitable for numerical evaluation.The master integrals in D = 4 − 2 dimensions are solved by deriving differential equations in Henn's canonical form [ 64,65] and subsequent integration to obtain the results as a Laurent expansion in in terms of generalized polylogarithms up to weight three.Furthermore, besides the corrections containing one-particle-irreducible two-loop (sub)diagrams the O(N f α s α) corrections contain reducible contributions which either involve a product of two one-loop subdiagrams or one-loop subdiagrams with an additional possibly unresolved QCD parton in the final state.We evaluate the O(N f α s α) corrections to single W/Z-boson production in a fully differential manner and study their effect on the (transverse) invariant-mass and transverse-momentum spectra of the W and Z boson, respectively.The calculation of virtual corrections of O(N f α s α) involves the issue of extending a gauge-invariant scheme for treating the W/Z resonance to this order.To solve this problem, we describe the generalization of the complex-mass scheme [ 66] (see also Ref. [ 45]), which is a standard method for a gauge-invariant treatment of resonances at NLO, for the application to W/Z resonances at O(α s α).Note that the consideration of N f -enhanced O(α s α) corrections is already sufficient for this step, since absorptive parts in the W/Z propagators necessarily involves closed fermion loops.
The paper is organized as follows: In Section 2 we briefly summarize the properties of the O(N f α s α) corrections, give explicit results of the O(α s α) contributions to the EW gauge-boson self-energies in terms of two-loop master integrals and discuss their renormalization and the generalization of the complex-mass scheme needed at O(N f α s α).Furthermore, we describe the reduction of the occurring two-loop diagrams to master integrals and the calculation of the integrals.The explicit results of the master integrals and the transformations needed to obtain Henn's canonical form of the differential equa-tions are provided in App. A. We discuss the phenomenological impact of O(N f α s α) corrections on transverse-momentum and invariant-mass distributions in Section 3, and Section 4 provides a short summary.
2 Details of the calculation 2.1 Survey of diagrams and structure of the calculation We consider the two types of DY-like pp scattering processes with ± denoting either e ± or µ ± .At leading order (LO), the charged-current process is entirely due to q q annihilation, but the neutral-current process receives contributions from both q q annihilation and γγ scattering.The γγ channel [ 16,18,19,20,67], however, delivers only a small fraction to the neutral-current cross section and does not develop a Z-boson resonance.Already the NLO EW corrections to this channel turn out to be phenomenologically irrelevant [ 19], so that we do not include the γγ channel in our calculation of O(N f α s α) corrections in the following, but restrict our calculation to q q( ) annihilation.NNLO corrections generically receive contributions from (i) "virtual-virtual" (vv-1PI) contributions involving one-particle-irreducible (1PI) twoloop (sub)diagrams, (ii) "virtual-virtual" (vv-red) contributions induced by diagrams containing reducible loop parts of the type (one-loop)×(one-loop), (iii) "real-virtual" (rv) contributions resulting from one-loop diagrams with one extra emission of a possibly unresolved particle (gluon, quark, photon), and (iv) "real-real" (rr) contributions induced by tree-level diagrams with two extra emissions of possibly unresolved particles.
Our focus on NNLO corrections of the order O(N f α s α) that are enhanced by the number N f of fermion flavours in the SM and on 2 → 2 scattering processes with four massless external fermions restricts the possible contributions to those categories considerably.
In order to produce the enhancement factor N f in loops, a closed fermion loop has to be present either in a one-or two-loop subdiagram.For the considered process class f1 f 2 → f3 f 4 with f i denoting the external massless fermions, those fermion loops only occur in gauge-boson self-energies. 1 This restricts the set of 1PI two-loop diagrams to the self-energy insertions shown in Fig. 1.To those EW self-energies, only contributions from closed quark loops contribute at O(N f α s α).The vv-red contributions are diagrammatically illustrated in Fig. 2; they combine the closed fermions loops (with either quarks or leptons in the loop) in the EW gauge-boson propagators with the NLO QCD loop Figure 1: One-particle-irreducible virtual-virtual (vv-1PI) two-loop contributions to DYlike processes at O(N f α s α).In the first diagram the two-loop O(N f α s α) self-energy insertions are shown, whereas the second and third diagrams show the finite gauge-boson fermion counterterms described in Section 2.3.
(a) Reducible virtual-virtual contributions within one diagram Figure 2: Different types of reducible virtual-virtual (vv-red) diagrams contributing at O(N f α s α) to DY-like processes, where the relative orders of the loop corrections are indicated in the vertex blobs.diagrams in all possible ways.The rv contributions similarly combine the closed fermions loops in the EW gauge-boson propagators with the real NLO QCD corrections.Figure 3 shows some of the corresponding diagrams for the gluon-emission channel, while their crossed counterparts from qg scattering are not depicted explicitly.Note that at O(N f α s α) there are no rr corrections with double real emission.Such contributions arise from g → q q and γ/Z → f f splittings at O(N f α 2 s ) and O(N f α 2 ), respectively, but at O(N f α s α) the corresponding contributions combine a gluon and a photon/Z splitting for a single spinor chain and, thus, vanish due to colour conservation.
In the following we describe in some detail the calculation and results of the twoloop contributions to the self-energies and the corresponding complex renormalization within the complex-mass scheme, which is employed for the gauge-invariant description of the gauge-boson resonances.The evaluation of the matrix elements including those selfenergies as well as the evaluation of the reducible vv and rv contributions proceeds fully analogously to the NLO QCD and EW calculations.Since there are no double-unresolved infrared-singular rr contributions, but only infrared singularities of NLO QCD type, we simply employ standard NLO QCD subtraction techniques to combine the vv-red and rv corrections; the vv-1PI corrections do not involve infrared singularities.
In total, we have performed two completely independent calculations, leading to two independent implementations, the results of which are in mutual numerical agreement.The first calculation builds on the Fortran program Rady, which is the basis for the NLO EW and QCD calculations described in Refs.[ 11,18,19].In order to generalize Rady to the calculation of O(N f α s α) corrections, we just had to dress all ingredients of the NLO QCD calculation with the EW gauge-boson self-energy contributions of O(α) and to add the relevant two-loop contributions to the EW gauge-boson self-energy corrections.Infrared singularities are handled with standard QCD dipole subtraction [ 68].The graphs and amplitudes for the two-loop self-energies were generated with FeynArts [ 69,70] and further algebraically reduced with inhouse Mathematica routines and KIRA [ 63,71].The genuine two-loop corrections of O(N f α s α) contain Goncharov Polylogarithms (GPLs) [ 72,73] up to weight three.In the first calculation the numerical evaluation of the necessary GPLs was performed in two steps.In the first step the GPLs were reduced by hand to Harmonic Polylogarithms (HPLs) [ 74] following the methods introduced in Ref. [ 75] and in the second step the HPLs were evaluated using the Fortran program CHAPLIN [ 76].The second, independent calculation of the corrected cross sections employs antenna subtraction [ 77] to handle infrared singularities present in the reducible vv-red and rv O(N f α s α) corrections, which were obtained analogously to the first calculation by dressing the NLO QCD calculation with EW gauge-boson self-energies of O(α).The twoloop self-energies were generated with QGraf [ 78] and algebraically reduced to scalar integrals via Matad [ 79] and FeynCalc [ 80,81].The reduction to master integrals was again performed with KIRA to get the final result in Mathematica.The GPLs contained in the genuine two-loop O(N f α s α) corrections were evaluated using the C++ library GiNaC [ 82].

Electroweak gauge-boson self-energies at O(α s α)
As explained above, the only 1PI two-loop building blocks required for the O(N f α s α) corrections are the EW gauge-boson self-energies at this order.More precisely, only the transverse parts Σ V V T (k 2 ) (V V = γγ, γZ, ZZ, WW) of those self-energies are needed, where k 2 denotes the virtuality of the gauge bosons V, V .For the precise relation between the two-point vertex functions Γ V V and the self-energies Σ V V we follow the conventions of Ref. [ 45] (identifying Σ W W ≡ Σ W and defining In the following, only the transverse self-energy parts Σ V V T will be considered, because the longitudinal parts Σ V V L are not relevant in our calculation.By definition, we do not include tadpole contributions in Σ V V T , since tadpoles fully cancel in the considered on-shell renormalization scheme, i.e. our results on Σ V V T correspond to the "parameterrenormalized tadpole scheme" (PRTS) as defined in Refs.[ 45,83].We decompose the O(α s α) contribution Σ V V T,(αsα) (k 2 ) to the self-energies according to where Σ V V T,(αsα),1PI comprises all genuine irreducible two-loop diagrams, as shown in Fig. 4, and Σ V V T,(αsα),δm represents all fermion loops with insertions of the quark-mass counterterms.In D = 4 − 2 dimensions, with µ denoting the arbitrary reference mass scale of dimensional regularization, the mass renormalization constants δm q in the on-shell scheme (see Fig. 4) is given by where C F = 4 3 denotes the quadratic Casimir factor of the fundamental representation of SU (3).Note that no other one-loop counterterm insertions in one-loop diagrams are relevant at O(N f α s α), because the only other potentially relevant renormalization constants of O(α s ) are the quark-field renormalization constants, but their contributions to Σ V V T,(αsα) fully cancel.
The gauge-boson self-energies are first expressed in terms of the two-loop two-point integrals 2 ) for m 1 = m 2 .Dotted lines represent a squared propagator.
where a graphical representation of these integrals is shown in Fig. 5.The prefactor in this definition is chosen in such a way that reducible integrals decompose into the product of the standard one-loop integrals defined in Refs.[ 83,45].The integral functions S abcde obey some obvious symmetries, which are exploited in the formulas below, Using Laporta's algorithm [ 62] as implemented in the program KIRA [ 63,71], we reduce the occurring two-loop integrals in terms of the minimal set of master integrals illustrated in Fig. 6.
For the transverse parts of the self-energies of the neutral EW gauge bosons we explicitly get where

The sums
q extend over all quark flavours q ∈ {u, d, c, s, t, b} with relative electric charges Q q and third components I 3 w,q = ± 1 2 of the weak isospin, and the vector and axial-vector couplings of quark q to the Z boson are denoted as Keeping the full dependence on D = 4 − 2 in order to facilitate the later specialization to specific mass patterns, the auxiliary functions f k (k = 1, . . ., 4) are given by ) ) S 02020 , (2.17) with suppressed arguments of the integral functions S abcde (s, m 2 , m 2 ).The O(α s α) contributions to the transverse part of the W-boson self-energy is given by where the sums j extend over the three generations of up-type and down-type quarks u j and d j , respectively.The auxiliary function f k (k = 5, 6) are given by with the Källen function and the arguments of the integral functions given by S abcde (s, m 2 1 , m 2 2 ).Note that the interchange (m d j ↔ m u j ) of the up-and down-type quark masses in (2.19) and (2.20) also concerns the arguments of the integral functions; this change of arguments can, however, be achieved by rearranging labels in S abcde using (2.7).
For massless fermions, only the functions f 1 and f 5 are relevant and given by to the relevant order in .
In addition to the self-energies Σ V V T (s) for non-vanishing s, we in particular need the W-boson self-energy at zero-momentum transfer in the application of the G µ inputparameter scheme below.In this limit the two-mass two-loop tadpole integrals T abc , defined by ] c , (2.25) are needed in addition.They obey the following symmetry relations (2.26) Since the numerical evaluation of Σ W T,(αsα) (0) is somewhat non-trivial in the above representation, we here give an explicit form for Σ W T,(αsα) (0) suitable for a numerical evaluation, which was obtained by explicitly expanding the master integrals about s = 0 with the help of the differential equations used to calculate them (see App. A), (2.28) The auxiliary functions f k (k = 7, 8) are given by where λ 0 is obtained by evaluating λ in (2.23) at s = 0, and the integrals S abcde have the arguments S abcde (0, m2 1 , m 2 2 ).Note that the appearing master integrals S 02020 and S 00202 are actually products of one-loop tadpole integrals and can be expressed in terms of T abc via The limits of the functions 2 ) (k = 7, 8) in which one of the two quark masses vanishes can be obtained by simply evaluating (2.29) and (2.30) with the corresponding mass set to zero.In the case in which both quark masses are zero the whole contribution of the corresponding massless quark generations to Σ W T,(αsα) (0) vanishes because of dimensional reasons.
The O(N f α s α) corrections to the EW gauge-boson self-energies have been calculated some time ago in Refs.[ 56,57,58,59,60,61].We have compared our results with the ones given in Ref. [ 56] and find full analytical agreement in the case of vanishing quark masses.For non-vanishing quark masses we find numerical agreement for Σ V V T,(αsα) (k 2 ) with those results 2 .

Renormalization and complex-mass scheme
In our calculation of O(N f α s α) corrections we employ straightforward generalizations of the on-shell renormalization schemes and their complexified versions used in the NLO EW calculations for DY-like processes described in Refs.[ 11,18,19].At NLO the real formulations and their complex generalizations are described in Refs.[ 83,45] and Refs.[ 66,45], respectively.
Since the reducible vv and rv contributions only involve one-loop subdiagrams, their calculation does not require any generalization beyond NLO.The only generalization to NNLO concerns the calculation of the counterterms required in the gauge-boson two-point functions and in the gauge-boson-fermion vertices for the EW gauge bosons.However, owing to our restriction to the N f -enhanced O(α s α) corrections, all relevant contributions to the needed counterterms originate from the contributions to the EW gauge-boson selfenergies Σ V V T considered above.In detail, we need the O(N f α s α) contributions to the following renormalization constants in the complex-mass scheme [ 66,45]: the gaugeboson mass renormalization constants δµ 2 W , δµ 2 Z , the gauge-boson field renormalization constants δZ V V , the renormalization constants δc w for the weak mixing angle, and the charge renormalization constant δZ e .
The W-and Z-boson mass parameters µ 2 V (V = W, Z) are defined as the locations of the poles in the complex k 2 plane of the W/Z propagators and are decomposed into real and imaginary parts according to where the real mass and width parameters M V and Γ V deviate from their counterparts M V,OS and Γ V,OS in the real on-shell (OS) scheme at the two-loop level.In good approximation, the connection is [ 45] M Since the on-shell mass and field renormalization of the EW gauge bosons is simply based on some momentum subtraction for the vertex two-point function, the perturbative contributions to the renormalization constants δµ 2 V and δZ V V are in one-to-one correspondence with the corresponding orders in the required self-energies Σ V V T .Denoting again the order of the contributions by some subscript "(α s α)" for O(α s α), we therefore get ) where Σ V V (k 2 ) ≡ (∂Σ V V /∂k 2 )(k 2 ).In quantities, in which the distinction between O(N f α s α) and O(α s α) is not necessary, we simply write (α s α) as subscript.
In order to avoid the evaluation of self-energies with complex k 2 , we follow the "simplified version" of the complex-mass scheme based on Taylor expanding Σ V up to the relevant order.This leads to Note that some care is required in order to catch all the relevant terms in the evaluation of δµ 2 V above.Firstly, there is no O(α s ) contribution to Σ V V T at NLO, and Γ V counts as O(α), so that no additional terms of O(α s α) arise from higher terms in the Taylor expansion (2.37) of Σ V V T at NLO.Secondly, we do not need to include any extra term like c W T as introduced in Refs.[ 66,45] that occurs at NLO EW as a consequence that W is rather a branch point than a pole of the W propagator, because this subtlety arises from an infrared singularity in on-shell diagrams with photon exchange of the W boson.However, at O(α s α), the self-energies Σ V V T do not involve infrared singularities, i.e.Σ W T,(αsα) is analytic at k 2 = µ 2 W , and no extra terms occur.The renormalization constants δc w and δs w for the (complex) cosine and sine of the weak mixing angle are fixed by the condition that the identity holds both for bare and renormalized quantities.Again, since Σ V V T does not receive O(α s ) contributions, we get for the contributions to δc w and δs w at O(α s α) δs w,(αsα) The determination of the charge renormalization constant δZ e beyond NLO deserves some care.It is derived from the condition that the renormalized fermion-photon vertex for on-shell fermions does not receive a correction in the "Thomson limit" of vanishing photon momentum.Using symmetry arguments similar to the arguments based on a Ward identity in quantum electrodynamics (QED), it is possible to express δZ e in terms of gauge-boson self-energies instead of vertex-correction formfactors.For the SM this derivation based on Lee identities is spelled out in App.C of Ref. [ 45] at NLO. Taking the fermion of the renormalization condition in the Thomson limit as a lepton, the only source for O(α s α) contributions in a generalization of this derivation are closed quark loops in the gauge-boson self-energies and related self-energies involving Goldstone bosons.Since those self-energies do not receive O(α s ) contributions, no reducible O(α s α) corrections occur in the derivation.Therefore, all identities that are given in App.C of Ref. [ 45] for O(α) corrections are valid for O(α s α) as well, with all corrections but the self-energies of the gauge-boson and Goldstone-boson sectors vanishing.The final result for δZ e then takes a form fully analogous to NLO, which holds as a consequence of Slavnov-Taylor (ST) identities.
The same result can be obtained more directly within the background-field method (BFM) [ 84,85,86,87,88], which is applied to the SM in Refs.[ 89,45].Owing to the gauge invariance of the background-field effective action, the Ward identities for the fermion-photon vertex takes the same simple form as in QED to all perturbative orders.In particular, Eq. (2.42) holds within the BFM to all orders. 3Consequently, the QEDlike result (2.42) for δZ e trivially carries over to the SM in its BFM formulation in each perturbative order.We note in passing that we have checked explicitly all ST identities for the O(α s α) contributions to the EW gauge-boson two-point functions Γ V V considered in the previous section.At O(α s α) these ST identities are formally identical to the BFM Ward identities given in Eqs. ( 59)-( 61) in Ref. [ 45].
The renormalization constants defined above enter the amplitudes for the O(N f α s α) corrections in two different ways.On the one hand, they are part of the renormalized gauge-boson self-energies Σ where we set µ A = 0. On the other hand, they enter the gauge-boson-fermion counterterms, where they change the LO coupling factors g σ V f f with chirality σ = ± by the factors δ ct,σ W f f ,(αsα) = δZ e,(αsα) − δs w,(αsα) where ) for a fermion f with relative electric charge Q f and third component I 3 w,f = ± 1 2 of weak isospin.All gauge-boson field renormalization constants cancel in the sum over all contributions, but in the above form the quantities Σ V V R,T,(αsα) and δ ct,σ V f f ,(αsα) are all ultraviolet finite individually.

Electroweak input-parameter scheme
In the following, we use the Fermi constant G µ as input for the EW coupling strength, instead of the fine-structure constant α(0) = e 2 /(4π), along with the gauge-boson masses µ W , µ Z , i.e. we work in the so-called "G µ -scheme", as e.g.described in Refs.[ 11,45].Formally, we derive the following value for α from G µ , i.e. we take α Gµ as a real quantity.The arguments given, e.g., in Sect.6.6.4 of Ref. [ 45] that this is a legal procedure in O(α) easily carry over to O(α s α).This reparametrization of α leads to the change in the charge renormalization constant, δZ e,(αsα) Gµ = δZ e,(αsα where ∆r quantifies the corrections to muon decay [ 90,91].The O(α s α) contribution ∆r (αsα) to ∆r is entirely given by the fermion-loop contributions to the gauge-boson self-energies and, thus, follows from the O(α) result [ 45,83,90,91] for ∆r with the corresponding substitution for the self-energies, where we have used (2.43).Similar to the situation at NLO, using the G µ scheme eliminates the sensitivity of the corrections to DY production to the light quark masses, since the mass-singular contribution Σ AA T,(αsα) (0) cancels in δZ e,(αsα) Gµ , and the universal corrections to the ρ-parameter are absorbed into the charged-current coupling α Gµ /s 2 w .Following the arguments of Sect.6.6.4 of Ref. [ 45], we can take the gauge-boson widths Γ V as independent input parameters, although they are strictly speaking not free parameters of the SM.We uniformly set them to their experimental values given below.Using different width parameters in LO predictions and corrections would unnecessarily obscure the impact of the calculated O(N f α s α) corrections we want to discuss.
3 Numerical results

Input parameters and event selection
The setup for the calculation is widely taken over from Refs.[ 43,44].The choice of input parameters closely follows Ref. [  We convert the on-shell (OS) masses and decay widths of the vector bosons to the corresponding pole masses according to (2.34).The electromagnetic coupling constant is set according to the G µ scheme.The masses of the light quark flavours (u,d,c,s) and of the leptons are neglected throughout.The CKM matrix is chosen diagonal in the third generation, and the mixing between the first two generations is parametrized by the following values for the entries of the quark-mixing matrix, For reference, in Tab. 1 we give numerical values for the gauge-boson-fermion renormalization constants δ ct,σ V f f ,(αsα) defined in Eqs.(2.45) and (2.46) for V = W, Z. 4 The numerical values are calculated using the complex-mass scheme and the G µ input-parameter scheme, as described above, using the input values of Eq. (3.1).Note that in the OS scheme diagrams containing the gauge-boson-fermion renormalization constants in Tab. 1 dictate the size of the vv-1PI O(N f α s α) corrections close to the resonance of the amplitude.Therefore, in the resonance regions the size of the vv-1PI O(N f α s α) corrections is at the permille level due to the smallness of δ ct,σ V f f ,(αsα) .For the PDFs we consistently use the NNPDF2.3set [ 93], i.e. the NLO and NNLO QCD-EW corrections are evaluated using the NNPDF2.3QEDNLO set [ 94], which also includes O(α) corrections.The value of the strong coupling α s (M Z ) quoted in Eq. (3.1) is dictated by the choice of these PDF sets.The renormalization and factorization scales are set equal, with a fixed value given by the respective gauge-boson mass, with V = W, Z for W and Z production, respectively.
For the experimental identification of the DY process we impose the following cuts on the transverse momenta and rapidities of the charged leptons, p T, ± > 25 GeV, |y ± | < 2.5, (3.4) and an additional cut on the missing transverse energy GeV, (3.5) in case of the charged-current process.For the neutral-current process we further require a cut on the invariant mass M of the lepton pair, M > 50 GeV, (3.6) in order to avoid the photon pole at M → 0. Since there is no photon emission involved in the corrections of O(N f α s α), the issue of dressed leptons and photon recombination is not relevant for the calculated corrections.

Corrections to differential distributions
Figure 7 shows the relative correction of O(N f α s α) to the distributions in the invariant mass M of the lepton pair + − ( = e, µ) for Z production and in the transverse invariant mass M T,ν of the pair ν + for W + production, where M T,ν is the invariant mass that is calculated by taking only the transverse components of the respective three-momenta into account.In the calculation of δ the O(N f α s α) contribution dσ (N f αsα) to the differential cross section is normalized to the LO cross section dσ LO bin by bin in the histograms, where both contributions are evaluated with the same PDF set, so that δ is practically independent of the factorization scale µ F .The correction δ mildly depends on the renormalization scale µ R via its proportionality to α s (µ R ).Apart from the full O(N f α s α) contribution (red curves) we show the part of the correction that is furnished by reducible diagrams only (green curves) and the contribution delivered by the first two fermion generations (blue curves).In Fig. 7 we depict the regions of low and high M and M T,ν separately, where the resonant contributions of the intermediate W/Z bosons is contained in the low-mass plots on the l.h.s..More precisely, the whole region with M T,ν < ∼ M W is dominated by resonant W bosons, while the Zboson resonance shows up only for M ∼ M Z .We only show the relative corrections δ to illustrate their impact; results on the absolute predictions for the shown spectra and their distinctive shapes are discussed in numerous papers (see, e.g., Refs.[ 18,19]).As already expected from the size of the renormalization constants given in Tab. 1, from the results on O(α s α) corrections for stable W/Z bosons, and from the results in pole approximation [ 44], the impact of O(N f α s α) corrections is at the level of permille, and thus phenomenologically unimportant, in all regions where resonant W/Z bosons dominate the cross section.Away from the resonance regions, the corrections grow to 1.5-2%, which is the typical size of the corrections for M and M T,ν values of 300−1000 GeV.Corrections of this size are in fact phenomenologically relevant in those off-shell tails, in particular in the search for traces of new physics, as potentially induced by Z or W bosons.
It is interesting to note that the contribution of reducible corrections dominates over the impact of irreducible diagrams whenever the O(N f α s α) correction is sizeable.Furthermore, we notice that the contributions of the individual fermion generations are generically of similar size, i.e. there is no suppression of the third generation (with massive quarks) w.r.t. to the other generations.In fact for Z production the impact of the third generation is even larger than the sum of the first two.We note in passing that the t t threshold is observable in the M spectrum at M ∼ 2m t ≈ 346 GeV (lower right plot in Fig. 7) in the full O(N f α s α) correction (red) and its reducible part (green), but of course not in the contribution of the first two fermion generations (blue).From the comparison of the   three different curves we conclude that neither a neglect of the third quark generation nor the approximation by setting m t and m b to zero provides a viable approximation for the corrections.Such approximations are often useful for QCD corrections at low or high energies; for EW corrections such approximations in general fail, since the EW gauge-boson masses M W ∼ M Z enter the renormalization conditions and provide an additional scale.
Figure 8 shows the relative O(N f α s α) correction δ to the leptonic transverse-momentum distributions in the low-and high-energy regions.At LO, the regions in which resonant W/Z bosons dominate the spectra are characterized by k T, < ∼ M V /2 (V = W, Z), i.e. at the left side of the Jacobian peaks at k T, = M V /2.Note, however, that jet emission from the initial-state partons transfers some transverse momentum to the W/Z bosons, so that in the presence of QCD corrections (and to a lesser extent also in the presence of photonic corrections which are not discussed here) the k T, regions above the Jacobian peaks receive contributions that are enhanced by a W/Z resonance.This well-known effect leads to extremely large QCD corrections for k T, > M V /2 which grow to some 100%.This does not mean that perturbation theory does not work in this region, but only that the LO prediction is not a good approximation for the differential cross section there.This enhancement mechanism of NLO QCD over LO contribution also leads to an enhancement of the O(α s α) correction, since it involves jet emission as well.In Fig. 8 the suppression of the LO cross section for k T, > M V /2 results in large values of δ that grow even to about 15% for k T, > ∼ 250 GeV.If we used the QCD-corrected cross section in the normalization of δ, we would get relative corrections in the few-% range, which is of the expected size of the O(α s α) corrections as observed in the invariant-mass spectra above.The dominance of the real QCD corrections via the described recoil mechanism is also the reason for the extreme dominance of the reducible contributions in the O(N f α s α) corrections, because the irreducible corrections do not involve real-emission effects.As already noticed in the discussion of the invariant-mass spectra above, the contribution of the third generation relative to each of the first two is higher for Z production than for W production; this feature is even more pronounced in the transverse-momentum spectra.
Finally, we mention that we observe only permille corrections of O(N f α s α) to the integrated cross section and to differential distributions that are entirely dominated by resonant W or Z bosons, such as distributions in the lepton rapidities.

Summary
Next-to-next-to-leading-order corrections of mixed QCD×EW type seem to be the largest component of the yet unknown radiative corrections to Drell-Yan-like W/Z production at fixed perturbative order, at least for off-shell W/Z bosons.In the vicinity of the W/Z resonances, the corrections are known in the form of a pole approximation up to the corrections that solely concern the initial state, which are supposed to be small.Recent evaluations of those initial-state corrections for on-shell Z bosons have confirmed this expectation.For off-shell W/Z production, several ingredients have been presented in recent years, including results for rather complex two-loop integrals, but no cross-section predictions have been presented yet.This paper takes a first step towards the numerical evaluation of the O(α s α) corrections by presenting results on the corrections of O(N f α s α), which are nominally enhanced by the number N f of fermion generations.These corrections comprise all diagrams with closed fermion loops and form a gauge-invariant part of the full O(α s α) corrections.
The genuine two-loop part of the calculation involves only self-energy complexity and was feasible by a straightforward application of current two-loop techniques.The twoloop integrals were reduced to master integrals with the help of Laporta's algorithm as implemented in the program KIRA, and the master integrals were evaluated via differential equations.We have successfully compared our results on the EW gauge-boson selfenergies to existing results in the literature and give explicit analytical results to allow for cross-checks with upcoming similar calculations.Generally, the description of resonance processes including higher-order corrections is delicate, in particular because of issues with gauge invariance.In order to guarantee a gauge-invariant description that is uniformly valid on resonance and in off-shell regions, we have generalized the complex-mass scheme, which is a standard procedure for treating resonances at NLO.It is interesting to note that the consideration of all O(N f α s α) corrections is already sufficient for the generalization of the complex-mass scheme for the full O(α s α) corrections, since the W/Z propagators that develop the resonance are affected at O(α s α) only by diagrams involving closed quark loops.
Concerning real corrections, focusing on O(N f α s α) leads to drastic simplifications in comparison to the full O(α s α) corrections as well.Since the O(N f α s α) corrections do not involve photon emission, but only up to a single emission of QCD partons, one-loop subtraction techniques are sufficient to treat infrared singularities.Specifically, we have applied dipole and alternatively antenna subtraction.
Our discussion on numerical results shows that O(N f α s α) corrections to observables that are dominated by resonant W/Z bosons, such as integrated cross sections or rapidity distributions, are at the permille level and thus phenomenologically negligible.This could be already concluded from the existing results on on-shell W/Z production or from results in pole approximation.Off-shell regions in differential distributions, however, receive sizeable corrections.For instance, the invariant-mass distribution for lepton pairs in Z production and the respective transverse-invariant-mass distribution in W production receive corrections at the level of 1.5-2% for (transverse) invariant masses of ∼ 300−1000 GeV.Nominally, transverse-momentum distributions of leptons even receive corrections of the order of 10% or more above the Jacobian peak at transverse momenta ∼ M V /2 (V = W/Z) if corrections are normalized to leading-order predictions.However, those corrections reduce to the few-% level after normalizing them to full predictions, since leading-order predictions systematically underestimate the distribution above the Jacobian peaks, which is a well-known phenomenon.
Considering the remaining theoretical uncertainty induced by missing higher-order corrections to W/Z production at hadron colliders, we have to keep in mind that the still unknown O(α s α) corrections without nominal N f enhancement are expected to be not smaller than the corrections of O(N f α s α).This is due to the enhancement of EW corrections at high energies originating from double (Sudakov) and single logarithms at NLO EW, which are known to factorize from QCD corrections in higher orders.With both NLO QCD and NLO EW corrections at the (known) level of some 10% in the TeV range of invariant masses, additional O(α s α) contributions at the few-% level can be expected.The presented results on O(N f α s α), thus, do not directly reduce the current theoretical uncertainty, but represent a relevant contribution to the full O(α s α) corrections and can serve as an estimate for the order of magnitude of missing corrections at this order.

A Calculation of the master integrals via differential equations
In this appendix we briefly describe the calculation of the two-loop master integrals via differential equations, which is based on transformations of the set of master integrals into Henn's canonical form [ 64,65] and subsequent integration of the new basis integrals in terms of a Laurent expansion in including terms up to order 1 , which involve Goncharov polylogarithms (GPLs) up to weight three.We start by describing the procedure for the general case of different non-vanishing masses and present some special cases with much simpler results afterwards.The most simple case, in which all masses are zero, has also been checked by direct integration with Feynman parameters.
Apart from the results outlined in the following, we have worked out an alternative solution for the master integrals, which is based on a more involved transformation to the canonical form, but leading to somewhat simpler expressions for the integrals.Numerically the two sets of obtained master integrals are in mutual agreement.

A.1 General case of different non-vanishing masses
We first change the basis of master integrals S abcde used to express the 1PI parts of the EW gauge-boson self-energies to the following set of basis functions, where we have used the shorthands Here and in the following, squared masses are always assumed to possess an infinitesimally small negative imaginary part, i.e. m 2 ≡ m 2 − i0.Moreover, we replace the kinematical variable s in favour of the dimensionless variable x, which rationalizes In terms of the kinematical input, the variable x is calculated according to where the two versions for λ > 0 are just distinguished to improve numerical stability.In order to ensure that s = 0, which will be our initial condition for solving the differential equation, corresponds to x = 0, we assume m 2 > m 1 in the following.The case m 2 < m 1 can be handled upon interchanging the mass values before the calculation of master integrals and appropriately interchanging the obtained master integrals using the symmetry relations (2.7).The transformation to the set of functions F , which is inspired by a corresponding but simpler transformation for the one-loop bubble integral, brings the differential equation of the master integrals for the evolution in the variable s (keeping the masses m 1 , m 2 constant) into the canonical form where f results from F by some rescaling, Schematically, the matrix A is given by where A 4 and A 6 are the 4 × 4 and 6 × 6 matrices which have the element A 44 = A 4,44 = A 6,11 = 0 in common and 0 m×n is the zero matrix of the indicated geometry.The explicit entries of A 4 and A 6 are given by with the auxiliary functions (A.9) Owing to the block structure (A.7) of the matrix A, the first four components of f and the last six components of f each define an independent system of linear differential equations, which can be solved independently; the fact that f 4 is part of either system does not disturb this feature.
As initial condition for the evolution of f (x, r) in x, we take the values f (0, r) corresponding to s = 0, where the functions S abcde reduce to vacuum integrals of the type defined in (2.25).For the functions in F this leads to the initial values (A.10) The integrals T 022 are just products of simple one-loop vacuum integrals, which are easy to calculate.The integral T 122 was first expressed in terms of T 111 with the help of KIRA, and T 111 was calculated by solving the corresponding Feynman parameter integral.The result for T 111 was also checked against the one published in Ref. [ 95].For the rescaled functions f , the initial values explicitly read With these initial values, the integration of the system (A.5) in terms of GPLs is straightforward.Since the functions f k (x, r) with k = 4, 5, 9 are constant in x, their solutions are trivially given by f k (x, r) ≡ f k (0, r), k = 4, 5, 9. (A.12) For the remaining functions f k (x, r), we give the results in terms of coefficients f (j) (x, r) of the Laurent series up to the relevant order in .Up to order 0 , those functions read For the evaluation of the self-energies given in Section 2.2, the functions f k (x, r), the results of which are getting more lengthy and untransparent, are needed as well; we provide those functions in an ancillary file supplementing the online version of this article.
To finally reconstruct the relevant master integrals S abcde in terms of a Laurent series in powers of , we first have to convert the coefficients f 2 ) of the components of F as defined in (A.6).By convention, we do not expand the global factor Γ(1 + ) 2 (4π) 2 contained in F and define k (s, m 2 1 , m 2 2 ) = f  k , however, get somewhat lengthy because of the explicit appearance of in the defining equations.Moreover, the basis set of master integrals used in the selfenergies in Section 2.2 is not identical with the one used in (A.1), i.e. a further change of basis has to be performed.Instead of reproducing unnecessarily lengthy formulas here, we provide the coefficients S (j) abcde needed for the self-energies in terms of the coefficients F (j) k given above in the mentioned ancillary file.

A.2 Two equal non-vanishing masses
In this appendix we consider the calculation of the master integrals S abcde for the special case m = m 1 = m 2 , in which the number of independent master integrals is reduced compared to the general case of the previous section owing to the symmetry relations (2.7).To solve the system of differential equations obeyed by those master integrals we consider the following basis of five functions, and have the element A 33 = A 3,33 = A 3,11 = 0 in common.Each of the matrices A 3 , A 3 defines a 3-dimensional system of linear ordinary differential equations that can be solved independently.

A.4 Massless case
The required master integrals for m 1 = m 2 = 0 can be obtained upon specializing the results from the previous section or via Feynman parameter integration in a straightforward way.The independent integrals are explicitly given by S 10110 = Γ(1 + ) The remaining ones follow from those via the symmetry relations (2.7).

Figure 3 :
Figure 3: Different types of real-virtual (rv) diagrams contributing at O(N f α s α) to DYlike processes, where the relative orders of the loop corrections are indicated in the vertex blobs.

Figure 4 :
Figure 4: Diagrams contributing to the EW gauge-boson self-energies at O(N f α s α), which all involve closed quark loops.In the first line the contributions to Σ V V T,1PI and in the second line the contributions to Σ V V T,δm are shown.

Figure 5 :
Figure 5: Two-loop sunset topology, corresponding to the self-energy integral S abcde defined in Eq. (2.6).

. 2 )
While b-quarks appearing in closed fermion loops have the mass m b given in Eq. (3.1), external b-quarks are taken as massless.

Figure 7 :
Figure 7: Relative O(N f α s α) corrections to distributions in the transverse invariant mass of the W bosons (upper plots) and in the invariant mass of the Z boson (lower plots), where the complete O(N f α s α) corrections are compared to the contribution originating from reducible graphs and to the contribution delivered by the first two fermion generations.

Figure 8 :
Figure 8: Relative O(N f α s α) to transverse-momentum distributions for Wboson (upper plots) and Z-boson production (lower plots), again with a comparison of full O(N f α s α) corrections to its reducible parts and to the contribution of the first two fermion generations.