Two-loop mixed QCD-EW corrections to neutral current Drell-Yan

We present the two-loop mixed strong-electroweak virtual corrections to the neutral current Drell-Yan process and we provide, in ancillary files, the explicit formulae of the infrared-subtracted finite remainder. The final state collinear singularities are regularised by the lepton mass. The evaluation of all the relevant Feynman integrals, including those with up to two internal massive lines, has been worked out relying on analytical and semi-analytical techniques, in the case of complex-valued masses.


Introduction
At hadron colliders, the production of a pair of leptons each with large transverse momentum, also known as Drell-Yan (DY) process [1], has played a central role in the development of the Standard Model (SM) of the strong and electroweak (EW) interactions. The same process that allowed the discovery of the W and Z bosons [2][3][4][5] provides today very precise information about the gauge sector of the SM and allows its test at the quantum level. The determination of the masses of the gauge bosons and of the effective weak mixing angle, with relative errors at the 0.01% and 0.1% level respectively [6][7][8][9], is possible thanks to the excellent quality of the data collected at the Fermilab Tevatron by the CDF and D0 experiments and at the CERN LHC by the ATLAS, CMS, and LHCb experiments. The kinematical distributions are fitted by using theoretical templates, keeping the EW parameters as fitting variables. The uncertainty in the evaluation of the templates enters in the total error of the experimental value as a theoretical systematic contribution [10][11][12].
Hence, it is of utmost importance to minimise the latter in order to make the comparison more significant. The direct search for signals of New Physics is carried on, in the DY process, with the precise measurement of the tail of the kinematical distributions. A deviation could then be interpreted as the distortion due to a new resonance or, in general, as the effect of a new higher-dimensional operator whose increment with the energy shows up in the tails. In order to appreciate the possibility of such a deviation, the most precise SM prediction is needed, with the smallest possible residual theoretical uncertainty. The inclusion of higher-order corrections makes the theoretical predictions more precise, but it requires to overcome challenging mathematical and computational problems.
The radiative corrections to the DY processes in the strong coupling α s were computed at next-to-leading-order (NLO) [13] and next-to-next-to-leading-order (NNLO) [14,15] in perturbative Quantum Chromodynamics (QCD). Recently, the calculations at next-tonext-to-next-to-leading-order (N 3 LO) QCD of the inclusive production of a virtual photon [16], a W boson [17], and for the total neutral current (NC) DY process [18] have been completed. Fully differential NNLO QCD computations including the leptonic decay of the vector boson have been presented in refs. [19][20][21][22][23], and the first estimates of the fiducial cross sections for the NC DY process at N 3 LO QCD have appeared [24]. A plethora of papers discusses the inclusion of higher-order corrections in the threshold region [25][26][27][28][29][30][31][32][33]. The complete NLO EW corrections in the EW coupling α have been computed for W production in refs. [34][35][36][37][38], while the Z production case has been considered in refs. [39][40][41][42][43]. The NLO QCD and NLO EW corrections are separately large, often reaching the O(10 − 20%) of the leading order (LO) when we consider differential kinematical distributions, raising the suspect that also mixed NNLO QCD-EW effects might be potentially large. The fact that in specific phase-space corners NLO QCD and NLO EW corrections have opposite signs leads to large cancellations [44]. When such cancellations take place, a reliable estimate of the residual uncertainties based on the previous orders is problematic, since in these regions we might expect non-negligible effects from the next perturbative orders. Only an explicit calculation can solve the problem.
The combination of NLO QCD and EW results, consistently matched with QCD and QED Parton Showers, has been presented in [45][46][47][48]. These Monte Carlo tools provide an approximation of the mixed QCD-EW corrections, valid in the soft and/or collinear limits for the radiation of multiple additional partons, but can not assess the size of the genuine NNLO QCD-EW terms. Instead, they provide a rough estimate of the residual uncertainties.
Since high precision is required by the studies aiming at the determination of the EW parameters and the searches for New Physics, an important effort of the theory community has been recently devoted to the evaluation of the mixed QCD-EW corrections. The NNLO mixed QCD-QED corrections were obtained in ref. [49] for the inclusive production of an on-shell Z boson and extended at differential level for the decay of the Z boson into a pair of neutrinos [50]. The production of an on-shell Z boson, including NNLO mixed QCD-QED initial state corrections, but also the factorised NLO QCD corrections to Z production and the NLO QED corrections to the leptonic Z decay, has been discussed in ref. [51]. The complete NNLO QCD-EW corrections at O(αα s ) in analytical form for the total cross section production of an on-shell Z bosons have been presented in refs. [52][53][54][55], whereas results differential in the final-state leptonic variables have been presented in refs. [56] and [57] for the on-shell Z and W production respectively. Beyond the on-shell approximation, results have been obtained in the pole approximation [58], which offers the possibility to move beyond the strictly on-shell limit and is based on a systematic expansion around the W or Z resonance. By definition, its range of validity is expected to be restricted to this phase-space region. The dominant part of the mixed QCD-EW corrections at the W or Z resonances has been evaluated with this technique in refs. [59,60]. Precise predictions at NNLO QCD-EW level in the whole phase-space have been obtained in ref. [61] for the subset of O(n F α s α). The charged-current DY process pp → ν + X, has been considered in ref. [62] including the mixed QCD-EW corrections, with the reweighted two-loop virtual corrections in the pole approximation and all the other contributions in exact form. Very recently the complete set of NNLO QCD-EW corrections to the NC DY process pp → + − + X has been presented in ref. [63], with the inclusion of the exact two-loop virtual contributions 1 .
The virtual contributions were one of the bottlenecks for a complete O(α s α) calculation, because the evaluation of the 2 → 2 two-loop Feynman diagrams with internal masses is at the frontier of current computational techniques. Progress on the evaluation of the corresponding two-loop master integrals has been reported in refs. [66][67][68][69] and, recently, the computation of the two-loop helicity amplitudes for NC massless lepton pair production was discussed in ref. [70]. In this article, we provide in detail the independent computation at O(αα s ) of the two-loop amplitude for the process qq → + − . The calculation considers a massive leptonic final state, and features collinear logarithms of the lepton mass. These results, combined with the remaining perturbative contributions, allow the numerical evaluation of the hadron-level cross-section at this perturbative order, as presented in [63].
In Section 2, we provide the framework of our computation, with the description of the treatment of γ 5 and of the ultraviolet (UV) renormalisation procedure. In Section 3, the details of the computation are presented, from the amplitude level until the renormalised matrix elements. In Section 4, the infrared (IR) subtraction procedure is discussed. In Section 5 we describe how we obtain the UV renormalised IR subtracted finite remainder and we present the details necessary for the numerical evaluation of the analytical expressions. In Section 6 we draw our conclusions. The unrenormalised matrix elements in terms of master integrals, together with a benchmark grid for the finite reminder, are provided as ancillary material. Appendix A gives the reader the information needed to make use of the files.

The process
We consider two-loop mixed QCD-EW corrections to lepton pair production in the quark annihilation channel, as given by (2.1) We have separately computed the results for the uū and dd channels, with u and d denoting generic massless up-and down-type quarks. In the paper we present explicit results for the uū channel. The bare amplitude 2 for this partonic process admits a double perturbative expansion in the two coupling constants In this paper, we present the calculation of the following interference terms: which contribute at O(αα s ) to the unpolarised squared matrix element of the process in eq. (2.1). These scalar quantities contain many one-and two-loop Feynman integrals. We apply the integration-by-parts (IBP) [71,72] and Lorentz invariance (LI) [73] identities to reduce those scalar Feynman integrals to a limited set, the so-called master integrals (MIs). The corrections in eq. (2.3) can then be expressed as the sum of the MIs, each multiplied by its corresponding coefficient. The latter are rational functions of the kinematical invariants and of the masses which appear in the process. Virtual corrections are affected by singularities both of UV and IR type. The renormalisation procedure allows us to cancel the UV divergences and to express the amplitude in terms of renormalised parameters, which are in turn related to measurable quantities. The presence of IR divergences in the virtual corrections, on the other hand, is the counterpart of those which appear upon phase-space integration in the processes with the emission of additional massless partons, in our case a photon and/or a gluon. The universality of the IR divergent structure of the scattering amplitudes arising from soft and/or collinear radiation has been widely discussed in the literature and allows us to build a subtraction term, independently of the details of the full two-loop calculation. We exploit this property to check the consistency of our computation, in particular for the prescription related to the Dirac γ 5 matrix. To regularise both UV and IR divergences, we work in d = 4 − 2ε spacetime dimensions. The renormalised amplitude can thus be written as a Laurent expansion in ε. At one-loop level we have: while at two-loop we have: In eq. 2.6 we have retained only terms up to ε 0 , as we are only interested in the finite cotributions at NNLO. However, this series contains higher powers of ε which will be relevant for higher-order calculations. In these expressions m l is the lepton mass, and µ W , µ Z are the complex masses of the W and Z boson defined as: The mass and decay width M V , Γ V are real parameters and are the pole quantities, i.e. correspond to the position of the pole in the complex plane of the gauge boson propagator. We also introduce the Mandelstam variables, defined as follows: The on-shell conditions of the external particles are given by The coefficients C i in eq. (2.6) constitute the IR-unsubtracted two-loop interference terms, while those in eqs. (2.4,2.5), computed at one-loop, are necessary to prepare the subtraction term. The IR-subtracted expressions that can be obtained by their combination constitutes the main result that we present in this paper.

Treatment of γ 5 in Dimensional Regularisation
The calculation of radiative corrections, including chiral quantities in dimensional regularisation, faces the problem of defining a generalisation of the inherently four-dimensional object γ 5 in the arbitrary space-time dimension d. In [74], 't Hooft and Veltman proposed to abandon, in d dimensions, the anticommutation relation of γ 5 {γ µ , γ 5 } = 0 (2.10) and to keep the cyclicity of the Dirac traces; the price is a significant computational load, due to the different treatment of the four-momenta components defined in the four spacetime dimensions and of those in the remaining d−4. In refs. [75,76] Kreimer et al. retained the anticommutation relation, with a substantial computational gain, and abandoned the cyclicity property of the trace operation. It was demonstrated that, at least in one-loop calculations, the choice of a unique point to start the evaluation of all the traces needed in one calculation leads to the systematic cancellation of the prescription dependent terms, along with that of the IR poles. It has recently been explicitly checked in ref. [70], for the NC DY process at O(αα s ) , that the implementations at two-loop level of the 't Hooft-Veltman-Breitenloner-Mason and of the Kreimer prescriptions yield, as expected, different results for the scattering amplitudes, but are in perfect agreement when one considers the finite correction after subtraction of all the IR divergences.
In this computation, we retain the anticommutation relation together with and keep a fixed point to write and then evaluate the Dirac traces. As a result, we end up with zero or one γ 5 at the rightmost position of the Dirac trace. In the single γ 5 case, we perform the replacement where, ε µνρσ is the Levi-Civita tensor. The contraction of two such tensors yields and all the Lorentz indices are considered to be d-dimensional.
The presence of a prescription-dependent term of O(ε) in the squared matrix element affects all the coefficients in the Laurent expansion, with the exception of the highest pole: in fact, the product of a term of O(ε) with a singular factor ε −k , with k > 0, generates a contribution of O(ε −k+1 ). Such prescription-dependent terms will be generated both in the unsubtracted squared matrix element and in the subtraction term. The cancellation of the IR singularities, expected on general grounds, requires that also the prescription-dependent terms cancel accordingly. In the present calculation, the IR subtraction term is computed by following the properties of universality of the radiation in the IR limits, combining the universal divergent structure with the Born and one-loop amplitudes. The construction of this subtraction term is completely independent with respect to the evaluation of the two-loop amplitude and it provides a non-trivial check of our algebraic manipulations. We observe the cancellation of all the lower order poles, when combining the full twoloop amplitude with the subtraction term, which hints in favour of the consistency of our approach.

Ultraviolet renormalisation
The renormalisation at O(αα s ) of the NC DY process has already been discussed in detail in ref. [61]. We report here the basic steps that we have implemented to obtain the complete two-loop renormalised amplitude.

Charge renormalisation
The bare gauge couplings g 0 , g 0 and the Higgs doublet vacuum expectation value v 0 are expressed in terms of their renormalised counterparts g, g , v via the introduction of appropriate counterterms. The relation of g, g , v to a set of three measurable quantities, like for instance G µ , µ W , µ Z (with G µ the Fermi constant) or α, µ W , µ Z (with α the fine structure constant) 3 , eventually allows the numerical evaluation of the amplitude. We define for convenience two additional bare quantities: the sine squared of the on-shell weak mixing angle, which we abbreviate as , and the electric charge e 0 = g 0 s W 0 , but we stress that only three of the above parameters are independent. We introduce the relation between bare and renormalised input parameters The on-shell electric charge counterterm δe at O(α 2 ) and O(αα s ) has been discussed in ref. [78], from the study of the Thomson scattering. The mass counterterms δµ 2 W , δµ 2 Z in the complex mass scheme [79] have been presented in ref. [61]. In terms of the transverse part of the unrenormalised V V gauge boson self-energies, they are defined as follows: at the pole in the complex plane q 2 = µ 2 V of the gauge boson propagator, with V = W, Z. From the study of the muon-decay amplitude, we derive the following relation The finite correction ∆r was introduced, with real gauge boson masses, in ref. [80] and its O(αα s ) corrections were presented in ref. [81,82]. We evaluate it here with complex-valued masses.
We consider now the bare couplings which appear at tree-level in the interaction of the photon and Z boson with fermions. The UV divergent correction factors δg α Z contributes to the charge renormalisation of the Zff vertex in the (α, µ W , µ Z ) input scheme Working out the explicit expression of δg Gµ Z , the dependence on the electric charge counterterm cancels out. Analogously, in the case of the γff vertex, the electric charge renormalisation is given by in the (α, µ W , µ Z ) scheme and by in the (G µ , µ W , µ Z ) scheme. The counterterm contributions to the renormalised amplitude are obtained by replacing the bare couplings in the lower order amplitudes with the expressions presented in eqs. (2.18-2.21) and expanding δg A,Z up to the relevant perturbative order. We show in the next Section how δg A,Z enter in the renormalisation of the gauge boson propagators.

Renormalisation of the gauge boson propagators
The renormalised gauge boson self-energies are obtained, at O(α), by combining the unrenormalised self-energy expressions with the mass and wave function counterterms. In the full calculation, we never introduce wave function counterterms on the internal lines, because they would systematically cancel against a corresponding factor stemming from the definition of the renormalised vertices. We exploit instead the relation in the SM between the wave function and charge counterterms [58] and we directly use the latter to define the renormalised self-energies. We obtain:

Checks
The physical scattering amplitude does not depend on the gauge choice used to quantise the theory. In our computation we use the background field gauge (BFG) [83], which restores some QED-like Ward identities in the full EW SM and guarantees in turn some properties of the Green's function which appear in the calculation. For instance, the UV finiteness of the vertex corrections, summed with the external wave function of the fermions entering in that vertex, is an important property which can provide an additional check of the calculation. This check has been performed explicitly at the one-loop level, based on the knowledge of the UV poles of the one-loop Feynman integrals. At two-loop level, we consider the UV finiteness of the box corrections, we explicitly renormalise the tree-level propagators, and we assume that the cancellation of the UV poles between vertex and external wave function corrections takes place. Since the propagator corrections are IR finite, we conclude that the sum of vertex, box, and external wave function corrections contains only IR poles. We have explicitly checked that the coefficients of these poles exactly match the ones predicted by the universal structure of the IR divergences in gauge theories [84]. This result hints in favour of the validity of the Ward identity also at two-loop level.
3 The UV-renormalised unsubtracted amplitude (a) include an initial-state QCD vertex correction; the second factor can be, alternatively: the final-state EW vertex corrections, the external lepton wave function correction, the one-loop self-energy corrections, the one-loop mass and charge renormalization counterterms, and the one-loop external quark wave function correction. They are schematically represented in Figures 2(a)-(f) respectively. Also in this case the BFG properties allow to identify UV-finite combinations of Feynman diagrams. In the interference term M (0) |M (1,1) we recognise different subsets associated to different combinations of the EW bosons exchanged in the loops. We consider the following cases: γγ, γZ, ZZ, W ; in the neutral cases there are either one boson in a loop and one in a tree-level line or both bosons in the loops, as depicted in a few representative cases in Figure 3 a-c for the γZ subset; in the charged case we find one or two W s always in the loop lines Figure 3 d-f. The calculation of the bare matrix elements up to two-loop follows a general procedure. For the generation of the Feynman diagrams we have used two completely independent approaches, one based on QGRAF [85] and a second one on FeynArts [86]. Two independent in-house sets of routines, written in FORM [87] and Mathematica [88] have been used to perform the Lorentz and Dirac algebra. After these first manipulations, the intermediate expressions contain a large number of scalar Feynman integrals, which can be reduced to a much smaller set of MIs by means of IBP and LI identities. We have executed the reduction algorithms with two independent in-house programs based on the public codes Kira [89], LiteRed [90,91] and Reduze2 [92,93]. We have cross-checked the corresponding expressions at one-and two-loop level, obtained in the two independent approaches, finding perfect agreement. The evaluation of the quark wave function corrections (Figure 1-b) deserves a comment because it proceeds in two steps. First, we generate all the Feynman diagrams contributing to the fermion self-energy and reduce their Feynman integrals to MIs. Second, since we have to compute the derivative of the self-energy with respect to the external invariant, we express the derivative of the MIs as a combination of MIs with a standard reduction procedure. We eventually evaluate this combination in the on-shell limit.
To simplify the calculation, we consider an approximation of the scattering amplitude in the small lepton mass limit, i.e. when m is negligible compared to the masses of the gauge bosons and to the energy scales of the process. We thus consider the ratio m / √ s and keep only terms enhanced by a factor log(m 2 /s), while we discard all those contributions vanishing in the m → 0 limit. This approximation is used in the reduction to MIs procedure, in order to distinguish the integrals where the dependence on m has to be considered from those where the massless lepton limit can be immediately taken, with a remarkable simplification in the latter case. After such procedure, the only contributions with an explicit leftover dependence on the lepton mass are the ones stemming from factorisable contributions with initial-state QCD and one-loop final-state vertex correc-tions. Among the box diagrams, those with one or two photons exchanges individually develop logarithms of the lepton mass, but the latter eventually cancel [94] in the physical amplitude.
The identification of all the MIs which contribute to the final result is given using the notion of integral family, which we present in the next Section.

Integral families and the representation of the result in terms of MIs
We use the word topology to identify the flow of the external momenta in the internal lines of a Feynman integral. An integral family is defined as a complete set of denominator factors, sufficient to reduce all the scalar products present in one scalar Feynman integral in terms of linear combinations of these quantities. In general, an integral family includes more denominator factors than those appearing in a specific Feynman diagram. In the present calculation we discuss two-loop four-point topologies and we find that the corresponding integral families contain 9 denominator factors, i.e. two more than the maximum number of internal lines of the Feynman diagram. The Feynman integrals with a smaller number of internal lines are called subtopologies, and belong to the same integral family of the integral under study.
By associating a scalar integral to a specific integral family, we can express all the scalar products in terms of the denominators of that particular family. We thus obtain, after all the algebraic simplifications, integrals with the following form: where the exponents α k can be either positive (0, 1 or 2) or negative. We introduce the symbols Ds and we use them to precisely define the various integral families. Those relevant for the full calculation are listed here below.
In the cases where the lepton masses are required, because of the presence of final-state collinear singularities, we find the following integral families which accommodate all the appearing scalar integrals: Additional families, where the lepton masses can be immediately discarded, are: where V can be either Z or W . Each MI is thus uniquely identified by its integral family and the exponents α j , e.g. B 20 [1, 1, 1, 0, 1, 1, 0, 1, 0] is the symbolic form of .

(3.5)
In the following, we list the appearing MIs for all the integral families. The massless two-loop form factor integrals were presented in [95]: The form factor type two-loop MIs with one or two massive lines were obtained in [98,99]: The two-loop box MIs with one and two massive lines have been computed in [66][67][68]. We relate our MIs as given below to the ones obtained in [66] through the IBP reduction and obtain the results for our MIs:  (Figures 1-a, 1-f), in terms of the MIs with the corresponding rational coefficients schematically described in eq. (3.10), as a useful intermediate representation.

Infrared Singularities and Universal Pole Structure
The renormalised matrix elements are UV-finite, but still contain IR divergences originating from soft and/or collinear massless partons. These singularities are guaranteed to cancel when combined with the NNLO QCD-EW real and real-virtual corrections and subtracted of the initial-state collinear poles, to obtain an IR-safe observable. Nevertheless, their presence at intermediate stages of an inclusive computation, as it is the case with the results provided in this paper, requires the development of methods to consistently handle them; they must be systematically subtracted, leaving a finite remnant. While at NLO computations the Catani-Seymour dipole subtraction [100][101][102] and FKS subtraction [103] are the two methods most widely used, at NNLO several techniques have been proposed (see e.g. [104] and references therein), each one with a different subtraction procedure for the various contributions to physical observables.
The subtraction of the IR poles from the two-loop matrix elements is in general achieved via a process-independent subtraction operator I, that can be constructed by using the universality of the IR singularity structure of the scattering amplitudes, known in the case of a massless gauge theory up to two-loop level [84,[105][106][107]]. An explicit study was performed in [108,109] to obtain the IR structure of the two-loop amplitudes for mixed QCD⊗QED corrections to NC Drell-Yan production considering massless leptons. The case of theories with massive particles has been studied in [110][111][112][113][114]. In particular, the IR structure for one-loop QCD corrections to top quark pair production has been examined in detail in ref. [115], while at two-loop level in refs. [116][117][118]. This can be appropriately abelianised [119] to obtain the IR structure in the present case by replacing the top quarks with massive leptons.
Different subtraction operator within different subtraction schemes can in principle differ from each other: anyhow the universality of the IR divergences implies that the subtracted virtual contribution at two loops may differ, at most, by a finite constant. For the sake of definiteness, we illustrate here the subtraction operator used in our calculation, which has been performed in the framework of the q T -subtraction formalism [120]. However, we stress that our results can be easily converted and combined to the description of real radiation in any other IR subtraction scheme, by adding the appropriate finite remainder.
The IR subtraction functions (I) at one-loop are given by Q l and Q u are the charges of the lepton and of the initial-state quark, and the Casimir of the fundamental representation of SU(N), C F , is given by Using the one-loop subtraction functions, we obtain the finite contributions to the one-loop QCD and EW amplitudes, respectively, as follows: The mixed two-loop subtraction operator 5 is given by In order to describe the contributions entering in I (1,1) , let us discuss first the IR singularities of the one-loop and two-loop boxes. The fact that the qq initial state has zero electric charge implies that there are no QED collinear singularities stemming from the exchange of a photon between an initial quark and a final lepton lines. Such terms appear in the individual diagrams, but cancel in the physical amplitude. This conclusion can be also seen as a confirmation of the universality of initial state QED collinear singularities. This property is observed at one loop and confirmed also when going to O(αα s ) by dressing the initial state quark line with all the possible gluon exchanges. In the sum of all the two-loop box diagrams we can thus expect that the photon exchange develops at most a soft divergence. The presence of a colorful initial state quark line but of a colorless final state lepton line leads to a natural separation between the gluon IR-divergent exchanges in the initial state and the QED interaction in the rest of the diagram. For this reason, all the IR singularities, soft and collinear, stemming from the two-loop box diagrams, can be subtracted by the product of I (1,0) and I (0,1) , times the Born amplitude. The divergences in the product of initial and final factorisable contributions are obviously canceled by the same product of one-loop subtraction operators. The initial state vertex corrections have a richer structure, including planar and non-planar topologies, both yielding IR divergent contributions. The former can be subtracted again by the product of I (1,0) and I (0,1) , while the latter require a genuinely two-loop subtraction operator. From the above comments we conclude that the operator I (1,1) is applied to the Born amplitude and receives contributions from the product of I (1,0) and I (0,1) and from the two-loop non-factorisable contributions appearing in the vertex corrections. We note that the application of I (1,1) can not make the complete amplitude IR finite. There are in fact additional configurations where e.g. the gluon yields a divergence and the EW interaction is finite. This is for instance the case in the two-loop boxes with the exchange of two W or Z bosons. They can be subtracted by applying the one-loop QCD subtraction operator I (1,0) to the one-loop EW amplitude subtracted of its one-loop QED divergences, as defined in eq. (4.6). An analogous subtraction can remove the photon divergences which multiply a finite QCD remnant, as defined in eq. (4.5).
Using eqs. (4.5-4.7), we obtain the finite and subtracted two-loop amplitude as follows withĨ (i,j) defined by dropping the term proportional to ζ 2 in I (i,j) . This choice is conventional and defines the finite part of our subtraction term. The approximation of the amplitude in the small lepton mass limit retains all the terms enhanced by log(m l ), divergent in the m → 0 limit. The structure of these corrections reflects the universality property of the final-state collinear divergences, and is given, normalised to the Born squared matrix element, by lim m →0 where K represents all the other terms in the interference, constant in the m → 0 limit.
The final results, expanded in powers of , are written in terms of GHPLs and of the symbols associated to the 5 MIs whose solution is given in terms of Chen integrals. In the next Sections, we explicitly discuss the numerical evaluation of the latter. We note that the individual contributions from these 5 MIs, to the single pole in of the full unrenormalised matrix element, contain the Chen iterated integrals, but the summed result is independent of those.

Semi-analytical solution of MIs via series expansion
In ref. [66], 5 MIs with two internal massive lines have been solved in terms of Chen iterated integrals, but their analytical continuation to the physical region is extremely difficult. We thus turn to a different approach and we obtain these MIs in semi-analytical form, by solving the relevant system of differential equations via series expansion as initially proposed in ref. [131]. This approach is implemented in the Mathematica code DiffExp [132]. However, since the public version of DiffExp does not support, for the time being, the possibility of complex masses in the internal lines of the loop integrals, we have written an independent implementation of this procedure, that consistently takes into account the presence of complex weights and the analytic continuation to the complex plane [133].
By approaching the problem in such semi-analytical way, we have to consider the complete system (36×36) of differential equations [66], which includes the 5 MIs of eq.  predictions. In the euclidean region, we have checked our results for these 5 MIs against the values computed with PySecDec [135], with perfect agreement within the numerical precision of PySecDec. In the physical region such comparison can not be performed, because of convergence issues observed with PySecDec, and we have adopted a different strategy. We have solved the same system of differential equations with a real value for the gauge-boson masses, by using both DiffExp and our private implementation, finding perfect agreement 8 , for different values of the centre-of-mass energy, below and above all the physical thresholds.
In Figure 4 we present the real and imaginary parts of the values of the 0 coefficient of the MIs numbered from 32 to 36 in ref. [66], in a phase-space region relevant for LHC phenomenology. As a technical benchmark we provide in eqs. (5.4) the values of these five integrals, at √ s = 120 GeV and cos θ = 0, using the Z boson complex mass with the values presented in the next section. This approach faces all the problems of the analytical continuation immediately from the analysis of the structure of the differential equations. If it is possible to expand the solution in a region around the boundary condition points, then the problem becomes of algebraic nature, i.e. we have to solve for all the coefficients of the various powers of the expansion variable. The system of differential equations of ref. [66] can be cast in dlog form and thus contains the complete information about the singularity structure of the problem, in particular the position of the branch points relevant in the discussion of the polylogarithmic functions present in the solution. This information can be exploited to control the analytic continuation when we extend the solution to an adjacent region, either non-physical or physical. We remark that complete control on the analytic continuation has been verified when letting vary one kinematical variable at the time, while the others are kept fixed.

Numerical results
In this Section, we present the numerical evaluation of our result, the finite IR-subtracted UV-renormalised hard function H (1,1) . For every practical application, the evaluation in one phase-space point of the correction factor stemming from the two-loop O(αα s ) virtual corrections requires a cumbersome procedure. We consider it more convenient to prepare a numerical grid for H (1,1) , covering all the phase-space values relevant for NC DY in a given fiducial volume and parameterised in terms of the partonic centre-of-mass energy √ s and scattering angle cos θ. We notice that the virtual correction, after the subtraction of all the IR enhanced factors, is a smooth, slowly varying function of √ s and cos θ. We have prepared a grid with a sampling based on the known behaviour of the NLO-EW distribution, with special care for the Z resonance region, where a finer binning is necessary. For illustration purposes we consider here a grid with (130x25) points in ( √ s, cos θ), covering the intervals s ∈ [40, 13000] GeV and cos θ ∈ [−1, 1] and present the values in the ancillary file gridsth.m . The extraction of any intermediate value of H (1,1) can be obtained via interpolation of the grid points with excellent accuracy thanks to the smoothness of the function. With the sampling provided as ancillary file, we reproduce the exact correction with an accuracy at the 10 −3 level, in the whole phase space. This accuracy is sufficient for phenomenological studies, but can be increased whenever needed.
The preparation, on a single core, of a grid of values for all the 36 two-masses MIs requires a fixed step, executed only once, to move from the boundary conditions fixed in the non-physical euclidean region to the physical region, which can take ∼ 20 minutes. Moving in the physical region is faster, with an average time of ∼ 5 minutes per point. Once this first grid is ready, the evaluation of the corresponding grid for the H (1,1) function requires ∼ 3 minutes per point. Once the H (1,1) grid is ready, its interpolation in any simulation code requires a negligible time.
We use the following parameters: The values of the gauge boson masses and decay widths are compatible with those reported in the PDG [136], extracted in a running-width fitting scheme. The choice of the complex mass scheme to define the masses of the unstable W and Z bosons implies that all the functions which depend on these parameters become functions of complex variables. We have exploited the rescaling properties of the GHPLs to set their argument equal to one, adjusting all the other weights, which become in general complex 9 . For the numerical evaluation of all the GHPLs appearing in the final expressions, we use GiNaC and handyG [138] in two independent C++ and Mathematica codes to evaluate each phase-space point.
We also use HarmonicSums [139,140], PolyLogTools [141] and LoopTools [142] for several cross-checks of the one-and two-loop integrals available in closed analytic form, with real and complex masses. For the numerical evaluation of the 5 MIs of eq. (5.3), we use the semi-analytical approach described in Section 5.2. We present in Figure 5 the numerical grid expressing the UV-renormalised IR-subtracted two-loop O(αα s ) virtual correction due to the two-loop vertex and box diagrams ( Figure  1 a,b,f and Figure 2 a,b,f). The plots show the correction factor normalised to the Born cross section and expressed in units α π αs π , as a function of √ s and cos θ. We consider the range −1 ≤ cos θ ≤ 1 and, for the partonic center-of-mass energy, two different intervals, namely 50 ≤ √ s ≤ 250 GeV and 50 ≤ √ s ≤ 5000 GeV. With the zoom in the 50 ≤ √ s ≤ 250 GeV interval, we can appreciate the non-trivial structure of the corrections starting from the Z resonance and up to 110 GeV. The description of the W W and ZZ thresholds is smooth, thanks to the implementation of the complex-mass scheme. Above 200 GeV the shape of the correction factor is negative, smooth, and slowly decreasing; we do not observe the appearance of additional structures. The non-trivial cos θ dependence is due to parity violating effects, in largest size stemming from the interference between the two-loop AA subset of diagrams and the tree-level amplitude with a Z-boson exchange.

Conclusions
In the NLO case, a complete automation of all the steps of the calculation is in a quite mature stage for both QCD and EW corrections. In the NNLO case instead, only parts of the procedure are automatic, because of various conceptual and computational problems. In this paper we have presented the details of the complete calculation of the exact O(αα s ) virtual corrections to the NC DY process, describing how we have overcome the problems appearing in the different stages of the procedure that brings us from the Feynman diagrams to the numerical evaluation of the squared matrix element.
All the algebraic manipulation processes face the difficulty of dealing with potentially very large expressions, so that a systematic simplification work, based on the search for recurring patterns, is crucial to keep the size of the formulae at a manageable level. While the reduction to MIs is completely standardised, the analytical subtraction of the IR divergences is not mature at the same level, and we have achieved this successful check only via a dedicated calculation. The presence of internal massive lines yields a rich analytical structure of the Feynman integrals; for the evaluation of the two-loop boxes with two massive lines we have adopted a semi-analytical approach, valid for arbitrary complex masses.
The final correction factor, defined in the q T -subtraction scheme, has been cast in a relatively compact form, which we provide in the ancillary files. This finite virtual factor can be used in any simulation code, provided that the appropriate finite correction for the IR subtraction scheme is applied. For benchmarking we provide the grid with the numerical evaluation of the correction factor in the phase-space region relevant at the LHC.
The successful completion of all these steps can help to devise new strategies for the automation of the calculation of two-loop EW and mixed QCD-EW radiative corrections, relevant for processes at the LHC and future lepton colliders.
of the NNLO QCD-EW corrections to NC DY and for a careful reading of the manuscript. We would like to thank G. Heinrich for help with pySecDec and R. Lee for help with

LiteRed.
S.D., N.R. and A.V. are supported by the Italian Ministero della Università e della Ricerca (grant PRIN201719AVICI 01). R. B. is partly supported by the italian Ministero della Università e della Ricerca (MIUR) under grant PRIN 20172LNEEZ.

A Description of the ancillary files
In this subsection we present a description of the content of the ancillary files attached to the paper. A complete listing of all the symbols used in the analytical expressions and their correspondence with the physical quantities can be found in the README file.
In the ancillary Mathematica file NCDY ME11.m, we present M (0) |M (1,1),f in i.e. the finite IR subtracted and UV renormalised matrix elements of all the vertex and box type Feynman diagrams. As discussed earlier, we present the result in the small lepton mass limit i.e. dropping all the terms which vanish as m l → 0. In the ancillary file uU totaldelta vertbox.dat, we provide the numerical evaluation, for up-type quark initiated channel, of the expressions contained in NCDY ME11.m, divided by Born, multiplied by 2, divided by 16 and also subtracted of the logarithms of the lepton mass, using the numerical value of the input parameters as described in Section 5.3, for all the points of the grid specified in gridsth.m. We also provide the exact result including the full dependence of the lepton mass in ME11 partml exact.m for the following subsets: all the vertex and box type contributions for the γγ subset and the contributions from the Feynman diagrams with NLO QED correction to the lepton vertex in the γZ subset. The explicit expressions of the two-loop self energy and one-loop self energy with NLO QCD contributions, are known in the literature (cfr. Refs. [61] and Refs. [143]). Hence, we refrain ourselves from providing these results.
The results are written as a combination of GHPLs and of the symbols associated to the MIs whose solution has been obtained via a semi-analytical method. Several variables have been introduced to define the weights and arguments of the GHPLs, which we report here for the sake of definiteness. Extending the notation of eq. (4.4), the Landau variables x m are defined with respect to the mass m as follows, where m ∈ {m l , µ Z , µ W }. The kinematical ratios are defined as The presentation of the variables appearing in the expression is as follows: xrl =x m l ; xl = x m l ; yl = y m l ; zl = z m l ; xrZ =x µ Z ; xZ = x µ Z ; yZ = y µ Z ; zZ = z µ Z ; vZ = v µ Z ; vpZ = v µ Z ; µ Z is the complex conjugate of µ Z .