Master integrals for the NNLO virtual corrections to $\mu e$ scattering in QED: the non-planar graphs

We evaluate the master integrals for the two-loop non-planar box-diagrams contributing to the elastic scattering of muons and electrons at next-to-next-to-leading order in QED. We adopt the method of differential equations and the Magnus exponential to determine a canonical set of integrals, finally expressed as a Taylor series around four space-time dimensions, with coefficients written as a combination of generalised polylogarithms. The electron is treated as massless, while we retain full dependence on the muon mass. The considered integrals are also relevant for crossing-related processes, such as di-muon production at $e^+e^-$ colliders, as well as for the QCD corrections to top-pair production at hadron colliders. In particular, our results, together with the planar master integrals recently computed, represent the complete set of functions needed for the evaluation of the photonic two-loop virtual next-to-next-to-leading order QED corrections to $\mu e \to \mu e$ and $e^+ e^-\to\mu^+\mu^-$.


Introduction
In a previous work [1], we began the investigation of the next-to-next-to-leading-order (NNLO) virtual corrections to the elastic scattering of muons and electrons in Quantum Electrodynamics (QED), by classifying and evaluating the planar two-loop integrals arising from Feynman diagrams at this order in perturbation theory [2].
The NNLO QED corrections to the process µe → µe are crucial to interpret the highprecision data of future experiments like MUonE, recently proposed at CERN, aiming at measuring the differential cross section of the elastic scattering of high-energy muons on atomic electrons as a function of the spacelike (negative) squared momentum transfer [3,4]. This measurement will provide the running of the effective electromagnetic coupling in the spacelike region and, as a result, a new and independent determination of the leading hadronic contribution to the muon g-2 [3,4]. In order for this new determination to be competitive with the present dispersive one, which is obtained via timelike data (see [5] for a review), the µe differential cross section must be measured with statistical and systematic uncertainties of the order of 10ppm. This high experimental precision demands an analogous accuracy in the theoretical prediction.
Moreover, the NNLO QED corrections for the crossing-related scattering process e + e − → µ + µ − are important for some of the high-precision studies planned at upcoming lowenergy e + e − experiments, like Belle-II and VEPP-2000. Two interesting applications would be the following. The forward-backward asymmetry in muon pair production could be exploited to constrain non-standard eeµµ interactions [6], and the current estimates suggest that the knowledge of the NNLO QED differential cross section is needed, as QED itself produces an asymmetry starting at NLO. The knowledge of the QED radiative corrections to the e + e − → µ + µ − cross section will also be needed for precise measurements of the ratio R(s) = σ (e + e − → hadrons) /σ (e + e − → µ + µ − ) [7,8].
In this work, we complete the task of determining all functions required by the NNLO QED virtual corrections to µe scattering, by evaluating the two-loop integrals coming from non-planar four-point Feynman diagrams. Given the hierarchy between the electron mass m e and the muon mass m, m e /m ∼ 5 · 10 −3 , as in our former study we consider the approximation m e = 0. For the non-planar topology, integration-by-parts identities (IBPs) [9][10][11] yield the identification of a set of 44 master integrals (MIs), which we compute analytically by means of the differential equations method [12][13][14]. The system-solving strategy [15,16] is based on a consolidated procedure, which has been proven to be very effective in the context of multi-loop integrals involving several scales [1,[16][17][18][19]. Starting from the choice of a set of MIs obeying a system of first-order differential equations (DEQs) in the kinematic variables s/m 2 and t/m 2 which is linear in the space-time dimension d, we derive by means of the Magnus exponential matrix [16], an equivalent system of equations in canonical form [15], where the d-dependence of the associated matrices is factorised from the kinematics. The matrix associated with the canonical system is a logarithmic differential form which, in appropriate variables, has a polynomial alphabet. The canonical MIs can be therefore cast as a Taylor series around d = 4, with coefficients written as combinations of generalised polylogarithms (GPLs) [20][21][22][23].
For certain classes of MIs, like the ones of µe → µe and crossing-related processes, the choice of the boundary conditions may also constitute a challenging problem. Here, we exploit either the regularity conditions at pseudo-thresholds or the expression of integrals which are obtained by solving simpler, auxiliary systems of DEQs. Therefore, we limit the use of direct integration only to a small subset of simpler integrals used as input functions.
The package Reduze [24] has been used throughout the calculations, for the IBPs decomposition and for generating the DEQs obeyed by the MIs. The analytic expressions of the MIs have been numerically evaluated with the help of GiNaC [25] and were successfully tested against the numerical values provided either by the computer code SecDec [26] or, for the most complicated two-loop non-planar topologies (with 6-and 7-denominators), by an in-house algorithm. For such topologies, we identified an alternative set of quasi-finite integrals [27], more suitable for numerical integration, also with the help of Reduze.
The non-planar four-point functions hereby presented, together with the planar ones presented earlier [1], and the three-point functions available in the literature [28][29][30], make the analytic evaluation of the virtual two-loop amplitudes for the µe scattering in QED, as well as for the crossing-related processes, within reach.
In addition, we remark that the MIs of the QED corrections to e + e − → µ + µ − are a subset of those needed for the QCD corrections to the tt-pair production at hadron colliders. The complete two-loop QCD corrections to pp → tt are currently known only numerically [31][32][33][34][35]. The analytic evaluation of the MIs appearing in the leading-colour corrections to pp → tt, due to planar diagrams only, were considered in refs. [36][37][38][39]. In our former work [1], we extended the set of available functions for considering also sub-leading colour contributions. Very recently, the analytic calculation of the MIs for the planar doublebox integral with a closed top loop appeared in [40]. The analytic result for a non-planar three-point function, which constitutes a sub-diagram of the non-planar double box with closed heavy quark loop, was presented in [41]. The non-planar graphs hereby considered would also contribute to subleading-colour terms to tt-pair production, and their analytic evaluation was never addressed before. The numerical evaluation of (one of the) non-planar integrals computed analytically in this work has been recently considered in [42], in the context of a novel, promising method that aims at the numerical solution of differential equations.
The paper is organised as follows. In section 2, we set our notation and conventions for the four-point topology relevant for µe scattering. In section 3, we describe the general features of the systems of DEQs satisfied by the MIs, cast in d log-form, and we present the results for the non-planar two-loop MIs. In section 4, we describe the numerical evaluation of the non-planar four-point integrals. The information provided in the text is complemented by two appendices: in appendix A, we discuss the computation of the auxiliary integrals which have been used to extract some of the boundary constants and, in appendix B, we give the matrices associated with the d log-form.
The analytic expressions of the considered MIs are given in the ancillary files accompanying the arXiv version of this publication.

The non-planar four-point topology
In this paper, we consider the µe scattering process µ + (p 1 ) + e − (p 2 ) → e − (p 3 ) + µ + (p 4 ) , (2.1) in the approximation of vanishing electron mass, m e = 0, i.e. with kinematics specified by where m is the muon mass. Representative Feynman diagrams of the 10 relevant twoloop four-point topologies T i that contribute to the process are depicted in figure 1. The computation of the MIs belonging to the topologies T 1,2,3,4,5,7,8,9,10 has been discussed in [1]. In this paper, we complete the evaluation of all MIs required for the two-loop virtual amplitude, by determining the analytic expression of the MIs that belong to the non-planar topology T 6 . The calculation involves the evaluation, in d dimensions, of Feynman integrals of the type  is carried over around d = 6. We set ≡ (d * − d)/2, where d * = 4 and d * = 6 according to the case considered, and define our integration measure where µ is the 't Hooft scale of dimensional regularisation and Γ ≡ Γ(1 + ). Notice that our integration measure, when evaluated at d = 4 − 2 , agrees with eq.(3.2) of [1].

System of differential equations
By means of IBPs, the two-loop integrals that belong to T 6 can be reduced to a basis of 44 distinct MIs. In order to determine the analytic expression of the latter, we derive their DEQs in the kinematic variables s and t. The evaluation of the MIs can be further facilitated by parametrising the Mandelstam invariants in term of two independent dimensionless variables, w and z, which are defined by where the constraint s + t + u = 2m 2 is understood. Such change of variables rationalises the canonical DEQs. 1 A canonical basis of MIs in d = 4 − 2 can be identified by making use of the algorithm described in [16,17]. Namely, we start by choosing an initial set of MIs F i that fulfill DEQs with linear dependence on the dimensional regularisation parameter , where the T i are the integrals depicted in figure 2.
1 At an earlier stage of the project, we found that the variables x and y, defined through, remove all irrational terms appearing in the system of DEQs, individually. However, as pointed out by Lorenzo Tancredi -whom we acknowledge for the suggestion -, it is sufficient to rationalise just those combinations of irrational terms that appear in the DEQs, by means of w and z defined through in eq. (3.1), to yield a polynomial alphabet.  Figure 2: The 44 MIs T 1,...,44 for the two-loop non-planar topology T 6 . Thin lines represent massless propagators and thick lines stand for massive ones. Each dot indicates an additional power of the corresponding propagator. Numerator insertions are indicated explicitly on top of each diagram.
Subsequently, we use the Magnus exponential in order rotate the integrals of eq. (3.2) into a new basis of MIs I i that satisfy canonical DEQs in both variables w and z (or, equivalently, in s and t), where we introduced the abbreviation λ t = √ −t √ 4m 2 − t.
By combining the two DEQs in w and z into a single total differential, we get dI = dAI , (3.4) where I is a vector that collects the 44 MIs and with the M i being constant matrices with rational entries. The arguments η i of this d logform, which define the so-called alphabet of the DEQs, are the following 12 letters: In the present work, we compute the MIs in the kinematic region where all the letters are real and positive, which corresponds to the unphysical region t < 0 ∧ s < 0 . with the n-th order coefficient given by where I (i) (w 0 , z 0 ) is a constant vector and ∆ (k) the weight-k operator that iterates k ordered integrations of the 1-form dA along a piecewise-smooth path γ in the wz-plane. Since the rational alphabet given in eq. (3.6) has only algebraic roots, we can directly express (by first integrating in w and then in z) the iterated integrals of eq. (3.11) in terms of GPLs, which are defined as G( a n ; x) ≡ G(a 1 , a n−1 ; G( a n−1 ; t), (3.12) The length n of the vector a n corresponds the transcendental weight of G( w n ; x) and it amounts to the number of iterated integrations that define the GPL. The GPLs in our solution are of two classes, namely GPLs in w, with weights drawn from the set and GPLs in z, with weights drawn from In the region defined by eq. (3.7), the imaginary part of our solution I( , w, z) only originates from the integration constants I (i) (w 0 , z 0 ).

Boundary conditions
The general solution of the system of DEQs in terms of GPLs, which is obtained from the integration of eq. (3.4), must be complemented by a suitable set of boundary conditions. These boundary conditions can be determined either from the knowledge of the analytic expression of the MIs in special kinematic configurations or by imposing their regularity at pseudo-thresholds of the DEQs. For the problem under consideration, regularity conditions express the boundary constant as combinations of GPLs of argument 1, with weights drawn from the set a i ∈ {−1, −i, 0, i, 1}, which arises from the kinematic limits imposed on the alphabet given in eq. (3.6). We used GiNaC to numerically verify that for each MI, at each order in , the corresponding combination of constant GPLs is proportional to a uniform combination of the transcendental constants π, ζ k and log 2.
In the following, we specify how the boundary constants of each integral have been obtained: • The integrals I 1,..., 5,8,9,10,13,...,18,20,...,24,29,30,41 are common to the two-loop topologies discussed in ref [1], to which we refer the reader for the discussion of the boundary fixing. Furthermore, the integrals I 25,26 are related to I 23,24 by s ↔ u crossing, so that their boundary constants can be inferred directly from the ones of I 23,24 .
• The integrals I 6,7 are regular in the limit u → 0, where they can be reduced, via IBPs, to a single two-loop vacuum diagram. From the analytic expression of the latter, we obtain the boundary values I 6 | u=0 = 0 , (3.16) • The integrals I 11,12 are regular in the limit u → 0. In particular, we observe that their boundary values can be obtained as the limits (3.17) Therefore, we can generate the DEQs for the analogous triangle integrals with p 2 4 = 0, and u = 0 and an off-shell leg p 2 1 , solve them by using as an integration base-point the regular point p 2 1 = 0, and finally extract the boundary values of I 11 and I 12 by means of eq.(3.17). The details of this computation are reported in appendix A. In this way, we obtain • The boundary constants of the integrals I 27,28 can be fixed by imposing the regularity of their DEQs as t → 0.
• The boundary constants of the integrals I 31,32,33 are obtained by demanding regularity at t → 0, as well the reality of the integrals in the region s ≤ 0, u ≤ 0.
• The boundary constants of the integrals I 19,34,35... 40,42,43,44 are obtained by demanding their finiteness in the limit s → The analytic expressions of the MIs are given in electronic form in the ancillary files attached to the arXiv version of the manuscript.

Numerical evaluation of the non-planar four point integrals
The analytic expression of our MIs have been numerically evaluated in the region s, t < 0 by means of the GiNac library, and successfully checked against independent calculations. In particular, the integrals I i with i = 1, . . . , 36, 41 were computed with the package SecDec.
For the most complex topologies, corresponding to the non-planar four-point integrals I i with i = 37, . . . , 40, 42, 43, 44, we adopted a different strategy. As the numerical evaluation of those integrals is challenging, we identified an alternative set of independent MIs that are quasi finite [27] in d = 6. The latter have been computed semi-numerically by means of an in-house algorithm: starting from the Feynman parametrisation of the integrals, we carried out as many analytic integration as possible, until we reached a form where the left over multivariate integral could be numerically evaluated by means of Gauss quadrature.
Dimension-shifting identities [43,44] and IBPs, implemented in LiteRed [45,46], establish analytical relations between this set of integrals and the MIs we computed around d = 4. The definition of the 7 non-planar MIs that are quasi finite in d = 6 dimensions, together with our results at the phase-space point s = −1/7, t = −1/3, m 2 = 1, are collected in table 1. We identified them through educated guesses or with the help of Reduze. In the next subsection, we use the first of those integrals as an example to describe our evaluation strategy.
graph  Table 1: Numerical results for our set of quasi-finite non-planar MIs belonging to the 6and 7-denominators topologies (m 2 = 1).

The non-planar box in d = 6 dimensions
As an example, we describe the numerical evaluation of the non-planar scalar integral carried out in two steps.

Analytic integrations
By using Feynman parametrisation, the integral can be written as After integrating over k 1 and k 2 , one finds

4)
where we used the notation We perform as many analytic integrations as possible. In particular, we integrate over x 3 eliminating the δ-function, and we make the changes of variables In this way, the polynomial ∆ becomes linear in x 4 and x 5 , so that eq. (4.4) becomes where The integral over x 4 in eq. (4.6) is finite for d → 6 and, in this limit, we get where and Finally, we integrate over x 5 , and reduce eq. (4.8) to where (4.12)

Numerical integrations
The four remaining integration variables in eq. (4.11) are rescaled, and mapped onto a four-dimensional hypercube of unit side, so that the new variables t i have to be integrated over [0,1]. At this point, we have to consider the branch points of the dilogarithms that appear in eq. (4.11), which correspond the hypersurfaces defined by the equations It is necessary to sample carefully the integrand near these branch points. Therefore, for the integration over t 4 , we split the integration interval at the N 4 (t 1 , t 2 , t 3 ) real solutions z 4j (t 1 , t 2 , t 3 ) of eq. (4.14) which are on the interval [0, 1], Analogously, for the integration over t 3 , we split the integration interval at the N 3 (t 1 , t 2 ) real zeros z 3j (t 1 , t 2 ) of the discriminants (polynomials in (t 1 , t 2 , t 3 )) that appear in the zeros z 4j . These are the points where the hypersurfaces of eq. (4.14) are tangent to the hyperplane t 4 = constant, Analogously, for the integration over t 2 , we split the integration interval at the N 2 (t 1 ) zeros z 2j (t 1 ) of the discriminants (polynomials in (t 1 , t 2 )) that appear in the zeros z 3j , We proceed in a similar way for the last integration, To carry out the integration over a generic interval [t a , t b ], we perform the change of variables t i → u i , with 19) in order to deal with possible singularities at the endpoints. The variable u i should be integrated in (−∞, ∞) but we actually truncate the integration domain to (−M, +M ), with M suitably large (typically M ∼ 4), and we use Gauss-Legendre integration over 16 points. Note that all the singularities in the integrands are logarithmic, and therefore integrable, so we can safely set a very small value of ω, like 10 −30 . By using 16 subdivisions in each interval and in every variable we find that our integral, in the phase space point s = −1/7, t = −1/3, m 2 = 1, amounts to A similar procedure is adopted for the other integrals in Tab. 1. Case-by-case, after the analytic integrations, the corresponding integrands, in the d → 6 limit, are found to be combinations of logarithms, so that the decomposition of the integration domain, and the numerical integration can be carried out along the same lines as for the non-planar scalar box integral.

Conclusions
In this work, we presented the analytic evaluation of the two-loop master integrals needed to compute the non-planar Feynman diagrams contributing to µe elastic scattering in QED at NNLO. We adopted the same computational strategy previously applied to planar diagrams, and presented in the companion article [1]. Namely, we employed the method of differential equations and of the Magnus exponential to identify a canonical set of master integrals and we derived boundary conditions either from the regularity requirements at pseudothresholds or from the knowledge of the integrals at special kinematic points, possibly evaluated by means of auxiliary, simpler systems of differential equations. The considered master integrals were expressed as a Taylor series around four space-time dimensions, whose coefficients are written as a combination of generalised polylogarithms. We worked in the massless electron approximation, while keeping full dependence on the muon mass.
The scattering of high-energy muons on atomic electrons has been recently proposed as an ideal framework to determine, in a novel way, the leading hadronic contribution to the anomalous magnetic moment of the muon. The ambitious experimental goal of the MUonE project, namely measuring the differential cross section of the µe → µe process with an accuracy of 10ppm, requires, on the theoretical side, the knowledge of the QED corrections at NNLO. The results of the planar and non-planar master integrals we obtained imply the feasibility of the evaluation of the virtual corrections at the required order.
By crossing symmetry, our results are also relevant for muon-pair production at e + e −colliders operating well below the Z-pole, such as Belle II and VEPP-2000, as well as for the QCD corrections to heavy-quark pair production at hadron colliders. The former application is particularly interesting, as a precise knowledge of the differential cross section in QED could be exploited to constrain non-standard eeµµ interactions via the measurement of a forward-backward asymmetry. The required input integrals belong to the integral family which is identified by the set of denominators and by the external momenta A representative 4-propagator integral of this family is depicted in figure 3. IBPs reduce the integral family of eq. (A.1) to a set of 5 MIs, whose dependence on p 2 3 is parametrised in terms of the dimensionless variable The integral basis The general solution of the DEQs can be expressed in terms of harmonic polylogarithms (HPLs), i.e. GPLs with weights a i ∈ {−1, 0, 1}. The integrals I 1,2 , which are independent of x, are determined by direct integration, whereas the boundary constants of I 3,4,5 are obtained by demanding their vanishing in the regular limit x → 0. In particular, for the two triangle integrals, we obtain (A.10) These solutions are real valued in the interval 0 < x < 1. By analytic continuation to the region x < 0, we can extract the values of the two integral at p 2 3 = m 2 , which can then be used in eq. (3.17) .

B d log forms
In this appendix we collect the coefficient matrices of the d log-form for the master integrals in the non-planar integral family, defined in eqs. (2.3) and (2.5):