On-shell Z boson production through ${\cal O}(\alpha \alpha_s)$

The analytical expressions of the mixed QCD-EW corrections to on-shell Z boson inclusive production cross section at hadron colliders are presented, together with computational details. The results are given in terms of polylogarithmic functions and elliptic integrals. The impact on the prediction of the Z boson production total cross section is discussed, comparing different proton parton density sets.


Introduction
The production in hadronic collisions of a pair of leptons, each with large transverse momentum, is known as Drell-Yan (DY) process and it plays a fundamental role for our understanding of Quantum Chromodynamics (QCD) as the theory of the strong interactions.The lepton pair acts as a probe of the initial-state proton structure.It allows for the measurement of the proton collinear parton density functions (PDFs) and for the study of the QCD dynamics from the analysis of the lepton-pair transverse momentum distribution.The kinematical distributions of the final-state leptons allow for precision tests of the electroweak (EW) Standard Model (SM), with the determination of the weak mixing angle [1,2], and of the masses m W,Z [3,4] of the W and Z bosons.
The production of an on-shell Z boson represents an approximation of the full neutralcurrent DY process, in that special kinematical configuration dominated by the Z-boson resonance.The total cross section of this process is an important theoretical benchmark constraining the proton PDFs, so that its prediction is a cornerstone in the precision physics program at hadron colliders.This process is described in lowest order (LO) by quark-antiquark annihilation into a Z boson via EW interaction.The evaluation of the next-to-leading order (NLO) [5], next-to-next-to-leading order (NNLO) [6][7][8][9], and nextto-next-to-next-to-leading order (N 3 LO) [10][11][12][13][14] QCD corrections to the production of an on-shell gauge boson, supplemented by the resummation of the logarithmically enhanced terms due to soft gluon emission [15][16][17][18][19][20][21][22], allows for the accurate estimate of the total cross section, the reduction of the impact of QCD uncertainty and the precise assessment of its actual size (cfr.also Ref. [23]).
The best available QCD prediction for the inclusive production of a virtual gauge boson includes up to N 3 LO QCD corrections [12][13][14], in the γ * or W cases respectively.It shows a dependence on the QCD renormalisation (µ R ) and factorisation (µ F ) scale choices at the sub-percent level for virtualities Q > 70 GeV and at the percent level for smaller Q values, with a stronger sensitivity to the choice of the factorisation scale.
In this high-precision QCD framework, the inclusion of EW effects becomes mandatory.The NLO-EW corrections to the DY process have been computed in [24][25][26][27][28] and are comparable in size to the NNLO-QCD effects.The theoretical uncertainty associated with missing higher-order EW corrections is formally at the NNLO-EW level and it is significantly reduced compared to the leading order (LO) case.Using different input parameters as a mean to estimate the size of missing higher-order EW effects, one finds that the LO variation is at the O(3.5%) level, whereas the NLO-EW one is reduced down to the O(0.5%) level.Despite of these significant progresses, it is possible to observe the presence of some residual sources of uncertainty in the results listed above.The higher-order QCD predictions are only LO accurate from the point of view of the EW interaction, and thus they suffer from the uncertainty associated with the different choices of input parameters.If we consider the NNLO-QCD prediction, supplemented with the NLO-EW one, we find a scheme uncertainty at the O(0.88%) level.This value is significant for any precision test, it is comparable to the residual QCD uncertainty, and in any case is slightly larger than the corresponding estimate based only on the LO+ NLO-EW results.A specular discussion applies to the NLO-EW corrections, which are only LO from the point of view of the strong interaction.They suffer from large uncertainties under variations of the factorisation scale.The canonical µ F variation by a factor 2 about its central value yields a change of the LO cross section by ±18% and, in turn, a change of the NLO-EW correction at the O(0.5%) level.In order to increase the control on the theoretical error, it is therefore mandatory to include in the analysis the full mixed QCD-EW corrections, since they stabilize both the dependence on the QCD scales of the higher-order EW corrections and the dependence on the EW input parameters of the higher-order QCD corrections.
The mixed corrections to the DY process in the resonance region have been studied in the so-called pole approximation in Refs.[29,30], where QCD and EW effects are factorised between production and decay of the vector boson.As far as the production of an on-shell Z boson is concerned, the exact QCD-QED corrections have been considered in Refs.[31][32][33], while in [34][35][36] we have computed the mixed QCD-EW effects fully analytically for the inclusive production.In Ref. [37], the fully differential QCD×EW effects have been presented, using a combination of analytical and numerical techniques.As for the full neutral-current DY, the QCD×QED corrections to the production of a pair of neutrinos have been discussed in Ref. [38], and the complete NNLO QCD-EW corrections for a charged lepton pair in the final state have been presented in [39].Considering pole-approximation only for the virtual contributions, the QCD-EW corrections to the charged current DY has been obtained in [40].
In this paper, we explicitly present the analytical expressions and all the computational details to obtain the O(αα s ) mixed QCD-EW corrections to on-shell Z boson inclusive production cross section at hadron colliders.The results have been fully computed in analytical form and expressed in terms of polylogarithmic functions and elliptic integrals, requiring the evaluation of new two-loop master integrals (MIs) that were not available in the literature.The inclusion of the O(αα s ) corrections increases the accuracy of the prediction and it reduces the impact of the residual theoretical uncertainties.The evaluation of the hadron-level cross section requires the convolution of the partonic results with proton PDFs.For consistency, the latter must satisfy DGLAP equations with also a QCD+QED evolution kernel.We comment on the accuracy of our predictions, compared to a standard analysis that includes only QCD corrections.

The Z boson production cross section
In this section, we discuss the theoretical framework, providing the details of perturbative contributions from various partonic channels which constitute the mixed QCD-EW corrections to the inclusive production cross section of the Z boson at hadron colliders.We also discuss here the necessary steps to perform ultraviolet (UV) renormalisation and mass factorisation.
The hadron-level cross section for the inclusive production of an on-shell Z boson can be written, according to the factorisation theorem, as 2) The bare PDFs, f h i (x), describe the partonic content of the hadron h, where i can be a quark (q), anti-quark (q), gluon (g) or photon (γ).The partonic cross section σij ≡ σ(ij → Z+X), describes the inclusive production of a Z boson in the scattering of partons i and j.
In the evaluation of the radiative corrections, the presence of UV and IR divergences is handled in dimensional regularisation, with d = 4 − 2ε being the number of space-time dimensions.The individual matrix elements contain in general poles in ε.The cancellation of the UV singularities takes place via UV renormalisation.The combination of virtual and soft real emission corrections to the same underlying process leads to a cancellation of the IR soft poles.According to the Kinoshita-Lee-Nauenberg (KLN) theorem [41,42], the combination of the different partonic cross sections contains only singularities due to the emission of collinear partons from the initial state.The latter are universal, they can be factorised and reabsorbed in the definition of the physical proton PDFs, , by means of the mass factorisation kernel Γ , defined at the factorisation scale µ F .As a consequence, Eq. (2.2) can be written as follows with the convolution represented by the symbol ⊗, The constant σ 0 is defined through the Born cross section σ (0,0) , In the latter, z ≡ m 2 Z /ŝ, √ ŝ is the partonic center of mass energy, N C is the number of colors, and c ) is the combination of charges of the coupling of the Z boson to a quark q.C v,q and C a,q are given by W and Q q the third component of the weak isospin and the electric charge in units of the positron charge, respectively.θ W is the weak mixing angle.
∆ ij is the UV-and IR-finite partonic cross section for the partonic channel ij, expressed in units σ 0 .In perturbation theory, ∆ ij is expanded in series of α s and α as The main results of this paper are the expressions of the corrections ∆ (1,1) ij with i, j = q, q, g, γ.In the study of the NNLO QCD-EW corrections, we encounter different gauge invariant subsets, characterised by the exchange of additional photons or massive weak bosons W s or Zs.We reorganize accordingly the total O(αα s ) correction, introducing   W .

The partonic subprocesses
The inclusive production cross section receives, order by order in perturbation theory, virtual as well as real emission corrections.The additional final-state partons present in the latter are completely integrated over their respective phase space.The first term of the expansion σ (0,0) in Eq. (2.1), receives contributions from the single partonic process (see Fig. 1 (a)) q + q → Z . (2.8) σ (1,0) in Eq. (2.1) receives contributions from the partonic processes q + q → Z , q + q → Z + g , q + g → Z + q . (2.9) The first process receives NLO-QCD virtual corrections (Fig. 1 (b)), while the other two are evaluated at tree level, with the phase-space of the emitted parton (g/q) integrated out (Fig. 1 (c),(d)).In case of σ (0,1) , the processes are (see Fig. 2) where the first process receives NLO-EW virtual corrections (Fig. 2 (a)).At O(αα s ) we have double-virtual, real-virtual and double-real contributions.Double-virtual corrections are two-loop contributions, with one gluon and one/two EW gauge bosons in the loop, to the partonic process (see Fig. 3) q + q → Z . (2.11) In the real-virtual contributions we find one virtual loop and one real-emitted particle.The processes which contribute to this group are (see Fig. 4) q + q → Z + g , q + g → Z + q , q + q → Z + γ , q + γ → Z + q . (2.12) For the first two processes, the loop integral is of O(α) (Fig. 4 (a)-(h)), while for the others it is of O(α s ) (Fig. 4 (i)-(p)).At O(αα s ) , the double-real contributions are with two real-emitted partons.Their amplitudes are evaluated at tree level.In the cases with two (anti)quarks in the initial and two (anti)quarks in the final state, in addition to the Z, the scattering is mediated by either a gluon or an EW boson, so that the respective contributions have a different proportionality in terms of α and α s .The interference between the two groups of contributions is of O(αα s ) with respect to the Born process and is relevant for the calculation of σ (1,1) .The complete list of processes is (see Fig. 5) q + q → Z + g + γ , q + q → Z + q + q , q + q → Z + q + q , q + g → Z + q + γ , q + γ → Z + g + q , g + γ → Z + q + q , q + q → Z + q + q , q + q → Z + q + q , q + q → Z + q + q . (2.13) We denote with q a different quark flavour.
Examples of Feynman Diagrams contributing to the real-virtual corrections to the production of a Z boson in hadronic collisions, Eqs.(2.12).
We analytically compute the contribution at O(αα s ) to the total cross section of each of the above listed processes, which can be presented as a Laurent expansion in ε.After the combination of all the degenerate states and mass factorisation, the singularities cancel and the remaining non-vanishing finite contributions are the factors ∆ (1,1) ij .We present in Section 4 their analytical expressions, and provide them also in a Mathematica ancillary file.

Classification of the radiative corrections
The O(α) and O(αα s ) corrections can be organised in a gauge invariant way into a photonic and two weak subsets.The photonic subset (∆ γ ) contains all contributions involving a photon i.e. a virtual photon in the loop or a real-emitted photon or a photon-initiated Examples of Feynman Diagrams contributing to the double real corrections to the production of a Z boson in hadronic collisions, Eqs.(2.13).channel; at O(αα s ) these are dubbed QCD-QED corrections.The first weak subset (∆ Z ) contains contributions from an additional Z boson, either in the loop or as the intermediate propagator in the tree diagrams of double-real contributions.The other weak subset (∆ W ) contains contributions from one or two W bosons and includes Feynman diagrams with a non-abelian trilinear gauge boson vertex.This subdivision makes the computation, especially the mass factorisation, more organised and is useful for several intermediate checks.

Ultraviolet renormalisation
The prediction of the hadron-level cross section requires to express the bare couplings and masses in terms of physical parameters via renormalisation.The choice of the background field gauge (BFG) [43] allows to restore the validity of U (1) em -like Ward identities between the vertex corrections and the external quark wave function corrections in the full EW model.We explicitly verify these Ward identities at O(α) and O(αα s ) .We observe at O(α) that the sum of the photonic correction in the vertex and external quark wave function corrections is UV finite.The same holds for the corrections with the exchange of one virtual Z boson and, separately, for those with one or two W bosons.An identical pattern takes place at O(αα s ) .We are thus left with the discussion of the charge and external gauge boson wave function renormalisation.
If we choose to express (g, g , v), the SU (2) L and U (1) Y couplings and the Higgs field vacuum expectation value, in terms of (G µ , m W , m Z ) (dubbed G µ -scheme), where G µ is the Fermi constant, then the weak charge renormalisation is achieved by the replacement We denote with a 0 subscript all the bare quantities.We abbreviate with c W = m W /m Z the cosinus of the weak mixing angle (s 2 W = 1 − c 2 W ), and we define δg where Z ZZ = 1 + δZ ZZ is the ZZ wave function renormalisation constant.∆r is a finite correction [44] expressing the relation between the Fermi constant and the muon decay amplitude, δs are the electric charge and gauge boson mass counterterms.The δg Z factor is, in the BFG, an UV finite correction and this fact yields a considerable simplification in the study of the impact of different input schemes.The ∆r parameter and the counterterms can be evaluated in perturbation theory and we keep terms of O(α) and O(αα s ) [45,46].For consistency, they have to be expressed in terms of (G µ , m W , m Z ).In addition to the redefinition of the overall weak coupling, a second renormalisation correction modifies the vector coupling v q of the Z boson to the quarks: , with δZ AZ the renormalisation constant of the γ − Z mixing.In the BFG also this shift of v q is UV finite.
If we instead choose to relate (g, g , v) to the (α, m W , m Z ) set of inputs (dubbed α(0)scheme), the replacement of the overall coupling is δg Z , while the redefinition of the vector coupling remains the same as in the other scheme 1 .
The α(0)-scheme choice is historically [44,48] one of the simplest EW renormalisation input schemes, but the low scale at which the fine structure constant is measured yields in turn the appearance of large logarithmic corrections in the perturbative expansion.The results depend on the value of the light-quark masses or, alternatively, on an experimental input ∆α had (m Z ) = 4π γγ (0) [49][50][51][52] needed to evaluate the hadronic contribution to the running of the electromagnetic coupling at low scales, where Π (5) γγ (q 2 ) indicates the contribution of the first five light quark flavors to the photon vacuum polarisation at a scale q 2 (cfr.Ref. [46]).

Infrared singularities and mass factorisation
Each UV-renormalised process, Eqs.(2.11-2.13), is in general IR divergent, because of the exchange of a soft and/or collinear gluon or photon.Thanks to the KLN theorem, after combining all the partonic subprocesses and summing over all degenerate states, only initial-state collinear singularities are left.The latter are absorbed in the definition of the physical proton PDFs, by means of the mass factorisation kernel Γ .The subtraction kernels at O(αα s ) are based on the splitting functions computed in Ref. [53].The photonic initial-state collinear singularities are reabsorbed in the definition of the physical proton PDFs and are resummed to all orders via the DGLAP evolution of the parton densities with a QED kernel.The consistent evaluation at O(αα s ) of the hadron-level cross section requires a proton PDF set featuring DGLAP QED evolution, and, as already presented in Eq. ( 2.3), the inclusion of photon-induced subprocesses.

Computational details
In this section we present the details of the calculation.We first introduce our general strategy, which is based on the conversion of all the phase-space integrals into loop integrals, using the so-called reverse unitarity approach.We then discuss in detail the challenges posed by some new MIs with internal massive lines.We eventually provide a detailed account of how we deal with the appearing elliptic integrals.

General strategy
The evaluation of the partonic cross-sections beyond LO requires the computation of Feynman loop integrals from virtual diagrams as well as two-and three-particle phase-space integrals arising from real emissions.In the very beginning of inclusive NNLO calculations, the phase-space integrals were performed using parametric and angular integration.However, to benefit from state-of-the-art techniques, developed for virtual integrals, we use the method of reverse unitarity [54,55] to convert the phase-space integrals into loop integrals.The reverse unitarity technique relies on the fact that the following replacement, Cutkosky rule, allows to convert the Dirac delta function in the phase-space measure of each final state particle into the difference of two propagators with opposite prescriptions for their imaginary part, where 0 + is an infinitesimal positive real number.The resulting integrals can then be studied with the help of standard techniques for virtual integrals.
The amplitudes that one needs to calculate can be classified as follows: two-loop 2 → 1 virtual amplitudes and two-loop 2 → 2 forward-scattering amplitudes with two-or three-particle cuts, stemming from the real-virtual and double-real corrections, respectively.There are up to two internal masses, which are in general complex valued.Once the Feynman diagrams contributing to these processes are generated using tools like Qgraf [56] or FeynArts [57], we compute the interference with the corresponding tree-level.We use in-house Form [58] or Mathematica codes for this algebraic part of the computation.
The interference term is expressed in terms of a large number of scalar integrals, that are not all independent, evaluated in d space-time dimensions.Dimensionally regularised scalar integrals satisfy integration-by-parts (IBP) identities [59][60][61].These identities link different scalar integrals with each other and make possible the reduction of a large number of terms to a small set of independent quantities, called the MIs.For the IBP reduction process, we have used the public programs Kira [62], LiteRed [63,64] and Reduze2 [65,66].
The EW corrections entail the presence of masses in the loop propagators, specifically m Z and m W . Accordingly, the integrals that have to be evaluated depend, in general, on three scales (or two dimensionless ratios).The fact that the two masses are numerically very close to each other can be used to reduce, effectively, the number of scales in the loop integrals, minimising the complexity of the calculation of the MIs.In fact, we can conveniently express the W boson squared mass as This makes it possible to expand the loop integrands, which contain m 2 W , as a Taylor series in the parameter ξ.Thus, loop integrals will only depend on ŝ and m 2 Z , hence z, and the dependence on m 2 W is formally confined into the coefficients of those integrals.The Taylor expansion in ξ will generate loop propagators with higher powers, the so-called dotted propagators, without affecting the topology of the diagrams.Therefore, the set of MIs will not change.The number of terms of such power expansion that have to be retained depends on the phenomenological accuracy we need for the evaluation of the corrections.After IBP reduction, the generic structure of the squared matrix element integrated over the inclusive phase space will then be where the c k,j are rational functions of the kinematical invariant z and the number d of space-time dimensions, while the I k (z, d) are the N MIs of the process.In this work, we have considered only up to the n = 2 order in the ξ expansion 2 and we will discuss the impact of this truncation in Section 5.
The following step is to evaluate the MIs as a function of z and d.We use the method of differential equations [67][68][69][70][71][72][73][74] to solve the MIs.We differentiate one MI with respect to z and we use the IBP identities on the differentiated output to express it as a combination of the MI itself and of other MIs.Applying this procedure on all the MIs of an integral family, we obtain a system of first-order linear differential equations, which can be solved given a set of boundary conditions (BCs).In the current calculation the differential equations for most of the MIs, can be solved expressing the solution in terms of harmonic polylogarithms (HPLs), generalised harmonic polylogarithms 3 (GPLs) [75][76][77][78] and cyclotomic HPLs [79].
Three MIs in the evaluation of the double-real corrections, need elliptic extensions of such functions, known as elliptic polylogarithms.

Evaluation of the full set of Master Integrals
The MIs which appear in the reduction of the squared matrix element of the different partonic subprocesses can be classified in three groups associated to the two-loop virtual corrections, the two-particle phase-space integrals of the real-virtual corrections and the three-particle phase-space integrals of the double-real processes.All the two-loop virtual MIs were already available in the literature [80][81][82][83][84][85][86][87].We recompute them using the method of differential equations, starting with a generic off-shell value of the Z boson virtuality and then taking the on-shell limit of the results.The latter are expressed in terms of multiple zeta values and cyclotomic constants [79], which can be reduced to a set of independent constants as introduced in [88].The off-shell virtual integrals and most of their on-shell limits have been independently checked using Fiesta [89].
In the real-virtual and double-real corrections, the MIs with only massless internal lines were available from Ref. [55].The new MIs with internal massive lines in real-virtual and double-real contributions are solved by employing again the differential equations technique.In these systems of differential equations, square root letters appear, making a direct solution cumbersome.To solve this problem, it is customary to perform a change of variables and to rationalize the square root factors.The presence, in some of the MIs under study, of multiple square roots makes it impossible to rationalize all the letters with one single variable transformation.In these cases, we exploit the linearity property of the integral operator and divide the MI into two (or more) subsystems.At this stage, we can choose for each subsystem a different transformation rule which rationalizes its specific letters.We trade the explicit dependence on z of the full solution with the one on several new kinematical variables, but we obtain a simpler representation of the solution in terms of a compact alphabet as presented later in Eq. (4.2).
Some of the MIs of the present study satisfy coupled differential equations and the problems of rationalising the letters and decoupling the equations are intertwined.We describe in the following example how the two issues can be simultaneously handled.

Example : A subsystem of two real-virtual master integrals
We consider two integrals J 1 ≡ J 1 (z) and J 2 ≡ J 2 (z) (see Fig. 6), which satisfy the following system of differential equations: The functions r 1,2 are the inhomogeneous part of the equations and contain rational functions and GPLs.With the transformation rule . Two MIs, J 1 and J 2 , for the real-virtual corrections.Thin plain lines represent massless particles, while thick plain lines represent massive particles.
and the relations the first differential equation of Eq. (3.4) can be rationalised and, with appropriate boundary conditions, can be solved to obtain up to O(ε −1 ) as This allows us to readily obtain the solution for J . However, the differential equation for J (−1) 2 has contributions to its non-homogeneous part from r as well as J (−1) 1 . The transformation rule in Eq. (3.5) does not rationalize the former contributions, but it rather converts the linear kernels to much more involved polynomials.At this point the difficulty of solving the system for J 1,2 is evident, as cumbersome letters would appear in the GPLs, making the next integrations impossible.We reorganize the problem with the working hypothesis that, at two-loop level, all the coefficients of the poles in ε are expressible in terms of simple GPLs.In other words we expect to be able to find a combination of such that their single pole can still be determined with elementary integrations.
With this goal, we define a new integral J 0 and we observe that J (−1) 0 satisfies a linear differential equation independent of both J (−1) 1 and J (−1) 2 . We solve it and obtain Consistently replacing each J (n) 2 using Eq.(3.7), both in the system of differential equations and in the matrix elements, we also remove the dependence of J (0) 1 .However, the contribution from J (−1) 1 , can not be eliminated, as expected.Hence, the issue of rationalisation remains in the full non-homogeneous part of J (0) 0 .The non-homogeneous part of J (0) 0 can be split in the sum of two terms, with or without a relation with the square root letter: we use the change of variables to w to linearise only for the former, while we keep the variable z for the latter.When we solve the non-homogeneous equation, we consider two separate integrals, with the integrand functions depending only on w or only on z, with straightforward results in terms of GPLs, in the following form The absence of the square root in the letters allows a smooth numerical evaluation.In the previous example, we find that a certain combination of MIs can be helpful to avoid the appearance of "complicated" GPLs in the intermediate steps of the solution of the system, under the assumption that they are not expected in the coefficients of the poles in ε.This remark simplifies in turn the solution of the system for the finite part of the MIs.For example, the contributions from J (0) 1 and J (0) 2 do individually contain the problematic GPLs, which however would eventually cancel in the final result.Our direct evaluation of J (0) 0 instead avoids these GPLs from the very beginning.
Apart from the square root letters, we also find three MIs in the double-real corrections for which the coupled sub-system of differential equations can not be factorised to first order, indicating that the MIs contain elliptic polylogarithms.In the next sub-section we present the details of the computation of these MIs.

Solving the master integrals with elliptic kernels
The double real corrections receive contributions from scattering processes with quark and anti-quark both in the initial and final state, whose matrix elements include contributions with the exchange of either an EW boson or a gluon.The interference of these two terms is of O(αα s ) and is relevant to our calculation.In these interferences, a non-trivial topology gives rise to three MIs {I 1 , I 2 , I 3 } (see Fig. 7) which are of elliptic kind.The 3 × 3 system of differential equations is not first-order factorizable.The homogeneous part of the system is given in each order of ε as The homogeneous part of this system of equations is the same as the one studied for the corresponding virtual diagrams in Refs.[85,90].In those papers, the results are obtained in terms of elliptic integrals of the first kind and eMPLs, respectively.The current system can be solved order-by-order in ε, and we expect to find only standard HPLs in the poles, while eMPLs will appear in the finite parts and higher orders in ε.
Formally, since the IBP reduction introduces a 1/ε factor in the coefficient of these integrals in the matrix elements, then we would need to evaluate I 1,2,3 up to O(ε), in order to compute all the finite corrections.For the same argument, we expect that individual contributions to the coefficient of the single pole of the matrix elements are expressed via eMPL.We formulate a physical Ansatz as in the previous Section, imposing the simpler polylogarithmic structure of the single pole of the matrix element, i.e. the absence of eMPLs in its final expression.We thus find the following combination of the elliptic MIs, (3.10) The differential equation of I 3 .We solve it to obtain the following solution: This new combination of MIs explicitly exhibits the absence of eMPL in their single pole and allows in turn a straightforward check of the pole cancellations in the cross section.
Clearly, since the combination I (1) 0 does not remain independent of 3 we find eMPLs contributing to the finite part of the matrix elements.For the numerical evaluation, a series representation of the integrals is equivalent to the formal solution via eMPLs.Hence, we solve I 0 instead contains logarithmically enhanced terms.The latter yield a singular behaviour at a point which does not correspond to any physical threshold.These logarithms are thus expected to cancel against analogous terms stemming from other MIs.In order to achieve an exact analytical cancellation, we further elaborate the solution of I ell , by splitting its expression in two parts: one with the closed form of the logarithmic dependence and one regular remainder given via a Taylor expansion.We obtain ell .
are obtained in the two series expansions as follows The series expanded δ (0,1) 2 and I (0) 3 , are all regular and admit a Taylor expansion with fast convergence.This procedure provides a series representation of the elliptic functions present in the calculation.In Fig. 8, we present the numerical evaluation of the elliptic MIs for the whole range of z.

Analytical results
In this section we present our results.First, we introduce all the variables which appear in our calculation and several other necessary definitions.We then present the analytic finite partonic coefficients

Preliminaries: variables, functions and abbreviations
The partonic cross-section for the production of a Z boson, integrated over the phase-space of the additional emitted partons, depends on the variables s, m Z and m W and on the top quark and Higgs boson masses m t and m H respectively only through UV renormalisation.As described in Sec.3.1, we expand the virtual EW amplitudes containing W bosons according to Eq. (3.2).For the purpose of clarity and compactness, we present here only results for the n = 0 term in the ξ expansion, i.e. m 2 W = m 2 Z .The Feynman integrals have a rich structure due to the presence of different internal thresholds.When the solution is given in terms of polylogarithmic functions4 , the internal structure of the integrals is displayed by the values of the weights and of the independent variable of the polylogarithms.The singular or branching points of the solutions are expressed in terms of a set of monomials, the letters, forming an alphabet.The presence in the letters of square-roots, also at multiple levels, can be avoided by appropriate transformations which rationalise the initial expressions.In the current calculation we adopt the following changes of variables: After these manipulations, the resulting alphabet of the problem is given by: 2) The set {−1, 0, 1} is the well-known alphabet for HPLs.{{3, 0}, {3, 1}, {4, 1}, {6, 0}, {6, 1}} denotes the third-, fourth-and sixth-root of unity, which defines the cyclotomic HPLs.The set {{4, 1}, i 1 , i 2 } appears in a group of diagrams expressed in terms of GPLs, where i 1 and i 2 are given by We enlist all our definitions below In order to present the results in compact form, we introduce the following abbreviations for all the polynomials which appear in the denominator of the expressions: . (4.8)

The virtual corrections
In this section we present the results for the sole virtual corrections.The two-loop Feynman diagrams contributing to the renormalised q qZ vertex are shown in Fig. 3.We consider the interference with the tree-level and we express it in perturbative expansion of α and α s as follows The photonic contributions are denoted by F (m,n) γ and they have been presented in [31,97].Below, we present the renormalised contributions with one Z boson or one/two W boson/s in the loop, denoted by , respectively.We present here the results for the zeroth order in ξ.We have also computed the completely general case with different masses and compared it against the results available in the literature in [87], finding agreement.We performed the calculation in the background field gauge.Therefore, we included the sole wave function renormalisation, which is sufficient to obtain UV finiteness.In our expressions, the following constant c 1 [98] appears The one-and two-loop results for , in the case of u-quark initial state, are in order:  where we defined the following combinations of vector and axial-vector couplings: with C v,q and C a,q defined in Eq. (2.5) and q ∈ (u, d).

The partonic coefficients for QED
The photonic part (∆ (1,1) γ ) of the total hadronic cross-section has been defined in Eq. (2.7) and receives contributions from several partonic channels which are convoluted with the physical proton PDFs as follows 1)  gγ . (4.17) The sums run over all flavours of quarks (q) and antiquarks (q).The following combination of vector and axial-vector couplings, along with electric charge, appear in all these partonic channels with a photonic correction: In Eq. (4.18) q can be an up-or a down-type quark.The partonic cross-section for all the channels are given below.We have set the factorisation scale µ F = m Z in presenting these results.

The partonic coefficients for W
In this Section, we provide the total partonic cross-sections stemming from all the channels with Feynman diagrams where one or two internal W bosons are exchanged.Like before, the hadronic cross-section receives contributions from several partonic cross-sections convoluted with the PDFs as: where ∆r (1,1) , δg and δs 2
can thus be taken as an estimate of the impact of the mixed QCD-EW corrections on this observable 6 .We observe in Table 1 that different PDF choices lead to a spread in the central value prediction; comparing the width of their envelope with the mean value of the three predictions, in the QCD-only case, we observe a dispersion of the results ∆ env ranging between the O(1%) and the O(2%) level at the LHC energies, while a smaller value is observed at the Tevatron.This spread is compatible with the estimate of the PDF uncertainty evaluated according to the Monte Carlo or Hessian recipes used by each PDF collaboration.We remark that the PDF uncertainties range between the O(1%) and the O(±2.5%)level, depending on the energy, collider and PDF fitting group; larger values are instead obtained for a collider with √ S = 100 TeV.
For each PDF set, the shift induced by the inclusion of the  QCD-EW corrections is almost independent of collider and energy and ranges between -0.4% and -1.4%.The size of this shift of the central value is significant if compared with the residual subpercent QCD scale uncertainty reported in Refs.[12,13] at N3LO-QCD, but also with the corresponding estimate made at NNLO-QCD.In Table 2 we present the width of the QCD scale uncertainty band, evaluated considering independent variations of the renormalisation scale µ R = ξ R m Z and factorisation scale µ F = ξ F m Z .We consider 9 combinations obtained varying ξ R,F ∈ [ 1 2 , 1, 2] and we discard the two cases ξ R = 1 2 , ξ F = 2 and ξ R = 2, ξ F = 1 2 .The percentage correction is computed with respect to the central scales choice ξ R = ξ F = 1.
Since the parameterisation of the proton must fulfil the sum rules expressing charge and momentum conservation, then the presence of the photon density in the proton implies a reduction of the total momentum fraction carried by quarks and gluons.In turn, we observe a reduction of the cross section of the partonic channels induced by quarks and gluons, which in general should be compensated by the positive cross section of the additional photon-induced channels.We remark that a similar O(1%) reduction of the quark-and gluon-induced DY cross section for the production of a lepton pair is balanced by the contribution of the γγ → + − process, of comparable size [25,27,109].In the specific case of on-shell Z production, the γγ initial state does not contribute at O(αα s ) , but only at O(α 2 ) ; its absence explains, at technical level, the size and sign of the observed effects.
The relevance of the total on-shell Z production cross section as a standard candle for benchmarking purposes is well established.In the present study, it allows to observe how different PDF collaborations have implemented the QCDxQED evolution and the photondensity, with a different impact with respect to the corresponding pure QCD analysis.

Conclusions
We have presented the analytical expression of the total cross section for the production of an on-shell Z boson at hadron colliders.The result has been expressed in terms of generalised polylogarithmic functions.The three elliptic functions which appear in the double-real contributions have been represented via a series expansion.
These corrections stabilise the prediction of this standard candle, by reducing the size of the uncertainty due to missing higher orders [36].The introduction of EW and mixed QCD-EW radiative corrections requires the consistent usage of proton PDFs determined in the same theoretical framework.The comparison with the best prediction obtained in a pure QCD-based model shows that QCD-EW effects up to O(−1%) have to be taken into account and are relevant for the precision determination of this cross section.

Figure 1 .
Figure 1.Feynman Diagrams contributing to the Born (a) and to the O(α S ) corrections (b)-(d) to the production of a Z boson in hadronic collisions.Crossed diagrams are not shown.

Figure 2 .
Figure 2. Feynman Diagrams contributing to the O(α) corrections to the production of a Z boson in hadronic collisions.Crossed diagrams are not shown.
Figure 3. Two-loop Feynman Diagrams contributing to the O(αα S ) corrections to the production of a Z boson in hadronic collisions, Eq. (2.11).Symmetric diagrams are not shown.

Figure 7 .
Figure 7.The three elliptic MIs, I 1 , I 2 and I 3 , contributing to the real-virtual corrections.Thin plain lines represent massless particles, while thick plain lines represent massive particles.A dot on a line represents the squared propagator.

Table 2 .
Dependence of the cross sections for on-shell Z production, in picobarns, on the renormalisation and factorisation scale choices.The upper and lower percentage variations, compared to the central scales choice, is computed among 7 scales combinations (cfr.text).The results have been computed with the central replica of the NNPDF31 nnlo as 0118 and NNPDF31 nnlo as 0118 luxqed PDF sets.