Two-Loop Master Integrals for the mixed EW-QCD virtual corrections to Drell-Yan scattering

We present the calculation of the master integrals needed for the two-loop QCDxEW corrections to $ q + \bar{q} \to l^- + l^+$ and $ q + \bar{q}' \to l^- + \overline{\nu} \, , $ for massless external particles. We treat W and Z bosons as degenerate in mass. We identify three types of diagrams, according to the presence of massive internal lines: the no-mass type, the one-mass type, and the two-mass type, where all massive propagators, when occurring, contain the same mass value. We find a basis of 49 master integrals and evaluate them with the method of the differential equations. The Magnus exponential is employed to choose a set of master integrals that obeys a canonical system of differential equations. Boundary conditions are found either by matching the solutions onto simpler integrals in special kinematic configurations, or by requiring the regularity of the solution at pseudo-thresholds. The canonical master integrals are finally given as Taylor series around d=4 space-time dimensions, up to order four, with coefficients given in terms of iterated integrals, respectively up to weight four.


Introduction
The Drell-Yan production of Z and W bosons [1] is one of the standard candles for physical studies at the LHC. Due to the big cross section and clean experimental signature, Drell-Yan processes can be measured with small experimental uncertainty and, therefore, allow for very precise tests of the Standard Model of fundamental interactions (SM). They give access to the determination of important parameters of the weak sector, as for instance the sine of the weak mixing angle and the W boson mass, that together with the top and the Higgs masses provides stringent constraints on the validity of the SM at the TeV energy scale. Furthermore, Drell-Yan processes constitute the SM background in searches of New Physics, involving for instance new vector boson resonances, Z and W , originating from GUT extensions of the SM. Finally, the Drell-Yan mechanism is used for constraining parton distribution functions, for detector calibration and determination of the collider luminosity. For all these reasons, an accurate and reliable experimental and theoretical control on Drell-Yan processes would be of the maximum importance for future physics studies at colliders.
The theoretical description of Drell-Yan processes currently includes NNLO QCD and NLO EW radiative corrections, implemented in flexible tools able to provide predictions for inclusive observables as well as kinematic distributions. Current theoretical predictions are in good agreement with the experimental measurements. However, higher theoretical accuracy is needed in order to match the future experimental requirements, in particular in view of the run II of the LHC. A consistent part of an increasing theoretical accuracy regards higher-order perturbative corrections.
Very recently, NNNLO QCD corrections were calculated for the Higgs total production cross section in gluon-gluon fusion [2,3]. The residual factorization/renormalization scales variation moved from about 10-15% of the NNLO calculation (supplemented by NNLL resummation) to about 5% of the current result. These results will be applied to Drell-Yan as well, since they involve the evaluation of the same topologies for the calculation of the corresponding Feynman diagrams [4][5][6][7][8][9].
At the same order of accuracy (one can roughly thing to exchange two powers of α S with one power of α), the mixed QCD-EW corrections have to be taken into account. As in the case of QCD NNLO with EW NLO perturbative corrections, the mixed QCD-EW corrections are expected to become of similar size with respect to QCD NNNLO at high leptonic invariant mass [10].
At LO, the partonic process in the SM is mediated by the exchange of a photon or a Z/W vector boson, in the s annihilation channel: qq → γ, Z → l − l + and qq → W → lν.
At higher orders in the coupling constants, we can distinguish between QCD and electroweak (EW) or mixed (EW-QCD) corrections to the LO process. In the first case, only the initial state receives quantum corrections, since the leptonic final state does not couple to gluons.
The NLO QCD corrections to the total cross section were calculated in [11,12] and revealed a sizable increase of the cross section with respect to the LO result. The NNLO QCD corrections [13,14] stabilized, then, the convergence of the perturbative series.
QCD fixed-order corrections to the total production cross section are supplemented by the resummation of soft-gluon logarithmically enhanced terms, up to NNNLL approximation [15][16][17][18].
EW quantum corrections allow exchanges of quanta between initial and final states. Therefore, already at the NLO, massive four-point functions have to be evaluated. Since the bulk of the corrections for inclusive observables comes from the resonant region, in which the exchanged vector boson is nearly on-shell, electroweak NLO corrections to the total cross section were calculated for the W [19] and Z [20] in narrow-width approximation.
More exclusive observables are known in the literature. The Z and W production at non-zero transverse momentum p T is known at the NLO in QCD [21][22][23][24][25][26] and in the full SM [27]. The two-loop QCD helicity amplitudes for the production of a Z or a W with a photon have also been calculated [28]. For small p T (p T m W , m Z ) the convergence of the fixed-order calculation is spoiled by the large logarithmic terms α n S log m (m 2 W /p 2 T ) that have to be resummed [29][30][31][32][33][34][35][36][37][38][39]. Finally, the rapidity distribution of a vector boson is known at the NNLO in QCD [40].
The NLO corrections are available in a fully differential description. They are implemented in flexible NLO Monte Carlo programs, and merged with QCD parton shower in MC@NLO [41] and POWHEG [42]. In [43], the NLO EW and the QED multiple photon corrections are combined with NLO QCD corrections and parton shower. Pure QED generators are also available [44][45][46][47]. Although these implementations provide an accurate description of the process and allow for realistic phenomenological studies at the hadronic level, they are not accurate enough for the performances of the run II at the LHC. The NNLO results mentioned above, however, are widely inclusive and they cannot provide realistic descriptions, that necessarily have to include experimental cuts. Therefore, a fully differential description of the Drell-Yan process at the NNLO is needed. With this respect, the state of the art is represented by the two programs FEWZ [48], that includes also EW NLO corrections [49], and DYNNLO [50,51]. In these two programs, the decay products of the vector boson, the spin correlations and the finite-width effects are also taken into account. . virtual corrections to the two-loop 2 → 2 processes: q +q → l − + l + , and q +q → l − + ν , for massless external particles. The masses of the W and Z bosons are numerically close to each other, in fact ∆m 2 ≡ m 2 Z − m 2 W m 2 Z . Therefore, in the diagrams containing both Z and W propagators at the same time, one can perform a series expansion in ξ ≡ ∆m 2 /m 2 Z ∼ 0.25. Within this approximation, all topologies appearing in the twoloop QCD×EW virtual corrections to Drell-Yan scattering shall contain either no internal massive line, or one massive propagator with mass m W , or two massive propagators with the same mass m W [56]. Should they be needed for achieving higher accuracy within the virtual amplitudes, the coefficients of the series in ξ correspond to scalar integrals with higher powers of the denominators.
Using the code Reduze 2 [57,58], the dimensionally regulated integrals involved in the calculation are reduced to a set of 49 MIs, which are later determined by means of the differential equations method [59][60][61], reviewed in [62,63]. Of those 49 MIs, 8 contain only massless internal lines, 24 involve one massive line and 17 involve two massive lines. The system of differential equations obeyed by the MIs is cast in a canonical form [64], following the algorithm based on the use of the Magnus exponential, introduced in [65, 66] 1 . Boundary conditions are retrieved either from the knowledge of simpler integrals emerging from the limiting kinematics, or by requiring the regularity of the solution at pseudothresholds.
Finally, the canonical MIs are given as Taylor series in (= (4 − d)/2), up to order 4 , being d the dimensional regularization parameter. The coefficients of the series are pure functions, represented as iterated integrals with rational and irrational kernels, up to weight four. The solution could be expressed in terms of Chen's iterated integrals. Alternatively, we adopt a mixed representation, where, when possible, we make explicit the presence of Goncharov polylogarithms (GPLs) [70,71], also within the nested integration structure. This representation is suitable for the numerical evaluation of our solution.
While the two-loop four-point integrals with massless internal lines are well known in the literature [61,[72][73][74], the four point integrals with one and two massive internal lines considered here are new and represent the main result of this communication. We verified the numerical agreement of the MIs in the unphysical region against the results of SecDec [75][76][77]. In particular, because of the presence of irreducible irrational weight functions, we found it convenient to cast 5 of the 17 MIs with two massive internal lines as one-dimensional integral formulas [78], involving GPLs in the integrands. The numerical evaluation of our solutions can, therefore, be performed with the help of the GiNaC library [79] for the evaluation of GPLs.
The article is structured as follows. Section 2 contains our notation and conventions. In Section 3, we discuss the solution of the canonical differential equations in terms of Chen's iterated integrals. In section 4, we explicitly present the system of differential equations and the solutions for the one-and two-loop MIs that contain one massive propagator. In Section 5, we give the system of differential equations for the one-and two-loop MIs containing two massive propagators. Conclusions are given In Section 6. In Appendix A, we discuss the kinematic domain of our analytic results. In Appendix B, we give the matrices of the system of differential equations in canonical form.
Our results are collected in ancillary files, that we include to the arXiv submission.

Notations and Conventions
In this paper we study the two-loop corrections to the following partonic scattering processes: The external particles are considered mass-less and they are on their mass-shell, p 2 1 = p 2 2 = p 2 3 = p 2 4 = 0. The scattering can be described in terms of the Mandelstam variables in such a way that, for momentum conservation, we have s + t + u = 0. The physical region is defined by s > 0 , t = − s 2 (1 − cos(θ)) , (2.4) where θ is the scattering angle in the partonic center of mass frame, lying in the range 0 < θ < π. Therefore, while s > 0, t is always negative and −s < t < 0. The quantum corrections to the processes (2.1) and (2.2) can be expanded in power series of the coupling constants. At one loop, the QCD corrections consist on the exchange of a virtual gluon between the initial-state quarks. The final state is not affected, and at most mass-less three-point functions have to be evaluated. The EW corrections, instead, consist on the exchange of photons, Z and W bosons. Moreover, these quanta can be exchanged between the quarks in the initial state as well as the leptons in the final state, but they can also be exchanged between a quark in the initial state and a lepton in the final state. Consequently, in the calculation of the one-loop corrections one has to evaluate massive box and vertex diagrams. In the process of qq → lν one has to evaluate diagrams in which a Z and a W bosons are exchanged simultaneously. In order to reduce the number of scales present in the calculation, we expand the Z propagators around m W : is the effective parameter of the expansion. The coefficients of the series in ξ are Feynman diagrams with the same masses, and eventually with increased powers in the expanded denominator. Such diagrams depend only on s, t, and one mass m = m W . However, this does not cause any problem in the calculation, since diagrams with higher powers of the propagators are in any case reduced to the same set of MIs. For phenomenological purposes the first order in ξ might be sufficient, but in principle any order in ξ can be calculated without effort, just relying on the reduction procedure. We apply the same approximation to the two-loop diagrams as well.
We calculate the quantum corrections to the processes (2.1,2.2) using a Feynman diagrams approach. After considering the interference with the leading order, and summing over the spins and colors, we express the squared absolute value of the amplitude in terms of dimensionally regularized scalar integrals. These integrals are reduced to a set of MIs by means of integration-by-parts identities [80,81] and Lorentz-invariance identities [61], implemented in the computer program 2 Reduze 2 [57,58].
At one-loop, the topologies involved in the QCD and EW corrections are shown in figure 1, where we distinguish: a) the mass-less case; b) the exchange of one massive particle; and c) the exchange of two massive particles.
At two-loop, the topologies required by the O(αα S ) corrections are only planar. They are shown in figure 2. As for the one-loop case, we consider three classes of diagrams, according to the presence of massive particles.
Topologies a 1 ) and a 2 ) belong to the same 9-denominators mass-less topology. They reduce to 8 MIs, that were already known in the literature [61,[72][73][74]. Topologies b 1 )-b 3 ) have one massive propagator. They reduce to 31 MIs out of which 24 contain one massive propagator and 7 are part of the MIs for topologies a 1 ) and a 2 ). The three-point functions were already known in the literature [88][89][90]. The four-point functions are calculated and presented here for the first time. Topologies c 1 ) and c 2 ) have two massive propagators and they reduce to 36 MIs, out of which 17 contain two massive propagators, 15 contain one massive propagator (and they are included in the set of MIs for topologies b 1 )-b 3 )) and 4 contain only massless propagators. The three-point functions were known in the literature [91,92] and the four-point functions are presented here for the first time.
The routings for one-and two-mass topologies, at the one-and two-loop level, can be defined in terms of the following sets of denominators D n , where k i (i = 1, 2) are the loop momenta, and p i (i = 1, . . . , 4) are the external momenta: • One-mass topologies. For the one-loop one-mass integrals (figure 1 b), we have: At two loops (figure 2 b 1 -b 3 ), instead, we have: (2.7) • Two-mass topologies. For the one-loop two-mass integrals (figure 1 c), we have: At two loops (figure 2 c 1 and c 2 ), instead, we have: (2.8) In the following we consider -loop Feynman integrals in d dimensions, built out of p of the above denominators, each raised to some integer power, of the form where the integration measure is defined as with µ the 't Hooft scale of dimensional regularization, and In eqs. (2.10), (2.11) we used = (4 − d)/2.

System of Differential equations for Master Integrals
In this section, we describe the general structure of the systems of differential equations obeyed by the MIs, and the corresponding solutions. Sections dedicated to the one-mass and two-mass MIs will follow, where the details of their complete determination will be provided.
The b-and c-type MIs are functions of the Mandelstam invariants defined in eq. (2.3) and of the mass m. For their evaluation it is convenient to define the dimensionless ratios The b-type and c-type MIs obey systems of partial differential equations in x and y, which can be combined into matrix equations for their total differentials. In general, the vector of MIs F is solution of the following system of differential equation, where the matrix K depends both on the kinematic variables and on the spacetime dimension d = 4 − 2 .
By means of a suitable basis transformation, built with the help of the Magnus exponential [65,93] following the procedure outlined in Sec. 2 of [66], we obtain a canonical set of MIs [64]. Such a basis obeys a system of differential equation where the dependence on is factorized from the kinematics. Moreover, the coefficient matrices can be assembled in a (logarithmic) differential form, referred to as canonical d log-form. Hence, the canonical basis I obeys the following system of equations, where dA is the d log matrix written in terms of differentials d log η i (that contain the kinematic dependence) and coefficient matrices M i (with rational-number entries). The integrability conditions for eq. (3.3) read

General solution
The general solution of the canonical system of differential equations (3.3) can be compactly written at a point x = (x 1 , x 2 ) = (x, y) as where I( , x 0 ) is a vector of arbitrary constants, depending on , while dA depends only on the kinematic variables. In the above expression, the path-ordered exponential is a short notation for the series in which the line integral of the product of k matrix-valued 1-forms dA is understood in the sense of Chen's iterated integrals [94] (see also [95] and the pedagogical lectures [96]) and γ is a piecewise-smooth path such that γ(0) = x 0 and γ(1) = x. It follows from Chen's theorem that the iterated integrals in eq. (3.7) do not depend on the actual choice of the path, provided the curve does not contain any singularity of dA and it does not cross any of its branch cuts, but only on the endpoints. In this sense, if one fixes x 0 and lets x vary, eq. (3.6) can be thought of as a function of x. In the limit x → x 0 , the line shrinks to a point and all the path integrals in eq.(3.7) vanish, so that I( , x) → I( , x 0 ), i.e. the integration constants have a natural interpretation as initial values, and the path-ordered exponential as evolution operator. We assume that the vector of MIs at any point I( x) is normalized in such a way that it admits a Taylor series in : The solution I( , x) is then in principle determined through (3.6) at any order of theexpansion, and reads (up to the coefficient of 4 ) 14) The problem of solving (3.3), given a set of initial conditions I( x 0 ), reduces therefore to the evaluation of matrices of the type γ dA . . . dA k times , (3.15) whose entries, due to (3.4), are linear combinations of Chen's iterated integrals of the form We refer to the number of iterated integrations k as the weight of the path-integral. The empty integral (eq. (3.16) for k = 0) is defined to be equal to 1. We stress that only the matrices (3.15) do not depend on the explicit choice of the path. The individual summands of the form in eq. (3.16), which contribute to their entries, in general depend on such a choice. To keep the notation compact, we define which also emphasizes that the iterated integrals in (3.16) are in general functionals of the path γ.

Properties of Chen's iterated integrals
The general theory of iterated path integrals was developed by Chen [94]. Chen's iterated integrals satisfy a number of properties that we summarize for completeness: • Invariance under path reparametrization. The integral C [γ] i k ,...,i 1 does not depend on how one chooses to parametrize the path γ.
• Reverse path formula. If the path γ −1 is the path γ traversed in the opposite direction, then • Recursive structure. From (3.16) and (3.17) it follows that the line integral of one d log is defined as usual 20) and only depends on the endpoints x 0 , x It is convenient to introduce the path integral "up to some point along γ": given a path γ and a parameter s ∈ [0, 1], one can define the 1-parameter family of paths .
• Path composition formula. If α, β : [0, 1] → M are such that α(0) = x 0 , α(1) = β(0), and β(1) = x, then the composed path γ ≡ αβ is obtained by first traversing α and then β. One can prove that the integral over such a composed path satisfies ip,...,i 1 . (3.27) • Integration-by-parts formula. In order to compute the path ordered integral of k d log forms using the definition, eq. (3.16) (or, equivalently, eq. (3.24)), in principle one would have to perform k nested integrations. When a fully analytic solution cannot be achieved, numerical integration can as well be employed. Therefore one can use an alternative form of the Chen iterated integral suitable for the combined use of analytic and numerical integrations. In fact, we observe that the innermost integration can always be performed analytically using (3.20), so that only k − 1 integrations are left. For instance, in the case k = 2, For k ≥ 3, one can proceed recursively using eq. (3.24), assuming that the numerical evaluation up to the first k − 1 iterations is a solved problem. Using integration by parts, one can show that the numerical integration over the outermost weight g k can actually be avoided, leaving only k − 2 integrations to be performed i k−2 ,...,i 1 dt . (3.29)

Mixed Chen-Goncharov representation
In principle eq. (3.6) completely determines the solution, which can be written in terms of Chen's iterated integrals along an arbitrary piecewise-smooth path (see the discussion below eq. (3.6)). The initial conditions I( x 0 ) can be computed analytically, if possible, or by means of numerical methods. The number of iterated integrals that have to be evaluated numerically can be minimized by the use the of algebraic identities relating them. According to the discussion in section 3.2, the evaluation of the solution up to weight 4 requires in general 2 nested numerical integrations. In order to obtain results that allow for an efficient numerical evaluation, we have chosen to give the solution in a mixed representation that involves GPLs and general Chen's iterated integrals. The representation in terms of GPLs is particularly convenient because public packages exist, like GiNaC, that implement their numerical evaluation in a fast and accurate way. Whenever the alphabet is rational in the kinematic variables x i , one can always choose a path that allows to express the Chen iterated integrals in terms of GPLs, namely the broken line such that, in each segment, only one of the x i is allowed to vary. Along each segment, by means of factorization over the complex numbers, one can obtain a linear alphabet and, therefore, the GPLs representation. This approach is equivalent to integrating the differential equations for x and y separately. By integrating, say, the equation in x one obtains the solution in terms of GPLs of argument x up to an unknown function H(y). By taking the derivative with respect to y and matching to the equation in y, one obtains a differential equation for H(y). The latter can be again integrated in terms of GPLs of argument y, up to a constant.
As we will discuss in section 5, the alphabet for our differential equations is not always rational in the kinematic variables we use and, in that case, a representation in terms of GPLs cannot be given for the complete solution. To reach the mixed representation, we have exploited the property of path-independence of the coefficients of the -expansion of the solution eq. (3.6). In particular, eqs. (3.11)-(3.14) can be written in an equivalent alternative form using eq. (3.24): where x t is the point (x(t), y(t)) along the curve identified by γ. We see that, in order to build the weight-k coefficient, one must perform a path integration over the weight-(k − 1) coefficient. The choice of such path is independent of the path used to compute the former because, as we have already discussed, each coefficient is a function of the sole endpoints. In other words, as far as the weight-k coefficient of the solution is concerned, we are free to choose the integration path independently for each of the k integrations (for each component of I( x)).
To see how this can be useful in our computation, we note that the letters η i (in suitable variables, say x) can be grouped in two classes. The first contains the letters that are rational in the components of x and happens to represent the alphabet for most of the MIs we need to compute. The second is the class of letters that are non-rational functions of the variables. The two classes together constitute the alphabet for the 5 most complicated MIs.
Starting from the weight-1 coefficient of the solution, we proceed as follows. As far as the involved η i 's belong to the first class of letters, we can express the solution in terms of GPLs. We keep integrating in this manner until, at some weight k, the solution begins to involve non-rational η i 's. At this point we proceed with the path integration as in (3.30). Within this approach, the weight k − 1 solution is not expressed in terms of Chen's iterated integrals over an arbitrary path, but in terms of GPLs. We introduce the following notation to keep our results compact: where G γ m (x) and G γ n (y) stand for the GPLs G m (x) and G n (y) evaluated at (x, y) = (γ 1 (t), γ 2 (t)).

Constant GPLs
In the determination of the boundary values of the MIs we encountered constant GPLs of argument 1 with weights drawn from three sets. For the one-mass MIs there is only one relevant set, with four weights, . For the two-mass MIs we encountered the following two sets, with seven weights each where the former includes the third roots of −1 and the latter involves a subset of the sixth roots of −1. With the help of GiNaC, we verified that, at order k , the Taylor coefficient of each MI I (k) i contains a combinations of constant GPLs that turns out to be proportional to ζ k , namely amounting to q i,k ζ k , with q i,k ∈ Q. The resulting identities were verified at high numerical accuracy. As examples, we show, where for simplicity we omitted the argument (x = 1) of the GPLs and we defined the weight r ≡ (−1) 1/3 . For related studies see also [97][98][99][100].

One-mass Master Integrals
In this section we describe the computation of the MIs with one internal massive line, namely topology (b) of figure 1 and topologies (b 1 )-(b 3 ) of figure 2.

One-loop
The following set of MIs for the one-loop one-mass box obeys a differential equation in x and y, defined in eq. (3.1), which is linear in : Figure 3: One-loop one-mass MIs T 1,...,5 . Thin lines represent massless external particles and propagators; thick lines stand for massive propagators; an horizontal (vertical) dashed external line represents an off-shell leg with squared momentum equal to s (t); dots indicate squared propagators.
where the T i are depicted in figure 3. By means of the Magnus exponential [65,93], according to the procedure outlined in Sec. 2 of [66], we obtain the canonical MIs The alphabet of the corresponding d log-form, eq (3.4), is and the coefficient matrices read If x > 0 and 0 < y < 1 all the letters η i are positive. Since the alphabet is linear in x and y, according to the discussion in section 3.3, the solution can be conveniently cast in terms of GPLs.
Instead of choosing a particular basepoint x 0 , the integration constants of I 2...5 can be easily fixed by demanding regularity at the pseudothresholds t → −m 2 , u → 0, s → 0 and their reality in the euclidean region. On the other hand, I 1 is a constant and must be determined by direct integration:

Two-loop
At the two-loop order, the following set of MIs admits -linear differential equations in x and y (defined in eq. (3.1)): where the T i are depicted in figure 4. Once again, by means of Magnus exponentials, we are able to obtain a canonical basis: where λ ± = m 2 ± s . After combining the two differential equations into one total differential, we find a d log-form (3.4) with the alphabet which includes the additional letter η 6 as compared to one-loop (4.3). If x > 0 and 0 < y < 1 all the letters η i are positive. The coefficient matrices are given in the appendix (B.1). Since the additional letter is multilinear in x and y, also at the two-loop order we are able to obtain the solution in terms of GPLs (see the discussion in section 3.3). We hereby list the conditions imposed to which integrals for determining their boundary constants; • regularity at t → −m 2 and u → 0 and imposing reality on the resulting boundary constants: I 2,...,5,7... 10,12,14...17,19...31 , • limit s → 0: I 11,13 , • limit t → 0: I 18 . Figure 5: One-loop two-mass MIs T 1,..., 6 . The conventions are as in figure 3.
This leaves us with I 1,6 , to be determined by direct integration: . (4.10) Owing to the explicit representation in terms of GPLs, all the one-mass MIs can be computed in the whole (s, t) domain (see appendix A). Our results have been successfully checked against SecDec.
The analytic expressions of all the MIs are explicitly given in electronic form in ancillary files that can be obtained from the arXiv version of this paper.

Two-mass Master Integrals
In this section we describe the computation of the MIs with two internal massive lines, namely topology (c) of figure 1 and topologies (c 1 )-(c 2 ) of figure 2.

One-loop
We choose the following set of MIs, admitting a differential equation linear in where the T i are shown in figure 5. After applying the Magnus transformation we obtain the following canonical basis 2) The alphabet of the corresponding canonical d log-form, (3.4), is non-rational in s, t and m. In particular four square roots appear The latter can be rationalized through the change of variables We note that the above mapping is not invertible at s = 4m 2 . In terms of w and z, the alphabet reads 5) and the coefficient matrices are In the region 0 < w < z < 1 all the letters η i are positive. For a detailed discussion of how the interesting regions in the s, t plane are mapped to the C × C space of the w, z variables, see appendix A. The alphabet in (5.5) is linear in z but contains letters quadratic in w. As the latter can be linearized by factorization over the complex numbers, we are once again able to express the solution in terms of GPLs (see the discussion in section 3.3).
The integration constants of I 4,5,6 can be fixed by requiring their regularity at the pseudothresholds s → 0, t → −m 2 and u → 0. The boundary constant of I 2 can be fixed by taking the s → 0 limit. This leaves us with two integrals, I 1,3 , to be determined by direct integration:

Two-loop
At the two-loop order we start with the set of MIs where the T i are shown in figure 6. The MIs F admit -linear differential equations, except for one of them. We have indeed where K 0 , K 1 and K 2 do not depend on , and K 2 is non-vanishing only in the inhomogeneous part of the differential equation for F 36 . In a first step we apply the Magnus algorithm on K 0 + K 1 in order to remove K 0 , and in a second step we apply an ad-hoc transformation in order to remove the remaining non-linear piece. The corresponding canonical basis reads , As compared to the one-loop case (5.3) we encounter one additional square root in the canonical d log-form which is not rationalized by the change of variables in eq. (5.4). In terms of w and z, the alphabet reads the argument of the square root entering η 13,...,17 is 15) and the four coefficients in η 16 are given by In the region 0 < w < z < 1 all the letters η i are positive. As already stressed, the alphabet is not rational in w and z. This prevents us from expressing the complete solution in terms of GPLs. In particular, the structure of the coefficient matrices M i is such that the solution for I  .13), (3.14), involves path integration over d log's with non rational arguments. Nevertheless, the MIs I 1,...,31 admit a representation in terms of GPLs which is convenient for their numerical evaluation. As for the remaining MIs, we followed the procedure outlined in section 3.3: we express the solution up to weight 2 for I 32 and up to weight 3 for I 33,...,36 in terms of GPLs and then obtain an 1-fold integral representation for the higher weights (for I (4) 32,...,36 we use eq. (3.29)). We hereby list the conditions imposed to integrals I 1,...,31 for determining their boundary constants; • independent input: I 1,4,..., 7,13,14,20,25 , • regularity at s → 0: I 24 , • regularity at t → −m 2 : I 12,21 , • regularity at u → 0: I 19,30,31 , • limit s → 0: I 2,3,6,...,10,15...18,29 , • limit t → −m 2 and s → 0: I 28 , • regularity at s → 0 and matching to independent input: I 22,23 .
For the MIs I 32,...,36 we observe that regularity at u = s = t = 0, corresponding to x 0 = (w 0 , z 0 ) = (1, −1), implies that we choose as initial condition of our solution in terms of iterated integrals.
The MIs I 1,...,31 are represented in terms of GPLs, and can be computed on the whole (s, t) plane (except for the line s = 4m 2 , see appendix A for further comments).
The explicit evaluation of I 32,...,36 requires a careful choice of the integration path, in such a way that no branch cuts are crossed. We successfully checked our results in the unphysical region s < 0 (see appendix A) against the numerical values obtained with SecDec. The evaluation of our analytic result relies on the use of GiNaC for the computation of the GPLs and on a one-dimensional integration for the cases where non-rational weights appear in the most external iteration, according to the eq. (3.29). As for the latter, we exploited the propriety of path-independence to choose simple paths (that avoid the singularities on the way from the basepoint to the chosen endpoints). Let us remark that in this work we did not focus on the the computing performances of the numerical evaluation of the mixed Chen-Goncharov iterated integrals appearing in our analytic expression. This aspect, together with a study of the analytic properties of our solutions in the whole phase-space, requires a dedicated future investigation.
The analytic expressions of all the MIs are explicitly given in electronic form in ancillary files that can be obtained from the arXiv version of this paper.

Conclusions
In this article, we presented the calculation of the master integrals (MIs) needed for the virtual QCD×EW two-loop corrections to the Drell-Yan scattering processes, for massless external particles. Besides the exchange of massless gauge bosons, such as gluons and photons, the relevant Feynman diagrams involve also the presence of W and Z propagators. Given the small difference between the masses of the W and Z bosons, in the diagrams containing both virtual particles at the same time, we performed a series expansion in the difference of the squared masses. Owing to this approximation, we distinguished three types of diagrams, according to the presence of massive internal lines: the no-mass type, the one-mass type, and the two-mass type, where all massive propagators, when occurring, contain the same mass value. The evaluations of the four point functions with one and two internal massive propagators are the main novel results of this communication.
To achieve it, we identified a basis of 49 MIs and evaluated them with the method of the differential equations. With the help of the Magnus exponential, the MIs were found to obey a canonical system of differential equations. Boundary conditions were imposed either by matching the solutions onto simpler integrals in special kinematic configurations, or by requiring the regularity of the solution at pseudo-thresholds. The canonical MIs were given as Taylor series around d = 4 space-time dimensions, up to order four, whose coefficients were given in terms of iterated integrals up to weight four. The solution could be expressed in terms of Chen's iterated integrals, yet, we adopted a mixed representation in terms of Chen-Goncharov iterated integrals, suitable for their numerical evaluation. Further studies concerning the analytic properties of the presented MIs in the whole phase-space, and the optimization of their numerical evaluation will be the subject of a forthcoming publication.

Acknowledgments
We would like to thank Valery Yundin for his contribution during the early stages of this project. We thank Lorenzo Tancredi for clarifying discussions and Matthias Kerner for technical support on SecDec. Some of the algebraic manipulations required in this work were carried out with FORM [101]. Some of the Feynman diagrams were generated by FeynArts [102,103] and drawn with Axodraw

A Variables for the one-mass and two-mass integrals
In this section we discuss the domain of the variables employed in the analytic expressions of the MIs for the Drell-Yan process, both in the case with one massive propagator and in the one with two massive propagators.

A.1 One-mass type
For the evaluation of the one-mass MIs we simply rescale by the squared mass the Mandelstam invariants. All the analytic results are given in terms of two-dimensional generalized polylogarithms, functions of the variables In the unphysical region s < 0, x is real and positive. Correspondingly, y can be either positive or negative. The analytic continuation to the physical region requires the Feynman prescription on the invariants. There, s becomes positive, with a positive vanishing imaginary part, s + i0 + . Accordingly, x is negative: On the other hand, t is negative (with a positive vanishing imaginary part) and ranges between 0 and −s, −s < t < 0. The numeric evaluation of the MIs expressed in terms of GPLs of the variables x and y can be done in the whole s, t plain using the routines in [79] expressing our analytic formulas in terms of GPLs evaluated in 1 and giving the explicit imaginary part to the Mandelstam variables (see for instance [105]).

A.2 Two-mass type
For the evaluation of the two-mass MIs, see section 5, we find it convenient to introduce the reduced variables w and z defined by We note that the above mapping allows the evaluation of our results everywhere in the (s, t) plane, with the exception of the value w = −1 (corresponding to s = 4m 2 ). For that specific value of w, the t dependence in z gets lost by construction, and z = −1 independently on t.
The evaluation of the solution at s = 4m 2 requires further investigations and it will be addressed in a forthcoming publication.
where we explicitly used the Feynman prescription s + i0 + .

A.2.2 Range of values for z
The variable z depends both on s and t. In order to study the different regimes, we define the following function of s where the second equality follows from eq. (A.5). We also define the ratio so that the second of eqs. (A.4) reads We choose the following root of the above equation Note that eq. (A.12) contains square-roots of K. Therefore, in order to compute z when K < 4, we have to keep track of the vanishing imaginary parts of the quantities entering eq. (A.10). Region by region in the (s, t) plane, the correct sign of the vanishing imaginary part (if present) is determined by the Feynman prescription on s, t, u, i.e. s + i0 + when s > 0, and likewise for t and u.
Depending on the value of K, we distinguish three cases (here we keep the prescription for the vanishing imaginary part of K arbitrary): All the square roots in eq. (A.12) are real, so z is real with 0 < z < 1.
Note that, since K is a function of s and t, each case can arise from multiple regions in the (s, t) plane. In table 1 we summarize the solution for z in the different regions of the (s, t) plane, by displaying also the appropriate sign for the i0 + prescription (if a vanishing imaginary part is present).

B Two-Loop d log-forms
In this appendix we give explicitly the coefficient matrices of the d log-forms, eq. (3.4), for the one-mass and the two-mass two-loop MIs, discussed respectively in sections 4 and 5.

B.1 One-mass
For the one-mass case at the two-loop order, the d log-form is