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 QCD × EW corrections to q+q¯→l−+l+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q+\overline{q}\to {l}^{-}+{l}^{+} $$\end{document} and q+q¯′→l−+ν¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q+\overline{q}^{\prime}\to {l}^{-}+\overline{\nu} $$\end{document}, for massless external particles. We treat the 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 pseudothresholds. 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. Furthermore,

JHEP09(2016)091
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 the detector calibration, and for the 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 [2]. 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 the effort needed to increase the 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 [3,4]. 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 [5][6][7][8][9][10].
At the same order of accuracy (one can roughly thing to exchange two powers of α S with one power of α), the mixed EW-QCD corrections have to be taken into account. As in the case of QCD NNLO with EW NLO perturbative corrections, the mixed EW-QCD corrections are expected to become of similar size with respect to QCD NNNLO at high leptonic invariant mass [11].
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 [12,13] and revealed a sizable increase of the cross section with respect to the LO result. The NNLO QCD corrections [14,15] 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 [16][17][18][19].
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 [20] and Z [21] 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 [22][23][24][25][26][27] and in the full JHEP09(2016)091 SM [28]. NLO corrections to the production of a vector boson with a jet were considered in [29][30][31][32]. The two-loop QCD helicity amplitudes for the production of a Z or a W with a photon have also been calculated [33]. 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 [34][35][36][37][38][39][40][41][42][43][44]. Finally, the rapidity distribution of a vector boson is known at the NNLO in QCD [45].
Pure QED generators are also available [59][60][61][62][63][64][65][66][67][68][69]. 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 MC programs FEWZ [70,71], that includes also EW NLO corrections [72], and DYNNLO [73,74]. In these two programs, the decay products of the vector boson, the spin correlations and the finite-width effects are also taken into account.
A sizable impact on the pp(p) → W → lν distributions, and therefore on the determination of the W mass, comes from the QCD initial state radiation (ISR) with QED final state radiation (FSR) or from the real-virtual (factorizable) corrections. However, at the level of precision required (∆m W ∼ 10 MeV), the complete set of mixed EW-QCD corrections may be important and has to be considered.
The NNLO mixed EW-QCD corrections to the production of a leptonic pair, i.e. order αα S corrections to the LO partonic amplitude, consist of two-loop 2 → 2 processes, in which the quark-antiquark initial state goes in the final leptonic pair (l + l − or lν), one-loop 2 → 3 processes, in which the final leptonic pair is produced together with an unresolved photon or gluon, and tree-level 2 → 4 processes in which the leptonic pair is produced together with an unresolved photon and an unresolved gluon.
The QCD×QED perturbative corrections were considered in [78]. In [79], the mixed two-loop corrections to the form factors for the production of a Z boson were calculated analytically, expressing the result in terms of harmonic polylogarithms and related generalizations. In [80], the authors calculated the mixed corrections in the pole approximation near the resonance region. In particular, they worked out in this region three different contributions. The first are the ones coming from the QCD corrections to the production process, which turned out to be suppressed below the percent level. Then, there is the photonic final state radiation, which is the dominant contribution. Finally, they considered the non-factorizable contributions, due to soft-photon exchange between the production and the decay processes. These last contributions are suppressed way below the percent level and they are negligible for the curent phenomenological purposes. In [81], the mixed αα S JHEP09(2016)091 contributions in the resonant region coming from factorizable initial-final corrections were fully exploited. These contributions represent the bulk of the corrections near resonance. The mixed αα S corrections worked out in the pole approximation are phenomenologically sufficient for a detailed description of the Drell-Yan processes near the resonance. However, in order both to further check this statement and to be able to treat the process also in physical regions different from the resonance, a complete diagrammatic calculation would be desirable.
In this article, we present the calculation of the master integrals (MIs) needed for the virtual corrections to the two-loop 2 → 2 processes: q +q → l − + l + , and q +q → l − + ν , for massless external particles. The Feynman diagrams contributing to the process q +q → l − + l + can involve the exchange of up to two virtual massive vector bosons. However, they never contain W and Z propagators at the same time. As a consequence of that, the master integrals involved in the calculation contain at most two equal internal masses and therefore they depend three scales.
The situation is different for the diagrams contributing to the process q +q → l − + ν. In this case, double massive exchanges can occur with the two different bosons in the same Feynman diagram, and therefore one should evaluate master integrals depending on four different scales. However, 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 reducing effectively the number of scales on which the corresponding masters depend to three. Within this approximation, all topologies appearing in the two-loop QCD×EW virtual corrections to Drell-Yan scattering shall contain either no internal massive line, or one massive propagator, or two massive propagators with the same mass [82]. 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 [83,84], 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 [85][86][87], reviewed in [88,89]. 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 [90], following the algorithm based on the use of the Magnus exponential, introduced in [91,92]. 1 Boundary conditions are retrieved either from the knowledge of simpler integrals emerging in specific kinematic limits, 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 , d being 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 in general be expressed in terms of Chen's iterated integrals. We

JHEP09(2016)091
actually adopt a mixed representation, where, when possible, we make explicit the presence of Goncharov polylogarithms (GPLs) [96,97], 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 [87,[98][99][100], the four point integrals with one and two massive internal lines considered here are new and represent the main result of this communication.
All the master integrals with the exception of 5 of them were cast in closed analytic form in terms of GPLs. Because of the presence of irreducible irrational weight functions, it has been necessary to cast 5 of the 17 MIs with two massive internal lines as one-dimensional integral formulas [101], involving GPLs in the integrands. The numerical evaluation of our solutions can, therefore, be performed with the help of the GiNaC library [102] 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 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, and we present the corresponding solutions. Conclusions are given in section 6. In appendix A, we discuss the kinematic domain of our analytic results. In appendix B, we give the coefficient matrices of the systems of differential equations in canonical form.
Our results are collected in ancillary files, available in the arXiv submission.

Notation and conventions
In this paper we study the two-loop corrections to the following partonic scattering processes: The external particles are considered massless and 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, due to momentum conservation, we have s + t + u = 0. The physical region is defined by 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. We calculate the quantum corrections to the processes in eqs. (2.1) and (2.2) using a Feynman-diagrammatic approach. After considering the interference with the leading order, and summing over the spins and colors, we express the squared absolute value of JHEP09(2016)091 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 [103,104] and Lorentzinvariance identities [87], implemented in the computer program 2 Reduze 2 [83,84]. 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 of the exchange of a virtual gluon between the initial-state quarks. The final state is not affected, and at most massless three-point functions have to be evaluated. The EW corrections, instead, consist of the exchange of a photon, a Z boson or a W boson. Moreover, these quanta can be exchanged between the quarks in the initial state as well as between 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 perform a series expansion in the difference of the two squared masses. Expanding for instance the Z propagators around m W , we find: is the effective parameter of the expansion. The coefficients of the series in ξ are Feynman diagrams with equal masses, that therefore depend only on s, t, and the W mass, m W . The expanded denominator will eventually appear with a power α > 1. 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. At the one-loop level, the topologies involved in the QCD and EW corrections are shown in figure 1, where we distinguish: a) the massless case; b) the exchange of one massive particle; and c) the exchange of two massive particles.

JHEP09(2016)091
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: 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) and (2.11) we used = (4 − d)/2.

Systems 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 differential equation,

JHEP09(2016)091
where the matrix K depends both on the kinematic variables and on the spacetime dimension. By means of a suitable basis transformation, built with the help of the Magnus exponential [91,116] following the procedure outlined in section 2 of [92], we obtain a canonical set of MIs [90]. 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 ∂ a ≡ ∂/∂a and a and b are such that a, b = x, y (i.e. ∂ a is the derivative with respect to one of the kinematic variables).

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 P exp 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 [117] (see also [118] and the pedagogical lectures [119]) 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

JHEP09(2016)091
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 ) 10) 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 [117]. 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 C • 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 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 .

JHEP09(2016)091
• 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 JHEP09(2016)091 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:

Constant Goncharov polylogarithms
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 contains a combination 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 [120][121][122][123].

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 : where the T i are depicted in figure 3. By means of the Magnus exponential [91,116], according to the procedure outlined in section 2 of [92], 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. As a consequence, the analytic continuation to arbitrary complex values of x and y is straightforward, allowing the solution to be evaluated for any real and complex values of s, t, and m 2 . See appendix A for further details.
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: JHEP09(2016)091

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 JHEP09(2016)091 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 list the conditions imposed to the MIs for the determination of 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 . This leaves us with I 1,6 , to be determined by direct integration:

JHEP09(2016)091
. (4.10) As in the one-loop case, owing to the explicit representation in terms of GPLs, the analytic continuation to arbitrary x and y is straightforward, so that all the one-mass MIs can be computed for any real and complex values of s, t, and m 2 (see appendix A). Our results have been successfully checked against SecDec, both in the Euclidean and in the physical regions.
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 alphabet 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. 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). As already discussed in section 5 for the one-mass case, owing to the explicit representation in terms of GPLs, the analytic continuation to any w and z is straightforward, so that all the one-loop two-mass MIs can be computed for arbitrary real and complex values of s, t, and m 2 (except for s = 4m 2 ). For a detailed discussion of how the regions in the real s, t plane are mapped to the C × C space of the w, z variables, see appendix A. 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 ,

JHEP09(2016)091
As compared to the one-loop case (5.3) we encounter one additional square root in the canonical d log-form which prevents the alphabet from being 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 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  .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 we use eq. (3.29)).
For the MIs I 32,...,36 we observe that regularity at u = s = t = 0, corresponding to 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. As already discussed for the oneloop case, for such MIs the analytic continuation to arbitrary w and z is straightforward.

JHEP09(2016)091
Therefore the MIs I 1,...,31 can be computed for any real and complex values of s, t, and m 2 (except for s = 4m 2 , see appendix A for further comments). We checked our results in the Euclidean and physical regions against SecDec finding complete agreement.
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 Euclidean 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 efficiency 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 canonical systems 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 pseudothresholds. 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. While the solution could be expressed, in full generality, in terms of Chen's iterated integrals, 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 JHEP09(2016)091 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 with SecDec. Some of the algebraic manipulations required in this work were carried out with FORM [124]. The Feynman diagrams were generated by FeynArts [125,126] 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 case 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 Euclidean region s, t < 0. Therefore, if m 2 > 0, both x and y are real and positive.
The analytic continuation to the other values of s and t, with m 2 > 0, requires the Feynman prescription on the invariants. In particular, in the physical region s becomes positive, with a positive vanishing imaginary part, s + i0 + . Accordingly, x is negative, with a negative vanishing imaginary part: On the other hand, t is negative and ranges between −s and 0, so that 0 < y < x . The extension to the case of complex mass, m 2 → m 2 − imΓ, is straightforward. 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, m 2 ) domain using the routines in [102] expressing our analytic formulas in terms of GPLs evaluated in 1 and giving the explicit imaginary part to the Mandelstam variables (see for instance [128]).

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 of t. The evaluation of the solution at s = 4m 2 requires further investigations and it will be addressed in a forthcoming publication.
The arguments in the following sections rely on the assumption that m 2 > 0. Nevertheless, the whole discussion can be straightforwardly extended to the case of complex mass, m 2 → m 2 − imΓ. For instance, it is easy to see that in the physical region, defined by s > 0 and −s < t < 0 (see section 2), the real and imaginary parts of (w, z) would differ from those in the zero-width case. Nevertheless, the sign of the imaginary parts will be preserved. Therefore the zero-width case would simply arise as a limit. Furthermore, in the presence of a non-zero width, the change of variables (A.4) is always well defined for any real s, thus circumventing the difficulty mentioned above.

A.2.1 Range of values for w
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 and K − 4. 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

B.2 Two-mass
For the two-mass case, at the two-loop order, the d log-form is where we used the abbreviations introduced below eq. (5.13). The coefficient matrices are JHEP09(2016)091 (B.10)
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.