$\pi-\pi$ scattering, QED and finite-volume quantization

Using the Coulomb gauge formulation of QED we present a lattice QCD procedure to calculate the $\pi^+\pi^+$ scattering phase shift including the effects of the Coulomb potential which appears in this formulation. The approach described here incorporates the effects of relativity and avoids finite-volume corrections that vanish as a power of the volume in which the lattice calculation is performed. This is the first step in developing a complete lattice QCD calculation of the electromagnetic and isospin-breaking light-quark mass contributions to $\varepsilon'$, the parameter describing direct CP violating effects in $K_L\to\pi\pi$ decay.


I. INTRODUCTION
The K → ππ decays are important to our understanding of CP violation from the Weak interaction. Indirect CP violation in the K L → ππ decay arises from the admixture of the CP-even combination of K 0 and K 0 mesons in the long-lived K L state. It is described by the parameter ε whose measured magnitude is 2.228(11) × 10 −3 . Direct CP violation arises from the CP-odd component of K L , which can also directly decay into two pions. It is described by the parameter ε whose value is 3 orders of magnitude smaller than ε. Due to its small size, direct CP violation in K → ππ decay, represented by the parameter ε , is very sensitive to new physics.
The experimental measurement determines the direct CP-violating ratio Re(ε /ε) = 16.6(2.3) × 10 −4 [1,2]. The current standard model prediction for this quantity computed using lattice QCD is given by Re(ε /ε) = 21.7(2.6)(6.2)(5.0)×10 −4 [3]. Here the first error is statistical, the second is systematic and the third arises from the neglect of electromagnetism and the effects of the isospin-violating mass difference, m u − m d (which will be referred to below as m u − m d effects) that are the topic of this paper.
For most processes, electromagnetic corrections are on the order of the fine structure constant α = 1/137. However, the quantity ε /ε involves with equal weight the amplitudes A 0 and A 2 for the decay of a K 0 meson into two pions in the isospin 0 and isospin 2 ππ states, respectively. Since A 2 is suppressed relative to A 0 by a factor of 22, a feature called the ∆I = 1/2 rule, the electromagnetic modifications to A 0 can in principle induce corrections to A 2 which are 22 times larger than this usual 1/137 scale. Such effects can then propagate into the standard model prediction for Re(ε /ε).
These electromagnetic effects have been extensively studied using chiral perturbation theory [4][5][6][7][8][9]. In fact, it was the recent estimate of these isospin-breaking effects given in Ref. [9] that was used in Ref. [3] and appears above as an estimate of the error resulting from the neglect of these effects. Given the difficulty of this calculation and the large size of these corrections to Re(ε /ε), an ab initio calculation of these effects using lattice QCD would be of value. In this paper we outline a possible strategy for including electromagnetic and m u − m d effects in a lattice QCD calculation of ε and work out in detail the first step in this strategy. The approach presented here builds upon that proposed in Ref. [10]. See also the related paper of Cai and Daviodi [11].
There are a number of important challenges that must be addressed if electromagnetic and m u − m d effects are to be computed using lattice methods: 1. The calculation of K → ππ depends heavily on the finite-volume methods of Lüscher [12] and Lellouch and Lüscher [13] which rely on the exponentially-localized finite-range interactions of QCD. Adding electricity and magnetism (E&M) introduces long-range interactions, inconsistent with the Lellouch-Lüscher strategy. This problem is dramatically illustrated by the fact that the scattering phase shifts which play a central role in the finite-volume treatment of Lellouch and Lüscher are not even defined when E&M effects are included, with the long-distance wave functions acquiring phases which grow logarithmically with distance, ∼ η ln(kr) where k is the center-of-mass (CoM) momentum of the scattering particles, r the CoM separation of the outgoing particles and η = me 2 /(4πk) is the Sommerfeld parameter with ±e their charge and m their mass.
2. The usual treatment of K → ππ decay relies on isospin symmetry to distinguish two independent ππ final states, one with I = 0 and the other with I = 2. Electromagnetism and m u − m d effects break isospin symmetry, allowing the ππ states with I = 0 and I = 2 to mix. As a result, the final state scattering that is part of the K → ππ decay becomes a coupled, two-channel problem, requiring a more complex treatment of both the finite-volume eigenstates and the infinite-volume decay processes. 3. A process such as K → ππ decay, which involves the acceleration of charge, will contain well-known infrared singularities [14][15][16] which are removed by a careful treatment of the possible, near-degenerate final states, which include the intended ππ state as well as states with one or more emitted photons in addition to the two pions. While the effects of such soft radiation can be computed using standard methods for the case of the infinite-volume decay, possible photon emission in a finite-volume lattice calculation may introduce additional complications making the already challenging two-channel problem described above into a problem with four channels, two of which are three-particle channels.
We intend to carry out this calculation of the first-order E&M and m u − m d contributions to ε by exploiting the linear character of such a first-order calculation. We will separate the problem into simpler pieces which can be computed independently and then added together to obtain the final result. The first step in this divide-and-conquer approach treats the effects of electromagnetism in the Coulomb or radiation gauge. As is reviewed in the next section, in the Coulomb gauge, the E&M vector potential A is required to be transverse, ∇ · A = 0, and the resulting E&M Hamiltonian divides into two separate terms. The first is the familiar 1/| r | Coulomb potential between two charge density operators evaluated at equal times and separated by the displacement r. The second, independent piece involves the coupling of transverse photons to the spatial components of the E&M current operator.
These two components describe different physical phenomena in the rest frame of the kaon and their effects can be treated separately. Since, in a lattice calculation, we are fixed into a frame and will naturally work in a finite volume in which the kaon is at rest, no added difficulties are introduced by the choice of such a Lorentz non-covariant gauge.
The Coulomb potential gives the largest effect of electromagnetism when the charged pions are moving at non-relativistic velocities; effects which are familiar from non-relativistic quantum mechanics. Its long-range, 1/r character complicates the usual finite-volume methods of lattice QCD and introduces singular effects resulting from the scattering of charged particles at long distances, including the logarithms of r present in the asymptotic behavior of scattering solutions. It is this Coulomb potential component of the E&M problem that we will discuss in this paper.
The second component of the Coulomb-gauge Hamiltonian is the photon-current interaction which emits and absorbs massless photons. This component is the source of the usual infrared divergences but, as with the Coulomb interaction, also includes singular shortdistance effects which require careful treatment. This part of the problem is not treated in the present paper. However, for both the Coulomb and radiation parts, the long-distance effects which may be incompatible with the finite volumes used in lattice QCD are fundamentally classical. This suggests that these long-distance effects can be separated and computed analytically, leaving the portion of the problem to be treated using lattice methods free of these long distance difficulties.
For the case of the Coulomb potential, this separation of short-and long-distance effects is easily achieved by writing the Coulomb potential as the sum of two pieces: where we refer to the fixed separation R as the truncation radius and the two terms on the right-hand side of Eq. (1) as the short-distance truncated Coulomb potential V TC (r) (the left-hand term) and the long-distance compliment of the truncated Coulomb potential V TC (r) (the right-hand term). Given the linearity of the first order correction that we wish to compute, we can treat these two terms separately. The effects of V TC (r) can be directly combined with those of QCD in a lattice calculation provided we choose 2R to be smaller than the linear size of the volume used in this QCD+QED calculation. As we will show, the second term can be treated analytically if R is sufficiently large that terms which decrease exponentially for large R can be neglected.
In the present paper, we treat a portion of the Coulomb potential problem, studying π + π + scattering and derive results which can be combined with a lattice QCD calculation to obtain the π + π + scattering phase shift accurate to first order in α, including both the effects of QCD and the full Coulomb potential. As we will show, this combined lattice QCD and analytic calculation can be carried out with the only finite-volume errors being those which vanish exponentially in the linear extent of the volume used for the lattice QCD calculation. In a practical use of our result, there will be power-law corrections that result from the usual neglect of phase shifts for those partial waves with angular momentum l larger than some minimum value. The remaining part of the Coulomb problem needed for the calculation of ε requires examining the two-channel π + π − and π 0 π 0 scattering and the effects of the Coulomb potential on the actual K → ππ decay. This problem has been treated in the non-relativistic case in our earlier proceedings [10] and we plan to provide a solution to this problem including relativistic effects in a later paper.
The second part of the problem of computing the E&M and m u − m d contributions to ε requires determining the effects of the transverse radiation. In a fashion similar to our approach to the Coulomb potential, we plan to divide the transverse photons into two groups whose energies lie above or below a boundary energy E B . Those photons with energy E γ < E B are treated classically using the usual Bloch-Nordsiek methods while those whose energy E γ > E B will be treated explicitly using lattice QCD in a finite volume with linear extent L. The correctness of the classical treatment requires that E B /Λ QCD 1 so that structure-dependent effects can be neglected. At the same time, the accurate treatment of the hard radiation using finite-volume lattice QCD requires 1/L E B . These combined requirements will result in the neglect of corrections behaving as a power of 1/(Λ QCD L).
This treatment of the radiation part of the E&M problem is our long term strategy and the focus of current study.
Our use of the truncated Coulomb potential which allows us to include QED in a finitevolume lattice calculation is different from the more conventional approach of called QED L [17]. If specialized to the Coulomb potential alone, the QED L approach would express the 1/r Coulomb potential in a finite spatial volume of size L 3 as a conventional periodic Fourier series over wave numbers k = 2π(n 1 , n 2 , n 3 )/L for integers {n i } 1≤i≤3 from which the ill-defined term with k = (0, 0, 0) has been omitted. When beginning this project, we examined this approach to QED in a finite volume and present some of our results in appendices to this paper.
Here our results are closely related to those of Beane and Savage [18] and Beane et al. [19]. However, the treatment presented in Appendix A may provide a useful compliment to this earlier work since it does not involve an effective range approximation to the energy dependence of the scattering phase shift. We choose to use the truncated Coulomb potential because this position-space approach appears easier to understand and the absence of new power-law finite-volume corrections may be an important advantage in a calculation in which power-law finite-volume effects are being exploited to determine the scattering phase shifts. In Appendix B, we investigate the size of the power-law finite-volume corrections to the scattering phase shifts determined numerically in the case of non-relativistic quantum mechanics with a simple scattering potential and find a large 1/L correction which gives an easy-to-correct energy shift and a small upper bound on the 1/L 3 corrections, suggesting that the QED L approach may also work well for the problem at hand.
The paper is organized as follows. In Sec. II, we recall the Coulomb gauge treatment of QED and present the proposed method to compute the scattering of two identical, charged spin-zero particles including the combined effects of QCD and the instantaneous Coulomb potential where the latter is included only to first order in α. The usual finite-volume quantization method that can be used to determine the O(α) contributions of V TC to the π + π + scattering phase shift are described in Sec. III while in Sec. IV we present the details of an analytic calculation of the O(α) contributions of V TC to this phase shift. The analytic results of this section are expressed in terms of on-shell properties of QCD specified by the π + π + scattering phase shift in the absence of electromagnetic corrections and the π + electromagnetic form factor. Finally, conclusions and future plans are described in Sec. V.
The paper has two appendices, A and B, whose contents are described above.

II. FORMULATION AND STRATEGY
As described above, we propose to compute the QED corrections to π + π + scattering by using the Coulomb-gauge formulation of QED. In that approach the Minkowski-space QED Lagrangian is written as a standard textbook result [20] for the quantum treatment of the electromagnetic field.
Here j( r) and ρ( r) are the current and charge density operators for the quarks to which the E&M field couples. This is the treatment of QED that is used in the lattice calculation described in Sec. III. However, in Sec. IV where we study the E&M interations of pions at large distances, we will use "scalar" QED and j( r) and ρ( r) will be written in terms of the pion field and its spatial and temporal derivatives. Since we plan to compute E&M corrections to first order in α, the two interaction terms on the right-hand side of Eq. (2) can be treated independently and the results simply added together in the end. This will allow us to consider separately the Coulomb interaction with its long-range distortion of the two-particle scattering problem and the interaction j · A of the transverse photons which requires the Bloch-Nordsiek treatment [14]. In this paper we will focus on the corrections arising from the second term, the Coulomb interaction.
In our next step this Coulomb potential is itself further divided into two terms V TC and V TC : where contributions to V TC come only from separations | r − r | < R while V TC involves only separations | r − r | > R. Note, we will distinguish the operators V TC and V TC defined in Eqs. (3) and (4), which contain a factor of 1 2 and do not depend on r from the functions V TC (r) and V TC (r), defined following Eq. (1), which do depend on r and do not contain the factor of 1 2 .
As is discussed in greater detail in Sec. III, the finite-range V TC term can be used directly in a lattice calculation to determine its contribution to the π + π + scattering phase shift. Both the QCD interactions and those implied by V TC have a finite range and can be used with Luscher's finite-volume quantization condition to determine their combined contribution to the π + π + scattering phase shift. We can then Taylor expand this quantization condition in α. The α 0 term in this expansion gives the usual finite-volume result, determining the QCD scattering phase shift in terms of the computed finite-volume energy. The first-order term in this expansion of Luscher's quantization condition will determine the first order contribution to that phase shift coming from V TC in terms of the first-order shift in the finite-volume energy arising from the V TC term in the QCD + QED Hamiltonian, a quantity that can be directly computed from lattice QCD as is derived in Sec. III.
The contribution of V TC to the π + π + scattering phase shift appears inaccessible to the methods of lattice QCD because of its infinite range. However, as a result of the minimum separation R that enters V TC , the effects of V TC come from large distances where we will show that they can be calculated analytically in terms of the QCD π + π + scattering phase shift and the pion electromagnetic form factor if the center-of-mass energy is below the four-pion threshold.

III. NUMERICAL TREATMENT OF V TC
Because of the spatial cutoff at the distance R, the truncated Coulomb potential V TC can be used in a standard finite-volume lattice QCD calculation to determine its contribution to the π + π + scattering phase shift. As in lattice QCD calculations that include QED L , we can add V TC to QCD either perturbatively to some finite order in α or non-perturbatively including all orders in α. Since we wish to preserve the separation of the Coulomb interaction and the transverse radiation, it is natural to carry out this portion of the calculation to first order in α. We begin with Lüscher's finite-volume quantization condition [12]: where for simplicity we focus on the case of practical interest, that of s-wave scattering. If E is the energy of a two-pion state in a finite volume with sides of length L with periodic boundary conditions, then p = (E/2) 2 − m 2 must obey Eq. (5) for integer n. The function where and q = pL/(2π).
We then expand the quantized energy E = E (0) + αE (1) + . . . in a power series in α. If we perform a similar expansion of the phase shift δ 0 (p) = δ (0) 0 (p) + αδ (1) 0 (p) + . . . then Eq. (5) can also be expanded to relate the lattice result for E (1) to the desired order-α contribution of V TC to the scattering phase shift: where The first-order energy E (1) can be determined directly from a lattice calculation in the spatial volume L 3 . If O ππ (t) is a suitable ππ interpolating operator localized at the time t which is invariant under the allowed lattice translations and rotations, then using first-order perturbation theory, E (1) can be determined from the ratio of correlation functions: provided t f −t V and t V −t i are each sufficiently large that Eq. (9) is actually determining the matrix element of V TC between the ground state of the π + π + system and itself, while excited states are suppressed. Here we have modified the argument of V TC to insure the translational symmetry supported by the periodic boundary conditions imposed on the lattice volume in which it is being used: One additional minor complication that must be addressed in this calculation is the renormalization of the π + mass that results from the Coulomb interaction. This mass shift can be computed from a three-point function very similar to that appearing in Eq. (9) in which the π + π + interpolating operator O ππ is replaced by an interpolating operator for a single pion. This shift in mass can be eliminated by adding a first-order shift αm (1) to the quark mass. The quantity m (1) would be chosen so that when this order-α mass term is combined with the Coulomb potential the π + mass is not changed from its original value.
This first-order mass term should then be included in the calculation of E (1) described in Eq. (9) by simply adding it to the Coulomb potential operator in that equation. The resulting value for E (1) could then be used in Eq. (8) to determine the first-order E&M contribution to the scattering phase shift for two π + particles, each with a fixed physical mass, the same mass that was used in the original order-α 0 QCD calculation.

IV. ANALYTIC TREATMENT OF V TC
In this section we will calculate the contribution of V TC to the π + π + scattering phase shift δ l to first order in α in terms of two quantities: i) the phase shift δ l without QED corrections and ii) the electromagnetic form factor of the pion. This infinite-volume, orderα V TC correction can be added to the correction determined numerically to first order in V TC using lattice QCD as described in the previous section to obtain the entire contribution of the instantaneous Coulomb potential to the π + π + scattering phase shift. This combined result should have only finite-volume errors that fall exponentially as the volume grows with the exception of the power-law corrections that come from omitting the scattering from higher partial waves in usual treatment of finite-volume quantization and whatever approximations are introduced when determining the π + form factor.

A. General formulation
This analytic calculation can be performed by working with the usual relativistic Lipmann-Schwinger equation which expresses the full π + π + scattering amplitude as a sum of products of two-particle irreducible kernels connected by pairs of pion propagators as shown in Fig. 1. We will refer to this sum as the Lipmann-Schwinger series. As is conventional in such discussions we will treat the composite pion in QCD as an elementary particle in a relativistic (φ † φ) 2 field theory working to arbitrary order in a perturbation expansion in the (φ † φ) 2 interaction. This will result in a somewhat physical discussion of the pion structure, here introduced by the (φ † φ) 2 coupling, whose space-time scale is set by the pion mass. As we explore the approximations necessary for the validity of our approach, we will be able propagators, shown as a pion with a lightly shaded bubble representing a sum of all one-particle reducible graphs. This is the standard expression for the complete scattering amplitude M that is useful to discuss two-particle scattering below the four-particle threshold can include the effects of the Coulomb interaction V TC .
to estimate their size and appreciate the bounds on the scattering energy that must be imposed. We will then assume that the relations established are universal and will also hold true in the physical QCD problem provided that the spatial scale R appearing in V TC is larger than the distance scale of QCD so that the differences in structure between QCD and this (φ † φ) 2 model become irrelevant. This requirement of universal behavior, independent of quark structure, implies that the truncation radius obey R Λ QCD The on-shell center-of-mass ππ scattering amplitude M l (E) for angular momentum l that is determined by the graphs shown in Fig. 1 can be defined by Here the pion three-momenta p = pp and p = pp are proportional to the two unit vectors over whose directions we are integrating. The pion energies are given by where m is the pion mass and the total energy E in the center-of-mass system is given by The amplitude M p 4 , p 3 , p 2 , p 1 can be obtained directly from the usual connected timeordered product: which must be evaluated on-shell with (p i ) 4 = ω p for 1 ≤ i ≤ 4. For E < 4m the unitarity of the scattering matrix implies that M l (E) can be written as Except for a change in overall sign of the metric, these are similar to the conventions of Ref. [21]. We will use a Minkowski metric (1, 1, 1, −1) so that analytic continuation between Euclidean and Minkowski formulations can be accomplished by changes in phase of the energy or time arguments.
In fact, in the discussion below we will begin with a Euclidean-space version of Eq. (12) which defines the amplitude M E (p 4 , p 3 , p 2 , p 1 ) for Euclidean arguments. The amplitude M E is defined by an equation identical to Eq. (12) except that the time ordered product is to be computed using a Euclidean path integral or a Hilbert-space formalism in which the timedependent fields are defined using exponentials of the Hamiltonian without the usual factor of i. In addition, the four-vector dot products that appear in the Fourier transforms from position to momentum space in the Euclidean version of Eq. (12) use the Euclidean metric (1, 1, 1, 1). In the discussion below we will explicitly analytically continue the amplitude M E to obtain M .
In order to discuss this analytic continuation concretely, it is useful to remove the momentum-conserving delta function from Eq. (12), to write the equation in terms of the total incoming Euclidean four-momentum P and two relative four-momenta p and p which are defined by the rewritten version of Eq. (12): where we have used momentum conservation and translational invariance to reduce the number of four-vectors on which M E depends from four to three and to evaluate the connected four-point function at the location x 4 + x 3 = 0. In the center-of-mass system among the three four-momenta p , p and P , only P has a non-zero energy component. In the standard approach to analytic continuation from Euclidean to Minkowski space the real Euclidean energy P 0 is changed in phase as: P 0 → e −iθ P 0 with θ increasing from zero to π/2. In order to avoid potential exponential divergence in the integrals over (x 1 ) 0 and (x 2 ) 0 , these integration variables are also rotated in phase in the opposite direction: This procedure then requires the time-dependent Euclidean Green's function to be analytically continued in the time.
Although this approach to analytic continuation will not be used below, this interpretation is useful to determine precisely what Minkowski amplitude results after this analytic continuation is performed. By design the resulting Minkkowski-space Fourier transform will contain the phase whose dependence on the total four-momentum P is With our conventions, this is the correct phase for an incoming, positive energy only if P 0 = −E where E = 2ω p . Thus, when performing this analytic continuation we will begin with P 0 = −E and carry out the phase rotation P 0 = −Ee −iθ , increasing θ from zero to π/2.
Since we wish to calculate the first order effects of V TC on the ππ phase shift, we will consider those terms in which V TC enters one of the components appearing in the Lipmann-Schwinger series shown in Fig. 1. The V TC interaction can enter in three ways as illustrated The second and third cases arise when V TC appears in the Bethe-Salpeter kernel and are distinguished by a property of the resulting kernel. Since the V TC interaction involves the fourth components of two electromagnetic currents at positions ( z 2 , t) and ( z 1 , t) multiplied by the function θ(| z 2 − z 1 | − R)/| z 2 − z 1 | we are representing this O(α) insertion graphically by two vertices at the positions ( z 2 , t) and ( z 1 , t) joined by a wavy line corresponding to this position-dependent function and referred to here as a "photon" line. If cutting this line separates the graph in the kernel K into two parts then the corresponding amplitude is of type (c). Otherwise it is of type (b). We will refer to these three cases as insertions of types (a), (b) and (c).
As we will see, the V TC contributions from graphs of type (a) and (b) are exponentially suppressed for large R. This will be demonstrated in the next section by an analytic continuation from Euclidean space in which the evaluations of the self-energy and Bethe-Salpeter kernel graphs remain in Eulidean space where the exponential suppression coming from the large space-like separation of the two vertices in V TC can be easily seen. Diagrams of type (c), which contain the Coulomb scattering of two far-separated pions, will therefore contain all of the relevant effects of V TC and are evaluated in Section IV C.
Corrections to the subdiagrams appearing in Fig. 1 that are first-order in V TC . Diagram (a) is the correction to the dressed pion propagator expressed as two full pion propagators connected by a sum of one-particle reducible graphs (with darker shading) which cannot be separated into two parts by cutting a single pion line. (Here two subgraphs joined by a photon line are viewed as connected.) Diagram (b) shows the correction to the two-particle irreducible kernel which cannot be divided into two disconnected parts if the V TC photon line is cut. Diagram (c) represents the class of two-particle irreducible diagrams that can be divided into two disconnected parts when the V TC photon line is cut.

B. Analytic continuation
In order to show that the contribution V TC to graphs of type (a) and (b) is exponentially suppressed for large R for the case of π + π + scattering at an energy below the four-pion threshold, we will first demonstrate that these self-energy and kernel subgraphs can be evaluated in Euclidean space even after the original Euclidean amplitude M E has been analytically continued to physical Minkowski energies. To show this, we propose a particular procedure to carry out this analytical continuation in which the external lines and the internal two-pion loop integrals shown explicitly in Fig. 1 are evaluated in momentum space and the loop integration contours distorted to avoid the singularities as this continuation is carried out. However, during this process the self-energy and kernel graphs will be evaluated as functions of position and those positions (and the amplitudes in which they appear) will remain in Euclidean space.
Anticipating this strategy we express both the self-energy and two-particle irreducible kernels as Fourier transforms of Euclidean, position-space functions: Both expressions are to be evaluated in Euclidean space and Fourier transformed with Euclidean four-momenta. Note we have used the translational invariance of both types of diagram to remove one of the space-time arguments from each Green's function.
The Green's functions appearing in Eqs. (15) and (16)  With this approach, it is straight-forward to perform the needed analytic continuation in the total incoming energy E from its Euclidean to its Minkowski value: replace E by e −iθ E and vary θ from zero to π/2. This complex variable E enters the amplitudes appearing in the Lipmann-Schwinger series in two places. First, it appears in the exponents written explicitly in Eqs. (15) and (16). Since the exponential function is analytic, the continuation can be easily performed. However, when E becomes imaginary, the Fourier integrals acquire an exponentially growing factor and we must demonstrate these integrals remain convergent.
Second, E appears in the internal propagators shown in Fig. 1. Here the variation in E from a real to an imaginary value will force the integration contour over the loop variable k 0 to be distorted to avoid the singularities which move as the phase of E changes. The motion of these singularities is shown by the dotted lines in Fig. 1. This distortion of the contours of the loop energies will require a further analytic continuation of the kernel and self-energy functions which can also be carried out since these energy variables also appear explicitly in analytic exponential functions. Of course, we must demonstrate that this further introduction of exponentially growing time dependence does not cause these internal Fourier transforms to diverge.
Thus, our demonstration that the contribution of V TC when appearing in diagrams of type (a) and (b) is exponentially suppressed as R increases proceeds in two steps. First we will show that the analytic continuation described above is possible, allowing the actual evaluations of the quantities D and K to be performed in Euclidean space. Second we examine the case when V TC appears in either of these Euclidean-space quantities and use their Euclidean-space character to show this exponential suppression as R increases.
To address the first step we examine the convergence of the integrals over time in Eqs. (15) and (16) Thus, the largest exponential growth than can result from the continued Fourier transform shown in Eq. (15) with the exponent P 0 2 + k 0 takes the form exp (E − m)x 0 . However, the Euclidean-space, one-particle-irreducible, self-energy graphs which enter Eq. (15) themselves decrease exponentially as the time separation x 0 increases. This decreasing behavior will be determined by the least massive intermediate state that can appear between the two interpolating operators φ(0) and φ † (x) in Eq. (15) which must be a three-pion state with energy no smaller than 3m. Thus, our analytic continuation strategy will fail when E − 4m > 0, reminding us of the known requirement that we stay below the four-pion threshold in our study of physical ππ scattering.
The behavior of the analytically continued expression given in Eq. (16) for the two-particle irreducible kernel is very similar. Now there are three time integrals which may diverge as the time separations among the arguments in Eq. (16) grow. In examining the behavior expected, we might distinguish the case where this kernel is the first or last factor in the Lipmann-Schwinger series from the case where it appears somewhere in the middle. Since we are working in the center-of-mass system, the variables k 0 or k 0 would be zero in the first case while they can have an imaginary part as large as E/2 − m in the second. Thus, we will consider the second case since it includes the first.
It is convenient to replace the variables x 2 and x 1 in Eq. (16) by the average and relative coordinates X = (x 2 +x 1 )/2 and x = (x 2 −x 1 ). Given the symmetry between the values ±y 0 and between the values ±x 0 we will chose both x 0 and y 0 to be positive. Further we will first consider the case where X 0 is negative and its magnitude is larger than (x 0 + y 0 )/2. This case is illustrated in Fig. (4). Because we are discussing the dependence of the two-particle irreducible kernal on its four time coordinates, such a one-dimensional figure is sufficient for a general discussion. An amputated pion line joins the two-particle-irreducible kernel K at each of the four times, ± y 0 2 and ± x 0 2 + X 0 which implies that at each of these times the number of pions much change by an odd number. If we include the fact the kernel must be both one-and two-particle irreducible, the allowed minimum number of pions appearing in the three intermediate states between these four times must be those shown in Fig. (4).
We conclude that the dependence of the integrand on the three times x 0 , y 0 and X 0 is given by The first three factors on the left-hand side come from the exponents in the Fourier transform appearing in Eq. (16) while the remaining three factors come from inserting three-pion states between φ( y 0 2 ) and φ(− y 0 2 ) and between φ † ( x 0 2 + X 0 ) and φ † (− x 0 2 + X 0 ) while a four-pion intermediate state must be inserted between φ(− y 0 2 ) and φ † ( x 0 2 + X 0 ). Inspecting the right-hand side of Eq. (17), we recognize that the integrand decreases if any of the three positive quantities y 0 , x 0 and −X 0 is increased provided E is below the four-pion threshold, E < 4m. If we treat y 0 and x 0 as fixed can then discuss the behavior of the integrand as X 0 increases from it assumed large negative value. As X 0 increases we reach the point where the lower end of the interval [ y 0 2 , − y 0 2 ] collides with the upper end of the interval [ x 0 2 + X 0 , − x 0 2 + X 0 ]. At this point the X 0 behavior of the integrand changes from exp{−X 0 (E − 4m)} to exp{−X 0 (E − 2m)}. Thus, the integrand now decreases as X 0 becomes more positive. Continuing to increase X 0 the next change in behavior occurs when one of the two intervals [ y 0 2 , − y 0 2 ] and [ x 0 2 + X 0 , − x 0 2 + X 0 ] lies within the other and the X 0 behavior becomes more rapidly decreasing exp{−X 0 E}. Next as X 0 increases further and these two intervals partially overlap and then separate the exponential decrease becomes steeper changing to exp{−X 0 (E + 2m)} and finally exp{−X 0 (E + 4m)}.
Thus the maximum is reached when − y 0 2 = x 0 2 + X 0 at which point the integrand behaves as exp{(y 0 + x 0 )(E − 4m)} insuring convergence for E < 4m and providing the analyticity needed for our explicit Euclidean-space evaluation at physical energies below the four-pion threshold. The case where − y 0 2 = x 0 2 + X 0 could have been anticipated as giving the largest contribution since it avoids the situation with a four-pion intermediate state, e.g. in Fig. 4 the case where the "4 pions" segment has zero length.
We conclude that the Minkowski-space ππ scattering amplitude can be determined from the Lipmann-Schwinger series in which the self-energy corrections and two-particle irreducible kernel are computed in Euclidean space. This implies that when V TC appears in those subgraphs, its contribution will be suppressed exponentially when R becomes large.
In each case, the two vertices introduced by V TC , separated by the distance R, must be joined by at least two Euclidean pion propagators so that the V TC contribution will decrease exponentially as R grows -the property that we wish to establish.
We should point out that an estimate of this exponential decrease given by exp{−2mR} coming simply from two sequences of single pion propagtors joining one V TC vertex to the other is an overestimate of the rate of decrease because of the exponential growth of the Fourier transform factors discussed above. For example, once the large spatial separation R has been introduced into the self-energy or kernel subdiagrams from V TC there will be an advantage to a diamond pattern, such as that shown

C. Lipman-Schwinger equation solution to first order in V TC
The result of the previous section shows that the contribution to π + π + scattering from V TC will come only from diagrams of type (c) in Fig. 2 once terms that are exponentially suppressed for large R have been neglected. In this section we will evaluate that contribution, in this same, exponentially-accurate approximation. This calculation receives contributions from the four diagrams shown in Fig. 6. While the diagrams shown in that figure involve the off-shell Bethe-Salpeter scattering kernel and off-shell pion electromagnetic vertex, the presence of V TC will be shown to restrict their evaluation to an on-shell, long-distance region allowing them to be evaluated directly in terms of infinite-volume scattering data and the on-shell pion form factor, both quantities that can be directly determined from lattice QCD or, if we choose, from a combination of experiment and dispersion theory/phenomenology.
The result will be an explicit formula determining the contribution of V TC to the π + π + scattering phase shift. This is the missing component that was omitted from the phase shift calculation using finite-volume lattice methods which include only V TC . We note that all of the variables and formulae in this and later sections are expressed using Minkowski space conventions.
This calculation can be cast into a suggestive form by expressing the first-order contribution of V TC to M l of Eq. (11) given by the sum of the four classes of diagram shown in Fig. 6 as a product of incoming and outgoing factors multiplied by the V TC scattering kernel: Here the V TC scattering kernel K TC (k , k, P ) is shown as diagram (c) in Fig. 2 and can be written in momentum space after the total momentum conserving delta function has been removed: FIG. 6. The four types of diagram that must be evaluated to determine the correction to the scattering phase shift arsing from V TC . Here the pion self-energy insertions on the internal lines are not shown. In this Section a self-energy insertion carrying four-momentum k is written as the product of a free-particle propagator and a function of k 2 which is absorbed into the off-shell scattering amplitude M .
The factor of two results from a combination of the four equivalent contractions contributing to K TC (k , k, P ) and the factor of 1 2 appearing in the operator V TC defined in Eq. (4). This expression is analogous to that appearing in Eq. (16) for the complete scattering kernel except we have used translation symmetry to fix the origin at the midpoint between the locations of the two charge density operators.
As suggested by Fig. 6, Ψ in/out lm (k, P ) can be separated into two components: The first component, ψ 0 lm (x, P ), is constructed from the plane-wave component of the incoming state which corresponds to no initial/final state scattering beyond that resulting from V TC . The second component includes all initial/final scattering and involves the full off-shell scattering amplitude M (E) defined in Eq. (12). We will now consider in turn the calculation of the contribution from ψ 0 lm (k, P ) and then ψ in/out lm (k, P ).
1. Plane-wave, ψ 0 lm (k, P ), contribution The plane wave part, ψ 0 lm (k, P ) can be deduced directly from Eqs. (11), (12) and (14): Here p = (E/2) 2 − m 2 and p = pp. Although ψ 0 lm (k, P ) can be easily determined as in Eq. (22), when we consider the second term, ψ in/out lm (k, P ), it will be simplest to directly evaluate its contribution to the product K TC Ψ in/out (k , P ). If this is done for ψ 0 lm (k, P ), we find: As in previous equations the two arguments of the pion vertex function Γ 0 (q , q) are fourvectors but here the incoming four-momenta are written out as an explicit combination ( q, q 0 ) of spatial, q, and temporal, q 0 , components.
An important simplifying step in the present discussion and that below considering ψ in/out lm (k, P ) involves the treatment of the vertex functions Γ 0 ( P 2 ± k , P 2 ± k). For both of these functions to be on-shell, four conditions must be obeyed: If these conditions are obeyed then each vertex function can be expressed in terms of a single covariant form factor F depending only on the momentum transfer (k − k) 2 = | k − k | 2 : We can then introduce a charge density function ρ(r) given by: Thus, if Eq. (25) is obeyed we can write Here the extra bar has been introduced to distinguish the real function ρ(r) from the operator ρ(r) introduced in Section II.
These on-shell conditions are obeyed for the case that ψ 0 appears in both the right and left factors. As we will see below for the case of ψ in/out , in the limit of large R the effect of the two single-pion poles in the integrand of the integrals over the four-vectors k and k in Eq. (18) is also to limit the contributing values of k and k to those obeying the conditions in Eq. (25). Thus, in the case where k and k are off-shell we will view Eq. (28) as an approximation where terms that vanish on-shell have been dropped. Since those unwanted terms will come with factors which cancel the poles which contribute to the large R limit, the neglect of these terms is justified, as we will see, if terms that fall exponentially for large R are neglected.
For the case at hand, if Eq. (28) is used to replace each of the two vertex functions in Eq. (24), that equation can be rewritten as: provided we assume that k can be treated as on-shell and substitute z = w − r 2 + r 1 This result for the plane-wave piece ψ 0 takes a recognizable form if we examine the case where the plane-wave factor ψ 0 also appears as the left-hand factor and use the standard decomposition of the plane wave e i w·( p− p ) in spherical harmonics and spherical Bessel functions. If we define M 00 TC,l as the resulting purely plane-wave contribution to M TC,l then following the conventions of Ref. [22] we find where to obtain Eq. (32) we have expanded Eq. (13) to first-order in α to relate the expression computed in Eq. (31) to the scattering phase shift δ TC0 l caused by V TC alone, assuming no strong interaction scattering (but including the effects of the pion form factor). If Eq. (32) is solved for δ TC0 l we obtain the standard relativistic Born approximation for the scattering of two identical mesons, each with a charge distribution ρ(| r|), by the potential V TC .
The second scattering term ψ in/out lm (k, P ) in Eq. (20) takes a similarly simple form for the case of K TC . We will examine the case of ψ in lm (k, P ), starting with the explicit definition of the Feynman amplitude implied by the left-hand factor in diagrams (a) or (d) of Fig. 6: where we begin with the four-momentum k and the scattering amplitude M E (k, p, P ), defined in Eq. (14), in Euclidean space. As in the previous examination of ψ 0 lm (k, P ) we study the combination of the kernel K TC (k , k, P ) with ψ in lm (k, P ) Note that the combinatoric factor of two that appeared as a prefactor in Eq. (23) is absent here because of an additional factor of 1 2 that must be introduced when K and M are combined. Since both K and M are defined as amputated Green's functions, they each contain a combinatoric factor of two that will appear only once for the contractions joining the K and M subdiagrams. The factor of 1 2 removes this extra factor. The next step is to perform the integral over k 0 along the contour C shown in Fig. 7.
Here we have reproduced Fig. 3 but adopted the more familiar Minkowski labeling where the abscissa is the real part of k 0 which is now a standard Minkowski energy while the ordinate is the imaginary part of k 0 . We have also shown the location of the k 0 branch points corresponding to the three-pion state with total spatial momentum k.

Re(k 0 )
Diagram showing the treatment of the k 0 contour of integration when evaluating the firstorder contribution of V TC in Eq. (34). This is a copy of Fig. 3 but with the integration variable k 0 shown as a Minkowski rather than a Euclidean energy. This change is accomplished by relabeling the axes used in Fig. 3. The familiar two-particle threshold occurs when as E increases, the pole at E 2 − ω k collides with the pole at − E 2 + ω k and when k = 0 at the endpoint of the integral over | k|. In this diagram we also show the location in the k 0 plane of the three-particle branch cuts that arise from the k 0 dependence of the self-energy and two-particle irreducible kernel. These are shown as horizontal dotted lines extending to the right from the points ± E 2 + 3ω k 3 and to the left from the points ± E 2 − 3ω k 3 , moved below or above the real line for clarity.
We will simplify the integral over k 0 in Eq. (35) by using Cauchy's theorem to move the contour in that figure to the right. As shown in Fig. 7 the new contour has two parts. The first is a closed, circular clockwise contour C 1 surrounding the pole at k 0 = ω k − E 2 . The second piece is a vertical contour C 2 parallel to the imaginary axis with a real part obeying We will evaluate the contribution of the contour C 1 and show that for E < 4m the contribution from the contour C 2 to Eq. (18) vanishes exponentially as R increases. The contribution from C 1 can be evaluated from Cauchy's theorem so that Eq. (35) becomes For the previous case of ψ 0 lm we were able to simplify this expression by introducing a charge density ρ( r) and expressing each of the vertex functions Γ 0 (±k + P 2 , ±k + P 2 ) as a Fourier transform of ρ( r). We can introduce the same simplification in this case as well but As a result, we will continue to work with Eq. (36) and use Cauchy's theorem a second time to evaluate the integral over the magnitude of k exploiting the limit that R (which provides a lower bound on | z|) is large. We begin by using polar coordinates to represent the vector k: k = kk, expressing k as the product of the magnitude k and the unit vectork which is determined by the polar coordinates θ and φ in the usual waŷ k = (sin θ cos φ, sin θ sin φ, cos θ).
Defined in this way the variable k enters the amplitudes Γ 0 and M E in the exponent of a factor that performs an exponentially convergent Fourier transform. Thus, we expect these amplitudes to be analytic functions of k in the region around the real line. We will express the function e i k· z in the usual series [22]: Here k = | k| and z = | z| where the context should prevent confusion with the four vectors k and z . Since R is large we will replace the spherical Bessel function j l (kz) by its asymptotic Next we recognize that the integral in the right hand term of this equation over the interval 0 ≤ k < ∞ is exactly the same as the integral of the left-hand term over the interval −∞ > k ≤ 0 allowing us to rewrite integral over k in Eq. (36) in polar coordinate but with the integral over k extending over the entire real line: The integrand in Eq. (40) contains poles at k = ± Neglecting the contribution from the shifted contour which falls exponentially with in-creasing R (and hence increasing | z|), we retain only the contribution from this pole: We should also observe that the contribution to the k 0 integral in Eq. (35) coming from the displaced contour C 2 which we have neglected has no singularity for real values of k if E < 4m and hence can contribute only to a term falling exponentially in R for large R, justifying its omission.
We have now established the result described in Section IV C 1 following Eq. (28): since the large R behavior comes from the on-shell form factors, we can neglect any off-shell contribution and replace Γ 0 (±k + P 2 , ±k + P 2 ) by the Fourier transform of a charge density as given in Eq. (28).
Once justified, this simplification is best made earlier in our derivation of the large R behavior of the K TC ψ in lm (k , P ). Thus, we return to Eq. (36) but now written in terms of the charge density: We then repeat the steps used to set | k| = | p| but will work with the combined plane wave factor that now appears in Eq. (42): expanding e i k· w in spherical harmonics, combining the terms in j l (kw) so that the k contour runs from −∞ to +∞, shifting the k contour into the positive imaginary half-plane and keeping the k = +| p| pole.
With this altered treatment the vector k appears only in the factor e i k· w and the scattering amplitude M E (k, p, P ). When we expand the factor e i k· w in spherical harmonics the angular integral over the direction of k then selects only the term where what was labeled l in Eq. (41) now must equal l, assuming that the scattering amplitude M E (k, p, P ) is rotationally symmetric. The resulting simpler expression becomes: We can now combine these results for K TC ψ (0) lm and K TC ψ in lm l m by adding the results given in Eqs. (29) and (44): The full result including the both the factors Ψ out (k , P ) on the left and Ψ in (k, P ) on the right as in Eq. (18) and choosing the first-order term in an expansion of Eq. (13) in powers of α is given by: Here we have replaced the leading asymptotic behavior of the spherical Bessel functions j l (pw) and n l (pw) by those functions without approximation to reinstate the complete result of our derivation, restoring terms which fall as powers of 1/R but which had been omitted for clarity when l > 0.
We have also introduced the factor exp{−µw} which depends on the screening mass µ into Eq. 48 to regulate the long-distance logarithmic singularity that would otherwise be present in the integral over w. This singularity is the manifestation in perturbation theory of the long-range character of the Coulomb potential which appears in exact solutions of the Schrodinger and Dirac equations as a common phase which depends logarithmically on the radial variable w. Since this is a common phase, equal for all partial waves, it does not enter most physical quantities which implies that for small µ, the µ-dependence introduced arbitrarily into Eq. 48 will also not appear in the final physical quantities that we wish to calculate. For example, the quantity η +− that measures CP violation in the neutral kaon decay into two charged pions [23] is a ratio of decay amplitudes into a common π + π − final state from which this phase will naturally cancel.
Equation (48) provides the analytic method to calculate the contribution of V TC to the π + π + scattering phase shift in terms of physical QCD properties: the ππ scattering phase shift δ l and the Fourier transform ρ(| r|) of the π + electromagnetic form factor F (q 2 ). The result applies to the fully relativistic case provided the center-of-mass energy E < 4m and that R is sufficiently large that omitted terms falling exponentially in R can be omitted. The ππ scattering phase shift δ l can be determined by the usual application of Lüscher's finitevolume quantization methods. The convolution of two factors of the charge density with the potential V TC can be most easily determined by expanding the function V TC ( w + r 2 − r 1 ) assuming | w| | r i |. The leading term is determined by the pion's charge and the 1/R 2 terms by its charge radius.

V. CONCLUSION
As a first step in calculating the electromagnetic and m u −m d contributions to the measure of direct CP violation ε we have presented in detail a method to determine the contribution of the Coulomb potential to the π + π + scattering phase shift. If the quantization of QED is carried out in Coulomb gauge with ∇ · A = 0, this is well-defined and will give the complete E&M contribution when the effects of transverse radiation are included.
The calculation of these Coulomb effects is itself separated into two parts: one in which the separation r of the two charge operators in the Coulomb energy is less than R (referred to as the truncated Coulomb potential, V TC ) and a second in which this separation is greater than R, V TC . The effects of V TC can be directly determined from a finite-volume, lattice QCD calculation while those of V TC can be obtained from an analytic expression which we derive.
Given the finite range of QCD and the truncated potential V TC their contributions to the π + π + scattering phase shift can be computed with the standard finite-volume methods In two appendices we discuss the use of the QED L approach to E&M corrections again considering the Coulomb potential but in the non-relativistic limit. We present a derivation of a QED L , finite-volume quantization condition which determines the E&M corrections to scattering phase shift that does not involve an effective range approximation and provides an alternative to the earlier treatment of Bean and Savage [18]. We also carry out a numerical study which suggests that this QED L approach gives accurate results in spite of the presence of uncontrolled 1/L 3 corrections.
Two important future steps are required for a complete lattice calculation of the electromagnetic and m u − m d contributions to ε . First the current Coulomb potential treatment must be generalized to the two-channel I = 0 and I = 2 ππ system and the K → ππ decay analyzed. This was outlined in Ref. [10] for the non-relativistic case. Second the contribution of transverse radiation must be included which requires the usual treatment of infrared divergences and the effects of three-particle ππγ states. Each is the subject of on-going study.

ACKNOWLEDGMENTS
We would like to thank our RBC and UKQCD collaborators for their ideas and assistance. The QED L finite-volume formulation of QED [17] is widely used in lattice QCD calculations to include the effects of E&M. In this approach the Coulomb potential or photon propagator is written as a Fourier series appropriate for a function defined in a finite volume and obeying periodic boundary conditions. The singular modes at zero wave number are discarded. The resulting Coulomb potential or photon propagator is convenient to use in such a volume but differs from the corresponding physical quantity by terms which behave as 1/L n for n ≥ 1. These power law imperfections in the QED L formulation lead to lattice results with finite-volume errors that also fall as powers of 1/L with the first few terms taking universal, point-charge values. However, those errors falling with higher powers of 1/L will be unknown and must be removed by studying multiple volumes and extrapolating L → ∞. For a recent reference on this topic see the paper of Davoudi, et al. [24].
For a problem of the sort being studied here in which the quantization of finite-volume energies plays an essential role, the need when using QED L to extrapolate L → ∞ creates potentially serious difficulties. This motivates the truncated Coulomb potential approach developed in this paper in which such new power-law finite-volume corrections are avoided.
Never-the-less, given the frequent use of QED L , we present in this appendix a non-relativistic derivation of a finite-volume quantization condition that might be used in a QED L calculation to determine the E&M corrections to the low-energy π + π + scattering phase shift.
In the following Appendix B, we investigate its accuracy in a numerical study of a simple potential model. For the example studied, we find accurate results for finite L without the need for an L → ∞ extrapolation.
The derivation given in this appendix can be viewed as an alternative to that presented by Beane and Savage [18] and used by Beane, et al. [19]. In contrast to Ref. [18], we express the quantization condition directly as an explicit formula for the resulting π + π + phase shift rather than as a condition on the scattering length and effective range, defined in a fashion appropriate for a Coulomb scattering problem, which would correspond to that phase shift.

Properties of the QED L approximation to the Coulomb potential
The QED L potential is simply the Fourier-transformed Coulomb potential in a finite volume with the zero-mode omitted: where p n = 2π L (n 1 , n 2 , n 3 ), with the integers n i , 1 ≤ i ≤ 3 giving the three components of the vector n and 0 = (0, 0, 0) is the zero mode.
We are interested in the l = 0 component of this potential near the origin and how it deviates from the Coulomb potential, e 2 /4πr, at short distances. This is most easily done by noting the similarities between the QED L potential and the periodic Helmholtz function G 0 , introduced in Lüscher's discussion [12] of finite-volume energy quantization: where (A3) Substituting the known spherical harmonic expansion of the periodic Helmholtz function and evaluating the limit of vanishing k, we obtain the standard result: where Z 00 is the standard zeta function: The difference between the Coulomb potential and the QED L potential is plotted in Fig. 8 for L = 4. Even when approaching the boundary of the box, about which the QED L potential is symmetric, the difference between the QED L and Coulomb potentials is well represented by the two leading correction terms in Eq. (A4) with the largest deviations seen in those directions where the greatest distance from the origin can be reached. However, small direction-dependent discrepancies can also be seen.
From Eq. (A4) we see that near the origin the QED L potential resembles the classical e 2 /4πr Coulomb potential plus additional power-law terms. The second term acts as an  energy shift whose effects in the problem at hand can be completely removed. This term corresponds to the universal 1/L term in QED self-energy calculations [24,25]. The third term can be recognized as the potential energy due to a uniform charge distribution making the total charge in the volume zero. This and higher-order unphysical terms will combine with the finite-range strong interactions to introduce "structure-dependent" errors which can be removed by an explicit L → ∞ extrapolation. These are related to the structuredependent 1/L 3 terms found in QED L self-energy calculations [24,25].
2. Lüscher finite-volume quantization extended to include QED L Following this brief description of the properties of the V QED L (r) we will use this potential and derive an approximate quantization condition which relates the energy eigenvalues of two-particle eigenstates of the combined QCD and QED L Hamiltonian in a periodic volume of side L to the infinite-volume scattering phase shift that results from their combined effects.
In contrast with the earlier sections of this paper, we consider the scattering of two scalar particles in the non-relativistic limit. In this case, our use of the Coulomb potential and neglect of transverse radiation can be justified as giving the leading order contribution in a non-relativistic expansion.
We follow the standard approach of treating an alternative non-relativistic problem in which the QCD interaction is replaced by a simple scattering potential V S (r) in nonrelativistic quantum mechanics. Since the relation obtained between finite-volume energies and the scattering phase shifts arises from the constraint that the large r behavior of the scattering solution obeys periodic boundary conditions in a finite volume, we expect this same relation to hold for QCD where the large r behavior of the scattering states depends on the scattering phase shifts in the same way.
The result described here is a generalization of that obtained by Lüscher to include the effects of the Coulomb potential and the approach used in its derivation follows closely that of Ref. [12]. This problem was studied earlier by Beane and Savage [18] who present a similar quantization condition specialized to the low energy region where an effective range approximation is valid. However, the method used by Beane and Savage is general and could be used to obtain the result presented here.
In fact, the derivation given here and that of Bean and Savage differ in the same way as do the treatment presented by Lüscher in Ref. [12] and that of van Baal [26]. While Lüscher's treatment matches the asymptotic form of the full "strong-interaction" solution to the behavior of the finite-volume Helmholtz equation solution outside the range of the strong-interaction potential, van Ball, following Huang and Yang [27], replaces that potential by a simpler point-like pseudo-potential creating a problem which can be solved explicitly in a finite, cubic volume. The resulting relationship between the quantized finite-volume energy in this psuedo-potential problem and the phase shift predicted by that pseudo-potential is precisely the usual one.
Adding the Coulomb potential significantly changes the quantized finite-volume energies because of the long range behavior of the Coulomb potential. Further, as discussed above, the Coulomb potential itself must be modified if it is to be consistent with the introduction of periodic boundary conditions. For simplicity, we will consider the case where only the s-wave phase shift that appears in the infinite-volume scattering problem defined by the potential V S (r) alone (without the Coulomb potential) is non-zero.
Specifically, we consider the finite-volume energy eigenstates ψ V E ( r) of the non-relativistic quantum mechanical Hamiltonian: where V S (r) is a finite-range potential which plays the role of the strong interaction between two identical spin-zero particles with reduced mass µ. We examine the behavior of the wave function ψ V E in a region where r lies outside of the region of radius r 0 in which V S (r) is nonzero, r > r 0 and, to reduce the (r/L) n errors introduced by using QED L , we also require r L where L is the length of the spatial side of the finite volume. Following Lüscher we consider two descriptions of the finite-volume energy eigenstate ψ V E ( r) in this region. The first is given by the Green's function G QED L ( r) which obeys a Helmholtz equation which now includes the QED L potential: If the function G QED L ( r) is expanded in spherical harmonics, those components with l > 0 will be regular in the entire volume V and should agree with corresponding terms in a partial wave expansion of ψ V E ( r) provided E = k 2 /2µ since we are assuming that for l > 0 the phase shifts δ l induced by V S (r) are zero. Thus, for l > 0 we expect that the contributions to G QED L ( r) with angular momentum l will solve Eq. (A6) throughout the volume V . The l = 0 component of G QED L ( r) will not solve Eq. (A6) with non-zero V S (r) and will be singular as r → 0. In this appendix, we will use G QED L evaluated to first order in α.
The second description of the finite-volume energy eigenstate ψ V E ( r) in the region r 0 ≤ r L is provided by the spherical Coulomb functions F l and G l which solve the Schrödinger equation including the physical Coulomb potential and replace the spherical Bessel functions j l (kr) and n l (kr) which appear in Lüscher's treatment. Here we assume that for r/L sufficiently small we can ignore the difference between the potential V QED L which appears in Eqs. (A6) and (A7) and the usual Coulomb potential which appears in the equations obeyed by F l and G l . These functions are known to all orders in a perturbation expansion in e 2 and can be easily simplified to a first-order form.
The spherical Coulomb functions appear in the r ≥ r 0 description of the infinite-volume scattering solutions ψ E (r) and obey Eq. (A6) if V QED L is replaced by e 2 /(4πr). Specifically if we consider an infinite-volume energy eigenstate ψ E,l (r) with energy E and orbital angular momentum l, for r ≥ r 0 it can be written: ψ E,l (r) = a l cosδ l F l (η, kr) + sinδ l G l (η, kr) .
where the Sommerfield parameter η = µe 2 /(4πk),δ l is the scattering phase shift for this combined Coulomb plus strong interaction problem and a l is a normalization factor. Recall that the large r behaviors of F l and G l are determined by their asymptotic forms: with the usual Γ function Thus, for large r ψ E,l (r) ∼ cosδ l F l (η, kr) + sinδ l G l (η, kr) = sin kr − η ln(2kr) − 1 2 πl + σ l +δ l .
The ln(2kr) dependence present in the large-r behavior of ψ E,l (r) appears in the functions F l (η, kr) and G l (η, kr) so thatδ l is a well-defined phase that is conventionally used to describe the remaining combined effects of the "strong" and Coulomb interactions. The phase shift δ l is the infinite-volume quantity that we would like to determine from a finite-volume calculation.
Our Coulomb finite-volume quantization condition is then obtained by requiring that In the remainder of the appendix, we will distinguish with primes the energy and momentum (E and k ) that appear in the infinite-volume, Coulomb-potential quantities from the energy and momentum (E and k) that appear in the corresponding finite-volume quantities.
The small r behavior of the spherical Coulomb functions F l and G l is well known [28] and for l = 0, ρ = k r and η = e 2 µ/(4πk ) can be written: ψ(z) = Γ (z)/Γ(z) is the digamma function and γ E = 0.57721 . . . is Euler's constant.
Next we examine the finite-volume Helmholtz equation solution G QED L ( r) and determine its behavior for small r. As a first step we rewrite the differential equation, Eq. (A7) obeyed by G QED L ( r) as the Lipman-Schwinger integral equation: where G 0 ( r) is the Helmholtz solution with e 2 = 0 given in Eq. (A3). Equation (A19) can be expanded through first order in e 2 giving Now we must determine the singular and regular parts of the right-hand side of Eq. (A20).
We will determine the singular parts first. This can be done in position space where we argue that for small r the finite-volume solution G 0 ( r, E) can be approximated by its infinitevolume counter part: The singular part of the zeroth-order term in Eq. (A20) is given immediately by Eq. (A22).
However, the singular part of the first-order term can also be found using this same expression for G ∞ 0 ( r) by evaluating: where the dots represents terms that are constants or of higher order as r → 0. Since we are interested in the singularity at r = 0, we need not be concerned about the large r region in the integral appearing above. This behavior can be controlled, for example, by giving k a small, positive imaginary part.
Finally we must determine the regular part, specifically the limit at r = 0 of the right hand side of Eq. (A20) after the singular part determined above has been removed. Now, the behavior of these functions in the entire periodic volume becomes important and the simple, infinite-volume formulae used in the above discussion of the logarithmic singularity are not adequate. Lüchser has done this for the G 0 ( r) term in Eq. (A20) and found: where q = kL/2π and the evaluation at s = 1 in Eq. (A5) is described in Ref. [12]. We must now make a similar evaluation of the right most term in Eq. (A20).
In analogy with Eq. (A30), we define the regular part of the e 2 correction to the function ln(1/kr).
The relatively mild logarithmic singularity allows a direct numerical evaluation Y (q).
Next, following Lüscher we introduce a Coulomb function φ C (E) such that the combination: contains the same ratio of regular to singular parts as we have found in G QED L (E). Equating the ratio of regular divided by singular parts in the limit of small r in Eqs. (A32) and (A20) we obtain: Recognizing that η with µe 2 /k we can simplify in Eq. (A33) by expanding to first order in η to obtain a final formula for φ C (E) to order e 2 : For simplicity, we have not expanded the relation between the primed and unprimed variables to first order in e 2 .
We can summarize the results of this appendix by stating the relation between the energy eigenvalue of a finite-volume energy eigenstate and the corresponding infinite-volume nonrelativisitic s-wave scattering phase shiftδ 0 (E ) given by the quantization condition: where E = E + ∆E where ∆E is given in Eq. (A15). Equation (A35) is accurate up to but not including terms of order (r 0 /L) 3 when the scattering potential V S vanishes outside the radius r 0 .

Appendix B: Numerical tests of QED L finite-volume quantization
After identifying the O( r 2 L 3 ) correction to the QED L potential, one could become concerned about the size of such corrections. In order to understand the efficacy of the QED Lmodified finite-volume quantization condition, it is prudent to work with a model system where the phase shift is known exactly. In such a system, one can attempt to discern the size of the effects of neglected e 4 terms and of the neglected inverse powers of L by comparing with the precisely known result. This additional step is useful before embarking on a costly lattice QCD calculation comparing the truncated Coulomb and QED L results in order to gauge expectations for accuracy. In this Appendix, we perform a numerical calculation of the phase shift in a model theory and determine the residual structure-dependent error introduced by the unphysical O(e 2 L −3 ) term.

δ-shell potential
Quantum mechanics with a δ-shell potential, is a sufficiently simple system to solve both analytically in infinite volume and numerically in a finite volume. The δ-shell potential has the feature necessary for application of Lüscher's quantization condition, as well as the QED L -modified finite-volume quantization condition, that there exists an exterior region of | r| > r 0 beyond which the interaction does not exist.
In the numerical finite-volume studies, the interaction range will be kept such that r 0 < L/2.
In infinite volume, within the regions interior and exterior to the δ-shell, the radial component of wave function will behave as a free particle and a solution to the spherical Bessel equation. In the exterior, the wave function will be a linear combination of the regular and irregular spherical Bessel functions ψ(r > r 0 ) = kAj l (kr) − kBn l (kr) , where phase shift is found through the ratio B/A = tan δ l . Since the wave function must remain finite, in the interior region the wave function will contain only the regular solution up to an arbitrary normalization.
For a δ-shell potential, the wave function must be continuous at the shell and its derivative must have a discontinuity proportional to the value of the wavefunction on the shell.
Applying these conditions, the phase shift is found to be where the first term in the denominator comes from the Wronskian of these functions For the s-wave phase shift, this equation reduces to the textbook result: . (B5)
The exterior and interior wavefunction are parameterized in the just as in (B2) and (B3) and are subject to the same boundary conditions at the radius of the δ-shell. Using the Wronskian of the Coulomb wavefunctions G l (z)F l (z) − G l (z)F l (z) = 1, the phase shift in the presence of the Coulomb interaction will be given by By comparing numerical results with this analytical result, the efficacy of the QED L -modified finite-volume quantization condition can be tested.

Computational strategy
In the remainder of this appendix we will calculate numerically the finite-volume energies for the combined δ-shell and QED L potentials, apply the QED L finite-volume quantization condition given in Eq. (A35) and then compare the results for the scattering phase shift with the exact infinite-volume values given by Eq. (B7). Our numerical approach is to determine the finite-volume energies and wave functions to order e 0 and then to use firstorder perturbation to determine the O(e 2 ) correction.
To solve the Schrödinger equation with δ-shell potential alone in finite volume, again the wave function is parameterized in the regions interior and exterior to the δ-shell. The wave function in the interior region is the same as given in Eq. (B3), but the exterior region must satisfy the periodic boundary conditions. The exterior wavefunction is proportional to the finite-volume Green's function given in Eq. (A3). The boundary conditions at the δ-shell leads to the quantization condition where G 0,0 (r) is the s-wave component of the zeroth-order Green's function G 0 ( r) defined in Eq. (A3): We will impose no conditions on higher partial waves with l ≥ 4 making the usual assumption that the δ-shell phase shifts for l ≥ 4 can be neglected so that the regular behavior of G 0,l (r) needs no modification. Thus, this approximation is made in both our quantization condition and our numerical work and hence is not tested in their comparison. The zeros of Eq. (B8) will be determined numerically. The numerical evaluation of G(r; k) is described in Appendix B 4 and uses the Ewald summation technique to dramatically improve the convergence of the sums in Eq. (A3). From those energy levels and the corresponding wavefunction, the QED L potential will be included using first-order perturbation theory by calculating ψ|V QED L |ψ numerically, as described in Sec. B 5.

Ewald summation
An efficient numerical evaluation of the Helmholtz Green's function in three dimensions is a non-trivial task. In its typical representation, Eq. (A3), as a sum in momentum space the series is slowly convergent. If one truncates the sum to only contain terms with p < p max , then the leading neglected terms will be proportional to p −2 max , but there will be approximately p 2 max of them. The large number of these terms will create substantial contributions even if their individual summands are small and will also increase the computational costs as p max is raised. A representation as a sum over the periodic images in position space can be found using Poisson's summation formula G 0 ( r; k) = 1 4π n e ik| r− nL| | r − nL| (B10) which is even less convergent.
Variants of the Helmholtz Green's function are necessary in a wide range of applications in physics and engineering, such as determining the wave function of electrons within a quasi-periodic crystal. The evaluation of these functions has been expedited with a method called "Ewald summation" [29]. Ewald summation breaks the series into two components which represent the "near" contributions, which come from contributions of periodic images close to r and are evaluated in the position space representation, and "far" contributions, which come from the more distant periodic images and are evaluated with the momentum space representation [29]. As will be seen, the rate of convergence will depend on the cutoff between these two regions. If the cutoff were taken to either extreme that all points were in one of these regions, the rate of convergence will blow up as anticipated.
To perform Ewald summation, the summand in Eq. (B10) is rewritten as an integral e ikr 4πr = 1 2π 3/2 Γ e −r 2 ξ 2 e k 2 /ξ 2 dξ , where Γ is a particular contour from 0 to ∞ with some constraints on how it approaches those limits [29]. First Γ must leave the origin with a negative angle φ in the limits 0 > φ > − π 4 , second it returns to the real axis at ξ 0 , and then finally, continues on the real axis to ∞.
This integral is then broken into two regions, from 0 to ξ 0 and from ξ 0 to ∞. The first of these regions will be evaluated in momentum space and the second in position space G 0 ( r; k) = G 1 ( r; k) + G 2 ( r; k) = 1 L 3 p e i p· r e k 2 −p 2 4ξ 0 k 2 − p 2 + −1 4π n Re e i| r+ nL|k erfc | r + nL|ξ 0 + ik 2ξ 0 .
These sums converge significantly faster than the original sum. The first converges as e −p 2 /4ξ 2 0 and the second as e −n 2 L 2 ξ 2 0 . The original momentum space sum is the special case of ξ 0 = 0 and the position space sum is reached as ξ 0 → ∞ which both demonstrate the slow convergence. Using the optimal value of ξ 0 = 3 1 4 π L [30], these sums can converge to sufficient accuracy with just a few terms.
These summation methods will also be necessary for numerically evaluating the QED L potential efficiently. Unlike in a lattice application, with a UV regulator, the standard summation definition QED L potential is also slowly convergent. For the numerical work, the QED L potential will be evaluated using Eq. (A2).

Numerical analysis of δ-shell potential
The numerical analysis of the δ-shell potential is performed for two case. In one case the constant V 0 = 10, making the effects of the potential large compared to the kinetic energies studied. For the second we choose V 0 = 1. The radius of the δ-shell is held fixed to 1. The size of the finite volume L is varied from 4 to 8 to generate many energy levels.
The energy levels are calculated by finding the lowest zeros of Eq. B8 over a range of volumes. The phase shift without Coulomb interactions, determined through Lüscher quantization, is shown in Fig. 9. In both cases, the finite-volume quantization condition reproduces the true phase shift, shown as the curve in that figure, to high accuracy over the range of k studied. The largest discrepancies are O(10 −11 ). For the case of V 0 = 10, k cot δ, being linear in k 2 can be well described by the first two terms of the effective range expansion, but for the weaker potential more terms would be required.
Perturbation theory is used to find the energy shifts caused by the inclusion of the QED L potential. The unnormalized wavefunction corresponding to a finite-volume energy eigenstate for the δ-shell potential alone is approximately and then performing a Monte-Carlo integration for the exterior region. The value of the coupling is set to e 2 = 1 2µ , such that the effective Bohr radius, 4π/(e 2 µ) is significantly larger than the size of the box. In the region of small k, the Sommerfeld parameter η will begin to grow large and the perturbation theory used in the quantization condition begins to break down. Fig. 10 shows the results of the phase shifts from O(e 2 ) energy levels using the finite-volume quantization condition of Eq. (A35). Even though the phase shift changes quite dramatically from the neutral case, the vast majority of the phase shifts determined from the finite-volume quantization condition reproduce the infinite-volume result with subpercent level accuracy. The lowest k 2 for the weak interaction case of V 0 = 1 differ by O(1%) from the true infinite-volume non-perturbative result due to this breakdown of perturbation theory at large η. This study of a simple quantum mechanics problem implies that the structure-dependent effects of the unphysical L −3 introduced by the QED L potential on the phase shift are small for this case. Their contribution to the perturbative energy shift from the region within the radius of the strong interaction is, at their largest, O(1.2%) for L = 4 of the contribution from the true Coulomb interaction in the same region. Naïvely, one could expect that to