Muon-electron scattering at NNLO: the hadronic corrections

The Standard Model prediction for muon-electron scattering beyond leading order requires the inclusion of QCD contributions which cannot be computed perturbatively. At next-to- and next-to-next-to-leading order, they arise from one- and two-loop diagrams with hadronic vacuum polarization insertions in the photon propagator. We present their evaluation using the dispersive approach with hadronic $e^+e^-$ annihilation data and estimate their uncertainty. We find that these corrections are crucial for the analysis of future high-precision muon-electron scattering data, like those of the recently proposed MUonE experiment at CERN.


INTRODUCTION
The elastic scattering of muons and electrons is one of the most basic processes in particle physics. It has been studied since the late 1930s, when measurements of the collisions of muons in cosmic rays with atomic electrons played a crucial role in the discovery of the muon. In spite of this long history, few measurements are available. In the 1960s, the µe elastic scattering cross section was measured at CERN and Brookhaven using acceleratorproduced muons [1][2][3], and in the late 1990s the scattering of muons off polarized electrons was used by the SMC collaboration at CERN as a polarimeter for high-energy muon beams [4].
Recently, a novel approach has been proposed to determine the leading hadronic contribution to the muon g-2, measuring the effective electromagnetic coupling in the spacelike region via scattering data [5]. The elastic scattering of high-energy muons on atomic electrons has been identified as an ideal process for this measurement, and a new experiment, MUonE, has been proposed at CERN to measure the differential cross section of µe elastic scattering as a function of the spacelike squared momentum transfer [6]. In order for this new determination of the leading hadronic corrections to the muon g-2 to be competitive, the µe differential cross section must be measured with statistical and systematic uncertainties of the order of 10ppm. An analogous precision is therefore required in the theoretical prediction.
Until recently, the Standard Model (SM) prediction of the µe → µe process had received little attention. Only the next-to-leading order (NLO) QED corrections to the differential cross section were computed (long time ago) in [7][8][9][10][11][12][13] and revisited more recently in [14]. As a first check, we recalculated these corrections and found agreement with [14] (see also [15,16]). An important step forward was taken very recently by the authors of [17], who calculated the full set of NLO QED corrections with-out any approximation and developed a fully differential Monte Carlo code. They also computed the full set of NLO electroweak corrections.
The QED corrections at next-to-next-to-leading order (NNLO), crucial to interpret the high-precision data of future experiments like MUonE, are not yet known, although some of the two-loop corrections which were computed for Bhabha scattering in QED [18,19] and for tt production in QCD [20,21] can be applied to µe scattering as well. A first step towards the calculation of the full NNLO QED corrections to µe scattering was taken in [22][23][24], where the master integrals for the two-loop planar and non-planar four-point Feynman diagrams where computed. These integrals were calculated setting the electron mass to zero, while retaining full dependence on the muon one. The extraction of the leading electron mass effects from the massless µe scattering amplitudes has been recently addressed in [25] (see also [26][27][28]).
In this letter we will study the hadronic corrections to µe scattering. While at NLO these corrections are simply proportional to the product of the LO QED predictions and the hadronic part of the vacuum polarization, at NNLO their evaluation is complicated by the presence of non-factorizable two-loop diagrams. We will present their calculation using the dispersive approach with hadronic e + e − annihilation (timelike) data [29]. This approach was used, for example, to calculate the hadronic corrections to Bhabha scattering [30,31]. We point out that, taking advantage of the hyperspherical integration method, it was recently shown that these nonfactorizable diagrams can also be calculated employing the hadronic vacuum polarization in the spacelike region, without using timelike data [32]. We will conclude that the hadronic NNLO corrections to µe scattering are important in comparison with the expected accuracy of future high-precision experiments like MUonE, and we will show how to properly include them in their analysis.

DETAILS OF THE CALCULATION
The SM prediction for the unpolarized differential cross section of the elastic scattering at leading order is where m (M ) is the electron (muon) mass, s, t, u are the Mandelstam variables satisfying s+t+u = 2m 2 +2M 2 , α is the fine-structure constant, and λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz is the Källen function.
In a fixed-target experiment where the electron is initially at rest, E µ is the energy of the incoming muons or antimuons, and E e is the electron recoil energy, the variables s and t are given by Both positive and negative muon beams will be available for the MUonE experiment. For E µ = 150 GeV, which is a typical energy available at the M2 beam line in CERN's North Area, s = 0.164 GeV 2 and −0.143 GeV 2 < t < 0.
To evaluate the hadronic corrections to µe scattering at NLO and NNLO, let us consider the SM vacuum polarization tensor with four-momentum q, where j µ is the electromagnetic current and the sum runs over fermions with charges Q f . The weak interactions will be ignored. The transverse part of the full photon propagator is where Π(q 2 ) is the renormalized vacuum polarization function satisfying the condition Π(0) = 0. It receives contributions from the charged leptons (l), the five light quarks u, d, s, c, b with the corresponding hadrons (h), and from the top quark: While the purely leptonic contribution Π l (q 2 ), as well as Π top (q 2 ), involving only the top quark, can be computed order by order in α and α s [33][34][35], the hadronic one cannot be calculated in perturbation theory for any value of q 2 , because of the non-perturbative nature of the strong interaction at low energy. Yet, the subtracted dispersion relation and the optical theorem ImΠ h (s) = (α/3)R(s), where and is the effective fine-structure constant, allow to express the hadronic vacuum polarization in terms of the measured cross section of the reaction e + e − → hadrons [36]. At NLO, the hadronic correction to the µe differential cross section, of order α 3 , is due to the diagram in Fig. 1a. It is the same for positive and negative muon beams: As the leading hadronic contribution to the running of (13) accounts for the leading hadronic effect of the running of the electromagnetic coupling constant in the spacelike region. The extraction of this quantity from µe scattering data is, therefore, the goal of the MUonE experiment. Alternatively, its numerical value can be obtained using the Fortran libraries alphaQEDc17 [37][38][39][40] and KNT18VP [40][41][42][43][44][45] based on hadronic e + e − annihilation (timelike) data. Let us define the ratio For E µ = 150 GeV, it reaches the maximum value of 2.1 × 10 −3 at t = t min = −0.143 GeV 2 . For the same incoming muon energy, the maximum value of the top quark contribution is a tiny 2.0 × 10 −9 . At NNLO, the hadron-induced corrections to µe scattering, of order α 4 , can be split into four classes of diagrams: I. Class I consists of tree level diagrams with one or two vacuum polarization insertions (Fig. 1b). Their contribution to the differential cross section is pro- , the reducible part of the second-order hadronic contribution to the running of α(t). We note that, in Fig. 1a, a virtual photon can be emitted and reabsorbed by the hadronic insertion. This irreducible part of the second-order hadronic contribution to the running of α(t) will not be considered as part of class I (although of the same order), because its effect is commonly included in the ratio R(s) as final-state radiation and, therefore, it is already incorporated in the NLO hadronic corrections in Eq. (13) [46,47].
II. QED one-loop diagrams in combination with one hadronic vacuum polarization insertion in the tchannel photon (Fig. 1c). Their contribution to the differential cross section is proportional to Π h (t) and a combination of one-loop QED corrections to µe scattering.
III. Real photon emission diagrams with a vacuum polarization insertion in the t-channel photon (Fig. 1d). They contain terms proportional either to Π h (t e ) or to Π h (t µ ), where t e (t µ ) is the square of the difference of the initial and final electron (muon) momenta. In general, t e = t µ = t because of the presence of the final-state photon.
All the diagrams in classes I-III are factorizable, since each of them can be reduced to the product of a QED amplitude multiplied by the function Π h (q 2 ) evaluated at some q 2 value fixed by the external kinematics. A fourth class of non-factorizable diagrams must also be considered: IV. One-loop QED amplitudes with a hadronic vacuum polarization insertion in the loop. They can be further subdivided into vertex and box corrections (Fig. 1e).
We remind the reader that there are no light-by-light contributions to the µe cross section at NNLO (order α 4 ). They appear at N 3 LO (order α 5 ). We calculated the amplitudes in class IV employing the dispersion relation in Eq. (10). The factor Π h (q 2 )/q 2 appearing in the loop -where q now stands for the loop momentum -is replaced by the r.h.s. of Eq. (10), where q appears only in the denominator of the term 1/(q 2 −z). Therefore, the dispersion relation effectively replaces the dressed propagator with a massive one, where z plays the role of a fictitious squared photon mass. This allows to interchange the integration order and evaluate, as a first step, the one-loop amplitudes with a "massive" photon. The results obtained for the z-dependent scattering amplitudes are then convoluted with the ratio R(s).
All four classes of diagrams were generated using FeynArts [48] with a modified version of the QED model that contains, besides leptons and photons, a fictitious massive gauge boson (the "massive" photon arising from the dispersion relation). The amplitudes were calculated and reduced to one-loop tensor integrals with Form [49] via the FormCalc [50] package, and exported as a Fortran code for the numerical evaluation of the dispersive and phase-space integrals. Two independent parametrizations of the 3-body phase space were employed to crosscheck the hard bremsstrahlung cross section. For the numerical evaluation of Π h (q 2 ) in the spacelike region, appearing in classes I-III, we relied on the native implementation available in the Fortran libraries alphaQEDc17 and KNT18VP. The one-loop tensor coefficients were computed with the library Collier [51], which features dedicated expansions for the evaluation in numerically unstable regions (small Gram or other kinematical determinants). We particularly benefited from this library when we convoluted the z-dependent amplitudes with the R(s) ratio provided by alphaQEDc17 or KNT18VP. Indeed, in performing the dispersive integrals in class IV diagrams, the squared photon "mass" z appearing inside the loop functions can acquire values which are orders of magnitude larger than the typical energy scale of the scattering process. Collier provides numerically stable results in this treacherous region and allows the numerical integration to converge. The dispersive integrals were performed with the subroutines in QUADPACK [52], while for the phase space integration we employed the VEGAS algorithm [53] in the Cuba library [54]. To check our results, we produced an independent Mathematica implementation using FeynCalc [55,56] and Package-X [57]. The results obtained by FeynCalc in terms of scalar one-loop functions were then evaluated numerically using analytic expressions provided by Package-X. The use of Mathematica's arbitraryprecision numbers, with a large number of digits, allowed us to keep track of precision at all steps and avoid instabilities during the numerical dispersive and phase-space integrations. We found perfect agreement between the two implementations.
The lepton masses were kept different from zero throughout the calculation, so that the matrix elements were free of collinear singularities. Ultraviolet singularities were regularized via conventional dimensional regularization and UV-finite results were obtained in the onshell renormalization scheme. The amplitudes of class II and the boxes of class IV develop IR poles which are cancelled by those arising from the phase space integration of the real emission diagrams of class III. We employed both the FKS subtraction scheme [58,59] as well as the traditional QED procedure to assign a vanishingly small mass to the photon to remove the soft singularities and to obtain an IR-finite cross section.

RESULTS
The ratio of the NNLO hadronic corrections to the µe differential cross section, with respect to the squared momentum transfer t e , and the LO prediction, is shown in Fig. 2 for the processes µ + e − → µ + e − (upper panel) and µ − e − → µ − e − (lower panel) for E µ = 150 GeV. The red lines indicate the total hadronic contribution arising from classes I-IV, while the blue ones show the sum of the contributions of classes II, III, and IV, but not I. The reason for this split is the following. The goal of the MUonE experiment is to determine ∆α h (t) = −Π h (t), the leading hadronic contribution to the running of the effective fine-structure constant in the spacelike region, from µe scattering data. In order to extract the NLO hadronic correction to the µe differential cross section, given by Eq. (13), which contains Π h (t), the experimental data will have to be subtracted, via a Monte Carlo event generator, of the total NNLO hadronic corrections (classes I-IV). If, instead of ∆α h (t), one wants to extract the hadronic corrections to the resummed photon propagator, then the corrections of class I should not be subtracted from data, as their contribution to the differential cross section accounts for the second-order reducible hadronic contribution to the running of α(t). Figure 2 shows that, when the muon/antimuon beam has an energy of 150 GeV, for most of the kinematic region scanned by the squared momentum transfer t e the factor K NNLO h (t e ) is of order 10 −4 -10 −5 . These corrections are therefore larger than the O(10 −5 ) precision expected at the MUonE experiment. Moreover, our Fortran code, available upon request, can calculate the NNLO hadronic corrections to any µe scattering differential distribution with arbitrary kinematical cuts and can therefore be implemented in future full NNLO µe scattering Monte Carlo codes.
At NLO, the tiny contribution of the top quark to the vacuum polarization can be separated from the hadronic one. At NNLO, these contributions mix with each other. The plots in Fig. 2 were obtained adding Π top (q 2 ) to Π h (q 2 ), so that the full top quark contribution has been included in the shown NNLO prediction. Its effect is however totally negligible.
As our calculation of the NNLO hadronic corrections to the µe differential cross section is based on the hadronic e + e − annihilation data, the precision of our prediction is limited by the experimental error on the R(s) ratio. We estimated the uncertainty of our results induced by this error by comparing the values obtained with the libraries alphaQEDc17 and KNT18VP. For each value of t e , we found that the relative difference between the two calculations of dσ NNLO h /dt e is about 1% or less. We therefore assigned to our dσ NNLO h /dt e predictions a relative uncertainty of 1%, which corresponds to an error in K NNLO h (t e ) of O(10 −6 ) or less, well below the precision expected at the MUonE experiment.
By employing the well-known one-loop analytic expression for Π l (q 2 ) instead of Π h (q 2 ), i.e. substituting the hadronic blob in Fig. 1e with an electron or a muon loop, we compared our results for the vertices in class IV with the two-loop analytic expressions for the QED form factors of Ref. [60]. We found perfect agreement between our numerical dispersive integrations and their explicit two-loop results (see also [32]).

CONCLUSIONS
In this letter we presented the NLO and NNLO hadronic corrections to the differential cross section for the processes µ ± e − → µ ± e − (+γ), where (+γ) indicates the possible emission of photons.
We showed that, in a fixed-target experiment where the electron is initially at rest and the energy of the incoming muons or antimuons is 150 GeV, the corrections to the differential scattering cross sections with respect to t e , the square of the difference of the initial and final electron momenta, are of order 10 −4 -10 −5 for most of the kinematic region spanned by t e . These corrections will therefore play a crucial role in the data analysis of future high-precision muon-electron scattering experiments like MUonE, whose goal is to reach a relative precision of order 10 −5 . The relative theoretical uncertainty of our predictions, induced by the experimental error of the hadronic e + e − annihilation data, is estimated to be about 1% or less. It is therefore well below the precision expected at the MUonE experiment.
Making use of crossing relations, our results are also relevant to muon-pair production at present and future e + e − colliders operating at high and low energies.