NNLO QCD$\times$EW corrections to Z production in the $q\bar{q}$ channel

We present the first results for the ${\cal O}(\alpha\alpha_s)$ corrections to the total partonic cross section of the process $q\bar q\to Z+X$, with the complete set of contributions, that include photonic and massive weak gauge boson effects. The results are relevant for the precise determination of the hadronic $Z$ boson production cross section. Virtual and real corrections are calculated analytically using the reduction to the master integrals and their evaluation through differential equations. Real corrections are dealt with using the reverse-unitarity method. They require the evaluation of a new set of two-loop master integrals, with up to three internal massive lines. In particular, three of them are expressed in terms of elliptic integrals. We verify the absence, at this perturbative order, of initial state mass singularities proportional to a weak massive virtual correction to the quark-gluon splitting.

The production of an electrically neutral gauge boson at hadron colliders is one of the historical processes for our understanding of Quantum Chromodynamics (QCD). The case of the decay of the Z boson into a pair of high transverse momentum leptons is known as Drell-Yan (DY) process and it is particularly important for the setting of several high-precision tests of the electroweak (EW) sector of the Standard Model (SM). It allows for instance a precise measurement of the weak mixing angle and of the properties of the Z boson. The Z boson DY production is one of the processes known with high perturbative accuracy. The pioneering calculations of the next-to-leading order (NLO) [1] and nextto-next-to-leading order (NNLO) [2] QCD corrections to the total inclusive cross section have been extended later to the fully differential description of the leptonic final state [3][4][5][6]. Finally, the evaluation of the next-to-nextto-next-to-leading order (N 3 LO) QCD corrections at the production threshold has been presented in Refs. [7][8][9][10] for the gauge boson total cross section and rapidity distribution. The impact of the NLO EW corrections, studied in Refs. [11][12][13][14], is at the O(1%) level as far as the total cross section is concerned and it is comparable to that of the NNLO QCD contributions. Kinematic distributions may receive additional enhancements in specific phasespace regions, yielding corrections at the O(10%) level or more. Since the high-precision determination of EW parameters requires control over the kinematic distributions in some cases at the per mille level (cf. Refs. [15][16][17] for a discussion on specific examples), the evaluation of the mixed QCD-EW corrections has emerged as necessary for both the study of the gauge boson resonances and of the high mass/momentum tails of the kinematic distributions [18,19]. First analytic results have been presented in Refs. [20][21][22][23][24], and compared with the approximations available via Monte Carlo simulation tools [25,26]: while the bulk of the leading effects, separately due to QCD and QED corrections, can be correctly evaluated for several observables, the remaining sub-leading QED effects and the genuine QCD-weak corrections are still missing in these tools. Furthermore, a realistic estimate of the theoretical uncertainties must account for several sources of ambiguity related to the recipes used in the matching of separate results for the QCD and EW contributions to the scattering amplitude. For these reasons an exact calculation of the full set of O(αα s ) corrections to the DY processes is desirable. In Refs. [27][28][29] the mixed QCD-QED corrections to the total cross section and transverse momentum spectrum of an on-shell Z boson have been discussed. The evaluation of all the Master Integrals (MIs) relevant to compute the full set of QCD-EW mixed corrections to DY process (including offresonance terms) has been documented in Refs. [30,31].
In this letter, we present the first results for the total inclusive cross section of production of an on-shell Z boson in the quark-antiquark partonic channel, including the complete set of QCD-EW corrections of O(αα s ). We retain the dependence on the massive states exchanged in the loops. As a consequence of that, the calculation involves a set of two-loop phase-space integrals, previously not available in the literature. Their analytic expression will be presented in a forthcoming paper. We also have the occasion to check the infrared structure of the corrections up to NNLO level, including the cases where a massive EW boson is exchanged. We verify the absence of initial state mass singularities proportional to a weak massive virtual correction to the quark-gluon splitting.
The calculation we are presenting in this letter is an important step towards the evaluation of the full set of QCD-EW corrections to the hadronic cross section. The inclusive production cross section σ tot of a Z boson at hadron colliders (pp → Z + X) can be written, using the factorization theorem, as s are the ratio of the squared Z boson mass, m Z , with S andŝ, the hadronic and partonic center of mass energy squared, respectively. S andŝ are related byŝ = x 1 x 2 S through the Bjorken momentum fractions x 1 , x 2 . The bare cross sectionσ ij of the partonic process ij → Z + X is convoluted with the bare parton densitiesf i (x). The sum over i, j includes quarks (q), antiquarks (q), gluons (g) and photons (γ). In the SM, we have a double expansion of the partonic cross sections in the electromagnetic and strong coupling constants, α and α s , respectively: whereσ (m,n) ij is the correction of O(α m s α n ) to the lowestorder inclusive total cross sectionσ (0,0) ij of the partonic scattering ij → Z. For a given initial state, the inclusive total cross section receives contributions from processes with different final state multiplicities, due to real parton emissions. In this letter we focus on the qq initiated scattering, and, for definiteness, we treat the case of an up-type quark: qq = uū. The full set of O(αα s ) corrections toσ uū stems from the evaluation of the following scattering processes: where d represents a down-type massless quark. Explicit expressions for the process (6) and QCD-QED contributions to process (3) have been presented in Ref. [32] and Ref. [20,33], respectively. The corresponding results for dd initiated subprocesses can be derived from our results with the replacements Q u ↔ Q d , I are the electric charge and the third component of the weak isospin, for a fermion f , respectively.
The process (3) receives contributions from two-loop 2 → 1 Feynman diagrams, that have to be interfered with the tree-level uū → Z and constitute the virtual corrections. The processes (4)-(5) receive contributions from one-loop 2 → 2 Feynman diagrams that have to be interfered with the corresponding tree-level. We refer to them as real-virtual corrections. The last three processes, (6)-(8), receive contributions from tree-level 2 → 3 Feynman diagrams interfered with themselves and we refer to them as double-real corrections.
The full set of O(αα s ) corrections can be organized in two gauge invariant subsets: QCD-QED and QCDweak contributions. Processes (3)-(7) contribute to the former, and processes (3), (4), (7) and (8) to the latter. While one gluon exchange, real or virtual, is always present, we identify three groups of contributions to the amplitudes depending on the presence of one real or virtual photon, of one virtual Z boson, or of one/two virtual W bosons. We further observe that the last two groups are separately gauge invariant. In our definition of total cross section we do not include the processes with the emission of one extra massive on-shell gauge boson, as their measurement depends on the details of the experimental event selection. Furthermore, these corrections do not contribute to the infrared structure of the process.
The amplitude of the two tree-level processes (7) and (8) has two components of O( √ αα s ) (an internal gluon exchange) and O( √ αα) (an internal weak boson exchange), respectively and their interference is, therefore, of O(α 2 α s ).

COMPUTATIONAL DETAILS
We follow a diagrammatic approach to obtain all the relevant contributions to the inclusive production cross section uū → Z + X. A detailed description of the computation will be presented in a dedicated publication. In this letter we sketch an outline of the procedure. We need to include contributions with two-loop virtual corrections, with one real emission and one loop (realvirtual), with two real emissions (double-real), and factorisable contributions stemming e.g. from the interference of two one-loop diagrams. We treat all the processes with the same algorithmic approach. Firstly, we compute all the Feynman diagrams contributing to a given amplitude with FeynArts [34] and QGRAF [35], we perform algebraic simplifications with FORM [36] and Mathematica; we use integration-by-parts (IBP) [37][38][39] and Lorentzinvariance (LI) identities [40] to reduce the Feynman integrals to MIs. The reduction to the MIs is carried out using the computer programs Kira [41], LiteRed [42,43], and Reduze 2 [44,45] The entire procedure is performed within dimensional regularization in D = 4 − 2ε spacetime dimensions. Then, we employ the method of differential equations [40,[46][47][48][49] to compute the MIs, for both the pure virtual and real emission corrections. In the latter case, the phase-space delta functions are dealt with via the reverse unitarity technique [50,51], which is based on the observation that the replacement known as Cutkowsky rule holds in terms of distributions: (9) It is thus possible to rewrite the phase-space measure of each final-state particle as the difference of two propagators with opposite prescriptions for their imaginary part (where η stands for an infinitesimal positive real number). We transform the integration over the full phase space of the additional parton/s for processes (4)- (8), into the evaluation of the cut two-loop integrals with an on-shell condition on the lines that correspond to the final-state particles.
The pure virtual MIs are already available in the literature [52][53][54][55][56][57], in the case of off-shell Z boson. Since in our case the Z boson is on-shell, we have computed these integrals taking the appropriate on-shell limit. The two-and three-body phase-space MIs with only gluon or photon lines are already available in the literature [51]. To validate our routines developed for the present calculation, however, we have recomputed them and found complete agreement with the known expressions. We have computed all the new MIs, with one or two internal massive lines, with the differential equations method. We have fixed the boundary conditions calculating the soft limit (z → 1) of the MIs.
After integration over the phase-space of the emitted real partons, the partonic total cross section depends solely on the variable z. The virtual contributions are proportional to δ(1 − z) and are therefore constants, which are found from the on-shell limit of the virtual MIs, i.e. evaluating the corresponding generalised harmonic polylogarithms (GPLs) [58][59][60][61] at z = 1. All the constants arising from this limit can be reduced to the basis introduced in Ref. [62]. The part that corresponds to processes (4)-(8) is expressed almost entirely in terms of δ(1 − z) and of GPLs, or cyclotomic Harmonic Polylogarithms [63], functions of z. In some parts of the calculation the package HarmonicSums [64] has been used. Three MIs appearing in processes (7) and (8) satisfy elliptic differential equations, whose homogeneous behaviour has already been studied in Ref. [55]. We have obtained their complete solution with a series expansion around z = 1 (see for instance [55,[65][66][67][68]).
In the calculation of the MIs, the masses of the W and Z bosons are set equal to m Z , to avoid the presence of an additional energy scale in the problem, which would make the analytical solution of the differential equations in terms of known functions more complicated. While for the virtual corrections this choice is not strictly necessary, since the knowledge of the MIs for off-shell Z would allow for a complete and exact calculation in the case of m W = m Z , the reduction of one mass scale in the computation of the real emission processes is in fact very effective and reduces the complication of the calcu-lation. Moreover, the equal-mass choice does not prevent us from obtaining an analytical solution with arbitrary precision for each of the affected MIs. In fact, we can perform an expansion of the integrand in powers of the ratio δ 2 m = (m 2 Z − m 2 W )/m 2 Z , and reduce all the terms of the series to a combination of the same basic equalmass MIs. We stress that the couplings of the Z boson to fermions are expressed in terms of the physical value of the weak mixing angle sin 2

Ultraviolet renormalization
The calculation is performed in the EW background field gauge (BFG) [69], which allows the identification of two sets of ultraviolet (UV) finite amplitudes. On the one hand, the combination of 1PI vertex and external quark wave function corrections, which satisfies, also in the EW SM, a QED-like Ward identity, with the consequent cancellation of the UV poles. On the other, the external Z boson wave function and the lowest-order coupling renormalization corrections, whose combination, orderby-order in perturbation theory, is also UV finite.
We need to perform the renormalization of the couplings and the fields up to O(αα s ) for process (3), while we need only the O(α) renormalization of process (4). One-loop QCD corrections to processes (3) and (5) are UV finite, after field renormalization, again because of a QED-like Ward identity. We remark that the Z boson field and the EW couplings do not receive O(α s ) renormalization corrections. The renormalization of the quark field receives EW corrections and we consider this in the on-shell scheme. The EW gauge sector of the SM Lagrangian depends on three parameters (g, g , v), the two gauge couplings and the Higgs-doublet vacuum expectation value. After the introduction of counterterms and renormalized parameters, we express the latter as a combination of (G µ , m W , m Z ) [73], respectively the Fermi constant, the W and Z boson masses.
A subset of the EW corrections can be reabsorbed in a redefinition of the weak mixing angle that appears in the vector coupling of the Z boson to fermions. These corrections are split, in the EW BFG, in two UV-finite groups, one due to vertex corrections, the other due to external γ − Z corrections and to the weak mixing angle counterterm (a shortcut for a combination of W and Z mass counterterms). In BFG the second group vanishes, because of a Ward identity [69] satisfied by the γ − Z wave-function correction.

Infrared singularities and mass factorization
The O(αα s ) corrections are organized in two gauge invariant subsets: QCD-QED and QCD-weak contributions. The former involve the exchange of two massless bosons, yielding the maximal degree of infrared singularity at the second perturbative order, i.e. ε −4 . The latter have only the poles due to a soft and/or collinear gluon. The cancellation of the soft singularities takes place separately in the two subsets, once the contribution of virtual corrections and of the corresponding soft real emissions are combined. To be more precise, for the QCD-QED subset, the process (7) does not yield soft singularities, so that the cancellation takes place when the processes (3)-(6) are combined. In the case of the QCD-weak subset, soft singularities appear only in processes (3) and (4) and cancel when the two are summed. When we consider the combination of the cross sections of the processes (3)-(8) we are thus left with initial state collinear singularities only. The processes (3)-(7) contribute to initial-state collinear singularities within the QCD-QED subset, while in the QCD-weak case only processes (3)-(4) have initial state collinear singularities of QCD origin. These singularities can be removed by mass factorization. The physical parton densities f i (x, µ F ) are defined, at the factorization scale µ F , by introducing the mass factorization kernel Γ ij , which subtracts the initial state collinear singularitieŝ The kernel can be expanded as a series in α and α s where Γ (1,0) ij is the QCD leading order (LO) splitting kernel, Γ (0,1) ij is its QED analogue and Γ (1,1) ij is the mixed QCD-QED contribution to the splitting kernels, recently presented in Ref. [70]. After the replacement of Eq. (10) in Eq. (1), we obtain the total cross section expressed in terms of subtracted, finite, partonic cross sections σ ij (z, µ F ): The σ ij admit a perturbative expansion in powers of α and α s , in analogy to Eq. (2). In this letter, we present the results for σ (1,1) uū . In processes (3) and (4) the weak virtual correction to the splitting vertex q → qg might induce an additional contribution to the subtraction kernel Γ (1,1) ij . However, we have checked that such a term vanishes, in the massless quark case, as a consequence of the conservation of the vector and axial-vector currents.

RESULTS
In order to discuss the size of the different sets of radiative corrections, we define: is the Born cross section of the process uū → Z, with N c the number of colours and C v/a,u the vector/axial-vector couplings of the Z boson to the up quark. ∆ (1,1) uū,K with K = γ, Z, W are the corrections due to the exchange of a photon, a Z boson, and of one or two W boson/s including the lowest order charge renormalization counterterms, respectively. For the sake of comparison, we introduce the NLO-QCD correction to the same partonic process, defined as α s σ In Figure 1 we present, as a function of the partonic  We observe that, in the high-energy limit (z → 0), the cross sections are damped by the incoming flux factor, proportional to z. The divergent behaviour for z → 1, due to the exchange of at least one massless boson, is also evident for all the contributions. The values of the EW charges, in the two subsets with one Z (red) or with one/two W s exchange (magenta), are responsible for the different size and for the opposite sign of the two contributions, visible in the z → 1 limit. We observe that in the case of the dd → Z + X process, the contributions with one/two W s exchange have similar size but opposite sign. The total contribution to the hadron-level cross section from this subset of diagrams of the two partonic processes is expected to undergo an important cancellation, modulated by the convolution with the proton PDFs. The QCD-QED corrections, shown in blue in Figure 1, are not monotonic, contrary to the NLO-QCD ones and have a maximum for z ∼ 0.85. They are smaller than the QCDweak contribution for z ∈ [0.8, 0.9], but become larger in absolute size when z → 1, because of the higher power of the threshold logarithms. The possibility of having a second Z boson in a resonant configuration yields the kink of the ∆ (1,1) uū,Z curve (red) at z = 1/4, as it can be observed in the inset of Figure 1.
In conclusion, we have presented the first results for the total inclusive partonic cross section for the process qq → Z+X, including the exact O(αα s ) corrections, with both photon and W/Z boson exchanges. The results are analytic and are expressed in terms of GPLs, but also contain three elliptic MIs, which have been computed with a series expansion around z = 1. The complete solution of the infrared structure of the process and the exact evaluation of all the relevant virtual corrections represent an important step towards the evaluation of the hadron-level cross section for Z production at this perturbative order.