A dispersive approach to two-photon exchange in elastic electron-proton scattering

TL;DR

This paper presents a dispersive method to accurately calculate two-photon exchange corrections in elastic electron-proton scattering, incorporating nucleon and Delta states, and compares results with recent experimental data.

Contribution

It introduces a dispersive approach that avoids off-shell uncertainties and ensures correct high-energy behavior, improving the theoretical understanding of two-photon exchange effects.

Findings

Dispersive calculations agree with experimental cross section ratios.

Method reduces uncertainties compared to traditional loop diagram evaluations.

Results align with polarization transfer measurements.

Abstract

We examine the two-photon exchange corrections to elastic electron-proton scattering within a dispersive approach, including contributions from both Nucleon and Delta intermediate states. The dispersive analysis avoids off-shell uncertainties inherent in traditional approaches based on direct evaluation of loop diagrams, and guarantees the correct unitary behavior in the high energy limit. Using empirical information on the electromagnetic nucleon elastic and N-Delta transition form factors, we compute the two-photon exchange corrections both algebraically and numerically. Results are compared with recent measurements of e+p to e-p cross section ratios from the CLAS, VEPP-3 and OLYMPUS experiments, as well as with polarization transfer observables.

Tables

Table 1. Table 1: Table of values for I i ​ j subscript 𝐼 𝑖 𝑗 I_{ij} in Eq. ( 42 ), keeping only logarithmic terms in the infinitesimal regulator λ 𝜆 \lambda . The dependence on s 𝑠 s and other kinematic variables is suppressed for clarity.

$N$	$i$	$j$	$I_{ij}$
0	0	0	$d_0(\lambda ,\lambda )-d_0({\Lambda }_1,\lambda )-d_0(\lambda ,{\Lambda }_2)+d_0({\Lambda }_1,{\Lambda }_2)$
1	1	0	${\Lambda }_1^2\left(d_0({\Lambda }_1,\lambda )-d_0({\Lambda }_1,{\Lambda }_2)\right)$
1	0	1	${\Lambda }_2^2\left(d_0(\lambda ,{\Lambda }_2)-d_0({\Lambda }_1,{\Lambda }_2)\right)$
2	2	0	$-{\Lambda }_1^2\left(c_0\left(\lambda \right)-c_0\left({\Lambda }_2\right)+{\Lambda }_1^2\left(d_0({\Lambda }_1,\lambda )-d_0({\Lambda }_1,{\Lambda }_2)\right)\right)$
2	0	2	$-{\Lambda }_2^2\left(c_0\left(\lambda \right)-c_0\left({\Lambda }_1\right)+{\Lambda }_2^2\left(d_0(\lambda ,{\Lambda }_2)-d_0({\Lambda }_1,{\Lambda }_2)\right)\right)$
2	1	1	${\Lambda }_1^2{\Lambda }_2^2{d}_0({\Lambda }_1,{\Lambda }_2)$

Equations

$$s$$

$$t$$

$$u$$

$$\epsilon$$

$$\tau$$

$$\nu$$

$$s$$

$$\epsilon$$

$$\tau$$

$${\mathcal{M}}_{\gamma }=-\frac{e^2}{q^2}{j}_{\mu }^{\gamma }{J}_{\gamma }^{\mu },$$

$$j_{\mu }^{\gamma }$$

$$J_{\gamma }^{\mu }$$

$${\Gamma }^{\mu }(q)=F_1(Q^2){\gamma }^{\mu }+F_2(Q^2)\frac{i{\sigma }^{\mu \nu }{q}_{\nu }}{2M},$$

$${\sigma }_R^{\mathrm{Born}}=\epsilon G_E^2(Q^2)+\tau G_M^2(Q^2),$$

$$G_E(Q^2)$$

$$G_M(Q^2)$$

$${\mathcal{M}}_{\gamma \gamma }^{\mathrm{box}}=-i{e}^4\int \frac{d^4{q}_1}{(2\pi {)}^4}\frac{L_{\mu \nu }{H}_R^{\mu \nu }}{(q_1^2-{\lambda }^2)(q_2^2-{\lambda }^2)},$$

$$L_{\mu \nu }$$

$$H_R^{\mu \nu }$$

$$S_F(k_1,m_e)=\frac{(/k_1+m_e)}{(k_1^2-m_e^2+i{0}^{+})}.$$

$$S_{\alpha \beta }(p_R,M_R)={\delta }_{\alpha \beta }{S}_F(p_R,M_R),$$

$$S_{\alpha \beta }(p_R,M_R)=-S_F(p_R,M_R){\mathcal{P}}_{\alpha \beta }^{3/2}(p_R),$$

$${\mathcal{P}}_{\alpha \beta }^{3/2}(p_R)=g_{\alpha \beta }-\frac{1}{3}{\gamma }_{\alpha }{\gamma }_{\beta }-\frac{1}{3p_R^2}\left(/p_R{\gamma }_{\alpha }(p_R{)}_{\beta }+(p_R{)}_{\alpha }{\gamma }_{\beta }/p_R\right),$$

$${\cal M}_{{\gamma\gamma}}^{\rm xbox}(u,t)=-{\cal M}_{{\gamma\gamma}}^{\rm box}(s,t)\Big{|}_{s\to u}\,.$$

$${\left|{\mathcal{M}}_{\gamma }+{\mathcal{M}}_{\gamma \gamma }\right|}^2$$

$$\equiv {\left|{\mathcal{M}}_{\gamma }\right|}^2\left(1+{\delta }_{\mathrm{TPE}}\right),$$

$${\delta }_{\mathrm{TPE}}=\frac{2\mathrm{Re}\left({\mathcal{M}}_{\gamma }^{†}{\mathcal{M}}_{\gamma \gamma }\right)}{{\left|{\mathcal{M}}_{\gamma }\right|}^2}.$$

$${\mathcal{M}}_{\gamma \gamma }⟶{\mathcal{M}}_{\gamma }{\Delta }_{\mathrm{IR}},$$

$${\delta }_{\gamma \gamma }\equiv {\delta }_{\mathrm{TPE}}-2\mathrm{Re}{\Delta }_{\mathrm{IR}}.$$

$${\Delta }_{\mathrm{IR}}(\mathrm{MTj})=-\frac{\alpha }{\pi }\log \left(\frac{M^2-s}{M^2-u}\right)\log \left(\frac{Q^2}{{\lambda }^2}\right).$$

$$\mathrm{Re}{\Delta }_{\mathrm{IR}}(\mathrm{MTj})$$

$$\mathrm{Im}{\Delta }_{\mathrm{IR}}(\mathrm{MTj})$$

$${\sigma }_R={\sigma }_R^{\mathrm{Born}}\left(1+{\delta }_{\gamma \gamma }\right).$$

$${\mathcal{M}}_{\gamma \gamma }$$

$$F_1^{\prime }$$

Full text

††preprint: JLAB-THY-17-2430

1. I Introduction

2. II Formalism

3. II.1 Elastic $ep$ scattering

4. II.1.1 One-photon exchange

5. II.1.2 Two-photon exchange

6. II.2 Dispersive approach

7. II.2.1 Analytic method

8. II.2.2 Numerical method

9. III Impact of two-photon exchange

10. III.1 Nucleon intermediate state

11. III.2 $\Delta$ intermediate state

12. III.3 Verification of TPE effects

13. III.3.1 $e^{+}p$ to $e^{-}p$ ratio

14. III.3.2 Polarization observables

15. IV Outlook

16. A Generalized form factors

17. B Analytic expressions for imaginary parts of Passarino-Veltman functions

18. C Form factor reparametrizations

I Introduction

The nucleon’s electroweak form factors are some of the cornerstone observables that characterize its extended spatial structure. Since the original observation Hofstadter and McAllister (1955) some 60 years ago that elastic scattering from the proton deviates from point-like behavior at large scattering angles, considerable information has been accumulated on the detailed structure of the proton’s and neutron’s electric and magnetic responses. Almost universally the underlying scattering reaction has been assumed to proceed through the exchange of a single photon between the lepton (typically electron) beam and nucleon target.

A major paradigm shift occurred around the turn of the last century with the observation of a significant discrepancy between the ratio of electric to magnetic form factors of the proton measured using the relatively new polarization transfer technique Jones et al. (2000); Gayou et al. (2002) and previous extractions of the same quantity from cross section measurements via Rosenbluth separation. It was soon realized Blunden et al. (2003); Guichon and Vanderhaeghen (2003) that a large part of the discrepancy could be understood in terms of additional, hadron structure-dependent two-photon exchange (TPE) contributions, which had not been included in the standard treatments of electromagnetic radiative corrections Tsai (1961); Mo and Tsai (1969).

A number of approaches have been adopted to computing the TPE corrections to elastic scattering, including direct calculation of the loop contributions in terms of hadronic degrees of freedom Blunden et al. (2003); Kondratyuk et al. (2005); Blunden et al. (2005); Kondratyuk and Blunden (2007); Nagata et al. (2009); Graczyk (2013); Lorenz et al. (2015); Zhou and Yang (2015); Zhou (2017), modeling the high energy behavior of box diagrams at the quark level through generalized parton distributions Chen et al. (2004); Afanasev et al. (2005), or more recently dispersion relations Gorchtein (2007); Borisyuk and Kobushkin (2008, 2015); Tomalak and Vanderhaeghen (2015). Each of these methods has its own advantages as well as limitations (for reviews, see Refs. Carlson and Vanderhaeghen (2007); Arrington et al. (2011); Afanasev et al. (2017)), and to date no single approach has been able to provide a universal description valid at all kinematics.

Most of the attention on the TPE corrections in recent experiments has been focussed on the region of small and intermediate values of the four-momentum transfer squared, $Q^{2}\lesssim$ few GeV2, where the expectation is that hadrons retain their identity sufficiently well that calculations in terms of physical degrees of freedom give reliable estimates. Traditionally, this approach has required direct evaluation of the real parts of the two-photon box and crossed-box diagrams, with nucleons or other excited state hadrons in the intermediate state parametrized through half off-shell form factors (with one nucleon on-shell and one off-shell). Because the off-shell dependence of these form factors is not known, usually one approximates the half off-shell form factors by their on-shell limits.

For nucleon intermediate states, the off-shell uncertainties are not expected to be severe. On the other hand, for transitions to excited state baryons described by effective interactions involving derivative couplings, such as for the $\Delta$ resonance, the off-shell dependence leads to divergences in the forward angle (or high energy) limit, and signals a violation of unitarity. Furthermore, from a more technical perspective, in order to evaluate the TPE corrections analytically in terms of Passarino-Veltman functions, the loop integration method requires the transition form factors to be parametrized as sums or products of monopole functions. This can prove cumbersome in some applications, since such parametrizations are usually only valid in a limited region of spacelike $Q^{2}$, and may be prone to roundoff errors in numerical evaluation. It would naturally be highly desirable to be able to compute the loop integrations with a more robust numerical method that is valid for form factor parametrizations based on more general classes of functional forms.

The limitations of the previous loop calculations are especially problematic in view of new measurements of ratios of $e^{+}p$ to $e^{-}p$ elastic scattering cross sections Rimal et al. (2016); Rachek et al. (2015); Nikolenko et al. (2014); Henderson et al. (1969), which have provided high precision data that are directly sensitive to TPE effects. Some of these data are in the small-angle region, where the off-shell ambiguities in the loop calculations make the calculations unreliable. To enable meaningful comparison between the data and TPE calculations over the full range of kinematics currently accessible, clearly a different approach to the problem is needed.

In this paper, we revisit the calculation of TPE corrections within the hadronic approach, but using dispersion relations to construct the real part of the TPE amplitude from its imaginary part. The dispersive method involves the exclusive use of on-shell transition form factors, thereby avoiding the problem of unphysical violation of unitarity in the high energy limit. The dispersive approach to TPE was developed at forward angles by Gorchtein Gorchtein (2007), and at non-forward angles by Borisyuk and Kobushkin Borisyuk and Kobushkin (2008, 2012, 2014, 2015), and more recently by Tomalak and Vanderhaeghen Tomalak and Vanderhaeghen (2015). A feature of the latter two analyses has been the use of monopole form factor parametrizations, which allowed the computations to be performed semi-analytically. In this work we extend the dispersion relation approach to allow for more general classes of transition form factors.

In Sec. II of this paper we review the formalism for elastic electron–nucleon scattering for both one-photon and two-photon exchange processes, and introduce the main elements of the dispersive approach. We describe analytical calculations of the imaginary part of the TPE corrections using the more restrictive monopole form factors, for which one can obtain analytic expressions in terms of elementary logarithms. We also describe the more general numerical method that allows standalone calculation of the imaginary part using a general class of transition form factors.

The results of the calculations are presented in Sec. III, where we critically examine the differences between the new dispersive method and the previous loop calculations with off-shell intermediate states. While the differences are relatively small for the nucleon elastic contributions, the effects for $\Delta$ intermediate states are dramatic at high energies and forward scattering angles. We also compare in Sec. III the new results with the recent data on $e^{+}p$ to $e^{-}p$ cross section ratios from the CLAS Rimal et al. (2016), VEPP-3 Rachek et al. (2015); Nikolenko et al. (2014) and OLYMPUS Henderson et al. (1969) experiments, as well as with polarization data sensitive to TPE contributions Meziane et al. (2011). Finally, in Sec. IV we summarize our results, and discuss possible future developments in theory and experiment. For completeness, in the appendices we give the full expressions for the generalized form factors in Appendix A, and analytic expressions for the imaginary parts of Passarino-Veltman functions in Appendix B. We also provide convenient reparametrizations of the nucleon and $\Delta$ vertex form factors in Appendix C that can be used in the analytic calculations.

II Formalism

In this section we present the formalism on which the electron–nucleon scattering analysis in this paper will be based. After summarizing the kinematics and main formulas for the elastic scattering amplitudes and cross sections at the Born and TPE level, we proceed to describe the new elements of the analysis that make use of dispersive methods, including both analytic and numerical evaluation of integrals.

II.1 Elastic $ep$ scattering

For the elastic scattering process $ep\to ep$ the four-momenta of the initial and final electrons (taken to be massless) are labeled by $k$ and $k^{\prime}$, and of the initial and final protons (mass $M$) by $p$ and $p^{\prime}$, respectively, as depicted in Fig. 1. The four-momentum transfer from the electron to the proton is given by $q=p^{\prime}-p=k-k^{\prime}$, with $Q^{2}\equiv-q^{2}>0$. One can express the elastic cross section in terms of any two of the Mandelstam variables $s$ (total electron–proton invariant mass squared), $t$, and $u$, where

[TABLE]

with the constraint $s+t+u=2M^{2}$.

The elastic scattering cross section can be defined in terms of any two of the dimensionless quantities

[TABLE]

The inverse relationships are also useful,

[TABLE]

In the target rest frame the variables are given by

[TABLE]

where $E$ $(E^{\prime})$ is the energy of the incident (scattered) electron, $\theta_{e}$ is the electron scattering angle, and $\varepsilon$ ($0<\varepsilon<1$) is identified with the relative flux of longitudinal virtual photons.

II.1.1 One-photon exchange

In the Born (OPE) approximation the electron–nucleon scattering invariant amplitude can be written as

[TABLE]

where $e$ is the electric charge, and the matrix elements of the electromagnetic leptonic and hadronic currents are given in terms of the lepton ($u_{e}$) and nucleon ($u_{N}$) spinors by

[TABLE]

The electromagnetic hadron current operator $\Gamma^{\mu}$ is parametrized by the Dirac ($F_{1}$) and Pauli ($F_{2}$) form factors as

[TABLE]

where the Born form factors are functions of a single variable, $Q^{2}$. In our convention, the reduced Born cross section $\sigma_{R}^{\rm Born}$ is given by

[TABLE]

where the Sachs electric and magnetic form factors $G_{E,M}(Q^{2})$ are defined in terms of the Dirac and Pauli form factors as

[TABLE]

II.1.2 Two-photon exchange

Using the kinematics illustrated in the box diagram in Fig. 1(b), the contribution to the TPE box amplitude from an intermediate hadronic state $R$ of invariant mass $M_{R}$ can be written in the general form Blunden et al. (2003); Kondratyuk et al. (2005)

[TABLE]

with $q_{2}=q-q_{1}$, and an infinitesimal photon mass $\lambda$ is introduced to regulate any infrared divergences. (In general the mass $M_{R}$ can have a distribution $W$ which can be integrated over, but here we specialize to the case of a narrow state $R$.) The leptonic and hadronic tensors here are given by

[TABLE]

with $k_{1}=k-q_{1}$, $p_{R}=p+q_{1}$, and the electron propagator is

[TABLE]

The hadronic transition current operator $\gamma N\to R$ is written in a general form $\Gamma_{\gamma N\to R}^{\alpha\mu}(p_{R},q)$ that allows for a possible dependence on the incoming momentum $q$ of the photon and the outgoing momentum $p_{R}$ of the hadron, while $\mu$ and $\alpha$ are Lorentz indices.

The hadronic state propagator $S_{\alpha\beta}(p_{R},M_{R})$ in this work will describe the propagation of a baryon with either spin-$\nicefrac{{1}}{{2}}$ or spin-$\nicefrac{{3}}{{2}}$. For spin-$\nicefrac{{1}}{{2}}$ intermediate states, such as the nucleon, this reduces to

[TABLE]

and the transition operator $\Gamma_{\gamma N\to R}$ involves one free Lorentz index. For spin-$\nicefrac{{3}}{{2}}$ intermediate states, such as the $\Delta$ baryon, the propagator can be written

[TABLE]

where the projection operator

[TABLE]

ensures the presence of only spin-$\nicefrac{{3}}{{2}}$ components. Unphysical spin-$\nicefrac{{1}}{{2}}$ contributions are suppressed by the condition on the vertex $p_{R\alpha}\Gamma_{\gamma N\to R}^{\alpha\mu}(p_{R},q)=0$.

One can obtain the crossed-box (“xbox”) contribution directly from the box term (10) by applying crossing symmetry. For example, in the unpolarized case, we have

[TABLE]

In general, ${\cal M}_{{\gamma\gamma}}^{\rm box}(s,t)$ has both real and imaginary parts, whereas ${\cal M}_{{\gamma\gamma}}^{\rm xbox}(u,t)$ is purely real. The total squared amplitude for the sum of the one- and two-photon exchange processes shown in Fig. 1 is then

[TABLE]

where the relative correction to the cross section due to the interference of the one- and two-photon exchange amplitudes is defined as

[TABLE]

Within the framework of the simplest hadronic models, analytic evaluation of $\delta_{\rm TPE}$ is made possible by writing the transition form factors at the $\gamma$-hadron vertices as a sum and/or product of monopole form factors Blunden et al. (2003, 2005), which are typically fit to empirical transition form factors over a suitable range in spacelike four-momentum transfer. Four-dimensional integrals over the momentum in the one-loop box diagram can then be expressed in terms of the Passarino-Veltman (PV) scalar functions $A_{0}$, $B_{0}$, $C_{0}$ and $D_{0}$ ’t Hooft and Veltman (1979); Passarino and Veltman (1979). This reduction to scalar integrals is automated by programs such as FeynCalc Mertig et al. (1991); Shtabovenko et al. (2016). The PV functions can then be evaluated numerically using packages such as LoopTools Hahn and Perez-Victoria (1999). In this paper we are interested only in the imaginary parts of these PV functions, which are considerably simpler than the full expressions. This will be discussed in detail in the next section.

Note that the expressions (10) and (19) contain infrared (IR) divergences arising from the elastic intermediate state when the momentum $q_{i}\ (i=1,2)$ of either photon vanishes. In analyzing the TPE corrections for $ep$ scattering, it is convenient to subtract off these conventional IR-divergent parts, which are independent of hadronic structure, and which are usually already included in experimental analyses using a specific prescription (e.g. Mo & Tsai Mo and Tsai (1969), Grammer & Yennie Grammer and Yennie (1973), or Maximon & Tjon Maximon and Tjon (2000)).

In general, the TPE amplitude at the IR poles ($q_{1}\to 0$ or $q_{2}\to 0$) has the form

[TABLE]

where $\Delta_{\rm IR}$ is a function containing all the IR divergences that is independent of hadronic structure. Its form depends on the particular IR prescription being used. This is discussed extensively in the TPE review by Arrington, Blunden, and Melnitchouk Arrington et al. (2011), and we defer to that paper for details. The hard-TPE correction of interest is then

[TABLE]

In this paper we follow the prescription used by Maximon and Tjon Maximon and Tjon (2000), which is to evaluate the contribution to the numerator of Eq. (10) arising from the poles $q_{i}\to 0$, while keeping the propagators in the denominator intact. In this prescription,

[TABLE]

This expression has both real and imaginary parts. In our convention, $\log(-x)=\log{x}-i\pi$ for $x>0$, so explicitly the real and imaginary parts are

[TABLE]

After accounting for conventional radiative corrections, the measured reduced cross section $\sigma_{R}$ is related to the Born cross section by

[TABLE]

In practice, most experimental cross section analyses use the IR-divergent expression of Mo and Tsai Mo and Tsai (1969), so that if one uses the Maximon and Tjon prescription Maximon and Tjon (2000) (as we do in this paper) then the difference should be accounted for when comparing to experimental data (see Ref. Arrington et al. (2011) for further discussion).

The total TPE amplitude ${\cal M}_{{\gamma\gamma}}$ can be rewritten in terms of “generalized form factors”, generalizing the expressions of Eqs. (5)–(7), as described by Guichon and Vanderhaeghen Guichon and Vanderhaeghen (2003). Although the decomposition is not unique, and different generalized form factor conventions have been used in the literature, in this paper we use the basis of form factors denoted by $F_{1}^{\prime}$, $F_{2}^{\prime}$ and $G_{a}^{\prime}$, defined via

[TABLE]

where the vector $F_{1}^{\prime}$ and $F_{2}^{\prime}$ generalized form factors are the TPE analogs of the Dirac and Pauli form factors, while the axial vector $G_{a}^{\prime}$ generalized form factor has no Born level analog.

Rather than construct $\delta_{\rm TPE}$ and subtract the IR-divergent terms, as in Eq. (21), it is convenient to incorporate the IR subtractions directly into the generalized form factors $F_{1}^{\prime}$ and $F_{2}^{\prime}$ ($G_{a}^{\prime}$ is not IR-divergent),

[TABLE]

where $F_{i,{\rm TPE}}^{\prime}$ refer to the unregulated expressions. In terms of these regulated generalized form factors, the relative TPE correction is given by

[TABLE]

II.2 Dispersive approach

As noted earlier, the TPE amplitude ${\cal M}_{\gamma\gamma}$ has both real and imaginary parts. The real and imaginary parts can be related through dispersion relations Gorchtein (2007); Borisyuk and Kobushkin (2008), which forms the basis of the dispersive method discussed in this section. Our discussion in this section follows the formalism of Tomalak and Vanderhaeghen Tomalak and Vanderhaeghen (2014, 2015). An alternative treatment by Borisyuk and Kobushkin Borisyuk and Kobushkin (2008) starts from the annihilation channel, $e^{-}+e^{+}\to p+\bar{p}$.

Using the parametrization of the TPE amplitude ${\cal M}_{\gamma\gamma}$ in terms of the generalized form factors $F_{1}^{\prime}$, $F_{2}^{\prime}$ and $G_{a}^{\prime}$, we note that these TPE amplitudes have the symmetry properties Gorchtein (2007); Borisyuk and Kobushkin (2008)

[TABLE]

and satisfy the fixed-$t$ dispersion relations

[TABLE]

Here ${\cal P}$ denotes the Cauchy principal value integral, and $\nu_{\rm th}=-\tau$ is the threshold for the elastic cut, corresponding to an electron of energy $E=0$. The physical threshold for electron scattering is at $\varepsilon=0$ (or $\cos\theta_{e}=-1$), which requires $E\geq M(\tau+\nu_{\rm ph})$, with $\nu_{\rm ph}\equiv\sqrt{\tau(1+\tau)}$. This integral therefore extends into an unphysical region of parameter space, which requires knowledge of the transition form factors in the timelike region of four-momentum transfer. The crossed-box terms in the real part of the TPE amplitudes are generated by incorporating the symmetry properties into the dispersive integrals, which is equivalent to the use of Eq. (17) in the loop calculation.

For the interaction of point particles, such as in elastic $e\mu$ scattering, the real parts generated in this way agree completely with those obtained directly from the four-dimensional loop integrals of Eq. (10) Tomalak and Vanderhaeghen (2014). In general, however, there may be momentum dependence in the $\gamma$-hadron interaction, such as for the $\gamma N\Delta$ vertex (see Sec. III.2 below). In fact, the momentum dependence in a transition vertex function allows one to construct different parametrizations of that vertex function, for example, by using the Dirac equation, that are equivalent on-shell but differ off-shell. The additional momentum-dependence associated with this freedom will affect one-loop integrals because the intermediate hadronic states are not on-shell. This ambiguity is not present in the dispersive method, for which all the intermediate states are on-shell. In the context of TPE, this means that for any momentum-dependent interactions one should not expect agreement between the real parts of the generalized form factors calculated from Eqs. (29) and those calculated using the loop integration method. We will quantify these differences for the cases of the nucleon and $\Delta$ intermediate states in Sec. III.

II.2.1 Analytic method

The analytic approach used in previous work Blunden et al. (2003); Kondratyuk et al. (2005); Blunden et al. (2005); Kondratyuk and Blunden (2007); Tjon and Melnitchouk (2008); Nagata et al. (2009); Tjon et al. (2009); Graczyk (2013); Lorenz et al. (2015); Zhou and Yang (2015) relies on a parameterization of the transition form factors as a sum and/or product of monopole form factors. The most basic relation is

[TABLE]

which is to be applied at each photon–hadron vertex ($i=1,2$). More complicated constructions are straightforward to generate by repeatedly applying the feature that the product of any two monopoles is proportional to their difference. The general expression for an amplitude with form factors will thus involve a sum of “primitive” integrals with different photon mass parameters $\Lambda_{1}$ and $\Lambda_{2}$ for each of the two photon propagators, modified according to Eq. (30). The primitive integrals may yield spurious ultraviolet or infrared divergences, but these divergences will cancel when taking the sum. We give details on these constructions in Appendix B.

By means of the PV reduction scheme, a one-loop integral for the box-diagram amplitude can be written in terms of a set of scalar PV functions $A_{0}$, $B_{0}$, $C_{0}$ and $D_{0}$, corresponding to one-, two-, three-, and four-point functions. This can be visualized as a “pinching” of the four various propagators in the box diagram due to cancellations of the terms in the numerator with the propagator terms in the denominator. The PV functions can be evaluated numerically, and there are various computer programs to do this Hahn and Perez-Victoria (1999); van Hameren (2011); Carrazza et al. (2016); Denner et al. (2017).

In general the scalar PV functions are complex-valued. The imaginary parts of the TPE amplitudes are contained entirely in these functions. For the box (and crossed-box) diagrams in elastic $ep$ scattering there are only four of the PV functions that have imaginary parts. These four functions are the ones that arise in the $s$-channel box diagram with the electron and intermediate hadronic states on-shell. This is illustrated in Fig. 2.

Recall that an amplitude becomes imaginary when the intermediate state particles become real, or on their mass shells. This is formalized by the well-known Cutkosky cutting rules Cutkosky (1958). Namely, as a consequence of unitarity, the imaginary part of a scattering amplitude can be obtained by summing all possible cuttings of the corresponding Feynman diagram, where a cut is across any two internal propagators separating the external states from the rest of the diagram. Cut propagators are then put on-shell according to the rule $1/(p^{2}-m^{2}+i0^{+})\ \to\ -2\pi i\,\theta(p_{0})\,\delta(p^{2}-m^{2})$.

For elastic scattering, the two functions $C_{0}(s)$ arising when either photon propagator is pinched are identical, so there are only three distinct PV functions. Other PV functions where the electron or hadronic intermediate state (or both) are pinched have no imaginary parts for $ep$ scattering, and the $u$-channel crossed-box diagram also has no imaginary part. (Recall that the crossed-box amplitude can be obtained by replacing $s\to u$, with an appropriate overall changed in sign given in Eq. (28).) We will denote these three functions as $B_{0}(s)$, $C_{0}(s;\Lambda_{1})$, and $D_{0}(s;\Lambda_{1},\Lambda_{2})$. The full expression for these functions is

[TABLE]

In addition to the explicit dependence on $s$ and $\Lambda_{i}^{2}$, there is also an implied dependence on $M$, $W$, and $Q^{2}$ that is suppressed for clarity of notation (see Appendix B for details).

We define the imaginary parts of the PV functions by $\{b_{0}(s),c_{0}(s),d_{0}(s)\}\equiv\{\operatorname{Im}B_{0}(s),\operatorname{Im}C_{0}(s),\operatorname{Im}D_{0}(s)\}$. According to the Cutkosky rules, the imaginary parts correspond to putting the electron and intermediate hadronic states on-shell, $k_{1}^{2}\equiv(k-q_{1})^{2}=m_{e}^{2}$ and $(p+q_{1})^{2}=W^{2}$. Working in the center-of-mass (CM) frame, we define the electron variables

[TABLE]

In this frame we have

[TABLE]

with the shorthand notation

[TABLE]

In the physical region, the CM scattering angle $\theta$ satisfies the constraint $-1\leq\cos\theta\leq 1$, requiring $s_{M}^{2}\geq sQ^{2}$. However, the dispersive integral of Eq. (29) only requires $s_{M}>0$, meaning there is an unphysical region of parameter space where $\cos\theta<-1$, and $\sin\theta$ is purely imaginary. Therefore, in the dispersive approach expressions for the imaginary parts of the TPE amplitudes need to be analytically continued into this unphysical region.

Recall that in terms of the electron energy in the laboratory frame, $E$, the $ep$-invariant mass squared is $s=M^{2}+2ME$. After changing the integration variable from $q_{1}$ to $k_{1}$, and using the on-shell conditions, we find after some algebra the expressions

[TABLE]

where $Q_{i}^{2}=-q_{i}^{2}$ are the squared four-momenta of the virtual photons ($i=1,2$), with

[TABLE]

and $Q_{0}^{2}=2E_{k}E_{k_{1}}=s_{M}s_{W}/(2s)$.

The $b_{0}$ integral in Eq. (35) is trivial. Through the use of Eq. (36), the other integrals can be brought into the form

[TABLE]

The integrand here has poles when $|b_{2}|>|a_{2}|$, which can arise in the unphysical region when $Q_{2}^{2}$ becomes timelike. A simpler version of this integral was considered by Mandelstam Mandelstam (1958) for the case where the target and scattering particles have equal masses. The general expression has been given by Beenakker and Denner Beenakker and Denner (1990),

[TABLE]

For the $c_{0}(s)$ function, we set $a_{1}=1+\Lambda_{1}^{2}/Q_{0}^{2}$, $b_{1}=1$, $a_{2}=1$, $b_{2}=0$, and $c_{2}=0$. For $d_{0}(s)$, we set $a_{i}=1+\Lambda_{i}^{2}/Q_{0}^{2}$, $b_{1}=1$, $b_{2}=\cos\theta$, and $c_{2}=\sin\theta=\sqrt{1-b_{2}^{2}}$. In the unphysical region, $\cos\theta<-1$, so that $\sin\theta$ is purely imaginary. However, we note that the combination $b_{2}^{2}+c_{2}^{2}=1$, and therefore $0\leq X\leq(a_{1}a_{2}-b_{1}b_{2})$, independent of the value of $\cos\theta$. Thus Eq. (38) for $J$ is the proper analytic continuation of the integral for $d_{0}(s)$ into the unphysical region. Explicit expressions for $b_{0}(s)$, $c_{0}(s;\Lambda)$ and $d_{0}(s;\Lambda_{1},\Lambda_{2})$, including the IR limits $\Lambda\to\lambda$, are given in Appendix B.

In previous work Blunden et al. (2003, 2005) the TPE amplitudes were obtained by numerical evaluation of the PV functions using the program LoopTools Hahn and Perez-Victoria (1999). The real parts were used directly, and the imaginary parts were not needed. Here, we have constructed analytic expressions for the imaginary parts in terms of elementary logarithms, thus allowing a completely analytic evaluation of the imaginary parts of the TPE amplitudes. The real parts are then constructed from a numerical evaluation of the dispersion integrals of Eq. (29). The imaginary parts obtained here are, of course, identical with those obtained numerically in the earlier work Blunden et al. (2003, 2005). The real parts are numerically identical for the elastic $F_{1}^{\prime}$ and $F_{2}^{\prime}$ TPE amplitudes, while the $G_{a}^{\prime}$ amplitude differs, but in a numerically insignificant way. For the inelastic $\Delta$ states there are significant differences, especially as $\varepsilon\to 1$. These differences will be discussed further in Sec. III.

As an alternative to using the PV reduction method implemented in FeynCalc Mertig et al. (1991); Shtabovenko et al. (2016), one can work entirely with on-shell quantities. Using the on-shell conditions, we find that the TPE amplitudes are sums of integrals of the general form

[TABLE]

where $f(Q_{1}^{2},Q_{2}^{2})$ is a polynomial function of combined degree $N$ in $Q_{1}^{2}$ and $Q_{2}^{2}$,

[TABLE]

The coefficients $f_{ij}$ are functions of $s$, $W$, and $Q^{2}$, and satisfy $f_{ij}=f_{ji}$ for elastic scattering due to the symmetry under $Q_{1}^{2}\leftrightarrow Q_{2}^{2}$. Thus we can write

[TABLE]

with the “primitive” integrals $I_{ij}$ defined as

[TABLE]

For nucleon intermediate states we find $N=2$, indicating that monopole form factors are sufficient to eliminate the UV divergences (there is one power of $q_{i}$ at each photon–nucleon vertex in the $F_{2}$ term of $\Gamma^{\mu}(q_{i})$). For $\Delta$ intermediate states, however, we find $N=3$, which implies that dipole form factors (or a product of monopoles) are needed to eliminate the UV divergences, as there are up to two powers of $q_{i}$ at each vertex in $\Gamma_{\gamma N\to\Delta}^{\alpha\mu}$ (see Secs. III.1 and III.2 below). The integrals of Eq. (42) up to $N=2$ are given in Table 1 of Appendix B, and can easily be extended to more complicated form factor constructions.

II.2.2 Numerical method

In analogy with Eq. (39), the TPE amplitudes of interest have the general form

[TABLE]

where $f(Q_{1}^{2},Q_{2}^{2})$ is a polynomial function of combined degree 2 (3) in $Q_{1,2}^{2}$ for $N~{}(\Delta)$ intermediate states, and $G_{i}(Q_{i}^{2})$ are form factors that are real-valued and finite for all spacelike values of $Q_{i}^{2}$. For elastic scattering, the total integral is symmetric under the interchange $Q_{1}^{2}\leftrightarrow Q_{2}^{2}$. For the nucleon intermediate state, it is convenient to bring the IR subtractions of Eq. (26) into this integral. This is consistent with the Maximon and Tjon IR regularization scheme whereby the numerator of Eq. (43) vanishes whenever $Q_{i}^{2}\to 0$. It also vanishes for excited states under these conditions. Therefore we could actually set $\lambda\to 0$ without encountering any singularities in the integrals. This is unlike the analytic expressions of the previous section, where only the sum of individual IR-divergent expressions is independent of $\lambda$.

In the physical region there are no singularities in the integrand of Eq. (43), so evaluation of the integral is a straightforward 2-dimensional numerical quadrature over the domain $-1\leq\cos\theta_{k_{1}}\leq 1$ and $0\leq\phi_{k_{1}}\leq 2\pi$, following Eq. (36). However, this approach fails in the unphysical region. To get around this, Tomalak and Vanderhaeghen used a contour integration method Tomalak and Vanderhaeghen (2015), and applied it to the calculation of TPE amplitudes with monopole form factors. By summing the residue at the poles enclosed by the contour they were able to obtain algebraic expressions for the TPE amplitudes. These expressions are equivalent to the ones we obtained in the previous section using algebraic expressions for the PV functions. In this section we will follow this method, with modifications, to implement a numerical contour integration of Eq. (43) that allows for a more general parametrization of the transition form factors than a sum and/or product of monopoles.

In the complex $Q_{i}^{2}$ plane, we define a timelike half with $\operatorname{Re}Q_{i}^{2}<0$, and spacelike half with $\operatorname{Re}Q_{i}^{2}>0$. In general, the allowed form factors can have poles in $Q_{i}^{2}$ anywhere in the timelike half of the complex plane. With certain restrictions, which we will state explicitly, the form factors can have poles in the spacelike half as well.

Without providing a rigorous mathematical proof, we can nonetheless specify certain restrictions on the type of allowed form factors. Namely, they should have a simple functional form in $Q_{i}^{2}$, such as exponentials, polynomials, or inverse polynomials, that can be analytically continued to the complex plane. There should be no branch cuts, and any poles should either lie along the negative, real axis (at timelike $Q_{i}^{2}$), or occur in complex conjugate pairs. This is the case for a commonly used form factor parametrization [see Eq. (51)] in terms of a ratio of polynomials Kelly (2004); Arrington et al. (2007); Venkat et al. (2011).

The area of integration in Eq. (43) can be visualized as an integral over the photon virtual momenta $Q_{1}^{2}$ and $Q_{2}^{2}$ of Eq. (36), which form a symmetric ellipse in $Q_{1}^{2}$ vs. $Q_{2}^{2}$, centered at $\{Q_{0}^{2},Q_{0}^{2}\}$ Gorchtein and Horowitz (2008); Tomalak and Vanderhaeghen (2015). The boundary of the ellipse is defined by $\cos\phi_{k_{1}}=1$. Following Ref. Tomalak and Vanderhaeghen (2015), we make a change of variables to elliptic coordinates,

[TABLE]

The contours of constant $\alpha$ represent concentric ellipses with radial parameter $r=\sqrt{1-\alpha^{2}}$. From Eq. (36), in elliptic coordinates we have

[TABLE]

In the physical region the integral over $\theta_{k_{1}}$ can be rewritten as a contour integral on the unit circle $z=\exp(i\theta_{k_{1}})$, with $Q_{i}^{2}(z)$ regarded as functions of $z$ (and $r$) using

[TABLE]

and

[TABLE]

In anticipation of extending this formalism to the unphysical region, Eq. (45) can be simplified into the form

[TABLE]

with

[TABLE]

Recall that $\cos\theta$ is given in Eq. (33) in terms of $s$ and $Q^{2}$. Expressed in this way, it is clear that $Q_{1}^{2}(z)=Q_{2}^{2}(\beta z)$. Because the integrands of the full TPE amplitudes are symmetric under the interchange $Q_{1}^{2}\leftrightarrow Q_{2}^{2}$, this means that for every value $z_{1}$ associated with the poles in $Q_{1}^{2}$ in the complex $z$ plane, there is a corresponding pole at $z_{2}=\beta z_{1}$. In addition, for every pole in $Q_{1}^{2}$, both $z_{\rm 1i}$ and $z_{\rm 1o}=1/z_{\rm 1i}$ are poles in the complex $z$ plane, where $|z_{\rm 1i}|\leq 1$ lies inside the unit circle, and $|z_{\rm 1o}|\geq 1$ lies outside the unit circle. The values $z_{\rm 1i}=z_{\rm 1o}=1$ are the IR-divergent poles, where $Q_{1}^{2}=0$ for $r=1$. These could be regulated by introducing the $\lambda$ photon mass parameter, but as the integrand vanishes in our regularization scheme when $Q_{1}^{2}=0$, there is no IR divergence, and one can set $\lambda=0$. Analogously, for $Q_{2}^{2}$, both $z_{\rm 2i}$ and $z_{\rm 2o}=\beta^{2}/z_{\rm 2i}$ are poles inside and outside a circle of radius $|\beta|$, respectively. The values $z_{\rm 2i}=z_{\rm 2o}=\beta$ are the IR-divergent poles, where $Q_{2}^{2}=0$ for $r=1$. These results are valid in either the physical region, where $\beta=\exp(i\theta)$, or the unphysical region, where $\beta<-1$.

As an illustration of these points, consider the monopole form factors as given in Eq. (39), for which there is a pole in $Q_{1}^{2}$ along the negative real axis at $-\Lambda_{1}^{2}$. This yields poles $z_{\rm 1i}$ and $z_{\rm 1o}$ along the positive real axis in the complex $z$ plane, with the interior pole $z_{\rm 1i}$ lying between 0 and 1. For the corresponding poles associated with $Q_{2}^{2}$, in the physical region, with $\beta=\exp(i\theta)$, these lie along a line at angle $\theta$. In the unphysical region, with $\beta<-1$, they lie along the negative real axis, with $z_{\rm 2i}$ lying between $\beta$ and 0. A graphical representation of these results in the complex $z$ plane is shown in Fig. 3, with the dots representing the “inside” points $z_{\rm 1i}$ and $z_{\rm 2i}$, and the crosses representing the “outside” points $z_{\rm 1o}$ and $z_{\rm 2o}$.

For more general form factors, we restrict ourselves in the first instance to those functions with poles in the timelike half of the complex $Q_{i}^{2}$ plane. In this case it can be shown that all poles in $z_{\rm 1i}$ will be clustered around the positive real axis, with $|z_{\rm 1i}|\leq 1$ and $\operatorname{Re}z_{\rm 1i}>0$. Recall that for every pole $z_{1}$ there is a corresponding pole $z_{2}=\beta z_{1}$. In the physical region, all poles in $z_{\rm 2i}$ will therefore be clustered around the line at angle $\theta$, with $|z_{\rm 2i}|\leq 1$. In the unphysical region, the $z_{\rm 2i}$ poles must therefore satisfy $|z_{\rm 2i}|<|\beta|$ and $\operatorname{Re}z_{\rm 2i}<0$.

By Cauchy’s theorem, a closed contour integral is equivalent to summing the residue at the poles of the integrand enclosed by the contour. In continuing the integral from the physical to the unphysical region, the contour $\Gamma$ must be deformed in such a way that no poles cross the boundary defined by the contour. There are many possible choices of an appropriate contour in the unphysical region. However, based on the symmetry of the poles in $z_{1}$ and $z_{2}$, as discussed above, the closed contour shown in Fig. 3(b) merits special consideration. It consists of two semicircles connected by two straight line segments, which we denote as $\Gamma_{\rm CD}$. It clearly satisfies the criterion that no poles $z_{\rm 2i}$ or $z_{\rm 2o}$ cross the boundary of the contour, provided that the form factors only have poles in the timelike region of the $Q_{i}^{2}$ plane. Moreover, because of the symmetry of the TPE amplitudes under the interchange of $Q_{1}^{2}$ and $Q_{2}^{2}$, the line integrals along the contours $\Gamma_{1}$ and $\Gamma_{3}$ are equal, as are the line integrals along the contours $\Gamma_{2}$ and $\Gamma_{4}$. Furthermore, the real part of the integrand is symmetric between the upper and lower half planes, while the imaginary part is antisymmetric (hence the imaginary part integrates to 0 in any closed contour). Thus one only needs to compute the real part of the line integral along the contour $\Gamma_{1}$ in the upper half plane, plus the contribution from $\Gamma_{2}$,

[TABLE]

Form factor parametrizations $G(Q_{i}^{2})$ are driven almost exclusively by fits to data at spacelike values of $Q_{i}^{2}$. Aside from requiring no poles along the positive real axis, typically there is little or no consideration given to the location of poles in the complex $Q_{i}^{2}$ plane. Therefore, the appearance of poles in the spacelike half of $Q_{i}^{2}$ should be regarded as a nuisance rather than a reflection of any underlying physics. Nevertheless, there are form factor parametrizations in the literature that do have such poles. For example, the ratio of polynomials,

[TABLE]

is commonly used Kelly (2004); Arrington et al. (2007); Venkat et al. (2011). Requiring all $b_{k}\geq 0$ is sufficient to eliminate zeros in the denominator for positive real $Q_{i}^{2}$, but zeros can still occur in complex conjugate pairs with $\operatorname{Re}Q_{i}^{2}>0$.

To handle such cases we modify the contour of integration in the unphysical region, $\Gamma_{\rm CD}$. We choose instead a circular contour centered at $z_{c}=\textstyle{\frac{1}{2}}(1+\beta)$ of radius $\textstyle{\frac{1}{2}}(1-\beta)$, as illustrated in Fig. 4(a). The new contour is shown in red, while the $\Gamma_{\rm CD}$ contour is in light gray. As we require no poles to cross the boundary defined by the new contour, there are two possible failures. Firstly, the poles $z_{\rm 1i}$ are no longer restricted to the positive half plane, and can lie anywhere inside a circle of unit radius. This gives the possibility of a pole $z_{\rm 1o}$ crossing into the interior of the new contour, since the only restriction is $|z_{\rm 1o}|>1$. Secondly, the corresponding poles $z_{\rm 2i}$ are no longer restricted to the negative half plane, and can lie outside the new contour, since the only restriction here is $|z_{\rm 2i}|<|\beta|$. By careful analysis of the location of the poles for arbitrary values of $\beta$ and $Q^{2}$, we have derived the following condition:

The validity of the contour defined by a circle of radius $\textstyle{\frac{1}{2}}(1-\beta)$, centered at $z_{c}=\textstyle{\frac{1}{2}}(1+\beta)$, is that the poles $Q_{p}^{2}$ in the form factor $G(Q_{i}^{2})$ must satisfy the condition

[TABLE]

The condition (52) is equivalent to

[TABLE]

This condition is satisfied for all $\cos\eta<0$ (poles in the timelike half of $Q_{i}^{2}$), and for all $R>1$. We can visualize the condition using a circle of radius $\textstyle{\frac{1}{2}}Q^{2}$, centered at $\{\textstyle{\frac{1}{2}}Q^{2},0\}$, in the complex $Q_{i}^{2}$ plane, where the poles $Q_{p}^{2}$ in the form factors must lie outside this circle. An example using $Q^{2}=4$ GeV2 is shown in Fig. 4(b). If a form factor has poles in the spacelike half ($\cos\eta>0$), then the maximum $Q^{2}$ for which we can use the new contour is given by $Q_{\rm max}^{2}=|Q_{p}|^{2}/\cos\eta$.

III Impact of two-photon exchange

In this section we present the results of the numerical calculations of the TPE contributions to elastic $ep$ scattering within the dispersive approach, and discuss the differences with the traditional loop calculations with off-shell intermediate states. Including the nucleon and $\Delta$ intermediate state contributions, we compare the results with measurements of $e^{+}p$ to $e^{-}p$ cross section ratios from the recent CLAS Rimal et al. (2016), VEPP-3 Rachek et al. (2015); Nikolenko et al. (2014) and OLYMPUS Henderson et al. (1969) experiments, as well as with polarization observables sensitive to effects beyond the Born approximation Meziane et al. (2011).

III.1 Nucleon intermediate state

For the contributions to the box diagram from nucleon intermediate states, the inputs into the calculation are the proton elastic electric and magnetic form factors, $G_{E}^{p}$ and $G_{M}^{p}$. In the numerical calculations in this analysis we use the recent fit by Venkat et al. Venkat et al. (2011), which has the form of Eq. (51) with $N=3$ for both $G_{E}^{p}$ and $G_{M}^{p}$. As Fig. 5 illustrates, this fit gives similar results to other parametrizations, such as the ones by Kelly Kelly (2004) and Arrington et al. (AMT) Arrington et al. (2007). In contrast, the recent parametrization by Bernauer et al. Bernauer et al. (2010), which is based on a spline with 8 knots, displays distinctive wiggles at low $Q^{2}$ for both the electric and magnetic form factors. In Ref. Bernauer et al. (2010) a number of other functional forms were considered in fits to the world’s elastic electron–proton scattering data.

In the application of the numerical contour integration method to the calculation of the box diagram, care must be taken to ensure that all relevant poles are included inside the contour, following the condition given in Eq. (53). For some commonly used fits in the literature, such as those by Bosted Bosted (1995) or Brash Brash et al. (2002) (in which the denominators of $G_{E}^{p}$ and $G_{M}^{p}$ are fifth-order polynomials in $\sqrt{Q^{2}}$), or the Kelly parametrization Kelly (2002) (third-order polynomials in $Q^{2}$), there are no upper limits on $Q^{2}$, since all poles occur in the timelike region. On the other hand, for the AMT Arrington et al. (2007) and Venkat et al. Venkat et al. (2011) parametrizations, which involve denominators with 5th-order polynomials in $Q^{2}$, poles in the spacelike region limit the range of applicability to $Q^{2}<4.5$ GeV2. For the purposes of the data analysis in this paper, all the above parametrizations are valid; however, for future applications at higher $Q^{2}$ values care must be taken to ensure the chosen parametrization has a suitable pole structure.

The TPE contributions from nucleon intermediate states to the imaginary parts of the generalized $F_{1}^{\prime}$, $F_{2}^{\prime}$ and $G_{a}^{\prime}$ form factors are illustrated in Fig. 6 as a function of energy $E$ (in the lab frame), at a representative $Q^{2}$ value of 3 GeV2. (The results at other $Q^{2}$ values are qualitatively similar.) In the low energy region, $E\lesssim 0.1$ GeV, in Fig. 6(a) the TPE amplitudes display a logarithmic divergence in $E$. Although the imaginary parts diverge, the dispersive integrals (29) for the real parts remain finite. To accommodate this in our numerical analysis, we fit the low-$E$ expressions for $\operatorname{Im}F_{1}^{\prime}$, $\operatorname{Im}F_{2}^{\prime}$ and $\operatorname{Im}G_{a}^{\prime}$ to functions of the form

[TABLE]

with the parameters $a$ and $b$ determined by a least-squares fit for $E$ between $0.001$ and $0.01$ GeV. At high energies, the imaginary parts of the $F_{1}^{\prime}$ and $F_{2}^{\prime}$ form factors become constant, as is apparent for $E\gtrsim 6$ GeV from Fig. 6(b), where the magnified scale more clearly illustrates the asymptotic behavior. The imaginary part of the axial $G_{a}^{\prime}$ form factor falls off as $1/E$ for $E\to\infty$. This high energy behavior is sufficient to ensure the convergence of the dispersive integrals of Eq. (29). Note that the corrections to the form factors are relative to the Maximon-Tjon result for the infrared part of the TPE Maximon and Tjon (2000).

For our numerical calculation, we compute the imaginary part of the TPE amplitudes on a logarithmic grid of 51 points in $E$, ranging up to 100 GeV. In the unphysical region, the two-dimensional numerical contour integral is computed using either of the contours discussed in the previous section, as appropriate to the poles of the form factors. In the physical region, the numerical contour integral is on the unit circle (although a direct numerical integration of (43) using Eq. (36) can also be used). We then interpolate between the grid points with a spline fit to obtain a continuous function of $E$. To obtain the real part on a grid of 20 equally spaced points in $\varepsilon$, the dispersion integral is evaluated using the $\log E$ fit for $E<0.01$ GeV, the spline fit for $0.01<E<100$ GeV, and an extrapolation beyond 100 GeV using the known asymptotic behavior. This approach can be tested against the analytic results obtained in the previous section, as well as the known analytic results for $e\mu$ scattering.

Interestingly, for the real part of the TPE amplitudes the dispersive integral is dominated by contributions from the unphysical region, $E<M(\tau+\sqrt{\tau(1+\tau)})$, which for $Q^{2}=3$ GeV2 is $E=1.97$ GeV. This results in the generally smoothly decaying functions for $E\gtrsim 2$ GeV observed in Fig. 7. The real parts of each of the form factors are negative in the region illustrated, with the $F_{1}^{\prime}$ form factor having the largest magnitude, and the $G_{a}^{\prime}$ form factor the smallest magnitude. Compared with the direct loop calculation in terms of the off-shell nucleon intermediate states, differences arise for the Pauli $F_{2}$ form factor term, which translates to a difference for the $G_{a}^{\prime}$ form factor. The results for the $F_{1}^{\prime}$ and $F_{2}^{\prime}$ form factors are the same for both calculations, as was observed previously in Refs. Borisyuk and Kobushkin (2008); Tomalak and Vanderhaeghen (2015). Numerically, the differences are relatively small, however, as Fig. 6 indicates, becoming notable only for $E\gtrsim 3$–4 GeV.

The total TPE correction (27) to the elastic cross section from nucleon intermediate states relative to the Mo-Tsai infrared prescription Mo and Tsai (1969), denoted $\delta_{N}$, is shown in Fig. 8 as a function of $\varepsilon$ for $Q^{2}$ between 0.5 GeV2 and 3 GeV2. As found in previous loop calculations Blunden et al. (2003, 2005), the corrections at these kinematics are negative, and increase in magnitude with increasing $Q^{2}$. For large $Q^{2}$ values we expect the reliability of the hadronic calculation to deteriorate, but indications from earlier work Blunden et al. (2003, 2005) suggest that it remains sizeable. In practice, since only the $G_{a}^{\prime}$ form factor is affected, and its magnitude is considerably smaller than that of $F_{1}^{\prime}$ and $F_{2}^{\prime}$, as illustrated in Fig. 7, the off-shell effects play a relatively minor role in $\delta_{N}$, with the dispersive and loop results almost indistinguishable.

III.2 $\Delta$ intermediate state

The contribution to the TPE amplitude from intermediate states involving the spin-$\nicefrac{{3}}{{2}}$, isospin-$\nicefrac{{3}}{{2}}$ $\Delta$ baryons is computed from the $\gamma^{*}N\to\Delta$ electromagnetic transition operator, $\Gamma_{\gamma N\to\Delta}^{\alpha\mu}$. Usually this is expressed in terms of three Jones-Scadron transition form factors, $G_{M}^{*}(Q^{2})$, $G_{E}^{*}(Q^{2})$ and $G_{C}^{*}(Q^{2})$, corresponding to magnetic, electric, and Coulomb multipole excitations, respectively Jones and Scadron (1973). Although the $\gamma^{*}N\to\Delta$ cross section is diagonal in these functions, they are cumbersome to work with in the transition vertex function, and other parametrizations have also been suggested in the literature Pascalutsa et al. (2007); Kondratyuk et al. (2005); Nagata et al. (2009); Pascalutsa and Tjon (2004); Caia et al. (2004). In this work we follow Ref. Kondratyuk et al. (2005) and use the on-shell equivalent parametrization of the $\gamma^{*}N\to\Delta$ vertex

[TABLE]

where $p_{\Delta}$ and $q$ are the momenta of the outgoing $\Delta$ and incoming photon, respectively. The $g_{i}\ (i=1,2,3)$ transition functions are related to the Jones-Scadron form factors by

[TABLE]

where

[TABLE]

Since in practice the magnetic multipole dominates the $\gamma^{*}N\to\Delta$ transition, the $g_{1}$ function is determined mostly by $G_{M}^{*}$. The electric form factor $G_{E}^{*}$ determines the difference $g_{2}-g_{1}$, while $g_{3}$ is sensitive to $G_{E}^{*}$ and the Coulomb form factor $G_{C}^{*}$.

In the present analysis, we take the Jones-Scadron form factors from the phenomenological parametrization by Aznauryan Aznauryan and Burkert (2012); Aznauryan (2016),

[TABLE]

where

[TABLE]

with $Q^{2}$ in units of GeV2. Empirical fits to data suggest that the E1/M1 multipole ratio $R_{\rm EM}$ and the S1/M1 multipole ratio $R_{\rm SM}$ can be well approximated by

[TABLE]

As this parametrization only has poles for timelike $Q^{2}$, we can use the contour $\Gamma_{\rm CD}$ of Eq. (50) in the unphysical region.

The $\gamma^{*}N\Delta$ transition form factors $g_{1}$, $g_{2}$ and $g_{3}$ from the Eqs. (58)–(60) are illustrated in Fig. 9(a) as a function of $Q^{2}$. At moderate and large $Q^{2}$ values, $Q^{2}\gtrsim 0.1$ GeV2, the $g_{1}$ and $g_{2}$ form factors dominate, with the $g_{3}$ form factor essentially zero. At very low $Q^{2}\lesssim 0.02$ GeV2, the $g_{3}$ form factor rises rapidly and becomes larger than the largest contribution; for the phenomenological applications relevant to this paper, however, its role is essentially negligible.

The quality of the fit to the magnetic transition form factor $G_{M}^{*}$ is shown in Fig. 9(b), compared with data from several experiments Ungaro et al. (2006); Frolov et al. (1999); Joo et al. (2002) for $Q^{2}$ up to $\approx 6$ GeV2. The parametrization in Eq. (58a) is compared with an alternative parametrization from the recent analysis by Zhou and Yang Zhou and Yang (2015), which agrees with the data in the intermediate $Q^{2}$ region, $Q^{2}\sim 3$–4 GeV2, but underestimates (overestimates) the data at lower (higher) $Q^{2}$ values. Similarly, good agreement is obtained for the $R_{\rm EM}$ and $R_{\rm SM}$ ratios in Figs. 9(c) and 9(d), respectively, for the parametrizations in Eqs. (60) over the full range of $Q^{2}$ ($\lesssim 6$ GeV2) where data are available. The parametrization Zhou and Yang (2015) also agrees with the data at low $Q^{2}$ values, $Q^{2}\lesssim 1$ GeV2, but discrepancies appear for larger $Q^{2}$.

Using the Aznauryan parametrization Aznauryan and Burkert (2012); Aznauryan (2016) of the $\gamma^{*}N\Delta$ form factors, the TPE contributions from $\Delta$ intermediate states to the generalized form factors $F_{1}^{\prime}$, $F_{2}^{\prime}$ and $G_{a}^{\prime}$ are illustrated in Fig. 10 for both the imaginary and real parts. For the imaginary parts of the amplitudes, Fig. 10(a) shows resonance-like structure appearing in the unphysical region, $0.34<E<2$ GeV. As for the nucleon case, the unphysical region accounts for most of the dispersive integral, giving rise to smoothly decaying real parts of the amplitudes for $E\gtrsim 3$ GeV, as Fig. 10(b) illustrates.

As was observed in previous calculations of loop corrections Kondratyuk et al. (2005); Drell and Fubini (1959); Campbell (1969); Greenhut (1969), the contribution from $\Delta$ intermediate states, $\delta_{\Delta}$, is generally of opposite sign to the nucleon contribution $\delta_{N}$ for $Q^{2}\gtrsim 1$ GeV2, and increases in magnitude with increasing $Q^{2}$, as Fig. 11 illustrates. Unlike previous loop calculations Kondratyuk et al. (2005); Zhou and Yang (2015), however, the correction $\delta_{\Delta}$ in the dispersive approach is well-behaved for all $\varepsilon$, vanishing in the $\varepsilon\to 1$ limit. In contrast, the correction $\delta_{\Delta}$ from the loop calculation with off-shell $\Delta$ states diverges as $\varepsilon\to 1$, as illustrated in Fig. 12(a). Note that the imaginary parts of the TPE amplitudes are identical for both calculations — only the real parts differ.

The nature of this divergence can be seen by plotting $\delta_{\Delta}$ versus electron energy $E$ instead of $\varepsilon$, as in Fig. 12(b). The linear divergence in $E$ indicates a violation of the Froissart bound Froissart (1961), and the breakdown of unitarity. This “pathological” behavior is not due to the inapplicability of hadronic models when $E\to\infty$ and $\varepsilon\to 1$, as suggested in Ref. Zhou and Yang (2015). Rather, it arises from an unphysical behavior of the off-shell contributions at high energies for interactions with derivative couplings.

As Fig. 12 illustrates, the loop calculations at large $\varepsilon$ are actually very sensitive to the shape of the $\gamma^{*}N\Delta$ form factors employed. Whether one use a dipole approximation or a more realistic parametrization, leads to significant differences with the dispersive approach already for $E\gtrsim 3$–4 GeV. (Here, for simplicity only the magnetic contribution to $\delta_{\Delta}$ is shown, but the effects are similar for the other $\gamma^{*}N\Delta$ form factors also.) In fact, the differences between the dispersive and loop calculations are significant not just near $\varepsilon\approx 1$, but also at lower $\varepsilon$ values. Generally, the magnitude of the dispersive $\Delta$ corrections is smaller than the loop results, resulting in less cancellation with the intermediate state nucleon contribution.

III.3 Verification of TPE effects

Having detailed the calculation of the TPE corrections from the nucleon and $\Delta$ intermediate states, we next compare the role of these corrections in observables that are particularly sensitive to effects beyond the Born approximation. These include the ratio of unpolarized $e^{+}p$ to $e^{-}p$ elastic scattering cross sections, and polarization transfer cross sections for longitudinally and transversely polarized electrons and protons.

III.3.1 $e^{+}p$ to $e^{-}p$ ratio

One of the observables that is most sensitive to the effects of TPE is the ratio of $e^{+}p$ to $e^{-}p$ elastic cross sections, which in the one-photon exchange approximation is unity. Since the TPE terms enter the $e^{+}p$ cross section with opposite sign to that in the $e^{-}p$ reaction, the ratio

[TABLE]

where $\sigma^{e^{\pm}}\equiv d\sigma({e^{\pm}p\to e^{\pm}p})/d\Omega$, provides a direct measure of effects beyond the Born approximation. Earlier data from elastic $e^{+}p$ and $e^{-}p$ experiments in the 1960s from SLAC Browman et al. (1965); Mar et al. (1968), Cornell Anderson et al. (1966), DESY Bartel et al. (1967) and Orsay Bouquet et al. (1968) gave some hints of a small enhancement of $R_{2\gamma}$ at forward angles and low $Q^{2}$, but were in the region (at large $\varepsilon$) where TPE is relatively small and were consistent within errors with $R_{2\gamma}=1$.

More recently, several dedicated $e^{+}p$ to $e^{-}p$ ratio experiments have been performed in CLAS at Jefferson Lab Rimal et al. (2016), VEPP-3 in Novosibirsk Rachek et al. (2015); Nikolenko et al. (2014) and OLYMPUS at DESY Henderson et al. (1969) aimed at providing measurements of $R_{2\gamma}$ over a larger range of $\varepsilon$ and $Q^{2}$ with significantly reduced uncertainties. In Fig. 13 the $R_{2\gamma}$ ratio from the CLAS experiment is shown as a function of $\varepsilon$ at averaged $Q^{2}$ values of $\langle Q^{2}\rangle=0.85$ GeV2 and $\langle Q^{2}\rangle=1.45$ GeV2 [Figs. 13(a) and (b), respectively], and as a function of $Q^{2}$ at averaged $\varepsilon$ values of $\langle\varepsilon\rangle=0.45$ and $\langle\varepsilon\rangle=0.88$ GeV2 [Figs. 13(c) and (d), respectively]. Most of the data at the larger $\varepsilon$ values are consistent with unity within the errors, but suggest a nonzero ratio, $\approx 2\%$ – 4% greater than unity, at the lowest $\varepsilon$ value for the higher-$Q^{2}$ set. The trend is consistent with the ratio calculated here, which shows a rising $R_{2\gamma}$ with decreasing $\varepsilon$. At these kinematics the calculated TPE correction is dominated by the nucleon elastic intermediate state, with the $\Delta$ contribution reducing the ratio slightly. Note that both the data and the calculated TPE corrections here (and elsewhere, unless otherwise stated) are shown relative to the Mo-Tsai infrared result.

The same trend is seen when the $R_{2\gamma}$ data are viewed as a function of $Q^{2}$ for fixed $\varepsilon$. At the larger average $\varepsilon$ value, $\langle\varepsilon\rangle=0.88$, the effects are consistent with zero as well as with the small predicted TPE correction. At the smaller value $\langle\varepsilon\rangle=0.45$, on the other hand, the larger predicted effect is consistent with the larger $R_{2\gamma}$ values with increasing $Q^{2}$. Again the effects of the $\Delta$ intermediate state are small at low $Q^{2}$ values, but become visible at larger $Q^{2}$, where they improve the agreement between the theory and experiment.

Data from the VEPP-3 experiment at Novosibirsk Rachek et al. (2015); Nikolenko et al. (2014), taken at energies $E\approx 1$ GeV and 1.6 GeV, are shown in Fig. 14 as a function of $\varepsilon$. The data correspond to a $Q^{2}$ range between $\approx 0.3$ GeV2 and $\approx 1.5$ GeV2, with $\varepsilon$ down to $\approx 0.3$. The ratio at the low $\varepsilon$ values shows an effect of magnitude 1% – 2%, slightly below but still consistent with the calculated TPE result at the $\approx 1\sigma$ level.

Most recently, the OLYMPUS experiment at DESY Henderson et al. (1969) measured the ratio $R_{2\gamma}$ at an energy $E=2.01$ GeV over a large range of $\varepsilon\sim 0.45-1$, corresponding to a $Q^{2}$ range from $\approx 0.2$ GeV2 to 2 GeV2. The results for the ratio $R_{2\gamma}$ are shown in Fig. 15. Interestingly, in contrast to the results from the CLAS and VEPP-3 experiments in Figs. 13 and 14, at large $\varepsilon$ values the trend in the data is towards values of the ratio slightly below unity, whereas the calculated dispersive TPE corrections give a ratio that has a small, $\lesssim 1\%$ enhancement above unity. At the lower $\varepsilon$ values, the trend is toward increasing values of $R_{2\gamma}$, consistent with the TPE calculation. Within the statistical and systematic uncertainties, including the overall normalization uncertainty of the OLYMPUS data, the theoretical result is consistent with the data over the entire $\varepsilon$ range. Note also that the $\lesssim 0.5\%$ correlated systematic (normalization) uncertainty quoted for the OLYMPUS data Henderson et al. (1969) is somewhat smaller than in the other $e^{+}p/e^{-}p$ experiments Rachek et al. (2015); Nikolenko et al. (2014).

While the possibility of unexpected effects in the high-$\varepsilon$ region is intriguing, we should note that the ratio $R_{2\gamma}$ defined in Ref. Henderson et al. (1969) is normalized by a ratio of $e^{+}$ to $e^{-}$ events obtained from a Monte Carlo (MC) simulation of the experiment designed to account for differences between electrons and positrons, $R_{2\gamma}\to(N_{\rm exp}^{e^{+}}/N_{\rm exp}^{e^{-}})/(N_{\rm MC}^{e^{+}}/N_{\rm MC}^{e^{-}})$. Here $N_{\rm exp}^{e^{\pm}}$ is the number of observed $e^{\pm}$ events, and $N_{\rm MC}^{e^{\pm}}$ is the number of simulated $e^{\pm}$ counts, taking into account radiative effects and various experimental settings. Of course, in order to simulate the elastic scattering cross sections, some input about the $e^{\pm}p$ interaction is needed for the MC, and it is possible that this may introduce additional model dependence into the procedure. Indeed, simulations using radiative corrections computed to order $\alpha^{3}$ versus those computed to all orders through exponentiation show that the latter can give $R_{2\gamma}$ values as much as 1% higher at the lowest $\varepsilon$ points Henderson et al. (1969). In Fig. 15 the results shown correspond to the ratio $R_{2\gamma}$ extracted with radiative corrections computed to all orders in $\alpha$.

The relatively large overall uncertainties on all of the currently available $R_{2\gamma}$ data unfortunately precludes any definitive conclusions about TPE effects that can be reached, other than that the effects are generally consistent with zero, as well as with the signs and magnitudes expected from the dispersive TPE calculations. This scenario calls for an urgent need for new measurements of $e^{+}p$ to $e^{-}p$ ratios at large $Q^{2}$, $Q^{2}\gtrsim 2$ GeV2, and over a range of $\varepsilon$ values below $\varepsilon\sim 0.5$, where the TPE effects are predicted to be large enough ($\sim 2\%$) to be more clearly identified experimentally. On the other hand, the negative values of the slope in $\varepsilon$ predicted by the TPE calculations are generally consistent with the data from each of the CLAS Rimal et al. (2016), VEPP-3 Rachek et al. (2015); Nikolenko et al. (2014) and OLYMPUS Henderson et al. (1969) $e^{+}p/e^{-}p$ experiments.

III.3.2 Polarization observables

A complementary set of observables that can provide information on TPE effects involves polarization transfer in the elastic scattering of longitudinally polarized electrons from (unpolarized) protons, with measurement of the polarization of the final state proton, $\vec{e}p\to e\vec{p}$. Defining $P_{T}$ and $P_{L}$ to be the polarizations of recoil protons in the transverse and longitudinal directions relative to the proton momentum in the scattering plane, one has Guichon and Vanderhaeghen (2003)

[TABLE]

where $\sigma_{R}=\sigma_{R}^{\rm Born}(1+\delta_{{\gamma\gamma}})$, with the reduced Born cross section $\sigma_{R}^{\rm Born}$ given in Eq. (8), and $\delta_{{\gamma\gamma}}$ in Eq. (27). Note that the IR subtractions we have made in $F_{1}^{\prime}$ and $F_{2}^{\prime}$ are such that terms in the numerator and denominator of $P_{L}$ and $P_{T}$ cancel exactly, independent of regularization scheme. Taking the ratio of the transverse to longitudinal polarizations, we define

[TABLE]

where $\mu_{p}$ is the proton’s magnetic moment. In the Born approximation, this reduces to a simple ratio of the electric to magnetic form factors, $R_{TL}\to\mu_{p}G_{E}/G_{M}$, which is a function only of $Q^{2}$ and is independent of $\varepsilon$.

The GEp2$\gamma$ experiment at Jefferson Lab Meziane et al. (2011) measured the ratios $R_{TL}$ and $P_{L}/P_{L}^{(0)}$, where $P_{L}^{(0)}$ is the Born level longitudinal polarization, at several values of $\varepsilon$ for fixed $Q^{2}=2.49$ GeV2. These are shown in Fig. 16 as a function of $\varepsilon$, compared with the dispersive TPE calculations including nucleon and $\Delta$ intermediate states. The TPE effect on the longitudinal polarization ratio is very small, with the $P_{L}/P_{L}^{(0)}$ ratio only marginally below unity for all $\varepsilon$ values. Although the trend of the data suggests an increasing effect at high $\varepsilon$, the data are consistent with the TPE calculation if (correlated and uncorrelated) systematic uncertainties are taken into account.

The $\varepsilon$ dependence of the $R_{TL}$ ratio in Fig. 16 is also very weak, and in good agreement with the dispersive TPE calculation, especially once the $\Delta$ intermediate states are included. As for the $R_{2\gamma}$ ratio, higher-precision measurements of the polarization observables at larger $Q^{2}$ and lower $\varepsilon$ values would be valuable in more definitively identifying effects beyond the Born approximation.

IV Outlook

In this paper we have presented a new approach to evaluating two-photon exchange effects in elastic electron–proton scattering, based on a dispersion relation analysis of the scattering amplitudes. We considered two methods for evaluating the imaginary parts of the loop diagrams, using analytic and numerical methods, and including the contributions from nucleon and $\Delta$ resonance intermediate states.

In contrast to previous calculations based on the direct evaluation of loop integrals, the dispersive analysis avoids uncertainties associated with off-shell intermediate states, and guarantees the correct behavior in the high energy limit. This problem is particularly egregious for the case of derivative interactions, such as for the $\Delta$ baryon, where the TPE amplitude in the off-shell calculation diverges linearly with energy in the forward limit. The dispersive approach, on the other hand, respects unitarity and is well-behaved at all energies.

The analytic dispersive method, which has been used recently in the literature Borisyuk and Kobushkin (2008, 2012, 2014, 2015); Tomalak and Vanderhaeghen (2015), has the advantage of allowing closed analytic expressions for the amplitudes in terms of simple logarithms, provided the vertex form factors can be parametrized in terms of sums or products of monopole functions. This has the virtue of increased speed of computation, but is limited by the accuracy of the monopole parametrization of the proton’s electric and magnetic form factors, which typically deteriorates markedly for $Q^{2}\gtrsim 5$–6 GeV2.

The numerical contour method, in contrast, allows for a wide range of form factor parametrizations, and is relatively straightforward to implement. We find that in practice a simple contour is valid for any parametrization which has poles in the timelike region of $Q^{2}$. For more elaborate parametrizations that have with poles in the spacelike region, a special choice of contour can be used up to some maximum $Q^{2}$ that depends on the exact location of the poles.

To verify the utility of the dispersive approach, we have compared the results of the numerical TPE calculations with the most recent data on the ratio of $e^{+}p$ to $e^{-}p$ elastic scattering cross sections, which is directly sensitive to electromagnetic effects beyond the Born approximation. We find good agreement with the data from the CLAS Rimal et al. (2016) and VEPP-3 Rachek et al. (2015); Nikolenko et al. (2014) experiments. The results are also consistent, within the experimental uncertainties, with the more recent OLYMPUS experiment Henderson et al. (1969), which suggests a trend in the opposite direction at near-forward angles compared with the TPE calculations and the other data sets Rimal et al. (2016); Rachek et al. (2015); Nikolenko et al. (2014).

For the future, it will be important to extend the present framework to inelastic non-resonant intermediate states, including the $\pi N$ continuum, and allowing for widths of resonances. Efforts in this direction have been made recently in the literature Borisyuk and Kobushkin (2015); Tomalak et al. (2016), and will be aided by better knowledge of the empirical virtual Compton scattering amplitudes at non-forward angles. Beyond this, a longer term challenge is to further generalize the dispersive approach to higher $Q^{2}$ and intermediate state masses, merging the phenomenological hadronic description with one that expresses the TPE amplitudes explicitly in terms of partonic degrees of freedom.

On the experimental front, new, higher precision data on $e^{+}p$ and $e^{-}p$ cross sections at larger $Q^{2}$ values, $Q^{2}\gtrsim 3$ GeV2, and lower $\varepsilon$ values, $\varepsilon\ll 1$, are needed to unambiguously identify TPE effects directly. The possibility of achieving this with a dedicated positron source at the 12 GeV Jefferson Lab facility remains an exciting prospect Arrington (2009). The technology described here for the TPE calculations can also be readily applied to the evaluation of $\gamma Z$ interference contributions in parity-violating electron–proton scattering, in the extraction of the strange electroweak form factors and the weak charge of the proton Arrington et al. (2011); Afanasev and Carlson (2005); Gorchtein and Horowitz (2009); Gorchtein et al. (2011); Rislow and Carlson (2011); Tjon et al. (2009); Sibirtsev et al. (2010); Hall et al. (2013, 2014, 2016).

Appendix A Generalized form factors

Beyond the Born approximation, the total amplitude for the elastic $ep\to ep$ scattering process can, for a massless electron, be decomposed into three independent amplitudes, or generalized form factors Guichon and Vanderhaeghen (2003). In this appendix we describe a practical method for extracting these generalized form factors from the total amplitude.

The objective is to map ${\cal M}_{{\gamma\gamma}}$ of Eq. (10) onto the generalized matrix element $\widehat{\cal M}_{\rm gen}(F_{1}^{\prime},F_{2}^{\prime},G_{a}^{\prime})$ given by Eq. (25),

[TABLE]

where the generalized current operator $\widehat{\Gamma}^{\mu}$ is defined as

[TABLE]

and the generalized form factors are functions of two variables, taken to be $Q^{2}$ and $\nu$.

The basic idea behind this method is to project $\widehat{\cal M}_{\rm gen}$ onto three linearly independent quantities (pseudo-observables) using Dirac trace techniques. The pseudo-observable projections are linear combinations of the amplitudes $F_{1}^{\prime}$, $F_{2}^{\prime}$, and $G_{a}^{\prime}$. The same projections are also made for ${\cal M}_{{\gamma\gamma}}$. One can then invert the transformation matrix to obtain $F_{1}^{\prime}$, $F_{2}^{\prime}$, and $G_{a}^{\prime}$ in terms of the pseudo-observable projections of ${\cal M}_{{\gamma\gamma}}$.

While any three linearly independent projections will suffice, it is convenient to use the same functional form given by $\widehat{\cal M}_{\rm gen}$. Consider the pseudo-observable

[TABLE]

where $x$, $y$, and $z$ are linear combinations of $F_{1}^{\prime}$, $F_{2}^{\prime}$, and $G_{a}^{\prime}$, and $A$, $B$, and $C$ are placeholder coefficients representing the three independent projections. The overall factor $(-e^{2}/q^{2})^{2}$ is irrelevant for this derivation, and can be absorbed into the coefficients $\{A,B,C\}$. The expression for $\sigma_{\rm gen}$ can be obtained from

[TABLE]

where the leptonic and hadronic tensors are given by

[TABLE]

respectively. Note that rather than keeping the vector and axial-vector terms separate, we have combined them into one compact expression. This is possible because the parity-violating terms in the combined amplitude vanish after taking the traces. Evaluating the traces, one finds

[TABLE]

where the transformation matrix $\mathbb{M}$ is given by

[TABLE]

Inverting the matrix $\mathbb{M}$ will obtain the relationships of interest,

[TABLE]

The same projections as in Eq. (66) can now be made for ${\cal M}_{{\gamma\gamma}}$,

[TABLE]

The set of projected functions $\left\{x,y,z\right\}$ are then combined to give the TPE amplitudes using Eq. (81). There is an apparent kinematic singularity in Eq. (81) at $\nu^{2}=\tau(\tau+1)=\nu_{\rm ph}^{2}$, which is the threshold between the physical and unphysical regions. However, the full expressions for the imaginary parts of $F_{1}^{\prime}(Q^{2},\nu)$, $F_{2}^{\prime}(Q^{2},\nu)$, and $G_{a}^{\prime}(Q^{2},\nu)$ are continuous, smooth, and finite across this boundary. Nevertheless, for numerical work we avoid directly using this kinematic point.

From Eq. (10), the pseudo-observable $\sigma_{\rm gen}$ has the form

[TABLE]

where the leptonic and hadronic tensors of Eq. (11)–(12) are

[TABLE]

The decomposition of the total amplitude into a basis of generalized form factors is not unique. Another convention in the literature Chen et al. (2004); Tomalak and Vanderhaeghen (2015) is to use the generalized matrix element

[TABLE]

The relationship between the set $\{F_{1}^{\prime},F_{2}^{\prime},G_{a}^{\prime}\}$ and the set $\{\widetilde{F}_{1},\widetilde{F}_{2},\widetilde{F}_{3}\}$ is $F_{1}^{\prime}=\widetilde{F}_{1}+\nu\widetilde{F}_{3}$, $F_{2}^{\prime}=\widetilde{F}_{2}$, and $G_{a}^{\prime}=-\tau\widetilde{F}_{3}$ Chen et al. (2004). This relationship can be easily derived by contracting $\widetilde{\cal M}_{\rm gen}$ with $\widehat{\cal M}_{\rm gen}^{\dagger}$, in analogy with Eq. (66), and comparing the transformation matrices.

Appendix B Analytic expressions for imaginary parts of Passarino-Veltman

functions

In the notation of the LoopTools package Hahn and Perez-Victoria (1999), the full dependence of the $s$-channel PV functions on kinematic variables is

[TABLE]

Following the discussion in Sec. II.2.1, the full expressions for the imaginary parts of the these functions are

[TABLE]

where $s_{M}$ and $s_{W}$ are given by Eq. (34), the quantity $J$ is given by

[TABLE]

and we have introduced the dimensionless variable $z=Z/Y$, with

[TABLE]

The expression for $d_{0}(s;\Lambda_{1},\Lambda_{2})$ in Eq. (88c) has been rewritten in a form that is numerically stable for very small values of $\Lambda_{i}$. By comparison, the form given in Eq. (37), with $X^{2}=Y^{2}-Z^{2}$, is susceptible to roundoff error in this limit. We also note that while the real part of $D_{0}$ has a logarithmic dependence on $m_{e}$, the imaginary part does not, and is therefore finite in the limit $m_{e}\to 0$, which we have used throughout this paper.

The PV functions $c_{0}(s;\Lambda)$ and $d_{0}(s;\Lambda_{1},\Lambda_{2})$ are IR-divergent in the limit $\Lambda^{2}\to 0$. Replacing $\Lambda^{2}\to\lambda^{2}$, and keeping only logarithmic terms in $\lambda^{2}$, we have the three IR-divergent combinations

[TABLE]

One can use these expressions to explicitly show that there is no residual dependence on $\lambda$ in the total TPE amplitudes, after subtracting the model-independent Maximon and Tjon IR-divergent terms. For numerical calculations, it is more convenient to use the full expressions (88), with a value for $\lambda$ satisfying the criterion $\lambda^{2}\ll m_{e}^{2}$.

For sums of monopole form factors, the integrals of interest are sums of primitive integrals $I_{ij}$ given in Eq. (42), with $N=i+j$. Table 1 gives the values for $I_{ij}$ up to $N=2$, which suffices for nucleon intermediate states. In general, for other states like the $\Delta$ one needs up to $N=3$, which requires the form factor to behave like $1/Q_{i}^{4}$ asymptotically, but the procedure follows analogously.

Appendix C Form factor reparametrizations

Here we present the reparametrizations of the nucleon and $\Delta$ vertex form factors in terms of sums and/or products of monopoles, suitable for use in the analytic expressions of Sec. II.2.1. Fits are over the range $0<Q^{2}<8$ GeV2, with $F_{1}(Q^{2})$ being a 5-parameter fit, while all others are 4-parameter fits.

The nucleon form factors are fit to the parametrization of Ref. Venkat et al. (2011),

[TABLE]

with $Q^{2}$ in GeV2. Defining

[TABLE]

the $\Delta$ transition form factors are fit to the parametrization of Ref. Aznauryan and Burkert (2012); Aznauryan (2016),

[TABLE]

The form factors $g_{1}(Q^{2})$, $g_{2}(Q^{2})$, and $g_{3}(Q^{2})$ can be obtained from simple combinations of these parametrizations, while still allowing for implementation in analytic form.

The use of these reparametrizations in the numerical integration method allows for a test of our codes against the analytic results. We were routinely able to obtain a relative agreement at the level of five significant digits. Similarly, the validity of using these reparametrizations in the analytic codes can be tested against the numerical results using the original functional forms. In general, we found that the relative differences between the analytic and numerical evaluations of the imaginary parts of the TPE amplitudes were comparable to the relative differences between the original and reparametrized forms over the relevant range of $Q^{2}$. As the original and reparametrized vertex form factors given in this appendix agree at roughly the 2% level when averaged over the range $0<Q^{2}<8$ GeV2, we find a 2% agreement in $F_{1}^{\prime}(Q^{2},\nu)$, $F_{2}^{\prime}(Q^{2},\nu)$, and $G_{a}^{\prime}(Q^{2},\nu)$ up to $Q^{2}=4$ GeV2. This suggests that the TPE results are not very sensitive to pole structure of the vertex form factor parametrizations in the complex plane.

Acknowledgements

We thank J. Bernauer, D. Higinbotham and A. Schmidt for helpful discussions and communications, and I. Aznauryan for providing a parametrization of the $N\Delta$ transition form factors. PGB thanks Jefferson Lab for support during a sabbatical leave, where part of this work was completed. This work was supported by NSERC (Canada) and DOE Contract No. DE-AC05-06OR23177, under which Jefferson Science Associates, LLC operates Jefferson Lab.