NNLO QED contribution to the $\mu e\to \mu e$ elastic scattering

We present the current status of the Next-to-Next-to-Leading Order QED contribution to the $\mu e$-scattering. Particular focus is given to the techniques involved to tackle the virtual amplitude and their automatic implementation. Renormalization of the amplitude will be also discuss in details.


Introduction
The current measurements of the muon anomalous magnetic moment a µ = (g µ − 2)/2 indicate a discrepancy of ∼ 3.5σ from the Standard Model prediction [1]. This fact might be a hint of new physics. New upcoming experiments at Fermilab and J-PARC will improve significantly the precision of a exp µ , hence higher accuracy for a th µ from the theoretical side is required.
A very clean way to achieve the aimed precision involves the determination of the leading hadronic contribution ∆α had (q 2 ) to the electromagnetic coupling in the space-like region [2,3]. The MUonE experiment, recently proposed at the CERN, is meant to obtain an estimation of ∆α had (q 2 ) through a measurement of the differential cross section of the µe-elastic scattering. In order to be competitive with the time-like datas, a precision of 10ppm is required, both on statistical and systematic errors.
QED contributions to µe-elastic scattering constitute an irreducible background for the measurement, and they represent a relevant source of systematic errors; results for Leading Order (LO) and Next-to-Leading Order (NLO) corrections are already available [4], as well as the Nextto-Next-to-Leading Order (NNLO) hadronic ones [5,6]. However, Next-to-Next-to-Leading Order (NNLO) corrections to µe-elastic scattering are mandatory to achieve the 10ppm precision.
In this work, we summarize the current status of the double-virtual contributions to the NNLO QED µe-elastic scattering. The review is organized as follows: Section 2 recaps the external kinematic definitions we use for characterizing the µe-elastic scattering; Section 3 states the basics of the calculation process, stressing the importance of the Feynman integrals. The next sessions will present the heart of this calculation: the decomposition of the amplitude and the evaluation of the master integrals. The decomposition strategies are presented in Section 4. We discuss the main features of the Adaptive Integrand Decomposition [7,8] and its implementation Aida [9,10], and the Integral Reduction algorithms involving the Integration-by-parts identities [11][12][13][14]. These method are exploit to obtain a reduction of the full amplitude in terms of a minimal set of master integrals. Section 5 is focused on the evaluation approaches, presenting the analytical and the numerical ways. Analytical strategy is based on solving systems of differential equations [15][16][17][18] which master integrals obey. The main results of the analytical method for the µe-elastic scattering are collected into [19,20], and cross-checked numerically. Alternatively, a numerical method which can be userful in these calculation is the so-called sector decomposition [21], which will be briefly introduced.
Section 6 is devoted to present the complete Mathematica implementation which embeds each step of the calculation, from the generation of the amplitude to evaluation of the master integrals. A flowchart of the global algorithm is presented in Fig. 2.
In Section 7 we define the renormalization procedure which will be adopted to obtain a UV finite µe-elastic Scattering Amplitude. Both MS and on-shell schemes are employed to renormalize respectively the QED coupling constant and the muon fermion field and mass.
Lastly, we present the main result achieved for the µeelastic scattering and we point out the next steps needed to complete the full NNLO contribution.

µe-elastic Scattering process
The process under investigation is the elastic scattering where p i are the on-shell momenta: The physical phase-space region is constrained by the following conditions The ratio between electron and muon mass is of the order 10 −5 . This enforce the validity of the massless electron approximation, which yields many advantages to the NNLO calculation. Therefore, from now on we set m e = 0 and m µ = m.

Double-virtual contributions
Within the perturbative QFT approach, the QED cross section σ of the eµ-scattering is a series expansion in terms of the coupling constant α, For 2 → 2 scattering process, σ is also related to the Feynman amplitude M through the well-known relation where C flux is the flux factor and dPS 2 is the infinitesimal phase-space element. Matching Eq. (4) and (5), the orderby-order expression of σ is manifest, where M 0 is the Born amplitude, M 1l and M 2l are respectively the one and two-loop Feynman amplitude, M r and M 2r are the real and the double-real emission amplitudes, with respectively one and two more particles emitted in the final state, and M 1l,r the real-virtual amplitudes, one-loop amplitude with one more particle emitted in the final start. The double-virtual contribution M 2l to σ NNLO (Fig. 1) is a combination of dimensionally regularized Feynman integrals where the coefficients A¯i depend on the kinematic variables s = (s, t, m 2 ). The Feynman integrals 1 I i (s) have the general form N¯i(k, s) and depend on the kinematic invariants s and the dimension d. The numerator N¯i is a general polynomial of the scalar products between loop and external momenta. The denominators D j are inverse scalar propagators, whose form is D j = q 2 j (k, p)−m 2 j . The multi-indexī = (i 1 , · · · , i n ) contains the exponents of the denominators.
Direct integration of Feynman integrals objects is usually not allowed. A viable way to calculate them efficiently involves decomposition methods: Feynman integrals can be expressed into combinations of "simpler" integrals. Such integrals can be evaluated either analytically of numerically. In the next session, the algorithm employed in the calculation of the amplitude (7) is presented.

Adaptive Integrand Decomposition
As extensively discussed in [7,8], the Adaptive Integrand Decomposition (AID) is an algorithm that works at integrand level. It exploits the (rational and) polynomial structure of the integrand in order to achieve an iterative decomposition of the numerator.
The basic ingredient is the splitting of the ddimensional space into parallel and transverse component [28] w.r.t. the space spanned by the independent external momenta: It can be shown that Eq. (9) makes the polynomial structure of the integrand manifest, and liable to perform polynomial division without passing by the Gröbner basis [29]. Therefore, the numerator can be written as (10) The polynomial division can be iterated for every N¯i −ē j (k, s). After the complete division of the numerator, the integrand is cast into the following form The residues ∆¯j(k, s) are polynomials depending of the irreducible scalar products, a minimal set of scalar product between loop and external momenta. Therefore, the total double-virtual amplitude can be cast into the form AID has been successfully applied to many amplitude [7,8,10]. An automated Mathematica implementation of this algorithm has been developed, called Aida [9,10]. It accepts the FeynArts [30] output (namely the complete unreduced amplitude) and perform the AID, exploiting the finite fields reconstruction made available through Finite-Flow [31]. Its output is finally translate into integral notation.

Integration-by-parts identities
As a consequence of the d-dimensionality and the shift/rotational invariance of the Feynman integrals, an entire class of new relations [11] can be found: where v µ can be either a loop or an external momenta. These identities between integrals are called Integration-By-Parts Identities (IBPs) [11,12,32]. Letting f (k) be the reduced integrands coming from AID, IBPs are employed to generate a large system of equalities, which are not all independent. As a consequence, there exists a minimal set of integrals which can not be reduced further: they are known as Master Integrals (MIs) [33]. The general expression for a MI is where S 1 . . . S m are the irreducible scalar products. IBPs allow the ultimate decomposition of the Feynman integrals in the basis of MIs, k ∆¯j(k, s) which amounts the complete double-virtual contribution in this simple form, The choice of the MIs is not unique. Laporta algorithm [12] is an automatable way to choose MIs which respect a EPJ Web of Conferences "simplicity" criterion. Although Laporta MIs look simpler, they might be not the convenient choice if one wants to evaluate them analytically [34][35][36].
There exist several public codes [37][38][39] which can provide IBPs and the Laporta MIs for any input topology, such as reduze [13] and kira [14]. It is important to underline that these codes may accept as input a custom MIs basis; in this case, their output will be the IBPs expressed in terms of the chosen basis.

Evaluating the Master Integrals
The -expansion of MIs can be achieved by either analytical or numerical evaluation. In this picture, d = 4 − 2 and the expansion reads where n p is the order of the higher order pole into the T i (s) series. In practice, the series expansion will be truncated to a finite order, such that every Feynman amplitude is expanded up to O( ). Here, the evaluation strategies employed in this calculation are presented.

Differential equations
The analytical approach exploits the IBPs to build a system of differential equation [15][16][17][18] (see also for review [40,41]): where the MIs basis is chosen in such a way that the matrix A( , s j ) depends linearly on . If there exist a rotation matrix R(s j ) such that Eq. (19) can be cast into the canonical form [41], simbolically expressed as the PDE admits a solution in terms of the Chen's iterated integral and its -expansion casts F i (s j ) to the form (18). The boundary condition F k (s j0 ) are chosen by exploiting regularity conditions into the phase space appropriate kinematic limits. The matrix R(s j ) can be found by means of the Magnus exponential method [36,42,43], that provides a close formula for the rotation matrix, s j0 ), , The complete order-by-order expression for Ω(s j , s j0 ) can be found in [36].

Sector decomposition
The alternative way to obtain (18) involves numerical methods.
A general algorithm valid for any multi-loop Feynman integral is the sector decomposition [21]. Briefly, Feynman parametrization applied T i (s) changes the integration space from the Minkowski (or Euclidean) space to the ndimensional unit cube 2 The unit cube C n can be iteratively decomposed and deformed into multiple sub-domains, until the pole structure of the integral manifests multiplicatively into the integrand, and subsequently isolated by means of the end-point subtraction method. For a j = b j = 1, such method reads Once the sector decomposition has been applied, MIs lie to the following form The integrals appearing in Eq. (27) are finite and they can be integrated using numerical algorithms. The integration on the physical phase-space region usually requires contour deformation procedures, such that the threshold singularities of the integrand can be avoided and the numerical stability is preserved.
The complete algorithm has been implemented into a code, SecDec [52]. It goes through the interation of the Sector Decomposition algorithm and performs the numerical evaluation using Monte Carlo integration methods available into the Cuba libraries [53], with the possibility of deform the integration contour.

Automation
A complete automation of the evaluating algorithm presented in Sections 4. and 5. has been developed. Due to the symbolic structure of the amplitude, Mathematica has been chosen to be the main workstation. The specific codes performing the decomposition and the evaluation of the amplitude have been embedded into a Mathematica script.
The generation of the double-virtual contribution is carried out by FeynArts and FeynCalc, which provide the raw input for the decomposition algorithms. In particular, FeynCalc performs the Dirac and tensor algebra and deals with the eventual Dirac traces.
The amplitude is now ready to be decomposed at integrand level by means of Aida, obtaining the structure expressed in (16). The new integrands are then converted into integrals, by a convenient notation.
Integrals belonging to the amplitude are given as input to Reduze or Kira, which perform the IBP reduction. An interface has been build, that automates the configuration of Reduze and reads its output into a Mathematica database file. At this stage, the double-virtual amplitude takes the form (17).
The last step is the evaluation of the MIs. Analytical values of the MIs is stored into Mathematica package [19,20]. It automatically replaces the symbolic integrals to its -pole Laurent expansion. Alternatively, an additional interface with SecDec has been developed. It automatically provides its input files and automate the numerical evaluation. The output of SecDec is then stored into a Mathematica package. To check the consistency of the result, both of the approached have been used.
The algorithm here presented is actually completely general, and this calculation represents a strong check of the validity of the method. A flowchart representing the data flow is given in Fig. 2.

Renormalization
QED is a renormalizable theory. As a consequence, UV divergencies arising from divergent loop integrals can be regularized and canceled order-by-order in the coupling expansion.
A diagrammatic approach to the renormalization is the most convenient way to obtain a UV finite amplitude out of the NNLO scattering amplitude; renormalized perturbation theory will be the ideal framework.
In order to renormalize QED, we consider the QED Lagrangian with N f = 1 + 1 active flavor, one massive (muon) and one massless (electron) [54]: In this approach, fields, couplings and masses are treated to be bare (denoted with the subscript 0). Bare quantities The bare Lagrangian can be expressed in terms of the renormalized quantites: where and The counterterm Lagrangian L ct contains the divergent behaviour arising from the bare quantities.
Ward Identities in QED constrain ratios between the normalizations Z i [54].
where κ ξ and κ e are arbitrary finite constants. By choosing κ i = 1, the normalizations become Note that the choice of the renormalization scheme for the photon field A µ leads automatically to fix the same prescription to the coupling e. This is a direct consequence of the Ward Identities (34) in QED .
It is important to stress that Eq. (33) constrains only the divergent part of the normalizations, letting their finite part unconstrained. Hence, two different renormalization prescription for Z e and Z A can be chosen, e.g. Z A in onshell scheme and Z e in MS scheme. This alternative choice brings the Ward identity to be which yields no contradiction with Eq. (33). Imposing the Ward identities (33), the counterterms Lagrangian becomes Finally, by setting The Ward identities have been chosen in such a way that only field and mass renormalization is needed. The coupling constant e is automatically renormalized thanks to the Ward identities.
The counterterm Lagrangian provides additional Feynman rules, which lead to UV finite amplitude. Physical results can only be obtained by choosing a renormalization scheme that fixes the residual unphysical degrees of freedom.
The choice of the on-shell renormalization scheme for the muon mass and fields sets δR ψ,i = 0 identically. Therefore, the renormalization scheme choice applied in the µe → µe scattering is: -MS scheme for the coupling, the photon and the electron fields -On-shell scheme for the muon field and mass The muon field is treated differently due to its nonvanishing mass. Despite of the complications it introduces at integration level, the opportunity to use on-shell scheme for the muon yields simplification at renormalization level.

Summary
There are 69 two-loop Feynman diagrams contributing to the QED NNLO µe-scattering process (Fig. 1). The complete amplitude gets contributions from ∼ 10 4 Feynman integrals with maximum rank equal to 4. AID and IBPs reduction with respectively Aida and Kira have decreased significantly the number of integrals [19,20], providing N MI = 120 with max rank equal to 2 (Fig. 3). These MIs are the irreducible set of integrals that one has to evaluate and further reduction of the rank is not possible. The presence of irreducible scalar products, typical of two-loop topology, does not allow a reduction to scalar (rank-0) integrals.