Theory for muon-electron scattering @ 10ppm: A report of the MUonE theory initiative

We review the current status of the theory predictions for elastic $\mu$-$e$ scattering, describing the recent activities and future plans of the theory initiative related to the proposed MUonE experiment.


Introduction
There is renewed interest in obtaining precise theoretical predictions for elastic muon-electron scattering. This is to be seen in the context of MUonE [1], a recent proposal to perform a very precise measurement of µ-e scattering [2]. A comparison of experimental data with perturbative calculations can be used to extract the hadronic vacuum polarisation (HVP) through its contribution to the running of the QED coupling α. This follows the original idea of using scattering data to extract the leading hadronic contribution to the muon (g − 2) from the effective electromagnetic coupling in the space-like region [3]. The measurement of the running of alpha in the space-like region from small-angle Bhabha scattering was proposed in [4] and done in [5].
The effect of the HVP changes the differential cross section of µ-e scattering by up to O(10 −3 ), depending on the The proposal of MUonE is to scatter a 150 GeV muon beam on a Beryllium fixed target. In order to obtain sufficient statistics and reduce multiple-scattering effects [6], the target (about 60 cm in total) is split into many (about 40) thin layers. The measurements are done in several stand-alone stations of about 1 m length and 10 × 10 cm 2 transverse dimension. The scattering angles of the electron θ e and the muon θ µ (in the lab frame) are measured very precisely, but no further kinematic information is assumed to be available.
From an idealised point of view we thus consider where the initial-state electron is at rest and X stands for any further radiation. With the energy of the incoming muon set to E 1 = 150 GeV, the centre-of-mass energy is fixed as s = m 2 + M 2 + 2mE 1 ≃ (400 MeV) 2 , where m and M denote the electron and muon mass, respectively. The momentum transfer t ranges from t min ≃ −(380 MeV) 2 to zero. Hence, there are two widely different scales entering the process with m 2 ≪ Q 2 , where Q 2 stands for the large scales M 2 ∼ s ∼ |t|.
The resulting large logarithms ln(m 2 /Q 2 ) will have to be properly accounted for in the theoretical treatment of the process.
The incoming muon beam consists of either positively or negatively charged muons and is about 80% polarised. Since the electrons in the target are unpolarised and QED is parity conserving, the only effect of the polarisation is due to the electroweak contributions coming from the Z-boson exchange. At tree level, the latter contributes at the level of 10 ppm and, hence, has to be included.
There are several effects that result in differences from the idealised process. First, the electrons are bound, and the impact of bound-state effects should be estimated. Second, there are nuclear background processes due to µ-N scattering. From our point of view, however, the most important aspect of (1) is the selection of elastic scattering. Since photons are not detected, there is no way of telling how much radiation is present in X. A contribution including n photons results in a suppression by α n relative to the leading order (LO), i.e. a N n LO contribution. Another relevant process is open lepton-pair production, i.e. X = e + e − or, albeit with a very small phase space, X = µ + µ − . This amounts to a NNLO QED contribution. Finally, there is also a background from pion production where the pion subsequently decays into two photons, i.e. X = π 0 → γ γ. As a last option we mention X = π + π − , again with a very small phase space.
In the absence of additional emission X in the final state, i.e. for the elastic 2 → 2 scattering process, we can derive a simple functional relation between θ e and θ µ that we call the elasticity curve. Thus, allowing only events within a small band around this curve effectively selects nearly elastic events. However, from a theoretical point of view this is problematic. Making a stringent cut on the phase space is a further source of large logarithms, beyond the ln(m 2 /Q 2 ) mentioned above, that might need to be resummed.
The precision expected at the MUonE experiment also raises the question whether possible new physics (NP) could affect its measurements. This issue was addressed in [7], studying possible NP signals in muon-electron collisions at MUonE due to heavy or light mediators, depending on whether their mass is higher or lower than O(1 GeV). The former were analysed in a model-independent way via an effective field theory approach, whereas for the latter the authors discussed scenarios with light spin-0 and spin-1 bosons. Using existing experimental bounds, they showed that possible NP effects in muon-electron collisions are expected to lie below MUonE's sensitivity, therefore concluding that it is very unlikely that NP contributions will contaminate MUonE's extraction of the HVP. The authors of [8] addressed the sensitivity of MUonE to new light scalar or vector mediators able to explain the muon g − 2 discrepancy. They concluded that the measurement of the HVP at MUonE is not vulnerable to these NP scenarios. Therefore, the analyses of [7] and [8] reach similar conclusions where they overlap. These results confirm and reinforce the physics case of the MUonE proposal.
In what follows we will discuss all issues related to obtaining a theoretical prediction for µ-e scattering at 10 ppm. We start in Section 2 by briefly revisiting the kinematics of µ-e scattering. Next, we discuss in Section 3 the fixed-order perturbative calculations in QED. This is followed by a discussion on how to include the HVP in Section 4. Possible strategies on how to deal with and estimate the importance of contributions beyond those included in the fixed-order calculations are considered in Section 5. Finally, in Section 6 we give an outlook on how the various pieces can be combined into a general purpose Monte Carlo code that provides a sufficiently accurate theoretical prediction, before we present our summary in Section 7.

Kinematics of µ-e scattering
Let us begin by reviewing for the elastic µ-e scattering process, the basic relations between angles, energies and momenta in the laboratory frame (LAB) and in the centre-of-mass system (CMS). In a fixed-target experiment, where the electron is initially at rest, the Mandelstam variables s and t are given by Here, E 1 is the energy of the incident muon, E 4 is the electron recoil energy and is the Källén function. The third Mandelstam variable u is related to s and t in the usual way as s + t + u = 2M 2 + 2m 2 . It is also convenient to define the variable x that is related to t as With t min ≃ −(380 MeV) 2 the range of x is 0 ≤ x 0.933 and x = 0 corresponds to t = 0.
The parameters for the Lorentz transformation between the LAB and the CMS are We define the scattering angles θ e,µ in the LAB and θ * e,µ in the CMS as the angles between the direction of the incident muon and the outgoing electron or muon. While in the CMS we trivially have θ * e = π − θ * µ , in the LAB frame the two angles are correlated by the elasticity condition where and β * µ is the muon velocity in the CMS. In the θ e -θ µ plane, (7) defines the elasticity curve depicted in Figure 1. This is the fundamental constraint for MUonE to discriminate elastic scattering events from the background of radiative events and inelastic processes. Since g * µ > 1, the outgoing muon is always emitted in the LAB forward direction at an angle smaller than θ max µ = 4.8 mrad (for E 1 = 150 GeV), where On the contrary, the recoiling electron can be emitted in the whole LAB forward hemisphere, i.e. 0 ≤ θ e ≤ π/2, since g * e = β/β * e = 1. Therefore, if both scattering angles are below 4.8 mrad there is an ambiguity between muon and electron that must be resolved by µ/e discrimination.
The energy and the scattering angle of the electron in the LAB can be obtained by solving the boost relation E * 4 = γE 4 − βγ p 4 cos θ e for E 4 . This yields Going beyond the elastic process (2), by allowing for additional emission in the final state as described by (1)  have to extend the definitions of the momentum transfer. The variables now have to be distinguished. Sometimes it is useful to express t e in terms of the electron scattering angle as which follows directly from (10) and (11).
Since the contribution of the HVP to µ-e scattering is of central importance, in Figure 2 we show its leading effect as a function of the electron scattering angle θ e . More precisely, we show the NLO K factor defined as where σ (0) is the Born cross section and σ NLO h the hadronic contributions at NLO. As can be seen from Figure 2 and will be discussed in more detail in Section 4, the contribution of the HVP to µ-e scattering is larger for small θ e , whereas for θ e 20 mrad (corresponding to x 0.4) it is strongly suppressed.
The determination of the HVP will be obtained by a template fit of the shape of the distribution. For a simplified discussion it is useful to think in terms of a split into a signal region (small θ e ) and a normalisation region (large θ e ). In the signal region the effect is of the order of 10 −3 whereas in the normalisation region the HVP contribution amounts to 10 −5 and its error is expected to be below the experimental systematic uncertainty.

Fixed-order calculations
In order to achieve our goal of a relative accuracy of 10 ppm, we need to calculate µ-e scattering at least up to NNLO in the perturbative expansion in the electromagnetic coupling α ∼ 1/137. In addition, we need a flexible setup that allows for the computation of arbitrary infrared safe observables, i.e. a parton level Monte Carlo (MC). The latter aspects will be discussed in Section 6. In this section we discuss the main features of the analytic fixed-order computations.

Leading order
Starting at LO in QED there is a single diagram with a tchannel exchange of a photon. It is precisely this feature that makes this process ideal to extract the HVP. The dominant contribution of the HVP is simply given by the insertion of the hadronic bubble Π had into the photon propagator, as shown in Figure 3b. It is precisely the effect of this contribution that is shown in Figure 2.
As indicated in Figure 3 we often (formally) distinguish the charges of the electron and muon. Denoting them by q and Q respectively, the LO amplitude can be written as where the superscript indicates the number of loops. The two integer subscripts of the last expression indicate the power of q and Q. The two-particle final state is indicated by the subscript n of the first expression, where n = 2 is implicitly understood. To obtain a (differential) LO cross section dσ (0) we simply integrate the (squared) matrix element M (0) n over the two-particle phase space where cuts applied by the experiment and the definition of the observable are understood. The leading-order differential cross section is given by Because M (0) n ∼ 1/t 2 and, hence, dσ (0) /(dt) ∼ 1/t 2 the total cross section is not well-defined. Therefore, we always have to apply cuts to the integration to avoid the region t ∼ 0.
At LO, effects due to the electron mass m are suppressed by Hence, they have to be taken into account at -and even beyond -LO to achieve a 10 ppm prediction. The contributions due to the exchange of a Z boson are strongly suppressed because of its large mass M Z . However, the interference between the Z-boson and photon-exchange diagrams is suppressed with respect to the LO QED contribution only by Q 2 /M 2 Z ≃ 10 −5 . Hence, this effect is relevant and needs to be taken into account in the calculation.

Next-to-leading order
Going to NLO, the separately divergent real and virtual contributions have to be combined to obtain a physical result. Following earlier efforts [9,10], recently a fully differential NLO code [11] has been used to perform a detailed phenomenological study, taking into account the full m dependence. The NLO contributions can be split into gauge invariant parts by separating the contributions into powers of q and Q. Thus we decompose the NLO amplitude A (1) n ≡ A (1) (µe → µe) as The leptonic vacuum polarisation contributions are part of the full NLO calculation. However, we sometimes treat them separately as they are somewhat closer connected to the signal extraction. The virtual corrections are then obtained by integrating over the n = 2 parton phase space the renormalised squared matrix element Typically, the masses and wave functions are renormalised in the on-shell scheme, whereas for the coupling, either the onshell scheme or the MS-scheme can be used. Similarly, the real corrections are obtained by integrating over the n +1 = 3 parton phase space the squared matrix element  where the amplitude A (0) n+1 = A (0) (µe → µeγ) is also decomposed according to The cross section is obtained as the sum As illustrated in Figure 4, the terms ∼ q 4 Q 2 (∼ q 2 Q 4 ) in M (1) n and M (0) n+1 correspond to the corrections due to photon emission from the electron (muon) and are hence called electronic (muonic) contribution. There are also mixed terms ∼ q 3 Q 3 . The latter flip the sign if µ − is changed to µ + . As will be discussed in Section 6, electronic effects are actually dominant [11]. Regarding the virtual corrections, the electronic contribution only requires A (1) 3,1 which is simpler to compute than A (1) 2,2 . Such considerations become more important when discussing NNLO contributions.
Keeping a finite m complicates the computation of the virtual corrections. On the other hand, it serves as a regulator for collinear singularities which are replaced by log(m 2 /Q 2 ) and only soft singularities are left. In [11] the latter are regularised using a photon mass. There are two additional independent parton level Monte Carlo codes [12,13] using dimensional regularisation for IR singularities. These codes have been compared to [11] and full agreement has been found.
Electroweak (EW) NLO corrections are not expected to be required at the 10 ppm level. This was explicitly verified in [11].
3.3 Next-to-next-to-leading order A complete result for NNLO QED corrections to µ-e scattering is not yet available. However, there are already several partial results and a large theoretical effort is under way to complete the full NNLO calculation.
Following the notation of (15), (19) and (22) and using A n+2 = A(µe → µeγγ), the required amplitudes for the NNLO corrections are Similarly, for the matrix elements we need for the double-virtual (vv), real-virtual (rv) and double-real (rr) corrections. They have to be integrated over the n = 2, n + 1 = 3 and n + 2 = 4 parton phase space, respectively, The interplay between these three parts is illustrated in Figure 5 where different cuts to the same diagram squared represent contributions to dσ (rr) , dσ (rv) and dσ (vv) , respectively. From a theory point of view, there is a choice whether to include the sub-process µe → µe + ee in M (0) n+2 . Assuming m > 0, this is a separate IR finite contribution.
The main bottleneck for a NNLO calculation keeping the full m dependence is the evaluation of the two-loop amplitude A (2) n (m). It is not clear if a complete NNLO calculation with full m dependence is feasible in the next years. Fortunately, this is also not really required. The electronic contributions can be computed with full m dependence. For the remaining contributions, an approximate treatment for the NNLO corrections, i.e. an expansion in z is expected to be sufficient to obtain 10 ppm precision in the theoretical prediction.
Usually we refrain from listing the dependencies of the amplitudes on the momenta and the masses, m and M. However, sometimes we will have to indicate how we treat the dependence on the electron mass. Either we keep it completely as in A (2) n (m), or we set it to zero as in A (2) n (0). A third option is to consider an expansion in m, using m 2 ≪ {M 2 , Q 2 }. We will indicate this in the notation by writing A (2) n (z) having in mind m = zM with z ≪ 1.
Whether or not we keep the electron mass will alter the form of how the IR singularities of M (2) n manifest themselves. Using dimensional regularisation with d = 4−2ǫ, the highest pole of M (2) n (0) is 1/ǫ 4 which corresponds to double-softcollinear poles. On the other hand, M (2) n (m) and M (2) n (z) will only have 1/ǫ 2 poles. The double-soft-collinear poles are now replaced by 1/ǫ 2 log 2 (z).
A similar change happens in the real-virtual and doublereal contribution. Introducing an electron mass regularises the collinear singularities in the phase-space integration, again transforming the corresponding 1/ǫ poles to log(z) terms. Of course, for a physical cross section, all final-state collinear (and soft) singularities cancel. Thus, as for all regularisation procedures, a cross section is independent of which regularisation is chosen for the collinear singularities. The difference between dσ (2) (m) and dσ (2) (0) or dσ (2) (z) is in terms of the form z p log l (z) that are finite (and actually vanish) for m → 0. An advantage of the regularisation with a finite electron mass is that the initial-state collinear logarithms are manifest in the fixed-order contributions.
In what follows we will now consider how to obtain a sufficiently precise approximation to a complete NNLO calculation, comparing different approaches on how to treat the electron mass.

Massive electron
We start by noting that also at NNLO the electronic emission is the dominant contribution. This corresponds to the terms ∼ q 6 Q 2 in dσ (2) , with an example shown in Figure 5a. This part can actually be computed with full m dependence. It is a problem with only one active mass scale, m, in the loops and the two-loop virtual corrections 2 Re[A (2) 5,1 ×(A (0) n ) * ] can be obtained from the heavy-quark (actually lepton) two-loop form factor [14,15]. For the moment we do not include the HVP insertion in the two-loop diagrams. This will be dealt with in Section 4.
In order to combine this with the double-real and realvirtual contributions using dimensional regularisation, a suitable NNLO subtraction scheme has to be implemented. One example of such a scheme is the FKS 2 [16] which extends the NLO FKS subtraction scheme [17,18] to NNLO in the case of massive QED where only soft singularities are present. FKS 2 was successfully tested for the muon decay with full electron mass. Preliminary results for the q 6 Q 2 terms have (rr)(rv) (vv) (a) Example for a contribution to the squared matrix element ∝ q 6 Q 2 (rr)(rv) (vv) (b) Example for a contribution to the squared matrix element ∝ q 5 Q 3 recently been presented in [19,20]. Of course, it is trivial to adapt these computations for the purely muonic emission, i.e. the terms ∼ q 2 Q 6 . But they are expected to be numerically much less important. A similar calculation has been done in the context of lepton-proton scattering in [21], where the electronic terms for e p → e p scattering have been computed using a phase-space slicing method.
Contributions where both emission from the electron and muon line are involved are technically much more challenging. An example is shown on Figure 5b. In fact, an analytic computation of the two-loop amplitude, keeping the full m dependence, is probably not feasible in the near future. However, the reduction to master integrals of M (2) n (m) is currently under investigation. This could be combined with a numerical evaluation of the master integrals. Another approach would be a completely numerical evaluation of the amplitude, even avoiding a reduction to master integrals. Even if a complete result for M (2) n (m), suitable for a Monte Carlo code, is not expected to be available within the next few years, these efforts are extremely useful as cross checks for other approaches (see below).
Apart from the two-loop amplitude also the full realvirtual corrections M (1) n+1 (m) need to be computed, i.e. the interference between one-loop and the Born amplitude of eµ → eµγ. Even though their integration over the n + 1 = 3 phase space entails an IR (soft) singularity, the order ǫ terms Fig. 6 Examples of Feynman diagrams ∼ q 3 Q 2 contributing to the real-virtual corrections to µ-e scattering at NNLO in QED. Related diagrams ∼ q 2 Q 3 , where the real photon is radiated from a muonic line, can be obtained from these by means of symmetries.
of M (1) n+1 (m) are not really required, if a suitable subtraction scheme such as FKS 2 is used.
The calculation of these interference terms was considered in [22], where both cases with massless and massive electrons were studied. The real-virtual contributions are encompassed in the A (1) n+1 term, and comprise 44 Feynman diagrams which are generated with the M packages F A and F C [23,24]. Four of these diagrams do not actually contribute since they automatically cancel out at integrand level because of Furry's theorem, while the remaining ones can be split into two sets: i) where the real photon emission occurs from a muonic internal or external line and ii) where the emission is from an electron line. When both leptons are massive the contributions from the two sets can be related via symmetries, namely exchanging the electron and muon masses and charges as well as the respective external momenta. This fact was exploited to halve the number of diagrams to be evaluated to just 20. Representative diagrams are depicted in Figure 6.
These unrenormalised amplitudes are then inserted in the real-virtual interference term (29), which is computed with massive electrons to assess the impact of the massless electron approximation enforced for the calculation of other contributions. The first steps towards a fully-analytic evaluation of this contribution were undertaken, using the same automatic framework employed for the calculation of M (2) n . The integrands depend on the two mass parameters plus five kinematic variables, which were parametrised using the Momentum Twistor parametrisation [25][26][27] in preparation for the adaptive integrand decomposition in A . Subsequently the amplitudes were simplified via integration-by-parts identities [28][29][30] generated with the package K [31], identifying 45 master integrals.
The interferences were then matched with counterterm amplitudes generated in F C , employing the on-shell renormalisation scheme. The cancellation of the leading ultraviolet (UV) poles in dimensional regularisation in the renormalised amplitudes was verified numerically.

Massless electron
Neglecting the electron mass reduces the difficulty of the problem from extreme to very high. Fortunately, there has been an impressive effort devoted to this computation, such that the evaluation of the two-loop amplitude for massless electrons, A (2) n (0), is close to completion. The amplitude A (2) n (0) receives contributions from 69 Feynman diagrams, which are generated with the help of the packages F A /F C [23,24] and its evaluation requires the calculation of O(10 4 ) integrals. Owing to the use of adaptive integrand decomposition [32,33], implemented in the in-house package A [34], and integration-by-parts identities [28][29][30], implemented in the public routines R - [35,36] and K [31], the amplitude -to be precise, the interference term 2 Re A (2) n ×(A (0) n ) * -has been simplified, and written as a linear combination of an integral basis formed by about 120 elements [37,38], dubbed master integrals. The latter have been successfully evaluated by means of the differential equation technique [39][40][41][42][43] in combination with the Magnus exponential method [44,45]. Originally evaluated in a non-physical region, where the mathematical complexity was found to be more limited, the analytic evaluation of the master integrals to the physical scattering region was recently obtained [46]. The analytic expressions of the master integrals, numerically evaluated with the help of G N C [47], were successfully validated against the numerical values provided either by S D [48] or, for the most complicated integrals, coming from the 7-propagator graphs, by an in-house algorithm. Representative diagrams are depicted in Figure 7.
The non-trivial evaluation of the (unrenormalised) twoloop amplitude for the µ-e scattering required the development of a high-level automated tool, exploiting the synergy of different packages embedded in a M framework, whose flowchart is depicted in Figure 8. UV divergences arising from divergent loop integrals are regularised within the dimensional regularisation scheme, to be later removed by means of a diagrammatic approach to the renormalisation. In particular, the counterterm Lagrangian provides additional Feynman rules, which can be adopted in our automatic frame- work for generating the additional diagrams, yielding a UV finite amplitude. Given the masslessness of the electron, we adopted an hybrid renormalisation scheme choice: -MS scheme for the coupling; on-shell scheme for the muon mass.
For a recent review on the state of the calculation, see also [49]. In principle, the full computation can be done with massless electrons. In the phase-space integration, this results in collinear singularities. Hence the subtraction (or any other) scheme used will need to be adapted to this case. Unfortunately this will destroy the simple divergence structure of massive QED that was exploited in FKS 2 .
There is one further subtlety when performing the calculation with massless electrons everywhere. While final-state singularities will cancel for m = 0 in any sufficiently inclusive observable [50], the same cannot be said about initialstate collinear singularities1. The corresponding ǫ poles will remain unless properly treated. There are multiple somewhat related ways to make these expressions well defined such as the Weizsacker-Williams approach, the structure function approach [51,52] or the QED parton distribution function approach [53]. These techniques where honed in the LEP era and will work at the required accuracy.
A final problem with a purely massless calculation is the restrictions imposed by the phrase 'sufficiently inclusive 1Such poles exist even though the initial-state electron is at rest as the total cross section is Lorentz invariant. observable' of the KLN theorem [50]. This will make a quantity such as θ e inaccessible without breaking IR safety or defining a jet-like observable.

Massified electron
Given that the electron mass is a natural cutoff for collinear emission, it seems to be natural to use m as a collinear regulator. Apart from reducing the complexity of the IR subtraction for the real integration, this will also facilitate the combination of a fixed-order result with parton-shower Monte Carlo codes and automatically produce the log(m) terms that are present in distributions.
In order to do this, we will have to massify M (2) n (0), i.e. transform it into M (2) n (z) that contains the leading logarithmic terms log(m). Initially, this problem has been addressed for Bhabha scattering [54,55] and then been generalised [56,57] using a factorisation approach. A further generalisation is needed if two different non-vanishing masses exist, as in our case M and m. This has been studied in the context of the muon decay [58] and will allow to obtain M (2) n (z) from M (2) n (0). To achieve this, an approach based on soft collinear effective theory (for an introduction to SCET see [59]) and the method of regions [60] is used. Loop integrals contributing to M (2) n (m) are expanded in z by taking into account all relevant scalings of the loop momenta k i (regions) and expand the integrand in all these regions. After expansion of the integrand, the integrations are simplified and adding up all contributions reproduce the expansion of the full integral. In our case, the relevant regions are hard k i ∼ (1, 1, 1), soft k i ∼ (z, z, z) and, ultrasoft k i ∼ (z 2 , z 2 , z 2 ). Further, we need collinear k i ∼ (z 2 , 1, z) for the in-coming and anti-collinear scaling k i ∼ (1, z 2 , z) for the out-going electron. Here we have used light-cone coordinates k i = (k + , k − , k ⊥ ). Each external electron defines a collinear direction (either the one of k − or of k + ) that has to be taken into account.
In principle, this can be done to any power in z. However, restricting ourselves to the leading power, the matrix element factorises as The hard contributions correspond to the massless matrix element M (2) n (0) the computation of which was discussed above. The soft part, S, is also process dependent and will have to be computed for µ-e scattering. However, the computation of the soft part is much simpler than the full amplitude. It obtains contributions from fermion-loop diagrams and can be tested against a fully massive computation of the fermion loops [61]. The collinear contributions are contained in the factor Z i (m) that is process independent and known [58] and has to be added for each external electron. Finally, ultrasoft contributions exist for individual integrals and diagrams, but they cancel for the amplitude in agreement with the SCET expectations.
This result can now be combined with a fully massive evaluation of dσ (rv) (m) and dσ (rr) (m). However, the fully massive real corrections contain poles that will not naively match the poles obtained through massification, instead causing a mismatch at O(z). This mismatch can be avoided by either expanding the analytic poles of the real corrections or calculating the fully massive poles of the two-loop amplitude from first principles [62,63]. For the phase-space integration only soft singularities have to be regularised with dimensional regularisation. Putting everything together results in a Monte Carlo code that provides results complete at NLO and includes all leading in z terms at NNLO. However, it does not systematically include non-leading terms at NNLO, i.e. terms that vanish in the limit z → 0, such as α 2 z log(z). It should also be mentioned that in the region t → −0 the counting used in this expansion breaks down. Of course the cross section is divergent in this region anyhow such that this problem can be avoided with an appropriate cut.
3.4 Beyond next-to-next-to-leading order A complete calculation at N 3 LO is a daunting task. However, keeping in mind that the dominant contribution to any loop order stem from emission of the electron, i.e. terms O(q 2+2n Q 2 ), at least a partial N 3 LO result might be achievable. Once more it is the remarkable simplicity of QED that allows us to extend the subtraction scheme proposed for NNLO, FKS 2 , to even higher loop orders [16]. The N 3 LO extension, aptly named FKS 3 , has already been worked out and shown to retain the simplicity of FKS.
The necessary ingredients for this endeavour are the A (3) 7,1 part of the three-loop A (3) n , the A (2) 6,1 part of the two-loop A (2) n+1 , the A (1) 5,1 part of the one-loop A (1) n+2 , and the A (1) 4,1 part of the tree-level A (0) n+3 . The latter two are, at least in principle, easy to obtain thanks to the advances made in the automation of one-loop calculations. The former two are more challenging. However, impressive progress has been made in calculating the heavy-quark form factors at threeloop [64,65] and an efficient tool to numerically evaluate generalised polylogarithms [66] is available. A big remaining problem is the two-loop real-double-virtual contribution. However, at least in principle, it should be possible to adapt and massify the calculations performed for γ * → qqg which are part of the NNLO calculations to three-jet production.

Next-to-leading order
The NLO and NNLO corrections to the muon-electron differential cross section involve non-perturbative QCD contributions given by diagrams with an HVP insertion in the photon propagator (see Figs. 3b and 9). Let us define the SM vacuum polarisation tensor with four-momentum q as where j µ em (x) = f Q fψ (x)γ µ ψ(x) is the electromagnetic current and the sum runs over fermions with charges Q f . The function Π(q 2 ) is the renormalised vacuum polarisation satisfying Π(0) = 0. It is commonly subdivided into the leptonic part Π lep , which receives contributions only from charged leptons, the hadronic part Π h , from hadrons containing the five light quarks u, d, s, c, b, and the contribution from the top quark Π top . The weak interaction will be ignored.
In perturbation theory Π lep and Π top can be computed order by order in α and the strong coupling α s [67][68][69][70]. On the contrary, the HVP cannot be calculated in perturbation theory for any value of q 2 because of the non-perturbative nature of strong interactions. Nevertheless, we can express Π h in terms of the measured cross section of the reaction e + e − → hadrons [71] thanks to the subtracted dispersion relation and the optical theorem where and is the effective fine-structure constant. The numerical value for the HVP can be obtained by using the Fortran libraries alphaQEDc17 [72][73][74][75], KNT18VP [75][76][77][78][79][80], and VPLITE [80,81] based on hadronic e + e − annihilation (timelike) data. The hadronic contribution to the µ-e cross section at NLO, due to the diagram in Figure 3b, is where ∆α h (t) = −Π h (t) is the leading hadronic contribution to the running of α(t). The goal of the MUonE experiment is the extraction of ∆α h (t) from µ-e scattering data. Note that the NLO hadronic corrections incorporate also the contribution from the diagram in Fig. 9a where a virtual photon is emitted and reabsorbed by the hadronic insertion. This irreducible part of the second-order hadronic contribution to the running of α(t) is not considered as part of the NNLO corrections because its effect is commonly included in the ratio R(s) as final-state radiation [82,83]. The impact of the hadronic contribution at NLO is shown in Figure 2 as a function of θ e . For later reference, in Figure 10 we show the same contribution as a function of t e . The factor K NLO h (t e ) depicted in Figure 10 is defined in analogy to (14). In accordance with Figure 2, the effect is larger for large values of |t e |.
Before we move on to the hadronic corrections at NNLO, following the analysis of [7] we will briefly discuss the impact of the SM corrections -and possibly NP -on the extraction of ∆α h (t) at MUonE. This experiment will extract ∆α h (t) from the shape of the differential µ-e scattering cross section by a template fit method [2]. The basic idea is that ∆α h (t) can be obtained measuring, bin by bin, the ratio N i /N n , where N i is the number of scattering events in a specific t-bin, labelled by the index i, and N n is the number of events in the normalization t-bin corresponding to x(t) ∼ 0.3 (for this value of x, ∆α h (t) is comparable to the experimental sensitivity expected at MUonE and its error is negligible). Therefore, this measurement will not rely on the absolute knowledge of the luminosity. To extract the leading hadronic corrections to the µ-e scattering cross section in the t-bin i, one can split the theoretical prediction into where σ is the LO QED prediction integrated in the t-bin i, 2∆α h,i is the leading hadronic correction obtained from (36), δ i is the remainder of the SM corrections, and δ NP,i is a possible NP contribution. The experimentally measured ratio N i /N n can then be equated with the ratio of the theoretical predictions, As ∆α had,n is known with negligible error, if (δ i − δ n ) is computed with sufficient precision, one can extract 2∆α h,i + δ NP,i − δ NP,n , bin by bin, from N i /N n . Equation (38) shows that the impact of the SM corrections on this extraction can only be established after subtracting their value in the normalization region, and that the MUonE experiment will not be sensitive to a NP signal constant in t relative to the LO QED one, i.e. such that δ NP,i = δ NP,n [7].

Next-to-next-to leading order
At NNLO, we split the hadron-induced corrections to µ-e scattering, of order α 4 , into four classes of diagrams. The first three classes contain factorisable contributions, i.e. amplitudes that can be written as the product of a QED amplitude times the function Π h (q 2 ) evaluated at some q 2 fixed by the external kinematics. They are: Class I: tree-level diagrams in combination with one or two vacuum-polarisation insertions (Figure 9b). Their contribution to the differential cross section is proportional to Π h (t)[Π h (t) + 2Π lep (t)], the reducible part of the second-order hadronic contribution to the running of α(t).
Class II: QED one-loop diagrams in combination with one HVP insertion in the t-channel photon (Figure 9c). Their contribution to the differential cross section is proportional to Π h (t) and a combination of one-loop QED corrections to µ-e scattering. Class III: real photon emission diagrams with a vacuumpolarisation insertion in the t-channel photon (Figure 9d). They contain terms proportional either to Π h (t e ) or to Π h (t µ ).
Moreover a fourth class of non-factorisable diagrams must be considered: Class IV: one-loop QED amplitudes with a hadronic vacuum polarisation insertion in the loop. They can be further subdivided into vertex and box corrections (Figure 9e).
There are no light-by-light contributions to the µe cross section at NNLO (order α 4 ) -they appear at N 3 LO (order α 5 ). In addition, the analysis of future µ-e scattering data will also require the study of µ-e scattering processes with final states containing hadrons, as for instance e − µ ± → e − µ ± π + π − and e − µ ± → e − µ ± π 0 . However, as √ s ≃ 400 MeV, the available phase space is quite small: √ s − M − m − 2m π 0 ≃ 20 MeV for the former and √ s − M − m − m π 0 ≃ 160 MeV for the latter process.
As the HVP per se is of non-perturbative nature, the hadronic NNLO corrections rely inevitably on some external data for their numerical evaluation. These inputs can be of two kinds: we can either use the R ratio and the traditional dispersive method, or we can dismiss the e + e − → hadron dataafter all, MUonE aims at measuring a HLO µ independently on R -and employ the very same space-like data measured by MUonE. The two approaches are the following.
To R: The traditional approach to calculate the amplitudes in class IV uses the dispersion relation to replace the dressed photon propagator inside the loopwhere q now stands for the loop momentumwith the r.h.s. of (33), where the momentum q appears only in 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 R ratio. Also for the amplitudes in classes I-III we rely on the dispersion relation (33) to compute the HVP in the space-like region. This method was employed, for example, to compute the hadronic corrections to muon decay [84,85] and Bhabha scattering [86][87][88]. The hadronic NNLO corrections to µ-e scattering based on the R ratio were presented in [61]. Two independent codes were developed. The first is a standard Monte Carlo which uses Collier [89] for the evaluation of the oneloop tensor integrals and employs the FKS subtraction scheme [17,18]. The second code is developed in M and takes advantage of the analytic expressions of the one-loop integrals from Package-X [90] and M 's arbitrary-precision numbers to check for numerical instabilities during the dispersive and phasespace integrations. Perfect agreement was found between the two implementations. The results are presented below. Not to R: This alternative approach was presented in [91].
The factorisable diagrams in classes I, II and III depend on Π h (t), which is the quantity extracted by MUonE from the diagram in Figure 3b. As discussed in [91], also the non-factorisable corrections in class IV, where the vacuum polarisation appears inside a loop, can be calculated employing the HVP in the space-like region, without making use of the R ratio. Indeed, the loop integrals containing Π h can be computed via the hyperspherical integration method. After introducing spherical coordinates for the loop momentum and continuing internal and external momenta to the Euclidean region, one can write the loop propagators as an expansion in Gegenbauer polynomials. Then, the integration over the angular variables is performed analytically thanks to the orthogonality properties of these polynomials, so that each diagram is eventually cast in the form of a residual radial integration, which is computed numerically once provided with the HVP in the space-like region. The expressions of the kernels f (Q 2 , s, t) were presented in [91]. Their implementation into a numerically stable code is necessary for future use in the Monte Carlo. The fact that the hadronic NNLO corrections can be obtained from Π h (q 2 ), with just q 2 < 0, suggests the possibility for MUonE to determine the HVP in an iterative way without making use of the R ratio. As a first step the hadronic NNLO corrections can be switched off in the Monte Carlo and a first approximation for Π h (q 2 ) extracted. Afterwards, the Monte Carlo can be supplied with such first approximation to compute the hadronic NNLO corrections, then a second approximation extracted and the process further iterated. Alternatively, if a functional form for Π h (q 2 ) is chosen to fit the HVP [92] the same ansatz can be employed at NNLO, under the assumption that it satisfies the correct asymptotic behaviour at infinity.
The dispersive and the hyperspherical methods are of course identical from the mathematical point of view; however the to R and the not to R approaches differ for the underlying theoretical assumptions. If we use the R ratio, we make a distinction between the HVP entering at NLO and at NNLO. On the one hand, at NLO we leave the HVP in a free form to be fitted from data. On the other hand, at NNLO we choose a different Π h (q 2 ) whose values are given by the R ratio via the dispersion relation.
On the contrary, by employing the hyperspherical method in the not to R approach we treat the HVP in a consistent way to all orders without making any a priori assumptions. Moreover, only in latter case, the MUonE determination of a HLO µ becomes truly independent and completely uncorrelated from time-like measurements.
Let us now discuss the size of these hadronic corrections. The ratio of the NNLO hadronic contribution to the µe differential cross section, with respect to the squared momentum transfer t e , and the LO prediction, is shown in Figure 11 for the processes µ + e − → µ + e − (left panel) and µ − e − → µ − e − (right panel) for E 1 = 150 GeV. These corrections were computed in [61] using the dispersive approach. The black 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. Figure 11 shows that for a muon beam with energy 150 GeV, most of the kinematic region scanned by the momentum transfer t e results in a factor K NNLO h (t e ) which is of order 10 −4 -10 −5 . These corrections are therefore larger than the O(10 −5 ) precision expected at the MUonE experiment.
The uncertainty on K NNLO h due to the error on R was estimated 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. Therefore a relative uncertainty of 1% was assigned on dσ NNLO h /dt e , 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.

Beyond fixed order
In Section 3 we have discussed the computation of the cross section in a strict expansion in the coupling α. At N n LO the cross section contains large logarithms of the form that potentially invalidate the perturbative expansion since αL m is not necessarily a good expansion parameter. The prefactor α 2 is from the Born cross section. In the total cross section the only logarithms that survive are due to initialstate collinear emission from the electron, as discussed in Section 3. However, for differential cross sections final-state collinear logarithms L m can be present. Going beyond the total cross section, as is required for MUonE, can result in additional large logarithms. In order to select elastic scattering, the emission of real radiation has to be restricted. Naively, this is done by vetoing photon emission with energy larger than a cutoff ∆. For the moment we ignore the fact that in practice another observable has to be chosen since photons are not detected. Restricting the emission of real radiation will result in additional logarithms of the form L ∆ = log(∆ 2 /Q 2 ) which again can be large if the cut is severe, i.e. ∆ is small. Thus for each order in α we obtain up to two powers of large logarithms. Of course, this is closely related to the 1/ǫ 2 (or the 1/ǫ log(m)) singularities for each perturbative order, discussed in Section 3. Therefore, the nth order correction to a differential cross section has the structure where the sum runs over 0 ≤ n 1 , n 2 ≤ n. Again, the prefactor α 2 is due to the Born term. In what follows we will omit this factor when discussing powers of couplings and logarithms.
Another potential source of large logarithms is related to the so-called factorisation (or collinear) anomaly [93][94][95]. This is related to the breaking of a scaling symmetry between collinear and soft modes in SCET and occurs due to the presence of two non-vanishing masses [58]. In practice it means that the separate factors of (31) might contain singularities that are not regularised through the usual dimensional regularisation. While these singularities cancel between the various factors on the r.h.s. of (31) the left-over of these cancellations corresponds to a logarithm of the form log(mM/Q 2 ).

Leading logarithm
The terms of (42) with n 1 = n 2 = n are the leading logarithms (LL). They can be resummed using a parton shower (PS). A PS has the advantage that the kinematics of the emitted photons is retained so that exclusive events can be generated. This makes sure that the resulting program remains fully-differential in all resolved particles which cannot be guaranteed in analytic calculations.
Roughly speaking, there are two avenues to numerically resum the LL contributions. The starting point is either soft emission or collinear emission. In the first case, the wellknown Yennie Frautschi Suura (YFS) exponentiation [62] of soft emission is used. This allows for a numerical implementation taking into account soft emission to all orders [96][97][98][99]. Such a resummation is well suited to be combined with fixed-order calculations performed with FKS 2 , the subtraction scheme suggested earlier. In fact, FKS 2 exploits the YFS structure of the matrix elements where M f (ℓ) n is free of IR singularities. The latter are all absorbed by the exponential of the integrated eikonalÊ that governs soft emission. For a precise definition of all quantities in (43) and more details see [16]. A recent example where a NNLO QED calculation is merged with a YFS resummation can be found e.g. in [100].
Taking collinear emission as a starting point, a QED parton shower can be constructed through subsequent collinear emission of photons governed by the e → e γ splitting kernel P(z) = (1 + z 2 )/(1 − z), where z is the momentum fraction of the electron after the split. This procedure has been used by the BabaYaga [101][102][103][104][105] event generator. It can be combined with fixed-order calculations and extended to next-to-leading collinear logarithms.
As both, YFS Monte Carlo and QED parton shower include the leading soft-collinear emissions, they agree at LL. Going beyond LL, the QED PS also includes hard (non soft) collinear radiation, i.e. it includes all leading collinear logs α n L n m . It can be further adapted to also include soft wideangle emission [106]. As we will discuss below and in Section 6, this difference beyond various implementations at LL can reveal useful information to assess the theoretical error. To exploit this, work is ongoing to implement µ-e scattering in both frameworks and compare.

Next-to-leading logarithm
The terms of (42) with n 1 + n 2 = 2n − 1 are the next-toleading logarithms (NLL). Whether or not their resummation is required and possible depends on the precise definition of the quantity that selects elastic scattering and on the value of the cut parameter ∆. A partial resummation of NLL terms can be done with improved Monte Carlo generators.
A complete resummation beyond LL requires a precise definition of the physical quantity for which the resummation is carried out. Given that for µ-e scattering we are not primarily interested in particular distributions of certain observables, but rather in a precise description of the fiducial cross section measured by MUonE, it is not possible to precisely match the observable to the measured quantity.
The most important cut that has to be made for the extraction of the HVP is to choose elastic events. In theory this can be achieved in several ways, all of which restrict the phase space for emission of photons. As a first example we mention a cut on the invariant mass m eγ of the electron-jet, i.e. the cluster of the outgoing electron plus all potentially emitted photons. This quantity is sensitive to large-angle soft emission and, contrary to e.g. m 2 µγ , also to small-angle hard emission. We thus define where according to (1) p X = i=1,n p iγ is the sum over (up to n) photon momenta. The elastic events can be chosen by making a cut m eγ − m ≤ ∆. Another well-studied option is to use the transverse momentum. In our case we have to take the transverse momentum p 4⊥ of the electron w.r.t. its tree-level direction, which can be determined from the muon scattering angle θ µ . In case of no photon emission, p 4⊥ → 0 and we can impose a cut |p 4⊥ | 2 < ∆ 2 to select elastic events. While these quantities are useful from a theoretical point of view and likely enable a resummation beyond LL, they are unfortunately not very useful from an experimental point of view. None of these quantities can actually be measured since neither are photons detected nor are the momenta or energies of the electron and muon measured. In practice, the experimental procedure to select elastic events has to rely solely on the scattering angles. Such a quantity is likely to be rather involved and not amenable to direct resummation. In Section 6 we will consider an acoplanarity cut as one example of a cut that can realistically be applied by MUonE to restrict radiative events.
A possible way forward is to use the analytic resummation of several different variables to construct approximations for perturbative coefficients beyond those included in the fixed-order approach. These results can then be implemented as approximate matrix elements in a fully differential parton-level Monte Carlo. Through comparisons of results obtained by using different versions of resummation it is possible to obtain a realistic error estimate of the approximations. This procedure has been used for example for top-quark pair production [107] and provided an improved prediction with a robust estimate of missing terms. A quantity that offers itself for resummation in the context of µ-e scattering is dσ/dt e where in the region t e → t min large logarithmic corrections are present.

Estimate of the theory error
As discussed earlier in the present section, higher-order contributions can be included according to different methods, such as a QED PS algorithm, YFS MC exponentiation or analytic resummation. Regardless of the approach used to account for the corrections beyond NNLO, the accuracy of the theoretical predictions due to missing perturbative contributions must be carefully estimated, as it represents a component of the total systematic error.
For this purpose, it seems advisable to evaluate the theoretical uncertainty step-by-step, as the different theoretical ingredients become available. In the following, a possible strategy for the theoretical uncertainty estimate is illustrated, at the level of differential distributions. It is assumed that, in addition to the fixed-order NLO QED calculation, also the NNLO QED matrix elements are implemented in a fully fledged Monte Carlo simulation tool. The technical accuracy, related to the details of the implementation of the fixed-order radiative corrections, can be controlled by means of two completely independent codes, which are assumed to exist.
In the relatively short term, a first assessment of the theoretical accuracy could be given as follows: by comparing the predictions for the photonic corrections at NLO and NNLO accuracy. This comparison can be performed for the full set of corrections but also separately for the gauge-invariant subsets of contributions due to electron radiation, muon radiation and electronmuon interference. The importance of this procedure is twofold as a) it would allow to settle the hierarchy of the different classes, which is a crucial prerequisite to identify the sources of corrections that need to be resummed at all orders and b) it would provide information about the convergence of the perturbative series (in particular, for those kinematical regions where the NLO corrections are particularly large). An estimate of the missing third order can be given by by computing the NNLO leptonic and hadronic corrections due to the combination of the two-loop vacuum polarisation contribution and real pair emission. It is known that these corrections give rise to large collinear logarithms but also that they are typically smaller than purely photonic corrections. In particular, the contribution due to electron loop and real e + e − radiation produce collinear logarithms L m . Taken separately, the LO cross section of the process eµ → eµ(e + e − ) results in contributions α 2 L 3 m . However, if combined with the virtual electron-loop contributions the α 2 L 3 m cancel and we are left with collinear logarithms α 2 L 2 m and α 2 L m . The computation of this class of corrections therefore would allow to probe the size of those NNLO logarithmicallyenhanced corrections that are of non-photonic nature; by comparing the finite-order expansion of a given resummation approach with the exact perturbative calculation at NNLO accuracy. Again, this comparison could be performed for the complete set or the gauge-invariant subsets of photonic corrections and would allow to quantify the size of the NNLO remainder beyond the LL approximation at O(α 2 ).
Over the longer term, assuming that different methods to account for the contribution of multiple photon radiation will be available and matched to the NNLO calculation, the theoretical uncertainty could be more reliably estimated as follows through a comparison of the exact NNLO calculation and the O(α 3 ) expansion of a given resummation procedure. From this comparison, one would get an evaluation of the whole set of higher-order contributions beyond NNLO; by comparing the all-order predictions of the different methods developed for the description of multiple photon radiation. As remarked earlier, approaches such the QED PS and YFS MC exponentiation provide the same LL structure but may differ in the partial resummation of NLL contributions. Moreover, this comparison could be extended, wherever possible, to include the results of analytic resummations possibly featuring a complete resummation beyond the LL approximation. As a whole, this procedure would provide a robust estimate of missing higher-order terms at NLL accuracy; under the assumption that two independent implementations of a MC code based on the matching of NNLO corrections with resummation will be available and different techniques for exponentiation will be used, a comparison between the predictions of the two codes could provide further important information about the theoretical accuracy. Actually, because of the reasons already emphasised, the two calculations are expected to differ for contributions dominated by terms of the order of α 3 L 2 m . Hence, this comparison would allow to probe the size of the most important NLL O(α 3 ) contributions, similarly to the previous point above but at the level of completely matched formulations.
Besides the contributions due to purely photonic corrections, the extreme accuracy of MUonE will presumably also demand for the inclusion of the dominant effects beyond NNLO from fermionic corrections. Among those, contributions due to electron pairs are the most important. The evaluation of NLO photonic corrections to the cross section eµ → eµ(e + e − ) combined with the corresponding virtual electron-loop contributions will exhibit α 3 L 3 m terms. The determination of these terms could be achieved, for example, by convoluting the NNLO cross section with a standard QED PS simulation or by means of an appropriate generalization of the basic ingredients of the PS algorithm. The resummation of pair production contributions can be also shown to take place to all orders of the perturbative expansion [108][109][110].
If necessary, many of the above estimates could be put on firmer ground by computing the full set of virtual and real photon corrections due to the radiation from a single leg at N 3 LO accuracy, as discussed in Section 3.4.
To summarize, the accuracy of NNLO calculations combined with the contributions due to multiple photon radiation will be limited by the approximate inclusion of NLL contributions at O(α 3 ). A careful estimate of their impact on the observables measured by MUonE will set the scale of the overall theoretical uncertainty.

Monte Carlo
In the MUonE experiment, the extraction of the HVP contribution to the effective electromagnetic coupling will be based on a template fitting method. In this procedure, differential cross sections are calculated according to a given theoretical input and compared to the data, as a function of the parameters entering the ∆α had (q 2 ) modelling. Inevitably, this requires the implementation of the theoretical predictions into a fully flexible MC code. The latter is also needed for a high-precision calculation of the normalisation cross section, as well as the evaluation of the detector efficiencies and the assessment of a number of experimental systematics. Thus, the MC is the experimentally-oriented completion of any theory calculation as it goes to the heart of the data analysis.
We describe here what is presently available in the sector of MC tools for simulations of the µ-e scattering process and the most important phenomenological results. We are also interested in providing a recipe on how to convert future theoretical achievements or theory MC into useful tools for the experimentalists and phenomenologists. A sketch of the ongoing efforts towards the realisation of MC codes with increased accuracy or the simulation of relevant contributions to µ-e scattering is also given.
Until now, feasibility studies and preliminary simulations by the MUonE collaboration have been performed using a MC event generator that includes NLO electroweak corrections. The theoretical content of the NLO MC is described in detail in [11] and will not be repeated here. Suffice it to say that the MC developed in [11] is based on a calculation of the full set of NLO electroweak corrections to µ-e scattering without any approximation, including finite mass contributions. More interesting facts are the main computational features of the NLO MC. They can be summarised as follows: the generated events are fully exclusive, i.e. all the momenta of the event particles can be stored in such a way that any observable can be studied and any further effect can be applied (experimental cuts, detector simulation, etc); both weighted and unweighted (constant weight) events can be generated. The use of weighted events speeds up event generation and, generally, reduces the statistical error due to MC integration; the incoming muon energy (beam momentum) can be spread by a Gaussian distribution around its nominal value, to match realistic beam preparation; the HVP contribution can be switched on and off, all the rest of the input parameters remaining unchanged. This gives the possibility of studying the contribution of ∆α had (q 2 ) to any observable at NLO accuracy, including experimental effects2; the generated events can be stored into R n-tuples for further analysis. The storage format includes all the relevant information for each run input and for each generated event. The flexible nature of the adopted format makes it suitable to facilitate the implementation of future theoretical developments, such as the inclusion of multiple photon emission; the code is equipped with a R interface for reading, analysing and manipulating the generated samples.
In the following, we show a sample of particularly interesting predictions obtained by means of the above NLO MC. Within the set of numerical results described in [11], we select those that are particularly relevant in the light of the efforts in the sector of NNLO corrections and resummation.
To that purpose, we are interested to address the following questions3: how the θ e -θ µ correlation of the elastic signal is affected by QED radiation at NLO and how the signal sensitivity can be recovered by applying suitable cuts; how the full NLO QED correction is shared among the different gauge-invariant subsets described in Section 3; how large finite electron-mass contributions are.
To answer the above questions, we provide numerical results for both the µ − e − → µ − e − and µ + e − → µ + e − process, since both options are relevant for the MUonE experiment and the mixed QED corrections ∼ q 3 Q 3 differ in the two cases.
We use the following input parameters: (46) where α(0) is the value used for the lepton-photon coupling. For the energy of the incoming muons, we assume E 1 = 150 GeV, which is the energy of the M2 beam line of the CERN SPS. Note that, under the fixed-target configuration of the MUonE experiment, the CMS energy corresponding to this muon energy is given by √ s ≃ 0.405541 GeV and that the Lorentz γ factor boosting from CMS to LAB is γ ≃ 370. Due to (11), a lower limit on E 4 implies an upper limit on t e . In this kinematical condition, the collinear logarithms L e = ln(|t max |/m 2 ) and L µ = ln(|t max |/M 2 ) amount to L e ≃ 13.4 and L µ ≃ 2.7, respectively.
To study the dependence of the radiative corrections on the applied cuts, we consider two different event selections defined by the following criteria: 1. θ e , θ µ < 100 mrad and E 4 > 0.2 GeV (i.e. t e −2.04 · 10 −4 GeV 2 ). The angular cuts model the typical acceptance conditions of MUonE and the electron energy threshold is imposed to guarantee the presence of two charged tracks in the detector; 2. the same criteria as in Setup 1, with an additional acoplanarity cut, applied to partially remove radiative events and thus enhancing the fraction of elastic events. We require acoplanarity π − |φ e − φ µ | lower than 3.5 mrad, for the sake of illustration.
The answer to the first question about the impact of NLO QED radiation on the θ e -θ µ elastic correlation is given by Figure 12. In that figure, we compare the correlation in the laboratory frame between the scattering angles of the outgoing electron and muon at LO and NLO, for the Setup 1 and Setup 2 defined above. It can be noticed that, in the absence of an acoplanarity cut (Setup 1), the correlation present at LO (elastic curve) is largely modified by the presence of events at relatively small muon angles, which originate from 3We focus on photonic corrections, as the contribution of purely weak NLO corrections is well below the 10 ppm level, as shown in [11].  Fig. 12 The correlation between the electron scattering angle θ e and muon scattering angle θ µ for the µ + e − → µ + e − process at LO (elastic curve) and NLO QED, for the selection criteria 1 and 2 defined in the text.
the bremsstrahlung process µ + e − → µ + e − γ. However, the tight acoplanarity cut (Setup 2) turns out to be effective in getting rid of most of these radiative events, thus isolating the elastic correlation curve. As shown in [11], in the presence of acceptance cuts only (Setup 1), the corrections to the electron scattering angle turn out to be quite sizeable at small angles, due to the emission of a hard photon in the radiative process µe → µeγ. However, this effect gets largely reduced when an elasticity cut is applied (Setup 2), yielding a correction in the 10-40% range for all the relevant distributions. In the presence of an elasticity cut that vetoes hard photon emission, the contribution of soft photons becomes enhanced and gives rise to large IR logarithms, as remarked in Section 5. Not to invalidate the perturbative expansion, those logarithms need to be resummed together with the contributions due to collinear emission. This can be achieved by means of exclusive MC techniques, such as YFS exponentiation or QED Parton Shower, or analytic resummation. As emphasised in Section 5, the latter method requires the identification of a kinematical quantity able to select elastic scattering. This poses the question how the elasticity band isolated by the cuts of Setup 2 can be approximated by a reasonably simple 'observable' suitable for resummation. To understand how the gauge-invariant subsets contribute to the overall NLO QED correction, we show in Figure 13 the impact of the different classes described Section 3 on the dσ/dt e (top plot) and dσ/dt µ (bottom plot) distributions. The upper (lower) panels refer to the µ + e − → µ + e − (µ − e − → µ − e − ) process. The squared momentum transfers t e and t µ are defined in (11) and (12), respectively. The main message that can be drawn from Figure 13 is that the NLO QED correction over the full range is, in general, the result of a subtle interplay between the various sources of radiation. A further general remark is that the mixed corrections due to (a) dσ/dt e (b) dσ/dt µ Fig. 13 The contribution of the QED gauge-invariant subsets to the cross section of the process µ + e − → µ + e − (upper panel) and of the process µ − e − → µ − e − (lower panel), as a function of the squared momentum transfer t e and t µ . The results refer to Setup 1 (solid lines) and Setup 2 (dotted lines) described in the text. electron-muon interference are of opposite sign for the two processes and particularly relevant for large |t µ,e | values. The latter behaviour has to be ascribed to the presence in the updown interference of logarithmic (and squared logarithmic) angular contributions of the type ln(u/t), which become potentially enhanced when either t or u are small. More in detail, one can see from Figure 13 that, in the presence of acceptance cuts only, the NLO correction is dominated by the contribution of electron radiation, the other effects being almost flat and much smaller over the full range. However, if an acoplanarity cut is applied, the contributions due to muon radiation and up-down interference corrections become visible for large |t| values, where they amount to some percent. Interestingly, the above contributions have the same sign in the µ + e − → µ + e − process (upper panel of the left plot) and sum up to contribute to the overall QED correction, whereas they  tend to cancel in the µ − e − → µ − e − process (lower panel of the left plot). Therefore, also in view of ongoing calculations at NNLO, these results indicate that all the gauge-invariant subsets have to be taken into account. The size of finite electron-mass contributions at NLO is illustrated in Figure 14, where the results are shown in percent of the fully massive LO differential cross sections. The predictions refer to both incoming µ + and µ − and are shown for the electron scattering angle and the squared momentum transfer t e , for the sake of illustration. Similar results hold for other distributions. A sensible assessment of mass contributions beyond logarithmic accuracy is in general a delicate issue for any fixed-order calculation and particularly tricky for µ-e scattering under MUonE conditions, where the limit of massless electron implies that its rest frame (i.e. the lab frame) can not be defined. To bypass this difficulty, we follow the procedure detailed in [11], that allows to get an estimate of finite electron-mass contributions for the sum of one-loop virtual and real soft-photon corrections. The one-loop virtual amplitude is split into a contribution that is proportional to logarithms of the (artificial photon mass) IR parameter and the remainder, M (1) n = M (1) n,IR + M (1) n,non IR . A similar split is done for the amplitude related to soft real emission. According to the notation of Section 3, the m → 0 limit of the NLO correction is then evaluated according to the chain formula that provides an IR-safe estimate of electron mass contributions, while keeping exact kinematics and phase space. As can be seen from Figure 14, the contribution of mdependent terms to the LO cross section is almost flat and below the 10 ppm level. The electron-mass corrections at NLO contribute to dσ/dθ e and dσ/dt e in the range from a few to some 10 −5 . We notice that the largest part of the finite m corrections is due to radiation from the electron line only, the full correction lying around it. The extra corrections w.r.t. electron line only are dominated by up-down interference and box diagrams. These results suggest that electron mass contributions beyond logarithmic accuracy can be neglected in a NNLO computation or, eventually, included at the level of electron line corrections only.
The presently available MC at NLO accuracy represents just the first step towards the realisation of a high-precision theoretical tool necessary for the data analysis of µ-e scattering by MUonE. The ultimate goal is the realisation of a MC code including NNLO corrections and resummation of QED contributions due to multiple photon radiation. However, over a relatively short term, a number of intermediate results could be obtained about some important contributions beyond NLO.
A first example is given by the matching of NLO corrections to a QED parton shower, following the formulation already applied to Bhabha scattering and e + e − annihilation processes in QED [101,102], Drell-Yan processes [103,104] and Higgs decay into four leptons [105]. This would allow to estimate the most relevant QED corrections beyond NLO under realistic event selection criteria.
A further prospect under consideration is the calculation of lepton pair corrections to µ-e scattering. These corrections appear at NNLO and are a combination of two-loop virtual lepton-loop corrections with the same-order contribution of real pair emission, i.e. µe → µe + (ℓ + ℓ − ), with ℓ = e, µ. One and two-loop diagrams with vacuum polarisation insertions in the photon propagator were considered some time ago for the case of the Bhabha scattering in the massless limit with the 0.1% accuracy [88,[112][113][114][115][116].
Such a calculation can be extended to the treatment of hadronic pair corrections, by combining the already available virtual hadronic contributions [61] with the process of pion pair production µe → µe +(π 0 π 0 , π + π − ). A further step that can be taken is the evaluation of the background process µe → µe + (π 0 → γγ), that could benefit, as for two-pion production, from the experience in the development of MC generators for the simulation of hadronic final states at flavor factories [80]. Finally, one should not forget that electrons are bound inside the target and the impact of bound-state effects should be evaluated. This will require considering the possibility of scattering of the incident muons off core valence electrons, for which off-shell effects due to the finite binding energy and momentum distribution must be considered.
All the above developments are under consideration and preliminary results are also available for most of them.

Summary
With this report we want to document that within the theory community there is sufficient interest, manpower, and expertise to provide the necessary theory support for the MUonE experiment.
The minimal goal that will be achieved in a first step is a fully differential parton-level Monte Carlo program containing the following contributions: (i) the fully massive NLO QED (and electroweak) contributions; (ii) the fully massive (dominant) electronic contributions at NNLO; (iii) the fully massive NNLO hadronic contributions; (iv) the remaining contributions at NNLO in a massified approach, i.e. neglecting finite electron mass terms. In addition, this fixed-order calculation will be matched to a parton shower taking into account multiple photon emission at leading logarithmic accuracy.
Given the current status described in this report, there are no further conceptual challenges that need to be overcome to achieve this. Needless to say that nevertheless there will still be numerous difficult issues to be sorted. Hence, ideally there will be at least two different implementations of such a code to facilitate debugging. To that end, the theory groups of Pavia and PSI are both committed to each produce an implementation.
In a second step, a detailed, realistic phenomenological analysis is required to investigate if this theory description is sufficient. A careful error estimate of the missing terms will be crucial. This analysis will be done in close collaboration with the experimental collaboration.
It is quite likely that a third step will be required, i.e. further improvements to the theory. Most probably, the nextto-leading logarithms will have to be addressed. A careful study of the logarithms related to the factorisation anomaly is also important. In connection to fixed-order calculations, it is not at all unrealistic to expect a fully differential N 3 LO computation of the dominant electronic contributions in time for the MUonE experiment. Also, a complete fully-massive NNLO calculation, possibly using numerical techniques, is a serious target for theory in the longer term.
In fact, all the theory questions that are to be addressed in connection with µ-e scattering are also of interest to a much wider community. The developments that are made in this -from a theory point of view -simple framework will undoubtedly lead to progress in related fields. Thus, apart from providing an alternative determination of the HVP, MUonE can also act as an icebreaker to free a path for further theory progress.