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.


I. 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 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 a 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 it was done in the electron-proton elastic scattering [4,6]. Such an elastic scattering experiment is presently being planned by the MUSE Collaboration [7]. Complementary, 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 previous 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 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 lepton-proton scattering in Sec. II, 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 Sec. III. We evaluate the TPE correction due to the subtraction function in the forward Compton amplitude T 1 , and provide an empirical determination of the subtraction function from data in Sec. IV. 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 Sec. V. Our conclusions are given in Sec. VI.

II. 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 Q 2 M 2 , M E, where E is the lepton beam energy (in the lab frame) and FIG. 1: TPE graph in elastic lepton-proton scattering.
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 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 loop-momentum 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) T Born 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 1 (0,Q 2 ), which in Eq. (12) is defined as 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 [12,13].
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 Section IV, and the dispersive contribution due to the structure functions F 1 and F 2 in Sec. V.

IV. 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.

A. Heavy-baryon ChPT subtraction function
First of all, we show the fit of Ref. [12] obtained by matching the heavy-baryon chiral perturbation theory (HBChPT) result to a dipole behavior: with the value of the magnetic polarizability β M = (2.5±0.4)×10 −4 fm 3 taken from PDG [14].

C. 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 T 1 (ν,Q 2 )−T R 1 (ν,Q 2 ), 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 By comparing Eq. (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 [ Using Eqs. (8,16,17), this yields an expression for T subt 1 (0,Q 2 ) which expressed in terms of 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 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 high-energy behavior of the quantity  (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. [12,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 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 small 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. [12,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 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. [13]. was speculated in Ref. [22] 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. [22] 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. [22] for µ − p scattering, however, it differs by an overall sign.
In Fig. 9 we compare the TPE correction to elastic muon-proton 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) betweenQ 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 though slightly smaller result. This can be understood as the empirically determined β(Q 2 ) changes sign as 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 aroundQ 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. [22] 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. [22] 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.

V. 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,23].
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, 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 scattering 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 %.    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 lepton-proton 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 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 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. [12]. We compare our result with the 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. [12]. It is compared with the Feshbach term for point-like 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].
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. 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. point particle results in the TPE correction δ QED 2γ given by 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 expansion, the leading Q 2 terms are given by where Li 2 (x) denotes the dilogarithm function. The leading IR divergent piece is given by Eq. (A8). In the limit Q 2 m 2 , the result of Eq. (A9) reduces to the expression of Eq.
Appendix B: Regge poles residues of the proton structure function F 1 from high-energy data The high-energy limit (ν very large at fixedQ 2 ) of the proton structure function F 1 is often parameterized through a Regge pole fit as where γ α 0 (Q 2 ) are the leading Regge poles residues, which can be extracted from the highenergy inclusive electron-proton scattering data. We will determine these residues from the Donnachie-Landshoff (DL) fit [19] to data for the proton structure function F 2 in the region of very small Bjorken variable x Bj ≡Q 2 / (2Mν): where f 0 (Q 2 ) = A 0 with parameters values (using GeV units for all mass scales) [19]: The F 1 structure function in the high-energy region is then obtained: where R ≡ σ γp L /σ γp T is the ratio of longitudinal to transverse virtual photon absorption cross sections on a proton. We will use the experimental result R 0 = 0.23 ± 0.04 atQ 2 > 1.5 GeV 2 from the H1 and ZEUS collaborations [24], and approximate R in our numerical estimates by the following expression, independent of W 2 ≡ 2Mν + M 2 −Q 2 : R = R(Q 2 ) = R 0 Θ Q 2 − 1.5 GeV 2 + R BC (Q 2 )Θ −Q 2 + 1.5 GeV 2 , where R BC (Q 2 ) is value obtained in the Christy and Bosted fit [20] evaluated at W 2 ≈ 2.63 GeV 2 . The latter corresponds with the W 2 value for which the ratio R from the BC fit R BC (Q 2 = 1.5 GeV 2 ) ≈ 0.23, and thus goes over into the H1/ZEUS value atQ 2 > 1.5 GeV 2 .
We use the relative uncertainties from the data of Ref. [24] in the wholeQ 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): 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].