On μe-scattering at NNLO in QED

We report on the current status of the analytic evaluation of the two-loop corrections to the μescattering in Quantum Electrodynamics, presenting state-of-the art techniques which have been developed to address this challenging task.


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 acceleratorproduced 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 [34], we took a first step towards the calculation of the full NNLO QED corrections to µe scattering. In particular, we considered 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 worked in the approximation m e = 0. In this case, integration-by-parts identities [35][36][37] yielded the identification of a set of 65 MIs, which we computed analytically by means of the differential equation method [38][39][40]. Elaborating on recent ideas to simplify the systemsolving strategy [41,42], we chose 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 [42], we derived an equivalent system of equations in canonical form [41], 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 [42][43][44][45]. We found that the matrices associated with the canonical systems admit a logarithmic-differential (dlog) 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) [46][47][48][49]. The final determination of the MIs was 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 were numerically evaluated with the help of GiNaC [50] and were successfully tested against the values provided by the computer code SecDec [51]. The package Reduze [52] was 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)computed in [34], independently. 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 [53,54]. Alternatively, our approach exploited either the regularity conditions at pseudothresholds or the expression of the integrals at wellbehaved kinematic points. The latter were 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 adopted 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, in [34], we showed its application for the determination of the MIs for the non-planar 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 [55][56][57][58][59].
The evaluation of the MIs is only a first step towards the complete evaluation of the two-loop amplitudes. By means of the adaptive integrand decomposition [60][61][62] and the integration-by-parts identities, we can decompose the whole amplitude in terms of MIs. To achieve this task, we have been developing a general framework for the automatic evaluation of two-loop amplitudes, called Aida (Adaptive Integrand Decomposition Algorithm), and we present its first application to the case of µe-scattering.
In the following, we report on our findings.
The LO QED prediction for the differential cross section of the scattering in (1) is 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
We focus on the evaluation of the master integrals (MIs) of the planar two-loop 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.
We consider -loop m-denominator Feynman integrals In our conventions, the integration measure is defined as with µ being the 't Hooft scale of dimensional regularization and We display the relevant planar four-point topologies at one-and two-loop in families: • the one-loop integral family, depicted in figure 2; • the first two-loop integral family, which includes the topologies T 1 , T 2 , T 3 , T 7 and T 8 of figure 1; • 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 identify 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 [42][43][44][45], by means of Magnus exponential matrix, we identify a set of MIs I obeying canonical systems of DEQs [41], 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 dlog-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: The MIs presented in this paper are computed in the kinematic region where all letters are real and positive, 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 , 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 piecewise-smooth path γ in the xyplane. Since the alphabet given in eq. (12) is rational and has only algebraic roots, the iterated integrals (16) 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 GPLs fulfill shuffle algebra relations of the form where the shuffle product m n denotes all possible merges of m and n while preserving their respective orderings.
The analytic continuation of the MIs to the physical region defined in sec. 2 can be obtained through by-now standard techniques.

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. We choose the following set of MIs, which satisfy an -linear DEQ,  Figure 3: One-loop MIs T 1,...,5 .
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   (0; y)) .
(5.8) This set of MIs satisfies a canonical DEQ of the form given in eq. (10). 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.(5), 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  Figure 4: Two-loop MIs T 1,...,34 for the first integral family.

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 -10 -

Two-loop master integrals
In this section we present the results for the planar twoloop MIs contributing to the NNLO virtual QED corrections to µe scattering.

The first integral family
For the two-loop first integral family associated to the topologies T 1 , T 2 , T 3 , T 7 and T 8 of figure 1, the follow-ing set of 34 MIs fulfill an -linear system of DEQs, F 10 = 3 T 10 , F 11 = 3 T 11 , F 12 = 3 T 12 , where the T i are depicted in figure 4. Through the Magnus exponential, we rotate this set of integrals to the canonical basis I 1 = F 1 , This set of MIs I satisfies a system of DEQ of the form given in eq. (10), which can easily be integrated in terms of GPLs.
To determine the solution of the DEQ for the MIs of both two-loop families, we choose proper boundary val-ues 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. (12). For more details, we forward the reader to Ref. [34].
All results have been numerically checked with the help of the computer codes GiNaC and SecDec.

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. We are confident that 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.

Adaptive Integrand Decomposition
The decomposition of multiloop scattering amplitudes in terms of independent functions, together with the subsequent determination of the latter, is a viable alternative to the direct integration which, for non-trivial processes, may require the calculation of a prohibitively large number of complicated Feynman integrals.
Decomposing multi-loop amplitudes in terms of independent integrals can become problematic when the number of the scales of the diagrams increases, due to the exchange or to the production of massive particles, or when a large number of external particles are scattered, or when the morphology of the contributing diagrams becomes involved. The integrand decomposition algorithm has the advantage of treating scattering amplitudes involving massive particles at the same price of amplitudes for massless scattering. The output of the reduction procedure is the partial fractioning of the original integrand, namely the determination of the remainders of the successive division between the numerator and (the partitions of the product of) the denominators. Upon integration, the partial fraction formula correspond to rewrite the original amplitudes as a combination of independent integrals. However, the result of the integrand decomposition represents an intermediate step towards the complete amplitude reduction. In fact, additional relations among those integrals, like integration-by-parts identities, can minimise the number of independent master integrals (MIs) which can appear in the final formulas.
The integrand decomposition algorithm [63][64][65][66][67][68][69] played a key role for the automation of one-loop corrections to high-multiplicity scattering processes. The extension of this approach at two-loop and beyond [70][71][72][73] has been under intense investigation. The recent developments on the integrand side have been accompanied by important developments for novel derivation of the integral relations needed to identify MIs [74][75][76][77], as well as by progress in the ability of computing the latter analytically [41,42,78] as well as numerically [51,79]. This vivid research has been largely due to the deeper understanding of the properties of the integrands of Feynman graph, and of the refined algebraic and differential calculus which control them.
In this section, we recall the main features of the Adaptive Integrand Decomposition (AID) [60][61][62], implemented in the package Aida, and show its first applications to the µe-scattering.
Witin AID, the space-time dimension, d = 4 − 2 , are decomposed diagram-by-diagram, as d = d + d ⊥ , into a space spanned by the external momenta flowing in the leg -with dimensions d -referred to as parallel (or longitudinal) space, and a space spanned by the external momenta flowing in the leg -with dimensions d ⊥ -referred to as orthogonal (or transverse) space. The latter is formed by the union of the four-dimensional complement of the longitudinal space, and of the extra dimensional −2 -space.
In the structure of the Feynman integrals, d-dimensional loop momenta are defined as with x ji e α j , j=d +1 x ji e α j + µ α i , l=d +1 x li x l j + µ i j .
In Eq. (33), l i is a vector of the d -dimensional space spanned by the external momenta, and λ i belongs the d ⊥dimensional orthogonal subspace. Moreover e α i i=1,...,4 is a 4-dimensional basis, while µ α i lie in the −2 -space. Within this parametrisation, denominators appear to depend on less variables than the numerators, yielding a simplification of the decomposition procedure.
Let us indicate with z the full set of ( +9)/2 variables where x i (x ⊥ i ) are the components of the loop momenta parallel (orthogonal) to the external kinematics, the denominators are reduced to polynomials in the subset of variables so that the general r-point integrand has the form Since numerator and denominators depend on different variables, the adaptive integrand decomposition can proceed along the following algorithm: 1. Divide: we divide the numerator N i 1 ...i r (τ, x ⊥ ) modulo the Gröbner basis G i 1 ···i r (τ) of the ideal J i 1 ···i r (τ) generated by the set of denominators. The polynomial division is performed by adopting the lexicographic ordering λ i j x , The Gröbner basis does not need to be explicitly computed, since, with the choice of variables and the ordering described here, the division is equivalent to applying the set of linear relations described above.

2.
Integrate: Since denominators do not depend on transverse variables, x ⊥ , we can integrate the residue ∆ i 1 ...i r over transverse directions. This integration is carried out by expressing ∆ i 1 ...i r in terms of Gegenbauer polynomials, i.e., Where ∆ int i 1 ...i r is a polynomial in τ whose coefficients depend on the space-time dimension d.

3.
Divide: the structure of the integrated residue suggests a second division. This can be seen from the dependence ∆ int i 1 ...i r has on the variables τ. In fact, after applying the division, similarly as in the first step of this algorithm, we get where the new residue ∆ i 1 ...i r (x ) can only depend on x .
At the end of the decomposition, the integrand is written as where the residue functions ∆ (k j ) are polynomials in the irreducible scalar products, whose coefficients may depend on the external kinematics as well as on d (after the first division, the polynomial ∆(k j ) have no ddependence). This algorithm has been succesfully applied to a few cases as shown in in [60], and also to the leading color contribution to the two-loop all-plus five-gluon amplitude [61,62,80].
The Mathematica package Aida takes as input the integrands of Feynman diagrams generated by FeynArts [81] and FeynCalc [82,83] and applies the AID algorithm to them. The output of Aida is the expression of the original amplitude, written as sum of Feynamn diagrams, in terms of integrals, which can be further reduced in terms of MIs by automatic software packages like Reduze.
For the µe-scattering, we succesfully applied Aida to the one-loop diagrams shown in fig.6 and to the two-loop planar and non-planar diagrams in fig.7. We considered also the automatic generation of the one-loop renormalisation countertem-diagrams, which can be processed together with the two-loop diagrams, to provide a result where the one-loop UV (sub)divergencies have been removed.

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 proceedings, we reported our investigations on the feasibility of the evaluation of the corrections at NNLO. In particular, we began by considering the twoloop 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 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. 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 nonplanar 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.
At the same time, we considered the generation of the two-loop scattering amplitudes within the novel software Aida, implementing the adaptive integrand decomposition algorithm, and the decomposition of its output in terms of a minimal number of master integrals, by interfacing Aida to existing softwares carrying out the integration-by-parts algebra.
Our preliminary studies show that the evaluation of the virtual NNLO corrections are feasable. We shall report on that in future communications.