Leading order corrections to the Bethe-Heitler process in the $\gamma p\rightarrow l^+l^-p$ reaction

The ratio of di-lepton production cross sections on a proton, using the $\gamma p\rightarrow l^+ l^- p$ process, above and below di-muon production threshold allows to extract the effective lepton-proton interaction, which is required to be identical for electrons and muons if lepton universality is exact. To test for a scenario of broken universality at the percent level, of the size which could explain the different proton charge radii extracted from electron scattering and from muonic hydrogen spectroscopy, we evaluate all one-loop QED corrections to this process, including the full lepton mass dependencies. We furthermore show that two-photon exchange processes with both photons attached to the proton line vanish after averaging over di-lepton angles, and estimate the relatively small radiation off the proton. We compare the full one-loop calculation with a soft-photon approximation of the same order, and present estimates for a planned experiment.

The proton radius puzzle, the discrepancy between extractions of the proton charge radius from electron scattering or electronic hydrogen spectroscopy on the one hand and from muonic hydrogen spectroscopy on the other hand, is not solved yet. The initial discrepancy amounted to around 5.6 σ when comparing both values: the extraction using elastic electron scattering, from which the A1@MAMI Collaboration reported the value R E = 0.879 (8) fm [1,2], and the muonic hydrogen spectroscopy, which reported the value R E = 0.84087 (39) fm [3,4], with more than an order of magnitude higher precision. This puzzle has spurred a lot of activity, resulting in a new round of experiments in the field, which are crucial to scrutinze and improve our understanding of systematic errors in such precision measurements. Recent measurements using electronic hydrogen spectroscopy [5,6] as well as new electron scattering experiments [7] have each reported support for both large and small values of the proton charge radius. Several further experiments [8,9] have reported preliminary results and are at the stage of final analysis.
The possible explanations of electron versus muon discrepancy by a new physics scenario require to give up lepton universality, since the Standard Model has the same tree-level couplings for all leptons. In this context new experiments have been proposed to measure the proton form factors using muon beams: the MUSE@PSI experiment [30,31], which is ongoing with the aim to compare low-energy elastic electron scattering versus elastic muon scattering off the proton, and the COMPASS@CERN [32] experiment, which plans to measure the elastic muon scattering off a proton using a 100 GeV muon beam.
In Ref. [34], the authors suggested to test lepton universality by comparing the cross sections for electron-and muon-pair production in the reaction γ p → l + l − p. According to the findings of Ref. [34], a measurement of the cross-section ratio below and slightly above the di-muon production threshold can test lepton universality without having to rely on the precision which can be achieved for an absolute cross section measurement. It was found that the difference between both values for the extracted proton charge radius amounts to an effect of 0.2 % on this cross section ratio. Such measurement would therefore allow to test lepton universality at the 3σ significance level, if one is able to measure such cross section ratio with a precision of around 7 × 10 −4 . An upcoming experiment at MAMI is presently conceived to perform such measurements [35].
To make a conclusive statement from such kind of measurements, it is clearly necessary to include higher-order corrections. In our previous work [36], we have estimated the one-loop and higher-order QED corrections in the soft-photon approximation. We have found that the effect of radiative corrections is of the order 1 % on the cross section ratio, and thus significantly larger than the 0.2 % effect between both values for the radius. Therefore, the knowledge of the full one-loop corrections becomes imperative for an interpretation of an upcoming experiment.
In this work, we extend our calculation of the QED corrections in the soft-photon approximation and present a full one-loop QED calculation keeping all terms in the lepton mass. Moreover, we estimate the size of radiative effects on the proton side and indicate the vanishing two-photon exchange effect after the integration over the lepton-pair angles.
The outline of the paper is as follows. In Sec. II, we introduce the Bethe-Heitler process at tree level and define the relevant kinematic variables. In Sec. III, we present details of the full one-loop QED calculation. We introduce the crossing relation to simplify the evaluation of the crossed diagrams from direct ones. We demonstrate that our full one-loop calculation reproduces the correct double-logarithmic behavior of the soft-photon approximation [36].
In Sec. IV, we present details of the numerical approach to this calculation, using the Mathematcia packages FeynArts and FormCalc, which serves to cross-check the analytic result. In Sec. V, we present the calculation of the real corrections for the emission of a soft photon, as done in our previous work [36]. We verify that the infrared divergences in these contributions cancel with the infrared divergent part of the virtual one-loop corrections. In Sec. VI, we give the expression for the full one-loop QED corrections and exponentiate parts associated with the soft photon contributions. In Sec. VII, we present an estimate of the radiative corrections on the proton side. We prove that two-photon-exchange diagrams do not contribute on the level of the cross section after integrating over the di-lepton phase space. In Sec. VIII, we present our numerical results, and show the effect on the absolute cross section as well as on the cross section ratio of muon-and electron-pair production. We conclude in Sec. IX. Several technical details are discussed in five appendices.

II. LEPTON-PAIR PRODUCTION AT TREE LEVEL
The Bethe-Heitler process at tree level is described by two graphs, see Fig. 1. We use p (p ) for the momenta of the initial (final) proton, and p 3 (p 4 ) for the momenta of leptons l − (l + ) respectively. The initial photon has momentum p 1 , and the virtual photon momentum in the one-photon exchange graphs of Fig. 1 is defined as p 2 = p − p . The Mandelstam variables for this process are defined as The on-shell condition for external particles implies: At leading order, the scattering amplitude M 0 is given by where the electromagnetic vertex Γ ν for the proton is expressed as with the proton's Dirac and Pauli form factors F D and F P , respectively.
The corresponding unpolarized differential cross section dσ 0 is given by dσ dt ds ll dΩ where E γ is the lab energy of the initial photon and Ω CM l + l − ll is the solid angle of the lepton pair in their center-of-mass frame, in which the lepton velocity is denoted by In Eq. (12), we average over all polarizations in the initial state and sum over the polarizations in the final state. We express the cross section as a product of hadronic and leptonic parts as dσ dt ds ll dΩ where the fine-structure constant is defined as α ≡ e 2 /4π ≈ 1/137. Furthermore, the unpolarized leptonic tensor at leading order L µν 0 (including the average over the initial photon polarization) is given by and the unpolarized hadronic tensor H µν by Using (11), the unpolarized hadronic tensor can be expressed as wherep ≡ (p + p )/2. We have defined with τ ≡ −t/(4M 2 ) and where the electric (G E ) and magnetic (G M ) form factors are defined as which are functions of the spacelike momentum transfer t.
A very compact expression for the contraction of the tree-level lepton tensor L µν 0 with the hadronic tensor H µν can be found in [33,34]. We show the expression in our notation in appendix A.
For the electric and magnetic proton form factors, which enter the total cross sections for lepton-pair production, we exploit the fit of Ref. [2], which is based on a global analysis of the electron-proton scattering data at Q 2 < 10 GeV 2 with an empirical account of TPE corrections.
In the experimental setup, when only the recoil proton is measured, one has to integrate (14) over the lepton angles: The kinematical invariant t is in one-to-one relation with the recoiling proton lab momentum p (or energy E ): whereas the invariant s ll is then determined from the recoiling proton lab scattering angle: where s can be expressed in terms of the initial photon-beam energy E γ as: In Ref. [34], the authors calculated the cross section ratio R between electron-and muon-pair production: which depends on the invariant mass of the lepton pair s ll , and a reference point s 0 ll to which the measurement is normalized.
The corresponding plot for the kinematical range accessible at MAMI is shown in Fig. 2.
The normalization is shown for the choice s 0 ll = s ll , i.e., at each point above the muon-pair production threshold the sum of the cross sections for muon-and electron-pair production is divided by the corresponding cross section for electron-pair production. In this plot, the blue curve describes the scenario, when lepton universality holds, i.e., G µ E = G e E , while the red curve corresponds to a case when lepton universality is broken by an amount of 1%, which would correspond with the difference in proton radii as extracted from muonic Hydrogen spectroscopy [3,4] and from electron scattering [1,2]. The blue band describes the 3σ deviation if this observable is measured with an absolute accuracy of 7 × 10 −4 . We will show in this work that radiative corrections shift this curve by more than 3σ, making their inclusion indispensable for a comparison with experiment.

III. CALCULATION OF THE ONE-LOOP CORRECTIONS
In this work, we calculate all one-loop leptonic corrections contributing to the γp → l + l − p in three independent setups: (i) an analytic calculation using the techniques of Integration-By-Parts (IBP) identities, (ii) a numerical calculation using the mathematica package FormCalc, which uses the Passarino-Veltman tensor reduction with numerical implementation of scalar integrals in LoopTools, (iii) a calculation of the self-energy and vertex diagrams with the help of projection techniques.

A. Crossing relations
One can find relations between a given loop diagram and its crossed counterpart. Let us assume that the first diagram (direct diagram) in Fig. 3 is given by Then the crossed diagram is given by where Γ c ( p 1 , p 2 , p 3 , p 4 , m) can be related to Γ d ( p 1 , p 2 , p 3 , p 4 , m) by reading Γ d ( p 1 , p 2 , p 3 , p 4 , m) backwards and changing the direction of the momentum flow, i.e.: Γ is always a product of an odd number of Fermion propagators, so we can pull out the minus sign and, instead, change m to −m: This relation also holds for the tree level. Interfering the loop-diagrams with the tree diagrams we get therefore: Summing over spins, we find: The interference with the crossed counter-part gives: and, summing over spins: where we used cyclicity and complex conjugation of the argument of the trace in the second step. Comparing Eqs. (33) and (35), one can easily verify that: Since an odd number of gamma matrices does not contribute to the trace, only terms with an even number of powers of m contribute. Therefore the substitution m → −m has no effect, and the crossing relation reads: .
Having calculated any diagram, one can therefore easily obtain the expression for its crossed counterpart by only exchanging p 3 and p 4 . Note that the complex conjugation in Eq. (36) only affects the algebraic, i.e. the trace over gamma matrices, and not the analytic part (the i prescription in the integrals) of the diagram. For the unpolarized cross section, which only involves the real part of such interference, this complex conjugation is of no relevance.
B. Lepton self-energy

Lepton self-energy at first order
We show the first-order lepton self-energy diagram in Fig. 4. In the following, we use for the calculation of all Feynman diagrams the Feynman gauge and dimensional regularization for ultra-violet (UV) divergences ( UV = 2 − d/2 > 0) and for infrared (IR) divergences ( IR = 2 − d/2 < 0). The self-energy is then given by: It can easily be reduced to where s ≡ (k ) 2 , and where the finite parts of the master integrals A 0 and B 0 are given by the expansions in Appendix B.
For k = m, we find: Note that we separate the UV-divergent part and imply only the regular part of the scalar integrals A 0 , B 0 and C 0 in all expressions in the main part of this paper.
The on-shell renormalization condition fixes the pole at (k ) 2 = m 2 with residue equal to one. This gives the renormalization constants Z 2 and Z m : The evaluation of Σ(k ) and its derivative, results in the renormalization constants: The renormalized self-energy is then given bỹ which reads as: +m

Self-energy diagram
In the on-shell scheme, only self-energy diagrams contribute with the virtual photon attached to an internal fermion line. In Fig. 5, we show the corresponding diagram, contributing to the Bethe-Heitler process. The amplitude is given by where the renormalized self-energyΣ is given by Eq. (46). The interference of the direct self-energy diagram with the lowest-order diagrams is given by: where with and a (1) Evaluating these two traces, we find for the direct self-energy diagram: where where L µν 0,d denotes the direct part of the tree-level lepton tensor, which is given by L µν 0,d = L 1;µν SE , and where T SE 1 is given by Eq. (A1).
In the limit for small lepton masses, the expression for the direct and crossed lepton selfenergy contributions keeping only terms with logarithmic terms in the lepton mass scaled by s ll reduces to:

Decomposition of the half-off-shell vertex
In order to calculate the one-loop vertex corrections on the lepton side, we need to evaluate the half-off-shell vertexΓ µ , shown in Fig. 6. It can be constructed with the following set of Lorentz structures: with Using the Dirac equation: we see, that k does not appear in the decomposition. All in all, we find 6 independent Lorentz structures in the decomposition of the vertex: with scalar form factors F i ± , s = (k ) 2 and projectors to the on-shell states: The structure proportional to F 3 ± does not contribute to the cross section, since q µ gets either contracted with the photon momenta or with the hadron tensor; both contributions give zero due to gauge invariance. However, for the proof of gauge invariance, we need to consider this structure.
One can construct projectors to extract the six form factors bŷ where the coefficients a For the form factors, we find: The corresponding expressions for the form factors F 3 + and F 3 − , which do not contribute to the observables due to gauge invariance, are given in Appendix D.

The on-shell limit and renormalization
In this section we show, that the form factors F 1 + , F 2 + and F 3 + reproduce the correct on-shell limit. The on-shell form factor can be decomposed into two structures, according For s = m 2 , we find for the on-shell expressions of the form factors: where we defined F 1 + (m 2 , q 2 ) has an infrared divergence, which arises from C 0 (m 2 , q 2 , m 2 , 0, m 2 , m 2 ).
We reproduce the on-shell Dirac form factor which is given by Moreover, we reproduce the Schwinger correction to the electron magnetic moment given by In the on-shell subtraction scheme, the vertex counter term is defined to fix the electron charge e at q 2 = 0. Only the Dirac form factor F 1 (q 2 ) is UV divergent, and one finds at q 2 = 0 the renormalization constant: This leads to the renormalized (on-shell) form factor:

Ward-Takahashi identity
As a further check of our one-loop expression for the half-off-shell vertex we check the Ward-Takahashi identity, which reads: where Σ(k) denotes the Fermion propagator. Contracting the vertex of Eq. (61) with the photon momentum q µ , we find: Plugging in the form factors (Eqs. (65)-(68), (D1), (D2)), we find: This is indeed equivalent to the right-hand side of Eq. (78): considering the one-loop self-energy of the fermion, as can be seen by using Eqs. (39) and (40).

Vertex diagrams
Having determined the half-off shell one-loop expressions for the vertex, we can now evaluate the the two vertex diagrams contributing to the Bethe-Heitler process, as shown in Fig. 7. The first diagram of Fig. 7 is given by Comparing this to the expression of M 0 , given by Eq. (10), 3 additional Lorentz structures appear. The interference of the first vertex diagram with the two tree level diagrams is given where a (1) a and Evaluating these 4 traces, we find for the first vertex diagram: where ll +M 2 s ll t 2 ll + 2ss ll t 2 ll + ts ll t 2 ll + us ll t 2 ll + 2tu 2 t ll + 2M 2 s 2 ll t ll + ss 2 ll t ll + us 2 ll t ll + 2M 4 tt ll −4M 2 tut ll + 4M 4 s ll t ll + 2u 2 s ll t ll − 2M 2 ss ll t ll − 3M 2 ts ll t ll − 6M 2 us ll t ll + 2sus ll t ll The second diagram of Fig. 7 is given by The interference of the second vertex diagram with the two Born diagrams is given bȳ where a (1) a (4) and These lepton tensors enter the calculation like before: where The contributions of the crossed vertex diagrams can be calculated by using the crossing relations derived in Sec. III A. Note that the replacement t ll → u ll also affects the scalar form factors.

Vacuum polarization at first order
We show the first-order vacuum polarization diagram in Fig. 8. The photon propagator can be written as: where D µν 0 is the leading order photon propagator and Π αβ is the vacuum polarization, which is given at first order in α by: Due to gauge-invariance, q µ Π µν = q ν Π µν = 0, and the vacuum polarization can be decomposed as: where Π(q 2 ) is given by [37]: where v is defined in Eq. (73).
The UV-divergence in Eq. (114) is removed by the renormalization constant Z 3 : which is fixed by requiring, that the renormalized vacuum polarizationΠ(q 2 ) has a pole with residue 1 at q 2 = 0: The renormalized vacuum polarization is then given by: The renormalized photon propagator is therefore given bỹ Note that, due to gauge invariance, only the term proportional to g µν contributes to a physical cross section.
We show the contributing vacuum polarization diagrams in Fig. 9. The amplitude for the sum of both diagrams is given by whereΠ(t) is given by Eq. (117) and M 0 by Eq. (10).
Note that the correction due to the vacuum polarization is the same for muon-and electron-pair production. It therefore drops out in the cross section ratio of lepton-pair production, considered in this paper. This remains also true when one considers hadronic vacuum polarization.

E. Lepton box diagrams and reproduction of soft-photon result
In this section, we will present the calculation of the lepton box diagrams, as shown in Fig 10. For the analytic calculation, we use QGRAF [40] to generate all diagrams in a representation, which can be processed with the algebra program Form [41]. To generate Integration-By-Part (IBP) identities, which are used to express a set of integrals in terms of a smaller set of so-called master integrals, we use the program Reduze 2 [42]. These master In the general case with a finite lepton mass, the calculation of the lepton box diagrams leads to very large expressions. Nevertheless, we were able to extract the analytical expressions with the setup described above.
Having calculated the leptonic tensor L µν , we were able to check gauge invariance by contracting with the external photon momenta, i.e., It turns out, that this identity is only fulfilled, once all diagrams (self-energy, vertex and box) and all counter-terms are taken into account. Therefore, only the renormalized lepton tensor is gauge invariant.
Here we give the leading contribution stemming from double logarithms, for the limit of small lepton masses m 2 → 0 and show that we reproduce the correct asymptotic behavior for s ll >> 4m 2 , which was derived in Ref. [36]. In this case, the contribution to the cross section can be written as where B and C are defined in Eq. (107), while the other coefficients are defined as: +4tu 2 s ll t ll + 6s 2 ts ll t ll − 3t 2 us ll t ll + 8stus ll t ll − C, andD is given by by making the replacements t ll → u ll and u → −u − s − t + 3M 2 in D.
By adding the contribution of self-energy, vertex and box diagrams, i.e Eqs. (57), (107) and (121), all mixed double-logarithmic terms involving the factor ln(m 2 /s ll ) drop out in the sum. Furthermore, the contribution proportional to ln(m 2 /s ll ) and ln 2 (m 2 /s ll ) factorizes in terms of the tree-level amplitude, such that we can write: where µ is the scale associated to the infrared divergence of the soft-photon loops. The dependence on µ cancels, once one takes real photon corrections into account.
We are thus able to reproduce the correct double-logarithmic behavior for the virtual corrections, which we derived in [36] in Eqs. (46,47), in the limit of β → 1, corresponding to s ll >> 4m 2 .
Furthermore, in Ref. [50] the radiative corrections to the dilepton pair photoproduction cross section were studied in the approximation of a small lepton mass. The correction is shown in Eq. (B1) of that paper. We compared our analytical result after expansion in m 2 and find exact agreement. To make the comparison, we had to make the identifications (kinematical quantities on left were introduced in Ref. [50]): Note that the L and H functions of Ref. [50] are given analytically by:

F. Integration over lepton angles
For an experiment which only measures the recoiling proton, we have to integrate the differential cross section over the lepton angles dΩ CM ll ll , as shown in Eq. (22). In order to perform the integration we have to express the kinematic invariants t ll and u in terms of the lepton angles θ ll and φ ll defined in the rest frame of the di-lepton pair: where The integration over the lepton angles θ ll and φ ll was done numerically in Mathematica.

IV. NUMERICAL CALCULATION OF VIRTUAL PHOTON CORRECTIONS
The numerical calculation of the virtual photon corrections were mostly done in Mathematica, using the packages FeynArts, FormCalc, FeynCalc and Looptools [43,44].
To obtain the leptonic tensor, we calculated the cross section of a truncated two-to-two scattering with an off-shell photon and an on-shell photon as incoming particles in dilepton production. This simplifies the calculation significantly since it can use the automated techniques to generate cross sections in FeynArts and FormCalc. Furthermore, one can directly apply the existing regularization and renormalization scheme from regular two-totwo scattering processes provided by FormCalc.
After obtaining the expression for the leptonic cross section, one can use the polarization vectors for the photons to construct the leptonic tensor structure. The renormalization constants are provided numerically within the FormCalc in the on-shell scheme, and are given in terms of one-and two-point functions and derivatives of these. After renormalization, the result can be tested to be UV finite by varying µ 2 .
After obtaining the expression for the leptonic cross section, one can use the polarization vectors for the photons to construct the leptonic tensor. We created a python script to replace all the photon polarization vectors into tensorial form in the format of FeynCalc notation. Note that after the contraction with the hadronic tensor, the scalar result is still infrared divergent. One has then to include the soft photon bremsstrahlung. For the regularization of infrared divergences, FormCalc introduces a photon mass. Cancellation of infrared divergences between the real and virtual corrections was tested by varying the photon mass.
Using the numerical result, we were able to check the analytic expression at a given phase-space point.

V. SOFT-PHOTON BREMSSTRAHLUNG
Having evaluated all one-loop virtual corrections, we have to take into account the radiation of real undetected photons to get an infrared finite cross section. We calculate real radiation only in the soft-photon limit and show the contributing diagrams in Fig. 11. Note that the diagram, where the photon is attached to the internal lepton line, does not con- tribute in the soft-photon limit. Denoting the momentum of the photon by k, the squared matrix element is given by: To calculate the contribution to the cross section, one then has to integrate over the undetected soft-photon energy up to a small value ∆E s , determined by the experimental resolution.
Due to the energy-momentum conserving δ-function, δ 4 (p 1 + p − p 3 − p 4 − p − k), the integration domain has a complicated shape in the lab system. The integration can be carried out in the rest frame S of the real (p 1 ) and virtual (p 2 ) photons, which is also the rest frame of the di-lepton pair and soft photon, defined by In such frame, the dependence of the integral with respect to the soft-photon momentum becomes isotropic. For the differential cross section, we then need to evaluate: where the integration is performed in the frame S.
The integrals are infrared divergent and can be carried out analytically after dimensional regularization. They have been worked out, e.g., in Ref. [38]. For the kinematics in system S, where the soft-photon momentum: with the lepton four-momenta: we obtain: where δ IR s;R is the infrared-divergent contribution due to real photon emission: and δ s;R is the corresponding finite part: The maximum value of the undetected soft-photon energy ∆E s is defined in the system S.
One can re-express it in terms of the detector resolutions. We consider the case of detecting the recoil proton only. The energy E and angle θ p of the scattered proton are measured in the lab frame. The missing mass M miss of the system is defined by where E s denotes the soft-photon energy.
The missing mass M miss is experimentally determined from the quantity: where τ is determined from the lab proton momentum by Eq. (23), and θ p is the experimentally measured recoil proton scattering angle in the laboratory frame.
For the process without radiation, this angle is given by Eq. (25), which can be equivalently obtained from Eq. (138) by the replacement M 2 miss → s ll : Combining Eqs. (138) and (139), we can express the soft-photon energy of Eq. (137) approximately as: Consequently, the experimental recoiling proton angular resolution, denoted as ∆θ p , determines the maximum value ∆E s of the undetected soft-photon energy, which enters the radiative correction of Eq. (135), as

VI. TOTAL RESULT AND EXPONENTIATION OF THE SOFT-PHOTON CONTRIBUTION
Adding the virtual one-loop corrections and the real soft-photon correction, we find a cancellation of all infrared divergences. The full result can be written as: where δ 1-loop is the finite part of the one-loop virtual corrections. In the limit s ll >> 4m 2 , one obtains from Eq. (123) the approximate result to double logarithmic accuracy: As shown in Ref. [36], we can exponentiate the double logarithmic part of δ exp , which in the soft-photon limit is given by: For s ll >> 4m 2 , this becomes By accounting for the exponentiation of large double-logarithms in m 2 /s ll , the full one-loop correction of the cross section is therefore given by: Note that, in Eq. (146) we have to subtract the leading double-logarithmic term −α/π(1 + β 2 )/(4β) ln 2 ((1 − β)/(1 + β)) from the virtual one-loop corrections and add the terms from the real soft-photon corrections that cannot be exponentiated, c.f. Ref. [36]. In the limit s ll >> 4m 2 , Eq. (146) is given by: in which the correct leading logarithm dependence from Eq. (123) is reproduced.
Furthermore, it was shown in Ref. [39] that the normalization factor F arises due to the physical assumption that in an experiment the sum of all soft-photon energies is smaller than ∆E s , instead of requiring that each soft-photon energy is individually smaller than ∆E s . Its leading correction from unity is given by: Although we account for the factor F explicitly, its deviation from unity is quite small: for s ll = 0.077 GeV 2 approximately −2.4 × 10 −3 for electrons and −8.5 × 10 −6 for muons.

VII. PROTON-LINE CORRECTIONS
Besides the one-loop leptonic corrections, we also estimate the one-loop hadronic corrections. As the latter are much smaller than their leptonic counterparts at low energies, due p 1 with the infrared-divergent contribution δ IR s;pR : wherev = 2 τ (1 + τ )/ (1 + 2τ ), and the corresponding finite part δ s;pR : The initial (β p ) and final (β p ) proton velocities in the lepton-pair c.m. frame are given by The soft-photon integral I (β p , β p ,v) can be expressed as with We provide an alternative expression for the integral I in the Appendix E.
Note the exact cancellation of the infrared divergences between the soft bremsstrahlung of Eq. (155) and the infrared part of the proton vertex correction coming from 2F 1 (t) with the renormalized form factor of Eq. (77), using the identification with v defined in Eq. (73).

VIII. RESULTS AND DISCUSSIONS
In Fig. 13 we show the radiative corrections to the cross section in the kinematical range of s ll between 0 and 0.08 GeV 2 and compare with our previous result in the soft-photon approximation. The muon threshold is at s ll = 4m 2 µ ≈ 0.045 GeV 2 (vertical dashed red line in Fig 13). We observe that the corrections for electrons are negative of order 10 percent, while the corrections for muons are positive of order 1 percent. The difference between the full one-loop calculation and the soft-photon approximation comes from terms which are not proportional to the double logarithmic form ln 2 (m 2 /s ll ) 2 . In Fig. 14 we show the effect of the full one-loop radiative corrections on the ratio defined in Eq. (27). The exponentiation has a considerably smaller effect on the full oneloop calculation, than on the soft-photon approximation. Furthermore, the soft-photon approximation clearly overestimates the effect of radiative corrections in this calculation. Taking radiative corrections into account, the ratio of Eq. (27) is now given by which depends on the measured invariant lepton mass s ll and the reference point s 0 ll , to which the cross section is normalized. δ e and δ µ are given by Eq. (146). One chooses s 0 ll < 4m 2 µ , such that the reference measurement is below the muon-pair-production threshold, and only electron pairs are created. In Fig. 15 we show the differential cross section ratio R of Eq. (27), including full oneloop QED corrections with ∆E s = 0.01 GeV. The radiative corrections to R are of the order of 1%. The red curve in Fig. 15 shows the scenario when lepton universality is violated by G µ E /G e E = 1.01, which is an effect of order 0.2%. Following Ref. [34], we use 3σ bands around the curves, with the experimental resolution σ = 7 × 10 −4 . One sees from this plot, that the inclusion of radiative corrections is indispensable, since the ratio of cross sections, defined in Eq. (27), is shifted to higher values by more than the 3σ band. The statement that lepton universality can be tested with a 3σ confidence level remains true if one adds radiative corrections as can be seen in Fig. 15.
We estimate the resulting correction on the proton side as a sum of soft bremsstrahlung the total hadronic radiative correction amounts to δ ≈ −(1.5 − 2) × 10 −4 in that range. The hadronic correction factor depends significantly on the momentum transfer t. We show this dependence for s ll = 0.077 GeV 2 and the photon beam energy E γ = 0.5 GeV in Fig. 17.
The resulting correction to the cross section is up to −3 × 10 −4 δ 0. We see from Consequently, our approximation of the proton as a point QED particle is reliable for the envisaged precision ≈ 7 × 10 −4 on the ratio of cross sections, which is required for a test of lepton universality.

IX. CONCLUSION
In this paper we calculated the first-order QED corrections to the Bethe-Heitler process in the γp → l + l − p reaction, keeping the full lepton mass dependence. This reaction may serve as a test for lepton universality as the authors in Ref. [34] pointed out, and such experiment is presently in the planning stage at MAMI.
The ratio of di-electron production cross sections above and below µ + µ − threshold shows a sensitivity of 0.2% when the difference between the larger proton charge radius from electron scattering is used versus the smaller proton radius which results from the muonic Hydrogen spectroscopy. Since the full one-loop radiative effects induce a correction around 1% on this same ratio, its inclusion is indispensable in this comparison.
The calculation was done in two independent setups, whose results were found to be in perfect agreement. Furthermore, our calculation reproduces the correct leading logarithms of the soft-photon approximation, and is in agreement with Ref. [50] in the limit of small lepton masses.
We also showed, that hadronic corrections are negligible at the required level of precision, due to the cancellation of box type graphs after integration over the lepton angles.
As a next step, we plan to extend our study to the reaction with an off-shell photon in the initial state, since the anticipated experiment at MAMI is designed to use a virtual photon.
Council (NSERC) of Canada. Our figures were generated using Jaxodraw [45], based on AxoDraw [46]. For the Mathematica plots we use the package MaTeX 2 .
Parts of this research were conducted using the supercomputer Mogon and/or advisory services offered by Johannes Gutenberg University Mainz (hpc.uni-mainz.de), which is a member of the AHRP (Alliance for High Performance Computing in Rhineland Palatinate, www.ahrp.info) and the Gauss Alliance e.V.