NNLO virtual and real leptonic corrections to muon-electron scattering

The recently proposed MUonE experiment at CERN aims at providing a novel determination of the leading order hadronic contribution to the muon anomalous magnetic moment through the study of elastic muon-electron scattering at relatively small momentum transfer. The anticipated accuracy of the order of 10ppm demands for high-precision predictions, including all the relevant radiative corrections. The fixed-order NNLO radiative corrections due to the emission of virtual and real leptonic pairs are described and their numerical impact is discussed for typical event selections of the MUonE experiment, by means of the upgraded Monte Carlo code MESMER.


Introduction
The value of the anomalous magnetic moment of the muon, a µ = (g−2) µ /2, is a fundamental quantity in particle physics, that has been very recently measured by the Fermilab Muon g − 2 Experiment (E989) with the unprecedented relative precision of 0.46ppm [1] for positive muons, improving on the 0.54ppm precision obtained by the E821 experiment at the Brookhaven National Laboratory [2].
The muon anomaly, i.e. the relative deviation of the magnetic moment from the value predicted by Dirac theory, is due to quantum loop corrections stemming from the QED, weak and strong sector of the Standard Model (SM) [3,4]. Hence, the comparison between theory and experiment provides a very stringent test of the SM and a deviation from the SM expectation is a monitor for the detection of possible New Physics signals.
There is a long-standing and puzzling muon g − 2 discrepancy between the measured value and the theoretical prediction, which presently exceeds the 3σ level. In particular, the new E989 result is 3.3σ greater than the theoretical prediction. This excess grows up to 4.2σ when the new experimental result is combined with previous measurements using both µ + and µ − beams [1]. The current status of the SM theoretical prediction for the muon g − 2 has been recently reviewed in ref. [5].
Future runs of the E989 experiment at Fermilab [6,7] and of the E34 experiment under development at J-PARC [8,9], are expected to bring the relative precision of the anomalous magnetic moment of the muon below the level of 0.2ppm. This calls for a major effort on the theory side in order to reduce the uncertainty in the SM prediction, which is dominated by non-perturbative strong interaction effects as given by the leading order hadronic correction, a HLO µ , and the hadronic light-by-light contribution. Actually, there is at present general consensus that the latter has a small impact on the theoretical uncertainty on a µ [5] (see also refs. [10][11][12] for very recent investigations through dispersion relations and Lattice QCD calculations).
Traditionally, a HLO µ has been computed via a dispersion integral of the hadron production cross section in electron-positron annihilation at low energies [5,[13][14][15]. More recent discussions on issues and constraints on data used in dispersion relations can be found in refs. [16,17]. Lattice QCD calculations are providing alternative evaluations of the leading order hadronic contribution [18][19][20][21][22][23][24][25][26][27][28][29][30][31]. Lately the BMW collaboration has presented a precise determination of a HLO µ with an uncertainty of 0.78% [32], a central value larger than the ones obtained via dispersive approaches, by an amount that ranges from 2 to 2.5σ depending on the reference value used for the dispersive approach, and in agreement with the BNL and FNAL experimental determinations [1,2]. To clarify this tension, alternative and independent methods for the evaluation of a HLO µ are therefore more than welcome and desirable.
Recently, a novel approach has been proposed in ref. [33] to derive a HLO µ from a measurement of the effective electromagnetic coupling constant in the space-like region via scattering data, making use of a relation between ∆α had (q 2 ) at negative squared momenta and a HLO µ (see also [34]). Shortly afterwards, the elastic scattering of high-energy muons on atomic electrons has been identified as an ideal process for such a measurement [35] 1 and a new experiment, MUonE, has been proposed at CERN to measure the differential cross section of µe scattering as a function of the space-like squared momentum transfer [39][40][41][42]. In order for this new determination of a HLO µ to be competitive with the traditional dispersive approach, the uncertainty in the measurement of the µe differential cross section must be of the order of 10ppm.
The MUonE challenging experimental target requires very high-precision predictions for µe scattering, including all the relevant radiative corrections, implemented in fully-fledged Monte Carlo (MC) tools. The necessary perturbative accuracy will be next-to-next-to leading order in QED (NNLO QED), combined with the leading (and, eventually, nextto-leading) logarithmic contributions due to multiple photon radiation. In recent years, a number of steps have been already taken to achieve this goal. A comprehensive review of the theoretical knowledge, as of the beginning of year 2020, of the µe scattering cross section for MUonE kinematical conditions has been published in ref. [43]. We list below the main milestones already reached in this theoretical endeavour. In ref. [44], the full set of NLO QED and one-loop weak corrections was computed without any approximation and implemented in a fully exclusive MC generator. It is presently being used for simula-tion studies of MUonE events in the presence of QED radiation. Important results were also obtained at NNLO accuracy in QED. The master integrals for the two-loop planar and non-planar four-point Feynman diagrams were computed in refs. [45,46], by setting the electron mass to zero while retaining full dependence on the muon mass. A general procedure to extract leading electron mass terms for processes with large masses, such as muon-electron scattering, from the corresponding massless amplitude was given in ref. [47] and supplemented with a subtraction scheme for QED calculations with massive fermions at NNLO accuracy [48]. More recently, the exact NNLO photonic corrections along the electron line, including all finite mass terms, have been implemented in two independent MC tools, Mesmer [49] and McMule [50]. In the former, the framework for the complete photonic NNLO calculation has been built and made available, including the full NLO calculation to the process µe → µeγ and the LO double radiative process µe → µeγγ.
Since the two-loop diagrams where at least two virtual photons connect the electron and muon lines are not completely known yet, their infrared (IR) part has been taken into account by means of the classical Yennie-Frautschi-Suura approach [51]. Very recently, the analytic evaluation of the two-loop corrections to the amplitude for the scattering of four fermions in QED has been carried out keeping the complete dependence on the mass of one fermionic current [52]. The two-loop hadronic corrections to µe scattering were computed in refs. [53,54]. Also possible contamination from New Physics effects has been studied in refs. [55,56] and shown to be below the MUonE sensitivity, thus reinforcing the robustness of the proposed approach for a reliable determination of a HLO µ through MUonE data. In this paper, we present the calculation of the complete fixed-order NNLO QED corrections which include at least one leptonic pair, with all the relevant virtual and real contributions, summing over all contributing lepton flavours. This class of corrections is named NNLO leptonic correction. In section 2 we describe the classification of the contributions and the methods used for their calculation. Where needed, the calculation of the virtual corrections is performed by means of the dispersion relation technique 2 already used for the hadronic NNLO corrections in refs. [53,54], taking into account all finite mass effects. The real-pair emission contributions µ ± e − → µ ± e − e + e − and µ ± e − → µ ± e − µ + µ − are calculated by means of exact tree-level matrix elements, including electron and muon finite mass effects. For completeness, we calculate also the NNLO virtual hadronic corrections, already calculated in refs. [53,54], and compare their size with the leptonic corrections on all relevant differential distributions. The whole class of corrections has been included in an upgraded version of the Mesmer 3 MC code. By means of the developed MC code we illustrate in section 5 numerical results relevant for typical running conditions and event selections of the MUonE experiment.
The effects induced by the insertion of a leptonic loop are discussed in sections 5.1, 5.2 and 5.3 and the impact of real lepton-pair emission is discussed in sections 5.4 and 5.5, with attention to the interplay between real and virtual pair radiation as well as to the role played by the "peripheral" diagrams within realistic MUonE event selection. A summary and future prospects are finally given in section 6. The work presented in this paper represents an additional step towards the implementation of a fully-fledged MC generator including the complete set of NNLO corrections matched to multiple photon emission.

Classification of NNLO leptonic contributions
The complete set of NNLO leptonic corrections to µ ± e − → µ ± e − scattering, which we denote with dσ α 2 N f , consists of three parts, with contributions from virtual and real leptonic corrections, which can be schematically indicated with the following equation For later use, we remind that the running of the QED coupling constant α through Dyson resummation of 1PI diagrams, in the on-shell renormalisation scheme which is adopted here, reads where ∆α(q 2 ) can be expanded in α as In the previous equation the sum of the perturbative bits runs over all leptons and top quark, while the non-perturbative part ∆α hadr (q 2 ), due to the exchange of hadronic states, is traditionally evaluated through dispersion relations and the optical theorem. We remark that in QED including the running of α in a scattering amplitude amounts to replace the photon propagator of virtuality q 2 , multiplied by the vertex couplings, according to .
We detail below the content of the three contributions of eq. (2.1): Figure 2: Left: NNLO reducible photon vacuum polarisation contribution. Right: QED NNLO irreducible vacuum polarisation contribution.  • dσ α 2 virt , the virtual two-loop contribution, consists of several classes, both factorisable and non-factorisable: 1. the squared absolute value of the leading-order (LO) amplitude A 0 with one-loop photon vacuum polarisation (VP) due to leptons (figure 1 (b)); 2. interference between LO diagram and two-loop diagrams with iterated photon VP ∆α LO as well as two-loop irreducible VP ∆α NLO (figure 2).
The sum of the two contributions in items 1 and 2 is factorisable, i.e. it is equal   b) box corrections, which are IR divergent ( figure 6). The sum of these corrections, the ones in items 3 and 4 cancel the IR divergences of the sub-set dσ α 2 γ , described below; • dσ α 2 γ includes the interplay between real photon radiation and leptonic loop insertions, in particular the interference between tree-level diagrams for the process µ ± e − → µ ± e − γ and the same class of diagrams with insertion of the leptonic loop on the virtual photon (figure 7). IR divergences stemming from this set of corrections are cancelled by a sub-set of the contributions considered in the previous item; • dσ α 2 real contains the O(α 2 ) tree-level amplitudes for the processes µ ± e − → µ ± e − + − , with = e, µ. The production of the τ + τ − pair is forbidden by the available centre of mass energy in the MUonE experiment. We have three different gauge-invariant classes: 1. lepton-pair production through photon radiation from the electron line ( The real-pair radiation contribution dσ α 2 real is particularly important for two reasons: on the pure phenomenological side, when two final state particles escape the MUonE acceptance, the 4-fermion final state mimics a signal event and therefore its contribution is a reducible background and has to be thoroughly investigated under typical MUonE running conditions. On the other side, the interference between different real pair radiation diagrams is expected to give a partial cancellation with its virtual counterpart. For instance, it is known from the literature on Bhabha scattering [60][61][62][63][64][65][66][67] that the interference between diagrams (a) and (b) of figure 8, once integrated over the phase space available to the electronic pair, develops a term proportional to α 2 log 3 (−t/m 2 ), which cancels when added to the virtual contribution originating from diagram (a) of figure 5. Similar partial cancellations are expected between virtual contributions and the interference between the diagrams corresponding to the cut of the closed leptonic loop. Moreover, a careful study of the impact of the peripheral diagrams is necessary in view of their potentially large contribution, as will be discussed in the following sections.

Calculation of virtual and real-virtual corrections
We build our NNLO virtual leptonic corrections by means of dispersion relation (DR) techniques, starting from the available NLO calculation implemented in Mesmer. According to refs. [57][58][59]68], and the seminal work of refs. [69,70], an amplitude involving a photon in a loop with a VP insertion, due to lepton , is obtained by replacing the photon propagator as follows −ig µν where the renormalised VP function Π(q 2 ) is obtained from its imaginary part by means of the subtracted DR m representing the mass of the lepton circulating in the loop. Since the q δ q λ term in eq. (3.1) does not contribute because of gauge invariance, the two-loop calculation can be performed starting from the one-loop amplitude with the replacement The above approach has been used in ref. [57][58][59][60] for the calculation of two-loop QED hadronic and leptonic corrections to Bhabha scattering. More recently the same technique has been used in refs. [53,54] to calculate the NNLO hadronic corrections to µe scattering in the context of the MUonE experiment. Partial preliminary results on NNLO virtual leptonic corrections to muon-electron scattering have been obtained in ref. [71]. Equation (3.5) leads in a natural way to perform a MC integration in dz, for the evaluation of NNLO virtual pair corrections, starting from the NLO calculation. Indeed, by exchanging the order of integration variables between the loop momentum and the effective photon mass √ z given by the DR, the one-loop amplitude A 1L (µ ± e − → µ ± e − ) can be convoluted with a kernel function to obtain the two-loop amplitude with the VP insertion in the following way where A 1L (µ ± e − → µ ± e − ; z) represents the one-loop amplitude with the mass of the virtual photon, where the VP is inserted, set to √ z.
In this way the CPU load for the two-loop calculation remains of the same order of the one-loop calculation, having only one additional dimension in the MC integration. An important point to be checked is the numerical stability of the algorithm. For this reason we performed several internal checks on the convergence of the MC error estimate and cross checks with available analytic expressions, which we document below.

Calculation of real corrections
The calculation of the differential and integrated cross sections for real pair emission, i.e. for the processes µ ± e − → µ ± e − + − , with = e, µ (the case = τ is not allowed at MUonE by energy conservation) is a standard tree-level calculation for 2 → 4 processes. Because of the presence of the small electron mass scale and the need of covering the full available phase space in the simulation, the matrix elements have to be computed in a numerically stable way and the numerical integration has to be performed with dedicated samplings, in particular for the final state with electronic pair. As will be discussed later, also the scattering amplitudes for the processes µ ± e − → µ ± e − τ + τ − with m τ = m e has been calculated.
The matrix elements (considering only QED interactions) have been calculated with the help of the symbolic manipulation program Form [72][73][74] and cross-checked with the automatic package Recola [75,76], finding perfect agreement (typically the numerical relative difference is at the 10 −12 level and never worse than 10 −8 ) 5 .
The phase space integration has been carried out with two independent implementations. For both of them the final state momenta are generated in the centre of mass frame and boosted in the laboratory frame. According to one implementation, the 4-body Lorentz invariant phase space In order to select subsets of gauge invariant classes of diagrams, for the calculation of the amplitude µ ± e − → µ ± e − τ + τ − with mτ = me, mµ with Recola, we used the rescaling factors λµ and λτ for the µµγ and τ τ γ couplings, respectively.
is decomposed into different parameterisations. The diagrams (a) and (b) of figure 8, with the lepton pair radiated from the electron line, require the following parameterisation: where P = p 1 + p 2 is the total initial state momentum (with p 1 and p 2 being the initial muon and electron momenta, respectively), dΦ n is the n-body phase space element and p 3 , p 4 , p 5 and p 6 are the final state muon, electron and + , − momenta, respectively. In order to obtain reliable numerical MC predictions, different sets of independent variable mappings are necessary to tailor the generation according to the peaking behaviour of diagrams (a) and (b) of figure 8. In particular, for diagram (a) the following set of independent variables is used: where θ * 3 and φ * 3 are generated in the rest-frame of the pair 56 ( q 5 + q 6 = 0), and θ 56 and φ 56 are the polar and azimuthal angles of the pair momentum (q 5 + q 6 ) in the centre of mass frame. For diagram (b), θ 56 and φ 56 are understood in the centre of mass frame with z-axis oriented along the outgoing electron p 4 . The diagrams (c) and (d) of figure 8, with the lepton pair radiated from the muon line, do not require dedicated parameterisation because m µ is of the same order of magnitude of √ s and therefore the logarithmic structure developed by the numerical integration is not leading. The peripheral diagrams (e) and (f) of figure 8 require a different phase space parameterisation [77,78]: The independent variables are the invariant masses and the angles where cos θ * 3 , φ * 3 are generated in the Q 356 rest-frame and cos θ * * 5 and φ * * 5 are generated in the Q 56 rest-frame.
When the + − pair is an electron-positron pair e + e − , all the necessary simmetrised channels with p 4 and p 6 exchanged are considered.
A second (and independent) phase-space decomposition has been implemented, which is the one used in the official release of Mesmer. In this case, in the principal integration channel, the independent variables are where Q 2 456 = (p 4 +p 5 +p 6 ) 2 , Q 2 56 = (p 5 +p 6 ) 2 , t 13 = (p 1 −p 3 ) 2 , φ 3 is the p 3 azimuthal angle in the p 1 + p 2 rest frame, θ † 56 and φ † 56 are the p 5 + p 6 angles in the p 1 − p 3 + p 2 rest frame, and θ 5 and φ 5 are the p 5 angles in the p 5 + p 6 rest frame. The importance sampling of the variables in eq. 4.8 follows as close as possible the peaking structure of the differential cross section, and multi-channel techniques are used to map different structures. For instance, in the case of an extra electronic pair and referring to figure 8, in the principal channel, while φ 3 , φ † 56 and φ 5 are sampled uniformly, Q 2 456 is sampled as 1/Q 2 456 to follow the singularity due to the internal fermion propagator of diagram (b), t 13 is sampled according to 1/t 13 to flatten the peaks of the photon propagator exchanged between the muon and electron lines and Q 2 56 is sampled following 1/Q 2 56 . The cosine cos θ † 56 is distributed according to in the p 1 − p 3 + p 2 rest frame, to flatten the peaks due to the internal electron propagator of diagram (a) or the photon propagator in diagram (e). Finally, cos θ 5 is sampled according in the p 5 + p 6 rest frame, to flatten the peak due to the electron propagator in diagrams (e) and (f). In the other channel, the roles of p 1 and p 3 are interchanged with p 2 and p 4 . For identical particles, a symmetrisation of the channels is performed by exchanging p 4 ↔ p 6 .

Numerical results
We are going to study the numerical impact of QED NNLO leptonic pair corrections to observables for virtual and real corrections separately, in order to understand the features and the size of each contribution.
The simulations are performed assuming a 150 GeV muon beam 6 and according to the following event selection criteria 7 : • θ e , θ µ < 100 mrad and E e > 1 GeV. The angular cuts model the typical acceptance conditions of the experiment and the electron energy threshold is chosen in order to single out a region where the signal of the experiment is not negligible.
We remark that this event selection is driven by the 2 → 2 kinematics and can be rephrased differently to describe the events more in detail. In the 2 → 2 kinematics, the outgoing muon is characterised by an angle θ µ <θ µ = arccos 1 − m 2 e /m 2 µ 4.84 mrad and, for a 150 GeV incident muon, an energy E µ >Ē µ 10.28 GeV 8 . When we consider real lepton pair radiation we can have more charged tracks in the detector 9 and therefore we need to introduce a more general event selection inspired by the one of the elastic 2 → 2 process. In particular, we can define a muon-like track, characterised by angle lower than 6 Lately, it became more likely that MUonE running conditions will use a 160 GeV muon beam. Here we still use a 150 GeV energy for direct comparison with past phenomenological results and because it does not change our conclusions.
7 It corresponds to Setup 1 in ref. [49] and to Setup 2 in ref. [44]. In the latter, also setups with Ee > 0.2 GeV were considered, which are disregarded here for the sake of brevity. 8 It corresponds to the laboratory energy of a muon scattered backward in the centre-of-mass frame. 9 In its present proposal, MUonE will not distinguish positive from negative charges. θ µ and energy larger thanĒ µ , and an electron-like track, with angle lower than 100 mrad and energy larger than 1 GeV. The generalised selection criteria can be defined as follows: • two tracks in the detector, one muon-like (θ µ <θ µ 4.84 mrad and E µ >Ē µ 10.28 GeV) and one electron-like (θ e < 100 mrad and E e > 1 GeV) This alternative definition will be used in the following when 4-body final states are studied, and we will name it basic acceptance cuts.
On top of it, we are going to consider also three extra selection cuts (and combinations of them), which are tailored to select elastic events in order to suppress reducible backgrounds and keep the sensitivity to the signal. The elasticity cuts we implement are: 1. cut 1 − an additional cut on the minimum electron-and muon-like scattering angles, θ e , θ µ > θ c , to reduce the impact of peripheral diagrams in events with additional electron pair emission. For illustrative purposes, we use the value θ c = 0.2 mrad, which appears to be a good compromise between the requirement of keeping sensitivity to ∆α(t) in the large |t| region 10 and effectively suppressing the reducible background; 2. cut 2 − acoplanarity cut ξ = |π − |φ e − φ µ || < ξ c = 3.5 mrad 11 , where φ e,µ are the azimuthal angles of the electron-and muon-like momenta; 3. cut 3 − elasticity distance δ < δ c = 0.2 mrad, where δ is defined as the minimal Cartesian distance of the (θ e , θ µ ) point in the θ e -θ µ plane from the elastic curve. The elastic curve can be parameterised by the function θ µ (θ e ) = arctan 2m e r cos θ e sin θ e E i µ − r rE i µ + 2m e cos 2 θ e , where r is defined in eq. (9) of ref. [35] and E i µ is the incident muon energy in the laboratory reference frame. The distance of the generic point (θ 0 e , θ 0 µ ) is defined as While the optimal value for δ c can be determined only with a realistic simulation that takes into account of all the experimental conditions and detector effects, we use for the present study the reference value δ c = 0.2 mrad.
When considering the emission of real leptonic pairs in the following, in order to study the effects from different gauge invariant contributions and their interplay with the virtual corrections from a purely theoretical point of view, neglecting the effect of identical particles, we will begin our analysis considering the processes µ ± e − → µ ± e − τ + τ − for the two cases m τ = m e and m τ = m µ . The three gauge invariant contributions are the radiation from the electron line, from the muon line and the two peripheral diagrams, according to the discussion of section 2 on the ingredients of dσ α 2 real . In this exercise we are completely inclusive over the additional "τ " pair and assume the cuts are applied to the final state electron and muon, as if they were distinguishable from the extra leptons in the final state (asymmetric cuts). The obtained numerical results must be considered "academic" since the event selection does not fully match the realistic MUonE experimental settings.
As a second step, we consider the true processes µ ± e − → µ ± e − e + e − and µ ± e − → µ ± e − µ + µ − , taking into account more realistic conditions. In MUonE the final state charged particles are detected through their tracks in the silicon detectors, without charge identification. We therefore require that exactly two tracks are reconstructed in the detector by asking that they have a threshold energy larger than E thr. = 0.2 GeV and an angle lower than θ thr. = 100 mrad. If more than two tracks are reconstructed, the event is rejected. Otherwise, of the two selected tracks, we require one to be muon-like and the other to be electron-like, as per definitions above. We will further show the effects of cutting through If p 1(2) is the incoming muon (electron) momentum and p 3(4) is the muon(electron)-like outgoing momentum 12 , the considered differential observables are the following ones: where t ee = (p 2 −p 4 ) 2 , t µµ = (p 1 −p 3 ) 2 , θ e and θ µ are the scattering angles, in the laboratory frame, of the outgoing electron-and muon-like tracks, respectively. We decide to represent the corrections as a differential K factor defined by i.e. we normalise the differential cross sections under study w.r.t. the LO ones. The choice is dictated by consistency and for direct comparison with results presented in refs. [44,49] and because the 10ppm accuracy goal is implicitly referred to LO results. We notice that in the following we show also the hadronic effects already considered in [54], for completeness and cross-check purposes.

NNLO virtual pair corrections: the factorised diagrams
In this section we show the impact of NNLO virtual pair corrections of figures 1 and 2 13 , which are factorised over the tree-level cross section and contribute to the running of the electromagnetic coupling constant. In figure 9, the red histograms show the contribution of 0  The same breakdown of the separate effects, displayed in incremental way, is shown for the t ee distribution in figure 10. For comparison purposes, in the inset we present the function 10 5 × ∆α NLO (t) for an electron or muon loop as coded in the KNT v3.0.1 package and as obtained through a DR, showing perfect agreement. The DR we use here is fully analogous to eq. (3.2) where the function R (z) is replaced by the NLO function R NLO (z) resulting from the QED limit of eq. (18) in ref. [79].

Vacuum polarisation on photonic vertex corrections
Among the two-loop irreducible contributions where the one-loop diagrams with photonic corrections include a VP insertion in the loop photon propagator, the diagrams of figure 5, together with their respective counterterms, are gauge-invariant and IR safe. Therefore   it is meaningful to investigate their numerical relevance. In figures 11 and 12 we show, separately, all the effects of different leptons (electron and muon) in the photonic vertex one-loop diagrams of figure 5 for the three observables θ e , θ µ , t ee . For each contribution the histogram refers to the calculation through DR and the points refer to the analytical calculation. For the case of a closed loop with a lepton with mass equal to the one of the current which the photon loop is attached to, the analytical calculation presented in Appendix A.2 of ref. [80] has been used. In all other cases (i.e. mass of the lepton flowing in the VP insertion different from the mass of the lepton of the external current), we use an in-house independent calculation which is described in Appendix A.
We would like to remark that for all cases excellent numerical agreement is found, which indicates also that the DR approach is implemented correctly.
From the phenomenological point of view, the largest contribution is given, as expected, by the electron loop on the electron line (red curves), which reaches the level of −0.037% at the kinematical boundaries corresponding to the maximum |t ee |. As known in the literature, see e.g. refs. [64,80,81], the leading logarithmic term of the electron loop in the vertex correction on the electron line is proportional to α 2 log 3 [m 2 e /(−t ee )] which is expected to be partially cancelled by the corrections corresponding to the real emission of an e + e − pair. This will be discussed in the following section with a numerical investigation.
The electron loop on the muon line is suppressed w.r.t. the previous contribution and reaches the level of −0.014% at the kinematical boundary. The smallest contribution is given by the muon loop on the muon line, which never reaches the 0.001% level.
Finally, the muon loop correction on the electron line is of the same size of the hadronic correction on the electron line, remaining at most at the level of a few ppm.   Figure 12: The same as figure 11 for the t ee distribution.

NNLO virtual pair corrections: the complete set of diagrams
In this subsection we present the numerical impact of the contributions due to the photonic one-loop amplitudes (virtual and real contributions) with an additional closed leptonic loop, namely the sum of the interference between the diagram in figure 1 (a) with the diagrams in figures 4, 5 and 6 plus the interference of the diagram in figure 1 (b) with the diagrams in figure 3. Also real radiation is added (i.e. the interference of the diagrams in figure 7 with the very same diagrams without loop insertion), as well as all the contributions described in sections 5.1 and 5.2. As we know from previous studies on NLO [44] and NNLO photonic [49] (and also hadronic [54]) corrections, the real radiation can give large contributions, in particular for the electron scattering angle distribution. For this reason in refs. [44,49] the additional acoplanarity cut (cut 2 as defined in the introduction to this section) was introduced, in order to partially remove hard radiation effects. e blob (e radiation) e + µ blob (e radiation) e + µ + had blob (e radiation) e blob (full radiation µ − ) e + µ blob (full radiation µ − ) e + µ + had blob (full radiation µ − ) e blob (full radiation µ + ) e + µ blob (full radiation µ + ) e + µ + had blob (full radiation µ + ) with acoplanarity cut (3.5 mrad) e blob (e radiation) e + µ blob (e radiation) e + µ + had blob (e radiation) e blob (full radiation µ − ) e + µ blob (full radiation µ − ) e + µ + had blob (full radiation µ − ) e blob (full radiation µ + ) e + µ blob (full radiation µ + ) e + µ + had blob (full radiation µ + ) Figure 13: Left: differential K NNLO factor on θ e distribution, including the complete set of NNLO virtual leptonic pair corrections. Right: the same observable including the acoplanarity cut. with acoplanarity cut (3.5 mrad) e blob (e radiation) e + µ blob (e radiation) e + µ + had blob (e radiation) e blob (full radiation µ − ) e + µ blob (full radiation µ − ) e + µ + had blob (full radiation µ − ) e blob (full radiation µ + ) e + µ blob (full radiation µ + ) e + µ + had blob (full radiation µ + ) Figure 14: Left: differential K NNLO factor on θ µ distribution, including the complete set of NNLO virtual leptonic pair corrections. Right: the same observable including the acoplanarity cut. e blob (e radiation) e + µ blob (e radiation) e + µ + had blob (e radiation) e blob (full radiation µ − ) e + µ blob (full radiation µ − ) e + µ + had blob (full radiation µ − ) e blob (full radiation µ + ) e + µ blob (full radiation µ + ) e + µ + had blob (full radiation µ + ) −120 e blob (e radiation) e + µ blob (e radiation) e + µ + had blob (e radiation) e blob (full radiation µ − ) e + µ blob (full radiation µ − ) e + µ + had blob (full radiation µ − ) e blob (full radiation µ + ) e + µ blob (full radiation µ + ) e + µ + had blob (full radiation µ + ) Figure 15: Left: differential K NNLO factor on t ee distribution, including the complete set of NNLO virtual leptonic pair corrections. Right: the same observable including the acoplanarity cut. Figure 13 shows the impact of the virtual pair corrections on the electron scattering angle distribution for different levels of approximation: radiation from the electron leg only or the complete calculation for µ − or µ + beam. For each of these cases, three different contributions are considered: only electron loop, electron plus muon loop and electron plus muon plus hadronic loop. On the left, in the absence of any acoplanarity cut, all the above described contributions are almost degenerate, with very large corrections up to about 2.5% for small scattering angles, driven by the electron loop on top of radiation from the electron line. On the right, the additional acoplanarity cut changes the sign of the corrections, which range from about −0.1% for large scattering angles to several −0.1% for θ e → 0. In particular, the contribution from the electron line radiation reaches the value of about −0.4%, driven by the electron closed loop. The muon and the hadronic loops are of the same order of a few 0.01%. The full corrections (solid lines for µ − beam and dot-dashed line for µ + beam) display a difference with respect to the radiation from the electron line at the few 0.1% level for small scattering angles. Figure 14 shows the impact of the virtual pair corrections on the muon angle distribution 15 . In this case, the corrections grow to dozens of 0.01% as θ µ reaches its maximum. We notice a flattening of the K NNLO factor when the acoplanarity cut is added (right). The effect is generally dominated by electron loop insertion, and in the full calculation corrections for an incident µ + or µ − are opposite with respect to corrections on the electron line only. Finally, figure 15 shows the impact on the momentum transfer t ee , where similar considerations apply.
We conclude the section by remarking that the application of extra elasticity cuts, in particular cut 1+3, does not change the overall effects and size of the corrections considered here. In contrast, they have a much more dramatic effect in reducing the reducible background due to real electron pair emission, as will be discussed in the following section 5.5.

NNLO real pair corrections, asymmetric event selection (academic case)
In this section we present the impact of the emission of a real leptonic pair, treated in an inclusive way (i.e. with asymmetric cuts). To exclude diagrams with the exchange of identical particles, we consider in this section the processes µ ± e − → µ ± e − τ + τ − , setting m τ = m e , and adopt the generalised event selection criteria described in the introduction of this section. At the end of this subsection, for completeness, we consider also the negligible contribution due to the emission of an extra muon pair (m τ = m µ ).
In the left panel of figure 16 we show the contribution of the "τ " pair (with m τ = m e ) on the t ee distribution. Three different gauge invariant classes of diagrams are considered: radiation from the electron line only (red line), radiation from the muon line only (lightblue line), radiation from both electron and muon lines (see figure 8 (a)-(d)), including their interferences but neglecting the peripheral diagrams (figure 8 (e) and (f)). For the latter two options, the blue and violet crosses refer to the µ + and µ − beam cases, respectively. The results obtained with the complete matrix elements, including also the peripheral diagrams, are displayed with the blue and violet line for the µ + and µ − beam, respectively. As a general comment we can see that the numerical impact of the real e + e − radiation, on the t ee distribution, is at the level of few 0.01%, with dominance of the radiation from the electron line w.r.t. the radiation from the muon line. The interference between the radiation from the electron and muon lines is of different sign when flipping the charge of the incoming muon, as already seen for the QED NLO corrections in ref. [44]. The effect of the peripheral diagrams is almost negligible.  no identical particles, asymmetric cuts Figure 16: Real e + e − radiation effects with academic event selection on t ee (left) and θ e (right) distribution.  no identical particles, asymmetric cuts ξ < 3.5 mrad Figure 17: Real e + e − radiation effects with academic event selection on t ee (left) and θ e (right) distribution, in the presence of the acoplanarity cut.
The right panel of figure 16 shows the electron scattering angle distribution, normalised to the LO prediction, where the large correction factor for the radiation from the electron line is clearly visible for θ e → 0, analogously to what happens for the NLO QED correction. As commented in ref. [44], the effect is in part due to the vanishing of the LO differential cross section as θ e → 0. For this reason, in ref. [44] the acoplanarity cut has been introduced, which allows to partially remove the enhancements. Applying the same acoplanarity cut on our signature, for the t ee observable we obtain the results displayed in the left panel of figure 17, where the maximum of the correction is reduced to about 0.033%, for the case of an incoming µ + beam. The corresponding corrections on the θ e observable are shown in the left panel of figure 17, and also here we remark the large suppression factor induced by the acoplanarity cut.
Within this case, where we academically suppose one final state electron to be distinguishable from those of the extra pairs, the contribution from the peripheral diagrams is almost negligible for all the considered observables.
We also notice that on the K NNLO factor for the t ee distribution, in the most inclusive case, without the acoplanarity cut, corrections due to the emission of an extra e + e − pair are positive and of the same order of magnitude as the negative corrections due to VP insertion on photonic vertex corrections, discussed in section 5.2. We see indeed a partial cancellation between real and virtual pair corrections, for this subset of contributions.
Before considering also identical particles and fully realistic event selection on real pair emission, we show in this academic case the impact of µ ± e − → µ ± e − µ + µ − events (i.e. µ ± e − → µ ± e − τ + τ − with m τ = m µ ) on the t ee distribution, for completeness. It is clear from figure 18, where the red curve shows the incoming µ + case and the blue one the incoming µ − case, that this process gives a totally negligible contribution, reaching at maximum a few units in 10 −8 w.r.t. the LO differential cross section. We also notice that below t ee −0.06 GeV 2 there are no events because of phase space constraints 16 .  Figure 18: Impact of the µ ± e − → µ ± e − µ + µ − process on the t ee distribution, for the academic selection cut.

NNLO real pair corrections, realistic event selection
The discussion on real electron pair radiation of section 5.4 shows the relevance of real pair emission in the NNLO complete calculation. However, those results are purely academic since the analysis is completely inclusive on the emitted electron pair and does not take into account the MUonE realistic event selection, where the four final state leptons can produce up to four tracks in the detector. According to the discussion of section 5, to mimic an elastic event, we can require that only two tracks (one muon-like, one electronlike) are seen in the detector. In this way, without charge identification, in principle every pair out of the six possible pairings can hit the detector while the remaining two particles are below energy threshold and/or out of the angular allowed range. This requires the use of the full matrix elements for the calculation of the cross sections of the processes µ ± e − → µ ± e − e + e − (with MUonE realistic event selection we are considering here, no events for the process µ ± e − → µ ± e − µ + µ − satisfy the criteria, because at least three charged visible particles are always present), including all the diagrams with the two indistinguishable electrons interchanged, which were neglected in the previous subsection, and the adoption of symmetric cuts that do not distinguish identical particles in the final state. This introduces a new important feature: there is no distinction between the electron of the underlying tree-level and the electron or positron of the emitted pair. As a consequence, in the photon propagators the exchanged momentum can become very close to zero, giving a dramatic enhancement in the cross section and in the distributions, especially at the corner of the phase space. In this enhancement the peripheral diagrams (figures 8 (e)-(f)) give by far the largest contribution. For instance, a sharp peak in a very small window of the scattering angles can be seen in figure 19, where it is evident that "switching-on" peripheral suppressed due to the further reduction of available phase space because of the larger pion mass. diagrams produces a huge effect.   Figure 19: Differential K NNLO factor for θ e (left) and θ µ (right) distributions for real e + e − radiation with realistic event selection, without acoplanarity cut. θ µ is plotted above 0.2 mrad for readability.
The kinematical configurations related to these enhancements correspond typically to small momentum transfer for t ee or t µµ (figure 20) and could spoil the sensitivity of MUonE to the hadronic component of the QED running coupling constant by distorting in a significant way the differential cross sections where the signal to be measured is more important. Even with the inclusion of the acoplanarity cut, the numerical impact of real pair radiation   Figure 20: Left: t ee differential K NNLO factor for real e + e − radiation with realistic event selection, without acoplanarity cut. Right: the same for t µµ .
remains at the level of 0.1%, with a large contribution for θ e → 0 reaching the 10% level, as can be seen in figure 21.
In order to reduce the impact of real pair radiation to a manageable level, we propose to further filter events through the cuts labelled as cut 1-3 and described in the introduction  with identical particles, symmetric cuts ξ < 3.5 mrad Figure 21: Left: Differential K NNLO factor on t ee distribution for real e + e − radiation with realistic event selection, with acoplanarity cut. Right: the same for θ e .
to section 5. The effectiveness of the additional cuts in suppressing real electron pair radiation events can be visualised with the scatter plot of the θ e -θ µ angles correlation displayed in figure 22. The upper panel shows, with red points, the scatter distribution of Figure 22: Scatter plot of (θ e , θ µ ) points for 5 · 10 5 µ − e − → µ − e − e + e − simulated events.
the coordinates (θ e , θ µ ) for a sample of 5 · 10 5 events for the process µ − e − → µ − e − e + e − with only one muon-like and one electron-like track. It is clearly seen that many points fall well outside the elastic correlation curve (black line), in particular they gather in the region of small θ e and θ µ . The blue crosses display the same sample after applying the additional cut θ e , θ µ > θ c (cut 1). It is clearly visible the removal of the events with θ µ close to 0. The background rejection efficiency of such a cut is of about 99.5%. In the second panel from the top, the effect of adding the acoplanarity cut (cut 1+2) is shown and the rejection efficiency increases to about 99.96%. The third panel from the top displays the effect of imposing the cut on the distance from the elasticity correlation curve, δ < δ c , on top of cut 1. In this case the rejection efficiency is at the level of 99.89%. Finally, the bottom panel shows the effects of combining the three cuts (cut 1+2+3). As can be seen, this combination is particularly effective in reducing the real electron pair production, since only a tiny fraction (as small as 0.007%) of background events passes through the selection. It is interesting to observe the complementarity between the ξ and δ cuts. The latter alone suppresses hard radiation contributions while allowing for a certain degree of acoplanarity between the azimuthal angles of the two final state tracks, due to soft transverse radiation. This analysis allows us to propose an improved event selection, cut 1-3, as defined in section 5 for the analysis of MUonE data. As can be seen in figures 23 and 24, with this event selection the real e + e − emission contribution gets reduced to the ppm level over the full kinematical range. It is worth noticing, for instance in the θ e K NNLO factor, a peak on the leftmost bin: it can be traced back to few events that, after applying cut 1+2, lie in the region θ e < δ c allowed by the elasticity cut.
It is understood that more complete simulations will be needed to select the optimum values for the cut parameters θ c , δ c and ξ c . Here, we have only demonstrated that by a sensible choice of the cut parameters the potentially large background due to electron-pair emission can be reduced at the 10 −5 level.   Figure 24: Left: Differential K NNLO factor on t ee distribution for real e + e − radiation with realistic event selection, with all cuts applied. Right: the same for t µµ .

Summary and prospects
In this work we have discussed the calculation of the NNLO leptonic corrections to µ ± e − → µ ± e − processes and their implementation into the Monte Carlo code Mesmer, used for simulations for the MUonE experimental proposal. Such a project aims at a high precision measurement of the effective QED coupling constant in the space-like region, which will allow for a determination of the leading order hadronic contribution to the muon anomaly, in a way completely independent from time-like dispersion relation and Lattice QCD approaches.
The NNLO leptonic corrections consist of virtual and real contributions. The former have been calculated with dispersion relation techniques including all finite mass effects. For the class of vertex corrections, the numerical results have been cross-checked with exact analytical formulae, valid for arbitrary internal and external masses, finding excellent agreement. Also the hadronic contributions have been included by means of two different parameterisations existing in the literature. The real corrections have been calculated by means of exact matrix elements, taking into account all finite mass effects.
On the phenomenological side, we analysed all different, gauge-invariant, classes of virtual contributions, investigating their numerical impact on the differential observables relevant for MUonE realistic running conditions. As expected, the leading contributions come from the closed electron loop diagrams and the overall size of the factorised and vertex contributions is of the order of few 0.01% (positive the former, negative the latter). The vacuum polarisation insertions on the NLO QED photonic corrections can reach the level of few 0.1%, with shapes driven by the NLO photonic corrections. The corrections are negative for all distributions, except for the electron scattering angle distribution, where they become positive, reaching the few % level for θ e → 0. Similarly to what happens to the large QED NLO correction to dσ/dθ e , the introduction of the acoplanarity cut mitigates the effects and brings the corrections to negative values. On the other hand the acoplanarity cut enhances the size of the corrections of IR origin.
The hadronic contributions are of the same order of magnitude of the muon loop insertions, as already noticed in previous studies on Bhabbha scattering at flavour factories. Particular attention has been paid to the real leptonic pair emission. Because of the low centre of mass energy combined with realistic event selection criteria, the µ + µ − pairs do not play almost any role at MUonE. On the other hand, the emission of e + e − pairs is particularly important. In general the processes µ ± e − → µ ± e − e + e − are distinguishable from the elastic µ ± e − → µ ± e − processes. However, under realistic MUonE event selections, the 2 → 4 processes can produce only two charged tracks in the final state and become, therefore, degenerate with the elastic signal events. In this case the real pair emission cross section has to be summed to the virtual pair emission contribution, to reach a reliable theoretical prediction. Because of the presence of two indistinguishable electrons and peripheral diagrams, we observed very large positive contributions for small θ e and θ µ as well as for small t ee and t µµ values, which are not tamed by the previously introduced acoplanarity cut. In order to control these effects we studied a combination of additional cuts, namely a small minimum scattering angle for the two observed tracks and a maximum distance from the elasticity correlation curve θ µ (θ e ). After considering all the above cuts, the size of the real pair radiation effects has been shown to remain below the 10 −5 level.
Therefore, the combined effects of virtual and real pair radiation under the simplified event selection we considered are dominated by virtual contributions, in particular the ones stemming from electron loops, which, depending on the observable, can reach the 1% level in units of the LO differential cross sections. It is true that there is a partial cancellation between a subset of the virtual and a subset of the real contributions, but the presence of peripheral diagrams in real pair emission spoils the expected cancellation and forces the introduction of elasticity cuts to keep under control a potentially large background to the signal.
The work presented in this paper represents an additional step towards the implementation of a fully fledged MC generator including the complete set of NNLO QED corrections matched to multiple photon emission, which will be ultimately needed for the analysis of MUonE data.  [82,83] that can be used. However, to the best of our knowledge, the complete analytic expression for the case m 1 = m 2 without any expansion in the mass ratio is not yet available. In this section we show a way to calculate these vertex form factors for m 1 = m 2 .
The kinematics of the process is given by The integrals that appear in the calculation are of the form I(n 1 , n 2 , · · · , n 7 ) ≡ C( ) d d k 1 d d k 2 1 D n 1 1 D n 2 2 · · · D n 7 The propagators are defined as: Seven master integrals, depicted in figure 26, are necessary to calculate the form factors. We solved these particular integrals using differential equation methods [84][85][86]: we chose a suitable basis of master integrals, such that the dimensional regularisation parameter factorises from the kinematics. The integrals were then encoded in a d log-form, also called canonical form [87], which we were able to obtain using the Magnus algorithm [88,89]. The kinematic variables were chosen in such a way that the arguments of the d log's were simple rational functions, which enabled us to express the integrals in terms of the generalised polylogarithms [90][91][92]. The master integrals were then calculated up to required expansion to determine the form factor: they are I(0, 2, 2, 0, 0, 0, 0), I(0, 2, 0, 2, 0, 0, 0), I(0, 2, 2, 1, 0, 0, 0), I(0, 2, 1, 2, 0, 0, 0), I(0, 2, 0, 2, 1, 0, 0), I(0, 1, 2, 1, 1, 0, 0) and I(0, 1, 2, 2, 1, 0, 0), which can be found in the ancillary file in the arXiv version of the paper. Here we explicitly write the first four master integrals up to finite term to show the used convention. The following equations show the solutions that are suitable for the region of interest: Throughout the form factor calculation, the Integration by Parts reduction [93][94][95] was carried out using the publicly available codes Reduze [96] and FIRE [97]. The differential equations were generated using Reduze and, to rationalise the arguments of the d log, we used the package RationalizeRoots [98]. The numerical validations were accomplished using SecDec [99] and the numerical evaluation of generalised polylogs was performed with GiNaC [92] and handyG [100]. A detailed result of the vertex form factor is presented in ref. [101].
The resulting expressions for the form factors have been used to cross-check with high accuracy the DR results of figures 11 and 12.