Elements of QED-NRQED Effective Field Theory: II. Matching of Contact Interactions

In 2010 the first extraction of the proton charge radius from muonic hydrogen was found be be five standard deviations away form the regular hydrogen value. Eight years later, this proton radius puzzle is still unresolved. One of the most promising avenues to resolve the puzzle is by a muon-proton scattering experiment called MUSE. The typical momenta of the muons in this experiment are of the order of the muon mass. In this energy regime the muons are relativistic but the protons are non-relativistic. The interaction between them can be described by QED-NRQED effective field theory. In a previous paper we have shown how QED-NRQED reproduces Rosenbluth scattering up to $1/M^2$, where $M$ is the proton mass, and relativistic scattering off a static potential at ${\cal O}(Z^2\alpha^2)$ and leading power in $M$. In this paper we determine the Wilson coefficients of the four-fermion contact interactions at ${\cal O}(Z^2\alpha^2)$ and power $1/M^2$. Surprisingly, we find that the coefficient of the spin-independent interaction vanishes, implying that MUSE will be sensitive mostly to the proton charge radius and not spin-independent two-photon exchange effects.


Introduction
The response of an on-shell proton to a one-photon electromagnetic probe can be parametrized in terms of two non-perturbative form factors, the so-called "Dirac" (F 1 (q 2 )) and "Pauli" (F 2 (q 2 )) form factors, see (2) below. Alternatively, one can define a different linear combination the "electric" (G E = F 1 + q 2 F 2 /4M 2 , where M is the proton mass) and "magnetic" (G M = F 1 +F 2 ) form factors. The slope of G E at q 2 = 0 defines the proton charge radius (squared), as (r p E ) 2 = 6 dG p E (q 2 )/dq 2 | q 2 =0 . The two main methods to extract r p E are lepton-proton scattering and spectroscopy of lepton-proton bound states, i.e. regular or muonic hydrogen. Until 2010, only values from electron-proton scattering and regular hydrogen were available. In 2010, the first extraction of r p E from muonic hydrogen was reported as r p E = 0.84184(67) fm [1]. Surprisingly, it was five standard deviations lower than the regular hydrogen value of r p E = 0.8768(69) fm [2]. The most recent update from muonic hydrogen is r p E = 0.84087(39) fm [3], while the 2014 CODATA value is r p E = 0.8751(61) fm [4]. This discrepancy is known as the "proton radius puzzle".
The most exciting interpretation of the puzzle is as a new interaction that distinguishes electrons and muons, see, e.g., the references in [5,6]. In such explanations the crucial questions are what assumptions one is making about this new interaction, in particular how does it couple to hadrons besides protons. At the same time, one would like to test more mundane explanations such as experimental or theoretical issues in the extractions of r p E from spectroscopy or scattering using both electrons and muons 1 . We briefly review those. Regular hydrogen spectroscopy: No issues were identified with the theoretical input. On the experimental side, new measurements were published in the last two years with error bars comparable to the 2014 CODATA value. In 2017 a group in Garching has published the value r p E = 0.8335(95) fm from 2S − 4P transition [7]. In 2018 a group in Paris has published the value r p E = 0.877 (13) fm from 1S − 3S transition [8]. More results are expected in the very near future [9].
Muonic hydrogen spectroscopy: No issues were raised concerning the muonic hydrogen experimental measurement. It should be noted that the extraction of r p E from muonic hydrogen spectroscopy is based on the work of one group. There are no plans by other groups to repeat the measurement. On the theoretical side there was much discussion in the literature. It is reflected in the different theoretical formula used in [1] and the 2013 update [3]. Due to its precision, the muonic hydrogen result involves a more complicated hadronic input, beyond a one-photon probe of the proton structure. In particular, as emphasized in [10], two-photon effects are a potential source of uncertainty. The imaginary part of the two-photon exchange amplitude is related to experimental data: form factors and structure functions, see section 4.2 below. The amplitude cannot be reconstructed from its imaginary part and the knowledge of a subtraction function W 1 (0, Q 2 ), where Q 2 = −q 2 , is required, see (30) below. The subtraction function W 1 (0, Q 2 ) is not known exactly. Its small Q 2 expression is calculable using Non-Relativistic QED (NRQED) [10]. Its large Q 2 expression can be calculated using the operator product expansion. The spin-0 contribution was calculated in [11] and corrected in [12]. The spin-2 contribution was calculated in [12]. In [12] the two limits were interpolated to give an estimate for the contribution of W 1 (0, Q 2 ) to two-photon exchange effects. The uncertainty on the interpolation is larger than in [3], but it is too small to explain the discrepancy. The estimate in [12] is consistent with the literature [13,14,15,16,17,18,19,20]. On the other hand, [21] finds a much larger uncertainty. Ultimately one would like to probe the muon-proton two-photon exchange effects by using a different method such as muon-proton scattering.
Electron-proton scattering: Extractions of r p E from the electron-proton cross section data or even from the form factor itself require an extrapolation to q 2 = 0. Since G E 's functional form is not known, such an extrapolation is not simple. In [22] it was first suggested to utilize the z expansion, used before for meson form factors, to perform the extrapolation. The z expansion is based on the known analytic properties of the form factor. The values obtained using z-expansion analyses [22,23] generally disfavor the muonic hydrogen result 2 . More recently, lattice extractions of the charge radius have used the z expansion, see e.g. [27,28,29], although currently the errors are typically too large to distinguish between the two values of r p E . On the experimental side a new low-Q 2 electron-proton scattering experiment called "PRad" [30] was recently performed at Jefferson Lab and its results are expected to be published in the near future.
Muon-proton scattering: This is the least studied method to extract r p E . To address that, a new muon-proton scattering experiment called MUSE is being built at the Paul Scherrer Institute [31]. It will start taking data in 2019. It is the first muon scattering measurement with the required precision to address the proton radius puzzle [9].
In muonic hydrogen the muon's typical momentum is mα ∼ 1 MeV, and both the muon and the proton can be treated non-relativistically. For the MUSE experiment the muon's typical momentum is of the order of the muon mass m ∼ 100 MeV, and the muon must be treated relativistically, while the proton can be treated non-relativistically. The appropriate effective field theory (EFT) for such kinematics was suggested in [32] and is called QED-NRQED 3 .
In a previous paper we have studied some aspects of this effective field theory [33]. Denoting by m(M ) the muon (proton) mass and using Z = 1 for a proton, we showed that one-photon exchange O(Zα) QED-NRQED scattering at power 1/M 2 reproduces Rosenbluth scattering [34], and the two-photon exchange O(Z 2 α 2 ) QED-NRQED scattering at leading power reproduces the scattering of a relativistic fermion off a static potential [35,36].
Two photon exchange effects sensitive to the proton structure start at O(Z 2 α 2 ) and power 1/M 2 . For QED-NRQED they appear as two four-fermion operators, one spin-independent and one spin-dependent, see (4) below. This paper's goal is to determine their Wilson coefficients in terms of the proton's hadronic tensor, i.e. performing a matching onto QED-NRQED.
The paper is structured as follows. In section 2 we present an overview of the general features of the matching calculation. In section 3 we preform the effective field theory calculation. In section 4 we preform the full theory calculation for the toy example of a non-relativistic point particle and for the proton. In section 5 we extract the matching coefficients for both full theories. We present our conclusions in section 6. The appendices describe the matching in Coulomb gauge, properties of the hadronic tensor, and list the NRQED Feynman rules.
2 General features of the matching calculation

Lagrangian
We consider the interaction between a relativistic spin-half field ("lepton") of mass m, denoted by , and a non-relativistic spin-half field ("non-relativistic proton") of mass M , denoted by ψ. For the calculation of the Wilson coefficients of the leading power four-fermion operators from the proton's hadronic tensor, we will need the interactions at O(Z 2 α 2 ) and power 1/M 2 , where Z = 1 for a proton. The QED-NRQED Lagrangian contains three parts that are relevant to the calculation.
The first part is the NRQED Lagrangian 4 , up to 1/M 2 : (1) where D t = ∂/∂t + ieZA 0 , D = ∇ − ieZA, σ are the Pauli matrices, and e is the electromagnetic coupling constant 5 . These are the components of D µ = ∂ µ + ieZA µ . The notation [∇·E] denotes that the derivative is acting only on E and not on ψ. For a review see [38]. The (hidden) Lorentz invariance of the Lagrangian, also known as reparameterization invariance, implies that c 2 = 1 [39,40,32]. The other Wilson coefficients can be related to the proton electromagnetic form factors defined by . The latter can also be determined by the hidden Lorentz invariance of the Lagrangian [39,40,32]. The NRQED Feynman rules are listed in appendix C. The second part is the usual QED Lagrangian that describes the lepton's interaction with the electromagnetic field: where Q = −1 for a muon or an electron 6 . The third part is the QED-NRQED contact interactions. At 1/M 2 we have two possible contact interactions, 4 At this power there are also interactions of four non-relativistic spin-half fields, see [33]. We will not need them in this paper. 5 We follow the conventions of [37], although in that paper the NRQED Lagrangian describes an electron. In other words, as in [33], we take e to be positive. 6 We do not include 1/M 2 operators of [32], since they have Wilson coefficients that start at O(α). At the lowest order in α these operators lead to terms of order O(Zα 2 ) but not O(Z 2 α 2 ) which are relevant to the QED-NRQED contact interactions.
where our notation follows that of [32]. The main goal of this paper is to express the Wilson coefficients b 1 and b 2 in terms of the components of the hadronic tensor. In the matching we will encounter two other operators that are explicitly suppressed by m/M 3 . These are b 3 mψ † ψ¯ /M 3 and 2b 6 mψ † σ i ψ¯ i 2 ijk γ j γ k /M 3 , in the notation of [32]. The other 1/M 3 contact interactions contain space-like covariant derivatives. Since we will set the external three-momenta to zero in the matching, we will not encounter them.

Infrared and ultraviolet singularities
We calculate the amplitude for forward off-shell + p → + p scattering at O(Z 2 α 2 ) and power 1/M 2 both in the full and the effective theory. The former is expressed in terms of the hadronic tensor defined below, and the latter in terms of the Wilson coefficients of the effective theory. Both amplitudes are IR singular. We regulate these by using a fictitious photon "mass" denoted by λ. In terms of powers of λ we will find singularities 7 that scale like 1/λ 3 , 1/λ 2 , 1/λ, and log λ.
The lepton mass m is also an IR quantity. Thus we find both in the full and effective theory terms that diverge in the limit m → 0. In terms of powers of m we will find singularities that scale like 1/m 2 , 1/m, and log m. The IR terms regulated by λ, m, or both, must cancel in the matching. This serves as a non-trivial check of the calculation.
We assume m, λ M , so we will expand the singular terms in powers of m/M and λ/M . This allows to simplify the matching calculation. This obviously holds for the muon-proton case where m/M ∼ 0.1.
The effective field theory is valid up to a cutoff scale Λ ∼ M . We use a hard momentum cutoff to regularize the UV singularities of the QED-NRQED integrals. In the matching we also need to expand the form factors in power of q 2 . This can introduce UV divergent terms in the full theory calculation. In practice this is not a problem for the matching calculation, since such term are regularized when using the full functional form of the form factors.

Gauge Invariance
Since the photon propagator in Coulomb gauge is different for space and time components, it is often used for NRQED calculations, see e.g. [37]. In the following we perform the matching both in Feynman and Coulomb gauges. Both the full and effective field theory amplitudes are off-shell and as a result each amplitude need not be gauge invariant by itself. As we will see, both the full and effective theory calculations are different for each gauge, but the Wilson coefficients b 1 and b 2 are the same, as they should be.
Interestingly, when we take the non-relativistic limit for the lepton, the effective field theory amplitude is the same in both gauges. In this limit one adds ψ † ψ¯ γ 0 and ψ † ψ¯ to form the spin-independent amplitude and ψ † σ i ψ¯ γ i γ 5 and ψ † σ i ψ¯ i 2 ijk γ j γ k to form the spin-dependent amplitude.
The difference between the gauges arises even for the simpler case of non-relativistic point particle interacting with a relativistic point particle. As we show in appendix A.2, the amplitude is the same for Feynman and general covariant gauge, but different in Coulomb gauge. Taking the non-relativistic limit for both particles results in a gauge invariant amplitude. In a nutshell, the reason is that for p = (m, 0) and χ a non-relativistic spinor, (/ p − m)(χ 0) T = (mγ 0 − m)(χ 0) T = 0 making the non-relativistic spinor effectively on-shell in this limit. The details are given in appendix A.2.
In the following we present results in Feynman gauge. The results in Coulomb gauge are presented in appendix A.

Effective Field Theory Calculation
We calculate the amplitude for forward off-shell +p → +p scattering at O(Z 2 α 2 ) and power 1/M 2 in QED-NRQED effective theory using Feynman gauge. The QED-NRQED calculation in Coulomb gauge is presented in appendix A. The interactions of the relativistic fermion ( ) with the photon arise from (3). The interactions of the non-relativistic fermion (p) with the photon arise from (1) and the Feynman rules are given in appendix C. Combining the two we find the diagrams shown in figure 1 that contribute at power 1/M 2 .

Example calculation
We illustrate the calculation of such QED-NRQED loop diagrams, by considering the leading power diagrams. These are the first two diagrams in figure 1. We will refer to the first diagram as "direct" and the second as "crossed". These diagrams arise from the NRQED field coupling in first term of (1). Applying the Feynman rules in appendix C we obtain where we have defined the integrals Our notation is such that the numerator m (m − l 0 ) corresponds to the superscript m (0) and the plus (minus) sign in the last term in the denominator correspond to the subscript D (C). The presence of both l 0 and l 2 in the denominator implies that the standard covariant Feynman integral technics cannot be used. Also, unlike pure NRQED loop diagrams, see e.g. [37], we also have relativistic fermion propagators. Furthermore, we need only terms that are singular in λ and O(λ 0 ) terms up to power 1/M 2 . We will present two methods to calculate these integrals.
In the first method we use contour integration for the l 0 integrals. Since we have no external 3-momenta in the problem, the angular integration is also immediate. We are left with the integral over | l|. It is convenient at this stage to rescale λ, m, and | l| by M so the integral depends only on λ/M and m/M . We use a simple version of the method of regions [43]. We divide the integration domain to a region where | l| ∼ λ and a region where | l| λ. The two  regions are separated by an arbitrary scale that cancels when we add the two together. For the lower region we expand in | l| ∼ λ. For the upper region we can set λ to zero. We add the two regions, confirm that the arbitrary cutoff cancels, and expand the O(λ 0 ) terms up to 1/M 2 . The results of the integrals in (6) are In the second method we expand the NRQED propagator before we calculate the integrals: We then use contour integration for l 0 integrals, perform the angular integration and integrate over | l|. Since l 2 appears only in the numerator after the expansion, the last integration is easier to do than in the previous method . Furthermore, we automatically generate the series in 1/M . The price we pay is that the IR singular terms in λ are expanded in 1/M too. This is not a problem in practice since such singular terms cancel in the matching. The presence of the l 2 in the numerator might make the integral UV divergent. We regularize those by a cutoff Λ ∼ M . Using the second method the integrals in (6) are It is easy to check that these expressions correspond the 1/M expansion of (7) by choosing the value of Λ appropriately. Due to its simplicity, we will use the second method in the results presented below. We note that most of the direct and crossed terms that contain an even inverse powers of M differ in sign, while those with odd power have the same sign. This is easily understood from (8) where such a pattern appears in the expanded propagator. For singular terms such relations may not hold. In particular, the direct diagram alone has a 1/λ 3 singularity.
Because of such relations, many terms cancel if we add the direct and crossed integrals. This simplifies the structure of the spin-independent parts of the amplitude that typically contain a sum of the direct and crossed integrals. Furthermore, since the amplitude scales as inverse mass squared, the terms proportional to 1/M will be multiplied by 1/m or 1/λ. Thus they are IR divergent and cancel in the matching. Terms proportional to 1/M 2 typically cancel in the sum and do not contribute to the matching coefficient of the spin-independent terms. A non-zero contribution will have to come from the full theory amplitude. This has important implications for what follows. The spin-dependent parts of the amplitude typically contain the difference of the direct and crossed integrals. Thus terms proportional to 1/M 2 would not cancel and may contribute to the matching coefficient of the spin-dependent terms.
The leading power amplitude, corresponding to the first two diagrams in figure 1, is .
We note that the entire leading power effective field theory contributions is given by terms that diverge when λ and/or m go to zero. Such IR terms will cancel in the matching.

Final result
Using similar methods we calculate the other diagrams in Figure 1. The sum of all the diagrams is Notice that the sum of all spin-independent terms diverges when λ and/or m go to zero. Such IR terms will cancel in the matching. The spin-dependent terms contain non-IR terms that can contribute to the matching coefficients, unless they match the full theory terms. The results of the QED-NRQED calculation in the Coulomb gauge are given in appendix A.

Toy example: non-relativistic point particle
For the full theory calculation we find it instructive to consider first the toy example of the point-particle "proton". While obviously this is not a realistic model, it exhibits some of the features we will encounter later in the rigorous description of the proton in terms of the hadronic tensor. We consider the interaction of a relativistic point particle, e.g. a lepton, with a nonrelativistic point-particle "proton". There are two diagrams that contribute to the amplitude, a "direct" diagram and a "crossed" diagram. The amplitude is where and v = (1, 0), k = mv, p = M v. From now on we will suppress the argument of u(k).
Since we will match onto QED-NRQED, we use non-relativistic normalization for the proton spinors u(p) † u(p) = 1. We also use the Dirac representation of the γ matrices where we have used the identity Neglecting terms linear in l that integrate to zero we find where for integrals over l i l j we replace l i l j → δ ij l 2 /3. We need the following integrals.
. (18) To calculate them we use the partial fractioning identity and define We can now express I, I 0 , I 00 ,Ĩ in terms of i, i 0 , i 00 ,ĩ: In calculating the i integrals it is convenient to combine denominators first via [44] 1 Calculating i, i 0 , i 00 ,ĩ we find Using the expressions above, the amplitude is We will see in appendix A, that in Coulomb gauge this amplitude (24) is different. In the next section we will use this result to extract b 1 and b 2 for the non-relativistic point particle case. As a check, we can take the non-relativistic limit for the lepton, i.e.ūu,ūγ 0 u → χ † χ , and which is a known result [41]. We will see in appendix A, that unlike (24), we get the same result for (25) in Coulomb gauge.

The proton
For the proton the forward lepton-proton O(Z 2 α 2 ) amplitude can be expressed in terms of the hadronic tensor that encodes the two-photon interaction of the proton. Here p is the proton momentum and s is its spin. As we review in appendix B, W µν (p, q) can be expressed as The four scalar functions W 1 , . As a result, W 1 , W 2 , H 1 are even functions of ν, and H 2 is odd, see appendix B. For a point particle, In terms of the hadronic tensor, the amplitude in Feynman gauge is We insert (27) into (28) and set the lepton 3-momentum to zero, i.e. k = (m, 0) and work in the rest frame of proton, i.e. p = (M, 0). To match onto NRQED, the proton Dirac spinors in (27) should be used with non-relativistic normalization u(p) † u(p) = 1 and we replace u p → (χ 0) T , where χ is a two-component spinor. After these simplifications we find In deriving this equation, we used that W 1 and W 2 are even functions of ν = 2M l 0 . As a check, we insert the point particle values of W 1 , W 2 , H 1 , H 2 into (29) and reproduce the result of (24).
To proceed, we use dispersion relations for the scalar functions W 1 , W 2 , H 1 , H 2 . For a given value of Q 2 the singularities of these function lie on the real ν axis and are symmetric under ν → −ν see, e.g., [44]. These consist of a pole at ν = ±Q 2 , corresponding to proton state, and a cut that starts at ν cut (Q 2 ) = Q 2 + 2m π M + m 2 π , the threshold for a production of a proton and pion pair. Using the symmetry of the scalar functions under ν → −ν we can express them as integrals over positive value of ν starting at ν 0 ≡ Q 2 .
W 1 satisfies a subtracted dispersions relation in ν, see [12] for a proof that utilizes the optical and low energy theorems. We assume that W 2 , H 1 , H 2 satisfy unsubtracted dispersions relations From the definition of W µν (26) one can obtain the identity [44] 2 Im W µν = X p, s|J µ |X X|J ν |p, s (2π) 4 δ (4) (q − p X + p) where p X is the four-momentum of the state X. The summation includes phase space integrals and sum over spins of the state X. This equation allows us to separate the proton and the continuum contribution in the imaginary part of W µν . The proton contribution is given in terms of form factors, and the continuum contribution in terms of inelastic structure functions.
Using the dispersion relations we can express the scalar functions as [10] where the superscript numbers denote the number of subtractions, the superscript p is the proton contribution, and the superscript c is the continuum contribution.
Using (2) we can express the proton contribution in (32) in terms of the form factors where . After simplifying the Dirac structure, as explained in appendix B, we find Inserting these expressions into (30) and (31) gives W p,1 1 (ν, Q 2 ) and W p,0 2 (ν, Q 2 ) agree with [10]. H p,0 1 (ν, Q 2 ) and H p,0 2 (ν, Q 2 ) agree, up to an overall normalization with the ones given by Drell and Sullivan 8 in [45].
The subtraction function W 1 (0, Q 2 ) is not known exactly. In the large Q 2 limit it can be calculated using the operator product expansion. The spin-0 contribution was calculated in [11] and corrected in [12]. The spin-2 contribution was calculated in [12]. For the matching we need the small Q 2 expression that can be calculated using NRQED [10] where m p is the proton mass, F 2 (0) = a p , r p E and r p M are the proton electric and magnetic charge radii, andβ is the magnetic polarizability. This equation assumes Z = 1. For the matching we will only need the first term W 1 (0, To match the QED-NRQED singularities we will need a contribution to H 1 (ν, Q 2 ) in the Q 2 → 0 limit. From the low energy theorems of Low [46] and Gell-Mann and Goldberger [47] we know that H 1 (ν, Q 2 ) includes also a term F 2 (0) 2 /(2M 2 ). We will follow Drell and Sullivan [45] and include it in the continuum contribution to H 1 (ν, Q 2 ), namely H c,0 1 (ν, Q 2 ) of (33). In the literature one often finds a different approach to calculate the proton contribution to W µν . In this approach the on-shell vertex of (2) is inserted into the off-shell amplitude of (26). The hadronic tensor is then calculated with Feynman rules that depend on F 1 and F 2 . This approach is often called "Born terms" or "Born approximation", see e.g. [45], although this is not a first term in a perturbative series which is the usual meaning of the Born approximation in quantum mechanics. In this approach one finds Compared to the expressions above, W 2 and H 2 are unchanged. The extra term for W 1 corresponds to assuming W 1 (0, Q 2 ) = 2F 2 (2F 1 + F 2 ). The extra term for H 1 is F 2 2 /2M 2 . Interestingly its value at Q 2 = 0 corresponds to the F 2 (0) 2 /(2M 2 ) term that we include in the continuum contribution. This was also noted in [45]. For matching the QED-NRQED IR singularities, the difference between the extra terms in (38) compared to the terms we have included above is irrelevant.
In summary, to match QED-NRQED IR singularities we need: (36), W 1 (0, Q 2 ) to zeroth order in Q 2 , namely 2F 2 (0) [2F 1 (0) + F 2 (0)], and a part of H c,0 1 (ν, Q 2 ): F 2 (0) 2 /(2M 2 ). In total The second line for W 1 and H 1 in (39) does not contribute to the IR singular terms. Using this information we return to the calculation of (29). First, to simplify the notation we rename l → q. The poles in (29) arise from the lepton propagator, the photon propagators (regularized by λ), and the poles from (36). The location of the poles is such that we can perform a Wick rotation l 0 In terms of the new variables we have Since ν = 2M iQx, W 1 , W 2 , H 1 (H 2 ) are even (odd) function(s) of x. This implies, Notice the factor of m in front of all the spin-independent terms. We now insert the expressions from (39). The IR singular terms arise from the Taylor expansion of the form factors in Q 2 . Expanding, integrating over x and Q, and expanding the result in inverse powers of M gives The log(Q/M ) in the last line of (42) is a UV divergence that arises from the Taylor expansion of the form factor. Such a term would be regulated when using the full functional form of F 2 . As a check, we can add the χ † χūu and χ † χūγ 0 u terms and compare the sum to the NRQED result on the left hand side of equation (7) of [10]. Setting F 1 (0) = 1 and expanding that result to order 1/M 2 , we find a complete agreement. As another check, the terms in (42) proportional to F 1 (0) 2 match the IR singular terms of (24) expanded to order 1/M 2 .
The IR divergences in (42) match exactly the IR divergence of (12), as they should. Furthermore, the χ † χūu and χ † χūγ 0 u terms include only IR divergent terms. The χ † σ i χū i 2 ijk γ j γ k u term is equal to corresponding term in (12) including the non-IR term. This would lead to a zero matching coefficient at order 1/M 2 as expected from chiral symmetry. In the next section we will use (42) to extract the matching coefficients.

Extraction of the matching coefficients at power 1/M 2
With the effective field theory and full theory calculations at hand, we can now extract the matching coefficients b 1 and b 2 of (4) which is the goal of this paper. We will calculate them for two cases. First, for the toy example of a non-relativistic point particle. In this case since the full theory calculation is known explicitly, we can find explicit expressions for b 1 and b 2 . Second, for the case of the proton. In this case the full theory is expressed in term of the components of the hadronic tensor and we will give an implicit expression for the matching coefficients in terms of an integral over these components.

Toy example: Extraction of matching coefficients for a nonrelativistic point particle
For this case the full theory amplitude is given by (24). To perform the matching we need to expand this equation to power 1/M 2 . We find To match onto the effective theory, we need to take (12) and set the values of the Wilson coefficients to that of a point particle, namely, c p.p.
The difference between the full and effective theory is As expected from chiral symmetry, the coefficient of the structure χ † χūu and that of χ † σ i χū i 2 ijk γ j γ k u both cancel in the matching. Surprisingly, also the coefficient of the structure χ † χūγ 0 u cancels in the matching. We conclude that for the toy example of a point particle, b p.p.
One may ask at which order would the structure χ † χūγ 0 u get a non-zero contribution. Looking at equation (24), the only term that is not divergent in the λ → 0 or m → 0 limit is 4m log M/[M (m 2 − M 2 )]. This term vanishes at power 1/M 2 , but at power 1/M 3 it gives −4m log M/M 3 . Since when expanding the propagator log M cannot be generated by the effective field theory calculation, we expect that this would give rise to a non-zero contribution to a matching coefficient. Since we have not calculated the QED-NRQED amplitude to power 1/M 3 , we cannot determine it at this stage. We emphasize again that the results in (46) are for a toy example that in no way represents the behavior of the real proton. We consider the case of the proton next.

Extraction of matching coefficients for the proton
For the proton the full theory amplitude is given by (41). To find b 1 and b 2 we only need the structures χ † χūγ 0 u and χ † σ i χūγ i γ 5 u. The QED-NRQED result is given in (12). The Wilson coefficient b 1 is determined by the relation The Wilson coefficient b 2 is determined by the relation Equations (47) and (48) are the main results of this paper. Given the expressions for the W 1 , W 2 , H 1 , H 2 one can find an explicit expression for b 1 and b 2 . In lieu of that, we can use (42) to find the contribution to b 1 and b 2 from the full theory that are proportional to F 1 (0), F 2 (0) and M 2 F 1 (0). Comparing (42) to (12) we find first that there is no contribution to matching coefficients for the structure χ † χūu and χ † σ i χū i 2 ijk γ j γ k u at power 1/M 2 as expected from chiral symmetry and similar to the point particle toy example. The contributions to b 1 and b 2 are As discussed in the previous section, the log(Q/M ) is a UV divergence that arises from the Taylor expansion of the form factor. Such a term would be regulated when using the full functional form of F 2 . The F 1 (0) 2 term agrees with the toy example of a point particle.
Surprisingly, there is no contribution to b 1 . Would that still hold if we include the full expression for W 1 and W 2 ? Comparing (47) and (48) we see that the former has an explicit factor of m in the left hand side. This extra factor of m can be understood by looking at (40). All the spin-independent terms in (40) contain an explicit factor of m with the exception of the term proportional to 2iQx. Since this term multiplies W 1 which is an even function of x, we must combine it with 2imQx from the lepton propagator. This introduced an extra factor of m also for this term, leading to the overall factor of m in (47). As a result, the entire spin-independent contribution has an extra factor of m. The only difference between the two structures is the coefficients of W 1 in (41), −3m for χ † χūu and m(1 − 4x 2 ) for χ † χūγ 0 u. We see no reason why the coefficient of the former would receive no contribution at power 1/M 2 , while the coefficient of the latter would, when including the full expression for W 1 .

Conclusions and outlook
Motivated by the proton radius puzzle, a new elastic muon-proton scattering experiment called MUSE will soon start taking data. In this experiment muons with energy close to the muon mass scatter of protons. In such an energy regime one can describe the proton using NRQED, but one has to use QED to describe the muon. An appropriate effective field theory for such kinematics is QED-NRQED.
QED-NRQED is an effective field theory that describes the electromagnetic interaction of a relativistic point particle, e.g. a lepton, with a non-relativistic spin-half, possibly composite, particle, e.g. a proton. The QED-NRQED Lagrangian can be organized in inverse powers of M , the non-relativistic particle mass. For addressing the proton radius puzzle the most important operators are the dimension six spin-independent ones, namely the Darwin term ψ † c D e[∇ · E]ψ/8M 2 , and the contact interaction b 1 ψ † ψ¯ γ 0 /M 2 . The Wilson coefficient c D is related to the proton charge radius, and b 1 encodes two-photon exchanges contribution from physics above the scale M .
In a pervious paper we have shown how this effective field theory reproduces the Rosenbluth scattering amplitude at O(Zα) and power 1/M 2 and a relativistic scattering off a static potential at O(Z 2 α 2 ) scattering at leading power in 1/M . In this paper we determined the Wilson coefficients of the contact interactions at O(Z 2 α 2 ) and power 1/M 2 .
To determine the coefficients we did a matching calculation where we match a full theory onto QED-NRQED. We have considered two cases of full theories: A toy example of a nonrelativistic point particle, and the real proton which is described by a hadronic tensor. For both cases we have done the matching calculation in Feynman as well as in Coulomb gauge.
The effective field theory calculation in Feynman gauge is described in section 3 and in Coulomb gauge in appendix A. We calculated the amplitude for forward off-shell + p → + p scattering at O(Z 2 α 2 ) and power 1/M 2 in QED-NRQED. The relevant Feynman diagrams are shown in figure 1. The typical loop integrals in QED-NRQED involve both relativistic and non-relativistic propagators. We presented two methods for such a calculation and explicitly demonstrated their equivalence for the leading power diagrams. The easier method to use involves an expansion of the NRQED propagator in power of 1/M before performing the loop integrals. The complete QED-NRQED amplitude is given in (12) for Feynman gauge and in (50) for Coulomb gauge. The two amplitudes are not the same. There are four structures in the amplitude χ † χūu, χ † σ i χū i 2 ijk γ j γ k u, χ † χūγ 0 u, and χ † σ i χūγ i γ 5 u. The first two involve an even number of Dirac gamma matrices and are expected to be chiraly suppressed. Only the last two are expected to lead to a non-zero Wilson coefficient at power 1/M 2 .
The full theory calculation in Feynman gauge is described in section 4 and in Coulomb gauge in appendix A. For the toy example of the non-relativistic point particle we presented an explicit expression for the amplitude in (24) for Feynman gauge and in (51) for Coulomb gauge. As in the effective field theory case, the two amplitudes are different in different gauges. Interestingly when we add in pairs the spin-independent terms and the spin-dependent terms we find that the two gauges agree. We explain the reason in appendix A.2. For the case of a real proton we present an implicit expression for the amplitude given in terms of integrals over the scalar functions that multiply the components of the hadronic tensor. The expression for the amplitude in Feynman gauge is given in (41) and in Coulomb gauge in (52). The two amplitudes are different in the different gauges.
We preformed the matching and the extraction of the Wilson coefficients b 1 and b 2 in section 5. Both the full and effective field theory amplitudes are IR divergent. We regulate these IR singularities by using a fictitious photon "mass" denoted by λ. The lepton mass m is also an IR quantity in the full and the effective theory. Each type of IR singularities, as well as mixed IR singularities, e.g. log(m/λ)/(mM ), cancel in the matching. For the toy example of a non-relativistic point particle the explicit expressions for b 1 and b 2 are given in (46) for Feynman gauge and in (54) for Coulomb gauge. Although both the full and effective theory amplitudes are different in each gauge, the Wilson coefficients are the same for both gauges. Surprisingly we find that b p.p.
. For the case of the proton we give implicit expressions for b 1 and b 2 in (47) and (48) for Feynman gauge. These are the main results of the paper. We show explicitly how the IR singularities cancel between the full and effective theory in both Feynman and Coulomb gauge based on the known low-energy behavior of the hadronic tensor. We find that b 1 = 0 at O(Z 2 α 2 ). Based on the structure of the full theory integral in Coulomb gauge (41) and in Feynman gauge (52) we argued that b 1 = 0 at O(Z 2 α 2 ) even if we include the full expression of the hadronic tensor. The reason is that the full theory expression includes an explicit factor of m in those integrals for all the spin-independent terms. The spin-independent part of the effective field theory expression contains only IR singular terms that cancel in the matching and do not contribute to b 1 .
The fact that b 1 = 0 at O(Z 2 α 2 ) is surprising. One might ask at which power two-photon exchange effects from scales above M enter. For the toy example of a point particle there is a term in the full theory amplitude that would give rise to such a spin-independent matching coefficient at power 1/M 3 . At this power we expect to find also other operators, some of which depend on the lepton momentum.
This implies that scattering amplitudes calculated within QED-NRQED are not sensitive to two-photon exchange effects from scales above M at O(Z 2 α 2 ) and power 1/M 2 . At power 1/M 2 the spin-independent part of the amplitude depends only on c D , or in other words, on the proton charge radius. For the case of m being the muon mass and M the proton mass, we conclude that two-photon exchange effects from scales above M are an order of magnitude smaller compared to proton charge radius effects. One of the suggested explanations for the proton radius puzzle involves an unusual behavior of W 1 (0, Q 2 ), since only its asymptotic low and high Q 2 terms are known. The vanishing of b 1 at O(Z 2 α 2 ) implies that the MUSE experiment will be less sensitive to such effects, but its extraction of the proton charge radius will be more robust.
The vanishing of b 1 at O(Z 2 α 2 ) arises from a combination of two phenomena. The first arises since the effective field theory amplitude contains only IR divergent terms for the spinindependent part of the amplitude. This happens since it is a sum of direct and crossed integrals that tend to cancel each other. Such IR terms must vanish in the matching and give no contribution to b 1 . The second effect arises from the full theory amplitude. The spinindependent parts of the hadronic tensor, namely W 1 and W 2 , are multiplied by an explicit factor of m or a term linear in the photon energy (measured in the proton rest frame). Since W 1 and W 2 are even functions of the photon energy such terms vanish unless they are multiplied by a product of the lepton mass and the photon energy from the lepton propagator, see (29) for Feynman gauge and (52) for Coulomb gauge. The combination of the two phenomena leads to the vanishing of b 1 at O(Z 2 α 2 ).
Beyond the proton radius puzzle, this result can be of interest in physics beyond the standard model, where generating hierarchies, even "little" ones, between the weak scale and the scale of new physics is an active topic of research. It would be interesting to see if the vanishing of b 1 at O(Z 2 α 2 ) can be used to generate such hierarchies.

Coulomb gauge is
Comparing to the QED-NRQED result in Feynman gauge (12) we notice several differences in each structure. Interestingly, if we add the terms multiplying χ † χūu and χ † χūγ 0 u we find the same answer as the corresponding sum in (12). The same is true if we add the terms multiplying χ † σ i χū i 2 ijk γ j γ k u and χ † σ i χūγ i γ 5 u. We turn now to the full theory calculation. For the toy example of a non-relativistic point particle we find in Coulomb gauge (39), performing the integrals, and expanding in inverse powers of M gives As for Feynman gauge, the log(Q/M ) in the last line is a UV divergence that arises from the Taylor expansion of the form factor. Such a term would be regulated when using the full functional form of F 2 .
As a check, we can add the χ † χūu and χ † χūγ 0 u terms and compare the sum to the NRQED result on the left hand side of equation (7) of [10]. Setting F 1 (0) = 1 and expanding that result to order 1/M 2 , we find a complete agreement. As another check, the terms in (53) proportional to F 1 (0) 2 match the IR singular terms of (51) expanded to order 1/M 2 .
We can use these results to extract the matching coefficients following the same procedure as in section 5. First, for the toy example of a non-relativistic point particle. Setting the Wilson coefficients to their point particle values in (50) and subtracting the result from the order 1/M 2 expansion of (51) we find that in Coulomb gauge b p.p.
This is the same result as in Feynman gauge, (46). Although the full and effective field theory differ between Feynman and Coulomb gauge, the Wilson coefficients are the same in both gauges, as one would expect. For the case of the proton the Coulomb gauge expressions analogous to (47) and (48) are longer and we will not give them here. To compare to Feynman gauge, we find the contribution to b 1 and b 2 from the full theory that are proportional to F 1 (0), F 2 (0) and M 2 F 1 (0). These are This is the same result as in Feynman gauge, (46). In particular there is no contribution to b 1 . Since in (52) the spin independent terms either include an explicit factor of m or a term linear in l 0 , we expect no contribution to b 1 even if we include the full expression for W 1 and W 2 .

A.2 Gauge invariance
We have found that the matching coefficients are equal for Feynman and Coulomb gauge as one would expect. The effective and full theory amplitudes, on the other hand, are different in each gauge. This can be expected. Since we want to distinguish, e.g.,ūγ 0 u fromūu, we cannot put the lepton on-shell. Since the amplitudes are off-shell, they are not guaranteed to be gauge independent. This is most prominent for the point particle case where we can compute the full theory amplitude explicitly. Surprisingly, when taking the non-relativistic limit for the lepton, the two gauges agree. We now explain this for the non-relativistic point particle case. For a non-relativistic point particle case we take u(p) → (χ 0) T . At its rest frame, p 0 = M , (/ p − M )u(p) → (p 0 γ 0 − M )(χ 0) T = 0. In other words, the non-relativistic particle is "effectively" on-shell. Consider now a general covariant gauge where the photon propagator is The (ξ − 1) terms add an / l on the non-relativistic point particle line. For example, one of the terms givesū(p)γ µ [/ l + M (1 + γ 0 )] / lu(p) in a the general covariant gauge analog of (14).
We now write / l = ±(/ p ± / l − M ) ∓ (/ p − M ). The first term cancels the non-relativistic point particle propagator and will vanish when we add the direct and crossed contribution to the amplitude. The second term vanishes as a result of the "effective" on-shell condition for the non-relativistic point particle. We conclude that a general covariant gauge amplitude is the same as that of a Feynman gauge. We can try and apply a similar procedure for Coulomb gauge. It turns out that the gauge dependent parts in Coulomb gauge arise from terms that include γ i l i . We can use the identity γ i l i = γ 0 l 0 − / l to generate terms analogous to covariant gauge. But unlike covariant gauge we have terms, e.g.,ū(k) generates terms with an / l only on the lepton line. Such terms do not vanish for a relativistic lepton thus generating Coulomb gauge dependent pieces. Once we take the lepton to be also non-relativistic, we can apply the same procedure as before and use that / k − m annihilates the lepton spinor in the non-relativistic limit. This explains why in the point particle case the relativistic-non-relativistic amplitude is different in Coulomb gauge while the non-relativisticnon-relativistic amplitude is the same in Coulomb and Feynman gauges.
We expect a similar reasoning to apply for the QED-NRQED amplitude and the proton amplitude, but we have not proven it in these cases. The former since it must reproduce the IR singularities of the point particle amplitude, and the latter since it is a generalization of the point particle case.

B Appendix: Properties of the Hadronic tensor
For completeness we derive several properties of the hadronic tensor W µν defined in (26). The structure of W µν is constrained by current conservation, by being a forward matrix element, and by the parity and time reversal symmetries of the electromagnetic and strong interactions. First, current conservation implies that q µ W µν = 0 and q ν W µν = 0.
Second, W µν is a forward matrix element so we can write it as where Γ µν is a general Dirac structure. The Dirac structure can be simplified by using that for on-shell spinors / vu(p) = u(p) where / v ≡ / p/m. In analogy to Heavy Quark Effective Theory (HQET) we can define the projector P + ≡ (1 + / v)/2, where P + u(p) = u(p). As was shown in [48], between two P + 's the Dirac basis reduces to four matrices: P + and s µ = P + γ µ γ 5 P + . The matrices s µ are a generalization of the Pauli spin matrices that satisfy v · s = 0.
Sinceū p (p, s)u p (p, s) is independent of the spin, the tensor multiplying the unit matrix in W µν must be symmetric under the interchange µ ↔ ν. The tensor multiplying s αβ in W µν is antisymmetric under the interchange µ ↔ ν. This is most easily seen in the rest frame of the proton, where s αβ are just the Pauli matrices, i.e. s ij = ijk σ k . For the Pauli matrices χ † ↑ σ i χ ↑ = −χ † ↓ σ i χ ↓ Consider now the identity matrix. Since it is symmetric under µ ↔ ν, it is multiplied by a linear combination of the tensors g µν , p µ p ν , q µ q ν and (p µ q ν + q µ p ν ). Using q µ W µν = 0 and q ν W µν = 0, we can reduce these to two linear combinations −g µν + q µ q ν /q 2 and (p µ − p · q q µ /q 2 )(p ν − p · q q ν /q 2 ).
All together we find equation (27). Since p 2 = M 2 , the scalar coefficients of these structures, W 1 , W 2 , H 1 , H 2 depend on q 2 and p · q, or alternatively on Q 2 and ν. Translation Combining W µν (p, q) = W νµ (p, −q) with (27) implies that W 1 , W 2 , H 1 are even functions of ν, while H 2 is an odd function of ν.

C NRQED Feynman rules
For QED-NRQED amplitudes we need the Feynman rules of both QED and NRQED. The QED Feynman rules are well-known and will not be listed here. Most of the NRQED Feynman rules used in this paper, but not all, are given in figure 3 of [37] by multiplying the vertices by −i and the propagators by i. That table does not include one of the Feynman rules arising from the term multiplying c S . Also, Coulomb gauge was used in [37]. In a non-Coulomb gauge, apart from the different photon propagator, there is also an additional interaction from the term multiplying c D . We list below all the Feynman rules up to order 1/M 2 for both Feynman and Coulomb gauges.

C.1 Propagators
• The Feynman rule for an NRQED fermion propagator is • In Feynman gauge the Feynman rule for the photon propagator is q ! µ ⌫ −ig µν q 2 − λ 2 + i In Coulomb gauge the photon propagator is [37] D µν (q) = • There are two Feynman rules for the photon propagator in Coulomb gauge.

C.2 Vertices
The NRQED operators can be classified by the power of M by which they are suppressed. We need Feynman rules from operators suppressed by up to and including two powers of M .
• The operator ψ † iD t ψ gives rise to a one-photon Feynman rulẽ pp 0 − iZe.
• The operator ψ † c D e [∇ · E] 8M 2 ψ gives rise to a one-photon Feynman rulẽ In a non-Coulomb gauge, where ∇ · A = 0, there is an additional one-photon Feynman rule q " • The operator ψ † ic S e σ · (D × E − E × D) 8M 2 ψ gives rise to the one-photon Feynman rules