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

Norman Christ  Physics Department, Columbia University, New York City, New
York 10027, USA  Xu Feng  School of Physics, Peking University, Beijing
100871, China  Collaborative Innovation Center of Quantum Matter, Beijing
100871, China  Center for High Energy Physics, Peking University, Beijing
100871, China  State Key Laboratory of Nuclear Physics and Technology, Peking
University, Beijing 100871, China  Joseph Karpie  Physics Department, Columbia
University, New York City, New York 10027, USA  Tuan Nguyen  Physics
Department, Columbia University, New York City, New York 10027, USA

(October 17, 2021)

###### Abstract

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^{\prime}$, the parameter describing
direct CP violating effects in $K_{L}\to\pi\pi$ decay.

###### Contents

  1. I Introduction
  2. II Formulation and strategy
  3. III Numerical treatment of $V_{\rm TC}$
  4. IV Analytic treatment of $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$
    1. IV.1 General formulation
    2. IV.2 Analytic continuation
    3. IV.3 Lipman-Schwinger equation solution to first order in $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$
      1. IV.3.1 Plane-wave, $\psi^{0}_{lm}(k,P)$, contribution
      2. IV.3.2 Interacting $\psi^{\mathrm{in/out}}_{lm}(k,P)$ contribution
  5. V Conclusion
  6. A ${\rm QED}_{L}$ finite-volume quantization
    1. A.1 Properties of the ${\rm QED}_{L}$ approximation to the Coulomb potential
    2. A.2 Lüscher finite-volume quantization extended to include ${\rm QED}_{L}$
  7. B Numerical tests of ${\rm QED}_{L}$ finite-volume quantization
    1. B.1 $\delta$-shell potential
    2. B.2 Combined $\delta$-shell and Coulomb potential
    3. B.3 Computational strategy
    4. B.4 Ewald summation
    5. B.5 Numerical analysis of $\delta$-shell potential

##  I Introduction

The $K\to\pi\pi$ decays are important to our understanding of CP violation
from the Weak interaction. Indirect CP violation in the $K_{L}\to\pi\pi$ decay
arises from the admixture of the CP-even combination of $K^{0}$ and
$\overline{K^{0}}$ mesons in the long-lived $K_{L}$ state. It is described by
the parameter $\varepsilon$ whose measured magnitude is $2.228(11)\times
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
$\varepsilon^{\prime}$ whose value is 3 orders of magnitude smaller than
$\varepsilon$. Due to its small size, direct CP violation in $K\to\pi\pi$
decay, represented by the parameter $\varepsilon^{\prime}$, is very sensitive
to new physics.

The experimental measurement determines the direct CP-violating ratio
$\operatorname{Re}(\varepsilon^{\prime}/\varepsilon)=16.6(2.3)\times 10^{-4}$
Batley  _et al._ (2002); Abouzaid  _et al._ (2011). The current standard model
prediction for this quantity computed using lattice QCD is given by
$\operatorname{Re}(\varepsilon^{\prime}/\varepsilon)=21.7(2.6)(6.2)(5.0)\times
10^{-4}$ Abbott  _et al._ (2020). 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 $\alpha=1/137$. However, the quantity
$\varepsilon^{\prime}/\varepsilon$ 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 $\pi\pi$ states, respectively. Since $A_{2}$ is
suppressed relative to $A_{0}$ by a factor of 22, a feature called the $\Delta
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 $\operatorname{Re}(\varepsilon^{\prime}/\varepsilon)$.

These electromagnetic effects have been extensively studied using chiral
perturbation theory Cirigliano  _et al._ (2000a, b, c); Wolfe and Maltman
(2001); Cirigliano  _et al._ (2004, 2020). In fact, it was the recent estimate
of these isospin-breaking effects given in Ref. Cirigliano  _et al._ (2020)
that was used in Ref. Abbott  _et al._ (2020) 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
$\operatorname{Re}(\varepsilon^{\prime}/\varepsilon)$, 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 $\varepsilon^{\prime}$
and work out in detail the first step in this strategy. The approach presented
here builds upon that proposed in Ref. Christ and Feng (2018). See also the
related paper of Cai and Daviodi Cai and Davoudi (2018).

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

The calculation of $K\to\pi\pi$ depends heavily on the finite-volume methods
of Lüscher Luscher (1991) and Lellouch and Lüscher Lellouch and Luscher (2001)
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, $\sim\eta\ln(kr)$
where $k$ is the center-of-mass (CoM) momentum of the scattering particles,
$r$ the CoM separation of the outgoing particles and $\eta=me^{2}/(4\pi k)$ is
the Sommerfeld parameter with $\pm e$ their charge and $m$ their mass.

  2. 2.

The usual treatment of $K\to\pi\pi$ decay relies on isospin symmetry to
distinguish two independent $\pi\pi$ final states, one with $I=0$ and the
other with $I=2$. Electromagnetism and $m_{u}-m_{d}$ effects break isospin
symmetry, allowing the $\pi\pi$ states with $I=0$ and $I=2$ to mix. As a
result, the final state scattering that is part of the $K\to\pi\pi$ 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. 3.

A process such as $K\to\pi\pi$ decay, which involves the acceleration of
charge, will contain well-known infrared singularities Bloch and Nordsieck
(1937); Yennie  _et al._ (1961); Lee and Nauenberg (1964) which are removed by
a careful treatment of the possible, near-degenerate final states, which
include the intended $\pi\pi$ 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 $\varepsilon^{\prime}$ 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 $\vec{A}$ is required to be transverse,
$\vec{\nabla}\cdot\vec{A}=0$, and the resulting E&M Hamiltonian divides into
two separate terms. The first is the familiar $1/|\vec{r}\,|$ Coulomb
potential between two charge density operators evaluated at equal times and
separated by the displacement $\vec{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 short-distance 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:

$\frac{e^{2}}{4\pi|\vec{r}|}=\frac{e^{2}}{4\pi|\vec{r}|}\theta(R-r)+\frac{e^{2}}{4\pi|\vec{r}|}\theta(r-R)$
(1)

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_{\rm TC}(r)$ (the left-hand term) and the long-distance
compliment of the truncated Coulomb potential
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\rm 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 $\pi^{+}\pi^{+}$ scattering and derive results which can be combined
with a lattice QCD calculation to obtain the $\pi^{+}\pi^{+}$ scattering phase
shift accurate to first order in $\alpha$, 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
$\varepsilon^{\prime}$ requires examining the two-channel $\pi^{+}\pi^{-}$ and
$\pi^{0}\pi^{0}$ scattering and the effects of the Coulomb potential on the
actual $K\to\pi\pi$ decay. This problem has been treated in the non-
relativistic case in our earlier proceedings Christ and Feng (2018) 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 $\varepsilon^{\prime}$ 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_{\gamma}<E_{B}$ are treated classically using the usual Bloch-
Nordsiek methods while those whose energy $E_{\gamma}>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}/\Lambda_{\mathrm{QCD}}\ll 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\ll E_{B}$. These combined
requirements will result in the neglect of corrections behaving as a power of
$1/(\Lambda_{\mathrm{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
finite-volume lattice calculation is different from the more conventional
approach of called ${\rm QED}_{L}$ Hayakawa and Uno (2008). If specialized to
the Coulomb potential alone, the ${\rm 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
$\vec{k}=2\pi(n_{1},n_{2},n_{3})/L$ for integers $\\{n_{i}\\}_{1\leq i\leq 3}$
from which the ill-defined term with $\vec{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 Beane and
Savage (2014) and Beane et al. Beane  _et al._ (2021). 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 ${\rm 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 $\alpha$. The usual finite-volume quantization method that
can be used to determine the $O(\alpha)$ contributions of $V_{\rm TC}$ to the
$\pi^{+}\pi^{+}$ scattering phase shift are described in Sec. III while in
Sec. IV we present the details of an analytic calculation of the $O(\alpha)$
contributions of $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ to this
phase shift. The analytic results of this section are expressed in terms of
on-shell properties of QCD specified by the $\pi^{+}\pi^{+}$ scattering phase
shift in the absence of electromagnetic corrections and the $\pi^{+}$
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
$\pi^{+}\pi^{+}$ scattering by using the Coulomb-gauge formulation of QED. In
that approach the Minkowski-space QED Lagrangian is written as

$\displaystyle L_{\mathrm{EM}}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\int
d^{3}r\left\\{\left(\partial_{t}\vec{A}(\vec{r})\right)^{2}-\left(\vec{\nabla}\times\vec{A}(\vec{r})\right)^{2}+\vec{j}(\vec{r})\cdot\vec{A}(\vec{r})\right\\}$
$\displaystyle\hskip 108.405pt-\frac{1}{2}\int
d^{3}rd^{3}r^{\prime}\rho(\vec{r})\frac{1}{4\pi|\vec{r}-\vec{r}\,^{\prime}|}\rho(\vec{r}\,^{\prime}),$

a standard textbook result Srednicki (2007) for the quantum treatment of the
electromagnetic field. Here $\vec{j}(\vec{r})$ and $\rho(\vec{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 $\vec{j}(\vec{r})$
and $\rho(\vec{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 $\alpha$, the two interaction terms on the right-hand side of Eq.
(II) 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 $\vec{j}\cdot\vec{A}$ of the transverse photons which requires the
Bloch-Nordsiek treatment Bloch and Nordsieck (1937). 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_{\rm TC}$ and $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$:

$\displaystyle V_{\rm TC}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\int
d^{3}rd^{3}r^{\prime}\rho(\vec{r})\,\frac{\theta\bigl{(}R-|\vec{r}-\vec{r}\,^{\prime}|\bigr{)}}{4\pi|\vec{r}-\vec{r}\,^{\prime}|}\rho(\vec{r}\,^{\prime})$
(3) $\displaystyle V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\int
d^{3}rd^{3}r^{\prime}\rho(\vec{r})\,\frac{\theta\bigl{(}|\vec{r}-\vec{r}\,^{\prime}|-R\bigr{)}}{4\pi|\vec{r}-\vec{r}\,^{\prime}|}\rho(\vec{r}\,^{\prime})$
(4)

where contributions to $V_{\rm TC}$ come only from separations
$|\vec{r}-\vec{r}\,^{\prime}|<R$ while
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ involves only separations
$|\vec{r}-\vec{r}\,^{\prime}|>R$. Note, we will distinguish the operators
$V_{\rm TC}$ and $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ defined in
Eqs. (3) and (4), which contain a factor of $\frac{1}{2}$ and do not depend on
$r$ from the functions $V_{\rm TC}(r)$ and
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(r)$, defined following Eq.
(1), which do depend on $r$ and do not contain the factor of $\frac{1}{2}$.

As is discussed in greater detail in Sec. III, the finite-range $V_{\rm TC}$
term can be used directly in a lattice calculation to determine its
contribution to the $\pi^{+}\pi^{+}$ scattering phase shift. Both the QCD
interactions and those implied by $V_{\rm TC}$ have a finite range and can be
used with Luscher’s finite-volume quantization condition to determine their
combined contribution to the $\pi^{+}\pi^{+}$ scattering phase shift. We can
then Taylor expand this quantization condition in $\alpha$. The $\alpha^{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_{\rm
TC}$ in terms of the first-order shift in the finite-volume energy arising
from the $V_{\rm 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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ to the
$\pi^{+}\pi^{+}$ 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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$,
the effects of $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ come from
large distances where we will show that they can be calculated analytically in
terms of the QCD $\pi^{+}\pi^{+}$ 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_{\rm TC}$

Because of the spatial cutoff at the distance $R$, the truncated Coulomb
potential $V_{\rm TC}$ can be used in a standard finite-volume lattice QCD
calculation to determine its contribution to the $\pi^{+}\pi^{+}$ scattering
phase shift. As in lattice QCD calculations that include ${\rm QED}_{L}$ , we
can add $V_{\rm TC}$ to QCD either perturbatively to some finite order in
$\alpha$ or non-perturbatively including all orders in $\alpha$. 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 $\alpha$. We begin with Lüscher’s finite-volume quantization
condition Luscher (1991):

$\delta_{0}(p)+\phi(q)=n\pi,$ (5)

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=\sqrt{(E/2)^{2}-m^{2}}$ must obey Eq. (5) for integer $n$. The function
$\phi(q)$ obeys

$\tan\phi(q)=-\frac{\pi^{3/2}q}{\mathcal{Z}_{00}(1,q)}$ (6)

where

$\mathcal{Z}_{00}(s,q)=\frac{1}{\sqrt{4\pi}}\sum_{\vec{n}\in
Z^{3}}(\vec{n}\,^{2}-q^{2})^{-s}$ (7)

and $q=pL/(2\pi)$.

We then expand the quantized energy $E=E^{(0)}+\alpha E^{(1)}+\ldots$ in a
power series in $\alpha$. If we perform a similar expansion of the phase shift
$\delta_{0}(p)=\delta_{0}^{(0)}(p)+\alpha\delta_{0}^{(1)}(p)+\ldots$ then Eq.
(5) can also be expanded to relate the lattice result for $E^{(1)}$ to the
desired order-$\alpha$ contribution of $V_{\rm TC}$ to the scattering phase
shift:

$\delta_{0}^{(1)}(p^{(0)})=-\left\\{\frac{d\delta_{0}(p)}{dp}+\frac{d\phi(q)}{dq}\frac{L}{2\pi}\right\\}_{p=p^{(0)}}\frac{E^{(0)}}{4p^{(0)}}E^{(1)}$
(8)

where $p^{(0)}=\sqrt{(E^{(0)}/2)^{2}-m^{2}}$.

The first-order energy $E^{(1)}$ can be determined directly from a lattice
calculation in the spatial volume $L^{3}$. If $O_{\pi\pi}(t)$ is a suitable
$\pi\pi$ 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:

$\displaystyle E^{(1)}=\frac{\left\langle
O_{\pi\pi}^{\dagger}(t_{f})\frac{1}{2}\int_{V}\\!\\!\int_{V}d^{3}r_{2}d^{3}r_{1}\rho(\vec{r}_{2},t_{V})V_{\rm
TC}(|\vec{r}_{2}-\vec{r}_{1}|_{L})\rho(\vec{r}_{1},t_{V})O_{\pi\pi}(t_{i})\right\rangle}{\left\langle
O_{\pi\pi}(t_{f})O_{\pi\pi}(t_{i})\right\rangle}$ (9)

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_{\rm TC}$ between the
ground state of the $\pi^{+}\pi^{+}$ system and itself, while excited states
are suppressed. Here we have modified the argument of $V_{\rm TC}$ to insure
the translational symmetry supported by the periodic boundary conditions
imposed on the lattice volume in which it is being used:

$|\vec{r}_{2}-\vec{r}_{1}|_{L}=\left\\{\sum_{i=1}^{3}\left(\min\Bigl{[}\bigl{|}(r_{2})_{i}-(r_{1})_{i}\bigr{|},L-\bigl{|}(r_{2})_{i}-(r_{1})_{i}\bigr{|}\Bigr{]}\right)^{2}\right\\}^{\frac{1}{2}}.$
(10)

One additional minor complication that must be addressed in this calculation
is the renormalization of the $\pi^{+}$ 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 $\pi^{+}\pi^{+}$
interpolating operator $O_{\pi\pi}$ is replaced by an interpolating operator
for a single pion. This shift in mass can be eliminated by adding a first-
order shift $\alpha m^{(1)}$ to the quark mass. The quantity $m^{(1)}$ would
be chosen so that when this order-$\alpha$ mass term is combined with the
Coulomb potential the $\pi^{+}$ 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 $\pi^{+}$ particles, each with a fixed physical
mass, the same mass that was used in the original order-$\alpha^{0}$ QCD
calculation.

##  IV Analytic treatment of
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$

In this section we will calculate the contribution of
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ to the $\pi^{+}\pi^{+}$
scattering phase shift $\delta_{l}$ to first order in $\alpha$ in terms of two
quantities: i) the phase shift $\delta_{l}$ without QED corrections and ii)
the electromagnetic form factor of the pion. This infinite-volume,
order-$\alpha$ $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ correction
can be added to the correction determined numerically to first order in
$V_{\rm TC}$ using lattice QCD as described in the previous section to obtain
the entire contribution of the instantaneous Coulomb potential to the
$\pi^{+}\pi^{+}$ 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 $\pi^{+}$ form
factor.

###  IV.1 General formulation

This analytic calculation can be performed by working with the usual
relativistic Lipmann-Schwinger equation which expresses the full
$\pi^{+}\pi^{+}$ 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 $(\phi^{\dagger}\phi)^{2}$ field theory working to
arbitrary order in a perturbation expansion in the $(\phi^{\dagger}\phi)^{2}$
interaction. This will result in a somewhat physical discussion of the pion
structure, here introduced by the $(\phi^{\dagger}\phi)^{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 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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ is larger than
the distance scale of QCD so that the differences in structure between QCD and
this $(\phi^{\dagger}\phi)^{2}$ model become irrelevant. This requirement of
universal behavior, independent of quark structure, implies that the
truncation radius obey $R\gg\Lambda_{\mathrm{QCD}}$

Figure 1: A graphical expansion of the full scattering amplitude $M$ (with the
label $M$) in terms of products factors of the two-particle irreducible kernel
(with the label $K$) joined by dressed pion 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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$.

The on-shell center-of-mass $\pi\pi$ scattering amplitude $M_{l}(E)$ for
angular momentum $l$ that is determined by the graphs shown in Fig. 1 can be
defined by

$iM_{l}(E)\delta_{l^{\prime}l}\delta_{m^{\prime}m}=\frac{1}{4\pi}\int\\!\\!\int
d\Omega_{\hat{p}^{\prime}}d\Omega_{\hat{p}}\,Y^{*}_{l^{\prime}m^{\prime}}(\hat{p}^{\prime})M\bigl{(}(\vec{p}\,^{\prime},\omega_{p}),(-\vec{p}\,^{\prime},\omega_{p}),(\vec{p},\omega_{p}),(-\vec{p},\omega_{p})\bigr{)}Y_{lm}(\hat{p}).$
(11)

Here the pion three-momenta $\vec{p}=p\hat{p}$ and
$\vec{p}\,^{\prime}=p\hat{p}\,^{\prime}$ are proportional to the two unit
vectors over whose directions we are integrating. The pion energies are given
by $\omega_{p}=+\sqrt{p^{2}+m^{2}}$ where $m$ is the pion mass and the total
energy $E$ in the center-of-mass system is given by $E=2\omega_{q}$.

The amplitude $M\bigl{(}p_{4},p_{3},p_{2},p_{1}\bigr{)}$ can be obtained
directly from the usual connected time-ordered product:

$\displaystyle\prod_{i=1}^{4}\left\\{\frac{p_{i}^{2}+m^{2}}{i}\right\\}\prod_{i=1}^{4}\int\left\\{\int
d^{4}x_{i}\right\\}e^{i(-p_{4}x_{4}-p_{3}x_{3}+p_{2}x_{2}+p_{1}x_{1})}\bigl{\langle}\phi(x_{4})\phi(x_{3})\phi^{\dagger}(x_{2})\phi^{\dagger}(x_{1})\bigr{\rangle}_{\mathrm{conn}}$
$\displaystyle\hskip
108.405pt=(2\pi)^{4}\delta^{4}(p_{4}+p_{3}-p_{2}-p_{1})M(p_{4},p_{3},p_{2},p_{1}),$
(12)

which must be evaluated on-shell with $(p_{i})_{4}=\omega_{p}$ for $1\leq
i\leq 4$. For $E<4m$ the unitarity of the scattering matrix implies that
$M_{l}(E)$ can be written as

$M_{l}(E)=32\pi\frac{\omega_{p}}{p}\frac{e^{2i\delta_{l}}-1}{2i}.$ (13)

Except for a change in overall sign of the metric, these are similar to the
conventions of Ref. Christ  _et al._ (2015). 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 time-
dependent 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^{\prime}$ which are defined by the rewritten
version of Eq. (12):

$M_{E}(p^{\prime},p,P)=\int\\!\\!\int\\!\\!\int\\!\\!d^{4}yd^{4}x_{2}d^{4}x_{1}e^{-ip^{\prime}y}\left\langle\phi(\frac{y}{2})\phi(-\frac{y}{2})\phi^{\dagger}(x_{2})\phi^{\dagger}(x_{1})\right\rangle_{\mathrm{amp}}e^{i(\frac{P}{2}+p)x_{2}}e^{i(\frac{P}{2}-p)x_{1}},$
(14)

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^{\prime}$, $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}\to
e^{-i\theta}P_{0}$ with $\theta$ increasing from zero to $\pi/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: $(x_{i})_{0}\to e^{i\theta}(x_{i})_{0}$, $i=1,2$. This
procedure then requires the time-dependent Euclidean Green’s function
$\left\langle\phi(\frac{y}{2})\phi(-\frac{y}{2})\phi^{\dagger}(x_{2})\phi^{\dagger}(x_{1})\right\rangle$
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
$\frac{\vec{P}}{2}\cdot(\vec{x}_{1}+\vec{x}_{2})+\frac{P_{0}}{2}\bigl{(}(x_{1})_{0}+(x_{2})_{0}\bigr{)}$.
With our conventions, this is the correct phase for an incoming, positive
energy only if $P_{0}=-E$ where $E=2\omega_{p}$. Thus, when performing this
analytic continuation we will begin with $P_{0}=-E$ and carry out the phase
rotation $P_{0}=-Ee^{-i\theta}$, increasing $\theta$ from zero to $\pi/2$.

Since we wish to calculate the first order effects of
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ on the $\pi\pi$ phase
shift, we will consider those terms in which
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ enters one of the
components appearing in the Lipmann-Schwinger series shown in Fig. 1. The
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ interaction can enter in
three ways as illustrated Fig. 2. In the first case
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ contributes to the self-
energy graph appearing in one of the pion propagators connecting the two-
particle irreducible kernels as shown in graph (a) of Fig. 2. The second and
third cases arise when $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$
appears in the Bethe-Salpeter kernel and are distinguished by a property of
the resulting kernel. Since the
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ interaction involves the
fourth components of two electromagnetic currents at positions
$(\vec{z}_{2},t)$ and $(\vec{z}_{1},t)$ multiplied by the function
$\theta(|\vec{z}_{2}-\vec{z}_{1}|-R)/|\vec{z}_{2}-\vec{z}_{1}|$ we are
representing this $O(\alpha)$ insertion graphically by two vertices at the
positions $(\vec{z}_{2},t)$ and $(\vec{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).

Figure 2: Corrections to the subdiagrams appearing in Fig. 1 that are first-
order in $V_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ photon line is cut.

As we will see, the $V_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ and are evaluated in
Section IV.3.

###  IV.2 Analytic continuation

In order to show that the contribution
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ to graphs of type (a) and
(b) is exponentially suppressed for large $R$ for the case of $\pi^{+}\pi^{+}$
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:

$\displaystyle D(P,k)$ $\displaystyle=$ $\displaystyle\int
d^{4}x\bigl{\langle}\phi(0)\phi^{\dagger}(x)\bigr{\rangle}_{\mathrm{1PI}}e^{i(k+\frac{P}{2})x}$
(15) $\displaystyle K(P,k^{\prime},k)$ $\displaystyle=$
$\displaystyle\int\\!\\!\int\\!\\!\int d^{4}y\,d^{4}x_{2}\,d^{4}x_{1}$
$\displaystyle\hskip
54.2025pte^{-ik^{\prime}y}\left\langle\phi(\frac{y}{2})\phi(-\frac{y}{2})\phi^{\dagger}(x_{2})\phi^{\dagger}(x_{1})\right\rangle_{\mathrm{2PI,amp}}e^{i(k+\frac{P}{2})x_{2}}e^{i(-k+\frac{P}{2})x_{1}}.$

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 (IV.2) are intended to be
entirely general including both arbitrary orders in the
$(\phi^{\dagger}\phi)^{2}$ interaction and zeroth or first order in
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$. While
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ results from
electromagnetism, since it contains only operators that are evaluated at the
same time, it does not change the structure of the intermediate states that
determine the time dependence of these Green’s functions. As a result, the
possible presence of $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ in
these Green’s functions will be ignored when discussing their Euclidean time
dependence.

Figure 3: Diagram showing the treatment of the contour of $k_{0}$ integration
needed in the two-pion loop integrals that appear explicitly in the Lipmann-
Schwinger series, using the momentum assignment shown in Fig. 1. The change in
the incoming Euclidan energy from $E\to-i|E|$ causes the distortion of the
$k_{0}$ contour from the vertical axis to the new contour $C$. This diagram
refers to Euclidean-space coordinates but follows an unconventional assignment
of directions for the real and imaginary $k_{0}$ axes to give it a more
familiar, Minkowski-space orientation.

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\theta}E$ and vary $\theta$ from zero to
$\pi/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 (IV.2). 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_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{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 (IV.2) when the Euclidean energy appearing in the time-
dependent exponent has become imaginary. For the self-energy case described by
Eq. (15), we can see from Fig 3 that the magnitude of the imaginary part of
$k_{0}$ will be non-zero when $E/2>\omega_{k}$ and can be limited to
$E/2-\omega_{k}<E/2-m$. Thus, the largest exponential growth than can result
from the continued Fourier transform shown in Eq. (15) with the exponent
$\frac{P_{0}}{2}+k_{0}$ takes the form $\exp\bigl{\\{}(E-m)x_{0}\bigr{\\}}$.
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 $\phi(0)$ and $\phi^{\dagger}(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 $\pi\pi$
scattering.

The behavior of the analytically continued expression given in Eq. (IV.2) 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.
(IV.2) 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^{\prime}_{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. (IV.2) by
the average and relative coordinates $X=(x_{2}+x_{1})/2$ and
$x=(x_{2}-x_{1})$. Given the symmetry between the values $\pm y_{0}$ and
between the values $\pm 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,
$\pm\frac{y_{0}}{2}$ and $\pm\frac{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

$\left[e^{y_{0}(\frac{E}{2}-m)}e^{x_{0}(\frac{E}{2}-m)}e^{-X_{0}E}\right]\left[e^{-3m|y_{0}|}e^{-4m|-\frac{y_{0}}{2}-(\frac{x_{0}}{2}+X_{0})|}e^{-3m|x_{0}|}\right]=e^{\frac{x_{0}+y_{0}}{2}(E-4m)}e^{-X_{0}(E-4m)}.$
(17)

The first three factors on the left-hand side come from the exponents in the
Fourier transform appearing in Eq. (IV.2) while the remaining three factors
come from inserting three-pion states between $\phi(\frac{y_{0}}{2})$ and
$\phi(-\frac{y_{0}}{2})$ and between $\phi^{\dagger}(\frac{x_{0}}{2}+X_{0})$
and $\phi^{\dagger}(-\frac{x_{0}}{2}+X_{0})$ while a four-pion intermediate
state must be inserted between $\phi(-\frac{y_{0}}{2})$ and
$\phi^{\dagger}(\frac{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
$[\frac{y_{0}}{2},-\frac{y_{0}}{2}]$ collides with the upper end of the
interval $[\frac{x_{0}}{2}+X_{0},-\frac{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 $[\frac{y_{0}}{2},-\frac{y_{0}}{2}]$ and
$[\frac{x_{0}}{2}+X_{0},-\frac{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 $-\frac{y_{0}}{2}=\frac{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 $-\frac{y_{0}}{2}=\frac{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.

Figure 4: Diagram showing the ordering of arguments used in Eq. (17) with time
increasing from right to left. The average time of the two incoming vertices
is $X_{0}$ which is chosen negative while the earliest outgoing vertex (at the
time $-y_{0}/2$) is assumed to be later than the latest incoming vertex (at
the time $x_{0}/2+X_{0}$). The times of the four vertices are marked by
crosses while the average incoming and average outgoing times are each shown
as vertical lines. The minimum number of intermediate pions which appear in
the amplitude $K(P,k^{\prime},k)$ in each interval is also indicated.

We conclude that the Minkowski-space $\pi\pi$ 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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ appears in those subgraphs,
its contribution will be suppressed exponentially when $R$ becomes large. In
each case, the two vertices introduced by
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$, separated by the distance
$R$, must be joined by at least two Euclidean pion propagators so that the
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$
there will be an advantage to a diamond pattern, such as that shown Fig. 5 for
the case of a self-energy subdiagram, in which the two spatially-separated
vertices in $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ are joined to
two temporally-separated vertices. With this arrangement the exponential
suppression associated with the time separation of the two additional vertices
is partially compensated by the exponentially growing Fourier factors. For the
case shown in Fig. 5, such a diamond geometry leads to a weaken exponential
suppression. The total exponent corresponding to the four pion propagators
shown in Fig. 5 can be minimized by the choice of $z_{0}=x_{0}/2$ and a value
of $x_{0}$ which grows as $E$ approaches $4m$ giving a reduced exponential
fall off $\propto\exp\\{-R\sqrt{E(4m-E)}\\}$ which varies from the expected
$\exp\\{-2mR\\}$ when $E=2m$ to zero as $E$ increases to the four-pion
threshold.

Figure 5: An example of the insertion of the potential
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ into a self-energy graph
that can be used to estimate the rate of exponential decrease of such an
amplitude as the truncation radius $R$ increases. The location in space and
time of each of the four vertices is shown. As is discussed in the text the
contribution of graphs with this topology in the Lipman-Schwinger series will
decrease for large $R$ as $\exp\\{-\sqrt{E(4m-E)}R\\}$.

###  IV.3 Lipman-Schwinger equation solution to first order in
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$

The result of the previous section shows that the contribution to
$\pi^{+}\pi^{+}$ scattering from
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ to the
$\pi^{+}\pi^{+}$ 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_{\rm TC}$. We note that all of the variables and
formulae in this and later sections are expressed using Minkowski space
conventions.

Figure 6: The four types of diagram that must be evaluated to determine the
correction to the scattering phase shift arsing from
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{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$.

This calculation can be cast into a suggestive form by expressing the first-
order contribution of $V_{\overline{\raisebox{6.54413pt}{}\mathrm{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_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ scattering kernel:

${M_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}}_{,\,l}=\frac{1}{2l+1}\sum_{m=-l}^{l}\int\\!\\!\int
d^{4}k^{\prime}d^{4}k\,\Psi^{\mathrm{out}}_{lm}(k^{\prime},P)^{*}K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(k^{\prime},k,P)\Psi^{\mathrm{in}}_{lm}(k,P).$
(18)

Here the $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ scattering kernel
$K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(k^{\prime},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:

$\displaystyle
K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(k^{\prime},k,P)$
$\displaystyle=$ $\displaystyle-2\int\\!\\!\int d^{4}y_{2}\,d^{4}y_{1}\int
d^{3}z\int\\!\\!\int
d^{4}x_{2}\,d^{4}x_{1}\;e^{-i(k^{\prime}+\frac{P}{2})y_{2}}e^{-i(-k^{\prime}+\frac{P}{2})y_{1}}$
$\displaystyle\hskip
7.22743pt\Bigl{\langle}\phi(y_{2})\rho(\frac{z}{2})\phi^{\dagger}(x_{2})\Bigr{\rangle}_{\mathrm{1PI}}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(z)\Bigl{\langle}\phi(y_{1})\rho(-\frac{z}{2})\phi^{\dagger}(x_{1})\Bigr{\rangle}_{\mathrm{1PI}}e^{+i(k+\frac{P}{2})x_{2}}e^{i(-k+\frac{P}{2})x_{1}}.$

The factor of two results from a combination of the four equivalent
contractions contributing to
$K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(k^{\prime},k,P)$ and the
factor of $\frac{1}{2}$ appearing in the operator
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ defined in Eq. (4). This
expression is analogous to that appearing in Eq. (IV.2) 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, $\Psi^{\mathrm{in/out}}_{lm}(k,P)$ can be separated
into two components:

$\Psi^{\mathrm{in/out}}_{lm}(k,P)=\psi^{0}_{lm}(k,P)+\psi^{\mathrm{in/out}}_{lm}(k,P).$
(20)

The first component, $\psi^{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_{\overline{\raisebox{6.54413pt}{}\mathrm{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 $\psi^{0}_{lm}(k,P)$ and then
$\psi^{\mathrm{in/out}}_{lm}(k,P)$.

####  IV.3.1 Plane-wave, $\psi^{0}_{lm}(k,P)$, contribution

The plane wave part, $\psi^{0}_{lm}(k,P)$ can be deduced directly from Eqs.
(11), (12) and (14):

$\displaystyle\psi^{0}_{lm}(k,P)$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{4}}\int
d^{4}x\;e^{-i\vec{k}\cdot\vec{x}}\left\\{\frac{1}{\sqrt{4\pi}}\int
d\Omega_{\hat{p}}e^{i\vec{p}\cdot\vec{x}}Y_{lm}(\hat{p})\right\\}$ (21)
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{4\pi}p^{2}}\delta(|\vec{k}|-p)\delta(k_{0})Y_{lm}(\hat{k}).$
(22)

Here $p=\sqrt{(E/2)^{2}-m^{2}}$ and $\vec{p}=p\hat{p}$. Although
$\psi^{0}_{lm}(k,P)$ can be easily determined as in Eq. (22), when we consider
the second term, $\psi^{\mathrm{in/out}}_{lm}(k,P)$, it will be simplest to
directly evaluate its contribution to the product
$\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\Psi^{\mathrm{in/out}}\right](k^{\prime},P)$.
If this is done for $\psi^{0}_{lm}(k,P)$, we find:

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{0}_{lm}\right](k^{\prime},P)$
$\displaystyle=$ $\displaystyle-2\int d^{3}z\int
d^{3}k\Biggl{\\{}\Gamma_{0}\Bigl{(}k^{\prime}+\frac{P}{2},\bigl{(}\vec{k},\frac{E}{2}\bigr{)}\Bigr{)}e^{i(\vec{k}-\vec{k}^{\prime})\cdot\vec{z}}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{z}|\bigr{)}$
$\displaystyle\hskip
86.72377pt\cdot\Gamma_{0}\Bigl{(}-k^{\prime}+\frac{P}{2},\bigl{(}-\vec{k},\frac{E}{2}\bigr{)}\Bigr{)}\frac{1}{\sqrt{4\pi}p^{2}}\delta(|\vec{k}|-p)Y_{lm}(\hat{k})\Biggr{\\}}$
$\displaystyle=$ $\displaystyle-\frac{1}{\sqrt{\pi}}\int d^{3}z\int
d\Omega_{\hat{p}}\;\Gamma_{0}\Bigl{(}k^{\prime}+\frac{P}{2},\bigl{(}\vec{p},\frac{E}{2}\bigr{)}\Bigr{)}e^{i(\vec{p}-\vec{k}^{\prime})\cdot\vec{z}}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{z}|\bigr{)}$
$\displaystyle\hskip
86.72377pt\cdot\Gamma_{0}\Bigl{(}-k^{\prime}+\frac{P}{2},\bigl{(}-\vec{p},\frac{E}{2}\bigr{)}\Bigr{)}Y_{lm}(\hat{p}).$

As in previous equations the two arguments of the pion vertex function
$\Gamma_{0}(q^{\prime},q)$ are four-vectors but here the incoming four-momenta
are written out as an explicit combination $(\vec{q},q_{0})$ of spatial,
$\vec{q}$, and temporal, $q_{0}$, components.

An important simplifying step in the present discussion and that below
considering $\psi^{\mathrm{in/out}}_{lm}(k,P)$ involves the treatment of the
vertex functions $\Gamma_{0}(\frac{P}{2}\pm k^{\prime},\frac{P}{2}\pm k)$. For
both of these functions to be on-shell, four conditions must be obeyed:

$\displaystyle
k_{0}=k^{\prime}_{0}=0\quad\quad|\vec{k}|=|\vec{k}\,^{\prime}|=\sqrt{\left(\frac{E}{2}\right)^{2}-m^{2}}.$
(25)

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^{\prime}-k)^{2}=|\vec{k}-\vec{k}\,^{\prime}|^{2}$:

$\Gamma_{0}\left(\pm k^{\prime}+\frac{P}{2},\pm
k+\frac{P}{2}\right)=EF\bigl{(}(\vec{k}^{\prime}-\vec{k})^{2}\bigr{)}.$ (26)

We can then introduce a charge density function $\overline{\rho}(r)$ given by:

$\overline{\rho}(|\vec{r}|)=\frac{1}{(2\pi)^{3}}\int
d^{3}qe^{i\vec{q}\cdot\vec{r}}F(\vec{q}\,^{2}).$ (27)

Thus, if Eq. (25) is obeyed we can write

$\Gamma_{0}\left(\pm k^{\prime}+\frac{P}{2},\pm k+\frac{P}{2}\right)=E\int
d^{3}r\,e^{\pm(\vec{k}-\vec{k}\,^{\prime})\cdot\vec{r}}\overline{\rho}(r).$
(28)

Here the extra bar has been introduced to distinguish the real function
$\overline{\rho}(r)$ from the operator $\rho(r)$ introduced in Section II.

These on-shell conditions are obeyed for the case that $\psi^{0}$ appears in
both the right and left factors. As we will see below for the case of
$\psi^{\mathrm{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^{\prime}$ in Eq. (18) is also to limit the contributing values of $k$
and $k^{\prime}$ to those obeying the conditions in Eq. (25). Thus, in the
case where $k$ and $k^{\prime}$ 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. (IV.3.1), that equation can be rewritten as:

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{0}_{lm}\right](k^{\prime},P)$
$\displaystyle=$
$\displaystyle-\frac{E^{2}}{\sqrt{\pi}}\int\\!\\!\int\\!\\!\int\\!\\!\int
d^{3}w\,d\Omega_{\hat{p}}\,d^{3}r_{2}\,d^{3}r_{1}$ $\displaystyle\hskip
72.26999pt\overline{\rho}\bigl{(}|\vec{r}_{2}|\bigr{)}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{w}-\vec{r}_{2}+\vec{r}_{1}|\bigr{)}\overline{\rho}\bigl{(}|\vec{r}_{1}|\bigr{)}e^{i\vec{w}\cdot(\vec{p}-\vec{k}\,^{\prime})}Y_{lm}(\hat{p}),$

provided we assume that $k^{\prime}$ can be treated as on-shell and substitute
$\vec{z}=\vec{w}-\vec{r}_{2}+\vec{r}_{1}$

This result for the plane-wave piece $\psi^{0}$ takes a recognizable form if
we examine the case where the plane-wave factor $\psi^{0}$ also appears as the
left-hand factor and use the standard decomposition of the plane wave
$e^{i\vec{w}\cdot(\vec{p}-\vec{p}\,^{\prime})}$ in spherical harmonics and
spherical Bessel functions. If we define
${M_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}^{00}}_{,l}$ as the
resulting purely plane-wave contribution to
${M_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}}_{,l}$ then following the
conventions of Ref. Jackson (1998) we find

$\displaystyle{M_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}^{00}}_{,l}$
$\displaystyle=$ $\displaystyle\frac{1}{2l+1}\sum_{m=-l}^{l}\int\\!\\!\int
d^{4}k^{\prime}\,d^{4}k\,{\psi^{0}_{lm}(k^{\prime},P)}^{*}K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(k^{\prime},k,P)\psi^{0}_{lm}(k,P)$
(30) $\displaystyle=$ $\displaystyle-2E^{2}\int\\!\\!\int\\!\\!\int
d^{3}w\,d^{3}r_{2}\,d^{3}r_{1}\,j_{l}^{2}(pw)$ $\displaystyle\hskip
72.26999pt\cdot\overline{\rho}\bigl{(}|\vec{r}_{2}|\bigr{)}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{w}-\vec{r}_{2}+\vec{r}_{1}|\bigr{)}\overline{\rho}\bigl{(}|\vec{r}_{1}|\bigr{)}$
$\displaystyle=$ $\displaystyle
32\pi\frac{\omega_{p}}{p}\delta^{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}0}_{l},$
(32)

where to obtain Eq. (32) we have expanded Eq. (13) to first-order in $\alpha$
to relate the expression computed in Eq. (IV.3.1) to the scattering phase
shift $\delta^{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}0}_{l}$ caused by
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ alone, assuming no strong
interaction scattering (but including the effects of the pion form factor). If
Eq. (32) is solved for
$\delta^{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}0}_{l}$ we obtain

$\delta^{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}0}_{l}=-\frac{p\,\omega_{p}}{4\pi}\int\\!\\!\int\\!\\!\int
d^{3}w\,d^{3}r_{2}\,d^{3}r_{1}\,j_{l}^{2}(pw)\overline{\rho}\bigl{(}|\vec{r}_{2}|\bigr{)}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{w}-\vec{r}_{2}+\vec{r}_{1}|\bigr{)}\overline{\rho}\bigl{(}|\vec{r}_{1}|\bigr{)}.$
(33)

the standard relativistic Born approximation for the scattering of two
identical mesons, each with a charge distribution
$\overline{\rho}(|\vec{r}|)$, by the potential
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$.

####  IV.3.2 Interacting $\psi^{\mathrm{in/out}}_{lm}(k,P)$ contribution

The second scattering term $\psi^{\mathrm{in/out}}_{lm}(k,P)$ in Eq. (20)
takes a similarly simple form for the case of
$K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$. We will examine the case
of $\psi^{\mathrm{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:

$\displaystyle\psi^{\mathrm{in}}_{lm}(k,P)=-\frac{1}{\sqrt{4\pi}(2\pi)^{4}}\frac{1}{\bigl{[}(\frac{P}{2}+k)^{2}+m^{2}\bigr{]}\bigl{[}(\frac{P}{2}-k)^{2}+m^{2}\bigr{]}}\int
d\Omega_{\hat{p}}M(k,p,P)Y_{lm}(\hat{p}),$ (34)

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 $\psi^{0}_{lm}(k,P)$ we study the combination of the kernel
$K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(k^{\prime},k,P)$ with
$\psi^{\mathrm{in}}_{lm}(k,P)$

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right](k^{\prime},P)$
$\displaystyle=$ $\displaystyle\int\\!\\!\\!\int
d^{3}zd^{4}k\;\Gamma_{0}\Bigl{(}k^{\prime}+\frac{P}{2},k+\frac{P}{2}\Bigr{)}e^{i(\vec{k}-\vec{k}^{\prime})\cdot\vec{z}}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{z}|\bigr{)}\Gamma_{0}\Bigl{(}-k^{\prime}+\frac{P}{2},-k+\frac{P}{2}\Bigr{)}$
(35)
$\displaystyle\hskip-14.45377pt\cdot\frac{1}{\sqrt{4\pi}(2\pi)^{4}}\frac{1}{\bigl{[}(\frac{P}{2}+k)^{2}+m^{2}\bigr{]}\bigl{[}(\frac{P}{2}-k)^{2}+m^{2}\bigr{]}}\int
d\Omega_{\hat{p}}M(k,p,P)Y_{lm}(\hat{p}).$

Note that the combinatoric factor of two that appeared as a prefactor in Eq.
(IV.3.1) is absent here because of an additional factor of $\frac{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 $\frac{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 $\vec{k}$.

Figure 7: Diagram showing the treatment of the $k_{0}$ contour of integration
when evaluating the first-order contribution of
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{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 $\frac{E}{2}-\omega_{k}$ collides with the pole at
$-\frac{E}{2}+\omega_{k}$ and when $k=0$ at the endpoint of the integral over
$|\vec{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
$\pm\frac{E}{2}+3\omega_{\frac{k}{3}}$ and to the left from the points
$\pm\frac{E}{2}-3\omega_{\frac{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}=\omega_{k}-\frac{E}{2}$. The
second piece is a vertical contour $C_{2}$ parallel to the imaginary axis with
a real part obeying
$|-\frac{E}{2}+\omega_{k},\frac{E}{2}-\omega_{k}|<\mathrm{Re}(k_{0})<3m-\frac{E}{2}$.
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

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right](k^{\prime},P)$
$\displaystyle=$ $\displaystyle\frac{i}{\sqrt{4\pi}}\int\\!\\!\int
d^{3}z\frac{d^{3}k}{(2\pi)^{3}}$ $\displaystyle\hskip
21.68121pt\cdot\Biggl{\\{}\Gamma_{0}\Bigl{(}k^{\prime}+\frac{P}{2},k+\frac{P}{2}\Bigr{)}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{z}|\bigr{)}\Gamma_{0}\Bigl{(}-k^{\prime}+\frac{P}{2},-k+\frac{P}{2}\Bigr{)}$
$\displaystyle\hskip
43.36243pt\cdot\frac{e^{i(\vec{k}-\vec{k}^{\prime})\cdot\vec{z}}}{4\omega_{k}E\bigl{(}\sqrt{k^{2}+m^{2}}-\frac{E}{2})}\int
d\Omega_{\hat{p}}M_{E}\Bigl{(}k,p,P\Bigr{)}Y_{lm}(\hat{p})\Biggr{\\}}_{k_{0}=\omega_{k}-\frac{E}{2}}.$

For the previous case of $\psi^{0}_{lm}$ we were able to simplify this
expression by introducing a charge density $\overline{\rho}(\vec{r})$ and
expressing each of the vertex functions $\Gamma_{0}(\pm
k^{\prime}+\frac{P}{2},\pm k+\frac{P}{2})$ as a Fourier transform of
$\overline{\rho}(\vec{r})$. We can introduce the same simplification in this
case as well but first we must show that for large $R$ both of the vertex
functions in Eq. (IV.3.2) can be replaced by their on-shell values. This
simplification is not immediately possible here because the integration over
$k_{0}$ carried out when deriving Eq. (IV.3.2) evaluated only the vertex
function $\Gamma_{0}(k^{\prime}+\frac{P}{2},k+\frac{P}{2})$ on-shell. For the
other vertex function the time component of the incoming four-momentum is
$-k_{0}+\frac{E}{2}=E-\omega_{k}$. Thus, some added steps are needed to show
that both vertex functions can be replaced by their on-shell values.

As a result, we will continue to work with Eq. (IV.3.2) and use Cauchy’s
theorem a second time to evaluate the integral over the magnitude of $\vec{k}$
exploiting the limit that R (which provides a lower bound on $|\vec{z}|$) is
large. We begin by using polar coordinates to represent the vector $\vec{k}$:
$\vec{k}=k\hat{k}$, expressing $\vec{k}$ as the product of the magnitude $k$
and the unit vector $\hat{k}$ which is determined by the polar coordinates
$\theta$ and $\phi$ in the usual way
$\hat{k}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$.

Defined in this way the variable $k$ enters the amplitudes $\Gamma_{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\vec{k}\cdot\vec{z}}$ in the usual series Jackson (1998):

$e^{i\vec{k}\cdot\vec{z}}=\sum_{l^{\prime},m^{\prime}}4\pi
i^{l^{\prime}}j_{l^{\prime}}(kz)Y_{l^{\prime}m^{\prime}}(\hat{z})Y_{l^{\prime}m^{\prime}}^{*}(\hat{k}).$
(37)

Here $k=|\vec{k}|$ and $z=|\vec{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^{\prime}}(kz)$ by its asymptotic
form

$j_{l^{\prime}}(kz)\asymp\frac{e^{i(kz-\frac{\pi
l^{\prime}}{2})}-e^{-i(kz-\frac{\pi l^{\prime}}{2})}}{2ikz}.$ (38)

Next we recognize that the integral in the right hand term of this equation
over the interval $0\leq k<\infty$ is exactly the same as the integral of the
left-hand term over the interval $-\infty>k\leq 0$ allowing us to rewrite
integral over $\vec{k}$ in Eq. (IV.3.2) in polar coordinate but with the
integral over $k$ extending over the entire real line:

$\displaystyle\int d^{3}k\frac{e^{i(kr-\frac{\pi l}{2})}-e^{-i(kr-\frac{\pi
l}{2})}}{2ikr}\to\int
d\Omega_{\hat{k}}\int_{-\infty}^{\infty}k^{2}dk\frac{e^{i(kr-\frac{\pi
l}{2})}}{2ikr}$ (39)

This transformation of the integral of the second term on the left-hand side
of this equation to an integral over negative values of $k$ can be justified
by performing a simultaneous change in the sign of the variable $k$ and the
sign of the direction of $\hat{k}$ by performing a parity transformation on
the integration variable $\hat{k}$. This has no effect on the variable
$\vec{k}=k\hat{k}$ appearing in Eq. (IV.3.2). However, in the product
$j_{l^{\prime}}(kz)Y_{l^{\prime}m^{\prime}}^{*}(\hat{k})$ it changes the
$e^{-(ikz-\frac{\pi l^{\prime}}{2})}Y_{l^{\prime}m^{\prime}}^{*}(\hat{k})$
term into $e^{+(ikz-\frac{\pi
l^{\prime}}{2})}Y_{l^{\prime}m^{\prime}}^{*}(\hat{k})$.111This argument is
most easily presented using the asymptotic form of $j_{l^{\prime}}$. However,
it is also true if used to transform all of the terms in $j_{l^{\prime}}$ with
one sign of the exponent into the terms with the other sign, i.e. a
transformation between the two terms that appear when $j_{l^{\prime}}$ is
written in terms of spherical Hankel functions of the first and second kind.

With this change of coordinates the contribution to Eq. (IV.3.2) from the term
with specific values of $l^{\prime}$ and $m^{\prime}$ becomes:

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right]_{l^{\prime}m^{\prime}}(k^{\prime},P)$
$\displaystyle=$ $\displaystyle\sqrt{4\pi}i\int
d^{3}z\int_{-\infty}^{\infty}\frac{kdkd\Omega_{\hat{k}}}{(2\pi)^{3}}4\pi
e^{-i\vec{k}\,^{\prime}\cdot
z}Y_{l^{\prime}m^{\prime}}(\hat{z})Y_{l^{\prime}m^{\prime}}(\hat{k})^{*}$
$\displaystyle\hskip-50.58878pt\cdot\Biggl{\\{}\Gamma_{0}\Bigl{(}k^{\prime}+\frac{P}{2},k+\frac{P}{2}\Bigr{)}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{z}|\bigr{)}\Gamma_{0}\Bigl{(}-k^{\prime}+\frac{P}{2},-k+\frac{P}{2}\Bigr{)}$
$\displaystyle\hskip-28.90755pt\cdot\frac{e^{ik|z|}}{2i|\vec{z}|}\frac{1}{4\omega_{k}E\bigl{(}\sqrt{k^{2}+m^{2}}-\frac{E}{2})\bigr{)}}\int
d\Omega_{\hat{p}}M_{E}\Bigl{(}k,p,P\Bigr{)}Y_{lm}(\hat{p})\Biggr{\\}}_{k_{0}=\omega_{k}-\frac{E}{2}}.$

The integrand in Eq. (IV.3.2) contains poles at
$k=\pm\sqrt{\left(\frac{E}{2}\right)^{2}-m^{2}}\pm i\epsilon$. Because of the
positive imaginary coefficient of the variable $k$ in the exponent of the
integrand, we can shift the integration contour, adding a positive imaginary
constant. We can then write the expression in Eq. (IV.3.2) as the integral
over this shifted horizontal contour with a positive imaginary part added to
the Cauchy-theorem contribution from the pole at
$k=\sqrt{\left(\frac{E}{2}\right)^{2}-m^{2}}+i\epsilon$ that must be crossed
as the integration contour is shifted.

Neglecting the contribution from the shifted contour which falls exponentially
with increasing $R$ (and hence increasing $|\vec{z}|$), we retain only the
contribution from this pole:

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right]_{l^{\prime}m^{\prime}}(k^{\prime},P)$
$\displaystyle=$ $\displaystyle\sqrt{4\pi}\frac{i}{(2\pi)^{2}16\omega_{p}}\int
d^{3}z\int d\Omega_{\hat{k}}e^{-i\vec{k}\,^{\prime}\cdot
z}Y_{l^{\prime}m^{\prime}}(\hat{z})Y_{l^{\prime}m^{\prime}}(\hat{k})^{*}$
$\displaystyle\hskip
21.68121pt\cdot\Biggl{\\{}\Gamma_{0}\Bigl{(}k^{\prime}+\frac{P}{2},k+\frac{P}{2}\Bigr{)}V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{z}|\bigr{)}\Gamma_{0}\Bigl{(}-k^{\prime}+\frac{P}{2},-k+\frac{P}{2}\Bigr{)}$
$\displaystyle\hskip 43.36243pt\cdot\frac{e^{ik|\vec{z}|}}{|\vec{z}|}\int
d\Omega_{\hat{p}}M_{E}\Bigl{(}k,p,P\Bigr{)}Y_{lm}(\hat{p})\Biggr{\\}}_{\begin{subarray}{c}k=\sqrt{\left(\frac{E}{2}\right)^{2}-m^{2}}\\\
k_{0}=0\hfill\end{subarray}}.$

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.3.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 $\Gamma_{0}(\pm
k^{\prime}+\frac{P}{2},\pm k+\frac{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
$\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right](k^{\prime},P)$.
Thus, we return to Eq. (IV.3.2) but now written in terms of the charge
density:

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right](k^{\prime},P)$
$\displaystyle=$
$\displaystyle\frac{iE^{2}}{\sqrt{4\pi}}\int\\!\\!\int\\!\\!\int\\!\\!\int
d^{3}z\,d^{3}r_{2}\,d^{3}r_{1}\frac{d^{3}k}{(2\pi)^{3}}$ $\displaystyle\hskip
21.68121pt\cdot\Biggl{\\{}e^{i(\vec{k}-\vec{k}\,^{\prime})\cdot\vec{r}_{2}}\overline{\rho}(\vec{r}_{2})V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\bigl{(}|\vec{z}|\bigr{)}e^{-i(\vec{k}-\vec{k}\,^{\prime})\cdot\vec{r}_{1}}\overline{\rho}(\vec{r}_{1})$
$\displaystyle\hskip
43.36243pt\cdot\frac{e^{i(\vec{k}-\vec{k}^{\prime})\cdot\vec{z}}}{2\omega_{k}E\bigl{(}\sqrt{k^{2}+m^{2}}-\frac{E}{2})}\int
d\Omega_{\hat{p}}M_{E}\Bigl{(}k,p,P\Bigr{)}Y_{lm}(\hat{p})\Biggr{\\}}_{k_{0}=\omega_{k}-\frac{E}{2}}.$

We then repeat the steps used to set $|\vec{k}|=|\vec{p}|$ but will work with
the combined plane wave factor that now appears in Eq. (IV.3.2):

$e^{i(\vec{k}-\vec{k}\,^{\prime})\cdot(\vec{r}_{2}+\vec{z}-\vec{r}_{1})}\equiv
e^{i(\vec{k}-\vec{k}\,^{\prime})\cdot\vec{w}},$ (43)

expanding $e^{i\vec{k}\cdot\vec{w}}$ in spherical harmonics, combining the
terms in $j_{l}(kw)$ so that the $k$ contour runs from $-\infty$ to $+\infty$,
shifting the $k$ contour into the positive imaginary half-plane and keeping
the $k=+|\vec{p}|$ pole.

With this altered treatment the vector $\vec{k}$ appears only in the factor
$e^{i\vec{k}\cdot\vec{w}}$ and the scattering amplitude $M_{E}(k,p,P)$. When
we expand the factor $e^{i\vec{k}\cdot\vec{w}}$ in spherical harmonics the
angular integral over the direction of $\vec{k}$ then selects only the term
where what was labeled $l^{\prime}$ in Eq. (IV.3.2) now must equal $l$,
assuming that the scattering amplitude $M_{E}(k,p,P)$ is rotationally
symmetric. The resulting simpler expression becomes:

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right](k^{\prime},P)$
$\displaystyle=$ $\displaystyle\frac{iE}{4\sqrt{\pi}}\int
d^{3}w\,e^{-i\vec{k}\,^{\prime}\cdot\vec{w}}Y_{lm}(\hat{w})$
$\displaystyle\hskip 21.68121pt\cdot\,\frac{e^{ipw}}{w}\left\\{\int\\!\\!\int
d^{3}r_{2}\,d^{3}r_{1}\overline{\rho}(r_{2})V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(\vec{w}-\vec{r}_{2}+\vec{r}_{1})\overline{\rho}(r_{1})\right\\}M_{l}(E).$

We can now combine these results for
$\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{(0)}_{lm}\right]$
and
$\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\psi^{\mathrm{in}}_{lm}\right]_{l^{\prime}m^{\prime}}$
by adding the results given in Eqs. (IV.3.1) and (IV.3.2):

$\displaystyle\left[K_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}\Psi^{\mathrm{in}}_{lm}\right]_{l^{\prime}m^{\prime}}(k^{\prime},P)$
$\displaystyle=$ $\displaystyle-\int d^{3}we^{-i\vec{k}\,^{\prime}\cdot
w}Y_{lm}(\hat{w})$ $\displaystyle\hskip 21.68121pt\cdot\left\\{\int\\!\\!\int
d^{3}r_{2}\,d^{3}r_{1}\overline{\rho}(r_{2})V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(\vec{w}-\vec{r}_{2}+\vec{r}_{1})\overline{\rho}(r_{1})\right\\}$
$\displaystyle\hskip
43.36243pt\cdot\sqrt{4\pi}E^{2}\left\\{2\frac{i^{l}\sin(p|\vec{w}|-\frac{\pi
l}{2})}{p|\vec{w}|}+\frac{1}{16\pi\omega_{p}}\frac{e^{ip|\vec{w}|}}{|\vec{w}|}\left[32\pi\frac{\omega_{p}}{p}\frac{e^{2i\delta_{l}}-1}{2i}\right]\right\\}.$
$\displaystyle\hskip-93.95122pt=-4\sqrt{\pi}E^{2}\int
d^{3}we^{-i\vec{k}\,^{\prime}\cdot w}Y_{lm}(\hat{w})$
$\displaystyle\hskip-72.26999pt\cdot\left\\{\int\\!\\!\int
d^{3}r_{2}\,d^{3}r_{1}\overline{\rho}(r_{2})V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(\vec{w}-\vec{r}_{2}+\vec{r}_{1})\overline{\rho}(r_{1})\right\\}\frac{e^{i\delta_{l}}i^{l}\sin\bigl{(}p|\vec{w}|-\frac{\pi
l}{2}+\delta_{l}\bigr{)}}{p|\vec{w}|}.$

The full result including the both the factors
$\Psi^{\mathrm{out}}(k^{\prime},P)$ on the left and $\Psi^{\mathrm{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 $\alpha$ is given by:

$\displaystyle\delta^{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}_{l}$
$\displaystyle=$
$\displaystyle\frac{p}{32\pi\omega_{p}}{M_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}}_{,l}e^{-2\delta_{l}}$
(48) $\displaystyle=$
$\displaystyle-p\omega_{p}\int_{0}^{\infty}w^{2}dw\left\\{\int\\!\\!\int
d^{3}r_{2}\,d^{3}r_{1}\overline{\rho}(r_{2})V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(\vec{w}-\vec{r}_{2}+\vec{r}_{1})\overline{\rho}(r_{1})\right\\}$
$\displaystyle\hskip
36.135pt\cdot\Bigl{[}\cos{\delta_{l}}j_{l}(pw)+\sin(\delta_{l})n_{l}(pw)\Bigr{]}^{2}e^{-\mu
w}$

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\\{-\mu w\\}$ which depends on the
screening mass $\mu$ 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 $\mu$, the $\mu$-dependence introduced arbitrarily into
Eq. 48 will also not appear in the final physical quantities that we wish to
calculate. For example, the quantity $\eta_{+-}$ that measures CP violation in
the neutral kaon decay into two charged pions Donoghue  _et al._ (2014) is a
ratio of decay amplitudes into a common $\pi^{+}\pi^{-}$ final state from
which this phase will naturally cancel.

Equation (48) provides the analytic method to calculate the contribution of
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ to the $\pi^{+}\pi^{+}$
scattering phase shift in terms of physical QCD properties: the $\pi\pi$
scattering phase shift $\delta_{l}$ and the Fourier transform
$\overline{\rho}(|\vec{r}|)$ of the $\pi^{+}$ 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 $\pi\pi$ scattering
phase shift $\delta_{l}$ can be determined by the usual application of
Lüscher’s finite-volume quantization methods. The convolution of two factors
of the charge density with the potential
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ can be most easily
determined by expanding the function
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}(\vec{w}+\vec{r}_{2}-\vec{r}_{1})$
assuming $|\vec{w}|\gg|\vec{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 $\varepsilon^{\prime}$ we
have presented in detail a method to determine the contribution of the Coulomb
potential to the $\pi^{+}\pi^{+}$ scattering phase shift. If the quantization
of QED is carried out in Coulomb gauge with $\vec{\nabla}\cdot\vec{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_{\rm TC}$) and a second in which this separation is greater than $R$,
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$. The effects of $V_{\rm
TC}$ can be directly determined from a finite-volume, lattice QCD calculation
while those of $V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ can be
obtained from an analytic expression which we derive.

Given the finite range of QCD and the truncated potential $V_{\rm TC}$ their
contributions to the $\pi^{+}\pi^{+}$ scattering phase shift can be computed
with the standard finite-volume methods giving a result with finite-volume
errors which fall exponential in the spatial extent $L$ of the volume, with
the exception of power-law corrections asring from the usual neglect of
scattering phase shifts for angular momenta $l\geq 4$. Similarly, the analytic
expression given in Eq. (48) determines the contribution of
$V_{\overline{\raisebox{6.54413pt}{}\mathrm{TC}}}$ to exponential accuracy in
the truncation radius $R$ in terms of the $\pi^{+}\pi^{+}$ scattering phase
shift without E&M corrections and the pion form factor. This control of the
power-law finite-volume corrections may be important when exploiting finite-
volume quantization in which the results themselves are power-law, finite-
volume effects.

In two appendices we discuss the use of the ${\rm QED}_{L}$ approach to E&M
corrections again considering the Coulomb potential but in the non-
relativistic limit. We present a derivation of a ${\rm 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 Beane
and Savage (2014). We also carry out a numerical study which suggests that
this ${\rm 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 $\varepsilon^{\prime}$.
First the current Coulomb potential treatment must be generalized to the two-
channel $I=0$ and $I=2$ $\pi\pi$ system and the $K\to\pi\pi$ decay analyzed.
This was outlined in Ref. Christ and Feng (2018) 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 $\pi\pi\gamma$ 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 work of N.H.C., J.M.K. and T.H.N. was supported in part by
U.S. DOE grant #DE-SC0011941 while X.F. was supported in part by NSFC of China
under Grant No. 11775002 and No. 12070131001 and National Key Research and
Development Program of China under Contract No. 2020YFA0406400.

##  Appendix A ${\rm QED}_{L}$ finite-volume quantization

The ${\rm QED}_{L}$ finite-volume formulation of QED Hayakawa and Uno (2008)
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\geq 1$. These power law
imperfections in the ${\rm 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\to\infty$. For a recent reference on
this topic see the paper of Davoudi, et al. Davoudi  _et al._ (2019).

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 ${\rm
QED}_{L}$ to extrapolate $L\to\infty$ 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 ${\rm QED}_{L}$, we
present in this appendix a non-relativistic derivation of a finite-volume
quantization condition that might be used in a ${\rm QED}_{L}$ calculation to
determine the E&M corrections to the low-energy $\pi^{+}\pi^{+}$ 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\to\infty$
extrapolation.

The derivation given in this appendix can be viewed as an alternative to that
presented by Beane and Savage Beane and Savage (2014) and used by Beane, et
al. Beane  _et al._ (2021). In contrast to Ref. Beane and Savage (2014), we
express the quantization condition directly as an explicit formula for the
resulting $\pi^{+}\pi^{+}$ 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.

###  A.1 Properties of the ${\rm QED}_{L}$ approximation to the Coulomb
potential

The ${\rm QED}_{L}$ potential is simply the Fourier-transformed Coulomb
potential in a finite volume with the zero-mode omitted:

$V_{{\rm QED}_{\rm
L}}({\vec{r}})=\frac{e^{2}}{L^{3}}\sum_{{\vec{n}}\neq{\vec{0}}}\frac{e^{i{\vec{p}_{\vec{n}}}\cdot{\vec{r}}}}{p_{\vec{n}}^{2}}\,,$
(49)

where $\vec{p}_{\vec{n}}=\frac{2\pi}{L}(n_{1},n_{2},n_{3})$, with the integers
$n_{i}$, $1\leq i\leq 3$ giving the three components of the vector $\vec{n}$
and $\vec{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\pi r$, at short
distances. This is most easily done by noting the similarities between the
${\rm QED}_{L}$ potential and the periodic Helmholtz function $G_{0}$,
introduced in Lüscher’s discussion Luscher (1991) of finite-volume energy
quantization:

$V_{{\rm QED}_{\rm L}}({\vec{r}})=\lim_{k\to
0}-e^{2}\left[G_{0}({\vec{r}};k)-\frac{1}{L^{3}k^{2}}\right]\,,$ (50)

where

$G_{0}({\vec{r}};k)=\frac{1}{L^{3}}\sum_{\vec{n}}\frac{e^{i{\vec{p}_{\vec{n}}}\cdot{\vec{r}}}}{k^{2}-p_{\vec{n}}^{2}}\,.$
(51)

Substituting the known spherical harmonic expansion of the periodic Helmholtz
function and evaluating the limit of vanishing $k$, we obtain the standard
result:

$V_{{\rm QED}_{\rm
L}}({\vec{r}})=\frac{e^{2}}{4\pi}\left[\frac{1}{r}+\frac{Z_{00}(1;0)}{\pi^{3/2}(L/2\pi)}+\frac{2\pi
r^{2}}{3L^{3}}+\mathcal{O}\left(\frac{r^{4}}{L^{5}}Y_{l=4}({\hat{r}})\right)\right]\,,$
(52)

where $Z_{00}$ is the standard zeta function:

$Z_{00}(s,q^{2})=\frac{1}{\sqrt{4\pi}}\sum_{\vec{n}}\frac{1}{(n^{2}-q^{2})^{s}}.$
(53)

The difference between the Coulomb potential and the ${\rm QED}_{L}$ potential
is plotted in Fig. 8 for $L=4$. Even when approaching the boundary of the box,
about which the ${\rm QED}_{L}$ potential is symmetric, the difference between
the ${\rm QED}_{L}$ and Coulomb potentials is well represented by the two
leading correction terms in Eq. (52) 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.

Figure 8:  Left: The ${\rm QED}_{L}$ potential (solid black line), its
difference from the Coulomb potential (solid blue line), and the two leading
correction terms of Eq. (52) (dashed black line) in a box of size $L=4$.
Right: The ${\rm QED}_{L}$ potential (solid black line), its difference from
the Coulomb potential (solid blue line), the two leading correction terms in
Eq. (52) (dashed blue line), and the $r^{2}$ correction term alone (dashed
black line), each divided the Coulomb potential $V_{C}$. For the upper, middle
and lower rows $r$ is the distance from the origin along a line parallel to an
edge, a face diagonal or a body diagonal, respectively.

From Eq. (52) we see that near the origin the ${\rm QED}_{L}$ potential
resembles the classical $e^{2}/4\pi 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 Borsanyi  _et al._
(2015); Davoudi  _et al._ (2019). 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\to\infty$ extrapolation. These are
related to the structure-dependent $1/L^{3}$ terms found in ${\rm QED}_{L}$
self-energy calculations Borsanyi  _et al._ (2015); Davoudi  _et al._ (2019).

###  A.2 Lüscher finite-volume quantization extended to include ${\rm
QED}_{L}$

Following this brief description of the properties of the
$V_{\mathrm{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 ${\rm 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 non-relativistic 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. Luscher (1991). This problem was
studied earlier by Beane and Savage Beane and Savage (2014) 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. Luscher (1991) and
that of van Baal van Baal (1991). 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 Huang and Yang
(1957), 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
$\psi^{V}_{E}(\vec{r})$ of the non-relativistic quantum mechanical
Hamiltonian:

$\left(\frac{-\nabla^{2}}{2\mu}+V_{S}(r)+V_{\mathrm{QED}_{L}}(r)\right)\psi^{V}_{E}(\vec{r})=E\psi^{V}_{E}(\vec{r}),$
(54)

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
$\mu$. We examine the behavior of the wave function $\psi^{V}_{E}$ in a region
where $r$ lies outside of the region of radius $r_{0}$ in which $V_{S}(r)$ is
non-zero, $r>r_{0}$ and, to reduce the $(r/L)^{n}$ errors introduced by using
${\rm QED}_{L}$, we also require $r\ll 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 $\psi^{V}_{E}(\vec{r})$ in
this region.

The first is given by the Green’s function $G_{{\rm QED}_{L}}(\vec{r})$ which
obeys a Helmholtz equation which now includes the ${\rm QED}_{L}$ potential:

$\left(\nabla^{2}+k^{2}-2\mu V_{\mathrm{QED}_{L}}(r)\right)G_{{\rm
QED}_{L}}(\vec{r})=\delta(\vec{r})$ (55)

If the function $G_{{\rm QED}_{L}}(\vec{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
$\psi^{V}_{E}(\vec{r})$ provided $E=k^{2}/2\mu$ since we are assuming that for
$l>0$ the phase shifts $\delta_{l}$ induced by $V_{S}(r)$ are zero. Thus, for
$l>0$ we expect that the contributions to $G_{{\rm QED}_{L}}(\vec{r})$ with
angular momentum $l$ will solve Eq. (54) throughout the volume $V$. The $l=0$
component of $G_{{\rm QED}_{L}}(\vec{r})$ will not solve Eq. (54) with non-
zero $V_{S}(r)$ and will be singular as $r\to 0$. In this appendix, we will
use $G_{{\rm QED}_{L}}$ evaluated to first order in $\alpha$.

The second description of the finite-volume energy eigenstate
$\psi^{V}_{E}(\vec{r})$ in the region $r_{0}\leq r\ll L$ is provided by the
spherical Coulomb functions $\mathcal{F}_{l}$ and $\mathcal{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_{\mathrm{QED}_{L}}$ which appears in Eqs. (54) and (55) and the usual
Coulomb potential which appears in the equations obeyed by $\mathcal{F}_{l}$
and $\mathcal{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\geq r_{0}$ description of the
infinite-volume scattering solutions $\psi_{E}(r)$ and obey Eq. (54) if
$V_{\mathrm{QED}_{L}}$ is replaced by $e^{2}/(4\pi r)$. Specifically if we
consider an infinite-volume energy eigenstate $\psi_{E,l}(r)$ with energy $E$
and orbital angular momentum $l$, for $r\geq r_{0}$ it can be written:

$\psi_{E,l}(r)=a_{l}\left\\{\cos\bar{\delta}_{l}\,\mathcal{F}_{l}(\eta,kr)+\sin\bar{\delta}_{l}\,\mathcal{G}_{l}(\eta,kr)\right\\}.$
(56)

where the Sommerfield parameter $\eta=\mu e^{2}/(4\pi k)$, $\bar{\delta}_{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 $\mathcal{F}_{l}$ and $\mathcal{G}_{l}$ are determined
by their asymptotic forms:

$\displaystyle\mathcal{F}_{l}(\eta,\rho)$ $\displaystyle\sim$
$\displaystyle\sin\left(\rho-\eta\ln(2\rho)-\frac{1}{2}\pi
l+\sigma_{l}\right)$ (57) $\displaystyle\mathcal{G}_{l}(\eta,\rho)$
$\displaystyle\sim$ $\displaystyle\cos\left(\rho-\eta\ln(2\rho)-\frac{1}{2}\pi
l+\sigma_{l}\right)$ (58)

where

$\sigma_{l}=\mathrm{arg}\Bigl{(}\Gamma(l+1+i\eta)\Bigr{)}.$ (59)

with the usual $\Gamma$ function

$\Gamma(z)=\int_{0}^{\infty}dte^{-t}t^{z-1}.$ (60)

Thus, for large $r$

$\displaystyle\psi_{E,l}(r)$ $\displaystyle\sim$
$\displaystyle\cos{\bar{\delta}_{l}}\mathcal{F}_{l}(\eta,kr)+\sin{\bar{\delta}_{l}}\mathcal{G}_{l}(\eta,kr)$
(61) $\displaystyle=$ $\displaystyle\sin\left(kr-\eta\ln(2kr)-\frac{1}{2}\pi
l+\sigma_{l}+\bar{\delta}_{l}\right).$ (62)

The $\ln(2kr)$ dependence present in the large-$r$ behavior of $\psi_{E,l}(r)$
appears in the functions $\mathcal{F}_{l}(\eta,kr)$ and
$\mathcal{G}_{l}(\eta,kr)$ so that $\bar{\delta}_{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 $\bar{\delta}_{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 these two descriptions of the $s$-wave component of finite-volume
solution $\psi^{V}_{E}(\vec{r})$ agree in the region where both are valid:
$r_{0}\leq r\ll L$. If we neglect the differences between
$V_{\mathrm{QED}_{L}}$ and $e^{2}/(4\pi r)$, then these two representations of
the $s$-wave part of the finite-volume solution $\psi^{V}_{E}(\vec{r})$ must
be the same and the most accessible way to compare them is to extend both into
the region near $r=0$ where, since they are identical, they must have the same
ratio of regular to singular parts. Thus, we must examine the small $r$
behavior of the $l=0$ component of the Green’s function $G_{{\rm
QED}_{L}}(\vec{r})$ determined by Eq. (55) and the spherical Coulomb functions
$\mathcal{F}_{0}$ and $\mathcal{G}_{0}$ in Eq. (56). Equating the ratios of
the regular and singular parts in each description will then provide the
quantization condition determining the phase shift $\bar{\delta}_{0}$.

In Eq. (52) we show the $1/L$ and $1/L^{3}$ terms that appear for small $r/L$
in the difference between the Coulomb potential and $V_{\mathrm{QED}_{L}}$.
Both terms are potential sources of error in the phase shifts determined from
this quantization condition. While we do not attempt to characterize the
errors introduced by the $r^{2}/L^{3}$ term in Eq. (52), we can easily remove
the $1/L$ error by equating the ratio of the singular and regular $s$-wave
parts of the ${\rm QED}_{L}$ Helmholtz solution $G_{E}(\vec{r})$ evaluated at
the finite-volume energy $E$ with those of the spherical Coulomb functions
$\mathcal{F}_{0}$ and $\mathcal{G}_{0}$ combined with the phase shift
$\bar{\delta}$ as in Eq. (56) evaluated at the energy $E^{\prime}=E+\Delta E$
where

$\Delta E=-\frac{e^{2}Z_{00}(1;0)}{2\pi^{3/2}L}\,.$ (63)

In the remainder of the appendix, we will distinguish with primes the energy
and momentum ($E^{\prime}$ and $k^{\prime}$) 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 $\mathcal{F}_{l}$
and $\mathcal{G}_{l}$ is well known DLMF  and for $l=0$, $\rho=k^{\prime}r$
and $\eta^{\prime}=e^{2}\mu/(4\pi k^{\prime})$ can be written:

$\displaystyle\mathcal{F}_{0}(\eta^{\prime},\rho)$ $\displaystyle\sim$
$\displaystyle c_{0}\rho\left[1+\eta^{\prime}\rho+\cdots\right]$ (64)
$\displaystyle\mathcal{G}_{0}(\eta^{\prime},\rho)$ $\displaystyle\sim$
$\displaystyle\frac{1}{c_{0}}\Bigl{[}1+2\eta^{\prime}\rho\bigl{\\{}\ln(2\rho)+\mathrm{Re}\left[\psi(i\eta^{\prime})\right]+2\gamma_{E}-1\bigr{\\}}+\cdots\Bigr{]}$
(65)

where

$\displaystyle c_{0}$ $\displaystyle=$
$\displaystyle\left(\frac{2\pi\eta^{\prime}}{e^{2\pi\eta^{\prime}}-1}\right)^{\frac{1}{2}},$
(66)

$\psi(z)=\Gamma^{\prime}(z)/\Gamma(z)$ is the digamma function and
$\gamma_{E}=0.57721\ldots$ is Euler’s constant.

Next we examine the finite-volume Helmholtz equation solution $G_{{\rm
QED}_{L}}(\vec{r})$ and determine its behavior for small $r$. As a first step
we rewrite the differential equation, Eq. (55) obeyed by $G_{{\rm
QED}_{L}}(\vec{r})$ as the Lipman-Schwinger integral equation:

$G_{{\rm
QED}_{L}}(\vec{r})=\int_{V}d^{3}r^{\prime}G_{0}(\vec{r}-\vec{r}\,^{\prime})\left[\delta(\vec{r}\,^{\prime})+2\mu
V_{\mathrm{QED}_{L}}(\vec{r}\,^{\prime})G_{{\rm
QED}_{L}}(\vec{r}\,^{\prime})\right]$ (67)

where $G_{0}(\vec{r})$ is the Helmholtz solution with $e^{2}=0$ given in Eq.
(51). Equation (67) can be expanded through first order in $e^{2}$ giving

$\displaystyle G_{{\rm QED}_{L}}(\vec{r})$ $\displaystyle\approx$
$\displaystyle
G_{0}(\vec{r})+\int_{V}d^{3}r^{\prime}G_{0}(\vec{r}-\vec{r}\,^{\prime})2\mu
V_{\mathrm{QED}_{L}}(\vec{r}\,^{\prime})G_{0}(\vec{r}\,^{\prime}).$ (68)

Now we must determine the singular and regular parts of the right-hand side of
Eq. (68). 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}(\vec{r},E)$ can be approximated by its infinite-volume counter part:

$\displaystyle G^{\infty}_{0}(\vec{r})$ $\displaystyle=$ $\displaystyle\int
d^{3}k\frac{e^{i\vec{q}\cdot\vec{r}}}{(2\pi)^{3}}\frac{1}{-\vec{q}\,^{2}+k^{2}}$
(69) $\displaystyle=$ $\displaystyle-\frac{1}{4\pi}\frac{e^{ikr}}{r}\,.$ (70)

The singular part of the zeroth-order term in Eq. (68) is given immediately by
Eq. (70). However, the singular part of the first-order term can also be found
using this same expression for $G^{\infty}_{0}(\vec{r})$ by evaluating:

$\displaystyle\int_{V}d^{3}r^{\prime}G_{0}(\vec{r}-\vec{r}\,^{\prime})V_{\mathrm{QED}_{L}}(\vec{r}\,^{\prime})G_{0}(\vec{r}\,^{\prime})$
(76)
$\displaystyle\hskip-144.54pt\approx\int_{V}d^{3}r^{\prime}G^{\infty}_{0}(\vec{r}-\vec{r}\,^{\prime})V_{C}(\vec{r}\,^{\prime})G^{\infty}_{0}(\vec{r}\,^{\prime})$
$\displaystyle\hskip-144.54pt=\frac{e^{2}}{(4\pi)^{3}}\int
d^{3}r^{\prime}\frac{e^{ik|\vec{r}-\vec{r}\,^{\prime}|}}{|\vec{r}-\vec{r}\,^{\prime}|}\frac{e^{ikr^{\prime}}}{(\vec{r}\,^{\prime})^{2}}$
$\displaystyle\hskip-144.54pt\approx\frac{e^{2}}{32\pi^{2}}\int_{0}^{\infty}dr^{\prime}e^{ikr^{\prime}}\int_{-1}^{1}\frac{d\cos(\theta)}{\sqrt{r^{2}+{r^{\prime}}^{2}-2rr^{\prime}\cos(\theta)}}$
$\displaystyle\hskip-144.54pt\approx\frac{e^{2}}{(4\pi)^{2}}\int_{0}^{\infty}dr^{\prime}e^{ikr^{\prime}}\left[\frac{1}{r}\theta(r-r^{\prime})+\frac{1}{r^{\prime}}\theta(r^{\prime}-r)\right]$
$\displaystyle\hskip-144.54pt\approx\frac{e^{2}}{(4\pi)^{2}}\left(\ln(1/kr)+\cdots\right)$

where the dots represents terms that are constants or of higher order as $r\to
0$. Since we are interested in the singularity at $r=0$, we need not be
concerned about the large $r^{\prime}$ 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. (68) 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}(\vec{r})$ term in Eq. (68) and found:

$\displaystyle\lim_{r\to 0}G_{0}(\vec{r})$ $\displaystyle\sim$
$\displaystyle-\frac{1}{4\pi
r}-\left.\frac{1}{4\pi^{2}L}\sum_{\vec{n}}\frac{1}{\left[(\vec{n})^{2}-\left(\frac{kL}{2\pi}\right)^{2}\right]^{s}}\right|_{s=1}$
(77) $\displaystyle=$ $\displaystyle-\frac{1}{4\pi
r}-\sqrt{4\pi}\frac{1}{4\pi^{2}L}Z_{00}\left(1,q^{2}\right),$ (78)

where $q=kL/2\pi$ and the evaluation at $s=1$ in Eq. (53) is described in Ref.
Luscher (1991). We must now make a similar evaluation of the right most term
in Eq. (68).

In analogy with Eq. (78), we define the regular part of the $e^{2}$ correction
to the function $G_{{\rm QED}_{L}}(\vec{r})$ as

$\displaystyle Y(q)$ $\displaystyle=$ $\displaystyle\lim_{r\to
0}\sum_{\vec{n}}\sum_{\vec{m}\neq
0}\frac{e^{i\vec{n}\cdot\frac{2\pi\vec{r}}{L}}}{(2\pi)^{4}\pi}\frac{1}{(\vec{n})^{2}-q^{2}}\frac{1}{\vec{m}\,^{2}}\frac{1}{(\vec{n}-\vec{m})^{2}-q^{2}}-\frac{1}{4\pi}\ln(1/kr).$
(79)

The relatively mild logarithmic singularity allows a direct numerical
evaluation $Y(q)$.

Next, following Lüscher we introduce a Coulomb function $\phi_{C}(E)$ such
that the combination:

$\cos\left(\phi_{C}(E^{\prime})\right)\mathcal{F}_{0}(k^{\prime}r)-\sin\left(\phi_{C}(E^{\prime})\right)\mathcal{G}_{0}(k^{\prime}r)$
(80)

contains the same ratio of regular to singular parts as we have found in
$G_{{\rm QED}_{L}}(E)$. Equating the ratio of regular divided by singular
parts in the limit of small $r$ in Eqs. (80) and (68) we obtain:

$\displaystyle\frac{(\frac{2\pi\eta}{e^{2\pi\eta}-1})^{\frac{1}{2}}\,k^{\prime}r\Bigl{[}\cos(\phi_{C})-\frac{e^{2\pi\eta}-1}{2\pi\eta}\sin(\phi_{C})2\eta^{\prime}\left\\{\mathrm{Re}\left[\psi(i\eta^{\prime})\Bigr{]}+2\gamma_{E}-1\right\\}\right]}{-\sin(\phi_{C})\left[(\frac{e^{2\pi\eta}-1}{2\pi\eta})^{\frac{1}{2}}\Bigl{(}1+2\eta^{\prime}k^{\prime}r\ln(k^{\prime}r)\Bigr{)}\right]}$
$\displaystyle\hskip-144.54pt=\frac{-\frac{1}{4\pi^{2}L}\sqrt{4\pi}Z_{00}(1,\frac{kL}{2\pi})+\frac{2\mu
e^{2}}{4\pi}Y\left(\frac{kL}{2\pi}\right)}{-\frac{1}{4\pi r}+\frac{2\mu
e^{2}}{(4\pi)^{2}}\ln(1/kr)}$

where to determine the left-hand side we have used Eqs. (64) and (65) while
the right-hand side has been determined from Eqs. (78) and (79). In contrast
with the zeroth order in $e^{2}$ case, the singular solution $\mathcal{G}_{0}$
at order $e^{2}$ contains a regular term which must be included. To allow a
direct comparison of Eq. (A.2) and Eqs. (64) and (65) we have not expanded
this equation to first order in $e^{2}$.

Recognizing that $\eta$ with $\mu e^{2}/k$ we can simplify in Eq. (A.2) by
expanding to first order in $\eta$ to obtain a final formula for $\phi_{C}(E)$
to order $e^{2}$:

$\displaystyle\cot\bigl{(}\phi_{C}(E^{\prime})\bigr{)}$ $\displaystyle=$
$\displaystyle-(1+\pi\eta)\frac{1}{\pi^{\frac{3}{2}}}\left(\frac{2\pi}{k^{\prime}L}\right)Z_{00}\left(1,\frac{kL}{2\pi}\right)$
(82) $\displaystyle\hskip 36.135pt+8\pi\eta
Y\left(\frac{kL}{2\pi}\right)+2\eta\Bigl{\\{}\mathrm{Re}\left[\psi(i\eta)\right]+2\gamma_{E}-1\Bigr{\\}}$

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 non-relativisitic $s$-wave scattering phase
shift $\bar{\delta}_{0}(E^{\prime})$ given by the quantization condition:

$\cot(\bar{\delta}_{0}(E^{\prime}))+\cot(\phi_{C}(E^{\prime}))=0\,,$ (83)

where $E^{\prime}=E+\Delta E$ where $\Delta E$ is given in Eq. (63). Equation
(83) 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 ${\rm QED}_{L}$ finite-volume quantization

After identifying the $O(\frac{r^{2}}{L^{3}})$ correction to the ${\rm
QED}_{L}$ potential, one could become concerned about the size of such
corrections. In order to understand the efficacy of the ${\rm
QED}_{L}$-modified 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 ${\rm 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.

###  B.1 $\delta$-shell potential

Quantum mechanics with a $\delta$-shell potential,

$V({\vec{r}})=V_{0}\,\delta(|{\vec{r}}|-r_{0})\,,$ (84)

is a sufficiently simple system to solve both analytically in infinite volume
and numerically in a finite volume. The $\delta$-shell potential has the
feature necessary for application of Lüscher’s quantization condition, as well
as the ${\rm QED}_{L}$-modified finite-volume quantization condition, that
there exists an exterior region of $|{\vec{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
$\delta$-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

$\psi(r>r_{0})=kAj_{l}(kr)-kBn_{l}(kr)\,,$ (85)

where phase shift is found through the ratio $B/A=\tan\delta_{l}$. Since the
wave function must remain finite, in the interior region the wave function
will contain only the regular solution

$\psi(r<r_{0})=kj_{l}(kr)\,,$ (86)

up to an arbitrary normalization.

For a $\delta$-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

$\tan\delta_{l}=\frac{\frac{V_{0}}{k}j_{l}^{2}(kr_{0})}{(kr_{0})^{-2}-\frac{V_{0}}{k}j_{l}(kr_{0})n_{l}(kr_{0})}\,,$
(87)

where the first term in the denominator comes from the Wronskian of these
functions $n_{l}^{\prime}(z)j_{l}(z)-n_{l}(z)j_{l}(z)=z^{-2}$. For the
$s$-wave phase shift, this equation reduces to the textbook result:

$\tan\delta_{0}=\frac{\frac{V_{0}}{k}\sin^{2}(kr_{0})}{1+\frac{V_{0}}{k}\sin(kr_{0})\cos(kr_{0})}\,.$
(88)

###  B.2 Combined $\delta$-shell and Coulomb potential

An identical derivation can be performed in the presence of the Coulomb
potential, where the spherical Bessel functions are replaced by the regular,
$F_{l}(\eta,z)$, and irregular, $G_{l}(\eta,z)$, Coulomb wavefunctions
replacing $j_{l}(z)$ and $-n_{l}(z)$ respectively. These wavefunctions are
related to the confluent hypergeometric function $U$ through their relation to
the Coulomb analogue to the Hankel functions

$\displaystyle H_{l}^{\pm}(\eta,kr)$ $\displaystyle=$ $\displaystyle
G_{l}(\eta,kr)\pm iF_{k}(\eta,kr)$ (89) $\displaystyle=$ $\displaystyle e^{\pm
i\Theta(kr)}(\mp 2ikr)^{1+l\pm i\eta}U(1+l\pm i\eta,2l+2,\mp 2ikr)\,,$

where $\Theta(z)=z-l\pi/2+\sigma_{l}(\eta)-\eta\ln(2z)$ and
$\sigma_{l}(\eta)={\rm arg}\,\Gamma(1+l+i\eta)$ is the Coulomb phase shift. .

The exterior and interior wavefunction are parameterized in the just as in
(85) and (86) and are subject to the same boundary conditions at the radius of
the $\delta$-shell. Using the Wronskian of the Coulomb wavefunctions
$G_{l}(z)F_{l}^{\prime}(z)-G_{l}^{\prime}(z)F_{l}(z)=1$, the phase shift in
the presence of the Coulomb interaction will be given by

$\tan\delta_{l}=\frac{\frac{V_{0}}{k}F_{l}^{2}(\eta,kr_{0})}{1+\frac{V_{0}}{k}F_{l}(\eta,kr_{0})G_{l}(\eta,kr_{0})}\,.$
(90)

By comparing numerical results with this analytical result, the efficacy of
the ${\rm QED}_{L}$-modified finite-volume quantization condition can be
tested.

###  B.3 Computational strategy

In the remainder of this appendix we will calculate numerically the finite-
volume energies for the combined $\delta$-shell and ${\rm QED}_{L}$
potentials, apply the ${\rm QED}_{L}$ finite-volume quantization condition
given in Eq. (83) and then compare the results for the scattering phase shift
with the exact infinite-volume values given by Eq. (90). Our numerical
approach is to determine the finite-volume energies and wave functions to
order $e^{0}$ and then to use first-order perturbation to determine the
$O(e^{2})$ correction.

To solve the Schrödinger equation with $\delta$-shell potential alone in
finite volume, again the wave function is parameterized in the regions
interior and exterior to the $\delta$-shell. The wave function in the interior
region is the same as given in Eq. (86), 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. (51). The boundary conditions
at the $\delta$-shell leads to the quantization condition

$\left[1-V_{0}r_{o}-kr_{0}\cot(kr_{0})\right]G_{0,0}(r_{0};k)+r_{0}\frac{\partial
G_{0,0}}{\partial r}(r_{0};k)=0\,.$ (91)

where $G_{0,0}(r)$ is the $s$-wave component of the zeroth-order Green’s
function $G_{0}(\vec{r})$ defined in Eq. (51):

$G_{0,0}(r)=\frac{1}{4\pi}\int d\Omega_{\hat{r}}G_{0}(r\hat{r}).$ (92)

We will impose no conditions on higher partial waves with $l\geq 4$ making the
usual assumption that the $\delta$-shell phase shifts for $l\geq 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.
(91) 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. (51). From those
energy levels and the corresponding wavefunction, the ${\rm QED}_{L}$
potential will be included using first-order perturbation theory by
calculating $\langle\psi|V_{{\rm QED}_{\rm L}}|\psi\rangle$ numerically, as
described in Sec. B.5.

###  B.4 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. (51), as
a sum in momentum space the series is slowly convergent. If one truncates the
sum to only contain terms with ${\vec{p}}<p_{\rm max}$, then the leading
neglected terms will be proportional to $p_{\rm max}^{-2}$, but there will be
approximately $p_{\rm max}^{2}$ 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_{\rm 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}({\vec{r}};k)=\frac{1}{4\pi}\sum_{\vec{n}}\frac{e^{ik|\vec{r}-\vec{n}L|}}{|\vec{r}-\vec{n}L|}$
(93)

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” Ham and
Segall (1961). Ewald summation breaks the series into two components which
represent the “near” contributions, which come from contributions of periodic
images close to $\vec{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 Ham
and Segall (1961). 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. (93) is rewritten as an
integral

$\frac{e^{ikr}}{4\pi
r}=\frac{1}{2\pi^{3/2}}\int_{\Gamma}e^{-r^{2}\xi^{2}}e^{k^{2}/\xi^{2}}d\xi\,,$
(94)

where $\Gamma$ is a particular contour from 0 to $\infty$ with some
constraints on how it approaches those limits Ham and Segall (1961). First
$\Gamma$ must leave the origin with a negative angle $\phi$ in the limits
$0>\phi>-\frac{\pi}{4}$, second it returns to the real axis at $\xi_{0}$, and
then finally, continues on the real axis to $\infty$. This integral is then
broken into two regions, from 0 to $\xi_{0}$ and from $\xi_{0}$ to $\infty$.
The first of these regions will be evaluated in momentum space and the second
in position space

$\displaystyle
G_{0}({\vec{r}};k)=G_{1}({\vec{r}};k)+G_{2}({\vec{r}};k)=\frac{1}{L^{3}}\sum_{p}\frac{e^{i{\vec{p}}\cdot{\vec{r}}}e^{\frac{k^{2}-p^{2}}{4\xi_{0}}}}{k^{2}-p^{2}}+$
$\displaystyle\frac{-1}{4\pi}\sum_{n}{\rm
Re}\left[e^{i|{\vec{r}}+{\vec{n}}L|k}{\rm
erfc}\left(|{\vec{r}}+{\vec{n}}L|\xi_{0}+\frac{ik}{2\xi_{0}}\right)\right]\,.$
(95)

These sums converge significantly faster than the original sum. The first
converges as $e^{-p^{2}/4\xi_{0}^{2}}$ and the second as
$e^{-n^{2}L^{2}\xi_{0}^{2}}$. The original momentum space sum is the special
case of $\xi_{0}=0$ and the position space sum is reached as
$\xi_{0}\to\infty$ which both demonstrate the slow convergence. Using the
optimal value of $\xi_{0}=3^{\frac{1}{4}}\sqrt{\frac{\pi}{L}}$ Kustepeli and
Martin (2000), these sums can converge to sufficient accuracy with just a few
terms.

These summation methods will also be necessary for numerically evaluating the
${\rm QED}_{L}$ potential efficiently. Unlike in a lattice application, with a
UV regulator, the standard summation definition ${\rm QED}_{L}$ potential is
also slowly convergent. For the numerical work, the ${\rm QED}_{L}$ potential
will be evaluated using Eq. (50).

###  B.5 Numerical analysis of $\delta$-shell potential

The numerical analysis of the $\delta$-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 $\delta$-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. 91 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\delta$, 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.

Figure 9:  The scattering phase shift $\delta$ (Left) and $k\cot\delta$
(Right) for the $\delta$-shell potential. The strength of the $\delta$-shell
potential is 1 for the upper plots and 10 for the lower plots. The black curve
represents the exact infinite-volume solution of Eq. 88.

Perturbation theory is used to find the energy shifts caused by the inclusion
of the ${\rm QED}_{L}$ potential. The unnormalized wavefunction corresponding
to a finite-volume energy eigenstate for the $\delta$-shell potential alone is
approximately

$\psi({\vec{r}};k)=\begin{cases}\frac{\sin(kr)}{r}&r<r_{0}\\\
\frac{\sin(kr_{0})}{r_{0}G_{0,0}(r_{0};k)}G_{0}({\vec{r}};k)&r\geq
r_{0}\end{cases}$ (96)

up to order $1/L^{2l}$ power-law corrections generated by ignoring
contributions with $l\geq 4$ in the interior region. Adding $V_{{\rm QED}_{\rm
L}}$ as a perturbation, the energy shift is calculated as $\Delta
E=\frac{\langle\psi|V_{{\rm QED}_{\rm
L}}|\psi\rangle}{\langle\psi|\psi\rangle}$ by analytically integrating the
component of the wavefunction for $r<r_{0}$ and then performing a Monte-Carlo
integration for the exterior region. The value of the coupling is set to
$e^{2}=\frac{1}{2\mu}$, such that the effective Bohr radius, $4\pi/(e^{2}\mu)$
is significantly larger than the size of the box. In the region of small $k$,
the Sommerfeld parameter $\eta$ 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. (83). 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 sub-percent 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 $\eta$.

Figure 10:  The scattering phase shift $\delta$ (Left) and $k\cot\delta$
(Right) for the combined $\delta$-shell and Coulomb potential. The strength of
the $\delta$-shell potential is 1 for the upper plots and 10 for the lower
plots. The continuous black curve represents the exact infinite-volume
solution of Eq. 90 while the points are calculated using the finite-volume
${\rm QED}_{L}$ quantization condition developed here. The dash black curve
represents the exact infinite-volume solution without the Coulomb interaction
given by Eq. 88.

This study of a simple quantum mechanics problem implies that the structure-
dependent effects of the unphysical $L^{-3}$ introduced by the ${\rm 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 be the size of the error in the phase shift. The relative error
in the resulting phase shift for this case is $O(0.1\%)$. In fact, it appears
the larger deviations from the infinite-volume non-perturbative result are
from the truncation of perturbation theory in the finite-volume quantization
condition of Eq. (83), which grow largest at small $k$ where $\eta$ goes above
1. Of course, there is no guarantee the the structure-dependent effects of the
$1/L^{3}$ term will be the same in an actual QCD calculation, but analyzing
this would require a full numerical lattice study comparing this ${\rm
QED}_{L}$ approach to a more accurate method such as the truncated Coulomb
potential approach presented in this paper.

## References

  * Batley  _et al._ (2002) J. R. Batley _et al._ (NA48), Phys. Lett. B544, 97 (2002), arXiv:hep-ex/0208009 [hep-ex] .