Leading-order QED radiative corrections to timelike Compton scattering on the proton

We evaluate the leading-order QED radiative corrections to the timelike Compton scattering (TCS) process $\gamma p \to l^- l^+ p$. We study these corrections in two energy regimes using different models for the TCS amplitude. In the low-energy regime we calculate the contribution due to the proton and its lowest-energy excitation, the $\Delta(1232)$ resonance. In the high-energy near-forward kinematical regime we calculate the TCS amplitude in a handbag approach in terms of Generalized Parton Distributions (GPDs). On the level of cross sections we find the QED radiative corrections to be in the $5 -10\%$ range in the low-energy regime and around $20\%$ in the high-energy regime. We show that in both the di-lepton forward-backward asymmetry as well as in the photon beam helicity asymmetry these corrections nearly cancel out, making them gold-plated observables to extract the real and imaginary parts of the TCS amplitude. We demonstrate in particular the sensitivity of these asymmetries on GPD parameterizations for a recent CLAS12@JLab TCS experiment.


I. INTRODUCTION
The virtual Compton scattering (VCS) process is a versatile tool to unravel the proton electromagnetic structure beyond the information contained in its elastic form factors. The e − p → e − pγ reaction which accesses the virtual Compton scattering process with an incoming photon with spacelike virtuality and outgoing real photon has been studied extensively both at low and high energies, see Ref. [1] for an early review. At low energies, it allows to access generalized polarizabilities of the proton [2], which encode the spatial deformations of the quark charge densities in a proton upon applying an external electromagnetic field. These observables have been extracted over the past two decades at the electron scattering facilities MIT-Bates, MAMI, and Jefferson Lab (JLab), see Ref. [3] for a recent review. At high energies and for near-forward kinematics, the e − p → e − pγ process is closely related to deepinelastic scattering. In this regime, pertubative Quantum Chromo Dynamics (QCD) allows to express the proton structure entering the deeply virtual Compton scattering (DVCS) process through Generalized Parton Distributions (GPDs), which access the correlation between the longitudinal momentum distribution of partons in a proton and their two-dimensional transverse spatial distributions. We refer the reader to Refs. [4][5][6][7] for the original articles on GPDs and to Refs. [8][9][10][11][12][13] 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 [14]. 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 [15].
The information accessed in the e − p → e − pγ reaction can be complemented through the γp → e − e + p reaction which accesses the timelike Compton scattering (TCS) process with incoming real photon and outgoing timelike photon, through production of a di-lepton pair. In the near-forward kinematical regime in which the timelike photon has a large virtuality, the non-perturbative information entering the TCS amplitude can also be expressed in terms of GPDs [16]. Furthermore, a combination of both DVCS and TCS observables allows for a stringent test of the applicability of the underlying QCD factorization theorem at these kinematics. Such measurement of the TCS process at large timelike virtuality has been proposed by CLAS12@JLab [17], and recently first data of this experiment have been reported [18,19].
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 spacelike photon and outgoing timelike photon. The DDVCS process is of particular interest as it allows to extend the beam spin asymmetry measurements of GPDs into the ERBL domain [20,21]. Recently 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 [22].
The measurement of the e − p → e − pl − l + process has recently also been proposed in the low-energy region in order to help reducing the theoretical model error in estimates of the hadronic correction to the muonic hydrogen Lamb shift. Over the past decade, measurements of the 2S-2P Lamb shift in muonic hydrogen (µH) have reported a proton charge radius with an order of magnitude improvement in its precision [23,24] as compared to the precision obtained in electron scattering. Initially these muonic atom extractions of the proton charge radius disagreed by around 5.6 standard deviations with the values obtained from energy level shifts in electronic hydrogen [25] or from electron-proton scattering [26,27], arXiv:2012.09565v1 [hep-ph] 17 Dec 2020 which triggered a lot of activity in recent years, see Refs. [28,29] for some reviews. A new Lamb shift measurement in electronic hydrogen [30] as well as a new experiment using electron scattering [31] are now both in support of the smaller proton radius value obtained by muonic measurements. To fully clarify the situation, further experiments with electron beams at MAMI and MESA [32], with muon beams at the Paul Scherer Institute [33] and at CERN [34], or through a direct comparison of γp → e − e + p versus γp → µ − µ + p [35] have been proposed or are underway. All these experiments aim at extracting the proton charge radius with improved precision. Experimentally so far, the most precise measurement is coming from the µH Lamb shift measurements, for which the proton form factors and polarizabilities are required as theoretical input to estimate the nextorder proton structure corrections. Those are at present the theoretical limitation on the proton radius extraction from Lamb shift measurements. It has been demonstrated however how these proton structure corrections can be empirically constrained through measurements of the forward-backward asymmetry in the e − p → e − pe + e − process [36].
To extract all of the above proton structure information from the single or double virtual Compton scattering process requires an estimate of the Quantum-Electrodynamic (QED) radiative corrections to these processes. For the e − p → e − pγ reaction, a detailed study of the radiative corrections has been performed [37] and has been applied to existing data [3]. In the present work, we report on the calculation of the leading-order QED corrections to the reaction γp → e − e + p in two different energy regimes, generalizing our calculation of Ref. [38,39] to include the TCS process. We study the radiative corrections on the cross section as well as the di-lepton forward-backward asymmetry and the photon beam helicity asymmetry. The present work will set the stage for a future study of the leading-order QED radiative corrections to the double VCS reaction e − p → e − pl − l + .
The outline of the present paper is as follows. In Section II we introduce the contributing Bethe-Heitler (BH) and TCS amplitudes to the γp → l − l + p reaction at tree level and define the relevant kinematic variables. In Section III we describe two different models for the double virtual Compton amplitude which are tailored for applications in two different energy regimes. In the lowenergy regime, motivated for applications to describe the hadronic structure in precision muonic atom spectroscopy, we calculate the contribution due to the proton and its lowest-energy excitation, the ∆(1232) resonance. In the high-energy near-forward kinematical regime in which at least one of the photons has a large virtuality, we use a QCD factorization theorem to describe the double virtual Compton amplitude on the proton in terms of a Compton amplitude on the quark convoluted with the non-perturbative structure of the proton encoded in the GPDs. In Sections IV, V and VI we calculate the virtual one-loop QED corrections due to vacuum polarization, as well as due to photons attached to the leptonic lines of the BH and TCS amplitudes. In both cases we calculate the corrections on the level of the amplitude, allowing for the calculation of polarized cross sections. In Section VII we calculate the soft-photon emission contributing to the process in which the di-lepton pair is measured. Including the soft-photon radiation gives infrared finite results for the observables. In Section VIII, we present our numerical results, and show the effect of the radiative corrections on the cross sections as well as on the forwardbackward and photon beam helicity asymmetries. We show results both in the ∆(1232) resonance region as well as in the kinematical regime of the CLAS12@JLab TCS experiment. For the latter, we show the sensitivity of the cross section and asymmetries on the underlying GPD parameterization. We conclude in Section IX.

II. BETHE-HEITLER AND TIMELIKE COMPTON SCATTERING PROCESSES AT TREE LEVEL
In this work we consider the process where the quantities in brackets denote the fourmomenta of the particles. We distinguish between two different contributions to (1), which are called the Bethe-Heitler (BH) process and the timelike Compton scattering process (TCS). We show the corresponding Feynman diagrams in  The blob on the proton line denotes the nucleon structure, whereas the blob on the lepton line denotes the QED amplitude, including radiative corrections.
The process (1) is defined in terms of three kinematic invariants: and two angles θ l and φ l , which are defined in the restframe of the di-lepton pair, with the polar angle θ l defined relative to the c.m. direction of q ≡ l + + l − . The matrix element of the BH process at leading order is given by: where m denotes the mass of the lepton and where the electromagnetic vertex Γ ν for the proton is expressed as with momentum transfer to the proton ∆ ≡ p − p, satisfying ∆ 2 = t, with M the proton's mass, and with F D (F P ) the Dirac (Pauli) proton form factors (FFs) respectively.
The general TCS matrix element is given by: where M µν is the Compton tensor which will be specified below, and which depends on the model used to describe the interaction with the proton. The unpolarized, fully differential cross section dσ 0 is given by where E γ is the lab energy of the initial photon, which is related to W as E γ = (W 2 − M 2 )/(2M ), and Ω * ll is the solid angle of the lepton pair in the l + l − center-of-mass frame, in which the lepton velocity is denoted by The tree-level amplitude M 0 is given by the sum of BH and TCS amplitudes: In Eq. (6), we average over all polarizations in the initial state and sum over the polarizations in the final state.

III. MODELS FOR THE DOUBLY VIRTUAL COMPTON AMPLITUDE
The doubly virtual Compton tensor M µν entering Eq. (5) is calculated from the process We show the Feynman diagram for this process in Fig.  2. The blob in this diagram represents the interaction of the incoming and outgoing photons with the nucleon. In the following we use the average photon (q) and proton (P ) momenta, q, µ q ′ , ν p p ′ Figure 2. Diagram representing the doubly virtual Compton process.
Although in this paper we study the process for a real initial photon, we will indicate below the extension to the general case of two off-shell photons.
The general doubly 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 [40]. Using gauge invariance it was shown that the number of independent amplitudes can be reduced from 34 to 18 [40]. The latter number corresponds with the minimal number of helicity amplitudes for a parity conserving process, which can be determined by accounting for the possible helicity states of the photons (3) and fermions (2). However it was realized in Ref. [40], 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 which do not have any kinematical constraints and are valid in the whole phase space. It was realized in Ref. [41] 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. 2 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 order to specify the doubly 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 highenergy 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 GPDs.
A. Born term and ∆-pole model at low energies At low photon energies, the leading Born (B) contribution to the Compton amplitude is described by two Feynman diagrams, shown in Fig. 3 (upper panel), in which a nucleon is propagating between both photon interactions. Its contribution to the tensor in Eq. (5) can be calculated as: where Γ µ i (Γ ν f ) are the initial (final) state proton vertices, given by analogous expressions as Eq. (4) in which ∆ is replaced by q (−q ) for Γ µ i (Γ ν f ) respectively. Note that the FFs entering Γ ν f correspond with a timelike virtuality. For the numerical evaluation of these FFs we use the paramaterization of Ref. [42]. This parameterization allows the analytical continuation based on dispersion relations into the unphysical part of the timelike region, 0 < q 2 < 4M 2 , in which no direct experimental extraction exists. In addition to the Born term, we also evaluate the matrix element of the leading contribution due to the ∆ resonance in the general case with two off-shell photons. The corresponding s-channel diagram is shown in Fig. 3 (lower panel) and its contribution to the tensor in Eq. (5) can be calculated as: In Eq. (12), Γ βµ γN ∆ andΓ αν γN ∆ refer to the γ * N ∆ vertex function and its adjoint respectively. They are shown in Fig. 4 and can be expressed in terms of three transition form factors as [43]: and its adjoint: where we defined Note that the FFs g M , g E , and g C appearing in Eq. (13) have spacelike virtuality (q 2 < 0), whereas the corresponding ones in the adjoint vertex of Eq. (14) have timelike virtuality (q 2 > 0). The form factors g M , g E , and g C can be expressed by the more conventional magnetic dipole (G * M ), electric quadrupole (G * E ), and Coulomb quadrupole (G * C ) transition FFs as: with the so-called Ash FFs parameterized, for spacelike virtuality Q 2 = −q 2 , through the MAID2007 analysis as [44,45]: with Q in GeV and the dipole FF G D (Q 2 ) = 1/(1 + Q 2 /0.71) 2 . For small timelike virtualities, 0 < q 2 < (M ∆ −M ) 2 , we can extrapolate the expressions for spacelike virtualities by the substitution Q 2 → −q 2 .

B. High-energy timelike Compton Scattering in terms of GPDs
For deeply-virtual Compton scattering, in which at least one of the photons has a large virtuality, we can express the doubly virtual Compton tensor using perturbative QCD in terms of GPDs. At leading order in the large virtuality, the deeply virtual Compton tensor can be calculated through the handbag diagrams, shown in Fig. 5. These handbag diagrams express the factorization of the process in terms of the Compton amplitude on the quark convoluted with the amplitude to find the quark in the nucleon, which is parameterized through the GPDs. To describe the deeply-virtual Compton process with two virtual photons, it is convenient to define the Lorentz invariants ξ and ξ as: To calculate the handbag diagrams, we first express the four-momenta P µ andq µ in terms of the lightlike four-vectorsp and n, withp · n = 1, as: . The variablesξ andξ are related to the invariants ξ and ξ introduced in Eqs. (17,18) as: We notice that the difference between the tilded quantitiesξ,ξ of Eqs. (21,22) and the quantities ξ, ξ of Eqs. (17,18) involve kinematical corrections due to the target mass M and momentum tranfer t. In the following, we will consider the Bjorken limitq 2 M 2 , in which:ξ With these kinematic definitions, the double deeply virtual Compton scattering (DDVCS) tensor at leading twist-2 can be expressed as [20]: where and where the lightlike four-vectorsp and n are obtained from Eqs. (19,20) as: Furthermore in Eq. (24), H, E,H, andẼ are the GPDs, which depend on the two quark momentum fractions x and ξ, and on the momentum transfer t. One can apply the above formula of Eq. (24) for the DDVCS tensor to two experimentally important limits. The first is the deeply-virtual Compton scattering (DVCS) process, in which the final photon is real (q 2 = 0), and the intial photon's virtuality is large, i.e.
for which one has: with Bjorken variable x B ≡ Q 2 /(2p · q). The second limit is the timelike Compton scattering (TCS) process with inital photon real (q 2 = 0), and large timelike virtuality, i.e. q 2 −t, for which one has: TCS : The TCS amplitude can then be obtained by using the expression for the DDVCS tensor of Eq. (24) in the TCS limit in Eq. (5).  In the following, we will consider the observables for the TCS process on an unpolarized nucleon at small momentum transfer −t q 2 . For these observables, the dominant contribution arises from the structure function H, which is normalized to the proton Dirac FF F D (t). When studying the influence of the radiative corrections on the TCS observables at small values of −t, we will therefore neglect the contribution of the GPDs E,H and E in our study below. For the numerical evaluation in this work, we will use the GPD parametrizations from the VGG model [8,[46][47][48], summarized in Ref. [12] as: (29) The parameterization is based on a double-distribution (DD) ansatz for the (x,ξ)-dependence of the up (down) for valence (sea) quarks respectively, and on a Reggeized ansatz for the t-distribution, which was found to give a global description of existing DVCS data [49,50]. Furthermore, we added an isoscalar so-called Dterm contribution in Eq. (29), which only depends on the two variables x/ξ and t, and which 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].

IV. VACUUM POLARIZATION AT FIRST ORDER
We start our discussion of the first order radiative corrections with vacuum polarization process which is shown to first order in Fig. 6. The photon propagator can be written as: where D µν 0 is the leading order photon propagator 1 The small s-quark GPD contribution is neglected in this work. and Π αβ is the vacuum polarization, which is given at first order in α ≡ e 2 /(4π) by: where m l is the mass of the lepton in the loop. Due to gauge-invariance, q µ Π µν = q ν Π µν = 0, and the vacuum polarization can be decomposed as: The scalar function Π(q 2 ) has an UV divergence. In dimensional regularisation it is given by [37]: where we defined v 2 ≡ 1 − 4m 2 l /q 2 . In the following we consider muons and electrons in the fermion loop.
The UV-divergence (in limit UV → 0+) in Eq. (34) is removed by the renormalization constant Z 3 : In the on-shell scheme this constant 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: and the renormalized photon propagator by: Note that, due to gauge invariance, only the term proportional to g µν contribute to the BH and TCS amplitudes, such that the corrections factorize as: For the evaluation ofΠ(s ll ) we need to perform an analytic continuation of Eq. (37): whereṽ ≡ iv = 4m 2 l /s ll − 1.

V. VERTEX CORRECTION TO THE TCS AMPLITUDE
We next consider the one-loop QED corrections to the process shown in Fig. 7. Figure 7. One-loop QED correction to the vertex.
The corresponding matrix element can be expressed in terms of two form factors, F e D and F e P called Dirac and Pauli electron form factors respectively: To first order in α, the Dirac form factor F e D is divergent, such that a regularisation procedure is needed. In dimensional regularisation it can be expressed as [37]: The Pauli form factor F e P is finite and is given by: where v is defined in the same way as below Eq. (34) by replacing the mass m l → m.
As can be seen from Eq. (43), F e D has an ultraviolet (UV) divergence (in limit UV → 0+), as well as an infrared divergence (in limit IR → 0−). The UV divergence gets removed by the on-shell subtraction scheme, in which the vertex counter term is defined to fix the electron charge e at q 2 = 0. One finds at q 2 = 0 the renormalization constant: Figure 8. Diagrams contributing to the one loop corrections for the process γ * (p1)γ * (p2) → l(p3)l(p4) This leads to the renormalized (on-shell) form factor: For the case when the momentum transfer q 2 becomes timelike, i.e. q 2 > 0, one has to perform an analytic continuation of both form factors.
To calculate the corrections to the dVCS matrix element, we just have to contract Γ µ with the dVCS tensor. Adding the vacuum polarization correction of Eq. (39), this yields the one-loop radiative correction to the TCS amplitude as: whereΓ ν denotes the renormalized vertex, and where we indicated that the momentum transfer which appears in the Compton tensor is given by t. The remaining IR divergence in Eq. (43) will be discussed in Section VII.

VI. ONE-LOOP CORRECTIONS TO THE BH PROCESS
In order to calculate the one-loop diagrams contributing to the BH process, we consider corrections to the process: We show the contributing diagrams in Fig. 8. We use the standard definition of Mandelstam variables: Using p µ 1 , p µ 2 , p µ 3 , g µν and γ µ as building blocks we find 34 independent tensors with two indices to expand the one-loop amplitude. As discussed in Sec. III A, we can reduce the number of amplitudes from 34 to 21 amplitudes, which are gauge invariant and free of any kinematic singularities, such that we can write: with τ µν i the tensor basis introduced in [40], and B i the corresponding invariant amplitudes.
We calculated all 21 contributing amplitudes up to the one-loop order using the same setup as in Ref. [39] in dimensional regularisation. We use QGRAF [52] to generate all contributing diagrams, FORM [53] to implement the Dirac algebra and Reduze 2 [54] to generate IBP identities. The amplitudes are then expressed as a sum over scalar integrals with coefficients which are rational functions in all external scales and the space-time dimension d = 4 − 2 . In total we need 7 master integrals, which are all listed in Ref. [39].
To subtract UV divergences we use the on-shell renormalization scheme. In addition to the vertex counter terms defined in Sec. V, we also need counter terms from the fermion self-energies Σ(p), where p denotes the incoming four momentum. The renormalization condition fixes the pole of Σ(p) at p 2 = m 2 with residue equal to one.
Besides the UV divergences, one also encounters IR divergences in the amplitude. We checked that the matrix element has the correct infrared structure after performing the UV renormalization. The infrared structure of the amplitude can be calculated using the soft-photon approximation, in which the loop momenta l of the photon is assumed to scale soft, i.e. l ∼ λ, where λ is small compared to all external scales. In this approximation the calculation is performed at leading order in λ. The contribution factorizes in terms of the born amplitude: with the three-point function: Note that the same formula applies for the TCS matrix element in the soft-photon limit. Therefore, in this approximation the correction also factorizes on the level of the cross section. It is given by: with the infrared-divergent part: and the finite part: As an additional check we were also able to reproduce the results of Ref. [39] numerically, where the corrections to the unpolarized BH cross section have been calculated.
Using Eq. (50) and adding the contribution from vacuum polarisation, the one-loop matrix element corresponding to the BH process with the incoming photon on-shell can be calculated by contracting with the photon polarization vector and the proton line: where we now can identify t ll = (q − l − ) 2 .

VII. SOFT-PHOTON BREMSSTRAHLUNG
As the virtual corrections have an IR divergence, we also need to account for the soft bremsstrahlung. These correspond to 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, it 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 photon emission.
To evaluate the soft bremsstrahlung contribution to the cross section, one has to integrate over the unobserved soft-photon phase space. It is easiest to perform this integral in a reference frame where the maximum soft-photon energy ∆E s is isotropic. Such a reference frame depends on the specific experimental conditions. It is in general given by the rest frame of the soft photon and the unobserved particle in the process.
In Ref. [38] we considered the case when the final proton is measured and the di-lepton pair remains undetected. Therefore, the integral has to be performed in the rest frame of the di-lepton pair and the soft photon, i.e. l + + l − + k = 0.
In this work, we consider the γp → l − l + p process where l − and l + are observed, while the recoiling nucleon remains unobserved. Therefore the bremsstrahlung contribution to the cross section is evaluated in the rest frame of the unobserved proton and soft photon. Defining a missing momentum p m ≡ p + k, this frame is defined by p m = 0. In such frame the bremsstrahlung contribution due to the soft-photon emission from the l − and l + is given by: where the maximal soft-photon energy in this frame is denoted by ∆E s . We can easily perform the integrations for the first two terms in Eq. (58), which yields the expression: whereβ − ,β + are the lepton velocities, defined in the frame p m = 0:β whereẼ ∓ are the corresponding lepton energies, which can be expressed as: Furthermore in Eq. (59), the integral I −+ is due to the interference between soft-photon emissions from the l − and l + lines. It has been worked out e.g. in Ref. [55] as: with β given in Eq. (7).
In general we can divide the soft-photon contribution of Eq. (59) in an IR divergent piece (δ IR s;r ) and a finite piece (δ s;R ) as: with IR divergent piece given by: and the finite part δ s;r expressed as: In the ultrarelativistic limit for the leptons, i.e. β ≈ 1,β ∓ ≈ 1, which is a very good approximation for the production of an e − e + pair as considered in the following, the above expression simplifies considerably. Using we obtain in this limit for the finite part δ s;r : Note that the IR divergent pieces from virtual and real corrections, ie. Eqs. (54) and (66), exactly cancel on the level of the cross section, thus giving an IR finte result.
For the purpose of estimating the di-lepton forwardbackward asymmetry in the following, it will be convenient to express the lepton angles in the l − l + rest frame as defined in [36], which are denoted by θ l , φ l for polar and azimuthal angles respectively. Eq. (61) then allows to express the di-lepton energiesẼ ∓ in terms of invariants and the angle θ l as: To evaluate the finite soft-photon cross section corrections of Eq. (67) or Eq. (69), we need to express the maximal soft-photon energy ∆E s , defined in the frame p m = 0 in terms of the experimental resolutions using: 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. For the case of detecting the dilepton pair, the lab energies (E ∓ ), the scattering angle between the pair (θ Lab ll ), and the angle between the incoming photon q and the virtual photon q (θ Lab γγ ) are measured. These four parameters are in correspondence with the kinematic quantities s ll , t, θ l , and φ l introduced in the cross section expression of Eq. (6).
In the Lab frame, the missing mass squared can be expressed, neglecting the lepton mass in the kinematics, as: Accounting for the finite resolutions in the Lab kinematical quantities, and adding their contributions quadratically, we can express the maximal soft-photon energy ∆E s as: As is evident from Eq. (73), the evaluation of the softphoton radiative correction depends on the specific experimental resolutions. For the purpose of illustrating the size of these corrections, we will provide predictions below where the soft-photon cut-off energy ∆E s is chosen to be in the one to few percent range of the beam energy, as a realistic value. For our predictions in the low-energy region (E γ 0.36 GeV), we will show the corrections for ∆E s = 0.01 GeV, whereas for the highenergy region (E γ 6.8 GeV), we will show the results for ∆E s = 0.05 GeV.

A. Observables
We use our setup to study the radiative corrections to the γp → e − e + p process in both low-and high-energy regimes. For both kinematical situations we study the effect of these corrections on the cross section, on the forward-backward asymmetry A F B , as well as on the beam helicity 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 (in the l − l + rest frame), and where dσ ± in A stand for the polarized cross sections for circular photon polarization ±1 respectively. Because of the opposite symmetry of the BH and TCS amplitudes (Fig. 1) under charge conjugation (odd versus even number of photon couplings to the lepton charge), the asymmetry A F B , which interchanges the kinematics for l − and l + , allows for a direct assessment of the interference term between the BH and the real part of the TCS amplitude. The BH and TCS processes separately yield a zero asymmetry. While A F B is proportional to the real part of the BH-TCS interference, A is proportional to the imaginary part of this interference. Note however that the complex TCS amplitude by itself also yields a contribution to A , which is usually very small unless the real part of the TCS amplitude becomes comparable in size to the BH amplitude. The observables A F B and A are thus complementary in accessing the complex TCS amplitude experimentally and in testing theoretical models.
In order to calculate the TCS observables we implemented all amplitudes in a C++ code in which the interference of different Feynman amplitudes can be evaluated numerically.

B. Results for TCS observables in the ∆(1232) region
We firstly study the importance of the radiative corrections on the γp → e − e + p observables in the low-energy kinematical region. The γp → e − e + p process was studied in Ref. [36] in the ∆(1232) resonance region at small values of s ll and −t. In this limit, the TCS amplitude approaches the forward real Compton scattering amplitude, for which a full dispersive calculation based on empirical structure functions exists [56]. It was found in [36] that around W = 1.25 GeV the full TCS cross section integrated over the di-lepton angles is reproduced by a Born + ∆-pole model for the TCS amplitude, as discussed in Sec. III A, within an accuracy of 5% or better. Therefore, we consider the Born + ∆(1232)-pole model to be realistic enough as a model for the TCS amplitude around the ∆(1232)-pole in the near forward direction in order to study the effect of the QED radiative corrections on this reaction. In the following, we will show the effect of the radiative corrections on the γp → e − e + p cross section as well as on the two asymmetries A F B and A for a c.m. energy of W = 1.25 GeV.
In Fig. 9 we show the s ll -dependence of the unpolarized γp → e − e + p observables for in-plane di-lepton kinematics (θ l = 0 • and φ l = 0 • ) in the ∆(1232) region and for a small value of −t. One notices that the TCS amplitude increases the cross section in this kinematical region by 50% or more, and that the forward-backward asymmetry A F B , which depends linearly on the TCS amplitude reaches values larger than 30%. We see that the first-order QED radiative corrections yield a 5 -8% correction on the γp → e − e + p cross section in this kinematical range. However, we also notice that in comparing the complementary forward (δ F for θ l = φ l = 0 • ) and backward (δ B for θ l = φ l = 180 • ) kinematics, these correction factors are nearly the same, such that the asymmetry A F B is to good approximation unaffected by the first-order QED radiative corrections. This makes A F B an ideal observable to extract the real part of the TCS amplitude.
Until recently, the only proof-of-principle experiment of the forward-backward asymmetry of the γp → e − e + p process was conducted by Alvensleben et al. at DESY [57]. It was performed at W 2.2 GeV in nearforward kinematics, and aimed at an experimental verification of the Kramers-Kronig dispersion relation for the forward Compton amplitude. The Kramers-Kronig relation encompasses a fundamental connection between the photon absorption and scattering based on analyticity and unitarity, and allows to evaluate the real part of the forward Compton scattering off protons through a dispersive integral over its imaginary part, which is evaluated using the empirical knowledge of the total photoabsorption cross sections. It was shown in a recent re-analysis of the Kramers-Kronig relation for forward Compton scattering [56] that the present database of the unpolarized photoabsorption cross section is not entirely consistent, and shows the largest discrepancies (in the 5 -10% range) in the region of the ∆(1232) peak and in its higher-energy tail region. It was furthermore shown in [56] that these discrepancies in the world database also directly limit the obtained precision of the extracted Baldin sum rule  value for the sum of proton electric and magnetic polarizabilities, α E1 + β M 1 . Being nearly unaffected by the first-order QED radiative corrections, our results show that the forward-backward asymmetry in the resonance region provides a direct measurement of the real part of the near-forward Compton amplitude and thus has the potential to settle the existing discrepancies in the Compton database.
In Fig. 10 we show the corresponding di-lepton mass dependence of the γp → e − e + p observables for a circularly polarized photon beam for out-of-plane di-lepton kinematics (θ l = φ l = 90 • ). As before we see that the first-order QED radiative corrections yield a 7 -9% correction on the γp → e − e + p cross section in this kinematical range. The small cusp in the curves of δ ± around s ll ≈ 0.044 GeV 2 is due to the muon threshold in the vacuum polarization. This gives a contribution with different sign for the cross sections dσ ± . For the corresponding photon helicity asymmetry A , which is directly proportional to the imaginary part of the BH-TCS interference, the first-order QED radiative corrections nearly drop out. For kinematics close to the ∆(1232)-pole position, where the imaginary part of the TCS amplitude is maximal, we find a large value of the asymmetry, varying between −65% and −45%.
In Fig. 11 we show the dependence of both asymmetries in the ∆(1232) region on the lepton angle θ l for fixed out-of-plane angle φ l (φ l = 0 • for A F B and φ l = 90 • for A ) and different values of t and s ll . We see that increasing the values of −t and s ll results in larger asymmetries due to the larger interference terms. We only show the radiatively corrected results in Fig. 11, as the difference between the tree-level result and the result including first-order QED radiative corrections is at the few per mille level on these asymmetries.

C. Results for high-energy TCS observables
In this section we study the effect of the first-order QED radiative corrections on the cross section as well as on the forward-backward and beam helicity asymmetries of the γp → e − e + p process in the high-energy, forward scattering regime in which the di-lepton pair is produced with a large virtuality. In this regime, a QCD factorization theorem allows to model the TCS amplitude in terms of GPDs as described in Sec. III B. We present here results in the kinematics of a recent CLAS12@JLab experiment [17]. The latter is the first TCS experiment where data have been reported [18,19] in the kinematical regime corresponding with an average c.m. energy of W = 3.69 GeV and a large di-lepton invariant mass squared of s ll = 3.24 GeV 2 .
In Fig. 12 we show the t-dependence of the unpolarized γp → e − e + p observables for in-plane di-lepton kinematics (θ l = 65 • and φ l = 0 • ), corresponding with the CLAS12@JLab experimental conditions. One notices that the cross section and the forward-backward asym- metry show a sizeable sensitivity on the D-term contribution to the GPD parameterization, which contributes to the real part of the TCS amplitude. Comparing the GPD double distribution parameterization (shown by blue dashed curves in upper and lower panels in Fig. 12), which was used in a previous global analysis of DVCS data in the valence region [49,50], we see that adding the dispersive estimate of Ref. [51] for the D-term contribution (dotted red curves) gives a sizeable contribution to the cross section. For the forward-backward asymmetry, which is directly proportional to the real part of the BH-TCS interference, it leads to a shift in the asymmetry by 10 -20%. We see that the first-order QED radiative corrections yield a nearly 20% correction on the γp → e − e + p cross section, which are important to account for in the extraction of the Compton form factors. However, we also notice that in comparing the complementary forward and backward kinematics, these correction factors are nearly the same, such that the asymmetry A F B is to good approximation unaffected by the first-order QED radiative corrections. This makes A F B a gold-plated observable to extract the real part of the Compton form factor, and test its sensitivity to the D-term contribution in GPD parameterizations.
In Fig. 13 we show the t-dependence of the γp → e − e + p observables for a circularly polarized photon beam for out-of-plane di-lepton kinematics (θ l = φ l = 90 • ), corresponding with the CLAS12@JLab experimental conditions. We notice that in out-of-plane kinematics the sensitivity to the D-term in the GPD parameterization is very small, while the radiative correction on the polarized cross sections is also around 20%. In the corresponding photon beam helicity asymmetry A , which is directly proportional to the imaginary part of the BH-TCS interference, the first-order QED radiative corrections again nearly drop out. As the D-term is purely real, its effect on A is very small, and contributes only through the squared modulus of the complex TCS amplitude. The imaginary part of the latter is rather well constrained by the precise electron beam spin asymmetry data for the corresponding DVCS process [49,50]. As the non-perturbative part in both processes is the same, a direct comparison between the GPD predictions for DVCS and TCS beam helicity asymmetries at the same values of −t and ξ, provides a stringent test of the applicability of the underlying QCD factorization theorem at these kinematics. Furthermore, as can be seen from the lower panel in Fig. 10, such test is basically unaffected by the first-order QED radiative corrections.

IX. CONCLUSIONS
In this paper we presented the first-order QED corrections on the lepton side contributing to the timelike Compton scattering process on a proton, γp → l − l + p. This reaction contains contributions from both the Bethe-Heitler amplitude, for which the dependence on the proton structure is parameterized through its spacelike elastic form factors, and from the TCS amplitude with intial real photon and final timelike photon. We calculated the first-order radiative corrections on the level of the amplitude, such that the individual subprocesses can be easily incorporated into different calculations. We studied the effect of radiative corrections in two different energy regimes. In the ∆(1232) resonance region the TCS amplitude was modeled as the sum of Born and ∆(1232)-pole contributions, which was found to give a very good description of the near-forward kinematical regime. In the high-energy, near-forward regime we calculated the TCS amplitude in terms of GPDs, in kinematics of a recent CLAS12@JLab experiment.
In the kinematics near the ∆(1232)-resonance position, we found that the first-order QED radiative corrections on the γp → e − e + p cross section are in the 5 -10% range, while in the high-energy kinematics of the CLAS12@JLab experiment, they are in the 20% range. Their inclusion is thus important to extract the low-energy constants or Compton form factors from the γp → e − e + p process.  Figure 12. Upper panel: t-dependence of the γp → e − e + p unpolarized cross section for in-plane kinematics of the dilepton pair corresponding with the CLAS12@JLab experiment [18,19]. We show our results for two GPD models and including the first-order QED radiative corrections, for a soft-photon energy of ∆Es = 0.05 GeV. The middle panel displays the radiative correction factor for the forward (δF for θ l = 65 • , φ l = 0 • ) and backward (δB for θ l = 115 • , φ l = 180 • ) cross sections. The lower panel shows the sensitivity of the forward-backward asymmetry AF B on the D-term contribution to the GPD, and the negligible radiative correction on this observable (solid and dotted curves nearly overlapping).  Figure 13. Upper panel: t-dependence of the γp → e − e + p unpolarized cross section for out-of-plane kinematics of the di-lepton pair corresponding with the CLAS12@JLab experiment [18,19]. We show our results for two GPD models and including the first-order QED radiative corrections, for a soft-photon energy of ∆Es = 0.05 GeV. The middle panel shows the cross section correction factors δ± for circular photon polarization ±1. The lower panel shows the sensitivity of the photon beam helicity asymmetry A on the GPD parameterization, and the negligible radiative correction on this observable (solid and dotted curves nearly overlapping).
Besides the corrections on the unpolarized cross section, we also studied in both kinematical regimes the effect of the first-order radiative corrections on the forwardbackward asymmetry A F B , obtained by interchanging the kinematics of the produced di-leptons, as well as on the photon beam helicity asymmetry A . While the asymmetry A F B accesses the real part of the BH-TCS interference, the asymmetry A accesses the imaginary part of the BH-TCS interference in the regime where the BH amplitude dominates. The asymmetries A F B and A are thus complementary in accessing experimentally the complex TCS amplitude.
We found that in both kinematical regimes the radiative corrections on these asymmetries are at the few per mille level only, although the radiative corrections on the γp → e − e + p cross sections are in the 10 -20% range. The reason for the near cancellation of radiative corrections is that the corrections on the forward and backward cross sections as well as on both beam helicity cross sections are of the same size, and thus nearly cancel out in the corresponding ratios. This makes the A F B and A gold-plated observables to extract the real and imaginary parts of the TCS amplitude respectively.
In the ∆(1232) region, we find a value around +30% for the forward-backward asymmetry, which provides a good opportunity for a direct measurement of the real part of the near-forward Compton amplitude. This will also allow for a comparison with existing dispersive extractions, in which the real part of the forward Compton amplitude has been obtained as a dispersive integral over its imaginary part. At present the latter extraction is limited in precision due to discrepancies in the world database for photoabsorption on a proton, especially in the ∆(1232) region, which also directly limits the obtained precision of the extracted Baldin sum rule value for the sum of proton electric and magnetic polarizabilities, α E1 + β M 1 .
In the kinematical regime of the CLAS12@JLab experiment, the measurement of both asymmetries A F B and A will allow to extract the Compton form factors from the TCS amplitude and compare them to the corresponding DVCS process with a spacelike initial photon and final real photon. As the non-perturbative part in both processes is the same, a direct comparison between the GPD predictions for DVCS and TCS observables, provides a stringent test of the applicability of the underlying QCD factorization theorem at these kinematics. In particular, the forward-backward asymmetry is a very sensitive observable to extract the D-term contribution to the GPD parameterization, as it only contributes to the real part of the TCS amplitude. We found that a dispersive estimate of the D-term contribution leads to a shift in A F B by 10 -20%. The latter sensitivity will be interesting to test by the CLAS12@JLab TCS experiment, and is complementary to the information one extracts from the DVCS process, when comparing cross sections for both electron and positron beams.
As a next step we plan to generalize our radiatively corrected calculations to the double virtual Compton scattering, in which the incoming photon has a nonzero spacelike virtuality while the outgoing photon has a non-zero timelike virtuality, accessed through the e − p → e − pl − l + reaction. In the ∆(1232) region, it was shown that this process allows for an empirical determination of the remaining unknown low-energy structure constant entering the hadronic correction to the muonic hydrogen Lamb shift [36]. In the high-energy near-forward region, double deeply virtual Compton scattering allows to extend the DVCS beam spin asymmetry measurements of GPDs into the ERBL domain [20,21]. Since we calculated all sub-parts of the amplitudes for off-shell particles, it will be straightforward to generalize our setup to the double virtual Compton observables, accessed through the e − p → e − pl − l + reaction.