The massless single off-shell scalar box integral -- branch cut structure and all-order epsilon expansion

We investigate the single off-shell scalar box integral with massless internal lines in dimensional regularization. A special emphasis is given to higher orders in the dimensional regularization parameter epsilon, its branch cut structure, and kinematic limits. Common representations of the box integral introduce superficial branch cuts, which we eliminate to all orders in the epsilon expansion. In the literature so far a satisfactory solution for this issue only exists up to finite order in the epsilon expansion. However, for calculations at NNLO in perturbation theory, higher orders in epsilon of this integral are required. In this paper, we present results to all orders in epsilon in terms of single-valued polylogarithms and explicitly determine the real and imaginary part of the box integral in all kinematic regions.


Introduction
With the upcoming electron-ion collider (EIC) [1], it is essential to evaluate the hardscattering coefficients for semi-inclusive deep-inelastic scattering (SIDIS) at next-to-nextto-leading order (NNLO) in QCD [2][3][4][5]. One type of contributions arising at NNLO are the one-loop corrections to the next-to-leading-order processes of real gluon emission and photon-gluon fusion. One-loop integrals such as these arise in virtually all higher-order perturbative calculations and were thoroughly studied in the past 40 years [6][7][8][9][10][11][12][13][14][15][16]. Given that tensorial loop integrals can be systematically reduced to scalar integrals using the Passarino-Veltman technique [17], it is sufficient to evaluate scalar integrals. The case of massless propagators is of particular interest for perturbative QCD calculations, since light quark masses can be neglected in high-energy scattering processes and gluons are massless. Loop integrals with vanishing internal masses generally have soft and collinear divergences, which must be regularized. The most commonly used method is dimensional regularization, which simultaneously serves to regularize ultraviolet and infrared divergences [18][19][20][21]. We will work in d = 4 − 2ε dimensions.
Both real and imaginary part of a loop integral do in general yield contributions to the absolute value squared of a matrix element, hence we must know both. The imaginary part of a loop integral arises in certain kinematic regions and will be governed by the causal +i0 from the propagators. For massless two-and three-propagator scalar loop integrals in general dimensions d, it is easy to keep track of this infinitesimal imaginary part. 1 Following the calculation given in ref. [16], one finds that the external momentum invariants are essentially replaced by (p i ± p j ) 2 → (p i ± p j ) 2 + i0 in the final results. This simple rule yields no unambiguously well-defined imaginary part for the four-propagator box integral in general dimension d, as we will discuss in section 2 below.
The scalar box integral with one off-shell external particle (i.e. with non-vanishing momentum squared) is especially important since it is needed for higher order (NNLO) perturbative calculations of processes such as SIDIS, Drell-Yan, or single inclusive annihilation (SIA), which all feature a single off-shell gauge boson. In ref. [22], this integral and its imaginary part were first given up to finite order ε 0 in dimensional regularization. The result for general dimensions d written in terms of three Gauss hypergeometric functions was obtained in refs. [23,24]. However, the simple prescription (p i ± p j ) 2 → (p i ± p j ) 2 + i0 given in ref. [24] for analytic continuation of the result to all kinematic regions only applies to the expansion up to finite order. Ref. [23, eq. (D.4)], which was used in the calculation of DIS structure functions at NNLO [25], does not touch on the question of analytic continuation. In [16], the box integral was expanded to all orders in ε, however the branch cut structure was also not discussed. The causal +i0 from the propagators was systematically kept in ref. [26], but the result was only given up to finite order in dimensional regularization. To regularize divergences occurring in kinematic limits, however, we need the full ε-dependence of any terms divergent in these kinematic limits.
In this paper, we independently recalculate the massless single off-shell scalar box integral for general dimensions d. We explicitly keep the causal +i0 from the propagators, such that our results are valid for arbitrary (real) values of the kinematic variables. This allows us to explicitly determine real and imaginary part to all orders of the ε-expansion, which goes beyond the finite O(ε 0 ) results of refs. [24,26] and the Mathematica implementation in Package-X [14]. In section 2, we calculate the internal-massless single-off shell scalar box integral using the well-known Feynman parametrization approach and an elegant factorization which allows for a trivial reduction to one-dimensional integrals. As a limiting case our result also includes the internal-massless box integral with on-shell massless external particles, which is needed for direct photon production and Compton-like scattering processes. In section 3, we perform the ε-expansion to all orders and the explicit determination of real and imaginary parts. The special cases of vanishing Mandelstam variables, as well as several identities for the Gauss hypergeometric function, and the definition of single-valued polylogarithms, are collected in the appendices. Additionally, we provide an example illustrating the need of higher orders in ε of the scalar box integral for phase space regularization in appendix D.

Calculating the massless single off-shell scalar box integral
The massless single off-shell scalar box integral is defined through the loop integral where the external momenta are labelled as indicated by the Feynman diagram depicted in figure 1. It is p 2 1 = p 2 2 = p 2 3 = 0 and q 2 = (p 1 + p 2 + p 3 ) 2 = 0. Note that we take all external momenta to be incoming. We have also explicitly kept the causal +i0 in the propagators, which is necessary to determine the physical side of the branch cuts that appear in the l l + p 2 Figure 1. The massless scalar box diagram with one off-shell external particle carrying momentum q. It is p 2 1 = p 2 2 = p 2 3 = 0 and q 2 = 0. All external momenta are taken to be incoming. Drawn with Tik Z-Feynman [27]. final result. 2 The prefactor is chosen as in ref. [28].
To evaluate this loop integral, we first combine the propagators into one generalized propagator of higher power using the well-known Feynman parametrization [29, eq. (A.39)] Next, we complete the square with respect to l and then shift the loop momentum in the usual manner, which leads to where we introduced the Mandelstam variables results in (2.8) Note that this expression is symmetric under exchanging s 1 ↔ s 2 . To reflect this, we write the scalar box integral as D 0 (s 1 , s 2 , q 2 ), where q 2 = s 1 + s 2 + s 3 is the squared momentum of the off-shell external particle. To decouple the Feynman parameter integrals, which are currently coupled via the delta function, we use the substitution The Jacobian of this substitution, which was inspired by ref. [30, page 44], is η 1 η 2 . Through the substitution, the η-and ξ-integrals factorize, (2.11) Here, the delta function was used to evaluate one of the η-integrals. The remaining η-integral thus became a representation of the Beta function, which was evaluated in terms of Gamma functions, For the resulting integral to safely converge, we must take ε < 0. 3 We will first evaluate the integral in this area of convergence and subsequently analytically continue the end result to larger ε. So far, our calculation did not differ much from the way the scalar box integral was calculated before (see e. g. ref. [26]) and was mainly presented for pedagogical reasons and introduction of notation. We will now deviate from the usual path by factoring the term s 1 s 2 s 3 out of the integrand. This is of course not sensible for s 1 s 2 s 3 = 0, ±∞, that is s 1 , s 2 , s 3 = 0. The special cases of one or several vanishing Mandelstam variables will be considered in appendix A. Using the identity (compare [26, eq. (19) which is not completely trivial since α is not an integer, we obtain (2.14) The i0 prescription determines on which side of the branch cut the respective factors should be evaluated if their argument is negative. The remaining integral on the second line only depends on the two dimensionless variables which is why factoring the term s 1 s 2 s 3 out of the integrand is advantageous. By substituting ξ 1 → ζ 2 = x 2 ξ 1 and ξ 2 → ζ 1 = x 1 ξ 2 , all dependence on the dimensionless x i is moved into the integration boundaries, where Note that we inverted the prefactor using Additionally, we introduced the abbreviation And furthermore, we have written the integrand in a form suggesting the substitution ζ i → 1 − ζ i . Subsequently substituting ζ 2 → ζ 12 = 1 − ζ 1 ζ 2 allows us to evaluate the resulting ζ 12 integral, While this integrand is finite at ζ 1 = 0, we must regularize the divergence at ζ 1 = 0 if we split the integral into two separate integrals. For this, we subtract and add 1 and subsequently substitute ζ 1 → ζ = 1 − ζ 1 (1 − x 2 ) in the first integral and ζ 1 → ζ = 1 − ζ 1 in the second integral, Splitting the first of these integrals in two, we find that the scalar box integral is given by Note that we have rewritten the Gamma functions in the prefactor such that they become unity in the limit ε → 0, using Γ(z + 1) = z Γ(z). Since the result for the internal-massless single off-shell scalar box integral is commonly given in terms of Gauss hypergeometric functions, we will identify them here, even though this will introduce spurious additional branch cuts. The three integrals appearing in eq. (2.22) are of the form This is trivially zero if χ = 0. All three integrals have vanishing χ for s 3 = 0. As mentioned above, this special case of vanishing Mandelstam variables will be considered separately in appendix A. The upper boundary of the third integral, also vanishes if q 2 = 0. We will account for this by setting the third integral to zero in case q 2 = 0. In case χ = 0, we can substitute ζ → χ −1 ζ, As ζ > 0, we could factor it out of the square brackets here. To split this integral into two separate integrals, we must regularize the divergence that occurs for ζ = 1/χ if χ > 1. We regulate it by introducing an infinitesimal imaginary part χ → χ + i0 in the denominator, with0 taken to be positive. 4 Note that the end result will not depend on the regulator i0, since the original integral I(χ) was not divergent for ζ = 1/χ. Proceeding as described, we obtain Here, the first integral became a representation of the Gauss hypergeometric function discussed in appendix B (see eq. (B.2)), while the second integral could be elementarily evaluated in terms of a logarithm. Both hypergeometric function and logarithm are evaluated on the branch cut for χ > 1. The regulator i0 tells us which side of these branch cuts to choose. Also note that the prefactor of the hypergeometric function vanishes in case χ = 0, since we assumed ε < 0 for D 0 to converge. Applying eq. (2.26) to the term in curly braces in eq. (2.22) and using a different regulator i0 i for each of the three integrals gives Writing the logarithms on the last line of this expression in terms of their real and imaginary parts, we see that only the imaginary parts yield a contribution, Here, Θ(x) is the Heaviside step function. Let us now choose the sign of the regulators i0 i such that the contribution of the logarithms completely cancels. The kinematic regions in the x 1 x 2 -plane in which the respective regulators yield a contribution are depicted in figure 2. For the right-hand side of eq. (2.28) to vanish, the regulators must satisfy In case both x 1,2 < 1 (white region), the right-hand side of eq. (2.28) vanishes independently of the choice of regulators. Note that since these conditions only determine the relative signs of the regulators, we can choose the absolute sign of one regulator. Let us select i0 3 ≡ +i0 everywhere. Then, both i0 1 and i0 2 must change their sign in the yellow-blue region in order to satisfy the conditions. To preserve symmetry, we let them change signs along the dashed diagonal depicted in figure 2, such that Combining our results, we find that the internal-massless single off-shell scalar box integral is given by Remember that the term in curly braces is a function of the two variables x 1 and x 2 defined in eq. (2.15). The global prefactor (. . .) ε and the prefactors [. . .] −ε of the hypergeometric functions can be recombined by first inverting the factor in square brackets using eq. (2.18) and subsequently applying eq. (2.13) from right to left, resulting in For convergence, we required ε < 0 during the calculation. For the expansion about the physical four dimensions, the result is analytically continued to larger ε.
The result was calculated for s 1 , s 2 , s 3 ∈ R\{0}. In appendix A, we will see that it is also valid for vanishing Mandelstam variables in the sense of appropriate limits. As discussed underneath eq. (2.24), the term on the last line is equal to zero for q 2 = 0, since ε < 0. Remember that in dimensional regularization kinematic limits need to be taken before analytic continuation in ε to assure consistent results.
For the kinematic variables, the causal i0 and the regulator i0 specify on which side of the occurring branch cuts the respective functions are evaluated. Remember that by construction this result does not depend on the regulator i0. Hence, the individually discontinuous imaginary parts of the hypergeometric functions must cancel each other, and the branch cuts of the hypergeometric functions are spurious. We will see that this is indeed the case in eq. (3.19) below.
In the region where s 1 , s 2 , q 2 < 0 and all three hypergeometric functions are away from their branch cut, our result agrees with those found in [ Expanding the hypergeometric functions in ε gives a result in terms of polylogarithms (see eq. (B.3)). As long as one avoids the branch cuts of the polylogarithms, the analytic continuation works via s 1 → s 1 + i0, s 2 → s 2 + i0, q 2 → q 2 + i0 as given in [24]. However, treating the branch cuts of the polylogarithms by this prescription leads to inconsistent results. Remarkably, the prescription works correctly for the result expanded up to O(ε 0 ) if the Abel identity [32] is applied to the dilogarithms first. Also the analytic continuation procedure introduced in [9] and put forward in [28] only applies to the dilogarithm as it is based on the Euler identity [32]. It is worth mentioning that all these prescriptions for analytic continuation make use of the causal +i0 and apply to each term individually. However from eq. (2.36) it is clear that the imaginary parts of the hypergeometric functions respectively polylogarithms are determined by a distinct regulator i0 and are interdependent.

Epsilon expansion of the scalar box integral
The scalar box integral given in eq. (2.35) can be expressed as where D 0 abbreviates the term in curly braces. By defining we can write where the upper sign choice is valid for x 1 ≥ x 2 and else the lower sign choice. To expand the scalar box integral around ε = 0, we need the ε-expansion of F ± (ε; x), which is performed in section 3.1. Using this key ingredient, we subsequently discuss the expansion of D 0 (ε; x 1 , x 2 ) in section 3.2 and of D 0 (s 1 , s 2 , q 2 ) in section 3.3. Tabulated results for the real and imaginary parts of exp iπε Θ − s 3 in the various kinematic regions are given in table 1 in section 3.3.

Epsilon expansion of F ± (ε; x)
To perform the ε-expansion of F ± (ε; x), which was defined in eq. (3.2), we observe that (x − i0 sgn 123 ) −ε develops an imaginary part for x < 0, while the hypergeometric function does so for x > 1. Hence, there are 3 regions of interest: where the ε-expansion of the hypergeometric function given in eq. (B.3) was used in the first step. Next, we substitute n → N = n + k and then interchange the N -and k-sum, For N ≥ 2, we now introduce the single-valued polylogarithms L n (x) discussed in appendix C by adding and subtracting Note that since the single-valued polylogarithms are bounded on R, any divergences in the kinematic variable x are explicitly contained in the logarithms.

Expansion for
develop an imaginary part. Following the same steps as before, we find Again, any divergences in the kinematic variable are contained in the logarithms since the single-valued polylogarithms are bounded on R.

Expansion for x ∈ (1, ∞)
In the area where x > 1, the hypergeometric function 2 F 1 1, −ε, 1 − ε, x ± i0 develops an imaginary part. To make it explicit, we employ the inversion formula given in eq. (B.11), By replacing ε → −ε, we can use the ε-expansion from eq. (B.3) for the resulting hypergeometric function. Additionally writing the last term in terms of its real and imaginary part, we obtain The ε-expansion of the cotangent term is given by [31, eqs where ζ n = ζ(n) = ∞ k=1 k −n are the integer values of the Riemann zeta function. Additionally expanding |x| −ε and collecting orders of ε yields Upon adding and subtracting the term (−1) n−1 n! ln n−1 1 x ln 1 − 1 x , we can identify the single-valued polylogarithms L n 1 x , as defined in eq. (C.6). Inverting them with the help of eq. (C.7), we find (3.12) Organizing this expression in terms of logarithms and single-valued polylogarithms yields This compact expression is well suited to be used for the ε-expansion of the complete box integral. The imaginary part is made explicit and all functions of the real variable x are manifestly real. The complex exponential stems from (x − i0 sgn 123 ) −ε , while the imaginary part πε Θ(x − 1) originates from the hypergeometric function. An important feature of the expansion of F as given in eq. (3.15) is that it is finite in the limit x → ±∞ in every order in ε. Using the limits of the single-valued polylogarithms given in appendix C, we obtain By collecting the ln |x| terms we can cast eq. (3.15) into the alternative form From this we can read off that F(ε; x) behaves like |x| −ε near x = 0 and it diverges logarithmically like ε ln |1 − x| near x = 1 . In the complete box integral the ln |1 − x| divergences always cancels.

Epsilon expansion of
By plugging the ε-expansion of F ± (ε; x) from eq. (3.14) into eq. (3.3), it is straightforward to obtain the ε-expansion of D 0 (ε; x 1 , x 2 ). We must consider seven areas corresponding to different kinematics, which are depicted in figure 3. Since D 0 (ε; Observing that we see that the imaginary parts originating from the hypergeometric functions cancel in all kinematic regions, as they should by construction. Hence, we obtain The ε-expansion of D 0 in the first region, where none of the three parts develops an imaginary part, is purely real. The logarithmic poles at x 1,2 = 1 cancel between the three F-functions, such that s 1 s 2 D 0 (ε, x 1 , x 2 ) in the areas specified by the signs of s 1 , s 2 and q 2 . For compactness we introduce the abbreviations F 1 ≡ F(ε; x 1 ), F 2 ≡ F(ε; x 2 ) and For the order ε 2 coefficient, we can use Abel's formula (C.16), in combination with the Euler identity (C.11) and inverted Landen identity (C.15), to reduce the number of appearing L 2 to just two, where we assumed x 1 , x 2 , (x 1 +x 2 −x 1 x 2 ) > 0. For the higher orders in ε there presumably exists no representation with fewer than three (single-valued) polylogarithms, since only for the dilogarithm there exists a suitable five term relation (see appendix C.3).

Epsilon expansion of
Hence, the imaginary part is determined by the signs of s 1 , s 2 and q 2 . Explicit results for the real and imaginary part in all areas are given in table 1. Note that the superficially divergent terms ln(1 − x 1 ) and ln(1 − x 2 ) cancel between F(ε; x 1,2 ) and F(ε; in the regions where x 1,2 = 1. Combining the results listed in table 1 and eq. (3.1) we can directly infer real and imaginary parts of the complete box integral D 0 (s 1 , s 2 , q 2 ). For completing the ε-expansion we can employ [31, eq. (6.1.33)] for all Gamma functions in the prefactor, where γ E is the Euler-Mascheroni constant. However, since the same factor of Gamma functions multiplies all one-loop integrals and a similar factor is present in corresponding phase-space integrals it is often not required to expand all of them explicitly. Expanding the term s 3 µ 2 s 1 s 2 ε introduces additional logarithmic terms in the ε-expansion. One might albeit need to keep this term unexpanded to regularize singular behavior in the kinematic limits s 1 , s 2 , s 3 → 0. This is exemplified in appendix D.
To sum up our result, we can cast D 0 (s 1 , s 2 , q 2 ) as given in eq. (2.35) into the form Plugging in the ε-expansion of F(ε; x) from eq. (3.15) we have an all-order ε-expansion of D 0 (s 1 , s 2 , q 2 ) with explicit real and imaginary parts valid in the entire kinematic domain of s 1 , s 2 , q 2 . Since the term in square brackets is finite for vanishing s 1,2 , D 0 (s 1 , s 2 , q 2 ) behaves like |s 1,2 | −1−ε in the limits s 1,2 → 0. For a comprehensive discussion of the behavior of F(ε; x) in different kinematic limiting cases we refer to section 3.1.4.
By bringing the factor s 3 µ 2 s 1 s 2 ε into the brackets, we obtain a form analogous to eq. (2.36),

Conclusion
We have revisited the massless scalar box integral with one external off-shell particle, keeping the causal +i0 throughout the calculation. The main new result is extending the explicit determination of real and imaginary parts beyond finite order in the dimensional regularization parameter epsilon in all kinematic regions. The epsilon expansion is expressed to all orders in terms of newly introduced single-valued polylogarithms. This makes our result free of the spurious branch cuts common to representations using ordinary polylogarithms. The single-valued polylogarithms are bounded, such that any divergent behavior in kinematic limits is contained in logarithmic terms. Since we have determined the box integral in all kinematic regions, our result is suitable for all kinds of calculations in higher orders of perturbation theory, where imaginary parts or higher orders in epsilon become important. The calculation method put forward in this work can be generalized to the scalar box integral with two non-adjacent end points off the light cone. A corresponding paper is in preparation [33].

Acknowledgments
We are grateful to Werner Vogelsang for valuable comments regarding this work. We thank the referee for raising the question if generalizing the result to more end points off the light cone is possible. This study was supported in part by Deutsche Forschungsgemeinschaft (DFG) through the Research Unit FOR 2926 (Project No. 409651613). J. H. is grateful to the Landesgraduiertenförderung Baden-Württemberg for supporting her research.

A Scalar box integral in case of vanishing Mandelstam variables
In this appendix, we evaluate the internal-massless single off-shell scalar box integral in case one or several Mandelstam variables vanish, since the result calculated in the main text was obtained for non-vanishing Mandelstam variables. For this, we go back to D 0 as given in eq. (2.11) and explicitly set the vanishing Mandelstam variables to zero. We will find that all cases agree with appropriate limits of the general result.
A.1 s 1 = 0 or s 2 = 0 Since the scalar box integral is symmetric under exchanging s 1 ↔ s 2 , we can without loss of generality consider s 1 = 0. In this case, the integrals in eq. (2.11) factorize, Here the ξ 2 -integration was performed elementarily. Note that ε < −1 was imposed for convergence. While the ξ 1 -integral can in principle be evaluated elementarily as well, first substituting ξ 1 → ζ = −s 2 − s 3 ξ 1 is useful for discussing the additional limits q 2 , s 2 = 0. Splitting the resulting integral into two separate integrals, we obtain Note that the second integral and thus the second term on the last line vanishes for q 2 = s 2 + s 3 = 0. The result can be analytically continued to values of ε beyond the original area of convergence. Let us compare this result to the limit s 1 → 0 of the general result as given eq. (2.36). The limit of the first and last term in the general result are of the form To evaluate this limit, we first apply the inversion formula for the hypergeometric function derived in eq. (B.11) in the appendix, which yields eq. (A.3) = lim .
Imposing the condition ε < −1, which is needed for convergence in case s 1 = 0, the second term vanishes. The remaining limit can be evaluated using L'Hôpital's rule, eq. (A.3) = lim The derivative of the hypergeometric function was read off the series expansion given in eq. (B.1). The second term in the general result also vanishes for s 1 → 0 since ε < −1.
Combining everything, we find that taking the limit s 1 → 0 of the general result given in eq. (2.36) is equivalent to setting s 1 = 0 in eq. (2.11).
A.2 s 1 = 0 and s 2 = 0 In this case, we set s 2 = 0 in eq. (A.1) and obtain This is identical to the expression obtained by taking the limit s 2 → 0 of eq. (A.2), i.e. taking the limits s 1 , s 2 → 0 of eq. (2.36).

A.3 s 3 = 0
In this case, eq. (2.11) reads where the ξ 2 -integration was performed elementarily. Note that ε < −1 was imposed for convergence. Next, we substitute ξ 1 → ζ = −s 1 ξ 1 in the first term and ξ 1 → ζ = −s 1 ξ 1 − s 2 in the second term, and subsequently split the latter integral into two separate integrals, Note that the last integral vanishes for q 2 = 0. Evaluating the remaining ζ-integrals yields which can be analytically continued beyond ε < −1 which was originally required for convergence. One can easily see that this result is equivalent to setting s 3 = 0 in the general result as given in eq. (2.36), since the resulting hypergeometric functions evaluate to 2 F 1 (1, −ε, 1 − ε; 0) = 1 according to eq. (B.1).
A.4 s 1 = s 3 = 0 or s 2 = s 3 = 0 As above, we can without loss of generality consider s 1 = s 3 = 0. In this case, eq. (2.11) yields (A.10) As in the previous cases, one has to impose ε < −1 for convergence before analytically continuing to larger ε. Let us compare this expression to the corresponding limit of the general result as given in eq. (2.36). For s 3 = 0, the hypergeometric functions evaluate to 2 F 1 (1, −ε, 1 − ε; 0) = 1.
The limit remaining to be taken is lim Imposing the condition ε < −1, which is needed for convergence in case s 1 = 0, the second term vanishes. Applying L'Hôpital's rule to the first limit, we find that the limit s 1 , s 3 → 0 of the general result given in eq. (2.36) is equivalent to setting s 1 = s 3 = 0 in eq. (2.11).
A.5 s 1 = s 2 = s 3 = 0 Setting all Mandelstam variables to 0 in eq. (2.11) and requiring ε < −2 for convergence yields This integral also has to vanish since the integrand has a mass dimension of −ε − 2, but does not depend on any mass scale. Taking the limit s 1,2,3 → 0 of the general result given in eq. (2.36) leads to the same conclusion.
Note that the integral may be ill-defined for z > 1, since in this case the term (1 − zt) changes sign within the range of integration. This change of sign can, in general, give rise to a discontinuity or branch cut. In the following, we discuss selected properties (including the branch cut structure) of the specific hypergeometric function 2 F 1 (1, −ε, 1 − ε; z) which appears in the internal-massless single off-shell scalar box integral.

B.1 Epsilon expansion of
The ε-expansion of this hypergeometric function was derived in ref. [16, eq. (B.28)], Since the polylogarithms Li n (z) have a branch cut along the positive real axis greater than 1 [32], the hypergeometric function also has a branch cut there. The branch cut structure will be discussed in section B.3.

B.2 Inversion relation for
An equation relating the hypergeometric function 2 F 1 (1, −ε, 1 − ε; z) to another hypergeometric function with argument 1 z is useful because it allows us to transform the argument away from the branch cut. To find this inversion relation, we write the hypergeometric function in terms of its ε-expansion and invert the appearing polylogarithms using the inversion relation for the polylogarithm given in ref. [32, eq. (7.20)], which holds for z ∈ C \ [0, ∞). 5 Plugging the second term on the left into eq. (B.3) leads to the closely related ε-expansion of 2 F 1 1, ε, 1 + ε; 1 z , while plugging the second term on the right into eq. (B.3) leads to the series expansion of an exponential function, such that To simplify the term on the second line of this expression, we interchange the two sums and subsequently shift the summation index to N = n − 2k, such that the N -sum becomes an exponential function, Thus, we obtain 5 For reference, the corresponding equation on the branch cut, i.e. for x ∈ (0, ∞), is Lin(x ± i0) + (−1) n Lin 1 The remaining polylogarithm constants can be written as [32, eq. (7.29)] where B 2k are the Bernoulli numbers, which are given by the coefficients in the following series expansion 6 Using this, we find that the term in square brackets is the series representation of πε csc(πε) = πε/ sin(πε) [31, eq. (4.3.68)], . (B.10) Hence the inversion formula for the hypergeometric function is The hypergeometric function 2 F 1 (1, −ε, 1 − ε; z) has a branch cut along the positive real axis starting at 1. This is explicitly seen through the inversion formula derived in eq. (B.11).
If we take x ∈ R >1 and i0 as an infinitesimal imaginary part, then we have Note that the imaginary part could be dropped in the hypergeometric function on the right-hand side since the function is well-defined for arguments 1/x < 1. Explicitly writing the factor (−x ∓ i0) ε in terms of its real and imaginary part, we obtain C Single-valued polylogarithms

C.1 Single-valued polylogarithms in the literature
There are several definitions of so called single-valued polylogarithms in the literature. Their common feature is that they do not have branch cuts opposed to regular polylogarithms. Depending on the intended use, different versions prove themselves to be most suitable [34].
Bloch and Wigner devised an important version of a single-valued dilogarithm [35] D 2 (z) = Im [Li 2 (x)] + arg(1 − z) ln |z| , (C.1) defined on C \ {0, 1}. This was generalized to polylogarithms by Ramakrishnan [36]. With a slight modification for odd n due to Wojtkowiak [37] and the notation R n = Re / Im if n is odd / even, his definition reads Similarly Zagier proposed [38, eq. (33)] Both are single-valued real-analytic functions from C \ {0, 1} to R. Zagier's version satisfies the functional equations of the polylogarithm in a clean fashion without product terms. However, both Ramakrishnan's and Zagier's functions are identically zero on R for even n. Hence, they are not suitable for our purpose. Also in an attempt to get rid of product terms in the functional equations, Lewin defined his version of polylogarithms by [39, eq. (3.19)] Kummer's functions are real-valued for real arguments, which makes them single-valued polylogarithms predating Bloch and Wigner's construction by well over a century. We observe from eq. (C.5) that L n (x) can be continued to a single-valued function for x > 1 by replacing ln(1 − x) → ln |1 − x| in definition (C.4).

C.2 Definition and basic properties of continuous single-valued polylogarithms
We define for n ∈ N ≥2 and real x a version of single-valued polylogarithms by For x > 1 the polylogarithms are taken on their principal branch, which is inherited from the principal branch of the logarithm, i.e. Li n (x) ≡ Li n (x − i0). As a result we have a function of a real variable x with the following properties: (a) L n (x) is single-valued, i.e. there is no branch cut and no imaginary part, in contrast to Li n (x).
(b) L n (x) is continuous for all x ∈ R.
(c) L n (x) is bounded on R, in contrast to Li n (x).
(d) L n (x) satisfies clean versions of the functional equations of Li n (x), i.e. without product terms, as for Lewin's and Zagier's functions.
These properties make the function L n (x) well-suited for various applications, where one is interested in having a well-behaved version of polylogarithms with real arguments. Plots of the first few single-valued polylogarithms are shown in figure 4. The qualitative behaviour of the functions is different for even and odd n. Here, ζ n = ζ(n) = ∞ k=1 k −n are the integer values of the Riemann zeta function.
For both odd and even n, L n (x) is real analytic on R\{0, 1}. At x = 0 its derivative diverges logarithmically, at x = 1 it is (n − 2)-times differentiable, the (n − 1)th derivative diverges logarithmically at this point.

C.3 Functional equations of continuous single-valued polylogarithms
The functional equations satisfied by single-valued polylogarithms are inherited from the corresponding polylogarithms [32]. Single-valued polylogarithms obey these equations in a cleaner fashion, where all product terms vanish. There are equations satisfied for all n and additional identities for small values of n. We explicitly give the identities for di-and trilogarithm, identities for higher order polylogarithms can be found in [32,39].

C.3.1 Equations of all single-valued polylogarithms
For all n there is an inversion relation and a duplication formula.
Alternatively we can write this in the suggestive form (ii) Duplication formula:

C.3.2 Equations of the single-valued dilogarithm
For the single-valued dilogarithm, defined as there is the richest structure. We have identities connecting the arguments x, 1 − x, x x−1 , 1 x , 1 1−x and 1 − 1 x . Additionally there is the famous 2-variable 5-term relation discovered independently by Spence [42, p.9] and Abel [43,XIV.(9)]. All discontinuities in these identities are induced by the mapping of the discontinuity at infinity.

C.3.3 Equations of the single-valued trilogarithm
Besides the inversion relation, the single-valued trilogarithm, defined as satisfies a single-variable identity connecting three terms and a three-variable 22 term relation discovered by Goncharov [44]. From the latter there follow a plethora of two variable identities as special cases. Since the single-valued trilogarithm is continuous at infinity, these identities are continuous.
(i) Inversion relation: (ii) Landen identity: (iii) Goncharov identity: Cyc(x,y,z) where cyclic summation for a function f (x, y, z) is defined as Cyc(x,y,z) f (x, y, z) = f (x, y, z) + f (y, z, x) + f (z, x, y) . (C. 21) D Example: Regularizing kinematic limits for Drell-Yan differential in rapidity In this appendix, we illustrate the regularization of kinematic divergences occurring in the box integral using the Drell-Yan process differential in momentum transfer Q 2 and rapidity y at O(αα 2 s ) as an example. This process was first calculated to NNLO in [45] using the reversed unitarity method. Instead, for the purpose of our example, we consider the NNLO calculation using the traditional treatment of the phase space integral as was done for the NLO calculation in [46]. We will explicitly see that orders beyond ε 0 of D 0 are required and the full ε-dependence of any terms divergent in the kinematic limit is needed here. To illustrate these points it is sufficient to look at the contribution to dσ dQ 2 dy from the quarkantiquark process. At NNLO the box integral appears through the interference term of one-loop and tree level diagrams for qq → γ * g, which are depicted in figure 5. Figure 5. NNLO contribution to qq → γ * from interference of box and tree-level amplitudes. There are three different boxes and two different tree-level diagrams for this process. In the limit y → 1 the arguments of all three F functions become zero. From eq. (3.18) we infer that forF(ε; x) = F(ε; x) − |x| ε the limit lim x→0F (ε;x) x is finite order by order in ε. Therefore, we can write D(s, u, Q 2 ) ∼ (1 − z) −ε + e iπε y ε (1 − z −ε ) + (1 − y) ε F (ε; (1 − z)(1 − y)) + e iπεF ε; − 1 − y y − e iπεF ε; − (1 − y)z y . (D.8) Note that the term in square brackets stays finite in the limit y → 1 even when divided by (1 − y) and hence D(s, u, Q 2 ) = A su (y) + (1 − y) 1+ε B su (y) , (D. 9) with both A su (y) and B su (y) finite in the limit y → 1 order by order in ε.
To convert the singular behavior for y → 1 into poles in ε, which is required to combine real and virtual corrections ensuring cancellation of infrared poles, we use the expansion (see e.g. eq. (110) in [46] or eq. (69) in [47]) where the plus distribution is defined by We see that the order ε 1 of the box integral contributes to the order ε 0 of the full result, since the 1 ε -poles generated by eq. (D.10) multiply the order ε 1 of the functions A st (y), A su (y), and A tu (y). Hence, it is required to perform the ε-expansion of the box integral to one order higher than needed if no regularization would be required. In the case of multiple regularizations, e.g. for y, z → 1 as in eq. (111) of [46], multiple poles multiply the box integral, requiring knowledge of the box integral to even higher orders in ε. Furthermore it was crucial to sum up any terms of the box integral singular in the kinematic limit to all orders, i.e. summing up terms of the form ε n ln n (1 − y) to (1 − y) −ε . This changed 1 ε to 1 2ε in (i) and (iv).