Hadronic corrections to $\mu$-$e$ scattering at NNLO with space-like data

The Standard Model prediction for $\mu$-$e$ scattering at Next-to-Next-to-Leading Order (NNLO) contains non-perturbative QCD contributions given by diagrams with a hadronic vacuum polarization insertion in the photon propagator. By taking advantage of the hyperspherical integration method, we show that the subset of hadronic NNLO corrections where the vacuum polarization appears inside a loop, the irreducible diagrams, can be calculated employing the hadronic vacuum polarization in the space-like region, without making use of the $R$ ratio and time-like data. We present the analytic expressions of the kernels necessary to evaluate numerically the two types of irreducible diagrams: the two-loop vertex and box corrections. As a cross check, we evaluate these corrections numerically and we compare them with the results given by the traditional dispersive approach and with analytic two-loop vertex results in QED.


Introduction
The goal of the new g−2 experiments at Fermilab in the United States and at J-PARC in Japan is to measure the muon anomalous magnetic moment with a precision of 1.6 × 10 −10 [1,2] -corresponding to 140 ppb -an improvement by a factor of four of the final BNL E821 experiment's uncertainty: δa exp µ = 6.3 × 10 −10 (540 ppb) [3]. By all means, the theoretical prediction must keep up with the experimental precision. At present, the Standard Model prediction of the muon g−2 is limited by the uncertainty of the Hadronic Leading Order (HLO) and Light-By-Light (HLBL) contributions, that cannot be computed in perturbative QCD.
Recently a new experiment, MUonE , has been proposed at CERN to determine a HLO µ by measuring the running fine-structure constant α, α(q 2 ) = α 1 + Re Π(q 2 ) , (1.2) in the space-like region (q 2 < 0) in µ-e scattering as a function of the squared momentum transfer t [22,23]. The function Π(q 2 ) is the renormalized photon vacuum polarization from which it is possible to extract the hadronic contribution Π had (q 2 ) by subtracting from Π(q 2 ) the leptonic part Π lep (q 2 ), calculated in perturbative QED up to four loops [18]. 1 The hadronic vacuum polarization Π had (q 2 ) in the space-like region can provide an independent determination of a HLO µ thanks to the formula [24] a HLO To determine a HLO µ with an error of about 2×10 −10 the MUonE experiment must measure the differential cross section with statistic and systematic uncertainties of the order of 10 ppm.
To this end, a Monte Carlo event generator that includes all relevant corrections needed to reach such level of precision must be developed. It must contain QED and QCD radiative corrections up to Next-to-Next-to-Leading-Order (NNLO) in α matched to leadinglogarithmic corrections resummed to all orders. A Next-to-Leading-Order (NLO) Monte Carlo generator based on the existing BabaYaga [25][26][27][28][29] framework is currently under development [30]. A first step towards the evaluation of NNLO corrections was presented in [31,32] where the QED two-loop master integrals were calculated for finite muon mass and vanishing electron mass.
In this paper we will focus on the hadronic contributions to the µ-e scattering cross section. These corrections are genuinely non-perturbative and cannot be calculated in perturbative QCD since the scattering process will take place at a center-of-mass energy of about 0.5 GeV. The hadronic contribution to the NLO cross section -order α 3 -comes from the diagram in figure 1a; it corresponds to the leading effect of the fine-structureconstant running in an expansion of α = α(0).
The hadronic corrections to the NNLO cross sections -order α 4 -can be divided into four classes of diagrams.
I. Tree-level diagrams with double vacuum polarization insertion (figure 1b), either two hadronic insertions or one Π had and one Π lep . These are the second order effects of the running of α.
II. QED one-loop diagrams in combination with one insertion of Π had in the t-channel photon (figure 1c). Their contribution to the cross section is linear in Π had (t) and can be obtained directly from the QED one-loop amplitudes.
III. Real photon emission with a dressed photon propagator in the t channel (figure 1d).
All the diagrams in class I-III are factorizable or reducible since they are given by the product of a QED amplitude times the function Π had (q 2 ) evaluated at q 2 = t. A fourth class of non-factorizable or irreducible diagrams must be also considered: IV. One-loop QED amplitudes with a dressed photon propagator inserted inside the loop. They can be further subdivided into vertex and box corrections (figure 1e).
Note that there is no LBL contribution to the cross section up to N 3 LO -order α 5 . 2 Moreover, we remind the reader that the analysis of future MUonE data will also require the study of µ-e scattering processes with final states containing hadrons. Final states of Bhabha scattering containing hadrons were studied in [33]. The traditional approach to calculate the amplitudes in class IV uses the dispersion relation, to replace the dressed photon propagator inside the loop -now q stands for the loop momentum -with the r.h.s. of eq. (1.4) where the momentum q appears only in the term 1/(q 2 − z). This allows us to interchange the integration order and calculate as a first step the one-loop integrals with the dressed photon propagator replaced by a massive 2 In figure 1a a virtual photon can be emitted and reabsorbed by the hadronic bubble. In the spirit of the common nomenclature, we do not consider this kind of two-loop diagrams as a part of the hadronic NNLO corrections. This effect is commonly included in R(s) as final state radiation, so no additional contribution has to be taken into account.
gauge boson of mass √ z. Later on the z-dependent scattering amplitudes are convoluted with the R ratio. The dispersive approach was employed for instance to calculate the hadronic corrections to Bhabha scattering [34,35]. A complete calculation of the hadronic corrections to µ-e scattering at NNLO with the dispersive approach will be presented soon [36].
The dispersive approach requires the R ratio as an input. Therefore the MUonE's determination of Π had (q 2 ) and a HLO µ from space-like data would still depend marginally on time-like data if the dispersive approach were employed in the evaluation of the hadronic NNLO corrections. The alternative is to use the very same space-like data measured by MUonE to calculate the hadronic NNLO corrections iteratively without dispersion relation and the R ratio. One could approximate the function Π had (t) through successive iterations: as a first step the hadronic NNLO corrections can be switched off in the Monte Carlo and a first approximation for Π had (t) extracted. Afterwards, the Monte Carlo can be supplied with this first approximation to evaluate the hadronic NNLO corrections, a second approximation calculated and the process further iterated. The factorizable diagrams in class I, II and III depend on Π had (t), so they are well suited to implement this iterative procedure. What about the contributions of the irreducible diagrams in class IV?
In this paper we will show that also the non-factorizable diagrams can be calculated using the hadronic vacuum polarization in the space-like region by making use of the hyperspherical integration method. This method was exploited for instance to evaluate parts of the QED three-loop contributions to the electron g−2 [37][38][39][40][41][42], and more recently to calculate the pion pole contribution to a HLBL µ [43,44] and in the dispersive approach to the HLBL [45,46].
The loop integrals containing Π had (q 2 ) can be calculated as follows: after analytic continuation of internal and external momenta to the Euclidean region, one introduces spherical coordinates for the loop momentum q. The angular dependence of the Feynman propagators can be made explicit by an expansion in Gegenbauer polynomials. Afterwards, the integration with respect to the angular variables can be carried out analytically by taking advantage of the orthogonality properties of these polynomials. In this way, the non-factorizable diagrams are left in the form of a residual radial integration, which is calculated numerically once provided with the hadronic vacuum polarization in the space-like region. One must note, however, that the integral (1.5) requires the knowledge of Π had (q 2 ) for any q 2 < 0, while Π had is experimentally accessible only in a finite range of t. Therefore the proposed iterative procedure will require in any case an extrapolation between the measured region and the high-energy tail close to infinity. Lattice data could also come to the aid in the intermediate region. Padé approximants could be used in this merging procedure. They have been employed, for example, to evaluate the QED vacuum polarization function at four loops and its contribution to the g−2 at five loops [47] knowing the first terms of Π lep (q 2 ) in an expansion around q 2 = 0, 4m 2 , +∞. It is beyond the scope of this paper to study the impact of such extrapolation and the error that it would introduce in the evaluation of a HLO µ . The paper is organized as follows. In section 2 we will review the hyperspherical integration method. The master formulae for the evaluation of the irreducible diagrams will be presented in section 3 for the vertex corrections and section 4 for the boxes. In section 5 we will make a comparison between the traditional dispersive method and the hyperspherical method. Conclusion are drawn in section 6. The appendix contains an example of a one-loop calculation with the hyperspherical method.

The Hyperspherical Integration Method
In this section we will give a short review of the hyperspherical integration method. Each of the diagrams in class IV contains an insertion of the SM vacuum polarization tensor with four momentum q, is the electromagnetic current and the sum runs over fermions with charges Q f . The weak interactions will be ignored. Throughout this paper we will always assume Π(q 2 ) to be the renormalized vacuum polarization fulfilling Π(q 2 = 0) = 0. In each loop diagram we choose the routing of the loop momentum q in such a way that the momentum flowing through the dressed photon propagator is exactly q, so that the loop integral has the following form: where n = 2 (3) for the vertex (the box) corrections, D i = (q + k i ) − m 2 i + iε are propagator denominators and k i are linear combination of the external momenta p 1 , · · · , p n . The numerator N (q, p 1 , . . . , p n ) is assumed to be a scalar function. We will work in D = 4 dimensions. This choice is dictated mainly from the fact that hyperspherical integration of three Feynman propagators, necessary for the boxes, are known only in four-dimensions [48]. More details will be given further on.
We begin with the analytic continuation of all external momenta into the Euclidean region and with a Wick rotation of the integration contour. Do we need at this point to add the propagator pole residues after the Wick rotation? When we compute a one-loop integral in the traditional way, after introducing Feynman parameters and shifting the loop momentum one is left with a denominator of the form 1/(q 2 − ∆ + iε) n , that has poles at q 0 = ± q 2 + ∆ ∓ iε. Since the poles are in the top-left and bottom-right quadrants of the complex q 0 plane, integrating over the real axis is equivalent to integrating over the imaginary axis. However in the hyperspherical approach, we cannot shift the loop momentum and therefore we are left at the denominator with the product of propagators of the form 1/[(q − p) 2 − m 2 ], that has two poles in the q 0 complex plane at The two poles are not centered at the origin. If p 0 is sufficiently large, q − 0 lies in the top-right quadrant of the q 0 plane. Therefore the integration over the real axis is different from the integration along the imaginary axis: the residue of the pole q − 0 must be taken into account also. Phrased differently, integration over the real axis (in blue in figure 2a) is equivalent to the red path in figure 2a, which proceeds along the imaginary axis but avoid the q − 0 pole turning around it. Analytic continuation of the external momenta into the Euclidean region moves then the location of these poles. For example we can let the energy p 0 in (2.3) to acquire a phase e iφ which is then varied from 0 to π/2. In this way, the pole q − 0 moves to the top-left quadrant, while the pole q + 0 to the top-right one. No pole should cross the integration contours. Therefore the blue path along the real axis must be deformed because the pole q + 0 moves to the upper side, as shown in figure 2b, while the path along the imaginary axis becomes straight, in red in figure 2b, since the pole q − 0 moves to the left. So after Wick rotation no pole residue must be included if both internal and external momenta become Euclidean.
Note in addition that Π had (q 2 ) has a branch point at the pion threshold q 2 = 4m 2 π . In the q 0 complex plane it corresponds to two branch points at q 0 = ± q 2 + 4m 2 π ∓ iε. They are unaffected by the analytic continuation because their position is independent on the external momenta. So they do not interfere with the Wick rotation. Furthermore, the vacuum polarization does not introduce any other isolated singularity since its poles in the q 2 complex plane, corresponding to unstable resonances, are hidden below the real axis in the unphysical sheet. Figure 2: Integration contour in the complex plane q 0 before (2a) and after (2b) analytic continuation of the external momenta into the Euclidean region. The Feynman propagator poles q ± 0 move in the complex plane and the blue and red integration paths must be deformed accordingly. Now four-dimensional hyperspherical coordinates can be introduced for the loop momentum q: We will denote with capital and lowercase letters the momenta in the Euclidean and in the Minkowski space, respectively. Note that after Wick rotation the vacuum polarization's argument becomes negative: Π had (−Q 2 < 0). The angular dependence of the propagators D i in (2.2) can be made explicit by the expansion whereQ andP are unit vectors along the direction of Q and P , n (x) are an orthogonal basis of functions over the interval [−1, 1] with respect to the weight function √ 1 − x 2 . This allows us to perform the angular integration using the orthogonality conditions Since after the angular integration the momenta are still space-like, i.e. p 2 i < 0, we need eventually to analytically continue back the results to the time-like region. An example of a one-loop integral calculation with the hyperspherical method is presented in the appendix, where its analytic continuation is also further discussed.
There is a caveat however: the integral of the product of three Gegenbauer polynomials evaluated at three differentQ ·P i is unknown and therefore the angular integration cannot be performed with the method described above. This occurs when we calculate the box diagrams: the product of three denominators -the fourth does not depend on the angles -would become the product of three Gegenbauer polynomials if the expansion (2.4) were employed. The angular integrals can be evaluated nevertheless by brute force integrating directly with respect to the three hyperspherical angles, avoiding the expansion (2.4). The general solution of an integral with three denominators in D = 4 was given long time ago by Laporta in a not very-well-known article [48]. This is the main reason why our calculation is carried out in D = 4 and not in dimensional regularization, even if the hyperspherical method can be applied to D dimensions as well (see e.g. [49]). After the hyperspherical integration, the residual Q 2 radial integral can be ill-defined because of the bad behaviour of the integrand at infinity -an ultraviolet (UV) divergence -or at some finite value of Q 2 , an infrared (IR) one. Vertex corrections are UV divergent but IR finite because Π had (0) = 0 regularizes the behaviour of the kernel at the origin. Thanks to the on-shell renormalization prescription, the radial integral can be regularized by calculating the vertex together with its counter-term, which is the vertex itself in the limit of zero momentum transfer t → 0. The counter term built in this way cancels the original UV divergence pointwise in momentum space [38,40].
Vice versa, the boxes contain soft IR divergences but are UV finite. A small photon mass λ can be introduced to regularize the integral. It is even possible to avoid an explicit calculation of IR divergent integrals containing Π had (q 2 ) by observing that in the soft limit the box diagrams are proportional to the "tree-level" amplitude, i.e. the Born amplitude with a dressed photon propagator in figure 1a. Indeed the soft pole arises when the momentum of the undressed photon goes to zero and the momentum of the dressed one is t. This suggests us the possibility to extract the IR poles with the following subtraction: The first integral on the r.h.s. of eq. (2.9) is now free of IR divergences and can be evaluated setting λ = 0. The soft pole appears only in the second integral of (2.9) that does not contain anymore Π had and can be calculated analytically with standard methods. The last technical ingredient to discuss is how to perform the angular integration when the loop momentum q appears also at the numerator. One occurrence of the scalar product k i · q can be always removed against one of the propagators D i by writing i . Additional k i · q in the numerator can be further simplified using the technique described in the appendix of ref. [40]. The one-denominator case is straightforward: given that C 1 (x) = 2x. The angular integral is performed by expanding the denominator and using the orthogonality condition (2.6). For the two-propagator case we write The term I µ must be a linear combination of the Euclidean vectors K j and K k that appear at the denominator with scalar coefficients. Introducing two orthonormal vectorsê µ 1 and e µ 2 in the two-dimensional space spanned by K j and K k , we can replace the loop momentum Q µ in (2.11) with its projection onto the space spanned by K j and K k : Now the integral I µ contains, through the scalar productê i · Q, terms like K j · Q/D j or K k · Q/D k , which can be simplified as before, leading to an integrand without Q in the numerator. Figure 3: The leading contribution of the hadronic vacuum polarization to the QED vertex.

The Vertex Corrections
Having introduced the hyperspherical method, we can now apply it to calculate the hadronic vacuum polarization contribution to the QED form factors which can be used in a second stage to calculate the irreducible vertex corrections in class IV. The 1PI amplitude Γ µ (k) describing the interaction between a photon and the initial and final states of an on-shell lepton , with four-momenta p 1 and p 2 , respectively, can be written in terms of the Dirac and Pauli form factors F 1 and F 2 : , m is the lepton mass and k = p 2 − p 1 is the incoming fourmomentum of the off-shell photon. Let us call F had 1 (k 2 ) and F had 2 (k 2 ) the leading contribution of the hadronic vacuum polarization to the form factors F 1 and F 2 given by the two-loop diagram in figure 3. The vertex corrections in class IV can be expressed in term of F had 1 (k 2 ) and F had 2 (k 2 ), with = e, µ. Since there the role of the off-shell photon with momentum k is played by the photon in t-channel, we will identify k 2 with the Mandelstam variable t < 0, and we will restrict the calculation to the region t < 0.
The form factors F had 1 (k 2 ) and F had 2 (k 2 ) are extracted from the amplitude with the projector technique [50] and the loop integral calculated with the hyperspherical method. The final expression for the form factors can be cast in the form where i = 1, 2 and x is related to the radial variable Q 2 by The angular integration gives a well-behaved integrand f 2 , while we need to renormalize F had 1 . We impose the on-shell renormalization condition F had 1 (0) = 0 by subtracting, as a counter term, the integrand itself in the limit k 2 → 0. The final expressions for the kernel functions appearing in the renormalized form factors are: valid in the scattering region t < 0. The inverse hyperbolic tangent appearing in (3.4) and (3.5) are always real-valued if 0 < x < 1 and y < 0. By taking the limit t → 0 in F had 2 , we recover the space-like formula for a HLO µ in eq. (1.3), which is usually derived by applying twice the dispersion relation. Our calculation shows that eq. (1.3) can be obtained directly, without making use of the dispersion relation (it was proven already in [51]). Moreover by substituting Π(q 2 ) →−1 in (3.2) and performing the integration analytically, we reproduce the Pauli form factor at one-loop [50]: where ξ is the Landau variable t/m 2 = −(1 − ξ) 2 /ξ. The same check cannot be done straightforwardly for the Dirac form factor because f 1 (x, y) is not integrable anymore in x = 0 if we set Π had (q 2 ) = −1, while we assumed Π had (0) = 0. To reproduce F 1 (k 2 ) at one loop, we can substitute in (3.2) which corresponds in figure 3 to the exchange of the dressed photon with an undressed one with fictitious mass λ. The integral (3.2) is now finite and the integration can be done analytically. Keeping terms that do not vanish in the limit λ → 0 we correctly recover the known result [50]:

The Boxes
We can now turn our attention to the box diagrams. There are two topologies to take into account: the uncrossed two photon exchange in figure 4a and the crossed one in figure 4b. Figure 4: Irreducible hadronic box diagrams contributing to µ-e scattering at NNLO. Muon and electron lines are depicted with thick and thin lines. Two additional diagrams, with the vacuum polarization in the other photon propagator, must be considered also.
They are related by the crossing s ↔ u plus an overall minus sign. Both photons can be dressed with the hadronic vacuum polarization. The diagram with the same topology but with dressed and undressed photon exchanged can be obtained by replacing the initial state momenta with the final state ones and vice versa. Therefore the contribution of the two diagrams to the unpolarized cross section is the same since no crossing of the Mandelstam variables occurs. Let us fix the notation for the process e − µ − → e − µ − . We choose the following set of propagators: where m 2 (M 2 ), p 1 and p 3 (p 2 and p 4 ) are the mass, the initial state and the final state momentum of the electron (the muon). The Mandelstam variables satisfy s + t + u = 2m 2 + 2M 2 , with the physical requirements We work at the level of interferences between the boxes and the Born amplitude which provide us with the scalar numerators N in eq. (2.2). After algebraic manipulation, these interferences are written as linear combinations of loop integrals that can be evaluated one by one via the hyperspherical method. Each of the box diagrams requires the evaluation of 14 master integrals that have the following form: The kernel functions denoted by arise from the angular integration of their arguments followed by analytic continuation of the external momenta back to the physical region. The solutions of the angular integration (for Euclidean momenta) are taken from the results in refs. [43,44,48]. The necessary angular integrals are the following: 1 1 1 where (4.20) (4.22) At this point we would like to comment on the analytic continuation that we performed in eqs. (4.10-4.18). The general expressions for the angular integrals in [43,44,48] are written in terms of squared Euclidean momenta fulfilling P 2 i > 0 and |P i ·P j | < 1. They must be continued to the on-shell conditions P 2 i = −m 2 i and (P 1 − P 2 ) 2 = −(s + iε), (P 1 − P 3 ) 2 = −t, (P 1 − P 4 ) 2 = −u. Analytic continuation affects the whole radial integral (4.8), not only the angular integration result. Indeed it is necessary to check if any singularity crosses, in the Q 2 complex plane, the integration path along the positive real axis when the P 2 i are continued from positive to negative values. This check is carried out explicitly for the one-loop example in the appendix.
Note that when making use of the results in ref. [43,44,48] one just need to identify the correct side of the branch cut of the functions L(z) after setting P 2 i = −m 2 i . This can be achieved by imposing the correct analyticity structure dictated by unitarity. Indeed after substitution of P 2 i = −m 2 i and (P 1 − P 2 ) 2 = −s etc., all the square roots in front of (4.13-4.18) become real and positive and all arguments z i are real too. The vacuum polarization function is real in the space-like region, therefore the imaginary part of L(z i ) is the only one that must be fixed. Let us consider for example the integrals (4.13) and (4.14) which are obtained by pinching the denominators D 3 or D 1 in the box. These integrals depend on t < 0 and therefore, since they are evaluated below the lepton-pair threshold, their imaginary part must be equal to zero. The imaginary parts of the L(z i ) can be chosen accordingly with this constraint. Similar arguments apply to the integrals for the crossed box diagrams that depend on the Mandelstam variable u.
In addition to the formulae in (4.10-4.18), angular integrals with scalar product q · p i at the numerator are necessary. Thanks to the reduction technique outlined at the end of section 2, they can be written as a linear combination of the integrals (4.10-4.16): Up to this point, we have given an account of the angular integration solutions. We can now introduce the explicit expressions of the radial master integrals that must be evaluated numerically once provided with the hadronic vacuum polarization at negative q 2 . The integrals are the following: with i, j, k = 1, 2, 3 (i, j, k = 1, 2, 4) for the uncrossed box in figure 4a (the crossed box in figure 4b) and i = j = k. The denominator D 0 = q 2 = −Q 2 does not depend on the angles, so we can assemble the kernel functions from the results in (4.10-4.18) straightforwardly.
The following integrals must be considered as well: with i = j, i, j = 1, 2, 3 (1, 2, 4) and k = 3 (4) for the uncrossed box (the crossed box). In eq. (4.35) and (4.36) the kernel functions contain two terms, each of them gives a UV divergent integral if taken alone. To avoid the introduction of an explicit UV regulator, which eventually cancels out in the final result, we take their difference to obtain a UV finite integral. In eqs. (4.33-4.36), the 1/|Q 2 + t| pole in the functions (4.17) and (4.18) yields a singular integral that corresponds to the soft IR divergence arising when the undressed photon becomes soft. Note on the contrary that the 1/Q 2 pole does not lead to a singular integral since the kernel behaviour is smoothed at Q 2 → 0 by the renormalized vacuum polarization. The IR singularity can be regularized by introducing a photon mass λ for the undressed photon. However if we perform the subtraction (2.9) for each integral, we can still employ the formulae presented before to build the kernel functions. Indeed, in eqs. (4.37-4.40) the first integral is now free of IR divergences because the factor Π had (−Q 2 ) − Π had (t) compensates the |Q 2 + t| at the denominator. Therefore we can set λ = 0 in this first term and use our results for the angular integrals. The soft pole appears only in the second terms where standard techniques can be employed for the evaluation of the integrals since Π had does not depend anymore on the loop momentum q.
Note that the simple subtraction (2.9) does not work for I 0ijk in (4.34). The kernel has a 1/Q 2 pole compensated at Q 2 = 0 by Π had (−Q 2 ) but not by a constant term like Π had (t), whereas the factor 2 Q 2 /(Q 2 + |t|) vanishes in the Q 2 → 0 limit while it gives one when Q 2 → |t|.

Dispersive vs Hyperspherical Method
With the formulae for the QED form factors and the boxes in our hand, we can now make a numerical comparison between the standard dispersive approach and the hyperspherical method. The comparison can be done not only with the hadronic vacuum polarization Π had (q 2 ), but also with the well-known analytic expression for Π lep (q 2 ) at one loop, which is a smooth function both for time-like and space-like q 2 .
Numerical integrations, either space-like or time-like, are performed with a Mathematica code employing machine precision numbers and without any symbolic manipulation of the integrand. This ensures that we can compare the two cases and use the same code for Π lep (q 2 ) as well as for Π had (q 2 ). The numerical values of Π had (q 2 ) and the R ratio are provided by the Fortran library alphaQED [52][53][54], and Rhad [55] for the regions where perturbative QCD applies, via a mathlink interface.
We make a comparison with the irreducible diagrams calculated with the dispersive method in [36]. The dispersion relation (1.4) effectively replaces the dressed photon propagator with the propagator of a massive gauge boson. These amplitudes are generated by FeynArts [56] with a modified version of the QED model that contains, besides leptons and photon fields, a massive gauge boson with squared mass equal to z. Later on, the amplitudes are reduced by FormCalc [57,58] to one-loop tensor coefficients which are calculated by the Fortran library Collier [59] via the CollierLink interface [60]. Collier features dedicated expansions in numerically dangerous regions (small Gram or other kinematical determinants). We particularly benefited from the use of this library because in the numerical evaluation of the dispersive integral (1.4) the photon mass z appearing inside the loop  can acquire values a few orders of magnitude larger than the typical energy scales of the scattering process. The numerically stable results provided by Collier in this treacherous region speeded up the convergence of the dispersive integrals.
We begin by comparing the form factors in eq. (3.2). By employing Π lep (q 2 ) instead of Π had (q 2 ), i.e. substituting the hadronic bubble in figure 3 with an electron or a muon loop, we can compare our numerical integration with the analytic results of ref. [61], where the QED form factors at two loops were presented. In [61] the vacuum polarization contribution was calculated with the lepton inside the bubble equal to the external one. The relative difference (F num i /F i ) − 1 between the form factors calculated numerically with the hyperspherical method, F num i = F hyp i , and the exact two-loop result F i is shown in figure 5a and 5b for the electron and the muon case, respectively. The comparison is done for values of the Mandelstam variable t accessible by the MUonE experiment at √ s = 0.405 GeV: −0.142 GeV 2 ≤ t ≤ 0 GeV 2 . The error bars are the uncertainty due to the numerical integration. Harmonic polylogarithms are evaluated with the HPL package [62,63]. In addition to that, we calculated the QED form factors by employing the dispersion relation (1.4) and the analytic expression of Im Π lep (q 2 ). The relative difference with F num  shown as well in figure 5a and 5b. Both methods are in very good agreement with the exact two-loop results, at the level of one part in 10 −8 .
The comparison with Π had (q 2 ) is shown in figure 6a and 6b for the electron and muon form factors, respectively. In this case lacking an "exact" two-loop expression, we choose as normalization F i the dispersive method's result. The values shown in figure 6 are obtained with the same code employed for Π lep (q 2 ), except for the use of Π had (q 2 ) and Im Π had (q 2 ) instead of the leptonic ones. We note that with the hadronic vacuum polarization there is a small systematic shift between the numerical values obtained with the two methods, a relative difference of about 10 −3 − 10 −4 . An improvement of the numerical integration error does not change the picture.
The source of this shift is the following. The function Π had (q 2 ) provided by the library alphaQED is not obtained from a direct integration of the R ratio via (1.4). It is actually calculated in a different way: first different experiments are integrated separately and then weighted averages of the integrals are taken. This procedure appears to be more reliable for error estimate, especially in the ππ channel. The imaginary part provided by alphaQED, i.e. the time-like R, is obtained by averaging data energy-bin-wise. Therefore the numerical integration of this "unified" R can slightly differ from the first procedure since integration and averaging do not commute in general [64]. This effect is shown in figure 7 where we compare, at space-like t, the difference between the hadronic vacuum polarization provided by alphaQED, denoted by Π FJ (q 2 ), and the values obtained by direct integration of the dispersion relation with R from the same library, denoted by Π DR (q 2 ). We note a small difference of the order of 10 −3 , compatible with the systematic shift appearing in figure 6. The two determinations of Π had (q 2 ) are nevertheless in very good agreement within the experimental uncertainty on Π FJ (q 2 ), which is also provided by alphaQED (shown by the orange band in figure 7).
Also the results of the dispersive and hyperspherical methods are in good agreement taking into account the experimental error from the R ratio. The muon and the electron form factor F 2 at t = 0 corresponds to a HLO µ and a HLO e . Their relative uncertainties are about 0.6% [7][8][9], much larger than the discrepancy appearing in figure 6. One should remind however that the kernel functions employed in the dispersive evaluation of F 1 and F 2 at t = 0 are different fromK(s) in the g−2 formula (1.1), so the integration procedure would give in principle a different relative error because the experimental data are weighted differently. The uncertainty on a HLO µ and a HLO e (0.6%) must be understood as an order of magnitude of the error at t = 0 and not as a precise estimate. A explicit calculation of the uncertainty for all t is beyond the scope this analysis: it would require the combination of systematic and statistical errors of the data together with their correlation matrices. However, even assuming in figure 6 that the relative error due to R is 0.1%, which is a factor of six smaller than the uncertainty on a HLO µ and a HLO e , the dispersive and the hyperspherical method would be still in agreement.
Let us move on to the box diagrams.   For each value of t, the dispersive method's result is obtained by performing only one numerical integration: the convolution between the z-dependent virtual-Born interference and the imaginary part of the vacuum polarization. On the contrary, to achieve good numerical stability with the hyperspherical method we had to evaluate separately for each box topology the 14 radial integrals I. Some of the kernel functions are very unstable around Q 2 = |t| and Q 2 = +∞ and therefore a dedicated series expansion must be employed in these regions. The IR divergent integrals in eqs. (4.37-4.40) with the constant term Π(t) in front of it were written in term of one-loop scalar functions and calculated with Collier.

Conclusions
The present error on the hadronic leading order contribution to the muon g−2 constitutes roughly 50% of the error budget in the SM prediction. The MUonE experiment proposed at CERN aims at measuring the running of the fine-structure-constant in the space-like region in µ-e scattering and to determine from it a HLO µ with an error of about 2 × 10 −10 . To reach such level of precision it will be necessary to measure the differential cross section with an uncertainty of the order of 10 ppm. To this end, a Monte Carlo generator with QED and QCD radiative corrections up to NNLO in α must be developed.
In this article we studied the hadronic contributions to the NNLO cross section and we presented a method to evaluate numerically the non-factorizable two-loop diagrams with space-like data for the hadronic vacuum polarization, without making use of the R ratio. In this way the same space-like data measured by MUonE, together perhaps with lattice data and QCD perturbative results, could be exploited to calculate these hadronic corrections. This would allow us to decouple the space-like determination of a HLO µ from any time-like input.
This work took advantage of the hyperspherical integration method, that was described in section 2, to express the irreducible vertex and box corrections as a convolution between the vacuum polarization evaluated at negative q 2 together with a kernel function obtained by analytic integration of the loop diagrams with respect to the hyperspherical angular variables. The vertex corrections were presented in section 3 in terms of QED form factors. In section 4 we showed that each of box contributions can be reduced to a linear combination of 14 integrals which are calculable with the hyperspherical method. Some of these integrals are IR divergent. By making a dedicated subtraction, we managed to remove the IR poles from the integrals explicitly containing the hadronic vacuum polarization and to isolate them in terms that are calculable analytically with standard methods.
Finally, in section 5 we showed that the numerical evaluation of these irreducible diagrams gives results in agreement with the standard dispersive approach and -when the analytic expression of Π lep (q 2 ) is employed -in agreement with analytic two-loop vertex results in QED. A complete calculation of the hadronic corrections to µ-e scattering at NNLO with the dispersive approach will be presented soon [36]. valuable suggestions. I wish to thank also R. Bonciani In this appendix we present an example of a one-loop calculation with the hyperspherical method and we discuss how to perform the analytic continuation between the Euclidean and the physical region. We consider, as an example, the loop integral in eq. (4.15): After continuation of external and internal momenta to the Euclidean region, Wick rotation and the introduction of hyperspherical coordinates, the loop integral is cast in the following form: We expand the propagators as series in Gegenbauer polynomials: We perform the angular integration by making use of the orthogonality property (2.6): (A.7) The series in the expression above can be calculated by defining z = (−Z 1 Z 2 ) and by taking the derivative w.r.t. z, that yields: where τ =P 1 ·P 2 . We then take the primitive and we impose the boundary condition n z n+1 n+1 C (1) n = 0 at z = 0. So the series is: where we used the addition formula arctan(x) − arctan(y) = arctan( x−y 1+xy ). Having performed the angular integrations, the loop integral takes the form: Since ultimately we are interested in the answer for time-like P 2 1 and P 2 2 we have to perform the analytic continuation before the Q 2 -integration. The most important point one has to check is whether any singularity crosses the integration path in the Q 2 complex plane when P 2 1 and P 2 2 are continued to negative values. Barring the poles coming from the divergences in the infinite sums in (A.3) and (A.4), which does not affect this analysis, the integrand is meromorphic in the variables Q 2 , P 2 1 and P 2 2 except for the square roots in the Z variables (A.5) and (A.6). Let's now study the behaviour of Z 1,2 when P 2 1,2 is continued from positive quantities to negative on-shell values P 2 1 = −m 2 and P 2 2 = −M 2 . In the Q 2 complex plane, Z 1 (Z 2 ) has branch points at Q 2 = (P 1 ± im) 2 (Q 2 = (P 2 ± iM ) 2 ). At the beginning P 1 and P 2 are real and positive and the integration is performed from 0 up to ∞. Figure 9 shows the path of the branch points as P 2 1 is varied to −m 23 . When P 2 1 = 0 they are located at Q 2 = −m 2 . When P 2 1 is continued to negative values, one of the branch point moves to the left and the other one reaches the origin at P 2 1 = −m 2 . The branch points of Z 2 behave in the same way. None of the singularities crosses the integration path and therefore we can continue P 2 1 (P 2 2 ) to −m 2 (−M 2 ) without distorting the Q 2 contour. Note however, if we had to continue P 2 1 to a value larger than m 2 , one of the branch point would have crossed the positive real axis and we would have needed to distort the contour to get the correct continuation of the integral (see also the discussion in ref. [40]).
The Euclidean result of the angular integration can be then continued to the on-shell configuration by setting P 2 1 = −m 2 , P 2 2 = −M 2 and (P 1 +P 2 ) 2 = P 2 1 +P 2 2 +2|P 1 ||P 2 |τ = −s, keeping in a first step (M − m) 2 < s < (M + m) 2 in order to leave the square roots in (A.9) real valued. In a second step, we continue s to the physical region s + iε > (M + m) 2 , giving to it a small (positive) imaginary part. As in (A.9) the square root √ 1 − τ 2 becomes i √ τ 2 − 1, we can rewrite the arctangent in terms of L(z) (the hyperbolic inverse tangent) via the identity: arctan(iz) = iL(z). Eventually the kernel function appearing in the Re Q 2 (P + im) 2 (iP + im) 2 (iP − im) 2 (P − im) 2 Im Q 2 Figure 9: Location of the branch points of Z 1 , in the Q 2 plane. The path shows how these branch points moves as P 2 1 is varied from a positive value to −m 2 .
integral (A.10) can be cast in the following form: Let's analyze this formula. As expected, the function L has an imaginary part for s > (M + m) 2 since s is continued above the physical threshold. For real z, the function L(z) acquires an imaginary part when |z| > 1, which happens in the bounded region 0 < Q 2 < λ(s, M 2 , m 2 )/s. Eq. (A.11) provides also the result for s < (M −m) 2 , which corresponds to a u-channel configuration with s substituted by u. In this case it gives the formula (4.16) for the crossed box. One can verify that no imaginary part is developed if s = u < (M −m) 2 , since s is below the physical threshold. Indeed the argument of L is monotonically increasing for Q 2 → +∞ and it is bounded between − (s − (M − m) 2 )/(s − (M + m) 2 ) (at Q 2 = 0) and zero (at Q 2 → +∞). Finally one can show that eq. (A.11) is equivalent to the expression in (4.15) and (4.16), that were derived from eq. 10 in ref. [48], by making use of the identities L(u) + L(w) = L( u+w 1+uw ) and L(z) = L(1/z) + iπ/2 (for |z| > 1). The formula (A.11) in the equal mass case, i.e. m 2 = M 2 , appears also in the calculation of the vertex form factors in eq. (3.4) and (3.5).