Soft-photon radiative corrections to the $e^- p \to e^- p l^- l^+$ process

We calculate the leading-order QED radiative corrections to the process $e^- p\rightarrow e^- p l^- l^+ $ in the soft-photon approximation, in two different energy regimes which are of relevance to extract nucleon structure information. In the low-energy region, this process is studied to better constrain the hadronic corrections to precision muonic Hydrogen spectroscopy. In the high-energy region, the beam-spin asymmetry for double virtual Compton scattering allows to directly access the Generalized Parton Distributions. We find that the soft-photon radiative corrections have a large impact on the cross sections and are therefore of paramount importance to extract the nucleon structure information from this process. For the forward-backward asymmetry the radiative corrections are found to affect the asymmetry only around or below the 1\% level, whereas the beam-spin asymmetry is not affected at all in the soft-photon approximation, which makes them gold-plated observables to extract nucleon structure information in both the low- and high-energy regimes.


I. INTRODUCTION
Double-virtual Compton scattering (dVCS) on a proton, the process γ * p → γ * p with initial and final virtual photons (γ * ), is a prime process to study and test models describing the electromagnetic structure of the nucleon beyond the information contained in the elastic form factors.
At low-energies, it allows to extract nucleon structure constants, which enter the expansion of the nucleon Compton amplitude. The real Compton scattering (RCS) limit, the process γp → γp with both photons real, has been used over many years as an experimental tool to access the nucleon electromagnetic polarizabilities, see e.g. Ref. [1][2][3][4] for reviews. The virtual Compton scattering (VCS) process, γ * p → γp with initial space-like virtual photon, which can be accessed as a subprocess of the e − p → e − pγ reaction, has also been studied extensively over the past three decades to access the generalized nucleon polarizabilities [1,[5][6][7]. These structure quantities allow to obtain, through a Fourier transform, a spatial representation of the deformation of the charge and magnetization distributions of the nucleon under the influence of an external static electromagnetic field [8].
The most general case of a double-virtual Compton process, with both initial and final virtual photons, has until now been studied only in special limits. The most useful extension is given by the forward doublevirtual Compton scattering (VVCS) process, where the initial and final photons have the same non-zero spacelike virtuality. In contrast to the processes discussed above, the forward VVCS process is not directly measurable. It enters however in the leading hadronic corrections to the muonic Hydrogen Lamb shift and hyperfine splitting. The interest in its improved estimate was spurred in 2010 by the ultra-precise determination of the proton charge radius from muonic Hydrogen Lamb shift measurements [9], which reported a radius value which was 4% smaller than the 2010 recommended value by the Committee on Data for Science and Technol-ogy (CODATA) [10] based on results from electron-proton scattering and ordinary Hydrogen spectroscopy measurements, and represents a 7σ difference. Over the past decade, major progress has been made in resolving this puzzle, see Refs. [11][12][13] for some recent reviews. The dominant theoretical model error in the extraction of the proton charge radius from muonic Hydrogen Lamb shift measurements to date results from the subtraction function entering the VVCS process [14][15][16]. At second order in the photon virtuality, this function is constrained by the magnetic polarizability, which is determined experimentally [17]. To fourth order in the photon virtuality, one low-energy constant in this subtraction function is at present empirically unconstrained [18], and one relies on chiral effective field theory calculations [15,19] or phenomenological estimates [20]. In Ref. [21] it was proposed to access this low-energy nucleon structure constant empirically through the forward-backward asymmetry in the e − p → e − p l − l + process, with l = e or l = µ. The dVCS amplitude contributing to that process, γ * p → γ * p, has an incoming photon with space-like virtuality and an outgoing photon with time-like virtuality.
A second kinematical region in which the virtual Compton processes are being used as a prime tool to study the partonic structure of the nucleon is at high energies, for near-forward kinematics, either through the e − p → e − pγ process with initial space-like photon with large virtuality, the deeply-virtual Compton scattering process (DVCS), or through the di-lepton photoproduction process γp → l − l + p process with outgoing time-like photon with large virtuality, the time-like Compton scattering (TCS). In such kinematical regime, pertubative Quantum Chromo Dynamics (QCD) allows to express the proton structure entering the DVCS and TCS processes through Generalized Parton Distributions (GPDs), which access the correlation between the longitudinal momentum distribution of partons in a proton and their twodimensional transverse spatial distributions. We refer the reader to Refs. [22][23][24][25] for the original articles on GPDs and to Refs. [26][27][28][29][30][31] for reviews of the field. Accessing the resulting three-dimensional momentum-spatial distributions of valence quarks in a nucleon through exclusive processes has been one of the driving motivations for the JLab 12 GeV upgrade [32]. Furthermore, accessing the sea-quark and gluonic structure of nucleons and nuclei through such processes is one of the main science questions that will be addressed at the future Electron-Ion Collider (EIC) machine [33].
A further extension of either the DVCS or TCS process in the high-energy near-forward region has been proposed through the e − p → e − pl − l + reaction (with l − either an e − or µ − ), which accesses the double deeply virtual Compton scattering (DDVCS) process with incoming space-like photon and outgoing time-like photon. The DDVCS process is of particular interest as it allows to extend the DVCS beam spin asymmetry measurements, which directly access GPDs, into the so-called ERBL domain [34,35]. A feasibility study of the DDVCS experiment has shown that the SoLID@JLab project with its high-luminosity and large acceptance is very promising to perform such measurements [36].
In order to use the e − p → e − pl − l + reaction as a tool of proton structure, it is imperative to quantitatively estimate the QED radiative corrections to this process, which is the main objective of the present work. Our work extends previous studies of radiative corrections for the VCS process [37], as well as more recently for the TCS process [38][39][40]. In Ref. [40] it was found for the γp → l + l − p process that the relevant asymmetries to extract the real and imaginary parts of the TCS amplitudes, the forward-backward and beam-helicity asymmetries, are nearly unaffected by the radiative corrections. In contrast the TCS cross sections receive sizeable corrections: in the low-energy region up to 10%, and in the highenergy kinematical region up to 20%. As for the single space-like or single time-like Compton scattering cases, it is crucial to have a good quantitative understanding of radiative corrections also in the double-virtual case in order to be able to extract relevant information about the proton structure from future experimental data. As a first estimate of the size of radiative corrections we will use the soft-photon approximation in this work. We distinguish between three different, gauge-invariant types of corrections, from which one contributes to the VCS case, a second one contributes to the TCS case, both of which are obtained as limits of our work, and a third type of correction which is new for the double-virtual case. We study the size of these corrections on the level of unpolarized cross sections as well as on the forward-backward and beam-spin asymmetries.
The outline of the present paper is as follows. In Sec. II we introduce the relevant Feynman diagrams at tree level. We distinguish between two different contributions: the Bethe-Heitler and the Compton scattering processes. In Sec. III we introduce the two different nucleon structure models which we use to describe the dVCS amplitude. In the low-energy regime we calculate the contribution from the Born process in terms of the protons form fac-tors as well as the ∆(1232) resonance excitation in combination with a low-energy expansion of the dVCS amplitude. In the high-energy regime we use the QCD factorization theorem to express the dVCS amplitude in terms of GPDs. In Sec. IV we calculate the virtual radiative corrections in the soft-photon approximation from the three gauge invariant types of contributions. We give analytic expressions for the finite and infrared divergent parts of all three contributions in terms of a factorizing contribution on the level of cross section. In Sec. V we calculate the contribution to the cross sections due to soft-photon bremsstrahlung. Taking real radiation into account, we cross-check analytically the cancellation with the infrared divergences from the virtual corrections. In Sec. VI we show our results for the observables in both the low-energy and high-energy kinematical regimes. We conclude in Sec. VII. Technical details are discussed in two appendices.

II. DI-LEPTON ELECTROPRODUCTION AT TREE LEVEL
In this work we study the di-lepton electroproduction process as a probe of proton (N ) structure, with l − either an e − or a µ − , where the quantities in brackets denote the particle four-momenta. At tree level, we distinguish between three different contributions, which we denote as the spacelike (SL) and timelike (TL) Bethe-Heitler (BH) processes, see Fig. 1, as well as the double virtual Compton process (dVCS), see Fig. 2. Figure 1. Tree level QED diagrams contributing to the e − p → e − pl − l + process. We distinguish between the spacelike (left) and the timelike (right) Bethe-Heitler processes. The crossed diagrams, for which in the spacelike process the order of the vertices on the produced di-lepton line are interchanged, and for which in the timelike process the order of the vertices on the electron beam line are interchanged, are not shown.
To specify the kinematics, it is useful to introduce the following four-momenta: (2) Figure 2. Tree level diagrams for the Compton scattering. The blob represents the (elastic and inelastic) interaction of the virtual photon with the nucleon. Figure 3. Planes defining the scattering angles which characterize the e − p → e − pl − l + process. The angles Φ and φ * l are defined with respect to the blue plane, which is the scattering plane of the virtual photons with four-momenta q and q .
The process (1) is defined by eight kinematical invariants, which we choose as: where Φ denotes the angle of the initial electron plane relative to the production plane. Furthermore, the angle θ * l (φ * l ) denotes the polar (azimuthal) angle respectively of the negative lepton in the rest frame of the l − l + lepton pair. In Fig. 3 we show the three different scattering planes defined by these angles.
We denote by m the mass of the electron, by m l the mass of the produced leptons, and by M the mass of the proton. The on-shell relations of the external particles are therefore: and the invariant s is obtained from the Lab electron beam energy E e as s = M 2 + m 2 + 2M E e .
The matrix element for the spacelike Bethe-Heitler (BH,SL) process (left diagram in Fig. 1) is given by: while the timelike Bethe-Heitler (BH,TL) process (right diagram in Fig. 1) is given by: In Eqs. (5,6), h(h ) denote the helicities of the intial (scattered) electrons, h − and h + are the helicities of the produced lepton pair, and s(s ) are the helicities of the intial (final) proton respectively. Furthermore, Γ α is the electromagnetic nucleon vertex given by: where F 1 (F 2 ) are the Dirac (Pauli) form factors of the nucleon respectively. The matrix element for the double virtual Compton scattering (dVCS) process ( Fig. 2) is expressed as: where M µν denotes the Compton tensor which depends on the model to describe the interaction of photons with the nucleon, and which will be specified below.
In the case of e − e + production we have to take into account that the electrons with momenta k and l − are indistinguishable. Thus for e − e + production, we have to consider, besides the direct (dir) contribution of Eqs. (5,6,8), also the contribution of all exchange (ex) diagrams where both electrons in the final state are interchanged. The Bethe-Heitler matrix elements corresponding with these exchange terms are given by (note that this only contributes in the case m l = m): and the exchange term corresponding with the dVCS matrix element is given by: To ensure the Pauli principle, one has to anti-symmetrize the amplitude under exchange of both electrons in the final state. Therefore, the full matrix elements for the e − p → e − pe − e + process is obtained as difference between the amplitudes for direct (dir) and exchange (ex) diagrams: while for µ − µ + production only the direct diagrams contribute. The fully differential cross section for the e − p → e − pl − l + process is given by where dΩ * l refers to the phase space of the produced lepton of the dilepton pair in the l − l + rest frame, and where β s ll is the lepton velocity in the l − l + rest frame:

III. MODELS FOR THE DOUBLE VIRTUAL COMPTON AMPLITUDE
The double virtual Compton tensor M µν entering Eq. (8) is calculated from the process We show the Feynman diagram for this process in Fig.  4. The blob in this diagram represents the interaction of the incoming and outgoing virtual photons with the proton. In the following we use the average photon (q) q, µ q ′ , ν p p ′ and proton (P ) momenta, The general double virtual Compton tensor M µν can be constructed using q µ , q µ , p µ , g µν and γ µ as building blocks. From these blocks, one finds 34 independent tensors with two indices [41]. Using gauge invariance it was shown that the number of independent amplitudes can be reduced from 34 to 18 [41]. However it was realized in Ref. [41], that there is in general a problem in such representation. For specific kinematical points the 18 tensors become linearly dependent and therefore do not form a basis at these specific points anymore. As a result the corresponding Compton amplitudes display kinematic singularities at these points. To bypass this problem, Tarrach introduced an overcomplete basis by introducing three additional tensors. Such an overcomplete basis does not have any kinematical constraints and is valid in the whole phase space. It was realized in Ref. [42] that the kinematic singularities and constraints of the Compton amplitude in a minimal basis are due to the Born terms, in which the intermediate state in the Compton process in Fig. 4 is a nucleon, and that for the non-Born contributions a minimal tensor basis consisting of 18 structures free of kinematical singularities and constraints exists.
In this work, we will only need the helicity-averaged amplitude, which is described by 5 independent tensors, and can be expressed as, following the notations of [42]: where T µν i are the spin-independent and gauge invariant tensors, symmetric under exchange of the two virtual photons, and are given by: T µν 3 = q 2 q 2 g µν + q · q q µ q ν − q 2 q µ q ν − q 2 q µ q ν , Furthermore, in Eq. (17), the invariant amplitudes B i are functions of four Lorentz invariants, with ν ≡ q · P/M .
In order to specify the double virtual Compton amplitude, we need to model the internal structure of the nucleon. In this work, we will consider two different models, which are tailored for applications in two different energy regimes. In a low-energy model, which is motivated for applications to describe the hadronic structure in precision atomic physics measurements such as the Lamb shift or hyperfine splitting in muonic Hydrogen, we consider the photons to interact with the nucleon and its lowest excitation, the ∆(1232) resonance. In a high-energy model, in which at least one of the photons is highly virtual, we use perturbative QCD which allows to factorize the Compton process on the nucleon in terms of a Compton amplitude on the quark convoluted with the amplitude to find the quarks inside the nucleon. The latter is parameterized through generalized parton distributions (GPDs).
A. Low-energy double virtual Compton amplitude

Born diagrams
In the low-energy regime, we describe the Compton tensor in terms of the leading Born (B) amplitude, given in terms of the proton form factors. The amplitude can be calculated from two Feynman diagrams shown in Fig. 5 (upper panel). Its contribution is given by where Γ µ i (Γ ν f ) are the initial (final) state proton vertices. Note that the FFs entering Γ ν f correspond with a timelike virtuality. For the numerical evaluation of these FFs we use the paramaterization of Ref. [43], which allows the analytical continuation based on dispersion relations into the unphysical part of the timelike region, 0 < q 2 < 4M 2 . Note that in this region no direct experimental extraction exists.

∆-pole model
In addition to the Born amplitude, we need a model for the non-Born contribution at low energies. The covariant baryon chiral perturbation theory (BChPT) provides a systematic framework for the calculation of the double virtual Compton scattering process, see Ref. [18]. The latter work has shown that BChPT is fully predictive at orders O(p 3 ) and O(p 4 /∆), in which p stands for a small momentum and with ∆ ≡ M ∆ − M the excitation energy of the ∆(1232) resonance. The O(p 3 ) contribution comes from the pion-nucleon (πN ) loops, and the O(p 4 /∆) contribution comes from the Delta-exchange (∆-pole) graph, which is shown in Fig. 5 (lower panel), and the pion-Delta (π∆) loops.
For the near-forward real Compton cross section (i.e. integrated over dilepton phase space), it was found that around W = 1.25 GeV the Born + ∆(1232)-pole contribution reproduces a full dispersive calculation based on empirical structure functions within an accuracy of 5% or better [21]. As we consider in this work kinematics around the ∆(1232) resonance, we will study as a first step the effect due to radiative corrections on the ∆(1232)-pole contribution.
The dominant contribution is coming from the magnetic dipole transition FF G * M . In the following we use only that dominant contribution, corresponding with the leading term in the so-called δ-expansion [46], to calculate observables, i.e. we set G * E = G * C = 0.

Low-energy expansion
The non-Born part of the dVCS amplitudes, denoted asB i , can be expanded for small values of ν, q 2 , q 2 , and q · q , with coefficients given by polarizabilities. The relations between these low-energy coefficients and the polarizabilities measured through real Compton scattering (γp → γp) and virtual Compton scattering (γ * p → γp) have been given in [18].
A special limit of the double virtual Compton process is given by its forward limit, denoted by VVCS, which corresponds with q = q and p = p. This limit is of particular importance as it enters the two-photon hadronic corrections to the electronic and muonic Hydrogen energy levels. The helicity averaged VVCS process is described by two invariant amplitudes, denoted by T 1 and T 2 , which are functions of the two kinematic invariants, Q 2 and ν, as: withĝ µν ≡ g µν − q µ q ν /q 2 ,p µ ≡ p µ − p · q/q 2 q µ , and where α em = e 2 /4π 1/137. The optical theorem allows to express the imaginary parts of T 1 and T 2 as: where F 1 , F 2 are the conventionally defined structure functions parameterizing inclusive electron-nucleon scattering, depending on Q 2 and x ≡ Q 2 /2M ν. The twophoton exchange correction to the µH Lamb shift can be expressed as a weighted double integral over Q 2 and ν of the forward amplitudes T 1 and T 2 [14]. Using the empirical input of F 1 and F 2 , the ν dependence of T 2 can be fully reconstructed using an unsubtracted dispersion relation, whereas the dispersion relation for T 1 requires one subtraction, which can be chosen at ν = 0 as T 1 (0, Q 2 ). The subtraction function is usually split in a Born part, corresponding with the nucleon intermediate state, and a remainder, so-called non-Born part, denoted byT 1 (0, Q 2 ). The Born part can be expressed in terms of elastic form factors and is well known, see e.g. [4] for the corresponding expressions. The non-Born part cannot be fixed empirically so far. In general, one can however write down a low Q 2 expansion ofT 1 (0, Q 2 ) as: where the term proportional to Q 2 is empirically determined by the magnetic dipole polarizability β M 1 [17].
Theoretical estimates for the subtraction term were given at order Q 4 in heavy-baryon chiral perturbation theory (HBChPT) [15], in BChPT, both at leading order (LO) due to πN loops, and at next-to-leading order (NLO), including both ∆(1232)-exchange and π∆ loops [18,19], as well as extracted from superconvergence sum rule (SR) relations [20]. The different estimates forT 1 (0) are compared in Table I. Even for these theoretically well motivated approaches, the spread among the different estimates is quite large. The resulting uncertainty due to this subtraction term constitutes at present the main uncertainty in the theoretical µH Lamb shift estimate. To reduce such model dependence, the dilepton electroproduction process on a proton has been proposed in [21] as an empirical way to determineT 1 (0).
As the forward VVCS process of Eq. (25) is a special case of Eq. (17), one can express the subtraction function entering the hadronic corrections to the µH energy levels as [18]:T where both non-Born amplitudesB 1 ,B 3 are understood in the forward limit (q = q ), i.e.B i (0, q 2 , q 2 , q 2 ) for i = 1, 3. In order to specifyT 1 (0, Q 2 ) up to the Q 4 term, we use the low-energy expansion in k ∈ {q, q } of the amplitudesB 1 ,B 3 [18]: where β M 2 is the magnetic quadrupole polarizability determined from real Compton scattering [47], and β M 1 (0) is the slope at Q 2 = 0 of the generalized magnetic dipole polarizability which is accessed through virtual Compton scattering, see Ref. [7] for a recent review. While the terms of O(k 0 ) and O(k 2 ) in the low-energy structure of the amplitudeB 1 at ν = 0 are empirically constrained from real or virtual Compton scattering, the low-energy constant b 3,0 is not determined empirically so far because the tensor structure T µν 3 in Eq. (18) decouples when either the initial or final photon is real. As such the lowenergy constant b 3,0 is the main unknown to date in the empirical determination ofT 1 (0). Below, we study the sensitivity of the e − p → e − pl − l + process, including the soft-photon radiative corrections, to this low-energy constant.

B. High-energy double virtual Compton amplitude in terms of GPDs
For the high-energy Compton scattering we calculate the Compton tensor in terms GPDs using pertubative QCD. This can be done by calculating the leading-order handbag diagrams shown in Fig. 7. For the evaluation of these diagrams, we will need the kinematic Lorentz invariants ξ and ξ , defined as: Furthermore, we introduce the two light-like vectorsp µ and n µ withp · n = 1, which are related to the four momenta P µ andq µ as: . The variablesξ andξ are related to the invariants ξ and ξ introduced in Eqs. (30,31) as: Although at leading-twist, corresponding with the kinematical regime for whichM 2 q 2 , one has we will keep in the following analysis the (small) difference in the kinematical quantities. The leading twist-2 double deeply virtual Compton scattering (DDVCS) amplitude on a proton is given by: where the coefficient functions C ± (x,ξ ) are defined as: and where the lightlike four-vectorsp and n are obtained from Eqs. (32,33) as Furthermore, for the purpose of studying the influence of the radiative corrections on the DDVCS observables at small values of −t, we will only consider the contribution of the dominant GPD H in our study below. For the numerical evaluation, we will use the GPD parametrizations from the VGG model [26,[48][49][50], summarized in Ref. [30], in terms of a double distribution (DD) and socalled D-term contribution (D), as: with the double distribution part for the proton given by the weigthed sum of the light quark flavor distributions: The isoscalar D-term contribution, is directly related to the subtraction function in a dispersive framework for the Compton amplitude. For its evaluation, we use the dispersive estimate of Ref. [51]. In order to satisfy exact electromagnetic gauge invariance for both incoming and outgoing virtual photons in the DDVCS process, we generalize the procedure introduced in Ref. [49] to add transversal correction terms which are formally of higher-twist as follows: where the transverse part ∆ ⊥ of the four-momentum transfer to the nucleon is defined as: Using the identities one immediately verifies that both q µ M µν DDVCS = 0 and q ν M µν DDVCS = 0. Using the parameterization of Eq. (41) for the GPD H in terms of a double distribution and a D-term part, the evaluation of the amplitude in Eq. (37) involves a principle-value integral which can be evaluated numerically, for the case 0 < ξ < ξ, as: with the singlet GPD defined as:

IV. VIRTUAL SOFT-PHOTON CORRECTIONS
In this work, we evaluate all one-loop virtual photon radiative corrections to the e − p → e − pl − l + process in the soft-photon approximation. This limit is defined by the scaling of the loop momenta: we only account for the regions of integration where the loop momentum l scales as: where λ is a small parameter compared to all external scales. We then calculate all contributions only up to order λ. The resulting corrections factorize in terms of the tree-level amplitude, which shows that this is a gaugeinvariant subset of the full one-loop corrections. From all soft-photon contributions, one can then further distinguish between three gauge invariant subsets: We give analytical expressions for the corrections of all three types. In order to regularize the infrared divergences coming from the integration over the soft-photon loop momentum l we use dimensional regularization [52]. We therefore perform the loop integration in D = 4 − 2 dimensions. Infrared (IR) divergences manifest themselves as 1/ IR poles in the regularized expressions. We are using the on-shell renormalization scheme. In addition to the diagrams with virtual soft photons we also have to consider infrared divergent counter terms. Those counter terms are introduced to regularize ultraviolet (UV) divergences (which manifest themselves as 1/ UV poles), which due to the on-shell renormalization condition can also carry IR divergences. Those need to be included in the calculation in order to get a finite result in the end.
We will subsequently discuss the virtual radiative corrections to the spacelike Bethe-Heitler process, the timelike Bethe-Heitler process, and the double virtual Compton process.
A. Corrections to the spacelike Bethe-Heitler process

Contributions of class (a)
In this section we calculate the soft-photon corrections for which the soft-photon is attached to the electron line. This corresponds to the left diagram in Fig. 8. In the following we suppress helicity states in all spinors to make the formulas more compact and better readable. The first diagram in Fig. 8 is given by, using Feynman gauge: which reduces in the soft-photon approximation to Here and in the following C 0 denotes the scalar one-loop three-point function. We give an analytic expression of that function for the two different cases we need in this work in Appendix B.
In addition to the contribution of Eq. (50), we also have to include the vertex counter term, which we show in Fig. 9. We are using the on-shell subtraction scheme, in which the counter term is defined to fix the electron charge e at q 2 = 0. In the soft-photon approximation one has to extract only the IR divergent piece of the full expression, as has been done in Ref. [38]. To calculate the vertex counter term we consider the decomposition into the two form factors F e D and F e P , where q = k − k . In this decomposition only F e D (q 2 ) is divergent. The renormalization constant Z 1 of the vertex is therefore given by yielding for the renormalized vertexΓ Since we work in the soft-photon approximation, we only extract the infrared divergent part of the full oneloop renormalized vertex which can be found for example in Ref. [37] and find After adding the vertex counter term to Eq. (50) and evaluating the three-point function, the infrared divergent part of the virtual correction to the cross section of the spacelike process is given by and the finite part by where β Q = 1 + 4m 2 Q 2 . Note that here and in the following, we define δ to be the correction on the level of the cross section, not the amplitude. This corresponds to taking twice the real part of the correction on the level of the amplitude.
For the crossed diagrams with l − and l + interchanged we find the same result as in Eqs. (55) and (56). In the limit of a small electron mass, i.e. m Q 2 , the correction simplifies to (57)

Contributions of class (b)
Here we calculate all contributions to the spacelike process, for which the soft-photon is attached to the dilepton line. The Feynman diagram corresponding to this correction is shown in Fig. 8 on the right. The matrix element is given by which in the soft-photon approximation reduces to As for the class (a) contribution we need to include counter terms. In addition to the infrared divergent piece of the vertex counter term as given by Eq. (54), we also need the counter term of the fermion self energy, which is shown in Fig. 10. To first order in α em , the self-energy of a fermion with mass m f and momentum k is calculated as .
(60) Eq. (60) has a UV divergence, which needs to be subtracted by an appropriate counter term. In the on-shell scheme this counter term is fixed by requiring that the fermion self-energy Σ(k ) has a pole at k 2 = m 2 f with residue equal to one. This fixes the wave-function renormalization constant Z 2 and the mass renormalization constant Z m f : The evaluation of Σ(k ) and its derivative, results in the renormalization constants The renormalized self-energy is then given bỹ and in the soft-photon limit, in which we only extract the IR divergence, we find (66) Adding the counter terms of the vertex and fermion self-energy to Eq. (59), we find for the total contribution the infrared divergent part and the finite part In the limit of small lepton masses, i.e. m 2 l s ll , we find (69)

Contributions of class (c)
In this section we calculate all diagrams, in which a soft-photon connects the electron line with the di-lepton line. We show the contributing diagrams of this class in Fig. 11. For the contribution of class (c) no counter term diagrams have to be considered. The first diagram in Fig. 11 is calculated as which in the soft-photon limit can be reduced to Evaluating the three-point function C 0 , we find that the infrared-divergent part is given by: and the finite part is given by: where We now consider two limits for this correction, in which the expressions simplify. The first limit corresponds to the case where the electron mass is small compared to all other scales. In this limit we find for the infrared divergent contribution and for the finite contribution If in addition to m 2 k·l − , also m l = m, i.e. considering electron pair production, we find The second diagram in the first row of Fig. 11 can be related to the previous one using Eq. (71) with the replacement l − → l + together with a sign change, Therefore the correction on the level of the cross section is given by The first diagram in the second row of Fig. 11 is given by which reduces to In this case, the second argument of the three-point function is positive. Therefore, an analytic continuation of this function to the timelike region has to be performed. This yields and We consider the two limits like before. In the limit of a small electron mass, we find for the infrared divergent contribution: and for the finite contribution: Considering electron production, m l = m, we find For the second diagram in the second row of Fig. 11 we can derive the correction in the soft-photon approximation from the previous result, leading to: Note that from Eqs. (77), (79), (87) and (88) we can see that the sum of all class (c) corrections is antisymmetric with respect to interchanging l + ↔ l − . This is in contrast to the contributions of class (a) and (b), which are symmetric with respect to the interchange of l + and l − .

B. Corrections to the timelike Bethe-Heitler process
In this section we calculate the soft-photon corrections for the timelike process. We show that they lead to exactly the same corrections as in the spacelike process.
The first diagram in Fig. 12 is given by: which in the soft-photon limit reduces to: (right) to the timelike BH process. One also has to consider the corresponding counter-term diagrams. For the correction, the crossed diagrams with ∆ and q interchanged yield the same result.
After adding the counter terms, on the level of the cross section the same correction as for the spacelike process is found: By the same argument, the second diagram in Fig. 12, including counter terms, yields the same correction as for the spacelike process from Eqs. (67) and (68): The same argument also applies to the four diagrams of class (c), shown in Fig. 13, which yield the same correction as for the spacelike process: (93)

C. Corrections to the Compton process
In this section we calculate the corrections for the Compton scattering, which in the soft-photon limit lead again to the same results as before for space-like and time-like Bethe-Heitler process. Therefore, on the level of the cross section, the correction in the soft-photon approximation can be factorized for the total process, and is given by: (94)  , and class (c) (lower two rows) to the dVCS process. One also has to consider the corresponding counter-term diagrams. The crossed diagrams with q and q interchanged yield the same result.

D. Sum of all virtual soft-photon corrections
Adding all contributions from classes (a), (b) and (c), we define the virtual soft-photon corrections on the cross section as: The correction can be separated in the IR divergent contribution: and a finite contribution: For convenience of the reader, we summarize all formulas derived in the previous sections: For the electroproduction of an e − e + pair, the formula simplify as: (104)

V. SOFT-PHOTON BREMSSTRAHLUNG
To cancel the IR divergences of the virtual soft-photon corrections, we need to include soft real radiation (soft bremstrahlung). On the level of cross section the IR divergences cancel, resulting in a finite physical result. The contribution due to soft bremstrahlung stems from Feynman diagrams in which an additional soft photon is emitted from an external fermion line. Denoting the momentum of the fermion line with l and the momentum of the soft photon by k γ , this corresponds to the amplitude: with the + sign, if the fermion is outgoing and the − sign if it is incoming, where Q f denotes the charge of the lepton, and where M 0 denotes the amplitude without soft-photon emission.
The evaluation of the bremstrahlung contribution requires integrating over the momentum of the unobserved soft photon up to an energy cutoff ∆E s . The integration has to be performed in a reference frame in which the dependence of the integral on the photon momentum is isotropic. The choice of this frame depends on the experimental condition. In the present work, we consider the process e − p → e − l − l + p where the dilepton momenta l − and l + are measured and the scattered proton with momentum p remains unobserved. Thus, the brem-strahlung integral has to be performed in a system in the rest frame of the unobserved proton and soft-photon.
Defining the missing momentum p m ≡ p + k γ , this frame is defined by the condition p m = 0. The bremstrahlung contribution to the cross section in this frame is given by: where the maximal soft-photon energy in that frame is denoted by ∆E s . The expression after performing the integration in Eq. (106) is lengthy and complicated in the general case. Here, we give explicit results only in the limit of a small electron mass, i.e. we only keep the logarithmic dependence on m. For the calculation, one considers the basic integral I ij corresponding to the interference of two terms like Eq. (105) from two fermion lines. This basic integral has been worked out in [53] for a generic lepton mass. In the limit m → 0, we find: Using this expression and the general one for finite lepton masses from [53], we can now easily perform the integration in Eq. (106). Let us stress again that we only keep the dependence of the electron mass m in the logarithms while for the lepton mass m l we keep the full dependence of the soft-photon integral. For the infrared divergent contribution we then find: while for the finite part we find: δ s;r ≡δ s;r a + δ s;r b + δ s;r c , E ± denotes the energy of the lepton with momentum l ± in the rest frame of the soft-photon and recoil proton, and whereẼ (Ẽ ) denotes the energy of the electron with momentum k (k ) in the same system. Adding Eqs. (96) and (108), we verify that the IR divergences from real and virtual soft-photon corrections cancel on the level of cross section.
As mentioned before, the integration of the soft-photon bremsstrahlung is performed up to a small energy cut-off ∆E s . This cut-off can be related to the experimental resolution of the detector. In the frame p m = 0 we find where to first order we have used p 2 m ≈ M 2 in the denominator, and where ∆p 2 m denotes the resolution in the missing mass squared. In order to express ∆E s in terms of Lab quantities, one needs to calculate the missing mass in that frame. Neglecting the lepton masses, we find where all quantities on the rhs have to be given in the Lab frame, where θ kk denotes the scattering angle between the incoming electron with momentum k and the outgoing with momentum k , and where θ ll denotes the Lab angle between the lepton pair momenta. Eqs. (115) and (116) allow one to express the maximal soft-photon energy ∆E s (defined in the system p m = 0) in terms of Lab quantities and detector resolutions.
In the following it will be convenient to express the energiesẼ ± ,Ẽ andẼ in terms of kinematic invariants. For the case of a large lepton mass, for which the formulas are lengthy and complicated, we use the formulas given in Appendix A and then boost to the rest frame of the recoiled proton and soft-photon to calculate the energies numerically. In the case of electron-pair production in which we can neglect the mass m compared to other quantities, the formulas become more compact. In that case, we also find more compact expressions for the bremstrahlung corrections. We find The energiesẼ ± ,Ẽ andẼ in the rest frame of the recoil proton + soft photon are given by: (120) VI. RESULTS

A. Observables
We use our setup to study the e − p → e − pl − l + process, including the first-order radiative corrections in both the low-and high-energy regimes. For both cases we study the effect of these corrections in the soft-photon approximation on the cross section, on the forward-backward asymmetry A F B , as well as on the beam-spin asymmetry A . These asymmetries are respectively defined as where dσ θ * l ,φ * l in A F B stands for the unpolarized cross section measured at lepton angles θ * l and φ * l (defined in the l − l + rest frame), and where dσ ± in A stand for the polarized cross sections for a polarized electron beam with helicity ±1/2 respectively. In the following we show plots ranging from θ * l = −180 • to θ * l = +180 • . This allows us to show forward and backward cross sections economically in one plot, since dσ(θ * l , φ * l +π) = dσ(−θ * l , φ * l ). The forward-backward asymmetry can therefore also be written as and, including radiative correction explicitly, it is given by: From Eq. (124) one can see that corrections that are symmetric under the interchange l − ↔ l + , corresponding with θ * l ↔ θ * l −π, drop out in the ratio. Therefore, to first order, only corrections of class (c) give a contribution to the asymmetry.
where A 0 F B denotes the uncorrected asymmetry. On the other hand the radiatively corrected beam-spin asymmetry (BSA), given by does not get modified in the soft-photon approximation, since the corrections are the same for both helicity cross sections, i.e. δ + = δ − , and therefore drop out in the ratio.
In the following, we show our numerical results for the e − p → e − pl − l + observables including the first order softphoton radiative corrections.

B. Results for dVCS observables in the ∆(1232) region
In this section we show our results in the low-energy regime in which we choose the dVCS center-of-mass energy W = 1.25 GeV. We model the dVCS amplitude in terms of the Born amplitude and the first proton excitation, the ∆(1232) resonance. As was found in Ref. [21], this model can reproduce the full calculation based on empirical structure functions from Ref. [54] with an accuracy in the few per-cent range for the process γp → e − e + p (i.e. for a real photon). Therefore, we can safely assume that the Born + ∆-pole model describes the dVCS amplitude sufficiently well also in the virtualphoton process for sufficiently small photon virtualities.
In order for the dVCS model to be also accurate for the e − p → e − pe − e + process, in which we need to antisymmetrize the full amplitude under exchange of both final electrons as given by Eq. (12), we choose the kinematics in such a way that also for the exchange dVCS amplitude the c.m. energy W ex remains in the ∆(1232) resonance region, and the photon virtualities entering the exchange process remain sufficiently small. As can be seen from Fig. 15 (upper panel), for the choice of an electron beam of 0.6 GeV, we find that W ex (blue dotted curve) is roughly of the same magnitude as W (dashed red curve), varying between 1.18 and 1.33 GeV as function of θ * e . Note that a larger electron beam energy leads to a larger value of W ex . From the lower panel of Fig. 15, we furthermore see that both photon virtualities in the exchange dVCS amplitude, denoted by Q ex (blue dotted curve) and s ll,ex (green dash-dashed curve), are both below 0.18 GeV 2 for the full range of the lepton angle θ * e . We are thus in a kinematic regime where we can study the sensitivity of the full amplitude to the low-energy constant b 3,0 described in Sec. III A.
Having studied the appropriate kinematics to describe both the direct and the exchange dVCS amplitude within the same model, we next explore the sensitivity of the e − p → e − pl − l + observables on the low-energy constant b 3,0 introduced in the low-energy expansion of Eq. (29). This low-energy constant is the main unknown in the determination of the O(Q 4 ) term of the subtraction func-tionT 1 (0, Q 2 ) entering the theoretical calculation of the µH Lamb shift.
In Fig. 16, we show the dependence on the lepton angle θ * l of the e − p → e − pl − l + differential cross section (upper panels), the forward-backward asymmetry (middle panels) as well as the beam-spin asymmetry (lower panels) for both e − e + and µ − µ + production (left and right panels respectively). We choose the kinematics as in Fig. 15. As can be seen from the upper panel, the interference between the dVCS process with the BH process amplifies the cross section for both e − e + and µ − µ + production by roughly a factor of two as compared with the BH process itself. Furthermore, the spread between the different theoretical estimates for the low-energy constant b 3,0 , as shown in Table I, increases the cross sections additionally by approximately 15% in both cases.
We also study the effect of the soft-photon radiative corrections on the cross section, as given by Eqs. (97) to (103), (109) to (112) and (117). For the real softphoton emission correction, we choose the soft-photon energy cut-off of ∆E s = 0.01 GeV, which corresponds to approximately 1.5% of the lepton beam energy. As can be seen from Fig. 16, the effect of the first-order radiative corrections is found to be quite sizeable on the level of cross sections. In the case of e − e + production the effect leads to a decrease of the cross section by around 30%, whereas for µ − µ + production it leads to a decrease of the order of 15%. Therefore, although the cross section by itself has a relatively high sensitivity on the low-energy constant b 3,0 , for an experimental extraction of b 3,0 the inclusion of the radiative corrections is imperative. A comparable importance of the radiative corrections was also found in the extraction of the proton generalized polarizabilities from the cross sections of the VCS process e − p → e − pγ [7,37].
The situation is different for the asymmetries. For the forward-backward asymmetry A F B , we find for the kinematics of Fig. 16 only a small sensitivity to the dVCS rad. corr. Figure 16. θ * l dependence of the e − p → e − pl − l + cross section (upper panels), forward-backward asymmetry (middle panels) and beam-spin asymmetry (lower panels) for e − e + production (left panels) and µ − µ + production (right panels) in the ∆(1232) region, for Φ = 90 • . The curves show the predictions for BH and BH + dVCS for two models showing the sensitivity to the low-energy constant b3,0. The black solid curves show the effect of the radiative corrections, for the hadronic model of the green dashed-dotted curves (these curves exactly coincide in the lower panels).  amplitude and its underlying hadronic model. However, this is mainly due to the choice of Φ = 90 • , for which the forward-backward asymmetry is completely dominated by the BH process. The sensitivity can be increased by varying Φ. In Fig. 17 we show the cross sections and forward-backward asymmetries for the same kinematics, but for smaller angles between the ( k, k ) and ( q, q ) scattering planes in Fig. 3: Φ = 30 • and Φ = 45 • . For these cases we find a 20% shift of the forward-backward asymmetry for the case including the ∆-resonance compared to the BH process by itself. Including the range of theoretical values for the dVCS low-energy constant b 3,0 , we find a further shift of the asymmetry of up to 5% on A F B , while the inclusion of radiative corrections is found to have a very small effect, around or below the 1% range on A F B .
For the beam-spin asymmetry A , we find a significantly higher sensitivity on b 3,0 than for the forwardbackward asymmetry, as shown in the lower panels of Fig. 16. Note that the result for Born + ∆-pole + b 3,0 (green dashed-dotted curves) and the result which in addition also includes the radiative corrections (black solid curves), coincide, since in the soft-photon approximation the radiative corrections drop out in the ratio of cross sections calculated for the BSA, as discussed above. Including the range of theoretical values for the dVCS lowenergy constant b 3,0 , leads to a shift in the BSA up to around 15% for e − e + production and up to around 10% for µ − µ + production.
As the BSA and A F B are basically not affected by the radiative corrections, a combined analysis of the cross section, the A F B , and the BSA holds promise to extract  Figure 18. Radiative corrections for the e − p → e − pe − e + process in the ∆(1232) region in the limit Q 2 → 0 for the comparable kinematic setup as was studied before for the γp → e − e + p process in Ref. [40]. The corrections are for soft-photon cut-off energy of ∆Es = 0.01 GeV.
In Fig. 18 we show in more detail how the radiative corrections to the e − p → e − pe − e + process vary when the intial photon approaches the real photon limit, i.e. Q 2 → 0. In this limit only the class (b) corrections contribute. We choose the kinematics comparable to Ref. [40], in which we studied the effect of radiative corrections for the process γp → e − e + p. In that study we found corrections of roughly 8% for the full set of one-loop QED corrections for the kinematics shown in Fig. 18. In the soft-photon approximation the corrections are somewhat over-estimated, as can be seen from the blue dotted curve corresponding with a correction of roughly 13%. Comparing the blue dotted curve with the red dashed-dotted one, we see that for quasi-real photons with a virtuality of 10 −4 to 10 −3 GeV 2 the inclusion of all corrections is important, and the description as a real photon underestimates the corrections by 10 − 15%. Note that in the region around Q 2 ≈ 0.03 GeV 2 the two outgoing electrons with momenta k and l − are becoming collinear. This explains the spiked behavior in the red dashed-dotted curve, since the logarithm with the argument proportional to the scalar product k · l − is becoming large.
In Fig. 19 (upper panels) we study in more detail the relative size of the radiative corrections due to the three different diagram classes (a), (b) and (c). While for e − e + production the corrections due to class (a) and (b) are dominant and negative, for µ − µ + production the main correction arises from class (a) as it involves the vertex correction on the beam electron, whereas the corrections between the produced µ − µ + pair are small and positive. Comparing left and right panel one can clearly see that the biggest difference between e − e + and µ − µ + production is due to the corrections of class (b), which correspond with the corrections from the produced dilepton pair. Furthermore Fig. 19 illustrates, as mentioned above, that the corrections of class (a) and (b) are symmetric under the interchange of l + and l − , corresponding to the angular shift θ * l → θ * l − π, while class (c) is anti-symmetric. As only class (c) contributes to the forward-backward asymmetry to lowest order, the smallness of the corrections of class (c) also explains why the A F B is largely unaffected by the radiative corrections.
Furthermore in Fig. 19 (lower panels), we show the sum of all three types of corrections also for a twice larger value of the soft-photon cut-off energy ∆E s . One notices the positive contribution to the cross section correction δ upon increasing the value of ∆E s .

C. Results for high-energy DDVCS observables
In this section we show our results for the e − p → e − pl − l + observables in the high-energy regime, in which we use GPDs to model the dVCS amplitude in the deeply virtual regime, the so-called DDVCS process, as described in Section III B. We will explore the sensitivity of this process to the modeling of the GPDs, in particular the D-term contribution, and quantify the effect of the radiative corrections in the soft-photon approximation.
As is conventional in the high-energy regime, we give the differential cross section with respect to the Bjorken scaling variable instead of the c.m. energy W describing the Compton process. Therefore, in this section we show differential cross sections with respect to the quantity ξ, which is related to the kinematical invariants through Eq. (30). The cross section differential w.r.t. ξ is related to the cross section differential w.r.t. W 2 as: In Fig. 20 we show a comparison of e − p → e − pl − l + cross sections for e − e + production (left panels) vs µ − µ + production (right panels). The cross sections are shown for an incoming electron beam energy of 11 GeV, which corresponds to the experimental setup of the CLAS12@JLab experiment and of the SoLID@JLab project. We show the cross section for ξ = 0.175, Q 2 = 2.75 GeV 2 , −t = 0.25 GeV 2 , Φ = 90 • and for three values of the di-lepton invariant mass s ll . While for e − e + production one observes a pronounced peak around θ * e ≈ −20 • for all three cross sections, the cross sections for µ − µ + production appear flatter. Furthermore, the cross sections for e − e + production are 4 to 15 times larger than the cross sections for µ − µ + production at the same value of s ll . This significant difference is due . Upper panels: θ * l dependence of the radiative corrections to the e − p → e − pl − l + cross section in the ∆(1232) region for e − e + production (left) and µ − µ + production (right), for the different classes of radiative corrections for ∆Es = 0.01 GeV. Lower panels: θ * l dependence of the total radiative correction for different values of ∆Es for both e − e + production (left) and µ − µ + production (right).
to the contribution of the exchange diagrams, which only contribute to e − e + production to satisfy the Pauli principle. Naively, one would expect the cross sections to be of roughly the same magnitude, since even for s ll = 0.5 GeV 2 the di-lepton invariant mass is 10 times larger than the µ − µ + production threshold, such that effects from the lepton mass don't play a crucial role. However, the anti-symmetrization of the final-state electrons for e − e + production yields a large contribution of the exchange diagrams, which increases with increasing values of s ll .
Furthermore, one can see from Fig. 20 that the BH process again serves as an amplifier of the dVCS process, and the BH+dVCS cross section is roughly 50% larger than the BH cross section. The cross section therefore has a strong sensitivity to the underlying GPD model. In particular through such interference, the D-term con-tribution to the GPD, using the dispersive estimate of Ref. [51], decreases the cross section by up to approximately 20% (green dashed-dotted curves vs red dashed curves).
In order to use the e − p → e − pl − l + as a tool to access the DDVCS amplitude, it is important to quantify the radiative corrections, which is an aim of this work. In Fig. 20 we show the impact of the radiative corrections on the cross section. For the real soft-photon emission correction, we choose the soft-photon energy cut-off of ∆E s = 0.05 GeV, which is roughly 1% of the e − p center-of-mass energy √ s = 4.64 GeV (corresponding to an electron beam of 11 GeV). We see from Fig. 20 that the soft-photon radiative corrections are very sizeable in the DDVCS regime, decreasing the cross sections by up to 50% for e − e + production and by up to 35% for  µ − µ + production (black solid curves vs green dasheddotted curves).
In Fig. 21 we show the θ * l dependence of the softphoton radiative correction factor on the e − p → e − pl − l + cross section in more detail, for the three values of s ll , corresponding to the cross sections shown in Fig 20. As mentioned above, using ∆E s = 0.05 GeV, the corrections in the DDVCS regime vary between −60% and −45% for e − e + production, while for µ − µ + production they are slightly smaller and vary between −45% and −30%.
In Fig. 22 (upper panels) we show the relative size of the corrections to the unpolarized cross sections stemming from the three different classes of diagrams for the central value of the squared di-lepton mass, s ll = 1.0 GeV 2 . We show both the virtual corrections (labeled "virt"), corresponding with Eqs. (97) to (103), as well as the virtual + real soft-photon corrections, with the real corrections given by Eqs. (109) to (112) and (117).
While the corrections of class (a) are just constant and negative, the corrections of class (b) and (c) have a nontrivial dependence on the di-lepton scattering angle θ * l . For class (b) this dependence is coming entirely from the real photon emission correction, since the di-lepton energies in the soft-photon frame depend on the lepton scattering angle θ * l , see Eq. (118). Comparing the soft-photon radiative corrections to e − e + and µ − µ + production we see that the largest difference is again coming from class (b), which is expected since it is most sensitive to the lepton mass m l . Let us note that, as discussed before for the low-energy case, classes (a) and (b) are symmetric with respect to the interchange l + and l − , while class (c) is anti-symmetric. Therefore to first order only class (c) contributes to the A F B .
In the lower panels of Fig. 22 we show the sum of all corrections for two different values of the soft-photon energy cutoff of ∆E s = 0.05 GeV (blue dotted curve) and ∆E s = 0.1 GeV (red dashed curve). The difference between both curves is a constant proportional to ln ∆E 2 s /m 2 . For the twice higher value of the softphoton energy we find in absolute values smaller corrections shifted by approximately 10%. Furthermore we see that the sum of all corrections is in absolute value smaller by more than 10% for µ − µ + production compared to e − e + production. As mentioned above the difference is coming mainly from corrections of class (b) which are (in absolute value) smaller for µ − µ + production.
In Fig. 23 we show the soft-photon radiative corrections for the e − p → e − pe − e + process, for roughly the same kinematics as in Ref. [40], in which we studied the effect of radiative corrections for the timelike Compton scattering (TCS) process, γp → e − e + p, with an on-shell incoming photon. In the soft-photon approximation the corrections to that process are equivalent to that of class (b) studied in the present work. Therefore, we find for corrections of class (b) the same order of magnitude of approximately -25 % as in [40] (blue dotted curve). Including also corrections of class (a) and (c) we find that the corrections increase with increasing Q 2 values, varying from -35 % to -50 % when varying Q 2 from 10 −5 GeV 2 to 10 −1 GeV 2 (red dashed curve). Furthermore, one observes around Q 2 = 1 GeV 2 the same spiked behavior as in the low-energy case in Fig. 18. The reason is again that in this kinematics the two outgoing electrons with momenta k and l − become collinear which leads to a large logarithm in the corrections of type (c). For an experiment such kinematic region should be avoided.
In Fig. 24 we show a comparison of both, beamspin and forward-backward asymmetries in the DDVCS regime, both for e − e + and µ − µ + production, and for a di-lepton invariant mass squared of s ll = 1.0 GeV 2 .
Studying the DDVCS process in this energy regime is . Upper panels: θ * l dependence of the radiative corrections to the e − p → e − pl − l + cross section in the DDVCS regime for e − e + production (left) and µ − µ + production (right), for the different classes of radiative corrections for ∆Es = 0.05 GeV. Lower panels: θ * l dependence of the total radiative correction for different values of ∆Es for both e − e + production (left) and µ − µ + production (right).
of particular interest as it allows to extend the DVCS beam spin asymmetry measurements of GPDs into the so-called ERBL domain [34,35]. The BSA is proportional to the imaginary part of the DDVCS amplitude of Eq. (37), and allows to access the GPDs directly unlike the real part of the amplitude which depends on a convolution integral over the GPDs. The numerator of the BSA directly yields for both the cases ξ > 0 (Q 2 > q 2 ) and ξ < 0 (Q 2 < q 2 ): with −ξ < ξ < ξ. In Eq. (128) c is a known factor, originating from the BH amplitude dependent on the nucleon elastic form factors, and the ellipses stand for the subdominant contribution of GPDs beyond H singlet . As H singlet is an odd function in its first argument, we thus see that the BSA for the DDVCS process changes sign when crossing the point ξ = 0. The BSA for the DVCS and TCS limits have the same magnitude but opposite signs, expressing the fact that the GPD information content in both limits is the same.
Given that the real BH process does not yield a BSA by itself, we see from Fig. 24 that the BSA has a significant sensitivity to the GPDs, yielding asymmetries . Radiative correction for the e − p → e − pe − e + process in the TCS limit Q 2 → 0 in a kinematic setup comparable to Ref. [40].
between -25% and +10% for e − e + production and between -25% and -5% for µ − µ + production. The difference between both cases is mainly due to the effect of anti-symmetrization in both outgoing electrons for e − e + production. Furthermore, unlike the DVCS and TCS cases, the BSA for DDVCS is also sensitive to the D-term contribution to the GPD, as it also yields a contribution to the imaginary part of the DDVCS amplitude. By comparing the red dashed and black solid curves in Fig. 24, we notice that the sensitivity to the D-term induces a change of the BSA by 5% or more over a large angular range. As noticed above, the radiative correction drop out of the BSA in the soft-photon approximation. We also show the forward-backward asymmetry in Fig. 24, and notice that the anti-symmetrization induces already a large effect for the BH process itself (blue dotted curves in Fig. 24). Adding the dVCS contribution changes the forward-backward asymmetry by up to 25% over a large angular range, while the effect of radiative corrections (black solid curves) is in the few percent range only. The sensitivity on the D-term for the forwardbackward asymmetry is much smaller, comparing the curve including the D-term (green dot-dashed line) and the curve excluding the D-term (red dashed line) we find a difference of up to approximately 5%.
For the calculation of the e − p → e − pe − e + cross sections shown in Fig. 20 (left panels) and the corresponding asymmetries shown in Fig. 24 (left panels) we have to ensure that the model used for the dVCS amplitude is applicable for both the direct and the exchange terms. In Fig. 25 we show the two photon virtualities entering the dVCS tensor for the exchange diagrams. One can see that in the kinematics considered one of the two virtu-alities is around or above 2 GeV 2 for nearly all lepton angles, which corresponds with the lower limit for which the QCD-factorization in terms of GPDs is expected to hold. This justifies the use of the handbag description of Sec. III B in terms of GPDs also for the exchange term. Furthermore in Fig. 25 (lower panel) we show the scaling variables entering the DDVCS tensor for the exchange diagrams. While ξ ex andξ ex both are roughly constant as function of the di-lepton angle, close to the value of ξ = 0.175 for the direct diagram, one notices a large angular variation for ξ ex andξ ex . Compared to the constant value ξ = 0.0758 entering the direct diagram, ξ ex varies between −0.08 and 0.08. Thus the e − p → e − pe − e + process has the unique feature, due to the anti-symmetrization in both outgoing electrons, that by varying the di-lepton angle θ * e one performs a systematic scan in the scaling variable ξ ex in the ERBL domain of the GPDs.
In Fig. 26 we also show the di-lepton momenta and scattering angles measured in the Lab frame as function of the di-lepton rest frame angle θ * l . From the upper panel of that figure it becomes clear, that most of the region is experimental accessible. Except for angles around 0 • and ±180 • , over most of the di-lepton angular range the lepton momenta are larger than 0.1 GeV, which makes it quite feasible to detect the particles. The scattering angles measured in the Lab system are shown in the lower panel in Fig. 26. In the region where the Lab momenta of the di-lepton pair are larger than 0.1 GeV, they range from −50 • to 50 • . This initial study seems promising for a measurement of the e − p → e − pe − e + process over a large range of scattering angles.

VII. CONCLUSIONS
In this paper we studied the soft-photon radiative corrections to the process e − p → e − pl − l + , where l = e or l = µ. The process contains two distinct contributions: firstly the space-and time-like Bethe-Heitler processes which only depend on the nucleon elastic form factors, and secondly the double virtual Compton scattering process. The latter is sensitive to the underlying hadronic model describing the virtual photon-nucleon interaction, and a measurement of e − p → e − pl − l + observables can therefore be used to test and study nucleon structure models for different energy regimes. In the present work we studied the e − p → e − pl − l + process in two different energy regimes.
In the low-energy regime, in which the center-of mass energy is close to the ∆(1232)-resonance, and in which both photon virtualities are typically below or around 0.1 GeV 2 , we described the interaction using a ∆-pole model together with a low-energy expansion of the dVCS amplitude. This regime is motivated to better constrain the hadronic corrections to precision atomic spectroscopy. In particular for the muonic Hydrogen Lamb shift, the main hadronic unknown to date results from a low-energy nu-  . θ * l dependence of the e − p → e − pl − l + beam-spin asymmetry A (upper panels) and forward-backward asymmetry AF B (lower panels) in the DDVCS regime, for e − e + production (left) and µ − µ + production (right). Curve conventions as in Fig. 20.
cleon structure constant, denoted by b 3,0 , which enters the empirical determination of the O(Q 4 ) term in the subtraction function T 1 (0, Q 4 ) of the forward double virtual Compton amplitude. We found that the spread between the different theoretical estimates for the lowenergy constant b 3,0 increases the e − p → e − pl − l + cross section by approximately 15% both for e − e + and µ − µ + production. Furthermore, we also found that the beamspin asymmetry and the forward-backward asymmetry, resulting from an interchange in the kinematics of the produced di-lepton pair, are sensitive to the low-energy constant b 3,0 . For the beam-spin asymmetry the range of theoretical values for this low-energy constant leads to a shift in the asymmetry up to 15% for e − e + production, and up to around 10% for µ − µ + production.
A measurement of the e − p → e − pl − l + observables in this low-energy regime is thus promising to extract the nucleon structure constant, which could help to reduce the main uncertainty in the theoretical µH Lamb shift estimate.
For the high-energy deeply-virtual regime, we modeled the dVCS amplitude in terms of GPDs. We studied the sensitivity of the e − p → e − pl − l + process to the modeling of the GPDs, in particular the so-called D-term contribution. In kinematics of future experiments at JLab, we found that dispersive estimates for the D-term contribution to the GPDs induce around 20% change in the e − p → e − pl − l + cross section. Furthermore, we also found a large sensitivity to the GPD model for the beam-spin as well as the forward-backward asymmetry. The beam-spin asymmetry is of particular interest as it does not involve any convolution over GPDs, but is directly proportional to the GPDs, mostly in a linear way, through interference with the Bethe-Heitler process. For the e − p → e − pe − e + process, the beam-spin asymmetry has the unique feature, due to the anti-symmetrization in both outgoing electrons, that by varying the di-lepton angle one performs a systematic scan in the average quark momentum fraction in the ERBL domain of the GPDs, due to the exchange term.
In order to use the e − p → e − pl − l + process in either the low-energy or high-energy regimes as a probe of nucleon structure, we also studied the QED radiative cor- rections on the observables, in the soft-photon approximation. We found that the radiative corrections have a large impact on the cross sections. In the low-energy regime we find that these corrections lead to a decrease of the cross section of up to 30% for e − e + production, and around 15% for µ − µ + production. In the high-energy deeply-virtual regime, the corrections even range up to 50% for e − e + production, and around 35% for µ − µ + production in JLab kinematics. For the forward-backward and beam-spin asymmetries the situation is different. For the A F B the radiative corrections were found to affect the asymmetry only around or below the 1% level, whereas the beam-spin asymmetry is not affected at all in the softphoton approximation. A combined analysis of the cross section and of both asymmetries thus holds promise to ac-cess the hadronic structure information in both regimes. A next step to interpret future measurements of e − p → e − pl − l + observables, would consist in performing a full one-loop radiative correction calculation, beyond the soft-photon approximation. Such calculation can build upon the work of Refs. [39,40], in which such study was performed for the related γp → l − l + p process. The latter study has shown that the soft-photon approximation can be expected to somewhat over-estimate the full one-loop corrections on the cross sections, while the beam-spin and forward-backward asymmetries remain nearly unaffected by the radiative corrections.