Two-photon exchange correction to muon–proton elastic scattering at low momentum transfer

We evaluate the two-photon exchange (TPE) correction to the muon–proton elastic scattering at small momentum transfer. Besides the elastic (nucleon) intermediate state contribution, which is calculated exactly, we account for the inelastic intermediate states by expressing the TPE process approximately through the forward doubly virtual Compton scattering. The input in our evaluation is given by the unpolarized proton structure functions and by one subtraction function. For the latter, we provide an explicit evaluation based on a Regge fit of high-energy proton structure function data. It is found that, for the kinematics of the forthcoming muon–proton elastic scattering data of the MUSE experiment, the elastic TPE contribution dominates, and the size of the inelastic TPE contributions is within the anticipated error of the forthcoming data.


Introduction
The present measurements of the proton charge radius from muonic hydrogen spectroscopy [1,2] differ by a puzzling 7σ from the radius value extracted from the hydrogen spectroscopy [3] and the elastic electron-proton scattering data [4]. This huge discrepancy, which has become known as the "proton radius puzzle", see Ref. [5] for a recent review, has led to intense theoretical and experimental activity in recent years. So far it has defied an explanation which could bring all three experimental techniques in agreement with each other. To shed further light on this puzzle, several new experiments involving muons are being planned. Their aim is to test the lepton universality in the interaction of a lepton with a proton. One can compare the elastic scattering of electrons and muons on the proton target and measure the proton charge radius in the muon-proton elastic scattering in a similar way as was done in the electron-proton elastic scattering [4,6]. Such an elastic scattering experiment is presently being planned by the MUSE collaboration [7]. Complementarily, one can also compare the electron and muon pair photoproduction on the proton as proposed in [8]. These experiments should be performed at the 1 % level or better of experimental accuracy in order to have an impact on the observed discrepancy in the proton charge radius. Such a level of precision therefore calls for studies of the higher order corrections in such processes, as the corrections to cross sections suppressed by one power in the fine-structure constant α = e 2 /(4π) ≈ 1/137 are also in the 1 % or few % range. In particular it requires studies of the two-photon exchange (TPE) correction to the unpolarized lepton-proton elastic scattering for the case when the mass of the lepton cannot be neglected relative to its momentum. In a previ-ous work [9], we have performed an estimate of the leading TPE contribution from the proton intermediate state, and provided estimates for the MUSE experiment. In this work we account for the inelastic intermediate states, i.e. all possible intermediate states in the TPE box graphs beyond the proton state. As our aim is an estimate of such corrections for the MUSE experiment, which corresponds with very low momentum transfers, we will estimate the inelastic TPE corrections through the near forward doubly virtual Compton scattering process [10]. The essential hadronic information is contained in the unpolarized proton structure functions and in one subtraction function. We clarify the TPE correction coming from the subtraction function, and we provide an empirical determination of the subtraction function based on the high-energy behavior of the forward Compton amplitude T 1 .
The plan of the present paper is as follows. We review the elastic TPE contribution to the unpolarized leptonproton scattering in Sect. 2, and derive a low momentum transfer expansion accounting for all terms due to the non-zero lepton mass. We subsequently discuss the forward unpolarized doubly virtual Compton scattering process which will serve as our starting point in the determination of the inelastic TPE corrections in Sect. 3. We evaluate the TPE correction due to the subtraction function in the forward Compton amplitude T 1 , and we provide an empirical determination of the subtraction function from data in Sect. 4. We provide the expressions for the inelastic TPE correction coming from the unpolarized proton structure functions and present the results of our numerical evaluation for the muon-proton elastic scattering in Sect. 5. Our conclusions are given in Sect. 6.

Elastic TPE contribution
In this work we study the TPE correction at low momentum transfer to unpolarized elastic scattering of a charged lepton with mass m and initial (final) momentum k (k ) on a proton with mass M and initial (final) momentum p ( p ); see Fig. 1 for the notations of kinematics and helicities. In this work, we consider the region of small squared momentum transfer where E is the lepton beam energy (in the lab frame) and We define the TPE correction δ 2γ from the difference between the expected cross section σ exp and the cross section in the one-photon exchange (OPE) approximation σ 1γ : corresponding with first order corrections in the fine structure constant α = e 2 / (4π ), with e the unit of electric charge. In the low Q 2 region, the dominant contribution to the TPE graph of Fig. 1 results from the proton intermediate state (elastic contribution). We have fully calculated this contribution in a previous work [9]. In this section, we start by studying the quality of some approximate expressions for the elastic contribution in the low Q 2 limit, for which analytical expressions can be provided.
The first TPE estimate is due to Feshbach and McKinley [11], who calculated the TPE contribution, corresponding with Coulomb photon couplings to the static proton (i.e. two γ 0 vertices). This so-called Feshbach term contribution to δ 2γ in Eq. (1) is denoted by δ F and can be expressed through the scattering angle in the laboratory frame θ and the lepton velocity v as As a next step, one may consider the TPE correction in the scattering of two point-like Dirac particles (corresponding with two γ μ couplings). In Appendix A we provide some analytical expressions of this contribution in the limit of small Q 2 M 2 , M 2 (E 2 − m 2 )/s, with the center-of-mass frame squared energy s = M 2 + m 2 + 2M E, both for the cases when Q 2 m 2 ; see Eqs. (A6)-(A8), and when Q 2 and m 2 are of similar size, see Eq. (A9).
We show in Fig. 2 (left panel) the comparison between the Feshbach term, the TPE contribution for point-like Dirac particles, and the TPE for a point-like proton, with inclusion of the magnetic moment contribution. It is seen that the Feshbach correction of Eq. (2) with account of the recoil correction factor (1 + m/M) describes the result for point-like Dirac particles quite well in the kinematics of the MUSE experiment.
We also show in Fig. 2 (right panel) the effect of the proton FFs, according to the full numerical calculation of Ref. [9]. In the low Q 2 kinematics of the MUSE experiment, the inclusion of the FFs provides a reduction of the TPE by around 40 % at Q 2 ≈ 0.025 GeV 2 , consequently one should use the

Forward unpolarized doubly virtual Compton scattering tensor
The TPE contribution, δ 2γ , for the muon-proton scattering process is in general given by the interference of the one-photon exchange amplitude (T OPE ) and the two-photon exchange amplitude (T TPE ) as The OPE expression in the denominator of Eq. (3) is given by where G E (G M ) are the proton electric (magnetic) form factors, respectively, and with kinematical quantities τ and ε 0 defined as in Eq. (A3). The interference between the OPE and TPE amplitudes in the numerator of Eq. (3) can be expressed as with all momenta defined as in Fig. 1, whereq is the loopmomentum over which one integrates, and where α denotes the on-shell proton electromagnetic vertex: with F D (F P ) the Dirac (Pauli) form factors of the proton, respectively. Furthermore in Eq. (5), M μν denotes the proton doubly virtual Compton scattering tensor.
The main aim of the present work is to quantitatively estimate the inelastic TPE contribution in the low Q 2 region, corresponding with the MUSE kinematics. For this purpose, we will approximate the hadronic tensor M μν in Eq. (5) by the forward doubly virtual scattering (VVCS) tensor. The unpolarized forward VVCS process γ * (q) + p(P) → γ * (q) + p(P) is described by two invariant amplitudes T 1 and T 2 , which are defined in this work through the tensor decomposition: with the photon energyν = (P ·q) /M and the squared photon virtualityQ 2 ≡ −q 2 . The absorptive parts of the amplitudes T 1 and T 2 are related to the proton structure functions F 1 and F 2 by In this work, we will approximate the unpolarized tensor M μν entering Eq. (5) for the process γ * (q 1 =q + q/2) + p(P − q/2) → γ * (q 2 =q − q/2) + p(P + q/2) in the low momentum transfer limit q → 0, i.e. q 1 ≈ q 2 , by [10]: By using the electromagnetic gauge invariance of the lepton tensor, the hadronic tensor of Eq. (9) can be expressed equivalently as which is the tensor form which we will use in our evaluations of δ 2γ . The real part of the amplitude T 1 can be expressed through a subtracted dispersion relation (DR) as integral over the invariant mass W 2 of the intermediate hadronic state as with the pion-proton inelastic threshold: W 2 thr = (M + m π ) 2 ≈ 1.15 GeV 2 , where m π denotes the pion mass, and where T subt 1 (0,Q 2 ) is the subtraction function atν = 0. The real part of the amplitude T 2 can be obtained from an unsubtracted DR: In Eqs. (12) and (13) . (15) Note that in the derivation of a DR as given e.g. in Eq. (12) the elastic (nucleon pole) term contribution, given by only the first term of Eq. (14), correctly appears. This pole contribution differs from the Born term by As this is an energy (ν) independent function, we have absorbed it in the definition of T subt The advantage of expressing the amplitude T 1 w.r.t. to its Born contribution, results from the fact that the non-Born amplitude in Eq. (17) starts atQ 2 and is usually parametrized in terms of polarizabilities, i.e. the function β(Q 2 ) atQ 2 = 0 is given by the magnetic polarizability β M : β(0) = β M [13,21]. In order to evaluate the inelastic TPE contribution using the forward non-Born VVCS amplitudes, we will need the information on the proton structure functions F 1 and F 2 as well as to specify the subtraction function in Eq. (17), which is parametrized through the function β(Q 2 ). We will evaluate the subtraction function contribution in Sect. 4, and the dispersive contribution due to the structure functions F 1 and F 2 in Sect. 5.

Subtraction function contribution to TPE correction
In the present section we will discuss the TPE correction due to the subtraction function T subt 1 (0,Q 2 ). For this purpose, we will compare three different estimates for β(Q 2 ) defined through Eq. (17). At low Q 2 we will use existing estimates from heavy-baryon and baryon chiral perturbation theory. Furthermore, we will provide an empirical determination of β(Q 2 ), based on the high-energy behavior of the Compton amplitude.

Heavy-baryon ChPT subtraction function
First of all, we show the fit of Ref. [13] obtained by matching the heavy-baryon chiral perturbation theory (HBChPT) result to a dipole behavior: (18) with the value of the magnetic polarizability β M = (2.5 ± 0.4) × 10 −4 fm 3 taken from PDG [14]. For the purpose of showing error bands in our numerical estimates, we choose the lower and upper edges of such bands to correspond with the values: = 530 MeV, β M = 2.1 × 10 −4 fm 3 and = 842 MeV, β M = 2.9 × 10 −4 fm 3 , respectively. The resulting bands for T subt 1 (0,Q 2 ) are shown in Fig. 3, and correspondingly for β(Q 2 ) in Fig. 4 (blue bands).

Empirical determination of the subtraction function
In this section, we discuss an empirical estimate of the function β(Q 2 ) at non-zeroQ 2 from experimental information on inelastic electron-proton scattering. Following the idea of Refs. [17,18], the subtraction function can be obtained from an unsubtracted dispersion relation for the amplitude , where T R 1 denotes a Regge amplitude which is chosen such as to match the high-energy behavior of the amplitude T 1 , i.e. T 1 − T R 1 → 0 forν → ∞. 1 The function T R 1 is chosen as a sum over the leading Regge trajectories: with the intercept α 0 > 0,ν 0 is a reference hadronic scale which is used as a free parameter and γ α 0 (Q 2 ) are the Regge residues. Using Eq. (8), the imaginary part of T R 1 yields the corresponding Regge structure: The Regge residues γ α 0 (Q 2 ) can be obtained by performing a fit to inclusive electroproduction data on a proton. In our work we use the Donnachie-Landshoff (DL) high-energy fit [19] to obtain the proton structure function F 1 as where the values of the Regge intercepts α 0 and the residue functions γ α 0 (Q 2 ) are detailed in Appendix B. By comparing Eqs. (21) and (22) we notice that the second term in Eq. (21) is chosen such that for the Regge trajectory with 1 < α 0 < 2 ("Pomeron"): whereas for the Regge trajectory with 0 < α 0 < 1 (Reggeon): This ensures that in all cases the quantity Consequently, one can write down an unsubtracted dispersion relation for T 1 − T R 1 at fixedQ 2 as: Using Eqs. (8), (16), and (17), this yields an expression for T subt 1 (0,Q 2 ) which, expressed in terms of W 2 ≡ 2Mν + M 2 −Q 2 , is given by where the lower integration limit in Eq. (26) s thr is given by corresponding with a branch cut of F 1 starting at W 2 thr and a branch cut of F R 1 starting at s 0 . Eq. (26) allows one to quantitatively estimate the subtraction function given the structure function F 1 , the Regge fit determining F R 1 of the form of Eq. (21), as well as the corresponding value of T R 1 (0,Q 2 ), which follows from Eq. (20) as and is also fully determined by the Regge fit. In our numerical evaluation of Eq. (26), we describe the proton structure function F 1 in the resonance region by the fit performed by Christy and Bosted (BC) [20]. This fit is valid in the following region of kinematic variables: 0 < Q 2 < 8 GeV 2 , and W 2 < 9.61 GeV 2 ≈ 10 GeV 2 . For the dispersion integral in Eq. (26) we connect the BC fit with the DL high-energy fit starting from W 2 = 10 GeV 2 . The latter fit is described in Appendix B. The resulting proton structure function F 1 is shown in Fig. 5 as it enters the integral of Eq. (26). We add a 3 % error band to the BC fit [20] and use the same error estimate for all Regge pole residues. We notice that at low values ofQ 2 , both fits either overlap or are very close around the matching point W 2 ≈ 10 GeV 2 .
With increasing values ofQ 2 there is a slight mismatch in both fits around W 2 = 10 GeV 2 , which is due to the fact that the BC fit has not accounted for the HERA high-energy data, and the DL fit has not accounted for the lower W data. Even though a combined fit of all data would be very worthwhile, or a smooth interpolating procedure between the BC and DL fits could easily be performed, for our purpose we will only need data at lower value ofQ 2 up to about 1 GeV 2 . For this purpose, we can just split the W 2 integral entering Eq. (26) in a region W 2 < 10 GeV 2 where we will use the BC fit and a region W 2 > 10 GeV 2 where we will use the DL fit.
In Fig. 6, we demonstrate explicitly the vanishing highenergy behavior of the quantity F 1 -F R 1 , which is the necessary condition for the unsubtracted DR of Eq. (26) to hold.
We furthermore provide another consistency check of our numerical implementation. As the Regge function T R 1 of Eq. (20) has an arbitrary scaleν 0 (or equivalently s 0 ), the total result should not depend on the specific choice of this parameter. We demonstrate this in Fig. 7, where we illustrate how the s 0 dependence of the individual contributions in Eq. (26) adds up to yield the total result which is independent of s 0 .
In Fig. 3 we present the empirically extracted subtraction function T subt 1 (0,Q 2 ) of Eq. (26) and compare it with the subtraction functions of Refs. [13,15]. The subtraction function T subt 1 (0,Q 2 ) should vanish linearly whenQ 2 → 0 according to Eq. (17). This general property therefore provides a quality check on the accuracy of an empirical determination as described above. One notices from Fig. 3 that the value of T subt 1 atQ 2 = 0 is compatible with zero within 1-1.5 σ . We would like to notice, however, that at present such empirical determination can unfortunately only give the correct order of magnitude of T subt 1 (0,Q 2 ). This is partly due to the non-perfect match between the proton F 1 fits for the resonance region and the large W region, as we have shown in Fig. 5. Despite this caveat, it seems, however, that with increasingQ 2 , T subt 1 (0,Q 2 ) changes sign in the range somewhere between 0.1 and 0.4 GeV 2 , which may be an indication of the range up to which the ChPT-based results can be used. To provide a more accurate determination of the functional dependence of T subt 1 (0,Q 2 ), a combined fit of all proton F 1 structure function data over the whole range of W , incorporating the Regge behavior at large W would be desirable. At intermediate values of Q 2 , below and around 1 GeV 2 , this will also require one to have more accurate data in the intermediate W range between 3-10 GeV. In the lower end of this range, such data can be provided by the JLab 12 GeV facility.
Using our empirical determination of T subt 1 (0,Q 2 ), we can extract β(Q 2 ) dividing T subt 1 (0,Q 2 ) byQ 2 according to Eq. (17). For the purpose of combining our empirical estimate of T subt 1 (0,Q 2 ) with the empirical value of β(0) as determined from RCS, we use the central curve in the empirically determined error band of T subt 1 (0,Q 2 ) (green band in Fig. 3) to extract β(Q 2 ) in the rangeQ 2 > 0.12 GeV 2 , and extrapolate it by a linear function to the PDG value of β M at Q 2 = 0. The resulting curve is displayed in Fig. 4. We will use the latter curve in the following to provide an empirical estimate for the subtraction function contribution to the TPE correction for the muon-proton elastic scattering at small momentum transfer.

TPE correction from the subtraction function
Using Eqs. (3)-(5), we can now estimate the TPE correction due to the T subt 1 (0,Q 2 ) contribution to the first term in the hadronic tensor of Eq. (7). Performing the traces in Eq. (5) explicitly, the subtraction function results in the following TPE correction in the region of low momentum transfers: where the lepton (photon) propagators ± K ( ± Q ) are defined as in Eq. (A5). The second term within the curly brackets of Eq. (29) can be simplified to yield the expression making explicit the overall proportionality of δ subt 2γ to the squared lepton mass m 2 .
The integration in Eq. (30) is performed through a Wick rotation, as detailed in Appendix C, and the resulting TPE correction is given by in terms of the dimensionless variable x = 4Q 2 /Q 2 with the weighting function f (x, a): At low momentum transfers the result of Eq. (31) starts from a term proportional to Q 2 . We show the x (orQ 2 ) dependence of the weighting function of Eq. (32) in Fig. 8. The TPE contribution due to the subtraction function also provides a correction to the 2S-2P muonic hydrogen Lamb shift, which is the largest hadronic uncertainty in this precise quantity [1,2]. Using the ChPT-based results for β(Q 2 ) as input, this TPE correction was estimated in Refs. [13,15] and found to be too small to resolve the proton radius puzzle. The correction from the above discussed empirically determined subtraction function to the 2S energy level in muonic hydrogen, after integration up toQ 2 = 1 GeV 2 , yields: 2 which is in fair agreement with the estimate of Birse et al. [13], though slightly smaller: E subt 2S ≈ 4.2 ± 1.0 μeV. Our result of Eq. (34) is also within errors of the analogous evaluation of Ref. [22], where the authors assumed the existence of a J = 0 fixed pole. It was speculated in Ref. [23] that to explain the proton radius puzzle would require a huge enhancement of β(Q 2 ) at largeQ 2 . In order to account for the experimentally observed discrepancy in E 2S of around 310 μeV [5], it would require an around two orders of magnitude larger TPE correction than the naturally expected result from the ChPT estimates. For this purpose, an ad hoc subtraction function, proposed to be added as an extra contribution on top of the ChPT-based subtraction functions discussed above, was conjectured in Ref. [23] with the following functional form: In such a scenario, the largeQ 2 region would also dominate the TPE correction to the muon-proton elastic scattering, and the integral of Eq. (31) would be approximated by where the last step gives the approximate expression in the limit Q 2 M 2 , M E, E 2 . This approximation corresponds in magnitude with the result of Ref. [23] for μ − p scattering, however, it differs by an overall sign.
In Fig. 9 we compare the TPE correction to elastic muonproton scattering (for MUSE kinematics) due to the above discussed ChPT as well as empirically determined subtraction functions. To estimate the size of uncertainties of the BChPT result [15], we plot a band corresponding with a variation of the upper integration limit in Eq. (31) betweeñ Q 2 = 0.9−5 GeV 2 . We notice that the HBChPT and BChPT results are in agreement within their uncertainties.
The TPE correction due to the empirically extracted subtraction function is also shown on Fig. 9, giving a similar 2 Note that for the total TPE correction to the muonic hydrogen 2S level one needs to add to the subtraction function contribution also the dispersive contribution, which was evaluated based on data in Ref. [21]. Blue band result for the HBChPT-based subtraction function [13]. Pink band result for the BChPT-based subtraction function [15].  . 10 The dependence of the integral of Eq. (31) on the upper integration limitQ 2 max for three different estimates of the subtraction function β(Q 2 ) as described in the text though slightly smaller result. This can be understood as the empirically determined β(Q 2 ) changes sign as a function ofQ 2 . The region ofQ 2 contributing to the above result is shown in Fig. 10. One sees that the TPE integral has largely converged for an upper integration limit value of around Q 2 max ∼ 1 GeV 2 . In Fig. 9, we furthermore also show the TPE correction to elastic muon-proton scattering resulting from the subtraction function conjectured in Ref. [23] to explain the proton radius puzzle through enhancing the TPE corrections by nearly two orders of magnitude. Even though the weighting functions entering the TPE corrections in the muonic hydrogen Lamb shift and the elastic muon-proton scattering are different, one notices from Fig. 9 that the subtraction function of Ref. [23] also yields a nearly two order of magnitude larger TPE correction for the elastic muon-proton scattering. To put this in perspective, we also display in Fig. 9 the model independent estimate of the elastic TPE contribution, which has to be added on top of the inelastic TPE contribution, and which is due to the Feshbach term of Eq. (2) corrected by the recoil factor (1 + m/M). One notices that the use of such large subtraction function would yield an inelastic TPE correction to elastic muon-proton scattering which in magnitude already would exceed the elastic Feshbach contribution around Q 2 = 0.02 GeV 2 , and would increase further with increasing Q 2 .

Inelastic contribution to TPE correction
Besides the subtraction function contribution, the inelastic TPE correction to elastic muon-proton scattering includes the contribution of the DR integrals in Eqs. (12) and (13). Using Eqs. (3) and (4) and working out the traces in Eq. (5), the corresponding contributions from the unpolarized proton structure functions F 1 and F 2 to δ 2γ are given by with definitions introduced in Eqs. (A2)-(A5) and the following notation: Our numerical studies of the inelastic TPE contribution indicate that, in the limit Q 2 m 2 , M 2 , M E, the momentum transfer expansion starts with a Q 2 term and contains no Q 2 ln Q 2 type of non-analyticity. This is unlike the elastic electron-proton scattering case, where in the limit m 2 Q 2 M 2 , M E a non-analytic behavior of the type Q 2 ln Q 2 is present at low Q 2 [10,12].
To provide numerical estimates of the inelastic TPE contribution to elastic muon-proton scattering due to the dispersion integrals, we express the corresponding integrals of Eqs. (37) and (38) in the form In Figs. 11, 12 we compare the W dependence of the integrand f (W ) in Eq. (40) for the μ − p and e − p elastic scatter-

Fig. 12
Same as Fig. 11, but for the lepton momentum k = 115 MeV ing processes. As input for the proton structure functions F 1 and F 2 , we use the fit performed by Christy and Bosted [20]. We find that theQ 2 integrations are well saturated when performed up toQ 2 = 8 GeV 2 , which is the largest value covered by the BC fit. As a test, we extended the BC fit beyond its fit region and found that the relative contribution from the region 8 GeV 2 <Q 2 < 12 GeV 2 is smaller than 0.015 %. Figures 11, 12 show results in different kinematics corresponding with the MUSE experiment. The TPE corrections to e − p are sizeably larger than for the μ − p case at low Q 2 . With increasing Q 2 , the μ − p TPE corrections increase, as is evident from the result at lower beam momentum in Fig. 12, where at Q 2 = 0.03 GeV 2 both corrections reach similar sizes. We furthermore notice in Fig. 11 that the integrand for the elastic e − p scattering displays a narrow peak corresponding with the quasi-real photon singularity (for both photons), see Ref. [10], which is absent for the μ − p case.
To estimate the inelastic TPE correction to elastic leptonproton scattering, we find that the W integration in Eq. (40) is well saturated when performed up to 3.1 GeV, which is the largest value covered by the BC fit. When again extending the BC fit beyond its fit range, for the purpose of a test, we checked that the relative contribution from the region The total TPE correction for μ − p elastic scattering is shown as sum of the elastic TPE, the TPE correction from the F 1 and F 2 proton structure functions and the TPE correction from the subtraction function of Ref. [13]. It is compared with the Feshbach term for pointlike particles, see Eq.
(2), corrected by the recoil factor (1 + m/M), the elastic contribution based on the box graph evaluation with dipole form factors [9] and the e − p total TPE correction [10] 3.1 GeV < W < 4 GeV to δ F 1 ,F 2 2γ is smaller than 1.5 %. We estimate the uncertainties of the numerical integration coming from the integration regions outside the BC fit and from the inaccuracies in the BC fit at 5-6 % level.
The resulting inelastic TPE corrections for the elastic μ − p scattering process are shown in Fig. 13 as a function of Q 2 for three values of muon beam momentum, corresponding with the MUSE kinematics. Note that the muon beam lab momenta k = 115 MeV, k = 153 MeV, and k = 210 MeV, correspond with the kinematically allowed regions of Q 2 < 0.039 GeV 2 , Q 2 < 0.066 GeV 2 , and Q 2 < 0.116 GeV 2 , respectively. We notice that for the small momentum transfers corresponding with the MUSE kinematics, the inelastic TPE corrections to elastic μ − p scattering are very small, in the range of δ 2γ ∼ 5 × 10 −4 . This is well below the anticipated cross section precision of around 1 % of the MUSE experiment. Furthermore, we notice that the TPE corrections due to the subtraction function and the dispersive F 1 , F 2 structure function integrals come with opposite signs, leading to a partial cancellation.
We present in Fig. 14 the total TPE correction as a sum of the Born TPE correction of Ref. [9], corresponding with a proton intermediate state, and the inelastic TPE of this work using the subtraction function of Birse et al. [13]. We compare our result with the Feshbach term of Eq. (2) for a point-like Dirac particle corrected by the recoil factor (1 + m/M), with the elastic TPE correction based on the box graph evaluation with proton form factors of the dipole form [9], and with the corresponding TPE correction for elastic e − p scattering of Ref. [10]. Contrary to the electron-proton scattering case, where the subtraction function contribution is negligible [10] as it is proportional to the lepton mass squared, in the case of muon-proton scattering the inelastic proton structure function contribution is partially canceled by the T 1 subtraction function resulting in a negligibly small inelastic TPE correction for the MUSE kinematics. Only with increased lepton beam energy or when going to larger Q 2 values one needs to start accounting for the inelastic TPE correction, which shifts the total correction a little closer to the Feshbach result.

Conclusions
In this work we have estimated the TPE correction to muonproton elastic scattering at low momentum transfer. For the elastic (proton) intermediate state contribution, we have derived a low-momentum transfer expansion accounting for all terms due to the non-zero lepton mass. Besides the elastic contribution, we have accounted for the inelastic intermediate states by expressing the TPE process at low momentum transfer approximately through the forward doubly virtual Compton scattering. The input in our evaluation of the inelastic TPE correction is given by the unpolarized proton structure functions and by one subtraction function, corresponding with the forward Compton amplitude T 1 at zero photon energy. For the latter, we have compared two estimates based on heavy-baryon and baryon chiral perturbation theory with an empirical determination. For the empirical determination, we have expressed the subtraction function through an unsubtracted dispersion relation for the amplitude T 1 −T R 1 . The function T R 1 is suitably defined through a Regge pole fit such that at high-energies (T 1 − T R 1 ) → 0, ensuring convergence and applicability of the unsubtracted dispersion relation. We have provided a numerical evaluation of the subtraction function based on a Regge fit of high-energy proton structure function data. It was found that the extracted subtraction function is compatible in magnitude with the chiral perturbation theory calculations, and thus cannot explain the proton radius puzzle through missing TPE corrections, which would have required a total TPE correction which is larger by around an order of magnitude compared with the empirical and chiral perturbation theory-based evaluations. Besides the subtraction function, the second part of the inelastic TPE contribution was obtained through dispersion integrals over the unpolarized proton structure functions. For the latter, we used a fit of the data in the proton resonance region. Using our formalism, we have provided estimates for the total TPE corrections in the kinematics of forthcoming muon-proton elastic scattering data of the MUSE experiment. We found that in the MUSE kinematics, the elastic TPE contribution largely dominates, and the size of the inelastic TPE contributions is within the anticipated error of the forthcoming data.
with the lepton momentum in the lab frame k ≡ |k|.
We also provide the more general expansion of Eq. (A1) in the low-Q 2 limit Q 2 M 2 , M 2 k 2 /s, where Q 2 needs not be very small relative to the squared lepton mass. For such an expansion, the leading Q 2 terms are given by  T . The experimental result R 0 = 0.23 ± 0.04 in the regionQ 2 > 1.5 GeV 2 from H1 and ZEUS [24] is connected with the ratio taken from the BC fit [20] (central curve). The error band reflects the experimental uncertainty in the value of R from the fit to the H1 and ZEUS data. The data points are from Refs. [25][26][27] into the H1/ZEUS value atQ 2 > 1.5 GeV 2 . We use the relative uncertainties from the data of Ref. [24] in the wholẽ Q 2 region. We show the resulting functional form of R(Q 2 ) in Fig. 15, and we compare its value with the data from Refs. [25][26][27] in the rangeQ 2 < 1.5 GeV 2 . We notice that our parameterization of R yields good agreement with the data.
Adopting the above Regge parameterization for F 2 , with the ratio R from Eq. (B8), we obtain from Eq. (B7) for F 1 the following Regge pole residues entering Eq. (B1):