Master integrals for mixed QCD-QED corrections to charged-current Drell-Yan production of a massive charged lepton

The master integrals for the mixed QCD-QED two-loop virtual corrections to the charged-current Drell-Yan process $q\bar{q}^{\prime} \rightarrow \ell \nu$ are computed analytically by using the differential equation method. A suitable choice of master integrals makes it successful to cast the differential equation system into the canonical form. We keep the dependence on charged lepton mass in the building of differential equations and then expand the system against the ratio of small charged lepton mass to large $W$-boson mass. In such a way the final results will contain large logarithms of the form $\log(m_{\ell}^2/m_W^2)$. Finally, all the canonical master integrals are given as Taylor series around $d = 4$ spacetime dimensions up to order four, with coefficients expressed in terms of Goncharov polylogarithms up to weight four.


Introduction
With decades of accumulation of diligent investigations from both experimental and theoretical communities, the Drell-Yan (DY) processes [1] have become one of the best understood physics subjects at the LHC. They provide an excellent environment to probe the inner structure of proton [2] and the nature of the electroweak (EW) sector of the Standard Model (SM) by precision measurements of many EW observables, such as W -boson mass and width [3,4], Weinberg weak mixing angle [5,6] and charge asymmetries [7,8]. These are built on the huge amount of event data collected at colliders thanks to the large cross sections and clean experimental signatures of DY processes. In addition to the precision test of the SM, DY processes could help to search for signals of new physics as well because the new vector bosons predicted by conjectured extensions of the SM can be produced by a similar mechanism [9][10][11][12][13][14][15]. Precision measurements of DY processes must be confronted with equally precise theoretical predictions, therefore the theoretical perturbative descriptions of DY processes should be extended to higher orders.
At the parton level, the DY processes at the lowest-order can be categorized into neutral-current processes (qq → γ/Z → + − ) and charged-current processes (qq → W → ν). The production rates of those DY processes were computed up to the QCD next-tonext-leading order (NNLO) decades ago [16][17][18][19], which turned out to be a great triumph in precision study of phenomenology of particle physics. In Refs. [20][21][22][23][24], W and Z productions in hadron collisions were studied fully differentially at NNLO in QCD. The threshold corrections to DY production at N 3 LO in QCD were presented in Refs. [25,26]. Recently, progress has been made toward N 3 LO at the inclusive level for both neutral-and chargedcurrent DY processes [27,28], and the lepton pair rapidity distribution in the photonmediated DY process at N 3 LO in α s was calculated in Ref. [29] for the first time. Due to the leptonic final state of DY processes, the QCD corrections are only involved in the initial state. However, the EW radiative corrections are more complicated; they cannot be fully factorized into the corrections to the production and decay of the intermediate gauge boson separately. The NLO EW corrections to the neutral-and charged-current DY processes have been studied in Refs. [30][31][32][33][34][35][36] and [37][38][39][40][41][42][43], respectively. A fully differential description of NLO corrections to DY processes has been implemented in automated Monte Carlo (MC) programs. For a review on QCD and EW radiative corrections to W -and Zboson observables, see Ref. [44]. The authors presented a systematic comparison of a dozen MC codes, which served as either a general MC framework or a dedicated program for (partial) DY processes. There are also some implementations at QCD NNLO, such as GENEVA [45,46], DYTurbo (improved reimplementation of DYqT, DYRes, DYNNLO) [47], MCFM [48], MATRIX [49] and MiNNLO PS [50]. A comparison of those NNLO MC generators was provided in Ref. [51] very recently.
To build an exhaustive comprehension of DY processes at NNLO, it is necessary to consider the mixed QCD-EW corrections which are not yet available. The mixed QCD-EW NNLO corrections to DY processes can be divided into two categories: factorizable and non-factorizable. For Z-boson production in hadron collisions, the mixed QCD-QED corrections have been studied at both inclusive and exclusive levels in Refs. [52][53][54], and the mixed QCD-EW corrections are given in Refs. [55][56][57]. As for W -boson production, the mixed QCD-EW corrections can be found in Ref. [58]. The radiative corrections of O(N f α s α) to off-shell W/Z production at the LHC are also provided in Ref. [59]. The non-factorizable initial-final corrections to DY processes in the resonance region can be estimated via the so-called pole approximation [60,61]. In the precision calculation of the non-factorizable mixed QCD-EW corrections, the primary technical problem is the evaluation of massive two-loop four-point master integrals (MIs) induced by the most complicated box-type Feynman diagrams. The two-loop box-type MIs have been calculated analytically and numerically in the approximation of m W = m Z and m = 0 [62][63][64]. Recently, the mixed QCD-EW two-loop scattering amplitude for qq → + − (m = 0) in dimensional regularization and its independence on γ 5 scheme were studied in detail in Ref. [65], and the complete mixed QCD-EW NNLO corrections to the DY production of a massless lepton pair were first computed in Ref. [66].
When lepton mass is taken into consideration, the collinear photon emission off the final-state charged lepton(s) gives rise to the radiative corrections that are enhanced by large logarithms of the form log(m 2 /Q 2 ), where Q stands for a characteristic scale of DY processes, like weak gauge boson mass or center-of-mass colliding energy. As is well known, those logarithms cancel exactly if collinear lepton-photon systems are treated fully inclusively, which is guaranteed by the Kinoshita-Lee-Nauenberg theorem [67,68]. In the presence of phase-space cuts and in kinematic distributions, in general, the masssingular contributions survive, but the collinear safety can be restored after performing the photon recombination procedure [38][39][40]43]. Those mass-singular contributions can be fully captured by considering only the QED part of the EW corrections, since they are only induced by collinear photon radiation from charged leptons. The analytic calculation of the MIs for the mixed QCD-QED corrections to DY production of a massive lepton pair has been accomplished [69]. In this paper we focus on the charged-current DY process qq → ν with non-zero lepton mass, and calculate the MIs for the mixed QCD-QED two-loop virtual corrections.
The rest of this paper is organized as follows. In Section 2 the notation and conventions are declared. In Section 3, the canonical differential equations are constructed, and the 46 canonical MIs for the mixed QCD-QED two-loop virtual corrections to qq → ν are evaluated analytically. Finally, a short summary is given in Section 4.

Notation and conventions
In this paper we study the charged-current DY process The light quarks are considered as massless while the charged lepton mass is kept nonzero throughout our calculation. All the external particles are on their mass-shell, i.e., p 2 1 = p 2 2 = p 2 4 = 0 and p 2 3 = m 2 , and the scattering amplitude can be expressed in terms of the Mandelstam invariants which satisfy s + t + u = m 2 due to momentum conservation. The mixed QCD-QED two-loop virtual corrections to qq → ν can be classified as follows: • Non-factorizable corrections All the non-factorizable mixed QCD-QED two-loop Feynman diagrams are box type and irreducible. The dimensionally regularized two-loop four-point Feynman integrals in d = 4 − 2 dimensions have the form as where l 1,2 are loop momenta and the integration measure is defined by For arbitrary n i ∈ Z, i = 1, ..., 9, the set of all such F (n 1 , ..., n 9 ) in Eq.(2.3) is called an integral family. For an integral F (n 1 , ..., n 9 ), we defined its sector as [s 1 , ..., s 9 ] with s i = Θ(n i −1/2) 1 . These non-factorizable box-type two-loop Feynman diagrams can be categorized into 8 topologies. The 4 topologies shown in Figure 1 belong to the integral family F, which is identified by the set of propagators while the other 4 topologies can be obtained from Figure 1 by exchanging p 1 and p 2 and thus belong to the family F * ≡ F p 1 ↔p 2 . The MIs for the non-factorizable mixed QCD-QED two-loop virtual corrections to qq → ν will be computed analytically in Section 3. Figure 1: Four topologies of Feynman diagrams for non-factorizable mixed QCD-QED two-loop virtual corrections to qq → ν. All the diagrams in this paper are drawn by using JaxoDraw [70].

• Factorizable corrections
There are totally 24 factorizable Feynman diagrams for the mixed QCD-QED twoloop virtual corrections to qq → ν. Four of them are reducible and can be evaluated analytically by using the techniques for one-loop integrals. The rest 20 Feynman diagrams are irreducible and can be regarded as the mixed QCD-QED two-loop virtual corrections to qq → W * . These two-loop three-point Feynman diagrams can be categorized into 8 topologies which belong to 3 different integral families. A detailed discussion on the MIs for the mixed QCD-QED two-loop virtual corrections to qq → W * is given in Appendix A.
3 Non-factorizable mixed QCD-QED virtual corrections to qq → ν As stated in Section 2, the 8 non-factorizable two-loop Feynman diagrams for the mixed QCD-QED virtual corrections to qq → ν are categorized into 8 topologies respectively, which belong to 2 integral families. Due to the fact that the 4 topologies/diagrams belonging to F * can be obtained from the corresponding ones belonging to F by exchanging p 1 and p 2 , only the 4 topologies in Figure 1 are considered in the following discussion.
The 4 topologies in Figure 1 belong to the integral family F, corresponding to the following 4 sectors of F, respectively, We compute the MIs of S in the unphysical region, namely s < 0, − m 2 W < t < 0, m 2 > 0 and m 2 W > 0, where all the MIs are real-valued. Then their values in the physical region can be obtained by analytic continuation.

Canonical differential equations
In this section, we elaborate on the construction of canonical differential equations for the MIs and present the solution of the canonical differential system. The scalar Feynman integrals of a given family are not independent: there are integration-by-parts (IBP) recurrence relations [71]. After performing IBP reduction procedure, a minimal number of irreducible integrals which are also called master integrals are obtained; any dimensionally regularized scalar integral in this family can be expressed as a linear combination of the MIs with coefficients being the rational functions of kinematic variables and spacetime dimension. This is quite important for using the differential equation method [72,73] to compute scalar Feynman integrals. In this work, we adopt Kira [74,75] based on Laporta's algorithm [76] to perform IBP reduction, and then obtain 46 MIs, f 1,...,46 of S, which are depicted in Figure 2.
The mass dimension of the scalar Feynman integral F in Eq.(2.3) is given by Thus, F can be nondimensionalized by multiplying a factor of Q −[F ] , i.e., where Q is a characteristic mass scale. In this work, we take Q = m W and introduce the following three dimensionless variables for convenience, Figure 2: A set of pre-canonical MIs. The thin and thick lines represent massless and massive external particles and propagators respectively (red for W boson and blue for charged lepton). One dot means the index of that propagator is raised to 2 and two dots go for 3.
Then the nondimensionalized scalar Feynman integrals are functions of x, y and z.
With the help of LiteRed [77,78], we obtain the partial derivatives of the 46 MIs with respect to x, y and z. By means of IBP identities, those derivatives are expressed as the linear combinations of the 46 MIs. This leads to the following set of partial differential equations for the vector of MIs F ≡ (f 1 , f 2 , ..., f 46 ) T , It was suggested in Ref. [79] that a good choice of MIs would simplify the calculation tremendously. By means of the following basis transformation T : F −→ G we obtain a canonical set of MIs G = (g 1 , g 2 , ..., g 46 ) T . Such a canonical basis G obeys the following system of differential equations, where the dependence on the dimensional regularization parameter is completely factorized from the kinematics. Since the canonical differential system in Eq.(3.9) is invariant under the scaling transformation G −→ a G, the canonical basis G defined in Eq.(3.8) has been required to be analytic and nonzero at = 0. Finally, it is worth mentioning that all the entries of the basis transformation matrix T are rational functions of x, y and z, except (3.10)

Integrations
The partial differential equations in Eq.(3.9) can be combined into a d log-form total differential equation, where C i (i = 1, ..., 15) are constant matrices with rational-number entries, and the arguments of the d log's, are called letters and together constitute the alphabet of our problem. Considering the mass hierarchy between W boson and charged lepton (m m W ), we may expand the differential system for small z. After the expansion for small z, the coefficient matrix dA can be expressed as where (3.14) and A n (x, y) (n ∈ N + ) are completely comprised of rational-function entries 2 . Therefore, the matrices A x,y (x, y, z) and A z (x, y, z) in Eq.(3.9) are analytic and singular at z = 0, respectively, and can be expanded as In the lowest-order approximation, the 7 letters α i (i = 1, ..., 7) constitute the alphabet of the differential system. As is well known, the leading terms in the z-expansion in Eq.(3.16) (A x,0 , A y,0 and A z,−1 ) are responsible for the logarithmic lepton-mass singularities of the canonical MIs G. Thus, the logarithmically divergent contributions of the form log n z can be fully captured by adopting the lowest-order approximation. In this work, we adopt the second-order approximation, i.e., (3.18) in order to get the analytic expression for G up to O(z log n z). In the unphysical region,

19)
α i (i = 1, ..., 7) are real and positive. To evaluate the MIs in the physical region, the same analytic continuation procedure as in Ref. [62] should be applied.
Since G is finite in the limit → 0, it can be expanded in a Taylor series around four dimensions, (3.20) From Eq.(3.9), one can establish the differential equations connecting G (n) and G (n−1) . These first-order differential equations can be integrated in an iterative manner, namely order by order in . Then the canonical MIs at O( n ) can be expressed as where C (n) is the integration constant vector at O( n ) and the matrices M (n) (n ∈ N) are defined recursively by with M (0) = 1. Actually, it is sufficient to expand G up to the order of 4 for the mixed QCD-QED two-loop virtual corrections to charged-current DY processes. Since the entries of A x, y, z are all rational functions of x, y and z after the expansion for small z, M (n) can be expressed in terms of rational functions and Goncharov polylogarithms (GPLs) [80,81]. A GPL with weights w i (i = 1, ..., n) and argument t is defined recursively by with G( ; t) = 1 and For the MIs G involved in the calculation of the non-factorizable mixed QCD-QED twoloop virtual corrections to qq → ν, the weights of the GPLs with arguments x, y and z are drawn from the following three sets, respectively, In the recursive definition for M (n) in Eq.(3.22), one encounters the integrals with integrands of the form r(t) G( w; t), where r(t) and G( w; t) represent rational functions and GPLs, respectively. Those integrals can be calculated recursively by using the following IBP recurrence relation, where w n = (w n , ..., w 1 ), w n−1 = (w n−1 , ..., w 1 ), and R(t) is an antiderivative of r(t). In this paper we employ the Mathematica package PolyLogTools [82][83][84] and C++ library GiNaC [85] for both symbolic computation and numerical evaluation [86] of GPLs.

Boundary conditions
To obtain the exact solution of the canonical differential system in Eq.(3.9), the integration constants (i.e., the boundary terms of the MIs g 1,...,46 ) have to be specified. For some of the MIs, we take their known analytic solutions from the literature or evaluate them by means of some much simpler differential systems. For all other MIs, the boundary terms are fixed by imposing a set of regularity conditions, i.e., demanding regularity of those MIs or their linear combinations in certain kinematic limits. It was found that the weight structure of integration constants in Ref. [69] is quite simple. The integration constants at O( n ) (n = 0, 1, 2, 3, 4) are proportional to 1, 0, ζ(2), ζ(3) and ζ(4), respectively. This property still holds in our case. For our problem, the conditions imposed to the MIs for the determination of their boundary constants are listed below.
• g 13,18 are determined by matching against the MIs with full lepton mass dependence, which are computed in Appendix B.
• The boundary constants of the rest 29 MIs are fixed by the following linearly independent regularity conditions 3 : regularity at y → 1: g 16,23 regularity at z → 0 and matching against their massless (m = 0) counterparts: The boundary conditions for h a,b,c are determined from the m → 0 behaviour of these integrals, where all the internal lepton propagators become massless. Around the z = 0 singularity, the differential equation for the z-evolution reduces to

∂G ∂z
To figure out the m → 0 behaviour of the 46 MIs, we perform a Jordan decomposition of the pole matrix A z,−1 [87], i.e., where J is the Jordan canonical form of A z,−1 , 30) and the similarity transformation S reads Then the z-evolution around the z = 0 singularity can be described in Jordan form, where H = (h 1 , h 2 , ..., h 46 ) T is the Jordan canonical basis defined by H = S G. Considering the Jordan block structure of J, the singular behaviour of H in the limit z → 0 can be simply factored out, and the solution of Eq.(3.32) can be expressed as and The regularity of RH at z → 0 can be applied to the following two aspects.
1. Boundary conditions:  (3.38) The matching of h a,b,c to their massless counterparts, i.e., Eq.(3.38), provides three boundary conditions for G.

Consistency checks:
R H: finite in the limit z → 0 (3.39) Even though z = 0 is a singularity of H, R H is finite when z → 0. Thus, the finiteness of the following 14 combinations of H in the zero lepton mass limit provides a consistency check for the analytic solution of the canonical differential system (3.9), z 2 h 1,2,3,5 , z h 6,7,8,9,10,11,12,13,14 ,

Solution and checks
We obtain the solution of the canonical differential system (3.9) Table 1. It shows that all these numerical results are in good agreement with each other within the calculation errors.
-As stated in Section 3.3, h i (i = 15, ..., 46) are finite in the limit z → 0, and their massless counterparts can be expressed as the linear combinations of I 1 , I 2 , ..., and I 31 . The correctness of Eq.(3.36) was confirmed by employing IBP recurrence relations and, consequently, the explicit expression of the 32 × 31 matrix K was obtained.
h i (i = 1, ..., 14) are divergent as m → 0, but the 14 combinations of them defined by Eq.(3.40) are finite in the zero lepton mass limit. The asymptotic behaviour of these combinations in the limit z → 0 was checked both analytically and numerically, which clearly demonstrate the regularity of these combinations at z = 0.

Summary
The DY production of dilepton via Z and single lepton via W in hadron-hadron collisions is one of the cutting-edge topics in the field of high energy physics. In this paper we present an analytic calculation of the MIs for the mixed QCD-QED two-loop virtual corrections to the charged-current DY process qq → ν. The scalar Feynman integrals involved in the non-factorizable corrections can be expressed in terms of the linear combinations of 46 MIs, which satisfy a set of first-order differential equations. By means of a suitable basis transformation, we obtain a canonical set of MIs. The -dependence of the canonical differential system is factorized from the kinematics. Therefore, the canonical MIs, as the solution of the canonical differential system, can be cast in a Taylor series around d = 4, with coefficients written in terms of Chen's iterated integrals. After an expansion for small z, all the entries of the coefficient matrices are simple rational functions, and thus the Chen's iterated integrals can be expressed in terms of GPLs. The boundary constants of these canonical MIs are fixed by using known or simpler integrals as independent input and by requiring the regularity of the solution at some special kinematic limits. We solve the canonical differential system in the second-order approximation and then obtain the Taylor series expansions of the 46 canonical MIs around d = 4 up to O( 4 ), which contain all the logarithmic terms of the form z m log n z (m = 0, 1, n = 0, 1, 2, 3, 4). As for the G Analytic O(z log n z) Analytic O(z 2 log n z) pySecDec   A Mixed QCD-QED virtual corrections to qq → W * The 20 two-loop Feynman diagrams for the mixed QCD-QED virtual corrections to qq → W * can be categorized into 8 topologies belonging to 3 vertex integral families F 1,2,3 . Some representative topologies of these Feynman diagrams are depicted in Figure 3.
• Family F 1 : Figure 3: Three representative topologies of Feynman diagrams belonging to F 1 , F 2 and F 3 respectively for mixed QCD-QED two-loop virtual corrections to qq → W * .
3 topologies ⊂ F 1 : can be reduced to 3 MIs f 5 , f a and f b , where f a and f b are depicted schematically in Figure 4.
• Family F 2 : • Family F 3 :  8) is a three-point non-planar topology, which is represented schematically by the non-planar graph in Figure 3. The two-loop scalar integrals induced by the non-planar Feynman diagrams via tensor reduction belong to S 3 , 1, 0, 1, 1, 1, 1] , (A.9) and can be reduced to a set of 3 MIs f 5 , f b and f c , where f c is a non-planar two-loop scalar integral defined by the last graph in Figure 4. In summary, there are totally 12 MIs for the mixed QCD-QED two-loop virtual corrections to qq → W * , which can be chosen as g 2 , g 3 , g 5 , g 9 , g 11 , g 12 , g 14 , g 15 , g 28 Figure 4), their analytic expressions in terms of GPLs up to weight 4 can be obtained from the literature [62,90].

B Auxiliary integrals
In this appendix we discuss the computation of the auxiliary vertex integrals used to extract the boundary constants of g 13 and g 18 in Section 3.3. The two-loop vertex integral family considered here is identified by the set of propagators with external momenta p 1 , p 2 and p satisfying Now we focus on the sector [1, 1, 1, 1, 1, 0, 0] (depicted in Figure 5) of this integral family.
It is obvious that f 13 ∈ [1, 0, 1, 1, 1, 0, 0] where The same strategy as in Section 3.2 is applied here to solve the canonical differential system in Eq.(B.7). To determine the unknown integration constants, the four simple integrals γ 1,2,3,4 serve as input to this system [62,69]. As for γ 5,6,7 , the integration constants are fixed by the regularity at κ → −1. Then we obtain the solution of the canonical differential system (B.7) in terms of GPLs of argument κ 6 , By the definition of G and Γ we can see that (g 1 , g 5 , g 10 , g 6 , g 18 , g 19 , g 13 ) T = x −2 Γ with κ = z/x . (B.10) Thus, the boundary constants of g 13 and g 18 are given by the constant components of γ 7 and γ 5 , respectively, (B.11) 6 We have also calculated Γ by using the Maple program HyperInt [91], and obtained the same results as in Eq.(B.9).

C Asymptotic behaviour of G
We analyze the asymptotic behaviour of G in the limit α i → 0 (i = 1, ..., 7) by using the code asy.m [92,93] included in FIESTA [94,95] and summarize it in Table 2.
W and m 2 → 0 correspond to α i → 0, i = 1, ..., 7, respectively. The chechmark means the MI is analytic, while the cross singular, in the corresponding limit.

D Explicit expressions of G
The analytic expressions of the 46 canonical MIs g i (i = 1, ..., 46) in terms of GPLs up to O( 3 ) in the lowest-order approximation are listed as follows: