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

We evaluate the master integrals for the two-loop, 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 series to determine a canonical set of integrals, finally expressed as a Taylor series around four space-time dimensions, with coefficients written as 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.


Introduction
The elastic scattering of muons and electrons is one of the simplest and cleanest processes in particle physics. In spite of this simplicity, µe scattering measurements are scarse. In the 60s, experiments at CERN and Brookhaven measured this scattering cross section using accelerator-produced muons [1][2][3][4]. At the same time, µe collisions were measured by cosmic-ray experiments [5][6][7][8]. The scattering of muons off polarized electrons was then proposed as a polarimeter for high-energy muon beams in the late 80s [9] and measured by the SMC collaboration at CERN a few years later [10].
Recently, a new experiment, MUonE, has been proposed at CERN to measure 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 [11]. 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 [11,12]. In order for this new determination to be competitive with the present dispersive one, which is obtained via timelike data, 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.
In this work we take a first step towards the calculation of the full NNLO QED corrections to µe scattering. In particular, we consider the evaluation of the master integrals (MIs) occurring in the decomposition of the genuine two-loop 2 → 2 planar box-diagrams, namely all the two-loop four-point topologies for µe scattering except for the crossed double box diagram. Given the small value of the electron mass m e when compared to the muon one m, we work in the approximation m e = 0. In this case, integration-by-parts identities [34][35][36] yield the identification of a set of 65 MIs, which we compute analytically by means of the differential equation method [37][38][39]. Elaborating on recent ideas to simplify the system-solving strategy [40,41], we choose a set of MIs obeying a system of first-order differential equations (DEQs) in the kinematical variables s/m 2 and t/m 2 which is linear in the space-time dimension d, and, by means of Magnus exponential matrix [41], we derive an equivalent system of equations in canonical form [40], where the d-dependence of the associated matrices is factorized from the kinematics. Let us emphasize that the use of Magnus exponential matrix to identify a canonical basis of master integrals turned out to be very effective in the context of multi-loop integrals involving several scales [41][42][43][44]. The matrices associated with the canonical systems admit a logarithmic-differential (d log) form, whose entries are rational functions of the kinematics; therefore, the canonical MIs can be cast in a Taylor series around d = 4, with coefficients written as combinations of generalised polylogarithms (GPLs) [45][46][47][48]. The final determination of the MIs is achieved after imposing the boundary conditions, implemented by requiring the regularity of the solutions at special kinematics points, and by using simpler integrals as independent input.
The analytic expressions of the MIs have been numerically evaluated with the help of GiNaC [49] and were successfully tested against the values provided by the computer code SecDec [50]. The package Reduze [51] has been used throughout the calculations.
It is important to observe that the MIs of the QED corrections to µe → µe scattering are related by crossing to the MIs of the QCD corrections to the tt-pair production at hadron colliders. The analytic evaluation of the MIs for the leading-color corrections to pp → tt, due to planar diagrams only, was already considered in refs. [30][31][32][33]. They correspond to the MIs appearing in the evaluation of the Feynman graphs associated to the topologies T i with i ∈ {1, 2, 3, 7, 8, 9, 10} in figure 1, which we (re)compute here indepen-dently. The MIs for the planar topology T 4 and T 5 , instead, would correspond to the MIs of subleading-color contributions to tt-pair production, and were not considered previously.
For certain classes of MIs, like the ones of the processes µe → µe and pp → tt, the choice of the boundary conditions may still constitute a challenging problem. In some cases considered in refs. [30][31][32][33], the direct integration of the MIs in special kinematic configurations was addressed by using techniques based on Mellin-Barnes representations [52,53]. Alternatively, here we exploit either the regularity conditions at pseudo-thresholds or the expression of the integrals at well-behaved kinematic points. The latter might be obtained by solving simpler auxiliary systems of differential equations, hence limiting the use of direct integration only to a simple set of input integrals. Our preliminary studies make us believe that the strategy we adopt for the determination of the considered integrals is not only limited to the planar contributions, but it can be applied to the non-planar graphs as well. In particular, we show its application for the determination of the MIs for the nonplanar vertex graph [25][26][27][28][29]. Moreover, due to the similarity of the cases, we are confident that it can be very helpful for completing the analytic evaluation of the MIs needed for the two-loop QCD corrections to pp → tt, which are currently known only numerically [54][55][56][57][58].
The paper is organized as follows. In sec. 2 we describe the kinematics of µe scattering and we give a brief review of the LO and NLO QED contributions to the cross section. In sec. 3 we fix our notation and conventions for the four-point topologies. In sec. 4 we discuss the general features of the systems of differential equations satisfied by the MIs and their general solution in terms of generalised polylogarithms. In sec. 5 we describe the computation of the one-loop MIs and in sec. 6 we present the results for the planar two-loop MIs. Finally, in sec. 7 we compute the MIs for the non-planar two-loop vertex. In sec. 8 we give our conclusions. 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 expressions of the dlog-form of the matrices associated to canonical systems.
The analytic expressions of the considered MIs are given in the ancillary files accompanying the arXiv version of this publication.
The LO QED prediction for the differential cross section of the scattering in (2.1) is , (2.3) where α is the fine-structure constant. The NLO QED corrections to this cross section were computed long time ago [13][14][15][16][17][18] and revisited more recently [19]. As a first check, we recalculated these corrections and found perfect agreement with ref. [19], both for the virtual corrections and the soft photon emissions. We note that some of the pioneering publications, like [14,16], contain typos or errors, so that they cannot be directly employed.
In the rest of this paper we will work in the approximation of vanishing electron mass, m e = 0, i.e. with the kinematics specified by p 2 1 = p 2 4 = m 2 and p 2 2 = p 2 3 = 0. The master integrals will be conveniently evaluated in the non-physical region s < 0, t < 0.

Four-point topologies
The main goal of this work is the evaluation of the master integrals (MIs) of the planar twoloop four-point functions contributing to µe scattering, drawn in figure 1. For completeness, we will discuss also the evaluation of the MIs of the one-loop four-point function in figure 2.
In our conventions, the integration measure is defined as with µ being the 't Hooft scale of dimensional regularization and We choose the following set of propagators for the relevant planar four-point topologies at one-and two-loop: • For the one-loop integral family, depicted in figure 2, (3.4) • For the first two-loop integral family, which includes the topologies T 1 , T 2 , T 3 , T 7 and T 8 of figure 1, (3.5) • For the second two-loop family, which contains topologies T 4 , T 5 , T 9 and T 10 shown in figure 1, For all families, k 1 and k 2 denote the loop momenta. In the following sections, MIs will be represented by diagrams where thick lines stand for massive particles (muon), whereas thin lines stand for massless ones (electron, photon).

System of differential equations
In order to determine all MIs appearing in the three integral families defined above, we initially derive their DEQs in the dimensionless variables −s/m 2 and −t/m 2 . Upon the change of variable, the coefficients of the DEQs are rational functions of x and y. According to our system solving strategy, by means of integration-by-parts identities (IBPs), we choose an initial set of MIs F that fulfills a system of DEQs where the matrices A x ( , x, y) and A y ( , x, y) are linear in the dimensional regularization parameter = (4 − d)/2, being d the number of space-time dimensions. According to the algorithm described in [41][42][43][44], by means of Magnus exponential matrix, we identify a set of MIs I obeying canonical systems of DEQs [40], where the dependence on is factorized from the kinematics, After combining both systems of DEQs into a single total differential, we arrive at the following canonical form where the generic form of the total differential matrix for the considered MIs reads as, with M i being constant matrices. The arguments η i of this d log-form, which contain all the dependence of the DEQ on the kinematics, are referred to as the alphabet and they consist in the following 9 letters: Let us observe that, currently, there is neither a proof of existence, nor any systematic algorithm to build a basis of integrals whose system of DEQs is linear in . Nevertheless, by trial and error, we have been always able to find it within the physical contexts we have so far studied [41][42][43][44], as well as for the µe scattering. We believe it is a very important property which could be considered a prerequisite for the existence a canonical basis: in fact, a system of DEQs whose matrix is linear in can be brought into canonical form by a rotation matrix built either by means of Magnus exponential, or equivalently by means of the Wronskian matrix (formed by the solutions of the associated homogenous equations, and their derivatives), as shown for the case of systems of DEQs involving elliptic solutions [59][60][61].
The MIs presented in this paper are computed in the kinematic region where all letters are real and positive, x > 0 , 0 < y < 1 , (4.7) which corresponds to the Euclidean region s < 0, t < 0. All MIs are chosen to be finite in the → 0 limit, in such a way that I(x, y) admits a Taylor expansion in , I( , x, y) = I (0) (x, y) + I (1) (x, y) + 2 I (2) (x, y) + . . . , (4.8) with the n-th order coefficient given by where I (i) (x 0 , y 0 ) is a vector of boundary constants and ∆ (k) the weight-k operator which iterates k ordered integrations of the matrix-valued 1-form dA along a piecewisesmooth path γ in the xy-plane. Since the alphabet given in eq. (4.6) is rational and has only algebraic roots, the iterated integrals (4.10) can be directly expressed in terms of GPLs, which are defined as with w n being a vector of n arguments. The number n is referred to as the weight of G( w n ; x) and amounts to the number of iterated integrations needed to define it. Equivalently one has Each set arises from a different kinematic limit imposed on the alphabet given in eq. (4.6). We used GiNaC to numerically verify that at each order in n (up to the order n = 4), the corresponding combination of constant GPLs is proportional to Riemann ζ n . In particular, ζ n functions are known to be primitive [62][63][64], i.e. they have irreducible coproducts, ∆(ζ n ) = 1 ⊗ ζ n + ζ n ⊗ 1. Therefore, they must have vanishing coproducts components ∆ p (ζ n ) = 0 with p ∈ Π(n), where is the set of the integer partitions of n with dimensions larger than one. We explicitly verified that the combinations of GPLs of argument 1, appearing in the boundary conditions, are also primitive, although the considered coproducts components do not necessarily vanish when acting separately on each GPL involved in those combinations. Therefore, we made a simple ansatz that these combinations could be proportional to ζ n , and we checked it by means of high-precision arithmetic. Below we show some examples for these identities, For related studies see also [65][66][67].

One-loop master integrals
In this section we briefly discuss the computation of the master integrals of the one-loop four-point graph shown in figure 2, corresponding to the integral family defined in eq. (3.4). We choose the following set of MIs, which satisfy an -linear DEQ,  Figure 3: One-loop MIs T 1,...,5 .
This set of MIs satisfies a canonical DEQ of the form given in eq. (4.4), whose coefficient matrices read (in this case, M 3 and M 9 vanish), The integration of the DEQ in terms of GPLs as well as the fixing of boundary constants is straightforward. I 1,3 are obtained by direct integration and, by using the normalization of eq.(3.2), are given by The boundary constants for I 2 , I 4 and I 5 can be fixed by respectively demanding regularity at pseudothresholds s → 0, at t → 4m 2 , and at s = −t → m 2 /2. The final expression of the other MIs are, with 5 (x, y) = − 5ζ 2 + 2G(−1; x) (2G(1; y) − G(0; y)) . (5.8)

Two-loop master integrals
In this section we present the results for the planar two-loop MIs contributing to the NNLO virtual QED corrections to µe scattering, which are the main results of this work. We first discuss the computation of the MIs belonging to the integral family defined in eq. (3.5), which is associated to the topologies T 1 , T 2 , T 3 , T 7 and T 8 of figure 1, and then the MIs belonging to the integral family defined by eq. (3.6), which groups the topologies T 4 , T 5 , T 9 and T 10 .

The first integral family
For the two-loop family defined in eq. (3.5), the following set of 34 MIs fulfill an -linear system of DEQs,  Figure 4: Two-loop MIs T 1,...,34 for the first integral family.
where the T i are depicted in figure 4. Through the Magnus exponential, we rotate this set of integrals to the canonical basis This set of MIs I satisfies a system of DEQ of the form given in eq.(4.4), which can easily be integrated in terms of GPLs. The 34×34 coefficient matrices are collected in appendix B.1.
To determine the solution of the DEQ, we need to choose proper boundary values for each master integral. The boundary fixing can be achieved either by knowing the integral at some special kinematic point or by demanding the absence of unphysical thresholds that appear in the alphabet of the generic solution, defined in eq. (4.6). Below we describe in detail how the boundary constants for each integral were chosen: • The boundary values of I 1,3,4,7,9 were obtained by direct integration, • The boundary constants of I 2,8,11,23 are fixed by demanding finiteness in the limit s → 0.
• The integrals I 17 and I 18 are regular in the s → 0 limit. By imposing the regularity on their DEQ we can only fix the constant of one of them, say I 18 . The boundary constants of I 17 must be then computed in an independent way. We observe that the value of I 17 ( , 0) can be obtained in the limit p 2 1 → m 2 of a similar vertex integral with off-shell momentum p 2 1 and s ≡ (p 1 + p 2 ) 2 = p 2 2 = 0, We discuss the computation of the auxiliary vertex integral in appendix A, where we show that the limit appearing in the r.h.s of eq. (6.12) is indeed smooth and gives, • The regularity of the four-point integrals I 21,22,25,28 ... ,31 in either s → 0 or t → 4m 2 provides two boundary conditions, which can be complemented with additional relations obtained by imposing the regularity of the integrals at s = −t = m 2 /2.
• The boundary constants of integral I 24 are determined by demanding regularity in the limit s → −m 2 and t → 4m 2 .
• The boundary constants of I 34 are found by demanding finiteness in the limit u → ∞.
All results have been numerically checked with the help of the computer codes GiNaC and SecDec, and the analytic expressions of the MIs are given in electronic form in the ancillary files attached to the arXiv version of this manuscript.
• The boundary values of I 9 can be obtained by direct integration and it is given by • The boundary constants of I 19,21 can be fixed by demanding regularity when s → 0.
• Finally, the boundary constants of I 36...42 can be all determined by demanding regularity in the simultaneous limits t → 9 2 m 2 and s → −2m 2 .
All results have been numerically checked with the help of the computer codes GiNaC and SecDec, and the analytic expressions of the MIs are given in electronic form in the ancillary files attached to the arXiv version of this manuscript. Figure 5: Two-loop MIs T 1,...,42 for the second integral family.

Towards the non-planar integrals
The complete computation of the NNLO virtual QED corrections to µe scattering requires the evaluation of one last missing four-point topology, which corresponds to the non-planar diagram T 6 of figure 1. In view of future studies dedicated to this last class of integrals, we hereby show how the previously adopted strategy, based on differential equations, Magnus exponential and regularity conditions, can be efficiently applied to compute the MIs of a simpler vertex integral belonging to same family.

Master integrals for the non-planar vertex
We consider the non-planar vertex depicted in fig. 6, whose integral family is defined as where the loop propagators are chosen to be Figure 6: Non-planar two-loop three-point topology The MIs belonging to this integral family, which will be part of the full set of MIs needed for the computation of T 6 , have been already considered in the literature [25][26][27][28][29]. In all previous computations, the determination of the boundary constants resorted either to the fitting of numerical values to trascendental constants [25][26][27] or to Mellin-Barnes techniques [29]. With the present calculation, we show that they can be fixed equivalently by imposing the regularity of the solution at specific kinematic pseudo-thresholds and by matching a particular linear combination of integrals to their massless counterpart.
In order to determine the MIs belonging to the integral family (7.1), we derive their DEQ in the dimensionless variable x. We identify a set of 14 MIs which fulfills an -linear system of DEQs, Figure 7: Two-loop MIs T 1,...,14 for the non-planar vertex.
where the T i are depicted in figure 7. By making use of the Magnus exponential, we can transform these MIs into the canonical basis which satisfies a system of DEQ of the form, where M i are the constant matrices The 3 letters are real and positive in the range x > 0, which corresponds to the Euclidean region s < 0. The general solution of the system can be written in terms of one-dimensional GPLs. In order to completely determine the solution of the DEQ, we fix the boundary constants as follows: • I 1,2,3,4,5 correspond, respectively, to I 1,2,5,6,4 of the first integral family, defined in eq. (3.5).
• The regularity at s → 0 of I 8,9,10 can be used to fix the boundary constants of one single master integrals, which we choose to be I 9 .
The boundary values I 8 ( , 0) and I 10 ( , 0) can be obtained in the limit p 2 4 → m 2 of similar vertex integrals with off-shell momentum p 2 4 and s ≡ (p 1 + p 2 ) 2 = p 2 3 = 0, The computation of the auxiliary vertex integrals, which is discussed in appendix A, leads to • The boundary constant of I 11 is determined by imposing regularity when s → 0.
• I 12 corresponds to I 17 of the integral family (3.5) and the boundary constant of I 13 can be fixed by demanding regularity when s → 0.
• The boundary condition for I 14 is determined from the m → 0, or equivalently s → ∞ behaviour of the solutions, where all the internal lines of the diagrams become massless. In this regime, we search for a combination of integrals behaving as, where F ( ) is finite as z → 0.
Following the ideas outlined in [68], we begin by performing a change of variables z = 1/x = (−m 2 /s), yielding a total differential equation of the form, where M i are the constant matrices Around the z = 0 singularity, the system reduces to We perform a Jordan decomposition of M 1 , identifying the matrices J and S, related by J ≡ SM 1 S −1 , The latter can be used to define a change of the integral basis, H ≡ SI, which, by construction, obeys a system of differential equations in Jordan form, In particular, the differential equation of the truly diagonal elements H i , say i = 2, 3, 4, 6, 10, 11, 12, 13, 14 (not belonging to any block-diagonal sector), obey a trivial first-order differential equation of the form, therefore, their expression is of the type, where H i,0 is a boundary constant which may still depend on . Among the possible choices of i, we look at the element i = 11, from which we infer the behaviour around z = 0 of the following combination of canonical integrals, involving the integral I 14 . On the other side, H 11 can be computed by taking the limit z → 0 on the r.h.s. of eq. (7.18) directly at the integrand level, for the integrals I must be evaluated in the limit m → 0 (or alternatively s → ∞). To this aim, we need to pull out the prefactor m 4 coming from the integration measure defined in eq. (3.2), and to consider the definition of the canonical integrals I in terms of the linear-basis F given in eq. (7.4), In the latter equation, we took into account the vanishing of the massless tadpole in dimensional regularisation and the symmetries arising from the massless limit of the F integrals. After applying the IBPs to the massless integrals, the contributions due to all subtopologies cancel and the contribution of the massless non-planar vertex F 14 m=0 [23] is the only one left, where Therefore, Finally the boundary constant of integral I 14 can be determined by demanding the equality of eq. (7.18) and eq. (7.23).
All results have been numerically checked with the help of the computer codes GiNaC and SecDec, and the analytic expressions of the MIs are given in electronic form in the ancillary files attached to the arXiv version of this manuscript.

Conclusions
The scattering of high-energy muons on atomic electrons has been recently proposed as an ideal framework to determine the leading hadronic contribution to the anomalous magnetic moment of the muon. The ambitious experimental goal of 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. In this work, after calculating the QED corrections at NLO, which were found to be in agreement with the latest results in the literature, we investigated the feasibility of the evaluation of the corrections at NNLO. In particular, we began by considering the two-loop planar box-diagrams contributing to this process. We employed the method of differential equations and of the Magnus exponential series to identify a canonical set of master integrals. Boundary conditions were derived from the regularity requirements at pseudothresholds, or from the knowledge of the integrals at special kinematic points, evaluated by means of auxiliary, simpler systems of differential equations.
The considered master integrals were expressed as a Taylor series around four spacetime 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. Besides µe scattering, our results are relevant also for crossing-related processes such as muon-pair production at e + e − -colliders, as well as for the QCD corrections to top-pair production at hadron colliders.
The evaluation of the missing contributions due to non-planar box graphs will be the subject of a dedicated, future work -we are confident that the techniques employed here can be systematically applied for that case as well.

A Auxiliary integrals
In this appendix we briefly discuss the solution of the system of differential equations for the vertex integrals which has been used in eqs. (6.12) and (7.7) as an input for the determination of the boundary constants of some of the MIs considered in this paper.
Auxiliary vertex integral for eq. (6.12) We consider the integral family identified by the set of denominators and by external momenta p 1 , p 2 and p 3 satisfying All integrals belonging to this family can be reduced to a set of 8 MIs, whose dependence on p 2 1 is parametrized in terms of the dimensionless variable The basis of integrals fulfills a canonical system of differential equations, In the Euclidean region x > 0 the general solution of the system of differential equations can be expressed in terms of harmonic polylogarithms (HPLs) [46], and the boundary constants of all master integrals, with the only exception of I 1 = 1 and can be fixed by demanding their regularity at x → 0. In particular, for the I 7 ( , x) we obtain This expression, when it is analytically continued to the region x < 0, has a smooth limit for x → −1 ( p 2 1 = m 2 ), which has been used in eq. (6.13).
Auxiliary vertex integral for eq. (7.7) We consider the integral family identified by the set of denominators and by external momenta p 1 , p 2 and p 3 satisfying (A.14) All integrals belonging to this family can be reduced to a set of 5 MIs, whose dependence on p 2 3 is parametrized in terms of the dimensionless variable The basis of integrals fulfils a canonical system of differential equations, In the Euclidean region x > 0 the general solution of the system of differential equations can be expressed in terms of HPLs. The boundary constants I 4,5 , which are the only MIs appearing for the first time in this computation, can be fixed by demanding their regularity at x → 0. In this way, we obtain 4 (x) = 0 , I The analytic continuation of these expressions to x → −1 ( p 2 3 = m 2 ) produces the smooth limits which have been used in eq. (7.8).

B d log-forms
In this appendix we collect the coefficient matrices of the d log-forms for the master integrals in the first and second integral family, respectively defined in eqs. (3.5,3.6).

B.1 First integral family
For the first integral family, given in eq. (3.5), we have (M 3 is vanishing for this integral family):