Radiative Corrections: From Medium to High Energy Experiments

Radiative corrections are crucial for modern high-precision physics experiments, and are an area of active research in the experimental and theoretical community. Here we provide an overview of the state of the field of radiative corrections with a focus on several topics: lepton-proton scattering, QED corrections in deep-inelastic scattering, and in radiative light-hadron decays. Particular emphasis is placed on the two-photon exchange, believed to be responsible for the proton form-factor discrepancy, and associated Monte-Carlo codes. We encourage the community to continue developing theoretical techniques to treat radiative corrections, and perform experimental tests of these corrections.


Contents
1 Introduction Radiative corrections are of central importance for modern high-precision, sub-atomic physics experiments. In order to extract physical observables from experiments with percent precision, radiative corrections on the experimental observables need to be known to sub-percent precision. Modern tools to calculate the size, and to estimate the uncertainty, of these corrections are a subject of active research in the field. Frameworks for determining these effects can be developed and must be checked experimentally.
In modern experiments the radiative corrections for QCD and QED interactions are of particular interest. QCD systems (e.g., nucleons or light nuclei) are often probed in scattering experiments, which allow for precision extractions of form factors and structure functions. Typical probes are leptons or photons, since they have no internal structure themselves. Traditionally, the electron would be used as the lepton in scattering experiments as it is stable, but its low mass allows for significant external radiation. This becomes a challenge as experiments move to ever higher energies in the deeply inelastic scattering regime. In modern experiments, the muon serves as an attractive alternative to the electron, as its higher mass typically leads to smaller radiative corrections. In the case of QCD systems being studied via meson decay channels, hadronic and QED corrections can be of similar and significant size. In e + e − interactions, special care must be taken to not neglect the mass of the electron in order to achieve the requisite experimental and theoretical precision.
The topic of radiative corrections is quite broad, and as such in this review article we focus on a few selected topics. In Secs. 2 to 4, we focus on the process of lepton-proton (ℓp) scattering where ℓ = e, µ can be an electron or a muon. In Sec. 2, we start with the theoretical base and the classification of radiative corrections. In Sec. 3, we study an important class of radiative corrections -the two-photon exchange (TPE) contributions, which for hadron targets require knowledge of hadronic structure. In Sec. 4, event generators for scattering experiments are discussed. In Sec. 5, we consider higher-order corrections to the purely QED process of e + e − annihilation in the deeply inelastic regime. In Sec. 6, we discuss the corrections for radiative decays of light mesons. Finally, we summarize the state of the field, and give recommendations to those involved for a path forward to push the global understanding of radiative corrections to higher precision than exists today in Sec. 7.

Lepton-proton scattering
Scattering experiments allow experimenters to probe the structure of the target. From a measurement of the unpolarized ℓp scattering cross section, one can for instance determine the proton electric and magnetic Sachs form factors (FFs), G E (Q 2 ) and G M (Q 2 ), functions of the momentum transfer t = −Q 2 . The Born cross section of elastic ℓp scattering, given by the tree-level one-photon exchange (OPE) diagram in Fig. 1 (right panel), can be written in a compact form as with the Mott cross section for scattering off a pointlike particle the reduced Born cross section and the photon polarization parameter Here, we introduced the dimensionless quantities τ = Q 2 /4M 2 , ε ℓ = 1/(2τ ℓ ) and τ ℓ = Q 2 /(4m 2 ℓ ), and the crossing-symmetric variable ν ≡ k 1 · (p 1 + p 2 )/2 = (s−u)/4, with M (m ℓ ) the proton (lepton) mass, and where s and u denote the usual Mandelstam variables of the elastic scattering process of Fig. 1. Furthermore, α ≃ 1/137 is the usual fine structure constant, ϵ 1 (ϵ 2 ) is the energy of the incoming (outgoing) lepton, and θ is the scattering angle in the laboratory frame, see Eq. 24.
The extraction of the FFs, cf. Eq. 3, is done via Rosenbluth separation [1], which requires measurements at the same momentum transfer Q 2 but for different energies and angles. To recover the Born cross section, the experimentally measured cross section needs to be corrected for radiative corrections, cf. Refs. [2][3][4][5][6][7]. These change not only the absolute value of the cross section but also its dependence on the relevant kinematical variables. In general, the O(α 2 ) Born cross section has to be corrected by additional virtual photons or inelastic scattering with emission of real bremsstrahlung.

Radiative corrections
In this subsection, we discuss the calculation of cross sections in a perturbative expansion of the electromagnetic coupling α. While many remarks are valid for arbitrary processes, we will illustrate the general procedure for elastic ℓp scattering where k 1 (k 2 ) and p 1 (p 2 ) are the incoming (outgoing) lepton and proton momenta, and h(h ′ ) and λ(λ ′ ) the respective helicities. The Born or leading order (LO) cross section for elastic ℓp scattering, given in Eq. 1, is of O(α 2 Z 2 ) and is a function of G E (Q 2 ) and G M (Q 2 ). These form factors are depicted as grey blobs in Fig. 2 and, for each appearance, a factor Z, the total charge, is included in the counting of the couplings. To improve the theory description, next-to-leading order (NLO) corrections to Eq. 1 have been computed in Refs. [2][3][4][5][6][7]. The NLO corrections can be decomposed into several gauge-invariant parts. The one-photon exchange (OPE) contributions are of the form O(α 3 Z 2 ) and also include vacuum polarisation (VP) corrections. They are illustrated in Figs Fig. 7. These corrections change not only the absolute value of the observables but also their dependence on the relevant kinematical variables, e.g., the square of the transferred momentum Q 2 and the photon polarization parameter ε, defined in Eq. 4.
If the radiative corrections are small, an NLO calculation might be sufficient. In particularly simple situations, these corrections can be included as a multiplicative factor. The error associated to this procedure is at the percent level. Radiative corrections can lead to contributions that are enhanced by large logarithms of the form L = ln (Q 2 /m 2 ℓ ). As we will discuss below, these are associated with hard collinear radiation. Another source of large logarithms are stringent experimental cuts that restrict the phase space for soft emission. As a consequence, modern high-precision experiments are required to introduce higher order corrections and a more sophisticated treatment of real radiation.
This leads us to consider fully differential computations of cross sections which allows one to obtain   precise predictions for fiducial cross sections. To keep full generality, the phase-space integration needs to be adaptable to the experimental setting. As a result, it is necessary to revert to numerical integrations, which are most efficiently carried out using Monte Carlo (MC) techniques. For the LO cross section dσ (0) ∼ α 2 Z 2 , the tree-level matrix element squared |A (0) | 2 , depicted in Fig. 2, is to be integrated over the two-body phase space Φ 2 . Thus, where the so-called measurement function S(k 2 , p 2 ) defines the quantity to be computed in terms of the   momenta of the final state particles. For simplicity, trivial elements such as the flux factor are also included in S. If a cut is to be applied, S is set to 0 if the event does not pass the cut. Going beyond LO, we encounter ultraviolet (UV) and infrared (IR) divergences. After renormalization, in a QED calculation where all fermion masses are kept at their finite values, only soft IR divergences remain. These soft singularities cancel between the virtual and real corrections. At NLO this is illustrated in Fig. 3, where some corrections due to photon emission from the lepton line are shown. The virtual corrections are obtained as dσ (1) v = dΦ 2 2 Re(A (1) A (0) * ) S(k 2 , p 2 ) , whereas the real corrections are given by Both terms are O(α 3 Z 2 ). Here A (1) is the oneloop amplitude of the process Eq. 5 and A (0) γ is the tree-level amplitude with an additional photon of momentum k γ in the final state. The integration in Eq. 8 has to be carried out over the three-particle final state Φ 3 . Also, the definition of q is more subtle now, since k 1 − k 2 ̸ = p 2 − p 1 .
To give mathematical meaning individually to dσ (1) v and dσ (1) r , the IR singularities need to be regularized. Many calculations in QED still use a (small) fictitious photon mass to do so. Then, the IR singularities manifest themselves as logarithms of this mass. Another approach that is more in line with the tremendous progress in computational techniques that has been made for QCD calculations is to use dimensional regularization (i.e. work in 4 − 2ϵ dimensions) also for IR singularities. In this case, IR singularites show up as poles 1/ϵ. See Sec. 3.2.1 for a detailed discussion. Independent of the chosen regularization, in the combination of dσ (1) v and dσ (1) r the remnants of the IR singularities (i.e. the photon mass logarithms or the 1/ϵ poles) cancel. This is ensured by the well-known limiting behavior of the real matrix element in the soft limit (S@) [8], where it can be written as an eikonal factor E times the Born matrix element S@|A (0) γ | 2 = E |A (0) | 2 .
(9) Hence, in the complete NLO corrections to Eq. 6, dσ (1) = dσ (1) v + dσ (1) r , the regularization can be removed (i.e. the photon mass or ϵ set to zero). For this cancellation to occur, the observable has to be IR safe. This means that the observable must not change whether or not an arbitrarily soft photon is emitted. At a technical level, we must require From a mathematical point of view, it is sufficient if Eq. 11 holds in the strict limit. However, if experimental cuts (directly or indirectly affecting real radiation) induce a difference between S(k 2 , p 2 , k γ ) and S(k 2 , p 2 ) for finite, but small photon energies ∆E γ , there can be a large logarithm ln ∆E 2 γ /Q 2 as a remnant. Such (soft) logarithms can lead to enhanced QED corrections. Another source of large corrections is collinear emission. While collinear singularites are regulated by the fermion masses, they still can also lead to enhanced logarithmic terms αL after combination of real and virtual corrections. These logarithms often form the dominant part of the higher-order corrections. This also entails that in fully differential QED calculations -contrary to QCD calculationsit is not possible to set the lepton masses to zero. Furthermore, the neglect of hard collinear emission potentially leads to a loss of accuracy. Thus, using the soft approximation for the real matrix elements can have severe implications on the accuracy of the results.
The extraction of the IR poles from Eq. 8 for an arbitrary IR-safe observable, i.e. an arbitrary S(k 2 , p 2 , k γ ) satisfying Eq. 11, is by now a standard procedure at NLO and the focus turns towards next-to-next-to leading order (NNLO) applications as discussed in Ref. [9]. Again, at NLO there are two widely used options, the slicing method and the subtraction method.
In the slicing method, the phase space is split into two parts, depending on the photon energy In the part, where the photon energy is larger than a chosen resolution parameter δ, the integration can be carried out for a massless photon without encountering singularities. In the part, where the photon energy is smaller than the resolution parameter, the soft (eikonal) approximation Eq. 9 is used for the matrix element. The integration over the photon phase space then simplifies and can be carried out analytically. In this term, the photon mass logarithms that cancel those of the virtual corrections are generated. This method relies on a suitable choice for δ. It should not be too large, to ensure that the soft approximation is good enough. It should also not be too small, otherwise numerical problems (large cancellations between the resolved and unresolved parts) occur. For more details, see Sec. 4. In the subtraction method a term (the soft limit of the integrand) is subtracted and added back The first term on the r.h.s. of Eq. 13 is very similar to the corresponding term in Eq. 12. The only differences are that the photon is treated as massless and the integration is done over the full phase space in Eq. 13. This results in IR 1/ϵ poles that cancel against the IR poles of the virtual contributions. The second term on the r.h.s. of Eq. 13 is finite and can be integrated in four dimensions. An efficient numerical evaluation of this term including the subtraction requires an adapted phase-space parameterization. Within the subtraction method, there is no motivation to split the real corrections into hard and soft parts. It is actually simpler to always include the full real matrix element and make no approximation whatsoever. The OPE corrections due to emission from the lepton line O(α 3 Z 2 ) and the VP corrections, illustrated in Figs. 3 and 4, are simple from a conceptual point of view as they involve standard QED pointlike interactions only. However, there are also corrections of O(α 3 Z 3 ) that involve multiple exchange of photons between the two fermion lines, in particular the notorious TPE corrections. The class of corrections depicted in Fig. 5 is more complicated than those of Fig. 3 for several reasons. Even if the proton was pointlike, these corrections now involve one-loop box diagrams at NLO. In addition, there is the complication that the photon-proton vertex is more involved than for a pointlike particle. Furthermore, for ℓp scattering there are also contributions where the intermediate state is not just a proton, as depicted in Fig. 6. These aspects will be discussed in Sec. 3.
A final class of NLO corrections are those where the emission of additional photons is restricted to the proton line. They are O(α 3 Z 4 ) and a sample contribution is depicted in Fig. 7. It is tempting to include these diagrams into the definition of the form factors. However, they contain IR divergences that, once more, cancel between the real and virtual parts. Thus, from a practical point of view it is more convenient to explicitly include these diagrams into the radiative corrections and define the form factors accordingly, i.e. exclude all QED effects from their definition. The numerical impact of these contributions is rather minor, as they do not generate collinear logarithms. Hence, it is often possible to simply neglect the O(α 3 Z 4 ) contribution.
As discussed above, at each order in α the QED corrections can be enhanced by a collinear logarithm L and possibly a soft logarithm ln ∆E γ 2 /Q 2 . Hence, the corrections can be considerably larger than naive scaling with α implies. Thus, for very precise predictions it is necessary to go beyond NLO. Regarding OPE corrections restricted to emission from the lepton line O(α 4 Z 2 ), the computations have been extended to NNLO. This involves doublevirtual, real-virtual, and double-real corrections, as illustrated in Fig. 8. Results have been obtained using photon-mass regularization and the slicing method in Refs. [10,11] as well as using dimensional regularization with the subtraction method in Ref. [12]. This class of corrections is a gauge invariant subset of the full corrections. The virtual part involves only the vector part of the two-loop heavy-lepton form factor, cf. Refs. [13][14][15]. For lepton masses much smaller than the proton mass, these corrections are typically also dominant. They involve up to two additional photon emissions and a precise description how these photons are treated is required. This amounts to a definition of S(k 2 , p 2 , k γ1 , k γ2 ). In particular, both photons can be collinear, leading to enhanced NNLO correction of the form (α L) 2 .  Corrections at NNLO that go beyond the OPE approximation lead to substantially more complicated calculations. If the proton is treated as pointlike, the NNLO results for electron-muon scattering in Refs. [16,17] can be adapted by simple replacements m µ → M and m e → m ℓ . This requires the computation of the two-loop amplitude A (2) (for double-virtual) involving (crossed) double-box diagrams [18]. For the real-virtual corrections A (1) γ is required where one-loop pentagon diagrams contribute. Efficient techniques for one-loop computations have been developed in connection with high multiplicity processes at hadron colliders, as reviewed in Ref. [19]. As part of this topical collection, the effect of pointlike NNLO corrections were compared to present uncertainties from hadronic correction [20]. It turns out that the foremer can be sizeable and need to be properly taken into account. Beyond pointlike photon-proton interaction make the formidable task of the two-loop amplitude computation even more daunting. At NNLO, depending on the experimental set up, it might also be necessary to include the process with emission of an additional lepton pair which also has been computed for electron-muon scattering in Ref. [21].
First steps to move beyond NNLO have already been taken. The form factor γ * → ℓ + ℓ − is known at three loops [22] and the subtraction method for massive QED has been extended beyond NNLO [23]. Hence, a next-to-next-to-next-to leading order (NNNLO) calculation of OPE contributions O(α 5 Z 2 ) seems to be feasible within a reasonable time frame. However, a full computation of 2 → 2 scattering processes with massive particles at NNNLO is currently out of reach.
As alluded to above, the corrections are often dominated by large logarithms associated with collinear or soft emission. A convenient method resumming all contributions of large logarithms at all orders has been proposed in Ref. [24]. All leading contributions of order (αL) n are taken into account. The accuracy can be further improved calculating nonleading contributions of order α(αL) n as a K-factor. The application to ep elastic scattering can be found in Ref. [25], and, including hard photon emission, in Ref. [26].
A more generic approach to include logarithmically enhanced corrections is to use QED parton showers. A recent review of this activity can be found in Ref. [27] and will be discussed in some more detail in Sec. 4. Parton showers allow for the inclusion of the leading logarithmic contributions in a process independent way. Combined with fixed-order calculations this is a powerful tool to improve further the precision for fully differential cross sections. Including next-to-leading logarithms from initial-state (or finalstate) collinear emission is also possible in a generic way [28]. However, in general the systematic inclusion of subleading logarithms is observable dependent and some examples will be discussed in Sec. 5. Contrary to QCD, in QED it is typically not required to actually resum the logarithms. The suppression by α is strong enough to ensure that α L is still reasonably small. Hence, the inclusion of the logarithms at the first few orders in α beyond the fixed-order results is sufficient.

Experimental observables
In the following, we discuss different experimental observables and how they are affected by the LO in α virtual corrections, cf. Eq. 7. Of particular interest are the unpolarized cross section, the longitudinal and transverse polarization transfer observables P t and P l , as well as the target and beam asymmetries A n and B n , introduced in Secs. 2.2.4 and 2.2.5 respectively.

Elastic lepton-proton scattering
Let us start by presenting the tensor decomposition of elastic ℓp scattering. Taking into account discrete symmetries, assuming parity and time reversal invariance of the electromagnetic interaction, this process can be described in terms of 6 independent complex scalar amplitudes: G M and F i (i = 2, . . . 6), which are generally functions of two independent kinematical variables, e.g., the Mandelstam variables t = (k 1 − k 2 ) 2 and s = (k 1 + p 1 ) 2 . The helicity amplitudes of the process can be split into a sum of helicityconserving and helicity-flip contributions, see, e.g., Refs. [29][30][31], where q = k 1 − k 2 = p 2 − p 1 , K = (k 1 + k 2 )/2, P = (p 1 + p 2 )/2, and γ · a ≡ γ µ a µ . The 6 complex amplitudes G M and F i , sometimes called generalized FFs, fully describe the spin structure of reaction Eq. 5 for any number of exchanged virtual photons. In the limit of m ℓ → 0, the contribution from the helicityflip amplitudes to observables vanishes, see Eq. 16. Relations between the helicity amplitudes and the generalized FFs can be found f.i. in Ref. [32].
In the Born approximation, i.e., considering only the tree-level OPE diagram with proton FFs in Fig. 1 (right panel), the number of amplitudes in Eq. 14 reduces to 2. To be more precise, the Born amplitudes read where F 2 (Q 2 ) is the Pauli FF of the proton. The fact that the FFs in the space-like region are real functions of the virtuality of the exchanged photon is a consequence of unitarity, the spin-1 nature of the virtual photon, parity conservation, and the identity of the initial and final states, see Refs. [33,34] for an explicit derivation.
To separate the Born contribution from effects of virtual radiative corrections to the elastic scattering, we define the following decomposition of the scalar amplitudes: (22) where G E has been introduced via The order of magnitude of these quantities is given by

Laboratory frame
The proton and lepton four-momenta in the laboratory frame can be written as where three momenta are denoted by bold symbols. The scattered lepton energy is written in terms of θ as while the Mandelstam variables s and t are expressed as The differential cross section in terms of the matrix element squared is given by with the invariant I 2 = (k 1 · p 1 ) 2 − m 2 ℓ M 2 and the energy of the recoil proton E 2 .
For the case where the scattered lepton is detected in the final state, one obtains where D = (M + ϵ 1 )|k 2 | − ϵ 2 |k 1 | cos θ, and dΩ is the differential solid angle of the scattered lepton. For the case where the recoil proton is detected in the final state, one obtains where θ p is the angle between the directions of the lepton beam and the recoil proton, and dΩ p is the differential solid angle of the scattered proton. Using the relation and ε T is the degree of linear polarization of transverse photons In this way, one can also describe the leading TPE correction, cf. Eq. 56.

Polarization transfer observables
The longitudinal and transverse polarization transfer asymmetries, P l and P t , are defined as (with, e.g., h = +) where S ′ = ±S ⊥ is the spin direction of the recoil proton in the scattering plane transverse to its momentum. In the Born approximation, their ratio is related to the ratio of electric to magnetic Sachs FFs [35,36] Equivalent information on the proton FFs can be obtained from measuring double-spin asymmetries on a polarized proton target [37]. Taking into account the leading virtual corrections to the asymmetries, the ratio modifies, up to terms of order O(α 2 ), as described in Refs. [29,38,39] as with where the shorthand notations for ratio of two-photon amplitudes relative to the magnetic FF have been used: Furthermore for the longitudinal polarization transfer P l separately, its expression relative to the 1γ-result P Born l is given in Refs. [29,38], up to terms of order O(α 2 ), by with P Born Note that in Eqs. 45 and 50 the lepton mass has been neglected, m ℓ = 0, which is not a valid approximation for low-energy muon scattering [40].

Beam and target normal single-spin asymmetries
The target (beam) normal single-spin asymmetries (SSA) A n (B n ) is defined for the scattering of an unpolarized beam (a beam polarized normal to the scattering plane, with s = ±s n the direction of incoming lepton,) off a target with normal to the scattering plane polarization S = ±S n (an unpolarized target) Both these asymmetries are zero in the Born approximation. Taking into account leading virtual corrections, the asymmetries can be expressed as [30,41,42] where Im denotes the imaginary part.
3 Two-photon exchange The TPE contribution to lepton scattering, shown in Fig. 11, is of particular importance for two main reasons: Firstly, hadronic corrections, cf. the blob in Fig. 11 (bottom panel), are notoriously difficult to calculate and often have a large relative uncertainty. Secondly, as will be explained in Sec. 3.1, the OPE may not be a good approximation for the extraction of FFs from the unpolarized cross section at large Q 2 . Therefore, the leading TPE correction following from the interference of the Born and TPE amplitudes, A (0) and A 2γ , has to be taken into account appropriately. A complete calculation should go beyond the soft-photon approximation and include the hard TPE with all possible intermediate states.
In the last decade, predictions of the TPE correction have been systematically improved and model dependence has been reduced by employing dispersion relations and effective field theories. This part of the paper is organized as follows. In Sec. 3.1, we illustrate the importance of the TPE by comparing Rosenbluth and polarization-transfer measurements of the proton FFs. In Sec. 3.2, we review theoretical predictions for the leading TPE correction, distinguishing between proton and inelastic intermediate states, as well as regions from small to large momentum transfer. In Sec. 3.3, we discuss effective field theory calculations. In Sec. 3.4, empirical extractions of TPE corrections and amplitudes are presented and updated based on new data. In Sec. 3.5, results of past experiments aimed at extracting the TPE are presented and compared to theoretical predictions. We finish with an outlook on future experiments and theory advances in Sec. 3.6. For further reading, we refer to the following reviews in Refs. [43][44][45][46] and the recent CFNS whitepaper, Ref. [47].

Rosenbluth vs. polarization transfer experiments
Polarization experiments in elastic ep scattering at Jefferson Lab Hall A [48][49][50][51], revealed that the ratios of proton FFs, G E (Q 2 )/G M (Q 2 ), extracted based on the Rosenbluth [1] or polarization transfer methods [35,36] in the OPE approximation deviate with increasing Q 2 . The FF ratio from Rosenbluth extractions is nearly constant as a function of Q 2 , while it decreases for polarization transfer measurements, as can be seen in Fig. 9. In the limit of OPE, one defines an unpolarized reduced cross section that is linear in ε and the corresponding ratio of measurements with transverse or longitudinal polarization of the recoiling proton is constant in ε as given in Eqs. 3 and 44. What is measured in experiment, however, will always include radiative corrections. Taking these into account and under some assumptions, in presence of TPE, the expressions modify according to Eqs. 36 and 45, respectively. One can see that after inclusion of the TPE correction, the interpretation of what is measured in Rosenbluth and polarization transfer experiments changes. Model independent considerations [52][53][54] show that it is still possible to recover experimentally the electric and magnetic proton FFs, even in presence of TPE, but this requires either the measurement of three time-odd or five time-even polarization observables (including triple spin observables, of the order of α) or, alternatively, the generalization of the Akhiezer-Rekalo recoil proton polarization method [35,36] with longitudinally polarized electrons and positron beams in identical kinematical conditions. In a series of papers in the early 2000's, it has been shown that the inclusion of "hard" TPE may reconcile the FF extractions from polarized and unpolarized ℓp scattering, within some assumptions (real-valued FFs and ε linearity of the TPE contribution). Guichon and Vanderhaeghen identified in Ref. [29] the experimental observables used to extract the G E /G M ratio as the ε-slope in the reduced cross section Neglecting lepton mass terms and using the shorthand notations for the TPE amplitudes of Eq. 49, δ 2γ has the form On the other hand, the polarization transfer experiments yield the ratio given by Eq. 45 Reference [29] showed that the TPE in the Rosenbluth method effectively corrects a small number, cf. the ε dependent term in Eq. 57 with R 2 EM (0) ∼ (1/2.79) 2 ∼ 0.13 which has an additional 1/τ suppression at large Q 2 . The FF ratio extracted from P t /P l is only mildly affected by radiative corrections as it involves asymmetries, whereas the cross sections themselves are strongly modified by radiative corrections at moderately and large momenta. Around the same time, Blunden [2] (no hard TPE was included), whereas for the polarized extraction no radiative corrections were taken into account.
Note that alternative explanations have been put forward, such as different calculations of radiative corrections, including also lepton structure functions [25,68], correlations between the Rosenbluth parameters [69], or acceptance problems in the analysis or experiment [70].
The electric FF is determined by the slope of the Rosenbluth plot, at fixed Q 2 , derived from the radiatively corrected cross section. Corrections at large Q 2 may reach 40% [59] and the ε-slope may essentially change. Even more dramatically, above 3 GeV 2 the slope of the measured cross section may even be negative before corrections, cf. Fig. 10 [71,72] where the effect of radiative corrections [3], including "soft" but no "hard" TPE, is shown. Note that G E must be real in the space-like region. This means that G 2 E is essentially determined by the ε-dependence of the applied radiative corrections. Different calculations show a difference of few percent in the slope, already at NLO [5,6].  [59] at momentum transfers of 1.75 (green squares), 3.25 (blue triangles) and 5 GeV 2 (red circles). Empty and full symbols show data before and after radiative corrections [3] (no hard TPE included), respectively. Figure taken from Ref. [72].
Concerning polarization observables, it is commonly assumed that the FF ratio given by polarization measurements is not (or less) affected by radiative corrections, giving therefore a more reliable result on this ratio than the one extracted by the Rosenbluth method. The polarization asymmetries P t and P l are ratios of cross sections, cf. Eqs. 42 and 43, in which the bulk of the radiative corrections that factorize (virtual corrections on lepton side and soft-photon emission corrections) drop out. Therefore, they are less sensitive to radiative corrections than the unpolarized cross section.
Already in the 1970's, a reason why the TPE correction could become important at large Q 2 was put forward [73][74][75]: when the transferred momentum is equally shared between the two photons the scaling in α may be compensated by the steep decrease of the FFs with Q 2 . The effect is therefore expected to increase with Q 2 and with the hadron mass. The advent of the high duty cycle and highly polarized electron beam at the Jefferson Lab, together with the availability of large solid angle spectrometers and hadron polarimetry in the GeV region, opened the way to high-precision experiments in elastic and inelastic electron-hadron scattering. Two experiments of elastic electron-deuteron (ed) scattering at large Q 2 in Hall A [76] and Hall C [77] claiming an error of 5%, showed a discrepancy of up to 15% in the elastic ed cross section at the same Q 2 , but different energies and angles. A possible explanation was brought up that the TPE contribution could be at the origin of these findings [78]. Finally the discrepancy was attributed to a systematic error of the Hall C spectrometer setting, as no Q 2 dependence was observed.

Theoretical predictions
In this subsection, we review theoretical predictions for the TPE contributions to ℓp scattering. We will refer to the (standard) TPE approximations by Maximon-Tjon [2] (or Mo-Tsai [3]), conventionally included in the experimental analysis, as "soft" TPE. (Non-standard) refinements of the TPE, that need to be included given the accuracy of modern scattering experiments, will be referred to below as "hard" TPE. Starting from the so-called elastic TPE with a proton intermediate state, Fig. 11 (top panel), we subsequently discuss TPE contributions with inelastic intermediate states, Fig. 11 (bottom panel), in regions of small, medium and large momentum transfer. Furthermore, we discuss the forward TPE as the limiting factor in the theoretical description of the energy spectra of light muonic atoms.

Infrared subtraction schemes
When plotting the elastic TPE contribution to δ 2γ or R 2γ , defined in Eqs. 56 and 90, one has to somehow treat the IR divergences they contain. Of course the only way that leads to a measurable quantity is to combine them with the corresponding real corrections as discussed in Sec. 2.1. However, this can be quite cumbersome in practise and makes comparing different TPE scenarios and parametrisations needlessly complicated. Instead, one defines an IR subtraction scheme that removes the IR divergence (and potentially some finite parts as well). The exact form of this subtraction term is somewhat arbitrary and solely determined by conventions since the result is no longer a physical quantity. A common way of doing this in the high-energy community is the definition of the Catani I operator [79] that simply removes the 1/ϵ poles in dimensional regularisation.
In the QED community it is not uncommon to add the first terms of Eqs. 12 or 13 This scheme, sometimes called eikonal subtraction [23], has the added advantage of giving the remnant a physical interpretation as long as the parameter δ is chosen small enough -it just approximately includes real corrections up to the cut-off δ.
The calculation of this integral is somewhat involved but the result is well known [23,80,81]. Since this physical interpretation exists, it also means that the result needs to be independent of how the TPE calculation was regularised, making comparisons between different methods easier. The TPE community usually follows either Mo-Tsai (MoT) [3] or Maximon-Tjon (MTj) [2]. Note that in this paper, we use the MTj convention for our Figs. 12-19. The latter defines the IR-subtracted TPE contribution as This integral can be evaluated both using photonmass regularisation or dimensional regularisation where we have defined λ s = m 4 + (M 2 − s) 2 − 2m 2 (M 2 + s). In the limit of m 2 ≪ s, M 2 , Q 2 , the expression for e − p scattering simplifies to To derive the MoT subtraction prescription, we instead start by replacing either the ℓ 2 propagator or the (ℓ − q) 2 propagator with 1/q 2 . The idea is to mimic the soft behaviour of the integral as either the ℓ 2 propagator goes to zero or the (ℓ − q) 2 propagator. This procedure results in two identical triangle functions Note that this is not yet the correct MoT prescription, hence MoT ′ . The loop integral is slightly more involved than the previous one but can still be trivially written in terms of the well-known Ellis-Zanderighi functions [82]. In this case, we need the integral I fin 6 (s; m 2 , M 2 ) commonly referred to as Triangle 6, If one chooses, I fin 6 could be expressed in terms of logarithms and dilogarithms. The resulting expression is not overly complicated but reproducing it here serves no practical purpose. However, in the limit of m 2 ≪ s, M 2 , Q 2 the expression is fairly compact where we have defined y = M 2 /(s−M 2 ) and used the usual definitions of the ζ function and dilogarithm.
To arrive at the correct MoT prescription [3] (cf. also [44]) we need to flip k 1 → −k 1 in the first term of (64), resulting in This operation is equivalent to setting s → s ′ = 2m 2 + 2M 2 − s, which spoils the crossing symmetry but eliminates the (iπ) 2 from the analytic continuation of the logarithm. The resulting expression can still be written in terms of I fin 6 but now in terms of s ′ rather than s The symmetry of the λ functions means that λ s ′ = λ s , which helps to simplify the result. The difference between δ MoT Crucially, the poles agree after taking the real part so that the IR cancellation in δ 2γ is unaffected.

Proton intermediate state
The first terms in the low-Q expansion of the TPE with proton intermediate state can be described model independently. The leading term is given by the well-known Feshbach correction [83] which stems from the interaction between a massless lepton and a structureless proton. Inelastic intermediate states start to contribute in the subleading Q 2 ln Q 2 term only [84], see discussion in Sec. 3.2.4. It follows that TPE could become larger in lepton scattering on heavy ions, inducing a large charge asymmetry even at small angles [85]. During the last two decades, predictions of the TPE contribution improved considerably. The simplest approach to evaluate the elastic TPE is with a hadronic model calculation of the box (and crossedbox) diagrams in Fig. 11 (top panel), assuming on-shell proton FFs, see for instance Refs. [55] and [32] for electron and muon scattering, respectively. A way to avoid model dependence is to instead use a dispersive approach to describe the scalar amplitudes G i and F i , introduced in Eqs. 14 and 37. Dispersion relations (DRs) allow one to express the real part of the amplitudes via their imaginary part, and the latter, using unitarity, can be related to physical observables. The s−channel cut in the elastic TPE diagram then permits the use of on-shell proton FFs. First dispersive evaluations of the elastic TPE in ep scattering, neglecting the electron mass, used unsubtraced DRs for the scalar amplitudes [86,87]. Considering the case of µp scattering, in which the lepton mass cannot be neglected, F 4 will require a once-subtracted DR [88].
The contribution from TPE with inelastic intermediate states, discussed in the following, is becoming important with increasing momentum transfer. It has been suggested that further subtractions, e.g., in the DR of the F 3 amplitude, can be introduced in order to minimize model-dependence from higher intermediate states [89]. Sum rules, which are exact in QED and approximated in QCD, give indications on the contribution of intermediate states [90].

Zero momentum transfer
In the limit of zero momentum transfer, Q 2 = 0, the TPE is not important as a radiative correction to the scattering process. δ 2γ , as defined by the interference of OPE and TPE amplitudes in Eq. 56, is vanishing in the forward limit (Q 2 → 0) at fixed ν, see discussion in Ref. [88]. At the next order in α, the contribution from |A 2γ | 2 to the cross section at zero momentum transfer is not vanishing, however, numerically suppressed. On the contrary, the TPE in forward kinematics is important in the description of the spectra of light muonic atoms, where its effect is enhanced as compared to ordinary atoms due to the heavier muon mass. Let us focus on (muonic) hydrogen. The forward TPE corresponds to a δ(⃗ r ) potential that gives an O(α 5 ) contribution to the energy levels. For comparison, the leading contribution from the Coulomb potential is of O(α 2 ) and the proton finite-size starts to contribute at O(α 4 ) through an OPE diagram with proton FFs. The off-forward TPE starts to contribute at O(α 6 ln α) through the so-called Coulomb distortion effect. See Ref. [91] for a comprehensive theory of the Lamb shift in light muonic atoms. Here, we want to focus on the dispersive formalism used to evaluate the O(α 5 ) effect of the forward TPE, since a similar approach is used to approximate the small momentum transfer corrections to the scattering process, cf. Sec. 3.2.4.
The forward TPE between a lepton and a proton can be expressed in terms of forward doubly-virtual Compton scattering (VVCS) off the proton. The TPE-induced shift of the nS-level in a hydrogen-like atom is given by the unpolarized VVCS amplitudes T 1 and T 2 [92] ∆E 2γ (nS) = 8παm ϕ 2 where ϕ 2 n = 1/(πn 3 a 3 ) is the wavefunction at the origin, a = (Zαm r ) −1 is the Bohr radius (in the following Z = 1 for hydrogen), m r is the reduced mass of the lepton-proton system and m is the lepton mass (m e or m µ , respectively, for hydrogen and muonic hydrogen). Furthermore,ν = q 0 andQ 2 = q 2 − q 2 0 are the energy and virtuality of the VVCS photons inside the TPE loop diagram. Similarly, the TPE contribution to the hyperfine splitting can be expressed through the spin-dependent VVCS amplitudes S 1 and S 2 . The TPE effect can then be evaluated with a datadriven dispersive approach [92][93][94][95][96][97][98][99][100][101], or predicted from variants of chiral perturbation theory [102][103][104][105][106][107][108], and in the future, from lattice QCD [109,110]. For recent reviews, discussing low-energy proton structure in the context of muonic-hydrogen spectroscopy and lepton scattering, including TPE and the proton radius, we refer to Refs. [43,[111][112][113][114][115].
Data-driven evaluations of the forward TPE make use of dispersion relations for the VVCS amplitudes, expressing them through proton structure functions measured in electron scattering where F 1 (x,Q 2 ) and F 2 (x,Q 2 ) are the unpolarized structure functions, with x =Q 2 /2Mν the Bjorken variable, andν el =Q 2 /2M . Similarly, the spindependent VVCS amplitudes S 1 and S 2 follow from dispersion integrals over the spin structure functions g 1 (x,Q 2 ) and g 2 (x,Q 2 ). Unfortunately, the T 1 amplitude in Eq. 72 requires a once-subtracted dispersion relation. Since the subtraction function T 1 (0,Q 2 ) cannot be fully constrained from experiment, the data-driven approach suffers from a model dependence.
The TPE contributions to the energy levels in muonic atoms are usually split into a so-called "elastic" contribution, corresponding to the TPE with proton intermediate state shown in Fig. 11 (top panel), and a so-called "polarizability" contribution, corresponding to the diagrams in Fig. 11 (bottom panel). The elastic part of the VVCS amplitudes is well-known and can be expressed through proton FFs. The uncertainty of the elastic TPE can therefore be reduced with improved input for the proton FFs. This has recently been done based on the FF descriptions from Refs. [116,117], which are both consistent with the small proton charge radius extracted from the muonic-hydrogen Lamb shift, see Refs. [118] and [91] for the elastic TPE in the muonic hydrogen Lamb shift and hyperfine splitting, respectively. Therefore, the dispersive approach is only used to reconstruct the inelastic part of the VVCS amplitudes. To this end, the dispersion integral is restricted to the inelastic structure functions, i.e., the region of x ∈ [0, x 0 ] with the inelastic threshold x 0 , e.g., the pion production threshold. In the following section, we will discuss the extension of this formalism to the TPE correction to ℓp scattering in the region of small momentum transfer.

Small momentum transfer
In the near-forward kinematics, i.e., in the small momentum transfer region, the TPE amplitude can be approximated with a modification of the dispersive approach presented above, see Eqs. 72 and 73, that reconstructs the VVCS amplitudes based on empirical input for the inelastic proton structure functions. In the following, we review the leading-Q 2 behaviour of the inelastic TPE corrections to the unpolarized cross section and the beam normal spin asymmetry as functions of the transverse cross section for photoabsorption off the proton. Recall that the leading terms in the Q 2 expansion of the TPE stem from the proton intermediate state [83,84], see discussion in Sec. 3.2.2.

Unpolarized cross section
The leading-Q 2 behaviour of the inelastic TPE correction to the unpolarized cross section reads [84,119] where W 2 is the invariant mass of the hadronic state in the internal Compton scattering process, andν π is the pion production threshold. This formula has been first published by Brown [84] in the early 1970's. Later, Eq. 74 has been re-evaluated by Gorchtein [119] based on the phenomenological fit of the total photoabsorption cross section in Ref. [120]. It is important to point out that the coefficient of the order Q 2 ln Q 2 term in δ 2γ is a model-independent result. Tomalak and Vanderhaeghen have extended this approach beyond the leading Q 2 ln Q 2 term and approximate the δ 2γ correction for ep scattering in the region up to Q 2 = 0.25 GeV 2 in Ref. [121] and for µp scattering in the region up to Q 2 = 0.08 GeV 2 in Ref. [122], see Fig. 19. As can be seen from Eqs. 72 and 73, in general, this requires input for the unpolarized proton structure functions, F 1 (x,Q 2 ) and F 2 (x,Q 2 ) (or equivalently, the transverse and longitudinal cross sections, σ T (ν,Q 2 ) and σ L (ν,Q 2 )), as well as the subtraction function T 1 (0,Q 2 ), where contrary to Eq. 74 theQ 2 dependence of the input is kept. Besides the limited applicability range of the near-forward approximation, the T 1 (0,Q 2 ) subtraction function introduces an additional model dependence for the scattering of massive leptons, while its contribution in the elastic ep scattering is suppressed by the electron mass and negligible [121].
Beam normal spin asymmetry The leading-Q 2 behaviour of the inelastic TPE contribution to the beam normal single spin asymmetry in the diffractive limit (high-energy and forward scattering) reads [123,124] where the total photoabsorption cross section is assumed to be independent of the photon virtuality and roughly constant as a function of the invariant mass of the photon-proton system. In the nuclear resonance region, where σ T strongly varies with energy, the cross section enters through an energy-weighted integral. Note that this result is not fully model independent, as it makes use of the Callan-Gross relation [125] between the longitudinal and the transverse cross section in the DIS region. The expression Eq. 75 has different dependence on momentum transfer Q and the beam energy ϵ 1 compared to contributions from the proton intermediate state. First, it does not fall off with the beam energy as m/ϵ 1 ; second, it has linear dependence on θ at a fixed energy ϵ 1 in the near-forward kinematics, whereas for elastic intermediate states (as well as for a celebrated Mott formula) the asymmetry falls as the third power, ∝ θ 3 . As a result, the contribution of inelastic intermediate states becomes dominant for energies above GeV [123,124,126]. In Ref. [127], subleading terms in the Q 2 expansion of B n were derived and the energy dependence of the photoabsorption cross section was taken into account. A comparison to Eq. 75 showed that the leading-term in the Q 2 expansion is indeed a good approximation. Note that while B n also has a double-logarithmic enhancements due to hard collinear quasi-real photons [123,128,129], it is highly suppressed for small scattering angles as it stems only from the helicity-flip Compton amplitude. However, double-logarithmic enhancement takes place at large scattering angles from a region of quasi-real Compton scattering [126], where both exchanged photons have low virtualities, while the mass of the excited hadronic state is the largest allowed by the incoming beam energy.
Below the threshold of two-pion production, the absorptive part of proton Compton amplitude that governs the beam asymmetry can be described, as required by unitarity, in terms of quadratic combinations of single-pion production amplitudes. This program was realized in Ref. [126] where, using MAID amplitudes for electroproduction of pions, the authors reached good agreement for the beam asymmetries measured at MAMI for a broad range of scattering angles.

Medium momentum transfer
The dispersion relation approach is an appropriate tool to improve our knowledge of inelastic contributions to the most uncertain TPE corrections for arbitrary scattering angles in elastic ep scattering at GeV energies and below. The total TPE correction is given as a sum over all possible intermediate states, proton, pion-nucleon, resonances, etc. As was mentioned in Sec. 3.2.2, the elastic (proton) intermediate state was included within unsubtracted dispersion relations in Ref. [86] and within subtracted dispersion relations in Refs. [88,89]. The first inelastic pion-nucleon intermediate state was evaluated from data-driven MAID parameterization for the pion electroproduction amplitudes [130,131] at low momentum transfers (Q 2 ≲ 0.064 GeV 2 ) in Ref. [132]. This straightforward, but numerically involved calculation, consistently accounts for the contribution of all resonances and non-resonant background. Going to higher momentum transfers requires analytical continuation of ep scattering amplitudes to the unphysical region of kinematics. Such a novel method of analytical continuation for multiparticle intermediate states was developed and tested in Ref. [133]. It allowed authors to evaluate pion-nucleon contributions to TPE corrections in elastic ep scattering at more than an order of magnitude larger momentum transfers (Q 2 ≲ 1 GeV 2 ) compared to their previous work [132].
For higher Q 2 , the dispersive approach for resonant intermediate states was developed in Ref. [134], and implemented for spin-1/2 and spin-3/2 resonances with mass below 1.8 GeV in Ref. [135]. They included a Breit-Wigner shape with a finite width for each resonance, and for the resonance electrocouplings at the hadronic vertices they used helicity amplitudes extracted from the analysis of CLAS electroproduction data [136].
3 GeV 2 the net excited state resonance contributions are small, and the total correction is dominated by the nucleon elastic intermediate state. The net effect of the higher mass resonances is to increase the magnitude of the TPE correction at Q 2 ≳ 3 GeV 2 , due primarily to the growth of the (negative) odd-parity N (1520) 3/2 − and N (1535) 1/2 − resonances, which overcompensates the (positive) contributions from the ∆(1232)3/2 + . At the highest Q 2 = 5 GeV 2 value, the total TPE correction δ tot reaches ≈ 6% − 7% at low ε. An estimate of the theoretical uncertainties on the TPE contributions can be made by propagating the uncertainties on the fitted values of the transition electrocouplings [136], which are dominated by the ∆(1232) 3/2 + and N (1520) 3/2 − intermediate states.
At low Q 2 , Q 2 ≲ 0.5 GeV 2 , the uncertainties are insignificant, but become more visible at higher Q 2 values, as illustrated by the shaded bands in Fig. 12, for Q 2 = 1 − 5 GeV 2 .

Large momentum transfer
For TPE calculations at large momentum transfers, quark degrees of freedom have to be included. In an approach based on generalized parton distributions and a handbag mechanism (with two photons coupling to the same quark) [41,42], TPE effects for the unpolarized cross section and for polarization asymmetries were calculated. Perturbative QCD (pQCD) approaches were considered in Refs. [38,137]. In particular, in Ref. [38], special consideration was given to the coupling of two photons to different quarks, noticing that a large momentum may be transferred with one less gluon compared to the case of pQCD description of elastic proton FFs.

Effective field theory calculations
For elastic ℓp scattering at lepton energies well below the proton mass, we can use low-energy effective field theories (EFTs). At such energies, we do not need the full knowledge of the non-perturbative properties of the proton. For OPE, we only need a few effective couplings, e.g., the proton charge, magnetic moment, and charge radius. Similar simplifications apply to TPE. In the following, we discuss two such effective theories: heavy baryon chiral perturbation theory (HBChPT) and QED-NRQED, and review the current status of their predictions for TPE. Before a detailed discussion, it is worth mentioning that radiative corrections in the elastic electron-proton scattering when Q 2 >> m 2 e formulated recently in the Soft-Collinear Effective Theory (SCET) approach within QED [138], which allows us also for a systematic resummation of large logarithms. Moreover, the two-photon exchange corrections in the hard momentum region were also evaluated in the SCET framework within QCD in Ref. [139].

Heavy baryon chiral perturbation theory
Within this topical collection, an exact analytical evaluation of the TPE contribution to elastic ℓp scattering in low-energy EFT, namely, HBChPT at NLO, is published [140]. The LO HBChPT Lagrangian contains one derivative only. The NLO Lagrangian will have two derivatives and terms of order 1/M , where M is the proton mass. In this framework, the proton being the heavy degree of freedom is treated nonrelativistically. The leptons, the electron and muon, along with the photons constituting the light degrees of freedom, are treated in standard covariant QED. Throughout, the masses of the relativistic leptons are kept non-zero. Previously, the NLO HBChPT prediction had been calculated using the soft-photon approximation (SPA) [141,142]. As before, Ref. [140] assumes that at low energies the dominant photon loop contributions arise only from the elastic proton intermediate state.
The inelastic proton intermediate states are beyond the intended accuracy of the evaluation. Notably, since the proton FFs, which parametrize the hadronic structure effects, enter in HBChPT at NNLO, the proton is point-like at NLO. All hard-and softphoton exchanges are now included in the two-photon loop corrections which are evaluated exactly. The methodology relies on the successive use of partial fractions and integration by parts facilitating the decomposition of the rather intricate TPE four-point loop functions into a system of two-and three-point scalar master integrals, which are straightforward to evaluate analytically.
To LO accuracy, our low-energy TPE result approximates the well-known McKinley-Feshbach two-photon correction in potential scattering theory, cf. Eq. 70. At NLO the finite analytical real parts of the TPE loop contributions to the elastic differential cross section are presented. The exact analytical evaluation of the NLO bremsstrahlung process is not included. This means that the systematic cancellations of the IR-singular terms that arise from the TPE loops are not discussed. While the details of the TPE evaluation are going to be presented in Ref. [140], a rigorous NLO evaluation of the bremsstrahlung subtracted IR-free TPE contribution to the radiative corrections shall be a subject matter of a later publication.

QED-NRQED
An EFT related to HBChPT, is QED-NRQED, suggested in [143][144][145][146] and further developed in [107,147,148]. In this EFT the proton is treated nonrelativistically, using NRQED, and the leptons and photons are treated using QED. Unlike HBChPT, the EFT does not include other hadronic degrees of freedom, such as pions. Such effects are encoded in the QED-NRQED coupling constants. Thus QED-NRQED is most useful for elastic ℓp scattering.
Up to dimension six, the independent NRQED Wilson coefficients are c F and c D corresponding to the proton magnetic moment and charge radius, respectively. There are also two dimension-six contact interaction terms: a spin-independent term and a spin-dependent term [146], where ψ (ℓ) is the proton (lepton) field and M is the proton mass, the cutoff scale of the EFT. In Ref. [147], it was shown that QED-NRQED ℓp scattering at O(Zα) and power 1/M 2 reproduces the known Rosenbluth scattering formula expanded to 1/M 2 . It requires just the Dirac Lagrangian and the NRQED Lagrangian up to 1/M 2 , implying that b 1 and b 2 are zero at O(Zα). QED-NRQED ℓp scattering at O(Z 2 α 2 ) and leading power in 1/M reproduces the O(Z 2 α 2 ) terms in the scattering of a lepton off a static 1/r potential [149]. Here Z denotes the proton charge in units of |e| and is used to denote the appropriate class of radiative corrections.
Since the Wilson coefficients b 1 and b 2 are zero at O(Zα), they are sensitive to TPE at scales above the proton mass. In Ref. [148], b 1 and b 2 were calculated at O(Z 2 α 2 ). They were extracted by calculating the ℓp → ℓp off-shell forward scattering amplitude at O(Z 2 α 2 ) and power 1/M 2 in the effective and full theory in both Feynman and Coulomb gauges. Two cases were considered: a toy example of a nonrelativistic point particle (p.p.), and the real proton described by a hadronic tensor, namely, the forward VVCS off the proton discussed in Sec. 3.2.3.
where Λ is the UV cutoff of QED-NRQED, and Q ℓ is the lepton charge in units of |e|. Surprisingly, b p.p.
. For the case of the real proton, Ref. [148] derived implicit expressions for b 1 and b 2 in terms of the components of the hadronic tensor. Considering only the contribution to the Wilson coefficients of F 1 (0), F 2 (0) and M 2 F ′ 1 (0), related to the proton charge, magnetic moment and charge radius respectively, Ref. [148] finds: The ellipsis denotes non F 1 (0), F 2 (0), M 2 F ′ 1 (0) terms. Surprisingly, again there is no contribution to b 1 . It does not follow from a symmetry of the EFT and it might be a one-loop "accident".
This implies that low-energy elastic ℓp scattering is much less sensitive to spin-independent TPE effects above the proton mass scale compared to spindependent ones. On the other hand, the proton charge radius extraction will be more robust.
With the calculation of the QED-NRQED dimension six couplings: c F , c D , b 1 , b 2 , at hand, the next logical step is the calculation of the differential cross section for elastic ℓp scattering. Such a calculation was not performed in the literature yet. Once performed, it would be interesting to compare its result to HBChPT.

Empirical extractions of TPE
Besides theoretical predictions, it is possible to extract TPE amplitudes, as well as the TPE contributions to the unpolarized cross section δ 2γ , from empirical cross section and polarization transfer measurements.

Phenomenological fits of δ 2γ
The size of the TPE contribution becomes more significant at larger Q 2 . A phenomenological fit of the Q 2 -dependent part, as an addition to the Feshbach correction, Eq. 70, is [150] where a = 0.069, b = 0.394 GeV −2 . This fit includes the world data set on unpolarized cross section measurements (updated to a common standard on radiative corrections) and measurements of the FF ratio using polarization and assumes that the whole discrepancy is driven by hard TPE. Similarly, Schmidt in Ref. [151] uses several existing parametrizations for FFs from Rosenbluth-type experiments and polarization measurements of the ratio to determine the hard TPE contribution, under the assumption that TPE preserves the linearity of the reduced cross section in ε.

Updated extraction of TPE amplitudes
A measurement of the three observables for elastic ep scattering, σ R , P l , and P t , as a function of ε for a fixed value of Q 2 allows one to extract in a modelindependent fashion the three TPE amplitudes, which conserve the lepton helicity, as can be seen from Eqs. 58, 45 and 50. Such an analysis has been performed in [38] at Q 2 = 2.5 GeV 2 , based on available cross section data as well as the ε-dependence of P t /P l and P l /P Born l as measured by the JLab/Hall C experiment [152]. The combination of both unpolarized and polarization experiments at a same Q 2 value provided the necessary three observables to extract the ε-dependence of the three TPE amplitudes Y M , Y E , and Y 3 , introduced in Eq. 49. A concise update of this analysis in view of more recent data is given here.
The data for P t /P l of the dedicated JLab/Hall C GEp-2γ experiment [66,152], as shown in the upper panel of Fig. 13, do not see any systematic TPE effect within their error bars of order 1%. One can therefore effectively fit the observable on the lhs of Eq. 59 assuming an ε-independent part, which equals its OPE limit, and extract the Q 2 -dependence of R EM , cf. Eq. 46, from this observable, using the parameterization given in Ref. [66] (Global fit II). For the longitudinal polarization observable P l /P Born l , shown in the lower panel of Fig. 13, the deviation from unity is 6.2 times the statistical uncertainty and 2.2 times the total uncertainty [66]. The TPE correction effect has been parameterized as [38] P l /P Born JLab/Hall C data are from the GEp-2γ experiment [152], using the updated analysis of [66]. The red curves are the fit to the data. Left: ε-independent fit according to Ref. [66]; right: updated fit from Ref. [38] according to Eq. (81). Note that the values for the polarization observables are shown without radiative corrections. The bulk of these corrections (virtual corrections on the lepton side and soft-photon emission corrections), which factorize in terms of the Born cross section, drop out of the asymmetries, and the hard TPE is therefore expected to be the leading correction. It has been checked in Refs. [66,152] using the MASCARAD program [153] that the standard radiative corrections yield a multiplicative correction around or less than 0.1% on the asymmetries.
The comparison of this fit function to the data is shown in the lower panel of Fig. 13 for a l = 0.09. Recently, an updated global analysis for σ R , including new high-Q 2 data from JLab, has been given in Ref. [63]. In Ref. [63], the data for the reduced cross section have been fitted by a linear in ε behavior with µ p the proton magnetic moment, and G D (Q 2 ) ≡ 1/(1 + Q 2 /0.71) 2 the standard dipole. As the TPE amplitudes vanish for ε → 1, the proton FF G M can be extracted from Eqs. 57 and 82 as The parameterization for σ R chosen in [63] does however not have the correct ε → 1 limit, in which the TPE amplitudes vanish. At a fixed value of Q 2 , the limit ε → 1 corresponds with the Regge (high-energy) limit, as Q/(2ϵ 1 ) ≃ (1 − ε) 1/2 / √ 2 for ε → 1, with ϵ 1 the electron lab energy. In an expansion around ε → 1, the leading TPE correction can be model independently expressed as [121] δ 2γ (ε → 1, In order to combine the empirical observation of a linearity of σ R over most of the ε region with the correct ε → 1 limit, the fit function of Eq. 82 for σ R can be improved as which does not change the fit value for small ε and has the correct limit for ε → 1: Using the fit function of Eq. 85, the TPE correction δ 2γ follows as In Fig. 14 (upper panel), we compare, for Q 2 = 2.5 GeV 2 , the fit functions for σ R of Ref. [63] (red dotted curve) with the modified fit function of Eq. 85 (green solid curve), as well as σ R in the OPE approximation, setting δ 2γ = 0, (blue dashed curve). The resulting εdependence of the TPE correction δ 2γ for Q 2 = 2.5 GeV 2 , based on the linear fit of Ref. [63], using the improved fit function of Eq. 85, using the spline and Padé fits of Ref. [150], are compared in Fig. 14   Based on the measured ε-dependence of the three elastic e − p observables (σ R , P l , and P t at Q 2 = 2.5 GeV 2 as shown in Figs. 13 and 14), an updated empirical extraction of the three TPE amplitudes, Y M , Y E , and Y 3 , at this Q 2 value was obtained, see Fig. 15.
The updated empirical extraction of the TPE correction δ 2γ for general values of Q 2 and ϵ, is compared to data from the OLYMPUS, VEPP-3 and CLAS experiments in Figs. 16, 17 and 18 (dot-dashed pink curves), respectively. Note that the empirical extraction displays a small oscillation in the low-Q region, cf. ε = 0.45, 0.88 panels in Fig. 18. This, however, is not a true feature but a remnant that propagated from the underlying fit of R EM (Q 2 ) that was used [66], and is enhanced due to two opposite sign contributions in the (1 − ε) term in Eq. 88.

Single-spin asymmetries
As stated above in Sec. 2.2.5, parity-conserving SSA for elastic ep scattering are zero in the first Born approximation. Therefore measurements of either target or beam SSA may provide useful information on TPE amplitudes. While beam asymmetries arising from TPE are of the order 10 −5 for electrons [30,123,126] and 10 −3 for muons [154], the target SSA are expected at percent level [42]. The beam asymmetries were measured with high accuracy at MIT [155], MAMI [156] and Jefferson Lab [157][158][159][160] using experimental setups designed to study parity-violating electron-helicity asymmetries. Measurements of target SSA were performed at Jefferson Lab on a polarized 3 He target [161]; results for Q 2 > 1 GeV 2 were in agreement with GPD-based calculations [42] for a neutron. The measurements of SSA in elastic electron scattering provided unambiguous evidence of TPE effects, while probing contributions arising from an absorptive part of the hadronic Compton amplitude entering the TPE mechanism.

Cross sections and charge asymmetries
The TPE correction, defined in Eq. 56, can be accessed experimentally by measuring the ratio of cross sections from elastic ℓ + p vs. ℓ − p scattering. Corrections that depend on an odd power of the lepton charge do not cancel in this ratio. This includes the leading TPE correction, δ 2γ . It also includes, however, the interference of lepton and proton bremsstrahlung radiation, δ b , which is comparable to δ 2γ . Their sum gives us the "odd" part of the radiative corrections, δ odd = δ b + δ 2γ . By contrast, the "even" part of the corrections, δ even , include logarithmically enchanced terms, log(Q 2 /m ℓ ), coming from lepton bremsstrahlung, vacuum polarization and vertex corrections. As a result, the even part is relatively large, implying |δ odd | ≪ |δ even |. The above mentioned cross section ratio can be written in the following way [45]: The TPE correction can be extracted from which is the measured ratio R exp.
± corrected for δ b and δ even .
Experimental possibilities are limited, as such a measurements requires both electron and positron beams of high quality and at relevant beam energies of a few GeV. In recent years, three experiments have published results, cf. Figs. 16-18: an experiment [162] at the VEPP-3 storage ring measured the cross section ratio with beam energies of 1.0 and 1.6 GeV with a non-magnetic spectrometer covering angles between 15 • and 105 • . Since the experiment lacks a relative luminosity determination of sufficient precision, the data is published relative to the most forward point set to a ratio of one. However, this point is sufficiently far away from ε = 1 that a proper comparison with curves requires a shift of the data to match the predictions at this normalization point. CLAS at Jefferson Lab [163,164] generated a wide-energyspread electron and positron beam by converting the CEBAF electron beam to a photon beam via a converter target, subsequent pair-production from that photon beam, and by cleaning and recombining the electron/positron beam via a chicane magnet system.
Using the detector setup of the CLAS collaboration in Hall B, the so produced wide-energy beam-proton scattering events are detected. By the curvature of the scattered lepton in the magnetic field, the charge of the particle can be determined. Because of the wideenergy spread of the beam, the resulting data is sorted in large bins in ε and Q 2 . The collaboration published data with several (non-independent) bin selections. The OLYMPUS experiment [165] at the DORIS ring at DESY, Hamburg, measured using a monoenergetic beam with an energy of around 2 GeV using an open-cell internal hydrogen target. Scattered leptons and protons were detected with a spectrometer upgraded from the former BLAST detector. OLYM-PUS reached the highest Q 2 of the three experiments of slightly more than 2 GeV 2 .

Comparison theory vs. experiment
The OLYMPUS data (Fig. 16) at smaller Q 2 / larger ε are below unity. This is only replicated by the (mostly) phenomenological models, while theoretical predictions produce a ratio below 1 only much closer to ε = 1. Overall, the phenomenological predictions compare well to the data, but both predict a somewhat larger effect at larger Q 2 / smaller ε, albeit within the uncertainty of the data. Note that, calculations prior to and motivating the experimental proposal [38] predicted a large and visible deviation of the ratio from unity already around Q 2 ≃ 2 GeV 2 .  The VEPP-3 (Fig. 17) and CLAS data (Fig. 18) also prefer the phenomenological extractions; however, the CLAS data are mostly compatible with Refs. [135] and [133] as well. The CLAS ε = 0.88 data set and the VEPP-3 E beam = 1.594 GeV data set show a ratio below 1 for small Q 2 / large ε, similar to the behavior seen in the OLYMPUS result.  A global comparison of the data from the three experiments and the calculation from [167] was done in Ref. [168]. The largest contribution to the charge asymmetry is due to the soft photon emission. A specific procedure was suggested to limit the effect of the different corrections applied to the data. The results do not favor an enhanced contribution of TPE in the considered kinematical range.

Outlook
The OLYMPUS data gives a hint that at higher Q 2 , TPE might not explain all of the discrepancy. Theoretical descriptions, which do not show a particular good agreement in the measured Q 2 range, are not at all tested at these larger Q 2 . Precision measurements of the FFs at these large Q 2 however depend on an accurate description of TPE. This makes experimental determinations of TPE at these kinematics highly desirable. Unfortunately, experimental activities are limited by the availability of suitable beams. The TPEX experiment [169] seeks to measure TPE at DESY using beams of 2 and 3 GeV and up to a Q 2 of 4.6 (GeV/c) 2 . This measurement would require a new extracted-beam line at the DESY test beam facility. Current plans for future construction at DESY make it unlikely that this experiment can be realized.
Jefferson Lab is planning to construct a positron source for CEBAF. Current timelines put the availability of beam beyond 2030, if funding can be secured. Such a capability would allow for a comprehensive TPE program, with possible measurements in Halls A, B and C [170][171][172][173].
At smaller Q 2 , MUSE [174,175] and AMBER [176] will take data for µp scattering with both beam charges. From this data, TPE can be determined at small Q 2 close to the validity of the Feshbach limit, and in the case of AMBER, at ε extremely close to 1. Theoretical predictions for the inelastic and total TPE in the kinematic region envisaged by the MUSE experiment are shown in Fig. 19.    [122] for the three different beam momenta envisaged by the MUSE experiment [177]. While the elastic TPE dominates, the main uncertainty stems from the inelastic part. The relative uncertainty for the total TPE is thus small and not displayed here.

Implementations and Generators for Lepton-Proton Scattering
Eventually the radiative corrections discussed above need to be made available to the experiments. For less precise experiments it is often sufficient to do this in terms of fiducial (differential) cross sections that include a (often simplistic) model of the kinematic and geometric acceptances of the detector. Doing this is very simple from a theoretical viewpoint as fiducial observables can easily be calculate by integrating over phase space with a measurement function S as defined, e.g., in Eq. 6. Even for very complicated measurement functions, this is a reasonably simple procedure. However, for precision experiments, such as the current and next-generation ℓp scattering experiments, fiducial cross section are insufficient; instead, corrections need to be included in the entire analysis pipeline in the form of an event generator. Using only elastic LO events in the experimental simulation would introduce a large systematic uncertainty that would be impossible to reduce with only fiducial calculations. Hence, the detector should be modelled using at least an NLO event generator. To account for these effects properly, however, one really needs an NNLO generator that is ideally matched to some kind of parton shower to further capture leading logarithms beyond NNLO. In other words, the experimental simulation needs to use the best theoretical calculation available. Similarly, it is important that the applied radiative corrections are well documented, so that future re-analyses can update the corrections with more advanced calculations.

Overview of experimental requirements
Experiments vary greatly in their requirements for kinematic range, supported event topology, and precision. Typically, the precision goals are informed by the size of other uncertainties, mainly the statistical precision, but also other unrelated systematic uncertainties, for example from geometry. But this can be misleading. While statistical uncertainty and often also other systematic sources have no or small correlation with kinematics, i.e. can be assumed to statistically average over many measurements, this is not true for radiative corrections, which have the potential to be highly correlated not only between data points of the same experiment, but even between multiple experiments if they use the same or similar formalism. Further, experiments might measure different event structures. In the case of ℓp scattering, for example, the experiment could detect outgoing lepton, proton, or both, affecting the effectively seen corrections. Additionally, experiments can vary widely in the energy resolution, possibly requiring to describe the radiative tails accurately over rather large photon momenta ranges. Table 1 gives an overview over contemporary experiments [63,150,[174][175][176][178][179][180][181][182][183][184][185][186]. We list beam energies, angle and Q 2 ranges, as well as detected particles. We also calculate an estimate of the effect of the proton radii on G E and G M as well as the fractional contribution of G E to the cross section, as an estimate for requirements on systematic uncertainties. Of the planned experiments, only AMBER will measure the scattered lepton in coincidence with the recoiling proton. AMBER, PRad and PRAD-II both have forward kinematics with ε very close to 1, while the other experiments tend to measure at larger angle ranges, ε < 1. AMBER and MUSE will measure with leptons of both charges. The lower energy experiments require proper treatment of the lepton mass, and in the case of ULQ2, MUSE and likely MESA, a treatment without ultra-relativistic approximations.
Note that in most experiments at small Q 2 , the absolute normalization of radiative corrections is less of an issue-it is generally left floating and determined by an extrapolating fit to Q 2 = 0. Typically experimental normalization uncertainties are significantly larger than radiative correction uncertainties, so that such a fit is unavoidable. On the other hand, experiments aimed to extract the radii typically measure at very small Q 2 to reduce systematical errors from the fit extrapolation to Q 2 = 0. Consequently, radiative corrections need to control relative, Q 2dependent, errors to a very high degree: The proton radius puzzle is a 4% difference in extracted radii, for recent reviews see Refs. [112][113][114]187]. A meaningful measurement would extract the form factor slope to better than 0.5-1% or so. A radius measurement at small Q 2 where the radius/slope effect is, say, 5%, would then need to control Q 2 -dependent systematic uncertainties over this Q 2 range to better than 0.05%.
For Rosenbluth separations at somewhat extreme kinematics, one of the form factor contributions to the measured cross section can be highly suppressed compared to the other one: G M is suppressed at small Q 2 , while G E is at large Q 2 . The uncertainty of radiative corrections on the cross section-in this case mainly the ε-dependence at constant Q 2 -then gets levered up by potentially large factors in the determination of the suppressed form factor.
A further consideration is the requirement for unweighted events, or variance-reduced weights. This strongly depends on the overall speed of the remaining MC chain. For example, the Mainz A1 MC simulation code does not employ particle path tracing through magnetic fields and can simulate thousands if not millions of events per second, enabling to use fully weighted events, possibly even those with negative weights. In contrast, OLYMPUS' simulation chain included path tracing, full digitization and tracking for every simulated event, resulting in speeds of few events per second. This requires either unweighted events or at least measures to reduce the variance of the event weights.

Overview of theoretical tools
For now we will focus on tailor-made event generators for ℓp experiments even though much effort has also been devoted to the development of generators at high energy (see [188] for a recent and comprehensive review). However, we will revisit a potential synergy in Sec. 4.2.3. For ℓp scattering, the most notable tools are ELRADGEN [189,190], ESEPP [5], and Simul++ [150,180,191]. Another tool in this context is McMule [12,192] (as part of this topical collection, a simple TPE model is compare to the effect of NNLO corrections using McMule [20]) which currently cannot generate events. However, work is underway to address this in a way that is convenient for experiments. Looking slightly beyond ℓp scattering, we have the code by Epstein and Milner (EM) [193] for Bhabha and Møller scattering; MERADGEN [194] for Møller scattering; MES-MER [11,21,195,196] for muon-electron scattering; and BabaYaga for Bhabha scattering and photon pair production [197][198][199][200][201]. While these generators are not directly applicable for ℓp scattering experiments, they could either be adapted into (trivially in the case of MESMER) or help inspire new ℓp generators. Hence, we will still include them in our comparison as shown in Tab. 2.
When we compare MC codes, we consider the order in perturbation theory they implement, whether they use resummation. Another important point is whether they use any further approximations such as as an ultra-relativistic approximation (which treats the electron as nearly massless E ≫ m) or the soft-photon approximation (which assumes E, √ s ≫ E γ , also known as eikonal or peaking approximation). We stress that the order in perturbation theory is always relative to the listed process, usually ep → ep. This means that an NLO calculation of ep → ep includes an LO calculation of ep → epγ as a subset (cf. Sec. 2.1). If one wishes higher accuracy, i.e. NLO, for the latter process, more complicated calculations are required which, if properly implemented, would allow the implementation of an NNLO code. As discussed in Sec. 2.1 there are two commonly used methods: slicing (cf. Eq. 12) and subtraction (cf. Eq. 13). Slicing tends to be more commonly used in NLO QED calculation as it is simpler to implement and use to generate events. While most MC generators discussed here are fixed order, either at NLO or NNLO, some also include a resummation. We will discuss the different strategies relevant here in Sec. 4.2.1. Since the ability to produce unweighted events is very important from an experimental point of view, we will discuss different strategies in Sec. 4.2.2. Finally, while Table 2 focuses on purpose-built tools for low-energy experiments, one should not forget about the generalist tools that were developed for the LHC and other highenergy experiments. In Sec. 4.2.3 we will discuss some potential synergies.

Parton shower and resummation
When calculating higher-order corrections, be it at NLO or NNLO, one often finds that the corrections are enhanced by the presence of large logarithms L between disparate scales such as ln(m 2 e /s), ln(m 2 e /m 2 p ), or ln(Q 2 /m 2 p ). These logarithms somewhat spoil the perturbative convergence since we usually get at least one such logarithm per order in α. This means that NLO, and sometimes even NNLO, is insufficient. One way out of this problem is just to compute even higher orders in α. However, this quickly becomes infeasible, with partial N 3 LO corrections just about feasible. Luckily, there is rarely a need to compute the full corrections since typically only the logarithmically enhanced matter. These terms follow a very predictable structure and can often be calculated to all orders in α. This process is called resummation and it introduces its own counting. Accounting for all terms of the form α n L n is referred to as leading logarithm (LL), the second "tower" α n L n−1 is called next-to-leading logarithm (NLL) and so on.
When combining a calculation at a given level of resummation, say LL, with a fixed order calculation at e.g. NLO, one must take care to avoid doublecounting. Both calculations include the same O(αL) term, albeit in different forms. Ensuring that this does not spoil the result requires the two calculations to be matched at the given order.
In QED we have two main sources of large logarithms: those related to the smallness of the electron mass (usually called collinear logarithms) and those related to restricting real radiation, either through explicit cuts on the phase space integration or by considering the endpoints of distributions (usually called soft logarithms).
Logarithms can either be resummed using a numerical algorithm such as a parton shower (see below) or analytically. A parton shower resummation contains the full LL but may contain partial results of the NLL as well, depending on method and implementation. The development of an exact NLL parton showers is an area of very active research at the high-energy frontier. Analytic resummation can be performed well beyond LL with state-of-the-art calculations reaching next-to-next-to-next-to-leading logarithm (N 3 LL). However, these calculations are observable dependent and cannot easily be used to generate events.

YFS exponentiation
The simplest way to resum soft logarithms is Yennie-Frautschi-Suura (YFS) exponentiation, a procedure that is based on the universality of soft singularities in QED [8]. In its simplest form, this accounts for large logarithms near the endpoints of differential distributions e.g. by exponentiating the singular behaviour of  [24,203,204] at LL. However, the same strategy has recently been used to calculate soft logarithms up to next-to-next-to-leading logarithm (NNLL) in the electron spectrum of muon decay.

YFS shower
This form of YFS exponentiation can be used to obtain very precise results for single distributions. However, without modification, it cannot be used to generate events. By using the YFS approximation on both real and virtual corrections, we can construct a locally finite prescription to generate an event with an arbitrary number of real but soft photons [205] This shower reproduces the soft LL for any observable as the resummation happens for each phase space point. Further, matching a YFS shower to a fixedorder calculation is relatively straightforward since one can relatively easily incorporate the fixed-order calculation in the YFS shower [205].

DGLAP-based parton showers
An alternative resummation approach is based on the structure function approach. Given an electron with a fixed momentum p, the structure function D(x, Q 2 ) is the probability density of observing an electron with momentum xp and virtuality Q 2 . Convoluting a structure functions for each external electron with the fixed-order cross section then resums the collinear logarithms to LL accuracy. The structure functions themselves are governed by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation [206][207][208] where P ee (x) is the (regularised) kernel that is known up to two-loop accuracy. A parton shower is an algorithm of solving the DGLAP equation numerically using MC methods [209]. This allows us not only to obtain an exact numerically solution for the structure functions but also to generate approximate photon momenta. The latter can be improved by using the YFS theorem [198]. Work is still ongoing to reach NLL accuracy with a general parton shower. However, analytic solution of the DGLAP equation at NLL are already possible in QED [28].

Unweighting procedures
So far, we have only discussed the generation of weighted events. This means that we generate a set of random momenta {p µ i } and use them to calculate an event weight w i that is comprised of the matrix element (squared) and a phase space mapping that accounts for the way we have generate the momenta. We then use this event to perform some sort of histogramming, potentially involving a detector simulation. The event weight itself may vary across many orders of magnitude and even be negative. This increases the number of events that need to be generated to obtain a reasonable uncertainty on the final histogram since negative events will lead to (potentially severe) cancellations.
If the histogramming step is comparatively fast, one would normally stop here and just increase the number of events generated. However, a realistic detector simulation can be extremely complicated and require hours of CPU time which is expensive even compared to the matrix element evaluation that usually takes less than a second. Hence, our goal needs to be the reduction of the number of events that need to be passed to the detector simulation. The easiest way to do this, is to generate events with uniform weight even if this requires the evaluation of many more events on the theory side. Even in the absence of negative weights, this is a non-trivial task: given a probability distribution function (PDF) that we in general only know numerically, we want to generate variables that follow this distribution.
In simple cases, it is possible to generate the event directly without weight; the energy of an emitted photon can be sampled as where 0 < # < 1 represents a uniformly distributed random number and E min some cut-off parameter. This is possible because soft photons are distributed following the eikonal approximation. However, this fails as soon as one considers higher-order corrections as the distribution becomes much more complicated. A more complicated method is called hit-andmiss: we generate a set of weighted events with weights w i and find the maximum weight W = max |w i |. A given event is rejected with probability p = 1 − |w i |/W , otherwise it is accepted with weight sign(w i )×W . At this stage we have significantly fewer events than we started with but their weights are all ±W .
This strategy can be somewhat improved upon by sampling the original distribution using an adaptive MC sampler rather than naively using uniformly distributed random numbers. The simplest such sampler is VEGAS [210] which will sample regions with large weights more often than those of low weight, improving efficiency. It does this by splitting each axis of the integration into (by default) 50 subdivisions that are each sampled the same. The sizes of divisions are than adjusted to sample regions of large integrand more finely. Note that this relies on the behaviour of the integrand being reasonably well aligned with the axis of the integration. This has been improved upon by cells instead of a grid in FOAM [211]. More recently, machine learning techniques were investigated for even more efficient sampling methods, e.g. [212][213][214].
Much effort has been dedicated to improve unweighting procedures as hit-and-miss is very inefficient, especially in the presence of many events with negative weights w i < 0. This can for example be done by optimising the event generation, i.e. ensure that fewer event with negative weights are generated in the first place by modifying the parton shower matching procedure [215,216]. Alternatively, one can carefully remove the events with negative weight without affecting physical observables [217][218][219].
Recently, this led to the development of cell resampling [220] which works process independently and preserves all physical observables. The idea here is to utilise the experiment's finite resolution that is anyway required to IR finiteness. We now create a small sphere around every event with negative weight, ensuring that the sphere is smaller than the experimental resolution and cannot be probed. Assuming enough event were generated, the sum of all weights in this sphere is positive since experiments always measure a positive cross section. We now replace all events by their absolute value, taking care not to disturb the sum of weights.

Synergy with high energy
Many MC generators were developed for the LHC physics program. While some of these are purposedbuilt for certain high-energy processes, other are more general. Some of the most mature, general-purpose tools are Sherpa [221,222], Herwig [223,224], Pythia [225], and MadGraph aMC@NLO [226]. While these tools were designed for hadron colliders, they could be extended for low-energy elastic e-p scattering. In the following, we briefly summarise the four main codes and what they could offer without major modifications.
While Sherpa usually generates its own NLO matrix elements, it can be fed the required expressions in through a shared library. This means that once provided with all matrix elements (including the TPE) as some kind of callable routine, it would be able to generate events at NLO accuracy out of the box. Further, it provides a full YFS shower [227] that can resum soft logarithm; collinear contributions are taken at fixed order.
Pythia itself is not a program but rather a library that can be used for event generation. As such, it is trivial to link to an external matrix element provider and it could then be used for the generation of NLO events. While naturally focused on QCD showers and hadronisation, it does include a QED shower that includes a fully coherent multiple treatment of photon radiation beyond the YFS approximation. Further, the resummation of initial-state collinear log is possible using a lepton parton distribution function.
MadGraph aMC@NLO is mainly designed for the automatic generation of NLO matrix elements which are then used in a QCD parton shower. The automatic generation is less helpful for our purposes but could be side-stepped with some effort. The code further provides a QED shower for final states and, similar to Pythia, can use a lepton parton distribution function for the inclusive resummation of ISR. Weight reduction to unit weights w = ±1 is only available when combined with the parton shower to avoid the introduction of an IR cut-off.
Herwig is a Monte Carlo event generator with QED showering. In normal operations for LHC physics, amplitudes are included up to NLO. However, for more complicated processes Herwig is already set up to obtain amplitudes from external libraries. This feature could be used to provide the generator with matrix elements relevant for ℓ-p scattering.
Facilitated by the development of many different generators, the high-energy community was forced to develop standardised interfaces for generators and event formats. The most commonly used ones are LHEF [228] and HepMC3 [229]. These standards are used in conjunction with the Rivet framework [230]. Rivet supports the development of generators by encouraging experimentalists to share their analyses in a common framework that can be applied both on data and Monte Carlo samples. New generators (at both the high-energy and low-energy frontier) will probably interface to these standards to cut down on complexity.

Conclusion and outlook
The precision requirements of modern ℓp scattering experiments require unprecedented theoretical calculations. A critical part of this is the development of event generators that can be used in experimental simulations. These requirements mean that we need to go beyond leading order for the radiative process ℓp → ℓpγ. This is also roughly equivalent to going beyond NLO for ℓp → ℓp. Currently, no event generator with these capabilities exist even though some codes (notably McMule and MESMER) could be adapted.
Unfortunately, the increased precision comes with a massive increase of the complexity of the distributions that need to be sampled. Luckily this is a known problem; the LHC community has worked on issues that are almost identical for decades. While not all ideas are immediately transferable to our process, many are. Additionally to the calculation of oneloop matrix elements (which is fully automatised for all standard theories), this includes methods to deal with IR singularities, negative weights, and numerical resummation using a parton shower.
Using these tried and tested techniques will allow us to construct event generators with (and beyond) NNLO accuracy for ℓp → ℓp in the coming years.

Higher-order corrections in QED inelastic scattering
The Rosenbluth formula [1] as well as the Akhiezer-Rekalo formalism [35,36], respectively, for unpolarized and polarized elastic ep scattering hold only in the OPE formalism where two real FFs, functions of Q 2 only, completely describe the reaction. In the case of TPE, instead of two real functions of Q 2 one has to deal with three amplitudes, neglecting the electron mass. Three additional charge-odd amplitudes, that are proportional to the electron mass, are of the same order, generally of complex nature and also depending on two kinematical variables. The latter are relevant for observables such as transverse polarization. Model-independent statements hold also in the annihilation region [231,232]. Assuming time and parity invariance and crossing symmetry, the reactions e + e − ↔ pp are described by the same electromagnetic FFs. However, in the time-like region of transferred momenta FFs are of complex nature, even in the case of OPE. Due to the large transferred momenta that are involved above the physical threshold one would therefore expect that TPE effects are enhanced.
In the time-like region, the angular distribution of one of the produced particles allows one to extract the moduli of the FFs. In principle such measurement is simpler than a Rosenbluth separation, as it requires only one setting of the collider equipped with a 4π detector. The ε(∼ cot 2 θ)-linearity of the Rosenbluth plot translates by crossing symmetry into a cos 2 θ dependence of the annihilation cross section. Then, a forward-backward asymmetry arises in the presence of charge-odd contributions, as hard TPE. The sum (difference) of the cross section at corresponding angles would cancel (enhance) the TPE contribution, taking advantage of its charge-odd nature.
The study of both annihilation reactions e + e − ↔ pp, related by time reversal is very instructive in this respect, as radiative corrections are in general different. In e + e − colliders, huge effort has been done to achieve a per-mille precision, based on the lepton structure function method [233]. This is made necessary by the use of the "initial state radiation" (ISR) process [234], that allows to scan a large region of transferred momenta, even in a fixed energy collision.
Radiative correction for the reactionpp → e + e − , induced by antiproton beams have been particularly studied [235] and are being implemented in a dedicated MC event generator (cf. Sec. 4), in the frame of the time-like FF program planned at FAIR [236]. The interest of the antiproton beam is also related to the possibility to measure the muon-antimuon final state [237].
The search of charge-odd contributions to the annihilation cross section has been looked for in the BaBar data, see Ref. [238], with the conclusion that they vanish in the limit of 2%, which is the order of magnitude of the interference between initial and final state radiation emission. The presence of TPE in the future data from PANDA at FAIR has been simulated in Ref. [236] pointing out that a contribution of ≥ 5% would be detectable.
The QED radiative corrections to both deepinelastic scattering and e + e − annihilation at large virtualities Q 2 and cms energies s ≫ m 2 e are very important for precision measurements since they can be differentially rather large in parts of the phase space. Furthermore, they strongly depend on the different possible measurements of the kinematic variables of the respective processes. Unlike the case in massless QCD, the emerging logarithms L = ln(Q 2 /m 2 e ) or L = ln(s/m 2 e ) are observable phasespace logarithms and cannot be resummed away as in the strictly massless case. In particular, they can be measured as has been done in many experimental analyses in deep-inelastic scattering and e + e − annihilation. In addition, the fermion mass corrections cannot be dealt with as purely massless corrections from the beginning, cf. e.g. [239]. They are, furthermore, not just given by splitting functions and massless Wilson coefficients, but one has to compute massive operator matrix elements (OMEs) instead of the splitting functions. Two renormalization group equations rule the corrections in the limit Q 2 ≫ m 2 e or s ≫ m 2 e . 1 These quantities are ruled by two renormalization group equations, cf. [240][241][242], containing also the charged lepton masses. They read 2 Here µ denotes the factorization and renormalization scale, g the QED coupling constant, β the QED β-function, γ m (g) the mass anomalous dimension, γ ij the QED anomalous dimensions, Γ li are the massive OMEs andσ lk are the massless sub-system scattering coefficients. Finally, s ′ (or (Q 2 ) ′ ) denote the sub-system cms energy (or virtuality). The renormalization group equations rule all corrections down to the constant term which does not depend on the scale logarithm. Here the renormalization group equations have been written in Mellin space. Since the radiatively corrected (differential) scattering cross sections are observables, the µ-dependence needs to be canceled, which is possible by a consistent expansion in the coupling constant g. Usually the consideration of the radiators Γ li is therefore not sufficient, with the exception of the leading order terms O(α k L k ), with α the fine structure constant, since there only the Born sub-system cross section contributes. They are leg-dependent as the other corrections and the relevant hard scales rescale with the collinear emission variable z accordingly as will be outlined below.
The representation of the radiative corrections using Eqs. 93 and 94 allows to calculate all but the power corrections (m 2 e /s) k , k ∈ N, k ≥ 1. The latter terms require the full phase-space integrals to be carried out and are believed to yield small contributions only in the limit of large scales.
In Sec. 5.1 we will consider the QED radiative corrections for the deep-inelastic scattering process. Section 5.2 deals with the initial state QED radiative corrections to e + e − annihilation and Sec. 5.3 contains the conclusions.

Higher-order corrections to deep-inelastic scattering
The QED corrections to deeply-inelastic scattering have been calculated in the context of different experiments starting with those at SLAC [3], CERN, with BCDMS and EMC [243], and with much more work for HERA, cf. Refs. [244,245]. These papers cover the O(α) corrections and we will mostly discuss higher-order corrections in the following. Furthermore, most of the calculations refer to unpolarizedlepton unpolarized-nucleon scattering. Later on we will also discuss the QED radiative corrections to polarized nucleon scattering. Early leading logarithmic O(α) approaches (LLA) are shown in Refs. [3,246,247], see also Refs. [24,248,249]. The O(α) electroweak corrections, also for νN scattering, were calculated in Refs. [250][251][252][253][254][255][256][257][258][259][260][261][262] and applied to data. This has been important for the measurement of sin 2 θ W from deep-inelastic data.
Complete neutral and charged current O(α) corrections to ep scattering were calculated in Refs. [263][264][265], calculating the bremsstrahlung contributions using MC integration. The neutral and charged current O(α) corrections to deep-inelastic ep scattering were calculated analytically for the differential cross sections d 2 σ DIS /dx/dQ 2 in Refs. [266][267][268]. The O(α) leading log corrections to the neutral and charged current process in leptonic variables were computed in Ref. [269], correcting earlier work in Ref. [270]. The LLA results agreed very well with the complete results in the neutral current case. For only leptonic contributions these corrections were calculated in Refs. [271,272] also. The radiative corrections are quite different for different choices of variables to define the differential cross sections, cf. Ref. [273]. The LLA corrections in the case of jet (or hadronic) variables were calculated in Ref. [274]. In Tab. 3 we display the z rescaling for some of the kinematic variables.
Here J (x, y, z) denotes the respective Jacobian and x = Q 2 /(sy), y = p.q/p.l, with p and l the nucleon and lepton 4-momentum and q 2 the virtuality of the process with Q 2 = −q 2 . z 0 denotes the hard bremsstrahlung threshold. The rescaling relations for other sets of variables, such as Jaquet-Blondel variables, the double-angle method, the θy and eΣ methods, see Refs. [245,273].
Let us illustrate the radiator method in lowest order LLA. One obtains Here L e denotes the radiation logarithm L e = ln(Q 2 /m 2 e ) − 1 and we have accounted for a nonlogarithmic correction in addition. P (0) ee denotes the O(α) QED splitting function for the ee transition [275]. The QED splitting functions and anomalous dimensions can be obtained from those of QCD [276][277][278][279][280][281] to three-loop order in the unpolarized and polarized cases by setting T F = 1, C F = 1, C A = 0. The photon-photon splitting function needs more care. The consistent construction of the different order radiative corrections including both the massive OMEs and the massless sub-system scattering cross sections is described in Ref. [282] in detail in Mellin space. Parts of the radiators depend also on generalized harmonic sums [283,284] while the others are given by harmonic sums [285,286]. However, their Mellin inversion to z-space can be written in terms of harmonic polylogarithms [287] at modified argument. Standard technologies, cf. Ref. [288], also allow one to calculate the radiators in z-space.
We also would like to mention that the transition matrix elements from e − ↔ e + were given in Ref. [241] to two-loop order. They are related to the different non-singlet anomalous dimensions γ N S,± [277,280] starting at two-loop order.
Different radiative correction packages such as TERAD, HERACLES, HELIOS and HECTOR were designed [243-245, 289, 290]. In addition to the O(αL) corrections, which are large in certain parts of the phase space, also the O(α 2 L 2 ) corrections were calculated for leptonic variables in Ref. [291] and for the whole set of kinematic variables in Ref. [292]. Related radiators were obtained in Refs. [293,294] and later corrections in Ref. [241].
Likewise, complete O(α) calculations have been performed for Jacquet-Blondel variables in Ref. [295] and for mixed variables in Ref. [296].
The O(α 2 L) corrections for the latter case were computed in Ref. [297], extending the results of Refs. [292,296]. All sub-leading corrections of this type are process dependent since they also depend on lower order hard scattering contributions. Therefore, resummations of non-leading radiators, cf. e.g. Ref. [28] and references therein, do not lead to a complete description. Furthermore, one has to deal with massive OMEs, cf. e.g. Ref. [241].
The QED corrections to the final state hadronic jet are dealt with as inclusive, since they usually cannot be resolved as electromagnetic showers do also occur in the hadronization process. Therefore, this part of the corrections is summed up according to the Kinoshita-Lee-Nauenberg theorem [298,299].
The radiators of the LLA corrections to higher order O(α k L k ) are process and radiation leg independent. Furthermore, the non-singlet corrections are the same in the unpolarized and polarized case. The radiators for the O(α 2 L 2 ) corrections were calculated in Refs. [291,292]. The corresponding higherorder corrections up to O(α 5 L 5 ) were computed in Refs. [203,[300][301][302][303][304][305][306], where the singlet corrections were calculated in Ref. [305,306]. One may also resum small-x contributions of O((α ln 2 (x)) k ), which has been performed in Ref. [307].
In the case of leptonic variables already the O(α) corrections do not only contain initial and final charged lepton radiation, but also the contributions due to the Compton peak, which was first dealt with in Ref. [3] and later considered in Refs. [266,269,271,308].
There are also hadronic initial state radiative corrections. In early studies [263,264,[266][267][268] some conceptional assumptions were still made, which are incompatible with the massless parton model, being generally assumed for the QCD corrections and also needs to be applied to the QED corrections. The correct treatment has been given in Refs. [274,309]. Later numerical studies were made in Refs. [310,311]. Table 3: Scaling behavior of various sets of kinematic variables for leptonic initial and final state radiation, cf. [245,273]. The shifted variablesQ 2 ,ŷ, the threshold z 0 and the Jacobian J depend on the choice of the external kinematic variables.
Initial State Radiationŝ = zs KinematicsQ 2ŷẑ 0 J (x, y, z) Finally, the radiative corrections to polarizedlepton polarized-nucleon scattering were also calculated at O(α) in Ref. [312] and for the singlet corrections to O(α 5 L 5 ) in Ref. [305], see also Refs. [190,313]. Finally, we would like to give some numerical illustrations on the leptonic QED radiative corrections. In Fig. 20 we show the corrections in leptonic variables. It is shown that the LO corrections deviate by some % from the O(α 2 ) corrections, which are therefore necessary. The soft exponentiation is adding much less. The reduction of the radiative corrections in deep-inelastic ep scattering is also possible by measuring the kinematic variables by other methods. Here we would like to mention, in particular, the double angle method, cf. [292], where the corrections are much smaller and do also behave rather flat as functions of the Bjorken variables x and y. In Figure 21 we present numerical results in the case of mixed variables [292,296,297] in O(α), corrected by the terms of O(α 2 L 2 ) and O(α 2 L). The dominant behaviour is given by the O(αL) corrections. However, several % of the corrections are added by the higher-order terms included. These are required for precision measurements.

Higher-order QED corrections to
e + e − annihilation Unlike the case in deep-inelastic ep scattering, the complete two-loop QED corrections at large cms energy have been calculated for e + e − annihilation [239], which are of central importance for precision measurements of the Z 0 peak, Z 0 H 0 production and tt production at threshold, cf. [314] and references therein. 3 In earlier attempts to calculate these corrections or parts of it in Refs. [294,317] the massless limit m e → 0 has been taken too early, which led to incorrect results. Using for the inclusive cross sections the factorization property [318], one may also calculate it referring to the massless sub-processes [319][320][321] and the corresponding massive OMEs calculated in Ref. [241]. Since the results of Ref. [241] did not agree with those from Refs. [294,317], there was also the possibility that factorization does not hold in the massive case. However, it could be proven by the full massive phase-space calculation in Ref. [239] that this is not the case and the reasons for the differences could be given in Refs. [239,314,321]. With the factorization relation at hand one can now calculate the three most leading logarithms in any order of α. This has been done to O(α 6 L 5 ) in Ref. [282].
In Tab. 4 we show the improvement of the peak position and the size of the width of the Z 0 boson by including higher and higher order QED initial state corrections, until we reach the level of ±10 keV as the theoretical error. We consider both the cases of an s-independent as well of the s-dependent width Γ Z .
The QED initial state corrections are also very important for the measurement of inclusive tt production in the threshold region, which will be used for a precision measurement of the top quark mass in the future. Evidently higher than the O(α 2 ) corrections are required to obtain a sufficient description, cf. Fig. 22. For a future precision measurement of the fine structure constant α at high energies [324] one needs also higher-order corrections to the forward-backward asymmetry. The corresponding LO corrections to O(α 6 L 6 ) have been calculated in Ref. [325] extending earlier results at O(α 2 L 2 ) [271]. The O(α) corrections were computed in Ref. [326].
The LLA initial state radiative corrected forwardbackward asymmetry is given by with σ and σ T (s) the radiatively corrected total cross section. The radiator is given bỹ where the parameter z 0 plays the role of an energy cut and z = s ′ /s. Here H LLA e is angular independent, and denotes the leading logarithmic contributions of the usual radiators, while H LLA FB , which is angular dependent and obtained by the following integral In z-space the radiator depends on cyclotomic harmonic polylogarithms [327].  In Tab. 5 we illustrate the improvement of the forward-backward asymmetry by including higher order leading-log initial state radiative corrections beyond the O(α 2 L 2 ) corrections [271]. The O(α 6 L 6 ) corrections stabilize the off-resonance forwardbackward asymmetries to six digits. Of course, also non-leading log corrections are needed to be calculated in the future.

Conclusions
We have summarized the status of the higher-order QED corrections to deep-inelastic ep scattering and for the QED initial state corrections to e + e − annihilation. In both cases the emphasis is on the nonpower corrections, i.e. on the logarithmic contributions down to the constant terms of O(α k L 0 ). These contributions can be computed using renormalization group techniques, which involve massive OMEs on one side and massless sub-system scattering processes on the other side. It is important in using this approach to cancel the dependence on the factorization scale by matching the different contributions. This is best done in Mellin-N space in a fully analytic way. The radiators are then given by Mellin convolutions of the respective pieces defining the OMEs and the sub-system scattering cross sections. The result can be Mellin-inverted analytically providing the correct radiators. The requested interplay between both types of contributions does not allow to study QED evolution only.
In the case of future experimental studies of e + e − annihilation it is planned to measure the width and the mass of the Z 0 -boson down to the 20 keV range. This requires initial state corrections of up to 6th order QED corrections by keeping three logarithmic orders.

Light-Hadron Decays and Reactions
The incorporation of radiative corrections is also essential for processes in a kinematical regime that is different from the previous discussions. Recently, it has been gradually recognized that in order to achieve desired precision in experiment, radiative corrections must be incorporated as a part of MC generators (see also Sec. 4 for more on MC concepts), and subsequently, the demand for such corrections has increased. Such an awakening resides in the fact that, typically, the QED NLO radiative corrections compete in size with the hadronic effects. And thus, once hadronic parameters are to be extracted from data, ignoring radiative corrections would lead to meaningless results, the examples of which we will see below.
Historically, the approach to radiative corrections varied. In some cases, they would be left out entirely. If an effort was done to include at least some of the relevant corrections, others that were not negligible were ignored. In some cases, approximate approaches (soft-photon limit or leading logarithms) were used, which often led to unreliable or misleading results. All the above-mentioned cases introduce artificial discrepancies between theory and experiment, as also discussed in Sec. 3.1. In other words, what is then considered to be "measured observables" or related "measured hadronic parameters" may include an unsubtracted QED part, which can turn out to be quite a significant contribution.
In the following, we will discuss the current status of QED radiative corrections in some sample decays of the lightest pseudoscalar mesons and baryons and their direct application in experiment.

Radiative corrections for π 0 decays
The NLO QED radiative corrections are now worked out for all the relevant neutral-pion decay channels: the π 0 Dalitz decay π 0 → e + e − γ, the π 0 double-Dalitz decay π 0 → 2(e + e − ), and the π 0 rare decay π 0 → e + e − . We will address these individually in the following subsections.
6.1.1 π 0 → e + e − Let us start with the last-mentioned process. Regarding the branching ratio, this channel is loop-and helicity-suppressed with respect to the radiative decay π 0 → γγ by eight orders of magnitude. It has thus been speculated that the π 0 → e + e − decay might be sensitive to possible effects of new physics. And it indeed seemed to be the case when the KTeV Collaboration published the results of their analysis on the precise measurement of the branching ratio in 2007 [328]. The direct subsequent comparison to the Standard Model (SM) prediction was then interpreted as a 3.3 σ discrepancy [329]. Many works followed trying to explain this difference both within (via modeling of the electromagnetic pion transition form factor) and outside the SM (by introducing various models including exotic particles). However, it is believed that the most important contribution in this regard was done by an explicit calculation of (two-loop) virtual radiative corrections [330] that brought the discrepancy down to the 2 σ level [331]: Until then, only leading-log approximation was available [332,333] that did not reproduce well accidental cancellations among various diagrams. New light can be shed on the branching ratio determination and the related remaining 2 σ discrepancy: A new measurement is under preparation by the NA62 Collaboration and the results should be known, conservatively estimating it, within a few months, i.e., by the end of 2023. In this analysis, all the known related radiative corrections [330,[334][335][336] are already included at the MC level.
Looking at the quantity that has been measured by KTeV [328], with x ≡ m 2 e + e − /M 2 π , it becomes immediately obvious why there is a need for precise knowledge of radiative corrections in order to compare this value with the predictions provided by theory, which can be written as B π 0 → e + e − ≈ 6.21 + 0.15 χ × 10 −8 , (102) with χ ≡ 2 χ (r) (770 MeV) − 5 2 . Theoretical predictions and models suggest (being rather conservative) χ (r) (770 MeV) ∼ 2-3 (see, e.g., Refs. [329,331,337,338]), so χ is expected to be within the interval (−1, 1). The result in Eq. 101 not only depends on the Dalitz-decay branching ratio serving here as the normalization channel, but the numerator allows for final states including soft photons with their energies limited by the requirement on the minimum of the normalized electron-positron invariant mass squared x. Eq. 101 thus needs to be further processed and it was done so in the KTeV analysis via a series of steps in order to provide a "no-radiation" value to be compared with theory. In other words, experimental results on the branching ratios of such processes will contain effects that need to be subtracted in some way, so they can be directly compared with theoretical predictions for branching ratios with no final-state radiation. Such a subtraction of radiative corrections can happen at different stages of the experimental analysis, but it seems that the earlier stage this happens, the better: Once the events accounting for photons in the final state are properly generated already at the MC level, comparing the resulting spectra with data helps to improve the quality of the analysis outcome.
We have already seen that in the KTeV rare-decay search described above, the Dalitz decay was used as the normalization channel. It is thus clear that it is essential to have the complete set of NLO QED corrections available also in this case. It then comes as a bonus that this process is also used for normalizing purposes in measurements in the kaon sector [336]. Furthermore, investigating the Dalitz decay itself leads to information about the singly off-shell (π 0 γγ * ) electromagnetic transition form factor (TFF) in the time-like region, i.e., direct access to the TFF slope without the need for model-dependent extrapolation from the space-like region, as used, e.g., in Ref. [339].
The radiative corrections for π 0 → e + e − γ were first estimated in terms of the total correction to the integral decay width [340], and only a decade later, the first corrections to the differential decay width appeared, which used the soft-photon approximation [341]. What can now be considered as a classical work on the Dalitz-decay radiative corrections came a year after [342]. The hard-photon corrections were included and the results were presented in terms of a table giving sample correction values throughout the Dalitz plot. One contribution was, however, missing: the two-photon-exchange contribution, or also known as the one-photon-irreducible (1γIR) contribution at one loop, which diagrammatically coincides with the bremsstrahlung contribution to the π 0 → e + e − [334]; see Fig. 23. It turned out, after a vivid + cross discussion in the literature [343][344][345], that this contribution was not negligible, and a direct subsequent calculation would show that omitting the 1γIR piece would introduce ∼ 15% discrepancy to the measured slope [335,346]. The most precise measurement of the slope in the time-like region to date was performed by the NA62 Collaboration with the results a π = 3.68(57)% [347]. It is worth mentioning that the slope of the radiative corrections amounts to circa −12.5%; see the low-x region of Fig. 24. This corresponds to the QED correc- tion of about −6% to the slope estimate in absolute numbers. This in turn means that without considering the QED corrections one would inevitably obtain a negative slope of similar size as the correct result a π stated above. In QED, in which hadronic effects of the π 0 → γ * γ * transition are included via the electromagnetic form factor, it turns out that one can precisely determine the Dalitz-decay branching ratio, directly following from the NLO calculation of the ratio R = B(π 0 → e + e − γ(γ))/B(π 0 → γγ) [336]. For this, it is sufficient to know the NLO radiative corrections and -from theoretical estimates confirmed by experiment -that the TFF slope is rather small. This then combines with the fact that the sensitivity to the only relevant (higher derivatives can be neglected since the pion mass lies well below the resonance region) hadronic parameter present (a π ) in such an equation (103) is suppressed because the region with high virtualities is (and the slope thus does not influence the final result much); this is represented by the (1 − x) 3 term above. Let us stress that such a precise determination is only possible since the complete set of NLO QED radiative corrections are now known. Believing the above determination, it becomes surprising that the KTeV Collaboration, soon after Ref. [336] was published, presented their results on the Dalitz-decay branching ratio [348] which is 3.6 σ away from the mentioned QED-based result. It is thus desirable to have another independent branchingratio measurement. There is ongoing activity in this direction in NA62.
A complete set of NLO radiative corrections in softphoton approximation for the neutral-pion double-Dalitz decay is computed in Ref. [349], in which the processes P →lℓl ′ ℓ ′ are addressed in general. This work also considers form-factor effects and improves upon the earlier work [350].

Radiative corrections for η (′)
Dalitz decays The knowledge of the radiative corrections for the neutral-pion Dalitz decay turns out to be very useful in different contexts, and it also serves as a starting point to the calculation of the corrections for the corresponding η (′) decay channels. These come with several additional challenges that one needs to overcome in order to be able to extract valuable information on the TFFs from data. Naive radiative corrections for η → e + e − γ [351] were indeed closely inspired by the solution for its neutral-pion counterpart, and only the mass of the decaying meson would be numerically addressed. However, one needs to take the situation a bit more seriously. The larger rest mass of the decaying η (′) has further consequences. The η-meson mass already lies above the muon-pair and pion-pair production threshold, which has some consequences for the treatment of virtual corrections. But what brings in a real complication is the fact that η ′ lies above the lowest-lying resonances ρ and ω. The decays and the radiative corrections then naturally become sensitive to the widths and shapes of the resonances. And this is a crucial difference compared to the π 0 case, where the TFF dependence was mild and could even be neglected in certain applications. Another technical complication then resides in the fact that η (′) mesons carry a nontrivial strangeness content that has to be addressed as part of the TFF treatment.
A more advanced approach to the η (′) → ℓ + ℓ − γ radiative corrections [352] then improves on the naive approach for η → e + e − γ in the following way. Regarding the vacuum polarization contribution, the muon-loop and hadronic contributions were added, the 1γIR contribution at one-loop level was included, and the form-factor effects (including the complete bremsstrahlung contribution) were calculated. Moreover, no assumptions regarding the final-state lepton masses were considered, so the muon channels could be consistently treated as well. An example of the resulting radiative corrections can be seen in Fig. 25.
Let us note that the developed exact treatment of the form-factor effects can also be used in the π 0 Dalitz decay. There it introduces circa 1% correction to the correction itself and is thus, as expected, negligible.

Radiative corrections for
The next technical complication in the calculation of radiative corrections to differential decay widths of Dalitz decays appears when the final-state photon is replaced by a massive particle. The first example of such a case that we discuss here appears outside the meson sector: the Dalitz decay of the neutral Sigma baryon. Radiative corrections within the soft-photonapproximation approach have been available for some time now [353]. The calculation of the bremsstrahlung contribution including the effects of hard-photon emission proceeds along similar lines as in the previous cases, although now the analytical integrals need to be generalized for a nonvanishing (Lambdahyperon) mass. It turns out that some of the integrals do not have analytical solutions anymore with such a generalization. The counterpart of the (until now triangular) 1γIR contribution of Fig. 23 are here twophoton-exchange box diagrams involving hadronic form factors; see Fig. 26. Although quite nontrivial to evaluate and contributing to the electron-positron asymmetry, they do not affect the one-fold differential decay width in terms of the invariant electronpositron mass squared, and thus the extraction of the Σ 0 → Λγ * transition form factor. Here lies a profound difference between the kinematical regimes of decay and scattering since, in the latter, the twophoton-exchange contribution matters, as discussed in Sec. 3.
Similarly to the π 0 -Dalitz-decay case, the (magnetic; electric is negligible) form-factor dependence is rather mild, and one can come up with precise predictions for Σ 0 → Λe + e − and Σ 0 → Λγ branching ratios [354]. Regarding the importance of the QED NLO radiative corrections (see Fig. 24), it was found, using again the simple consideration that the effect on the TFF-slope extraction is about −3.5% in absolute figures [354]. This is quite a significant value while extracting a quantity that is estimated to be of a size 1.8(3)%: When ignoring the QED effects, one would thus end up with a negative magnetic TFF slope, although again of a reasonable magnitude.

Radiative corrections for
Once the bremsstrahlung contribution to the Σ 0 Dalitz decay is established, the results can be readily used for the ω → π 0 ℓ + ℓ − decay. Here, various extractions of the ωπV correlator [355][356][357][358] are in tension with the simplest theoretical models, and one could wonder if applying the previously omitted QED radiative corrections (they were not available until very recently) could improve the situation. The corrections are, as in the previous cases, negative over a large region of the Dalitz plot (see Fig. 24), and one can thus expect the data points to be pushed upwards accordingly after the QED effects are subtracted. Needless to say, these corrections are more significant in the electron channel. A proper analysis is, of course, in place and should be part of the future experimental analyses of more precise data. The one-photon-inclusive NLO QED radiative corrections beyond the soft-photon approximation were calculated as a bachelor project in 2021 [359].

Radiative corrections for
The K + → π + ℓ + ℓ − decays allow us to access the K + → π + γ * TFF. Extracted both in the electron and muon channels, it allows us to test the leptonflavor universality. Most importantly, however, studying the TFF provides us with valuable information about weak transitions in the low-energy QCD sector. Extracting the relevant related hadronic parameters from experiment then requires, as expected, the knowledge of radiative corrections in order to subtract the QED effects that are present. Recently, following a number of previous extractions [360][361][362], the TFF parameters traditionally called a + and b + were measured by the NA62 Collaboration [363]; it is worth mentioning that radiative corrections were part of this analysis at the MC level [364,365].
Regarding the radiative corrections related to the lepton part of the K + → π + ℓ + ℓ − process, there are no additional complications, and one can use the same approach as in the previous cases once the form factor is suitably expanded to account for hadronic effects in the bremsstrahlung contribution beyond the softphoton approximation. On the other hand, the novel nontrivial obstacle is the fact that the involved kaon and pion are charged and can thus both also radiate a bremsstrahlung photon. Moreover, the transition K + → π + γ * γ ( * ) needs to be considered. Finally, virtual corrections that complement the bremsstrahlung from the meson part are somewhat nontrivial. Everything then should be ideally expressed in terms of a single form factor, the one that is extracted from data. This then allows for factorizing it out or for further iterations.
6.6 Radiative decay of K + → e + νγ with TREK/E36 at J-PARC The TREK collaboration (Time Reversal Experiment with Kaons) was formed in 2005 to pursue a program to search for T-violation in stopped K + decays [366]. When by 2010 it became clear that the required kaon beam intensity would not be achievable in a timely manner, the TREK collaboration still pursued a program on sensitive beyond-the-Standard-Model physics, with experiment E36, which was proposed in 2010 using only 50 kW beam power [367][368][369].
The focus of E36 was on testing lepton universality in the SM, by a precision measurement of the ratio of two-body decay widths, R K = Γ(K e2 )/Γ(K µ2 ). The experiment was approved in 2013 and constructed at the J-PARC K1.1BR beamline in 2014, reusing the toroidal spectrometer setup from KEK-PS experiment E246 [370] and the CsI(Tl) barrel [371]. A new scintillating fiber target, a spiral fiber tracker [372], and redundant PID systems (an aerogel counter array [373], a time-of-flight system, and a leadglass shower calorimeter in each sector [374]) were implemented. Data taking was limited to less than three months in 2015. First E36 results have been published on radiative kaon decays [375,376], with the result for R K being expected to become available in 2023. The E36 setup was also used for an in-situ search for dark photons, or more generally for light neutral bosons decaying into e + e − pairs over a broad mass range of 10-200 MeV/c 2 . The analysis, which is still ongoing, combines a charged particle (µ + , π + , or e + ) in one toroidal sector with two clusters from e + e − pair detection in the segmented CsI(Tl) calorimeter. Preliminary results from the light boson search have been reported [377], and the final sensitivity is expected to be reached in 2023.
The study of radiative kaon decay allowed to extract the branching ratio for the structuredependent (SD) process resulting in a hard photon radiated nearly in back-to-back kinematics with the positron. The process is accompanied by the internal bremsstrahlung (IB), predominantly of soft photons nearly in the direction of the positron. In the overlapping regions both processes appear, and the IB process is calculated and simulated as a background, and conventionally included in the two-body branching ratio for K e2 . The latest analysis also considers internal bremsstrahlung in the SD process K SD e2γ . The implementation of the IB radiative process in the MC simulation follows the scheme introduced by Gatti [378], allowing for analytic expressions of the differential branching ratio distributions to be used in an event generator. The experiment data were analyzed in terms of the lepton momentum in a narrow range. The radiative correction affecting the K e2 and K e2γ branching ratios were mostly due to the narrow acceptance chosen for the lepton momentum. In the ratio R γ = Γ(K SD e2γ )/Γ(K e2 ), the IB corrections partially cancel in the acceptance ratio. The first publication [375] accounted for the IB acceptance effect only in the monochromatic two-body decay K e2 . In the most recent result [376], the IB scheme was also implemented to correct the radiative yield in the three-body decay K e2γ . The obtained ratio is consistent with a recent lattice-QCD calculation [379], a gauged nonlocal chiral quark model (NLχQM) [380], however at tension with ChPT at order O(p 4 ) [381,382] and O(p 6 ) [383], and a previous measurement by KLOE [384] at the 4-σ level.
The latter data have been used by NA62 [385] to correct the universality ratio R K for the SD background, which should be revisited. 6.7 Radiative decay width ρ − → π − γ with COMPASS at CERN The COMPASS collaboration (COmmon Muon and Proton Apparatus for Structure and Spectroscopy) operated a multi-purpose, high-resolution and highrate capable spectrometer [386] at the CERN SPS M2 beamline with muon and hadron beams in the energy range 100-280 GeV. A program pursued since the very beginning some 20 years ago focuses on ultra-peripheral reactions of pion and kaon beams on atomic nuclei, in the domain of momentum transfers where the electromagnetic interaction dominates -so-called Primakoff reactions. It is intended to continue the program with a focus on kaon-induced reaction with AMBER [176], operating at the same site as COMPASS and as its successor experiment.  Fig. 27: Contributions to the cross section for the reaction π − γ → π − π 0 from the kinematic threshold, where the chiral anomaly dominates, to the region of the ρ(770) resonance. Adapted from the respective figure in [387]. The simplified model is used here to display the conceptual structure; for the data analysis to extract the chiral anomaly and the radiative width of the resonance, the dispersive framework [388] is used.
One of the channels of interest is π − γ ( * ) → π − π 0 , the photon from the exchange with the nucleus being quasi-real. According to the equivalent-photon formalism of Weizsäcker and Williams, for a 190 GeV pion beam on nickel, these photons cover a wide range of momenta of 0.1 to 2.0 GeV in the rest frame of the incoming pion. Thus the full low-momentum range for the reaction becomes accessible and covers the kinematic region in which the ρ(770) resonance dominates in the s-channel, π − γ ( * ) → ρ − → π − π 0 . The cross section in the low-energy region, governed by the chiral anomaly up to the resonance region, is depicted in Fig. 27. The strength of the resonance, as it appears in the reaction, depends on the radiative coupling ρ → πγ, entering in inverted kinematics. The measurement of the cross section thus allows for a precision determination of the radiative width of the ρ(770).
The reaction is affected by higher-order electromagnetic processes that have to be corrected for. A dominant contribution comes from a process with a t-channel exchange photon, coupling with the incoming photon to the outgoing π 0 (thus π 0 → γγ in inverted kinematics). Corrections due to photon loops and real-photon emission have been calculated and found to be on the level of a few percent, thus are to be taken into account for a precision analysis. The formulae have been given in [389] and revised [390]. The full implementation in terms of an event generator including bremsstrahlung photons has been started.

Summary and outlook
In this section, we have mainly presented results regarding the so-called inclusive radiative corrections. These are suitable to estimate the effect of QED corrections on the differential decay widths (and consequently on the extraction of TFFs or their parameters), assuming that the emission of additional photons has been ignored during the analysis. Any additional assumptions or cuts that are made in the experiment should be reflected in the calculation, so a close experiment-theory collaboration is essential for preparing the suitable tailored approach.
As pointed out earlier, the most desirable way to implement radiative corrections is as early as at the Monte Carlo level. Already at that stage, one might typically observe that there are events missing in the MC, and the agreement with data is not satisfactory. This leads to larger uncertainties and extracted parameters containing unsubtracted QED parts.
The corrections in Figs. 24 and 25 may seem large at first. But one needs to realize that the Dalitz decays are dominated by a one-photon exchange (i.e., by the low-x region due to a pole at vanishing photon virtuality). Integrating over the Dalitz plot then leads to the fact that the QED corrections at the dilepton production threshold will be receiving the largest weights leading to the determination of the corrections to the integrated decay width that are ∼ 1% as one would expect. These resulting values can then be independently checked by employing different methods (see, e.g., Ref. [341]), and it was done so for π 0 and Σ 0 Dalitz decays arriving at an exact agreement.
Note that using the soft-photon approximation, the corrections are typically negative all over the Dalitz plot, and the effects of the hard-photon corrections are instrumental for obtaining the observed positive total corrections (see, e.g., Ref. [353]).
We have stated examples of processes for which the radiative corrections have been worked out or used in the experiment. From the mentioned examples, it is straightforward to see the importance of radiative corrections and that with improving sensitivity and higher statistics of upcoming experimental set-ups, these will become increasingly indispensable in order to extract meaningful results from data. The list of light-hadron decays for which the exact NLO QED radiative corrections are available can thus be expected to get longer.

Summary and Conclusions
In this document we have outlined the state of the field for radiative corrections in ℓp scattering, QED corrections to deep-inelastic scattering, and in meson decays. In the case of ℓp scattering, the formalism of these corrections are laid out in detail at LO, NLO, and NNLO from a theoretical perspective. The experimental observables in unpolarized and polarized cross-section measurements are described in detail, in particular polarization transfer observables and SSAs.
We particularly emphasize the case of the TPE, as this is an active area of work for both theoretical and experimental efforts. While the TPE is believed to be the cause of the proton form factor discrepancy, and there has been significant theoretical effort to explain the discrepancy, this has not been conclusively demonstrated in experiment. It is of paramount importance that new TPE measurements using charge asymmetric beams be performed to verify existing calculations. Simultaneously, theoretical efforts using a variety of techniques must continue to be developed to inform experimental results.
As part of the theoretical effort to treat TPE, it is incumbent upon the community to develop event generators that can be easily integrated into analysis pipelines. These generators should be ideally capable of NNLO calculations that can capture the leading logarithms beyond NNLO. Such future generators, or improvements upon existing generators, are necessary for modern precision experiments. Generators are an active area of research in the community and efforts to continue this development are a welcome contribution to the community.
Beyond the scope of ℓp scattering experiments, there is more work to be done in e + e − annihilation experiments. While the complete two-loop QED corrections have been calculated, there are future experimental measurements of the Z 0 mass planned which require calculations of QED initial state corrections of up to 6th order.
For mentioned light-hadron decays, detailed calculations of radiative corrections involving twophoton-exchange contributions and hard-photon corrections were briefly presented. Future high-precision, high-statistics experiments will require further developments in calculations of higher-order corrections. It is evident that future work will result in exact NLO calculations of other meson decays.
While there is still much work to be done, there has been significant progress made in the area of precision, percent-level experiments, and NNLO calculations of their associated radiative corrections. Several future experiments are planned that will test calculations of radiative corrections, and all future experiments require modern treatments of their corrections. Continued experimental and theoretical efforts in the field are crucial as we progress into the next generation of precision physics measurements.

Acknowledgments
We