Multiple polylogarithms with algebraic arguments and the two-loop EW-QCD Drell-Yan master integrals

We consider Feynman integrals with algebraic leading singularities and total differentials in $\epsilon\,\mathrm{d}\ln$ form. We show for the first time that it is possible to evaluate integrals with singularities involving unrationalizable roots in terms of conventional multiple polylogarithms, by either parametric integration or matching the symbol. As our main application, we evaluate the two-loop master integrals relevant to the $\alpha \alpha_s$ corrections to Drell-Yan lepton pair production at hadron colliders. We optimize our functional basis to allow for fast and stable numerical evaluations in the physical region of phase space.


Introduction
The Drell-Yan process [1] is one of the most important and basic processes measured at the Large Hadron Collider at CERN. It is used for precision measurements of the W ± mass, the Z mass, and the weak mixing angle, as well as for new physics searches. In order to interpret the increasingly precise data, it is essential to include higher-order corrections in the theory predictions. In fixed-order perturbation theory, the pure Quantum Chromodynamic (QCD) corrections to the cross section are exactly known through to next-to-next-to-leading order [2][3][4][5], whereas the electroweak (EW) corrections are exactly known at just next-to-leading order [6][7][8]. Since these corrections turn out to have an important impact on the cross section, higher-order corrections should be incorporated. The next logical step is the inclusion of mixed EW-QCD corrections, which up to now have been studied in the pole approximation [9,10]. For the subset of mixed Quantum Electrodynamic (QED)-QCD corrections, the two-loop amplitude and its infrared structure have been calculated in [11].
The two-loop master integrals relevant to the virtual part of the mixed EW-QCD corrections to the Drell-Yan process have been studied in [12,13]. As it turns out, these integrals fall into a certain class, for which it is possible to derive -decoupled and, subsequently, d ln differential equations [14] after introducing various algebraic prefactors and forming appropriate linear combinations. For many of the integrals, no intrinsically nonrational prefactors appear in the -decoupled differential equations; these prefactors are either rational from the get-go or they become so in a particular parametrization of the kinematics. It has come to light, however, that the two most complicated two-loop integral topologies (right panels of Figures 1 and 2) involve three square roots at once which can not be rationalized simultaneously [15]. This leads to genuinely non-rational symbol letters [16] in the d ln differential equations, which makes the standard approach to the integration of d ln differential equations in terms of multiple polylogarithms impossible.
In the literature, other processes have been studied for which no simultaneous rationalization of the square roots which appear is forthcoming, e.g. a subset of the planar two-loop master integrals for Bhabba electron-positron scattering [17] and a subset of the planar two-loop master integrals for the next-to-leading order QCD corrections to Higgs plus jet production with full heavy top quark mass dependence [18]. So far, master integrals satisfying d ln differential equations which involve non-rational symbol letters have typically been treated as generic Chen iterated integrals [19], with a focus on the Euclidean region of phase space. Given our experience from other processes [20,21], one may expect significantly faster and also more stable numerical evaluations for results directly expressed in terms of suitably chosen multiple polylogarithms. A constructive algorithm to obtain such a solution in the presence of root-valued symbol letters was used in [22] to calculate Feynman integrals. However, the latter case is special since the symbol alphabet is univariate and simple, and a rational parametrization does in fact exist. In other words, it was unclear until now whether it should generically be possible to find linear combinations of multiple polylogarithms which solve d ln differential equations with symbol letters involving unrationalizable algebraic functions.
In this article, we show for the first time that it is indeed possible to integrate Feynman integrals of current phenomenological interest with unrationalizable roots in their symbol letters in terms of multiple polylogarithms, focusing primarily on the most complicated two-loop mixed EW-QCD Drell-Yan master integrals. Our method allows us to derive results written in terms of multiple polylogarithms well-suited for the physical region of phase space. The two-loop mixed EW-QCD Drell-Yan master integrals with two massive internal lines exhibit a rich structure of thresholds and pseudo-thresholds, which means that one must think carefully about the analytic structure of the involved functions. We discuss in detail a procedure to systematically filter the functions in a given region of phase space, such that no explicit Feynman +i 0 prescriptions are required for the kinematic variables and spurious singularities of individual multiple polylogarithms at pseudo-thresholds of the Feynman integrals are avoided. The outline of the article is as follows. In Section 2 we discuss the linear reducibility of examples which are well-known in the literature to frustrate the standard machinery of Feynman integral calculus. Specifically, we treat the integral topologies of Figure 2: a five-line master integral for Bhabba scattering and a six-line master integral for Drell-Yan production. Linear reducibility is a technical criterion which is of interest here because linearly reducible Feynman integrals are guaranteed to be integrable in terms of multiple polylogarithms to all orders in . In Section 3, we define a normal form basis with d ln differential equations for the two-loop master integrals for the mixed EW-QCD corrections to Drell-Yan production with two massive internal lines (see the right panel of Figure 1) and discuss a partial rationalization of the roots appearing in our integral basis definition. In Section 4, we show how to integrate the differential equations directly in terms of multiple polylogarithms even in the presence of root-valued symbol letters. In Section 5, we review the analytic continuation of multiple polylogarithms and outline our procedure to filter multiple polylogarithms with undesirable analytic properties out of our ansätze. In Section 6, we present results for the most complicated two-loop mixed EW-QCD Drell-Yan master integrals. In particular, we highlight the notable analytic features of our solution for the six-line integral from the right panel of Figure 2. We conclude in Section 7 and, for clarity, we give the set of complete d ln differential equations for the two-loop master integrals considered in this paper in Appendix A.

Linear reducibility for algebraic symbol letters
In this section, we discuss the direct integration of a five-line master integral for the two-loop QED corrections to massive Bhabha scattering (left panel of Figure 2). This particular five-line integral is of special interest because its symbol letters are believed to not be simultaneously rationalizable [23]. It has been known for quite some time that the planar master integrals for the two-loop QED corrections to Bhabha scattering (see the left panel of Figure 1) satisfy -decoupled differential equations [17], but, even at leading order in the expansion, it is not at all clear that integral (2.1) below may be expressed as a linear combination of standard multiple polylogarithms [24][25][26]. In fact, it was suggested in [27] that elliptic multiple polylogarithms [28] might actually be required for such cases. It is therefore of some importance to demonstrate that one can actually integrate it to all orders in in terms of Goncharov or multiple polylogarithms using a HyperInt-like approach [29][30][31][32][33][34][35][36][37][38][39].
In fact, this may be readily achieved by studying the polynomials which appear in the polynomial reduction [29,30,33] of the Symanzik polynomials at intermediate stages. The idea is to make changes of variables which factorize totally quadratic polynomials to allow for further integrations without producing square root-valued functions of the remaining Feynman parameters. Ultimately, the goal is to find a complete sequence of integration variables or, in the language introduced in [30] and refined in [31], to show that (2.1) is a linearly reducible Feynman integral in the chosen variables. Practically speaking, a Feynman integral is linearly reducible if it can be evaluated in terms of Goncharov polylogarithms, starting from the Feynman parametric representation but allowing for arbitrary variable changes along the way. In this case, the Feynman parametric representation is are the integral's Symanzik polynomials. We implement the delta function constraint at the very beginning by setting α 5 = 1 inside U and F. Our normalization for the integral is chosen to facilitate comparison with FIESTA 4 [40]. At the outset, the polynomial reduction of U and F indicates that just one integration can safely be performed, with respect to either α 1 or α 3 . The irreducible polynomials in the remaining Feynman parameters which could potentially appear in the result after integrating out α 1 are [31]: The final totally quadratic polynomial in L α 1 is conveniently dealt with by using the first non-trivial variable change from the discussion of the period P 7,11 in [33]. It is readily apparent from the form of the polynomial that two powers of α 2 factorize from it once the variable change is applied: Quite remarkably, the other polynomials in L α 1 which depend on α 4 still give rise exclusively to irreducible polynomials linear in α 2 after changing variables: From Eqs. (2.6), (2.7), and (2.8), it is clear that we can now safely integrate out α 2 .
Rerunning the polynomial reduction with HyperInt after changing variables shows that the irreducible polynomials could appear in the result after integrating out α 1 followed by α 2 . L α 1 α 2 is more nontrivial to treat than L α 1 and requires us to think carefully, both about what sort of result we expect to find and how the Mathematica-based private direct integration code written by one of us actually operates in practice. As should by now be clear, our goal is to delay the appearance of square roots until the final variable of integration is moved into the arguments of the Goncharov polylogarithms which depend on it. With the final integration of the most complicated weight-three functions in mind, it is useful to aim for a final integration domain of [0, 1]. This will naturally produce Goncharov polylogarithms of argument 1 which contain a complicated square root in the weights. Actually, due to the way that our integration script is written, it is most natural for us to integrate complicated finite Feynman parameter integrals such as the one considered in this section over the unit hypercube. This is because it is easier for our integration code to take a definite integral on [0, 1] than on [0, ∞), due to the fact that no complicated argument inversion formulae for non-trivial Goncharov polylogarithms need to be derived if the upper integration endpoint is always set to 1.
In this spirit, we initiate our treatment of L α 1 α 2 by trivially mapping one of the remaining integration domains onto [0, 1] via the transformation It is also useful to note that two of the quadratic polynomials in Eq. (2.9) depend on just a single Feynman parameter each, a situation which strongly suggests that a change of kinematic variables would be advantageous. Indeed, with Euclidean s and t in mind, we see that the variable changes rationalize two square roots which would otherwise arise from the third irreducible polynomial of L α 1 α 2 , 12) and the fourth irreducible polynomial of L α 1 α 2 , Curiously, transformation (2.11) was also found to be useful in the context of the multi-loop QED corrections to light-by-light scattering [41]. At this stage, it turns out that just a single non-trivial algebraic function of x 4 is generated if one attempts to naively continue the calculation by integrating out x 3 : (2.14) To our knowledge, a situation such as this was first discussed in the Feynman integral literature in the context of the fully-massless (three-loop) K 4 integral [42]. Using the method of parametrization by lines [42,43], it is completely straightforward to derive a final variable change for x 4 which rationalizes (2.14) and produces an integration domain of [0, 1]. As usual, one begins by identifying a suitable rational point of the algebraic variety For our purposes, the point is ideal. Eq. (2.15) defines a hyperbola which can be rationally parametrized by a family of lines passing through (2.16), We combine Eqs. (2.15) and (2.17) to determine x 4 as a function of y 4 : .
As desired, we see from the above that the integration domain for y 4 is the unit interval. By combining together the various variable changes given above, we obtain the complete sequence of integration variables {α 1 , α 2 , x 3 , y 4 } and establish that, up to an overall normalization factor, integral (2.1) may be evaluated as a rational linear combination of Goncharov polylogarithms to all orders in . As a sanity check, we explicitly evaluate all Feynman parameter integrals analytically using our direct integration code at O 0 and then evaluate the result obtained numerically to high precision using GiNaC [44,45] at the Euclidean phase space point v 1 , v 2 , m 2 = (1/7, 1/5, 2) . (2.19) We find which we could rapidly confirm agrees to five significant digits with an independent FIESTA 4 evaluation of the integral. Let us stress once more that the allegedly unrationalizable square root where GiNaC. If desired, a representation written strictly in terms of Goncharov polylogarithms may also be derived by fully factoring all of the generalized weights. Of course, the direct integration approach described in this section is not without its limitations. For one thing, it would at first seem quite non-trivial to find an explicit integration order for e.g. the two-loop double box with two massive internal lines from the two-loop mixed EW-QCD corrections to Drell-Yan lepton production (right panel of Figure  1). Actually, with the help of HyperInt, it is comparatively easy to find a complete sequence of integration variables; for a fixed number of loops and legs, one gets the impression that the number of internal lines is a less important factor in determining the difficulty of a linear reducibility analysis than the total number of massive lines. Although we could directly integrate through to weight four the six-line basic scalar integral (right panel of Figure 2) which first inherits unrationalizable symbol letters from its unrationalizable leading singularity, we find the direct integration of the top-level two-massive-line integrals to be very cumbersome at the technical level. Furthermore, the produced Goncharov polylgarithms are very complicated, not minimal, and suboptimal for numerical evaluations with GiNaC.
An alternative approach which can avoid these issues is to apply the differential equation method. We shall see that the differential equation method has the very appealing feature that one can fit an ansatz of multiple polylogarithmic functions appropriate for each physical kinematic region separately. In the following sections, we discuss in some detail how this is achieved for the mixed EW-QCD Drell-Yan master integrals with two massive lines, including, in particular, the six-and seven-line integrals with an unrationalizable square root in their symbol letters.

An -basis for the Drell-Yan master integrals with two massive lines
In this work, we are primarily concerned with neutral-current lepton pair production 1 in quark-antiquark annihilation, where all external particles are taken massless. The most complicated master integrals for the two-loop mixed EW-QCD corrections to above process are those with two internal lines of mass m, where, depending on the parent Feynman diagram, m may refer to either m W or m Z . In particular, as mentioned above, it has been known to us for some time that the integral topologies from the right panels of Figures 1 and 2 actually contain master integrals with unrationalizable symbol letters in their -decoupled differential equations. To index the seventeen master integrals with two massive internal lines, it suffices to consider a single integral family based on the two-loop planar double box with two massive internal lines from Figure 1 (Family E from Table 1). For completeness, we will also give definitions for those normal form integrals with zero or one massive internal lines which appear on the right-hand sides of the -decoupled differential equations for the master integrals with two massive internal lines. To index these auxiliary integrals, we reintroduce two additional integral families which were already studied in the physical region by two of us some time ago [13] (Family A and Family C from Table 1). In order to obtain a closed system of differential equations for the masters with two massive internal lines, we need to consider thirty-six integrals in total. Our notation for Feynman integrals in this section is exactly that of [13] (i.e. dots for doubled propagators, heavy lines for massive propagators, numerator insertions written in square brackets, and F:x for the crossed version of sector x from family F). For the kinematic invariants we use In the following, we keep the dependence on the internal mass parameter m implicit for the sake of brevity and since it is anyway clear from the thick-line notation. We build up 1 The two-loop mixed EW-QCD corrections to the charged-current Drell-Yan process could be accessed by expanding in 1 − m 2 W /m 2 Z , as this would allow one to make effective use of the equal-mass integrals relevant to the neutral-current process. Table 1. Integral families for the master integrals which appear in the differential equations for the two-loop mixed EW-QCD Drell-Yan master integrals with two massive internal lines.
our normal form basis out of the following thirty-six Feynman integrals: From the general principles discussed in [37], one can readily cast (3.3) above into a normal form basis for the integrals of interest. 2 Abbreviating the three square roots which appear as we find: Due to the presence of unrationalizable square roots, the available public software packages for the construction of a normal form basis of integrals [46][47][48] are not applicable to the problem at hand. For m 1 , . . . , m 36 , we employ the integration measure , which allows us to consider our integrals to be functions of two dimensionless kinematic variables. 3 For instance, One achieves a rationalization of two of our three roots, r 1 and r 2 , with the parametrization [12] and we hereafter work primarily with the variables w and z. It is immediate from (3.7) that, in the (w, z) representation, the only normal form integral which involves a square root in its definition is m 32 (i.e. the integral from the right panel of Figure 2). Indeed, we see from (3.10) and Appendix A of [12] that in the region s > 4 m 2 (see also Section 5.1 for a detailed discussion of the (w, z) parametrization).
In the following, we will replace r 1 , r 2 and r 3 according to (3.11)-(3.13) in the normal form definitions m 1 , . . . , m 36 and use these partially rationalized expressions for the entire physical region of phase space. That is, only the definition of m 32 involves a root in the prefactor, which requires a non-trivial analytic continuation from s > 4 m 2 to other regions of phase space on its own. Using these definitions, it is straight-forward to obtain the differential equations in the -decoupled form discussed in the next section. We used Reduze 2 [45,[49][50][51] to compute the integration-by-parts identities required to derive the differential equations.

Integrating root-valued symbols in terms of multiple polylogarithms
Our normal form basis for the mixed EW-QCD Drell-Yan integrals with two massive internal lines, (3.7), is chosen to bring the associated differential equations into an -decoupled form [14,52]: where l k are the symbol letters and A (k) are matrices of rational numbers. It was demonstrated already in [12] that a d ln form does in fact exist, but their choice of the symbol letters is not optimal for our purposes. For now, we proceed with the understanding that some algebraic symbol letters appear in Eq. (4.1), but leave their precise form to be determined by the analysis below. For x = (w, z) and m = (m i ), one can give a formal solution of Eq. (4.1) in terms of Chen iterated integrals: m( , x 0 ) are the boundary constants of the master integrals at the point x = x 0 , γ is a piece-wise smooth path connecting x 0 and x, and the path-ordered exponential in Eq. (4.2) is defined in the usual way as an infinite series of integral operators acting to the right, In (4.4), the product of two or more d ln terms is understood as an instruction to take the corresponding iterated integral of the kernel along the path γ. As a concrete example, let us consider a straight-line path τ on the real axis from 0 to x. One can write, e.g.: Note that, in the above example, one could equally well identify the iterated integral as a Goncharov polylogarithm, If all symbol letters are linear, it is always possible to integrate Chen iterated integrals in closed form using G functions. This is no longer true, however, if one encounters nonlinear or non-rational letters. In the former case, it is often possible to choose an appropriate integration order for which the non-linear letters appear only in the final integration kernel. At this stage, one may employ the generalized weights used in Section 2 to obtain concise results. If one encounters non-rational symbol letters, the standard integration algorithms for G functions cannot be applied. In many cases, a transformation can be found which simultaneously rationalizes all letters of the alphabet. Nevertheless, it can be proven that the symbol alphabet of the five most complicated integrals from (3.7), {m 32 , . . . , m 36 }, cannot be rationalized [15]. Of course, it is a priori not obvious whether the functional basis for {m 32 , . . . , m 36 } consists solely of multiple polylogarithms; at the outset, it is certainly possible that one could have to deal with a more involved space of functions.
Fortunately, our linear reducibility analysis from Section 2 guarantees that standard multiple polylogarithms suffice; what we need is a better way to integrate d ln differential equations. A clear alternative is to proceed by matching the symbol of the Chen iterated integrals to a suitable ansatz built out of logarithms and Li functions. For rational alphabets, an explicit algorithm was provided in [53] by Duhr, Gangl, and Rhodes to construct suitable Li function arguments, such that the functional basis contains no spurious letters. At weight one, it is clear that logarithms of the letters are admissable functions. The non-trivial step is to find a sufficiently large set of admissible Li n arguments; once these arguments have been found, it is straightforward to construct all admissible arguments for the Li functions of depth greater than one. By considering the symbol of these functions, we see that it is desirable to admit only those function arguments f which have the property that both f and 1 − f may also be written as a power product of the symbol letters. In practice, one therefore forms power products f out of the letters and tests if 1−f factorizes over the alphabet. The symbols of higher-depth Li functions are more complicated and lead to additional constraints. To treat Li 2,1 , Li 3,1 , and Li 2,2 , let F be the union of the set of admissible Li n function arguments and the set {1}. Then, the symbol dictates that a possible pair of arguments for Li n 1 , factorizes over the alphabet. It has proven useful for practical applications to impose further constraints on the functional basis to ensure real-valuedness and good numerical performance. An implementation of this method written by one of us in Mathematica has been applied successfully to various processes [13,21,34,35,37,39].
In the presence of square roots, we use a heuristic factorization algorithm to detect admissable function arguments. For a given expression g we are interested in factorizations of the form with a rational number c and a n ∈ Q. It is non-trivial to find such factorizations using standard computer algebra systems due to the presence of the root r in the symbol letters. We observe that the factorization (4.8) implies Replacing the variables by numerical samples allows us to find these relations using heuristic integer relation finders. To find the required factorizations, we employ the Lenstra-Lenstra-Lovász (LLL) algorithm [54] implemented in PARI/GP [55] for a parallelized C++ code written by one of us. We would like to stress that the definition of the symbol letters is not unique. One can replace a letter by power products of letters and it is a priori unclear which choice is optimal for practical purposes. For example, we find that out of the 17 letters presented in [12], only 16 combinations are actually required for the integration of the Drell-Yan integrals. Furthermore, since we consider actual derivatives of Feynman integrals, we are not sensitive to numerical letters like 2. In principle, one needs to include also ad-hoc letters like −1, 2, etc. in the construction of function arguments f . This is a problem occurring also in the purely rational case and does not seem to be a major obstacle in practice. In the applications we have considered so far, including −1 and 2 was sufficient. For the case of the Drell-Yan integrals, we were able to absorb the letter 2 by a redefinition of our symbol alphabet, as will be explained below.
In practice, we encounter two main problems specific to the case of non-rational symbol alphabets: (i) In general, one needs to allow for non-integer powers, e.g. 1/2, 1/4, etc. when forming power products f . Consequently, one may have to test many more expressions than in the rational case in order to construct enough function arguments to successfully "integrate the symbol." In fact, without additional constraints, the inherent combinatorial complexity makes this extension of the Duhr-Gangl-Rhodes algorithm too costly for two-loop calculations of current phenomenological interest.
(ii) Factorization over algebraic functions is not unambiguously defined; in principle, it is not clear where to "stop" factorizing. Consider, for instance, the set of letters { √ x, √ y, x + y} and note that it is a priori unclear whether one ought to factor the third letter further, i.e. x + y = ( In the case of the Drell-Yan integrals, both of these problems can be tackled by the following observations: given an alphabet, we introduce the subset of rational letters, L R , and the subset of intrinsically non-rational letters, L A . Let us assume that in L A we encounter a square root, r. We find it natural to take r itself to be an element of L A . For a given algebraic letter with a non-trivial rational part, l a , we define the conjugated letter,l a , by making the replacement r → −r in l a . For each l a ∈ L A we observe thatl a can be written as a power product of other symbol letters; that is to say, one could exchange any letter for its conjugate without affecting the singularity structure of the alphabet. Furthermore, we observe that the product l ala always factorizes over the rational part of the alphabet, L R . We find these observations to be quite compelling and therefore conjecture that these structural constraints hold also for other symbol alphabets of this type. Our conjecture is a strong one, since it implies that one can predict the form of the remaining unrationalizable symbol letters which appear in the presence of the square root r: knowing only the rational part of the alphabet, L R , and the unrationalizable square root which appears, it allows one to construct the algebraic part of the alphabet, L A , unambiguously. As we shall see, this insight allows for a drastic simplification of the algebraic part of the symbol alphabet relative to what was presented in Eqs. (5.13)-(5.19) of [12] and what we were able to derive ourselves initially using the above-mentioned heuristic factorization code.
We want to show this in some detail for the EW-QCD Drell-Yan master integrals with two massive internal lines. At the outset, after carefully considering various possibilities, we find the rational alphabet and the intrinsically algebraic alphabet where r = (1 + w 2 z 2 )(w + z) 2 + 2w z(w − z) 2 + 4w z 2 (1 + w 2 ), using our heuristic factorization code. Note that r already appears as a letter in Eq. (4.11). However, there are also two rather complicated-looking symbol letters with terms involving r 2 .
As stated above, knowing only r and L R , we can also construct an improved representation of the algebraic part of the alphabet by making an ansatz of the form where q a is a rational function in w and z, and then requiring that l ala factorize over L R . This produces a much simpler representation of the algebraic part of the alphabet, which reads:L A = r, The overall factors of 1/2 in Eq. (4.13) are judiciously chosen after the fact to prevent the appearance of explicit factors of 2 in our Li function arguments; as alluded to above, our original construction of the functional basis appended 2 as an auxiliary letter.
Using our heuristic approach to algebraic function factorization, we immediately find that the non-trivial elements of L A can indeed be expressed as power products of letters drawn from the improved alphabet: In summary, we define the full symbol alphabet for the two-massive-line Drell-Yan integrals, L = L R ∪L A , in terms of positive definite letters for physical phase space points which satisfy −1 < w < 0 and w < z < −w 2 (this component of (w, z) space corresponds to part of the s > 4 m 2 region, see Section 5.1 for details). From Eqs. (4.10) and (4.13), we have: Most importantly, our improved representation of the alphabet effectively solves the problem of finding power products of high degree, since we observe in practice that we do not need to consider square roots of any letter in L. For the integration of the EW-QCD Drell-Yan master integrals with two massive internal lines through to weight four, it was enough to consider power products of symbol letters up to total degree 9.

Analytic continuation and optimization of the functional bases
In this section, we review the salient features of the (w, z) representation for the mixed EW-QCD corrections to Drell-Yan production introduced in Section 3 above, as well as subtleties one encounters when analytically continuing multiple polylogarithms. Our primary goal in Section 5.1 is to motivate the analysis of Section 5.2, where we show how we avoid all explicit analytic continuations and +i 0 prescriptions by partitioning the physical phase space and finding several solutions to the differential equations in terms of well-behaved Li functions valid in judiciously chosen regions. It does seem to us, however, that a lucid and detailed discussion of the relevant issues is not always provided in the Feynman integral literature. We therefore give a rather detailed discussion of the analytic continuation of Feynman integrals in Section 5.1 below and hope that it may be of interest to some readers.

Analytic continuation
Feynman's +i 0 prescription for the propagators determines the value of a given Feynman integral in a specific region of phase space unambiguously. In principle, one could imagine to solve the Feynman integral in each region of phase space separately. Alternatively, one can try to solve the integral in one region and then use the solution to obtain a result for it in a neighboring region by analytic continuation. The latter method is of particular interest for the method of differential equations, since it typically involves regularity conditions in some regions of phase space, and one needs to transport this knowledge to the region of interest.
If we want to continue a specific representation of the solution based on just the solution itself (without reference to the original Feynman integral), we need to make sure Feynman's +i 0 prescription is maintained by appropriate complex values of the kinematic parameters. It is is essential to observe that the analytic continuation is along a path in complex phase space and that the +i 0 prescriptions must be respected for all points along this path, not just for the start and endpoint. This is in general non-trivial and needs to be checked for the representation at hand.
What is commonly referred to as "analytic continuation" in the physics literature should really be regarded as a two-step procedure in general: (i) The actual analytic continuation in the mathematical sense: given a solution for one region, derive a solution for a connected region in some representation.
(ii) A possible change of functional representation in the new region such that no explicit +i 0 prescription is necessary.
We will now work out the details for our current application. From the second Symanzik polynomials of our Drell-Yan master integrals, we see that their Euclidean region is given by s < 0, t < 0 and m 2 > 0. 4 Considering other regions of phase space, we observe that the Feynman propagator prescription can effectively be implemented by the replacements s → s + i 0, t → t + i 0, since these are external scales, and m 2 → m 2 − i 0, since this is an internal scale. The +i 0 prescription is relevant for s (t) whenever s (t) is positive. For the discussion of w and z, we will see that it is sufficient to view m 2 as being normalized to 1 without any imaginary part.
Let us focus on the (w, z) parametrization of our integrals in the physical region of phase space, s > 0, t < 0 and m 2 > 0. It was noted in [12] that the (w, z) representation has a point of non-analyticity at the physical two-mass threshold s = 4 m 2 and a rather different character depending on whether s is above or below this threshold. We find it natural to adopt the definitions We can see from the above that w and z are real variables which satisfy −1 < w < 0 and −1 < z < −w 2 in the region s > 4 m 2 . With our choices (5.1) and (5.2), we find that the Feynman prescription implies w → w + i 0 and z → z + i 0 in this region. On the other hand, both w and z become pure phases in the region 0 < s < 4 m 2 . Therefore, it makes sense to explicitly extract their real and imaginary parts, to emphasize that, in the below-threshold region, the imaginary parts of w and z are fixed in terms of the real parts of w and z. In particular, we can deduce from Eqs. representation is given in Figure 3.
In the following, we will study how to analytically continue solutions between different regions in s. In practice, we work with simple straight-line paths in the complex (w, z) space, checking after the fact that the chosen paths of analytic continuation always preserve the +i 0 prescription for the original kinematic variables. In what follows, we carefully go through typical elementary examples of analytic continuation in order to clearly illustrate the subtleties for our Li functions which one must be wary of to avoid introducing errors.
First, let us illustrate the importance of taking into account the complete path of analytic continuation rather than just the start and endpoint. Consider the analytic con-  tinuation of ln w 2 along a straight-line path from w i = w (0) + i δ to w e = −w (0) + i δ, for real w (0) and δ such that 0 < w (0) < 1 and 0 < δ < w (0) (see Figure 4). Such a path allows one to connect solutions for s < 0 with solutions for s > 4 m 2 . Naively, one might be tempted to erroneously implement the analytic continuation as ln w 2 i = ln w 2 (0) = ln (−w (0) ) 2 = ln w 2 e and incorrectly conclude that the analytic continuation is trivial along the chosen path. The problem here is that the logarithm is a multi-valued function and one must therefore carefully check whether or not the specified path of analytic continuation forces the polynomial function argument to cross the branch cut of the logarithm on the negative real axis, (−∞, 0). In the absence of any branch cut crossings, knowledge of the endpoint of the path is sufficient. However, in the presence of one or more branch cut crossings, the function leaves its principal Riemann sheet and one must add an appropriate monodromy contribution, taking into account the details of the path of analytic continuation.
Let us spell out in detail how to correctly analytically continue solutions valid for different values of s. We parametrize our chosen path in w space as with a parameter v ∈ [0, 1]. We then have for the argument of ln w 2 . As depicted in Figure 5, our path in w space takes w 2 (v) from just above the point w 2 (0) − δ 2 on the positive real axis to just below the point w 2 (0) − δ 2 and, crucially, it passes through the negative real axis at w(1/2) = δ i. The monodromy contribution to ln w 2 in this case is well-known to be simply 2πi, due to the fact that our path induces just one counter-clockwise branch cut crossing. Therefore, the analytically continued function is ln w 2 + 2πi in a neighborhood of the endpoint of our chosen path. This conclusion may be quickly checked by rewriting the function ln w 2 as 2 ln(w) in the Euclidean region where 0 < [w] < 1 before carrying out the analytic continuation along our chosen path. When the function is viewed in this alternative way, no branch cut crossing is induced and one can simply rewrite the function to explicitly extract its imaginary part: for [w] > 0, 2 ln(w) = 2 ln(−w) + πi = 2 ln(−w) + 2πi = ln w 2 + 2πi. In the two-step picture mentioned above we thus have 2 ln(−w) + 2πi = ln w 2 + 2πi (5.8) Depending on the type of representation we consider, either the analytic continuation (i) or the rewriting of the function to be independent of +i 0 prescriptions (ii) is trivial in this example. The relevant monodromy contributions become more complicated for higher-weight Li functions. Already for Li 2 , a new feature emerges: moving across its branch cut on the positive real axis, (1, ∞), onto a Riemann sheet other than the principal one actually exposes the existence of a hidden branch point at 1. To see this, consider Euler's identity for the dilogarithm of 1 − w 2 , The representation furnished by the right-hand side of (5.9) and the above discussion of ln w 2 make it clear that Li 2 1−w 2 has non-trivial monodromy for the path w(v) defined above in Eq. (5.5). For the path w(v), we find: This implies that Li 2 1 − w 2 picks up a monodromy contribution of −2πi ln 1 − w 2 due to the branch point at w = 1, despite the fact the function is continuous at that point on the principal sheet. When considering general analytic continuations of Li n 1 ,...,n k functions, one must take into account all function arguments to obtain the monodromy contributions as one moves along the path of analytic continuation. Building on the work of Goncharov [57], it was shown in [58] how one can easily compute the monodromy of an arbitrary Li function in terms of monodromies of simple logarithms using the coproduct. In particular, the presence of a monodromy contribution can be detected by studying the first entry of the symbol. For our master integrals, the sheer number of distinct branch cut crossings which can arise from the arguments of the various Li functions renders the analytic continuation between different regions rather involved in practice. While our final goal is to obtain representations in terms of well-behaved Li functions for kinematic regions, we find it convenient to also employ auxiliary representations for the purpose of analytic continuations.
If available, a representation in terms of G functions with the kinematic variables in the argument avoids many of the subtleties involved in the continuation described above. For the integrals with unrationalizable alphabets, however, this is not an option. For such cases, we find it useful to use expansions around regular and singular points [59] for the continuation. Using high-precision numerical evaluations for two expansion points with an overlapping region of convergence and the PSLQ algorithm [60], one can effectively transport analytic integration constants.
What has been discussed so far allows us to construct a functional basis for a given region of phase space, integrate the symbol in terms of these functions, and relate solutions for different regions by analytic continuation in order to fix the integration constants. We will now discuss how to construct domain-restricted but well-behaved ansätze of multiple polylogarithms, which don't require an explicit +i 0 prescription and perform well numerically.

Optimizing the bases of multiple polylogarithms
Starting from a Duhr-Gangl-Rhodes basis of multiple polylogarithms for the integrals m, we wish to remove Li functions which have suboptimal analytic properties for physical kinematics (s > 0). In fact, one finds that even the purely rational function arguments allowed by the letters (4.10) are non-trivial to treat systematically and, we will therefore restrict the discussion in this section to this rational subset. The main idea is to make a partition of the physical phase space into regions D i such that, inside each region, a solution to the differential equations may be constructed out of Li functions which never diverge or move off of their principal Riemann sheets for arbitrary phase space trajectories contained in D i . Although our primary goal is to show how to one can avoid supplying explicit +i 0 prescriptions for w and z, we also find it convenient to impose further aesthetic criteria on our polylogarithmic bases in order to simplify our results. For example, we find it useful in each ansatz to give precedence to those functions which do not involve the symbol letter 1 − w + w 2 , since this ensures that the master integrals which do not depend on 1 − w + w 2 are manifestly free of 1 − w + w 2 at the level of functions. 5 The first step of our analysis is to study the letters of the rational alphabet, L R , above and below the physical two-mass threshold at s = 4 m 2 . To proceed, we must determine under what conditions the logarithms of the letters either diverge or move off of their principal Riemann sheets. We find that, in practice, it is simplest to use Eqs. Reduce to work out how various constraints on polynomials of w and z map back to conditions on the original and more familiar kinematic variables, s, t, and m 2 . We find that Reduce works in an efficient way when the system of inequalities to be reduced is formulated in terms of real-valued parameters and the solution to the system does not involve roots of high-degree polynomials. Fortunately, these assumptions are always satisfied for the Li functions which have symbols built out of letters from L R . 6 To check whether it is possible to live without a +i 0 prescription, we must first understand for what values of s, t, and m 2 both the real and imaginary parts of the letters vanish simultaneously and for what values of s, t, and m 2 the imaginary parts of the letters vanish while their real parts happen to be negative: • ln(l 1 ) = ln(1 − w) diverges at the phase space boundary point s = 0.
• ln(l 2 ) = ln(−w) has no issues for m 2 > 0, except at the phase space boundary point s = 0 where it becomes ill-defined.
• ln(l 6 ) = ln(−z) has no issues for m 2 > 0, except at the phase space boundary point s = −t = 2 m 2 where it becomes ill-defined.
• ln(l 7 ) = ln(1 + z) diverges at the phase space boundary where t = 0 and at the two-mass threshold s = 4 m 2 .
• ln(l 8 ) = ln(1 − w z) diverges at the two-mass threshold s = 4 m 2 . 5 The presence or absence of the letter 1 − w + w 2 for particular master integrals is linked to the presence or absence of a one-mass threshold at s = m 2 . 6 In Section 6, the method described in this section is also used to prepare an ansatz involving Li functions which have symbols built out of letters from both LR and LA.
• ln(l 11 ) = ln(z − w) diverges when t = −m 2 for s ≥ m 2 and, in particular, at s = 4 m 2 . ln(l 11 ) also moves off of its principal Riemann sheet at s = − 4 m 2 t m 2 −t for −m 2 ≤ t ≤ 0. The next step is to determine from the above data how many regions of the physical phase space it makes sense to consider separately. In a sense, this step is the most nontrivial because different choices besides the one we ultimately make are possible. Our choices are guided by imposing analytic properties of the integrals onto the basis functions themselves. Due to the absence of u dependence in our planar double boxes, we expect them to have a regular limit as u approaches 0. We see from the above that, depending on whether s is less than or greater than 2 m 2 , regularity as u → 0 suggests the absence, respectively, of ln(l 9 ) or ln(l 10 ) at weight one. Of course, we wish to achieve that also the higher weight basis functions lack logarithmic singularities as l 9 → 0 for 0 < s < 2 m 2 and l 10 → 0 for s > 2 m 2 . This can be achieved by choosing the point s = 2 m 2 as a region boundary and then imposing an appropriate first-entry condition [58,[61][62][63] on the symbol to select suitable functions. In other words, we remove from consideration all Li functions which have the letter l 9 in the first entry of their symbols in 0 < s < 2 m 2 regions and all functions which have the letter l 10 in the first entry of their symbols in s > 2 m 2 regions. All of the other letters exhibit problems only at the physical one-and two-mass thresholds or at other exceptional points on the boundary of the physical phase space like s = 0. The only truly sick logarithm is ln(l 11 ), which can be dealt with by imposing an additional first-entry condition on l 11 in all regions.
Ultimately, we find it natural to partition the physical phase space into four regions: and where, as suggested by Eqs. Even for the twelve two-mass master integrals which have rational symbols in the (w, z) representation, it is not obvious that it suffices to consider a partition of the physical phase space into just four separate regions. Fortunately, no further subdivisions are necessary 7 and it is even possible to consistently impose stronger constraints in each region on the subset of Li functions which survive our first-entry cuts.
Our basis of Li functions in region D a may be further refined along the lines described in [13], where the EW-QCD Drell-Yan master integrals with a single massive line were evaluated in the physical region for the first time. Consider, for example, the set of 192 Li n arguments consistent with our first-entry conditions in this region: l 1 , l 2 , −l 2 , l 3 , l 4 , l 5 , l 6 , −l 6 , l 7 , l 8 , l 9 , for our final set of preferred Li n function arguments above the two-mass threshold. We can proceed analogously to filter the set of Li 2,1 , Li 3,1 , and Li 2,2 function argument pairs which survive our first-entry cuts. For the sake of brevity, let us simply state that we find 2496 preferred pairs of function arguments (f i , f j ) after using Reduce to throw away all pairs for which f i > 1 or f i f j > 1 for (w, z) in D a . As explained in [13], we can further improve the numerical performance of our functional basis in GiNaC by reordering the lists of function arguments obtained to ensure that those which lead to convergent power series expansions are given precedence. In region D b 1 , the first-entry conditions are the same as for D a and again lead to (5.16). As mentioned above, we can use Reduce to efficiently check which Li n basis functions have suitable analytic properties below the two-mass threshold as well; despite the fact that w and z become complex functions of the original kinematic variables in D b 1 and the remaining below-threshold regions, we can proceed as before provided that we first eliminate [w] and [z] via Eq. (5.15) and work only with [w] and [z]. It turns out that a surprisingly large number of exceptional phase space points and trajectories cause problems for particular basis functions in the (w, z) representation. Specifically, we find that certain Li n arguments assume values on the real axis greater than one • when s = 3 m 2 for −3 m 2 < t < 0.
We also observe that a number of function arguments have imaginary parts which vanish identically below threshold. Any such Li n function with an argument which can attain values greater than one on D b 1 must be discarded because it will be ill-defined unless small imaginary parts are assigned to w and z.
In D b 2 , besides those functions whose arguments have imaginary parts which vanish identically, we find that certain Li n arguments assume values on the real axis greater than one for − m 2 3 < t < 0. After discarding all of the offending Li n function arguments, we find that 184 preferred ones remain: l 1 , l 2 , −l 2 , l 3 , l 4 , l 5 , l 6 , −l 6 , l 7 , l 9 , In much the same way, we find 3420 preferred pairs of function arguments for Li 2,1 , Li 3,1 , and Li 2,2 in D b 2 . In D b 3 , besides those functions whose arguments have imaginary parts which vanish identically, we find that certain Li n arguments assume values on the real axis greater than one • when s = − 4 m 2 t m 2 −t for − m 2 3 < t < 0.

Outlook
In this article, we considered irreducible Feynman integrals satisfying d ln differential equations with non-rational symbol letters and described methods for their evaluation in terms of standard multiple polylogarithms. We show for the first time that a complete solution in terms of standard multiple polylogarithms can be obtained even in the presence of unrationalizable symbol letters. In particular, new techniques for the construction of an ansatz which matches the symbol in a particular region of phase space allowed us to calculate the two-loop master integrals for the mixed EW-QCD corrections to Drell-Yan lepton pair production in the physical region using Li functions through to weight four. We discussed in detail how to optimize the functional basis to allow for fast and stable numerical evaluations, systematically avoiding arguments on branch cuts which would require an explicit +i 0 prescription. The presented techniques and results may be applied directly to the calculation of the virtual amplitudes for the full EW-QCD corrections to Drell-Yan production of relative order αα s . As motivation, we also considered a rather non-trivial master integral for massive Bhabha electron-positron scattering and showed that it may be directly integrated from Feynman parameters to all orders in in terms of multiple polylogarithms despite the presence of root-valued symbol letters. Besides the applications discussed in this paper, d ln differential equations with root-valued symbol letters also appear in other interesting contexts, such as the next-to-leading order QCD corrections to H+jet production with full heavy top quark mass dependence [18]. One may therefore expect that our techniques could be fruitfully applied also to other problems of current phenomenological interest.