NNLO EW⊗QCD corrections to on-shell Z production

In this proceeding, we outline the computational details and present the results for NNLO mixed EW⊗QCD corrections to on-shell production of Z boson at the LHC. We use the method of differential equation to obtain all the master integrals analytically, appearing in virtual and real-emission contributions. These results increase the accuracy of the predictions and contribute to reduce the theoretical uncertainties. Copyright N. Rana. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 27-11-2021 Accepted 31-03-2022 Published 20-06-2022 Check for updates doi:10.21468/ciPostPhysProc.7.004


Introduction
The Drell-Yan (DY) process of the production of a lepton-pair in hadronic collisions is very important for our understanding of Quantum Chromodynamics (QCD) as the theory of the strong interactions, because of its significance in the precise measurement of the parton density functions (PDFs) and for the precise determination of the weak mixing angle and of the masses m W , m Z and decay widths Γ W , Γ Z of the W and Z bosons, respectively.
The on-shell Z boson production is a particular kinematical configuration of the full neutralcurrent DY process. The next-to-leading order (NLO) [1], next-to-next-to-leading order (NNLO) [2][3][4][5], and next-to-next-to-next-to-leading order (N 3 LO) [6][7][8][9] QCD corrections to the production of an on-shell gauge boson along with the resummation of the logarithmically enhanced terms due to soft gluon emission [10][11][12][13][14][15][16], allow to obtain a precise estimate of the total cross section. These computations also reduce the impact of QCD uncertainty and precisely assess its actual size.
The electroweak (EW) effects become very much necessary to include in such high-precision QCD framework. In [17][18][19][20][21], the NLO-EW corrections to the DY process have been obtained and the size of these corrections are comparable to the NNLO-QCD effects. From the perspective of the EW interaction, the higher-order QCD corrections are only leading order (LO) accurate. Hence, the different choices of the input parameters associate large uncertainties to the 004.1 higher-order QCD corrections. For example, we may choose to express the SU(2) L and U(1) Y couplings and the Higgs field vacuum expectation value, in terms of (G µ , m W , m Z ) (the G µscheme) where G µ is the Fermi constant, or in terms of (α, m W , m Z ) (the α(0)-scheme) [22,23] where α is the electromagnetic coupling constant. The LO differs by O(3.5%) for these two different input schemes and higher-order QCD corrections do not change the variation, whereas inclusion of the NLO-EW corrections reduces it to O(0.5%). If we consider the NNLO-QCD corrections, supplemented with the NLO-EW one, the input scheme uncertainty is at the O(0.88%) level, which is significantly large for any precision test, and hence, must be reduced. On the other hand, the NLO-EW corrections are only LO from the QCD viewpoint. They suffer from large uncertainties under variations of the factorization (µ F ) and renormalization (µ R ) scales. Hence, in order to increase the control on the theoretical uncertainties, it is necessary to include the mixed QCD-EW corrections, since they stabilize both the dependence on the unphysical scales and the dependence on the EW input parameters.
In Refs. [22][23][24][25], the charge renormalization and the corresponding counter-terms, necessary for the mixed QCD-EW corrections, were discussed. In Refs. [26,27], the mixed corrections to the DY process in the resonance region have been studied in the so-called pole approximation. The exact QCD-QED corrections for the production of an on-shell Z boson have been considered in Refs. [28][29][30]. In [31][32][33][34], the mixed QCD-EW effects has been obtained fully analytically for the inclusive production. In Ref. [35], the fully differential QCD-EW effects have been studied numerically. The QCD-QED corrections to the production of a pair of neutrinos have been presented in Ref. [36]. For the full neutral-current DY, results for the O(n F αα s ) contributions were presented in [37] and the complete NNLO QCD-EW corrections have been presented in [38]. Considering the reweighted two-loop virtual contributions in the pole-approximation, the QCD-EW corrections to the charged-current DY has been obtained in [39].
In this proceeding, we present the results obtained in [32][33][34]. We discuss the computational details to obtain the mixed QCD-EW corrections to on-shell Z boson inclusive production cross section in hadron colliders. The inclusion of these corrections increases the accuracy of the prediction and reduces the impact of the residual theoretical uncertainties. The fully analytic results have been expressed mostly in terms of generalized harmonic polylogarithms.

Theoretical framework
Following the parton model, the hadronic cross section for the inclusive production of an onshell Z boson is The bare PDFs (f h i (x)) describe the partonic (i) content of the hadron h, where i can be a quark, gluon or photon. The bare partonic cross sectionσ i j ≡σ(i j → Z + X ) contains initial-state collinear singularities and describes the inclusive production of a Z boson in the scattering of partons i and j. The initial-state collinear singularities are factorised universally and reabsorbed in the definition of the finite PDFs, f i (x, µ F ) = Γ ik (x, µ F )f k (x), through the mass factorisation kernel Γ ik . Following the mass factorisation, Eq. 1 can be written as follows

004.2
where the convolution is represented by the symbol ⊗, and the constant σ 0 is defined through the Born cross section σ (0,0) , z ≡ m 2 Z /ŝ, ŝ is the partonic center of mass energy, N C is the number of colors. For quark (q), c (0) q is given by where I (q) W and Q q are the third component of the weak isospin and the electric charge of the quark, respectively. θ W is the weak mixing angle.
The ∆ i j , with the constant σ 0 , is the ultraviolet (UV)-and infrared (IR)-finite partonic cross section for the partonic channel i j. In perturbation theory, ∆ i j is expanded in series in α s and α as The sum of the ∆ (m,n) i j from all the partonic channels, weighted by the corresponding PDFs, yields the total correction at a given perturbative order α m s α n . In [32,33], we have computed and studied the numerical effect of ∆ (1,1) i j . Recently in [34], we have presented the analytic expressions for ∆ (1,1) i j .

Computational details
In this section, we present the details of the calculation. We follow the generic loop computation procedure as our general strategy. The NNLO computation requires the evaluation of the Feynman loop integrals from the two-loop virtual diagrams as well as the two-and threeparticle phase-space integrals from real emissions. To be benefited from the state-of-the-art loop computation techniques, we use the method of reverse unitarity [40,41] 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 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 the standard techniques for the loop integrals. The Feynman diagrams of both the virtual and real emission processes are generated using tools like QGRAF [42] or FEYNARTS [43]. We use in-house FORM [44] and MATHEMATICA codes to interfere with conjugate amplitudes and to perform the algebraic part of the computation. With appropriate transformations of the loop momenta, we can associate each Feynman integral to an integral family, where an integral family is defined by a complete set of independent propagators needed to express all the scalar products among the loop and external momenta. The interference term is expressed in terms of a large number of scalar integrals. Integration-by-parts (IBP) identities [45][46][47] link different scalar integrals with each other linearly. By solving this algebraic system of linear relations, any scalar integral can be reduced to a small set of independent integrals, called the master integrals (MIs). To perform this reduction, we use public implementations of the IBP identities, KIRA [48], LITERED [49,50] and REDUZE2 [51,52]. The generic structure of the squared matrix element for a scattering process after the IBP reduction is where c k are the rational functions of the kinematical invariants and particle masses while the I k are the relevant MIs. The presence of two different masses (m W and m Z ) in the MIs, introduces two dimensionless scales. To minimise the complexity of the computation, we perform Taylor expansion of the W boson mass around the Z boson mass. We consider The two mass values, being numerically close to each other, allow the expansion in ζ. The required order of this expansion is determined by the phenomenological accuracy needed. After performing the expansion, the MIs depend on only one dimensionless variable z.
We use the method of differential equations [53][54][55][56][57][58][59] to solve the MIs. The MIs are functions of the variable z and the space-time dimensions d = 4 − 2 . We differentiate the MIs with respect to z followed by using the IBP identities on the differentiated output. Applying this procedure on all the MIs of an integral family, we obtain a coupled system of first-order linear differential equations, which can be solved given the boundary conditions are provided. 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 (GPLs) [60][61][62][63] and cyclotomic HPLs [64]. Three MIs in the evaluation of the double-real corrections are elliptic.
All the two-loop virtual MIs were already available in the literature [65][66][67][68][69][70][71][72]. We use the method of differential equations to recompute them. To use the method of differential equations, we compute the two-loop virtual integrals in off-shell limit of the Z boson and then take the on-shell limit of the results. The on-shell results are in terms of multiple zeta values (MZVs) and cyclotomic constants [64] which can be reduced to a set of independent constants as introduced in [73]. The off-shell virtual integrals and most of their on-shell limits have been cross-checked numerically using FIESTA [74]. The corresponding form factors were also presented in [28,72,75], with which we find perfect agreement.
In the real-emission corrections, the MIs with only massless internal lines were available from Refs. [41]. The new MIs with internal massive lines in the real-emission contributions are solved by employing the method of differential equations. However, in these systems, several square root letters appear. Generally in such scenario, we look for a variable transformation to rationalize the square root. However, in present scenario, one single transformation is not able to rationalize all the letters. We divide such system into two subsystems such that each subsystem is with a different transformation rule which rationalizes all the letters within the subsystem. By doing so, we obtain the result of an integral in terms of two different variables which are inter-related to each other. However, as all the letters are rationalized with respect to the variable of the subsystem, we obtain a compact and good alphabet. Also we observe a combination of some MIs at leading order in to provide a simple differential equation. We note that the MIs are independent in d-dimensions and hence in the higher order in , the combination does not hold true any more. However, the fact that it is true for the leading order, allows us to compute these MIs one less order in . In the following example, we explain it in detail.

004.4
We consider two integrals i , which satisfy a system of coupled differential equations.
With a transformation z → w 1 − w 2 (10) the differential equation for J 1 can be rationalized and solved to obtain up to O( −1 ) as This immediately provides us the solution for J . We obtain Consistently replacing each J (n) 2 using Eq. 11, both in the system of differential equations and in the matrix elements, we also remove the dependence of J 0 , which can be split in the sum of two terms, with or without a relation with the square root letter. We use the variable transformation to w to rationalize only for the former. We solve the non-homogeneous equation with two separate integrals, with respect to w or z, with straightforward results in terms of HPLs, in the following form This allows us the result without any square root in the letters. Thus, a certain combination of MIs can be helpful to avoid the appearance of more involved GPLs or at least postponing the appearance to the finite part of the matrix elements.

004.5
Additionally in double-real corrections, we also find three MIs (I 1 , I 2 , I 3 ) for which the coupled sub-system of differential equations is not first-order factorizable, i.e. the MIs contain elliptic polylogarithms. In the next, we present the details of the computation of these MIs.
The homogeneous part of the sub-system is the same as the one studied for the corresponding virtual diagram in Refs. [70,76], where, the results are obtained in terms of elliptic integrals of the first kind and elliptic multiple polylogarithms (eMPLs), respectively. In the present case also, the subsystem can be solved in terms of GPLs for the poles in ε and eMPLs for the finite parts. However, because of the IBP reduction, an additional factor 1/ε appears in the coefficient of these MIs, which requires the computation of these MIs up to O(ε), indicating an iteration over the eMPLs. To avoid this scenario, and also to avoid the eMPLs in the single pole of the matrix elements, we follow the similar procedure as earlier. We find the following combination of the elliptic MIs, The differential equations of I

Ultraviolet renormalization and Mass factorization
The prediction of the hadronic cross section requires to express the bare couplings and masses in terms of physical parameters via UV renormalization. The choice of the background field gauge [83] allows two separately UV finite contributions. The vertex correction with the quark wave function renormalization is UV finite. While the charge renormalization [24,25] with the renormalization of the external gauge boson creates another subset of UV finite contributions.
Thanks to the Kinoshita-Lee-Nauenberg (KLN) theorem [84,85], the sum of all the degenerate states is free of soft singularities arising from soft gluons and photons, and final-state collinear singularities arising from collinear partons. The initial-state collinear singularities are universally factorised and reabsorbed in the definition of the physical proton PDFs, as noted earlier, by means of the mass factorisation kernel Γ .

Results
Following the UV renormalization and mass factorisation, we obtain the finite partonic cross sections ∆ (1,1) i j , the analytical expressions of which are presented in [34]. To obtain the hadronic cross section, the finite partonic cross sections have further been convoluted with the physical PDFs. We have used extensively the packages GiNaC [63,86], handyG [87] and HarmonicSums [88][89][90], for numerical evaluation of all the polylogarithmic functions in the where m H and m t are the Higgs boson and the top quark masses. We write the total hadronic cross section σ H as follows: where σ i j indicates the contribution from the relative perturbative order O(α i s α j ) with respect to the Born. The NNLO QCD corrections are defined by B 1 = σ LO + σ 10 + σ 20 . We define B 2 as the NNLO QCD corrections (B 1 ) supplemented with the NLO-EW i.e. B 2 = B 1 + σ 01 . The mixed QCD-EW corrections are defined by B 3 = σ LO + σ 10 + σ 01 + σ 11 + σ 20 . We use NNPDF31_nnlo_as_0118_luxqed_nf_4 [91] proton PDF set, which is evolved with both QCD and QED kernels, to compute the B i s. We define a quantity A 1 , similar to B 1 but evaluated with the NNPDF31_nnlo_as_0118_nf_4 PDF set which is evolved with only-QCD DGLAP evolution. We use both the LHAPDF-6 [92] and HOPPET [93] to interface the PDFs. The relevant value of the strong coupling constant is α s (m Z ) = 0.1127. In the G µ -scheme, we find B 1 = 55651 pb , B 2 = 55501 pb , B 3 = 55469 pb .
The mixed QCD-EW corrections evidently reduces the uncertainty arising from the different choices of input schemes to 0.23%. The reduction in this scheme uncertainty explicitly shows the importance of the mixed QCD-EW corrections. In the G µ -scheme, which differs by −0.57% from B 3 . This essentially brings up the question of the definition of best prediction as both models are legitimate. The proton PDFs, evolved with the full QCD-EW model seems the logical choice.

Conclusion
In this proceeding, we have presented the results obtained in [32][33][34] for the computation of the mixed QCD-EW corrections to the on-shell Z boson inclusive total production cross section. We use the method of differential equation to obtain analytically all the new MIs, appearing in virtual and real-emission contributions. The proton PDFs, which are evolved using both the QCD and QED effects, are required to evaluate the hadronic cross section. The mixed QCD-EW corrections reduce the impact of two sources of theoretical uncertainties : the EW input parameters and the factorization scale variation uncertainties. The increased precision in the prediction of the inclusive Z production cross section will have an impact on the determination of the proton PDFs and on the hadron collider luminosity studies.