M{\o}ller scattering at NNLO

We present a calculation of the full set of next-to-next-to-leading-order QED corrections to unpolarised M{\o}ller scattering. This encompasses photonic, leptonic, and non-perturbative hadronic corrections and includes electron mass effects as well as hard photon radiation. The corresponding matrix elements are implemented in the Monte Carlo framework McMule allowing for the computation of fully-differential observables. As a first application we show results tailored to the kinematics and detector design of the PRad II experiment where a high-precision theory prediction for M{\o}ller scattering is required to achieve the targeted precision. We observe that the corrections become essential to reliably calculate the corresponding differential distributions especially in regions where the leading-order contribution is absent.


Introduction
Recent developments at the low-energy precision frontier renewed interest in high-precision calculations in QED. In order to meet the experimental accuracy when using state-of-the-art lepton beams, next-to-leading-order (NLO) or even next-to-next-to-leading-order (NNLO) corrections in QED have become highly desirable.
Electron-electron or Møller scattering is a prime example where a high-precision theory prediction has become important. As it is a ubiquitous process at electron beam lines, it has been investigated in connection with luminosity measurements [1,2] and background studies [3], including a recent dedicated measurement at very low energies to study electron mass effects [4]. Precise knowledge of QED effects might also be useful for the MOLLER experiment [5] that searches for parity violation as an indication towards new physics. The main motivation, however, is given by PRad [6,7] to which we turn below.
All of these experiments rely on Standard Model theory predictions for Møller scattering. Since at low energies the corresponding radiative corrections are QEDdominated, a high-precision theory calculation is feasible. To account for non-trivial detector geometries and acceptances these experiments rely on Monte Carlo event generators that combine matrix elements to physical observables. At NLO accuracy in QED such Monte Carlo codes have been developed [8][9][10]. While for most experiments this level of precision is sufficient, the situation is different for the planned PRad II experiment [7], the upgraded version of its predecessor PRad [6].
Both experiments measure the elastic scattering of electrons and protons to extract the proton charge radius. These efforts were triggered after a recent high-precision spectroscopy experiment using muonic hydrogen at the Paul Scherrer Institute [11,12] had measured a signifi-cantly smaller value for the proton charge radius compared to earlier results [13]. The PRad and PRad II experiments therefore play an important role in the resolution of this 'proton radius puzzle'.
A key feature of the PRad experiments is the suppression of systematic uncertainties by normalising to a simultaneous measurement of Møller scattering. To be precise, the differential cross section for electron-proton scattering (dσ/dθ) exp ep is extracted from the measured events N exp and the theoretical prediction for the Møller cross section (dσ/dθ) th ee via [14] dσ dθ It is therefore essential that the theory prediction for Møller scattering matches the experimental precision. As alluded to before, the existing NLO Monte Carlo event generators were sufficient for PRad. This, however, is not the case for PRad II that aims at improving the experimental precision significantly. It is therefore expected that missing higher-order QED corrections become the dominant source of systematic uncertainty [7].
We have therefore calculated the full set of NNLO QED corrections to Møller scattering. This includes not only photonic and leptonic corrections but also nonperturbative hadronic contributions. The corresponding matrix elements were implemented in the McMule framework, a Monte Carlo for MUons and other LEptons [15]. This allows for the computation of arbitrary fully-differential observables. This paper is organised as follows: we begin by briefly summarising our calculational methods in Section 2. Next, we present results tailored to the kinematics and detector design of the PRad II experiment in Section 3 and conclude in Section 4.

Calculation
We have calculated the full set of NNLO QED corrections to the process including effects due to the non-vanishing electron mass. These corrections are composed of three parts, as illustrated in Figure 1. First, double-real corrections are obtained by integrating the tree-level matrix element (squared) with two additional photons in the final sate over the photon phase space. Second, real-virtual corrections require the interference of the one-loop amplitude and tree-level amplitude with one additional photon. Finally, the double-virtual corrections involve the interference of the two-loop amplitude for e e → e e with its tree-level amplitude as well as the corresponding oneloop amplitude squared. These three individual parts are physically indistinguishable for soft photons and infrared divergent. Only the combination of all contributions leads to a finite, physical result. The corrections can be split into purely photonic contributions and fermionic ones that are due to vacuum polarisation (VP). The fermionic part takes into account all three lepton flavours as well as non-perturbative hadronic effects, indicated by the grey blob in Figure 1b. The corresponding method will be presented in the second part of this section after discussing our approach towards photonic corrections first. The double-real contribution of the fermionic part, shown top left in Figure 1b corresponds to the process with additional leptons (or hadrons) in the final state. Since this is a final state with a measurable difference it can be disentangled from Møller scattering. It is separately finite if fermion masses are not neglected. Hence, we have chosen not to include the double-real fermionic part in our results presented below. The matrix elements are regularised in d = 4 − 2ǫ dimensions and renormalised in the on-shell scheme. They are implemented in the parton-level integrator McMule that is based on FKS ℓ [16], a QED extension of the FKS subtraction scheme [17,18] beyond NLO. McMule allows to consistently perform the phase-space integration in the presence of soft singularities and to calculate physical observables in a fully-differential way.

Photonic corrections
All tree-level and one-loop photonic matrix elements were calculated analytically with the full dependence on the electron mass. The diagrams were generated with QGraf [19] and calculated with the Mathematica code Package-X [20]. While for most of the obtained expressions a sufficiently stable and efficient implementation was possible, a different approach had to be taken in the case of the numerically delicate real-virtual matrix element. In this case next-to-soft stabilisation, previously developed in the context of Bhabha scattering [21], was used. This method is based on the observation that the numerical instabilities are strongly enhanced when the emitted photon becomes soft. We have therefore expanded the real-virtual matrix element for small photon energies including the non-universal next-to-soft contribution. This allows to rely on OpenLoops [22,23] in its standard mode for the bulk of the phase space and otherwise to switch to the next-to-soft approximation. This approach therefore yields a stable and fast implementation of the problematic real-virtual contribution.
The two-loop matrix element for Møller scattering can be derived from corresponding Bhabha results via the crossing relation p 2 ↔ −p 4 . However, even though Bhabha scattering is one of the best studied processes in particle physics, the exact two-loop contribution is currently not known. The approximation where the electron mass m is set to zero was calculated some time ago [24]. Building upon this result it was later possible to determine the leading mass effects based on the universal factorising structure of collinear divergences [25][26][27]. In this approximation power suppressed terms of O(m 2 /q 2 ) are neglected with q 2 corresponding to any of the Mandelstam invariants s = (p 1 + p 2 ) 2 , t = (p 1 − p 3 ) 2 , and u = (p 1 − p 4 ) 2 . Even for low-energy observables this approximation is well justified due to the small electron mass. This massification therefore allows to straightforwardly calculate both the logarithmically enhanced and the constant mass terms based on the massless result. Recently, this procedure was extended to processes that include a heavy mass [28]. In this more general approach a factorisation anomaly [29,30] was discovered in fermionic contributions resulting in the breaking of naive factorisation. This in turn explains the occurrence of additional powers of large logarithms that result in a larger massification error. We have therefore only used the massification procedure for the photonic two-loop contribution and followed a different approach in the fermionic case.
The analytic continuation that has to be performed after crossing is non-trivial. Fortunately, this step was de facto already performed in the case of the massless twoloop matrix element [24] as it was calculated for all three possible kinematic channels for 2 → 2 processes. It was therefore straightforward to construct the massless twoloop matrix element for Møller scattering based on this. We then massified the corresponding expression ourselves arriving at the desired result.
We have verified the non-radiative part of our calculation at NLO by comparing to corresponding results in [31] and found exact agreement. We have further compared to full NLO results provided by [32] using the Monte Carlo code from [10]. On account of this code being an event generator and not a full dedicated integrator, we found agreement, up to Monte Carlo errors, at the level of precision achieved. To still allow for a robust check of our calculation we have implemented all photonic matrix elements in McMule such that they can directly be crossed to Bhabha scattering. This then allows us to verify our calculation indirectly by comparing to the much richer literature of Bhabha scattering. In particular, we have compared to the state-of-the-art Monte Carlo generator BABAYAGA that is based on a parton shower algorithm matched to the exact NLO result [33]. As presented in [21], we have found exact agreement at NLO and a deviation at the level of 17% for the NNLO correction, consistent with the expected precision of the logarithmic approximation from the parton shower.

Fermionic corrections
As mentioned above it is not sensible to apply the methods used for photonic corrections to the fermionic contributions because of the occurrence of the factorisation anomaly in the massification procedure. Moreover, at low energies non-perturbative hadronic contributions become relevant. Fortunately, all leptonic and hadronic contributions to Møller scattering are due to VP insertions. We can therefore calculate all of them simultaneously by including all contributions in the VP The one-loop contribution to the leptonic VP can be calculated trivially and the two-loop result can be extracted from [34]. For the non-perturbative hadronic VP, on the other hand, we rely on the Fortran library alphaQED [35][36][37] that calculates Π had based on experimental data. The full photon propagator −ig µν q 2 1 − Π(q 2 ) = −ig µν q 2 1 + Π(q 2 ) + Π 2 (q 2 ) + ... (5) automatically resums VP insertions to all orders in perturbation theory. In this paper, however, we take a strict fixed-order approach for all contributions. In what follows we therefore use the expanded version of the photon propagator given on the r.h.s. of (5). The impact of missing higher-order effects -together with resummation of soft and collinear photon emission -is left for future studies.
For a large class of diagrams, such as the two diagrams on the right of Figure 1b, the VP contributions factorise and can thus be calculated straightforwardly based on photonic NLO matrix elements. The verification of the photonic corrections described in Section 2.1 therefore also presents a strong check in this case. The calculation of the non-factorisable vertex and box diagrams, however, is more involved. Traditionally these contributions are calculated dispersively, see e.g. [38]. In our case we have instead used the hyperspherical formalism [39,40] where it was possible to reuse many of the results of the analogous calculation in muon-electron scattering [41]. In this approach the non-factorisable two-loop contribution is written as an integral over the radial part of the Wickrotated loop momentum which can be performed numerically. The kernel function K results from the analytic integration over the solid angle of the loop momentum. In the case of triangles this integration can be elegantly performed after writing the propagators in terms of Gegenbauer polynomials and making use of the corresponding orthogonality relation. For box integrals, however, this method does not work and the corresponding calculation is therefore more complicated. In this case the general result can be found in [42].
The main advantage of the hyperspherical approach compared to the dispersive calculation is that the VP is integrated over space-like momenta avoiding hadronic resonances. Furthermore, the threshold singularities are transformed into logarithmic divergences in the kernel functions through the angular integration. Nevertheless, the numerical integration over Q 2 has to be performed with care. Since large cancellations can occur in the kernel functions, they have to be expanded in these problematic regions. In particular, this is the case for Q 2 → 0, Q 2 → ∞, and around threshold singularities.
We have verified our calculation of the non-factorisable fermionic contributions by comparing to corresponding crossed results for Bhabha scattering. First of all, we have obtained complete agreement in the case of closed electron loops with the exact analytic calculation in [43]. Furthermore, we have compared the muon and the tau contributions to approximate results from [44] that neglect terms of O(m 2 f /q 2 ) with m f the mass of any of the fermions and q 2 ∈ {s, t, u}. As expected we do not find exact agreement but instead observe the correct converging behaviour when the energy scale q 2 is increased.

Results
With the contributions discussed in the previous section implemented in McMule we can compute any infrared safe observable of Møller scattering at NNLO in QED.
To illustrate this, as a first application we use a kinematic set up tailored to the PRad II experiment. We consider an electron beam with energy E b = 1.4 GeV incident on a target electron at rest. We refer to the outgoing electron with the smaller (larger) scattering angle as the 'narrow' ('wide') electron. The corresponding scattering angles are denoted by θ n and θ w , respectively. We further define the inelasticity η = E b + m − E n − E w and the coplanarity ζ = 180 • − |φ n − φ w | with E i and φ i the energy and azimuth of the narrow and wide electrons. Both quantities are zero for elastic events and can therefore be used to restrict hard photon emission. We approximate the experimental design with the cuts 0.5 • < θ n , θ w < 6.5 • , η < 3.5σ E , ζ < 3.5σ φ , with σ E = 37.7 MeV and σ φ = 2.1 • the expected detector resolutions for the considered beam energy [45]. In addition to this fiducial observable we also consider a simpler version that does not restrict hard photon radiation.
To this end we convert the angular cuts in (7) to corresponding restrictions for the t-and u-channel momentum transfers using the kinematic relations for elastic scattering. The corresponding cuts for this simplified scenario then read The order-by-order contributions, σ (i) , to the integrated cross section, σ 2 = σ (0) + σ (1) + σ (2) , for both the fiducial and the simplified cuts are presented in Table 1. The photonic and fermionic corrections are given separately and are denoted by σ (i) γ and σ (i) VP , respectively. Additionally, we show the corresponding K factors defined as As expected, the results for the fiducial and the simple cuts agree for the elastic corrections σ (0) and σ VP . We further observe large NLO contributions compared to the naive counting of the expansion parameter α/π. The NNLO K factors, on the other hand, are much smaller and in agreement with the naive prediction. For the total cross section we can thus conclude that the NNLO result yields a robust prediction for the considered observable with missing higher-order corrections well under control. Comparing the fiducial and simple results reveals that the inelasticity and coplanarity cut have a noticeable impact on the higher-order corrections. If hard photon emission was even more severely restricted, the NNLO corrections, too, would become large and the resummation of corresponding large logarithms would be required.
For the more realistic fiducial observable we present differential distributions in Figure 2. Figure 2a shows the corresponding results with respect to the energy of the 'narrow' electron. The angular distributions for both the 'narrow' and 'wide' electron are given in Figure 2b. For both figures the differential cross section is displayed in the upper panel. In addition, the middle panel shows the differential version of the K factor defined in (9). In regions where the LO cross section is zero the NNLO contribution effectively corresponds to an NLO correction resulting in a large K factor. For this reason the lower panel zooms into the other region where the NNLO K factor is small. In this case we observe a similar behaviour as for the total cross section with large NLO and small NNLO corrections. Contrary to the integrated case, however, the NNLO calculation is indispensable for distributions in order to obtain a reliable prediction also in the region where the LO contribution is zero.
We have presented integrated and differential cross sections for only one possible scenario for PRad II. Other configurations of interest to PRad and PRad II will be studied in the future and made available on https://mule-tools.gitlab.io/user-library The data that enter the results presented in this section can also be found there.

Conclusion
We have calculated the full set of NNLO QED corrections to Møller scattering including photonic, leptonic, and non-perturbative hadronic corrections. The fermionic part was calculated with a semi-numerical approach using the hyperspherical formalism without any approximation. The photonic two-loop matrix element, on the other hand, was computed based on the crossed massless Bhabha result via massification resulting in a parametrically small error of O(α 2 m 2 /q 2 ) relative to LO. Apart from the calculation of the two-loop matrix element, the main challenge in the computation of the photonic corrections is the numerical stability of the real-virtual matrix element. We have obtained a fast and stable implementation of this delicate contribution using a combination of OpenLoops and the next-to-soft stabilisation method recently developed in the context of Bhabha scattering.
All matrix elements were implemented in the partonlevel Monte Carlo integrator McMule. While not an event generator, this framework allows to compute arbi-   (7). trary fully-differential observables by consistently removing soft divergences originating in the phase-space integration. As a first application of our NNLO calculation we have generated integrated and differential cross sections relevant for the planned PRad II experiment. For the total cross section we obtain large NLO and small NNLO corrections ensuring a high-precision prediction. In the case of differential distributions the situation is more subtle. In the absence of a LO contribution the NNLO correction becomes essential to achieve a moderate precision of about 10%. To improve this theory accuracy further missing higher-order corrections could be approximated with a parton shower.