Charged multi-hadron systems in lattice QCD+QED
S. R. Beane W. Detmold R. Horsley M. Illa M. Jafry D. J. Murphy Y. Nakamura H. Perlt P. E. L. Rakow G. Schierholz P. E. Shanahan H. Stüben M. L. Wagman F. Winter R. D. Young J. M. Zanotti
hep-lat hep-ph nucl-th


DESY 20-28
FERMILAB-PUB-20-123-T
ICCUB-20-007
MIT-CTP/5183
NT@UW-20-03


Department of Physics, University of Washington,Box 351560, Seattle, WA 98195, U.S.A.
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A.
 School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD, UK
 Departament de Física Quàntica i Astrofísica, Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona, Martí Franquès 1, E08028 Barcelona, Spain
Department of Physics, University of Washington,Box 351560, Seattle, WA 98195, U.S.A.
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A.
 RIKEN Center for Computational Science, Kobe, Hyogo 650-0047, Japan
 Institut für Theoretische Physik, Universität Leipzig, 04109 Leipzig, Germany
 Theoretical Physics Division, Department of Mathematical Sciences,University of Liverpool, Liverpool L69 3BX, UK
 Deutsches Elektronen-Synchrotron DESY, 22603 Hamburg, Germany
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A.
 Regionales Rechenzentrum, Universität Hamburg, 20146 Hamburg, Germany
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A.
Fermi National Accelerator Laboratory, Batavia, IL 60510, U.S.A.
 Jefferson Laboratory, 12000 Jefferson Avenue, Newport News, VA 23606, USA
CSSM, Department of Physics, University of Adelaide, Adelaide SA 5005, Australia
CSSM, Department of Physics, University of Adelaide, Adelaide SA 5005, Australia

NPLQCD and QCDSF collaborations

  Systems with the quantum numbers of up to twelve charged and neutral pseudoscalar mesons, as well as one-, two-, and three-nucleon systems, are studied using dynamical lattice quantum chromodynamics and quantum electrodynamics (QCD+QED) calculations and effective field theory.
QED effects on hadronic interactions are determined by comparing systems of charged and neutral hadrons after tuning the quark masses to remove strong isospin breaking effects.
A non-relativistic effective field theory, which perturbatively includes finite-volume Coulomb effects, is analyzed for systems of multiple charged hadrons and found to accurately reproduce the lattice QCD+QED results.
QED effects on charged multi-hadron systems beyond Coulomb photon exchange are determined by comparing the two- and three-body interaction parameters extracted from the lattice QCD+QED results for charged and neutral multi-hadron systems.







Charged multi-hadron systems in lattice QCD+QED
    J. M. Zanotti
    Received 29 November 2019 / Accepted 07 March 2020
======================================================






§ INTRODUCTION



The interplay between the strong and electromagnetic (EM) interactions of the Standard Model is subtle but central to the complexity of visible matter. While EM interactions are typically much weaker than the strong interactions, it is their competition with strong isospin-breaking effects that leads to the observed proton–neutron mass difference and determines the stability of atomic nuclei.
Moreover, the hierarchy between the length scales where strong and EM interactions are important leads to the existence of chemistry and all of the complexity that it entails.
Progress towards understanding the combined effects of the strong and EM interactions from first principles has been made in recent years, and the EM–strong interaction decomposition of the neutron-proton mass difference, M_n - M_p, has been calculated from the Standard Model for the first time using  the numerical  lattice formulation of quantum chromodynamics  (QCD) and quantum electrodynamics (QED)   <cit.>.
Coupled lattice QCD+QED calculations have also been performed for meson and baryon masses <cit.>, leptonic decay rates <cit.>, and the anomalous magnetic moment of the muon <cit.>, as both QCD and QED must be accounted for in precision tests of the Standard Model. 

Systems of multiple charged hadrons exhibit rich phenomenology in which the interplay between QCD and QED effects is less well understood than for single-hadron systems.
The differences between pp, np and nn two-nucleon scattering channels are much more poorly constrained than their isospin average.
Lattice QCD+QED can potentially provide phenomenologically useful information on such isospin-breaking effects in nucleon-nucleon scattering;
furthermore, calculations with a range of quark masses and values of the QED fine structure constant can provide insight into the fine-tuning of QED effects and strong isospin breaking interactions, providing complementary information to experiment.
In large nuclei with charge Z∼1/α, where α is the EM fine structure constant, relativistic QED effects are expected to become nonperturbative and give rise to electron-positron pair-production and interesting consequences for the structure of elements <cit.>.


QCD+QED studies of systems involving two or more protons are of great phenomenological interest, but calculations of large nuclei are difficult using current lattice QCD+QED techniques <cit.>.
It is also interesting to consider systems of multiple charged mesons, which allow an investigation of similar QED effects without a number of the numerical complexities of multi-proton calculations.
In particular, charged multi-hadron systems include QED effects that become non-perturbative for hadron pairs with sufficiently small relative velocity v ≪α.
All calculations of charged hadrons must also reconcile the presence of net charge in a finite volume (FV) with Gauss's law, which can be accomplished in several ways <cit.>.
The current study makes use of the QED_L formulation <cit.>, in which the photon spatial zero-mode is removed, as has been recently studied in detail for single-hadron systems in Refs. <cit.>.
Subtleties related to the nonlocality of zero-mode subtraction in QED_L have been understood for single-hadron systems <cit.>, and similar methods are used to understand nonlocality in multi-hadron systems in this framework as described below. 


In lattice QCD calculations, hadronic scattering phase shifts are determined by relating them to the FV energies of two-hadron systems using FV relations first derived by Lüscher <cit.> and generalized in recent years (see Ref. <cit.> for a review).
Understanding scattering in lattice QCD+QED is complicated by the lack of a FV formalism that includes nonperturbative Coulomb effects, and thereby relates FV energy shifts calculable in lattice QCD+QED to infinite-volume (Coulomb-corrected) scattering phase shifts.
In the QED_L formulation, Ref. <cit.> argues that Coulomb effects can be treated perturbatively for sufficiently small volumes where L ≪ 1/(α M) and therefore the Bohr radius of the charged-particle system is larger than the infrared cutoff provided by the FV. 
However, the volume must also satisfy L ≫ 1/M in order for relativistic FV effects not accounted for in Lüscher's analysis to be negligible.
For some systems there may be a window of intermediate-sized volumes where relativistic FV effects and nonperturbative Coulomb effects can both be neglected.
Non-relativistic effective field theory (EFT) can be used to investigate this issue and to relate FV energy shifts for systems with more than two particles to two-, three-, and higher-body EFT low-energy constants (LECs).

This work studies QCD+QED effects for systems of multiple hadrons including up to twelve charged or neutral mesons and up to three nucleons in multiple finite volumes. 
A larger than physical value of the fine-structure constant, α≃ 0.1, is used in order to increase the size of QED effects and permit the study of systems with total charge Z satisfying Z α≳ 1.
As described in Ref. <cit.> and summarized below, the quark masses are tuned to remove strong isospin breaking effects so that splittings between particles that are degenerate in QCD can be identified as pure QED effects.
Calculations are performed in both lattice QCD+QED_L (LQCD+QED_L) and in an EFT appropriate for non-relativistic charged hadrons, referred to as non-relativistic QED_L (NRQED_L) below.
Two-body scattering lengths and three-body interactions of charged and neutral mesons are determined by tuning the hadronic interaction parameters appearing in the NRQED_L Lagrangian to reproduce the LQCD+QED_L FV energy levels.
QED effects on hadron interactions beyond Coulomb photon exchange are determined by comparing the two-body and three-body interaction parameters determined for charged and neutral mesons.
Studies of multi-baryon systems are used to probe volumes where L > 1/(α M) and Coulomb effects may not be perturbative according to the criteria of Ref. <cit.>. 
It is noteworthy that no significant nonperturbative Coulomb effects or relativistic effects are seen at the level of precision of the results obtained here, even for systems with  L > 1/(α M) or Z α > 1.

The practical window of L where relativistic effects can be neglected and Coulomb effects can be treated perturbatively is explored by comparing LQCD+QED_L and NRQED_L results for a variety of FV systems.
From these studies, it is argued that this window requires the spatial extent L of a cubic FV to satisfy

    1/ML≪ 1, α M L/4π≪ 1.

The combination α M L / (4π) in Eq. (<ref>) that quantifies the size of FV Coulomb effects is equivalent to the infinite-volume Coulomb expansion parameter α / v for a pair of hadrons moving back-to-back with one unit of FV momentum, i.e. ± 2π / L.
Eq. (<ref>) indicates that Coulomb effects can be treated perturbatively in current and future LQCD+QED_L calculations over a wide range of volumes and in particular for systems with L ∼ 1/ (α M), where the Bohr radius is commensurate with the volume.
This scaling differs from that suggested in Ref. <cit.> by the factor of 4π in the denominator of the α M L / (4π).

This manuscript is structured as follows. Section <ref> presents the details and methodology of the LQCD+QED_L calculations that are performed.  Section <ref> discusses the construction of NRQED_L for charged multi-hadron systems and provides the formalism needed to extract hadronic scattering lengths and other parameters from the results of these LQCD+QED_L calculations. Section <ref> describes the determination of hadronic scattering parameters from LQCD+QED_L and presents results for multi-meson and multi-nucleon systems. Conclusions are presented in Section <ref>. Appendix <ref> contains further technical details on the matching between QED_L and NRQED_L, and Appendix <ref> contains details on the fitting procedure used to extract energy levels from Euclidean correlation functions computed in LQCD+QED_L.




§ LATTICE QCD + QED



Lattice QCD+QED is a nonperturbative approach to QCD+QED that at intermediate stages implements an ultraviolet regulator defined by a lattice spacing a (where 1/a is assumed to be much smaller than the QED Landau pole).
Calculations are performed in Euclidean spacetime with a cubic spatial volume of extent L× L × L and a finite temporal extent T;
the quark, gluon, and photon fields satisfy periodic boundary conditions (PBCs) in all spatial directions.
Here, the QED_L formalism <cit.> is used to define charged-particle correlation functions as detailed below, which defines the LQCD+QED_L formalism.
The LQCD+QED_L gauge field configurations used were generated by the QCDSF collaboration. 
The full details of the generation of these ensembles are presented in Refs. <cit.>. 
For completeness, relevant aspects of the ensemble generation are described below.




 §.§ Lattice action and parameters




 
The LQCD+QED_L action used in this study is given by

    S=S_G+S_A+S_F,

where S_G is the tree-level Symanzik-improved SU(3) gauge action as defined in Refs. <cit.>.
The quark dynamics are encoded by an O(a)-improved stout link
non-perturbative clover (SLiNC) action:

    S_F =   ∑_q ∈{u,d,s}∑_x {1/2∑_μ[
    q̅(x)(γ_μ-1)e^-ie_q A_μ(x)Ũ_μ(x)q(x+μ̂) . . 
       . .- q̅(x)(γ_μ+1) e^ie_qA_μ(x-μ̂)Ũ_μ^†(x- μ̂)q(x- μ̂)]
    . 
       .  + 1/2 κ_qq̅(x)q(x)-1/4c_SW∑_μνq̅(x)σ_μν F_μν q(x) },

where Ũ_μ is a single-iterated stout-smeared
SU(3) link <cit.>, A_μ is a non-compact U(1) gauge field, and the U(1) quark charges are given by
e_q.
The field-strength F_μν appearing in the Sheikholeslami-Wohlert or “clover term” <cit.> involves the unsmeared SU(3) gauge-field as in Ref. <cit.>.
The clover coefficient c_SW was non-perturbatively determined for pure QCD in Ref. <cit.>.
Electromagnetic gauge fields are not included in the clover term,
however with these simulation parameters, the O(α a) effects are no larger than the
residual O(a^2) effects of pure QCD.
The photon action is described by the non-compact form:

    S_A = 1/2e^2∑_x,μ<ν(A_μ(x)+A_ν(x+μ)-A_μ(x+ν)-A_ν(x) )^2

where e is the U(1) gauge coupling corresponding to α = e^2 / (4π). Gauge fixing and the treatment of zero modes are discussed in Sec. <ref>.

The parameters of the lattice action
were tuned by identifying a point of approximate SU(3) flavor symmetry, where the average
light-quark mass takes its physical value—see
Ref. <cit.> for further discussion.
With dynamical QED, this is complicated by the fact that the quark
charges explicitly break the SU(3) flavor symmetry.
An approximate SU(3) flavor symmetry is then realized by tuning the quark mass
parameters such that the connected flavour-neutral pseudoscalar mesons[Since κ_s = κ_d and e_s = e_d, this theory exhibits exact
  U-spin symmetry, and therefore the connected
  dd̅ and ss̅ correlation functions are identical to the ds̅ correlation functions
The connected part of the uu̅ correction function can be interpreted in a partially quenched theory, see Ref. <cit.> for further discussion.]
are 
degenerate.
As it is inspired by Dashen's theorem <cit.>,
this prescription for separating strong and electromagnetic effects is known as
the “Dashen scheme” <cit.>.
This scheme preserves Dashen's theorem, whereby the neutral pseudoscalar mesons 
are
protected from receiving an electromagnetic self-energy correction in the chiral limit.
To reach the physical quark masses and charges, the SU(3) flavor symmetric point can then be lowered to the physical value before the symmetry is broken to fix the individual quark masses.
In this exploratory work a single set of approximately SU(3) flavor symmetric parameters is used.


The Dashen scheme provides a natural framework for separating QED
and QCD effects; at the SU(3) symmetric point, any splittings among pure-QCD multiplets are
identified as pure QED effects.
The explicit breaking of the quark masses from this point then allows
one to independently isolate the effects of strong (or quark mass)
flavor symmetry breaking and QED effects.
The action parameters used in this work correspond to the SU(3) symmetric point and are presented in Table <ref>. 




 §.§ Gauge fixing and the U(1) zero mode


The correlation functions of charged particles are not gauge invariant and hence ensemble-averaged quantities are only meaningful within some gauge fixing prescription.
In this work, Landau gauge is adopted by enforcing the condition

    ∑_μ(A_μ(x)-A_μ(x-μ̂))=0.

This condition can be imposed after generation of the gauge fields.
However, the Landau gauge condition leaves the 4-dimensional zero mode unconstrained.
These uniform background fields do not contribute to the gauge action, however they do couple to the fermionic action, Eq. (<ref>).
The action remains invariant under discrete shifts of the zero mode in units of 2π/L_μ.
This redundant gauge degree of freedom can be eliminated by mapping the 4-dimensional zero modes onto the finite interval -π/L_μ<A_μ(k=0)≤π/L_μ <cit.>, where A_μ(k) is the Fourier transform of A_μ(x) defined explicitly in Eq. (<ref>).
The leading effect of these zero modes is to induce a charge-dependent twist of the single-hadron momenta.
This small energy shift can be corrected in single-hadron states <cit.>, however this would introduce an unnecessary complication in the analysis of many-body interactions.
Instead, this work adopts the QED_L prescription <cit.> of setting the spatial (3-dimensional) zero modes of the gauge potential to zero on every timeslice.

With respect to the action used to generate the gauge configurations, the elimination of the 3-dimensional zero modes before computing quark propagators is not a gauge symmetry.
As a consequence, there is a partial-quenching effect whereby the valence quarks experience a different zero mode to the quarks in the sea.
The zero modes cannot affect closed fermion loops, and their only contribution to the determinant will be associated with fermion lines wrapping around the boundary of the lattice.
Consequently, this partial-quenching effect is exponentially suppressed in m_π L and is negligible in comparison to the power-law FV effects studied in this work.









 §.§ Correlation functions


The particular  gauge field ensembles used in this work are described in Table <ref>; along with the parameters used in their generation, the number of configurations and the average number of correlation function source locations that are used per configuration are also reported. 
Up and down/strange (equivalent for the masses used here) quark propagators are computed from each of the randomly chosen source locations  using three-dimensionally Jacobi smeared sources <cit.> (100 iterations with ρ=0.21) using a solver tolerance of 10^-12. Local and smeared fields are used in the sink interpolating operator to construct smeared-point (SP) and smeared-smeared (SS) correlation functions with the former providing cleaner signals in all cases.



FV energy levels are determined by analyzing two-point correlation functions

    G_h(t;  x_0) =⟨∑_ xχ̃_h( x,t) χ_h^†( x_0,0)⟩ ,

where χ_h (χ̃_h) is a source (sink) interpolating operator with the quantum numbers of the state being considered and x_0 is the source location. Two-point correlation functions are constructed for n∈{1,…,12} charged pions (ud̅) and neutral kaons (sd̅) 
using the techniques developed in Refs. <cit.>. At the sink, each color-singlet meson bilinear is separately projected to zero momentum; for example, for a system of n pions (h = nπ^+),

    G_nπ^+(t;  x_0) =⟨(∑_ xu̅( x,t) γ_5 d( x,t) )^n (u( x_0,0)γ_5 d̅( x_0,0))^n ⟩ .



The correlation functions for the single proton and neutron make use of standard local  interpolating operators χ_p;α =ϵ^ijk (u^i C γ_5 d^j) u^k_α  and χ_n;α =ϵ^ijk (d^i C γ_5 u^j) d^k_α where the parentheses indicate  contraction of spin indices. For two-baryon and three-baryon systems, the contraction techniques of Refs. <cit.> are used, again with each baryon separately projected to zero momentum at the sink. 



 §.§ Finite-volume energy level determinations


Finite-volume energy levels are extracted from the correlation functions for each system using correlated fits to their time dependence. For the multi-meson systems, which factorize easily into color-singlet components, thermal contributions where one component propagates forward in time and another propagates backwards in time are particularly relevant. The corresponding functional forms that are used to fit these correlation functions are 

    f_n M(t,Ě,Ž) = ∑_m=0^n Z_n M; m e^-E_m M t e^-E_(n-m)M(T-t)  + ∑_e=1^N_ states-1 Z_n M^(e) e^-E_n M^(e) t ,

where M ∈{π^+, K^0} labels the type of meson, N_ states is the total number of non-thermal states included in the fit, E_n M is the ground-state energy of the system with the quantum numbers of n mesons of type M, Z_n M; m is the overlap factor describing the amplitude of thermal contributions with m forwards propagating mesons and n-m backwards propagating mesons, and non-thermal excited states are also included with energies E_n M^(e) and overlap factors Z_n M^(e).
In practice, fits with N_ states∈{ 1,2,3 } are used in this work.
The vectors Ě and Ž indicate the energy and overlap factor parameters to be constrained in the fit. 
In order to determine the many parameters of these fitting functions, fits are performed sequentially for increasing n, with the energies E_m M for m<n used as input for f_n M as in Ref. <cit.>. 
Thermal effects arising from excited states are not found to be significant.



For baryon systems, statistical noise grows rapidly with the temporal separation between the source and the sink, and the contributions of thermal states are negligible with respect to the statistical uncertainties. For these systems a simpler fit function is used:

    f_b(t,Ě,Ž)= Z_b e^-E_b t + ∑_e=1^N_ states - 1 Z_b^(e)e^-E_b^(e) t ,

where b labels the type of baryon system and, as for meson systems, the second term corresponds to  a sum over the excited states included in the fitting model.
This work investigates single-nucleon systems with b ∈{ p, n } as well as two-nucleon systems with b ∈{ nn, np(^1S_0), np(^3S_1), pp } and three-nucleon systems with b ∈{^3H, ^3He}.

Best-fit parameters are determined from minimization of the correlated χ^2 function

    χ^2 (Ě,Ž) = ∑_t,t' = t_ min^t_ max (G(t) - f(t,Ě,Ž))(𝒮(λ^*)^-1)_t,t'(G(t') - f(t',Ě,Ž))

for the appropriate correlation functions, G(t), and fit function, f ∈{ f_nM, f_b }, [t_ min, t_ max] is the range of times included in the fit, and  𝒮(λ^*) is an estimate of the covariance matrix described below.
Finite sample-size fluctuations may make the sample covariance matrix ill-conditioned, and shrinkage techniques <cit.> are used to obtain a numerical stable inverse covariance matrix.
Following the application of shrinkage to LQCD in Ref. <cit.>, the covariance matrix including shrinkage is defined as 𝒮(λ) = (1-λ)𝒞 + λ𝒯, where 𝒞 is the bootstrap covariance matrix and 𝒯 = diagonal(𝒞), and therefore χ^2-minimization with 𝒮(λ) interpolates between correlated χ^2-minimization for λ = 0 and uncorrelated χ^2-minimization for λ = 1.
The optimal shrinkage parameter λ^* appearing in Eq. (<ref>) is chosen to maximize the expected closeness to the true covariance matrix and defined in Eq. (<ref>), see Appendix <ref> and Ref. <cit.> for further discussion.


Systematic uncertainties arise from the dependence of the fits on the functional forms and time ranges that are used. To address these, the time ranges are systematically sampled, and fits with and without excited states are attempted, with an information criterion used to select the appropriate number of excited states to include for each fit range choice.
A weighted average of the results from all successful fits passing various reliability checks is used to determine the final results.
Further details are presented in Appendix <ref>.


















For one- and two-nucleon systems, combined fits to the SS and SP correlation functions are performed using common energies but different overlap factors for local and smeared sources.
For mesons and three-nucleon systems, combined fits performed in this way are only marginally more precise than fits using SP correlation functions alone, and fits using only SP correlation functions are therefore used in what follows for simplicity.


An effective energy plot that removes thermal effects from backwards propagating states for n=1 systems, as well as constant contributions from thermal effects on n ≥ 2 correlation functions, is defined as

    E(t) = ArcCosh( G(t) - G(t+4)/2 [G(t+1) - G(t+3) ] ).

Note that this effective ground-state energy function differs from the form ArcCosh[ (G(t+1) + G(t-1))/(2G(t)) ] commonly used to remove thermal effects from single-pion correlation functions. The advantage of Eq. (<ref>) is that constant terms present in the spectral representation Eq. (<ref>) for n ≥ 2 mesons exactly cancel in the correlation function differences appearing in Eq. (<ref>), and in particular Eq. (<ref>) is able to exactly isolate the ground-state energy from a two-meson correlation function including thermal effects of the form e^-E t + e^-E(T-t) + const.
Additional terms with more complicated t dependence arise for n ≥ 3 meson correlation functions as shown in Eq. (<ref>). Effective energies for n ≥ 3 systems will therefore include contamination from thermal effects and will only plateau to the true ground-state energy in the region t ≪ T where thermal effects can be neglected.

Figs. <ref>-<ref> show effective energy plots, including correlation function results and fit results for E(t) for each case of n∈{1,…,12}π^+ on the L/a=48 lattice volume as well as the ground-state energy results from all successful fit range choices.
Analogous results for the L/a=32 volume and for n = {1,…,12}K^0 systems on the L/a ∈{32,48} volumes are shown in Appendix <ref>.
The resulting ground-state energies for n∈{1,…,12} meson systems on both lattice volumes are shown in Table <ref>.




 §.§ Hadron mass results


Single-meson masses computed with the gauge field ensembles used here have already been presented in Ref. <cit.>; the meson mass results obtained using the fitting produce employed here are shown for completeness in Fig. <ref>.
The K^0 energy is equal within statistical uncertainties for the L/a = 48 and L/a=32 lattice volumes. 
Volume dependence in E_π^+(L) and differences in the FV single-particle energies E_π^+(L) - E_K^0(L), however, are clearly visible.
Because of the quark mass tuning described in Sec. <ref>, which is designed to remove strong isospin breaking effects, these energy differences are ascribed to QED effects.
In order to interpret these QED effects, results for E_π^+(L) can be compared to the 𝒪(α, L^-2) prediction of QED_L <cit.> or equivalently NRQED_L <cit.>,

    E_π^+(L) = m_π + α/2 L c_1 ( 1 + 2/m_π L) + 𝒪( α^2, L^-3),

where c_1 = -0.266596 is a geometric constant that does not depend on the system under consideration.
Nonlocal effects from zero mode subtraction and charge radius effects introducing a new parameter both arise at next-to-next-to-next-to-leading-order (N^3LO) and are neglected here.
Fitting the L/a∈{32,48} results to Eq. (<ref>) gives 

    a m_π^+ = 0.15419(29),  m_π^+ = 449(1)(13) MeV,

where the lattice spacing a=0.068(2)  fm is determined in Ref. <cit.>. The first uncertainty in each term of Eq. (<ref>) includes statistical and systematic fitting uncertainties as detailed in Appendix <ref>, and the second uncertainty in the expression for m_π^+ arises from the uncertainty in a.
A constant fit to the L/a∈{32,48} results for E_K^0(L) gives 

    a m_K^0 = 0.13918(25),  m_K^0 =  404(1)(12) MeV.

For the values of the quark masses and α≃ 0.1 used in this work, strong isospin breaking effects approximately vanish and the difference between the charged and neutral pseudoscalar meson masses m_π^+ - m_K^0 = 45(2) MeV is attributed entirely to QED effects.
It is noteworthy that QED effects account for approximately 10% of the π^+ mass with these parameters. 









Results for the proton and neutron ground-state energies are shown in Table <ref> and Fig. <ref>.
Volume dependence of the neutron mass is expected to be exponentially suppressed and is found to be consistent with zero within statistical uncertainties.
The proton mass includes power-law FV corrections from QED effects identical to those shown for the π^+ in Eq. (<ref>).
These relativistic FV effects are suppressed by 𝒪((M_p L)^-1) and are therefore suppressed in E_p(L) compared to E_π^+(L) by m_π^+ / M_p ≪ 1.
The proton energy is found to have mild FV effects consistent with zero and with Eq. (<ref>).
Fits to the NLO QED_L prediction for E_p(L) and to a constant for E_n(L) give the results
   
    a M_p    = 0.4037(23),  M_p =  1171(7)(34)  MeV, 
    
          a M_n    = 0.3971(25),  M_n =  1152(7)(34)  MeV,

where the uncertainties are as defined for the π^+ after Eq. (<ref>).
Combining these results gives M_p - M_n = 20(10) MeV.
This result is approximately ten times larger than the QED contribution to the proton-neutron mass difference at the physical values of the quark masses and α <cit.>,
which is consistent with the expected linear dependence of the proton-neutron mass difference on α, given the value α≃ 0.1 that is used here.



§ FINITE-VOLUME NON-RELATIVISTIC QED



Hadronic EFT results relating FV energy levels to the LECs appearing in EFT Lagrangians are useful for interpreting LQCD+QED_L results, as evidenced by the use of fits to Eq. (<ref>) to describe the volume dependence of single-hadron energies and to extract m_π^+ and M_p at L=∞.
Analogous EFT results for charged multi-hadron systems are needed to extract hadronic scattering information from LQCD+QED_L results for multi-hadron FV energy levels, but EFT for charged multi-hadron systems is complicated by the presence of Coulomb interactions that are nonperturbative for hadron pairs with sufficiently small relative momentum as discussed below.
Further complications arise for non-relativistic EFTs for QCD+QED_L because of the nonlocality inherent in zero-mode subtraction.
However, non-relativistic EFTs have the advantage that FV energy shifts for systems of arbitrary particle number can be computed more simply than in relativistic EFTs, where the relation between scattering parameters and FV energy levels is not yet known for systems of four or more hadrons. 
This work therefore pursues the application of non-relativistic EFT to charged multi-hadron FV systems, and this section develops the formalism necessary for studying multi-hadron systems in NRQED_L including nonlocal effects arising from zero-mode subtraction.



 §.§ Finite-volume formalism for two charged hadrons




Interactions of two electrically neutral particles with mass M and relative momentum 2 p̌ are described at low energies by a scattering phase shift δ(p), where p = |p̌|.
The phase shift is an analytic function of the center-of-mass energy E^* = 2 √(p^2 + M^2) for energies below both the t-channel cut and the s-channel inelastic particle production threshold.
For p ≪ M, the energy can be described by the non-relativistic expansion E^*= 2M + p^2/M + …, and below the t-channel cut and inelastic threshold the phase shift admits a convergent effective range expansion p δ(p) = -1/a + r/2p^2 + …, where a is the scattering length (not to be confused with the lattice spacing a appearing elsewhere) and r is the effective range.
For neutral particles, this expansion is straightforwardly reproduced in terms of pionless EFT <cit.>, a theory of hadrons interacting via contact interactions organized in powers of derivatives.

Interactions of electrically charged particles are complicated by the fact that the t-channel cut associated with one-photon-exchange and the inelastic photon production threshold start at p=0.
This leads to p=0 singularities in contributions to the scattering amplitude from Coulomb ladder diagrams describing iterated one-photon-exchange, shown in Fig. <ref>.
Increasingly higher-loop Coulomb ladder diagrams are suppressed by powers of α, but the p=0 singularity becomes more severe.
The loop expansion for Coulomb ladder diagrams corresponds to a perturbative expansion in powers of

    η(p) ≡α/v(p) = α M/2 p,

where v is the relative velocity of the two charged particles and p≪ M is assumed throughout this section.
For η(p) ≳ 1, QED becomes nonperturbative and Coulomb ladder diagrams must be resummed to all orders in α.
As shown in nonrelativistic quantum mechanics by Bethe <cit.>, and in EFT by Kong and Ravndal <cit.>, the resummed scattering amplitude is nonanalytic in η and the effective range expansion is modified as

    (2πη/e^2πη - 1) p δ(p) + 2 p η[ Reψ(i η) - ln(η) ] = -1/a_C + 1/2r_C p^2 + …,

where ψ(x) = Γ^'(x)/Γ(x) is the digamma function, a_C is the Coulomb-corrected scattering length, and r_C is the effective range of the charged particle system.
To leading order in α, and at all orders in η and in the four-particle contact interactions as in pionless EFT, the charged particle effective range is unaffected by QED interactions and r_C = r <cit.>.
Conversely, the contact interaction associated with the scattering length must be renormalized to absorb divergences from Coulomb ladder diagrams.
In the  scheme, the running coupling that would be identified with the scattering length in the absence of QED is related to the Coulomb-corrected scattering length by <cit.>

    1/a^(μ) = 1/a_C + α M [ ln(μ√(π)/α M) + 1 - 3/2γ_E - 1/2μ r ],

and can be understood as the physical scattering length with Coulomb effects from infrared length scales > 1/μ removed.

The FV two-particle energy spectrum can be related to the scattering phase shift in non-relativistic quantum mechanics <cit.> and in quantum field theory <cit.>.
In a finite spatial volume of size L^3 with PBCs, the system exhibits reduced spatial symmetries characterized by covariance under the cubic group, and the momentum carried by a free particle is quantized as p̌ = 2 πň / L with ň∈ℤ_3.
In QED_L, zero-mode subtraction implies that the momentum carried by a photon propagator is restricted to p ≥ 2π / L.
The expansion parameter η(p) in Eq. (<ref>) is therefore restricted to

    η(p) ≤η_L ≡α M L/4 π

in QED_L Coulomb ladder diagrams.
For sufficiently small volumes, η≤η_L ≪ 1 and Coulomb photon exchange can be treated perturbatively.

The 𝒪(α) quantization condition relating the s-wave scattering amplitude to the two-particle spectrum in the A_1^+ representation of the cubic group in non-relativistic EFT in the approximation of negligible partial wave mixing was derived in Ref. <cit.>,

    (2πη/e^2πη - 1) p δ(p) + 2 p η[ Reψ(i η) - ln(η) ] = 1/π L𝒮^C(p) + α M [ ln( 4π/α M L) - γ_E ],

where η(p) is defined in Eq. (<ref>) and

   
    𝒮^C(x)    = 𝒮(x) - α M L/4π^3𝒮_2(x) + α M a_C^2 r_C/π^2 L^2ℐ 𝒮(x)^2 + O(α^2), 
    𝒮(x)    = lim_Λ_n →∞∑_ň∈ℤ^3^|ň|<Λ_n1/|ň|^2 - x^2 - 4πΛ_n, 
    𝒮_2(x)    = lim_Λ_n →∞∑_ň∈ℤ^3^|ň|<Λ_n∑_m̌∈ℤ^3 ∖{ň}1/|ň|^2 - x^21/|m̌|^2 - x^21/|ň - m̌|^2 - 4π^4 lnΛ_n ,


where ℐ≈ -8.9136 is a geometric constant whose evaluation is detailed in  Refs. <cit.> and sums over integer triplets ň are restricted to |ň| ≤Λ_n where Λ_n is a cutoff and the Λ_n →∞ limit should be taken as indicated.
LQCD+QED_L results for the p associated with FV energy levels below three-particle thresholds can be related to charged particle scattering phase shifts by Eq. (<ref>) and used to constrain parametrizations of the phase shift such as Eq. (<ref>).
Since Eq. (<ref>) neglects exponentially small FV effects present in small volumes, as well as (α M L)^n effects from Coulomb ladder diagrams with n photon propagators that must be resummed for sufficiently large volumes, Eq. (<ref>) is necessarily valid only for an intermediate range of L.
An important goal of this work is to test the range of volumes over which Eq. (<ref>) can be reliably used to extract Coulomb-corrected scattering parameters from LQCD+QED_L results.



 §.§ Charged two-hadron systems in NRQED_L


Eq. (<ref>) can be perturbatively expanded in powers of a_C / L and other higher-order effective range expansion coefficients.
Ref. <cit.> determined this expansion to 𝒪(α  M L) and 𝒪(a_C / L)^3 under the assumption that α M L ≪ a_C / L ≪ 1.
Following  Refs. <cit.>, the same result can be derived in non-relativistic quantum mechanics using a Hamiltonian that perturbatively includes the effects of relativity and allows straightforward generalization to many-particle systems.
The NRQED_L Lagrangian is given by

    ℒ   = ψ^†( i D_0 - Ď^2/2M) ψ -  1/2( 4π a/M) (ψ^†ψ)^2  - η_3(μ)/3! (ψ^†ψ)^3 + ℒ_γ^ξ + ℒ_r.

In this expression ψ is a non-relativistic hadron field, D_μ = ∂_μ + i e Q A_μ where Q is the electric charge operator, the gauge-fixed photon Lagrangian ℒ_γ^ξ is given in Eq. (<ref>), η_3(μ) is a renormalization-scale-dependent coupling associated with short-range three-body interactions, and four- and higher-body interactions are neglected.
ℒ_r includes effective range contributions and relativistic corrections involving two derivatives, and is given by

    ℒ_r    = ψ^†(  Ď^4/8M^3) ψ -  1/2( π a/M) ( a r - 1/M^2) (ψ^†ψ)(ψ^†Ď^2 ψ + ψĎ^2 ψ^†).

The coefficient of the (ψ^†ψ)(ψ^†𝐃^2 ψ + ψ𝐃^2 ψ^†) operator in Eq. (<ref>) can be fixed by demanding that the strong interaction EFT given by replacing D_μ with ∂_μ reproduces Eq. (<ref>) with α = 0 when both are expanded perturbatively in powers of L^-1 to[The Lagrangian in Eq. (<ref>) can also be obtained by studying the nonrelativistic limit of relativistic scalar field theory as in Refs. <cit.> .] 𝒪(L^-6).
This 𝒪(L^-6) threshold expansion of Lüscher's quantization condition is given in Ref. <cit.> and is verified below to be reproduced by the non-relativistic EFT defined by Eqs. (<ref>)-(<ref>)[Refs. <cit.> include the operators in Eq. (<ref>), but a factor of 2 discrepancy in the coefficient of the last term leads to a difference in the 𝒪(L^-6) threshold expansion result.].
Operators in ℒ_r lead to contributions to the threshold expansion suppressed by ar/L^2 as well as relativistic effects suppressed by 𝒪( (ML)^-2) that will be neglected in the power counting schemes discussed below. 
Additional relativistic corrections to Eq. (<ref>) arise from photon loops and operators with four and more derivatives, but these give rise to FV effects suppressed by powers of 𝒪( (ML)^-1) that will be neglected  as discussed below and detailed in Appendix <ref>. 

Introducing Fourier transformed fields

    ψ_ǩ = 1/L^3/2∫ d^3x  e^i ǩ·x̌ψ(x̌), ψ(x̌) = 1/L^3/2∑_ǩ e^-i ǩ·x̌ψ_ǩ,

the FV Hamiltonian for NRQED_L is given by

    H     = ∑_ǩψ^†_ǩ( ǩ^2/2M) ψ_ǩ + H_γ^ξ + H_ int,

where H_γ^ξ is the gauge-fixed photon Hamiltonian, whose explicit form will not be used below, and

    H_ int   =  - ∑_ǩψ^†_ǩ( ǩ^4/8M^3) ψ_ǩ + 1/2 L^3∑_p̌^', p̌, Q̌ V(p̌,p̌^') ψ^†_Q̌-p̌^'ψ^†_Q̌+p̌^'ψ_Q̌-p̌ψ_Q̌+p̌
        + η_3(μ)/(3!) L^6∑_Q̌,p̌,p̌^',q̌,q̌^'ψ_Q̌+p̌^'^†ψ_Q̌+q̌^'^†ψ_Q̌-p̌^'-q̌^'^†ψ_Q̌+p̌ψ_Q̌+q̌ψ_Q̌-p̌-q̌,

where the two-body potential includes strong interaction and Coulomb terms,

    V(p̌^',p̌)     = 4π/M(a + a/4(a r - 1/M^2) (p̌^2 + p̌^' 2) + …) + 4πα/|p̌^' - p̌|^2( 1 - δ_p̌,p̌^').

Relativistic corrections to the potential, including photon loop effects as well as nonlocal effects of zero-mode removal, can in principle be calculated by carrying out the FV analog of higher-order matching between QED, NRQED, and potential NRQED  (pNRQED) <cit.>, but these terms give rise to 𝒪(ML)^-1 effects that are neglected below.

Specializing to the case of identical bosons, operators are associated with the fields in Eq. (<ref>) and are defined to satisfy commutation relations [ψ^†_p̌^' ,ψ_p̌ ] = δ_p̌,p̌^'.
FV two-particle states defined by

    |p̌_1, p̌_2⟩ = 1/√(2)ψ_p̌_1^†ψ_p̌_2^†|0⟩

satisfy ⟨p̌_1,p̌_2|p̌_1,p̌_2⟩ = 1.
The ground state of the two-particle system in its center-of-mass frame is |p̌_1 = 0̌, p̌_2 = 0̌⟩.
Splitting the Hamiltonian into kinetic energy terms and an interaction term H_ int that will be treated perturbatively, the ground-state FV energy shift at leading order in Rayleigh-Schrödinger perturbation theory is then given by

    Δ E^LO = 0̌,0̌H_ int0̌,0̌ = 1/L^3 V(0̌,0̌)  = 4π a/ML^3 .

Standard perturbation theory techniques can be used to extended this result to higher orders in a/L and α.
At next-to-next-to-leading order (NNLO) in H_ int, the result is given by

    Δ E^NNLO,PC1   = 4π a_C/ML^3{ 1  -  ( a_C/π L)   ℐ + ( a_C/π L)^2  [ ℐ^2    -   𝒥]  }
         + α/π L{ - ( a_C/π L) 2𝒥 + ( a_C/π L)^2 { 2 ℐ𝒥 + ℛ_22 - 2 𝒦 + 4π^4[ ln( 4π/α M L) - γ_E ] }. 
       . - ( α M L/4π^3) 𝒦  + ( a_C/π L)( α M L/4π^3)[ 2 ℛ_24 + 𝒥^2 - ℒ] + ( α M L/4π^3)^2 ℛ_44} ,


where PC1 is a label for the power counting scheme discussed below and higher-order corrections in a_C/L ∼ r/L, as well as relativistic effects suppressed by 𝒪( (ML)^-1), have been neglected.
All sum-integral differences ℐ, 𝒥, 𝒦, ℒ, and ℛ_nm which appear in this expression are evaluated[See also Refs. <cit.> for more details on evaluating FV sums.] in Ref. <cit.>, which includes an evaluation of Eq. (<ref>) in EFT to leading order in α M L and to higher order in a/L, except for ℛ_44 which is a convergent sum given by ℛ_44≈ 55.47.

Eq. (<ref>) can be expressed as an expansion in a_C / L and the FV Coulomb parameter η_L defined in Eq. (<ref>),

    Δ E^NNLO,PC1   = 4π a_C/ML^3{ 1  -  ( a_C/π L)   ℐ + ( a_C/π L)^2  [ ℐ^2    -   𝒥]  }
         + 4 η_L/ML^2{ - ( a_C/π L) 2𝒥 + ( a_C/π L)^2 [ 2 ℐ𝒥 + ℛ_22 - 2 𝒦 - 4π^4[ ln( η_L ) + γ_E ] ] . 
       . - ( η_L/π^2) 𝒦  + ( a_C/π L)( η_L/π^2)[ 2 ℛ_24 + 𝒥^2 - ℒ] + ( η_L/π^2)^2 ℛ_44} .

This resembles a double power series expansion in the parameters η_L and a_C/L, suggesting that the NRQED_L threshold expansion should provide a good approximation of FV energy shifts for

    PC1: η_L ∼a_C/L≪ 1  .

It is noteworthy that η_L ≪ 1 only requires α M L ≪ 4π and is less restrictive than the condition α M L ≪ 1 discussed in Ref. <cit.>.
In the matching to the LQCD+QED_L simulations discussed below, the power counting η_L ∼ a_C / L numerically overestimates the size of QED effects, particularly on the smaller volumes studied, and the alternative power counting[The description of FV effects on charged hadron masses in QED_L as a dual expansion in α and (ML)^-1 and the possibility of using alternative power countings in LQCD+QED_L calculations depending on the values of these parameters is explored in Ref. <cit.>.]

    PC2: η_L^1/2∼a_C/L≪ 1,

will also be used in fits to LQCD+QED_L results.
Higher-order results for the threshold expansion for short-range contact interactions without QED <cit.> can be used to extend Eq. (<ref>) from NNLO in Eq. (<ref>) to N^3LO in Eq. (<ref>), 

    Δ E^N^3LO,PC2   = 4πa_C/ML^3{ 1  -  ( a_C/π L)   ℐ + ( a_C/π L)^2  [ ℐ^2    -   𝒥]  - ( a_C/π L)^3 [ ℐ^3 - 3 ℐ𝒥 + 𝒦] }
         + 4η_L/ML^2{ - ( a_C/π L) 2𝒥 - ( η_L/π^2) 𝒦. 
       . + ( a_C/π L)^2 [ 2 ℐ𝒥 + ℛ_22 - 2 𝒦 - 4π^4[ ln( η_L ) + γ_E ] ] }.

The parameter a_C(L) is equal to a_C plus a 1/L^3 suppressed correction arising from the interaction terms in Eq. (<ref>), 

    a_C = a_C(L) - 2πa_C(L)^2/L^3( a_C(L) r - 1/2M^2),

where (ML)^-2 suppressed effects are shown for completeness and lead to agreement with the 𝒪(L^-6) strong interaction relativistic threshold expansion of Ref. <cit.> after taking[The operators in Eq. (<ref>) lead to additional effective range and relativistic corrections to the right-hand-side of Eq. (<ref>), namely an additional term of the form -2a_C^2 r η_L ℐ / (ML^5) - η_L^2 𝒥 / (4π^2 M^3 L^4), but these are higher order than the effective range term in Eq. (<ref>) according to the power counting of Eq. (<ref>).]  α = 0.

In principle, LQCD+QED_L results for π^+π^+ FV energy shifts on multiple lattice volumes could be used to extract both a_C and r by constraining the 𝒪(L^-6) difference between a_C and a_C.
In the LQCD+QED_L calculations discussed below, (a_C(L) - a_C)/a_C can be estimated at LO in chiral perturbation theory (χPT) to be 2% and 8% for the L/a=48 and L/a=32 lattice volumes respectively.
To see whether a_C - a_C can be reliably determined, this estimate must be compared with an estimate of relativistic effects neglected in Eq. (<ref>), which as discussed in Sec. <ref> below modify the dominant QED FV effects by 𝒪( (ML)^-1).
For π^+π^+ systems on the L/a =48 lattice volume, the dominant (NLO) QED effect in Eq. (<ref>) amounts to a shift a_C →a_C - 2𝒥 (η_L / π^2) a_C, and radiation photon effects can be estimated to lead to a 𝒪( (ML)^-1) suppressed shift in (a_C(L) - a_C)/a_C of order   ∼ 2𝒥 (η_L / π^2)/(ML) = 𝒥α/(2π^3) ∼ 3% for α≃ 0.1.
Relativistic effects are therefore comparable to a_C - a_C for π^+π^+ systems and prevent the effective range contribution to a_C(L) from being disentangled from other FV effects neglected in Eq. (<ref>).
Therefore, a_C(L) - a_C will be neglected when fitting π^+π^+ LQCD+QED_L results to Eq. (<ref>).







 §.§ Zero-mode effects



The derivation of Eq. (<ref>) in Ref. <cit.> and the form of the NRQED_L Lagrangian in Eq. (<ref>) assume that charged particle NRQED_L contact interactions in FV are equal to their infinite-volume counterparts up to exponentially suppressed corrections.
More recently, however, it has been shown in Refs. <cit.> that this assumption is violated in the single-particle sector of NRQED_L because of the inherent nonlocality of zero-mode subtraction.
NRQED_L parameters can be obtained by calculating masses, scattering amplitudes, or other observables in both NRQED_L and QED_L and tuning the parameters of the NRQED_L Lagrangian to reproduce QED_L results.
In QED_L, on-shell photon exchange leads to power-law FV effects on charged particle masses suppressed by powers of α and 1/(ML) that have been studied in Refs <cit.>.
The 𝒪( α / (ML) ) and 𝒪( α / (ML)^2) corrections are independent of the structure of the charged particle and are described by one-loop diagrams in both QED_L and NRQED_L.
At order 𝒪(α / (ML)^3), structure-dependent effects involving magnetic moments and charge radii arise. 
Nonlocal effects from zero-mode subtraction also enter at order 𝒪( α / (ML)^3) because zero-mode subtraction leads to power-law FV effects from off-shell antiparticle modes in QED_L that are not reproduced by NRQED_L loop diagrams.
These effects can be included in NRQED_L by adjusting the Lagrangian to include additional particle-antiparticle interactions <cit.>, or more simply by adjusting the coefficients of mass operators in the NRQED_L Lagrangian by factors proportional to α/(ML)^3 <cit.>.
For charged scalars, although not for charged fermions, these effects vanish in the charged particle rest frame <cit.>.
This non-decoupling of antiparticle modes is a consequence of the nonlocality of NRQED_L, and for EFTs with breakdown scale Λ generically produces 𝒪( α / (Λ L)^3 ) corrections to LECs <cit.>.

Nonlocal effects from zero-mode subtraction could lead to power-law FV effects that modify four-hadron contact interaction couplings proportional to a_C in NRQED_L.
Considering the 𝒪(ML) enhancement of FV effects associated with Coulomb ladder diagrams, it is necessary to analyze nonlocal effects of zero-mode subtraction on a_C in order to determine whether Eq. (<ref>) is modified within the order of approximation considered.
The effects of zero-mode subtraction on a_C can be determined by matching any QED_L and NRQED_L correlation functions sensitive to four-hadron contact interactions.
One-particle FV self-energies only receive contributions from four-hadron contact interactions in 𝒪(α^2) diagrams containing closed loops of particle-antiparticle pairs.
Two-particle FV Green's functions receive contributions from four-hadron contact interactions at 𝒪(α^0), and it is convenient to calculate nonlocal FV effects on four-hadron contact interactions in NRQED_L by directly matching two-particle FV Green's functions in QED_L and NRQED_L.
This matching is detailed in Appendix <ref> and summarized below.






The 𝒪(ML) enhancement of Coulomb ladder diagrams arises when one intermediate particle propagator is placed on-shell and another intermediate particle propagator is nearly on-shell with virtuality k^2 / M, as compared with virtuality k when a photon propagator is placed on-shell.
This can be explicitly seen by comparing the integrands of the “box” diagram (rightmost top row in Fig. <ref>) with the “crossed-box” diagram (leftmost bottom row in Fig. <ref>).
In NRQED_L, the box diagram involves the integral


    i∫dk^0/2π( 1/k^0 - k^2/2M + i ϵ) ( 1/-k^0 - k^2/2M + i ϵ) (1/(k^0)^2 - k^2 + iϵ)^2 
       = -( M/k^2) (1/(k^2/2M)^2 - k^2) + d/dk^0[ (1/k^0 - |k|)^2 ( 1/k^0 - k^2/2M) ( 1/-k^0 - k^2/2M) ]_k^0 → -k
       = -M/k^6[ 1 + 𝒪(k^2/M^3) ] - 3/4k^5[ 1 + 𝒪(k^2/M^2) ].


In this expression, the first contribution involving the particle pole places the second particle propagator nearly on-shell with kinetic energy k^2/2M. 
The second contribution from the photon double pole gives both particle propagators off-shell kinetic energies k.
When FV effects are computed, k is replaced by the quantized values 2π n /L with n∈ℤ_3, and (after adding all necessary UV counterterms) amplitude suppression by powers of k/M implies suppression of FV effects by the corresponding power of (ML)^-1.
The NRQED_L crossed-box diagram involves the integral


    i∫dk^0/2π( 1/k^0 - k^2/2M + i ϵ)^2 ( 1/(k^0)^2 - k^2 + iϵ)^2 
       = d/dk^0[ (1/k^0 - |k|)^2 (1/k^0 - k^2/2M)^2 ]_k^0→ -k
       = 3/4k^5[ 1 + 𝒪(k^2/M^2) ] ,


where only the photon pole contributes and leads to particle propagators with off-shell kinetic energies k. FV effects associated with the crossed-box diagram are therefore suppressed by 𝒪( (ML)^-1) compared to the dominant contribution of the box diagram.
 


The 𝒪(ML) power enhancement of the box diagram is only present in the diagrams involving repeated s-channel Coulomb interactions shown in Fig. <ref> but not in the diagrams in Figs. <ref>-<ref> or other diagrams involving particle-antiparticle pair creation that vanish non-relativistically[More details on the power-counting of diagrams involving FV photon exchange are given in Ref. <cit.> and discussions of analogous power counting arguments for radiation pions are given in Ref. <cit.>.].
Furthermore, the 𝒪(ML) enhancement only occurs when the intermediate-state charged particles are both nearly on-shell.
As detailed in Appendix <ref>, power-law FV effects in QED_L that are not reproduced by loop diagrams with the leading order NRQED_L Lagrangian arise from antiparticle poles where intermediate states have a large virtuality of order 2M.
These do not receive the 𝒪(ML) enhancement of particle pole contributions, and zero-mode effects are found to be suppressed by 𝒪(α / (ML)^3).
Matching between QED_L and NRQED_L is explicitly performed for charged scalars in Appendix <ref>, and zero-mode effects are found to modify the four-scalar coupling in NRQED_L at order 𝒪( α/(ML)^3) for boosted systems but not to modify the coupling for scalars at rest at this order.
This shows that Eq. (<ref>) is valid for charged scalars in NRQED_L up to 𝒪(α^2) effects and 𝒪(ML)^-1 relativistic effects.
The vanishing of these effects in the center-of-mass frame arises from cancellations due to the specific form of the scalar-photon vertex functions.
It is possible that for charged fermions nonlocal effects also arise at 𝒪( α / (ML)^3) in the charged fermion rest frame.
The 𝒪( α/(ML)^3 ) suppression factor is consistent with the physical arguments of Ref. <cit.> that nonlocal effects arise from interactions between the subtracted zero-mode, which can be interpreted as a uniform background charge density that ensures Gauss's law is satisfied for the FV system <cit.>, and high-energy modes that have been integrated out of the EFT, and therefore that FV effects from zero-mode subtraction are suppressed by α times the inverse volume.











 §.§ Charged many-hadron systems in NRQED_L




The two-particle energy shifts in Eqs. (<ref>)-(<ref>) can be extended to a threshold expansion for FV effects on systems of n non-relativistic particles using Rayleigh-Schrödinger perturbation theory.
Unit-normalized many-particle states are given by


    |p̌_1,…,p̌_n⟩ = 1/√(n)ψ^†_p̌_1×…×ψ^†_p̌_2|0⟩,


and the leading order FV energy shift for the ground state of n identical bosons in the center-of-mass frame is given by


    0̌,…0̌H_ int0̌,…0̌ = 1/L^3n  2V(0̌,0̌) + 1/L^6n 3η_3(μ) .


Working to NNLO in the power counting of Eq. (<ref>), the energy shift of an n-hadron state is equal to n  2 times the two-body energy shift of Eq. (<ref>), plus additional contributions from induced three-body and four-body forces shown in Fig. <ref>.
Additional diagrams associated with radiation photon exchange shown in Fig. <ref> are suppressed by 𝒪( (ML)^-1).
The resulting threshold expansion for the FV energy shift for a system of n like-charged hadrons at rest is given by


    Δ E_n^NNLO,PC1   = 4π a_C/ML^3n  2{ 1  -  ( a_C/π L)   ℐ + ( a_C/π L)^2  [ ℐ^2    + (2n-5)   𝒥]  }
         + 4 η_L/ML^2n  2{ - ( a_C/π L) 2𝒥    - ( η_L/π^2) 𝒦 + ( η_L/π^2)^2 ℛ_44. 
       .  + ( a_C/π L)^2 [ 2 ℐ𝒥 + ℛ_22 + (4n - 10) 𝒦 - 4π^4[ ln( η_L ) + γ_E ] ] . 
       .   + ( a_C/π L)( η_L/π^2)[ 2 ℛ_24 + 𝒥^2 + (2n - 5) ℒ] } ,


where omitted terms are: quartic or higher in η_L ∼ a_C / L; relativistic effects suppressed by 𝒪( (ML)^-1); or three-body contact interactions where η_3 ∼ a_C^4/M is assumed so that 𝒪(L^-6) terms of the strong interaction threshold expansion appear at the same order.




At N^3LO in the power counting of Eq. (<ref>) there are contributions from short-distance three-body interactions well as the induced few-body interactions discussed above.
These introduce a new free parameter η_3(μ) that parametrizes the strength of six-particle operators in the NRQED_L Hamiltonian, Eq. (<ref>).
Combining Eq. (<ref>) with the N^3LO results for the QCD threshold expansion from Ref. <cit.>, the N^3LO FV energy shift for a system of n like-charged hadrons in the power counting of Eq. (<ref>) is given by


    Δ E_n^N^3LO,PC2   = 4πa_C/ML^3n  2{ 1  -  ( a_C/π L)   ℐ + ( a_C/π L)^2  [ ℐ^2    + (2n-5)   𝒥]  . 
       . - ( a_C/π L)^3 [ ℐ^3 + (2n - 7) ℐ𝒥 + (5n^2 - 41n + 63)𝒦] }
         + 4η_L/ML^2n  2{ - ( a_C/π L) 2𝒥  - ( η_L/π^2) 𝒦. 
       . + ( a_C/π L)^2 [ 2 ℐ𝒥 + ℛ_22 + (4n - 10) 𝒦 - 4π^4[ ln( η_L ) + γ_E ] ] }
        + n  31/L^6[ η_3(μ) + 64πa_C^4/M (3√(3) - 4π) ln(μ L) - 96 a_C^4/π^2 M𝒮_ MS].



The non-relativistic EFT three-body coupling η_3(μ) is renormalization scheme  and scale dependent.
The scale-dependence of η_3(μ) cancels the explicit ln(μ L) scale-dependence shown in Eq. (<ref>), and the scheme-dependence is compensated by scheme-dependence in the finite term 𝒮_ MS.
This scale-dependence arises from the ambiguity in separating short-distance three-body interactions described by contact operators in NRQED_L from long-distance two-body rescattering effects.
In relativistic theories this ambiguity does not arise, and the particle mass plays the role of the scale μ in relativistic descriptions of the L^-6ln L term in the three-body threshold expansion derived for generic relativistic field theories in Ref. <cit.>.
Equating the 𝒪(L^-6) term in the relativistic threshold expansion of the three-particle threshold amplitude ℳ_ 3,th in Ref. <cit.> with the corresponding 𝒪(L^-6) nonrelativistic threshold expansion provides a relation between the non-relativistic coupling η_3(μ) and the scale-independent ℳ_ 3,th in the absence of QED,


    η_3(μ)    = -ℳ_ 3,th/48 M^3 + 64 π a^4/M(3√(3) - 4π)ln(M/2πμ)  + 48 a^2 π^2/M^3 + 48 a^3 π^2 r/M
       + 12 a^4/π^2 M𝒮_3 + 768 π^3 a^3/M^2𝒞_3,


where a is the scattering length for a neutral two-particle system, 𝒞_3 = -0.05806 is a FV sum evaluated in Ref. <cit.>, and 𝒮_3 = 571.398 is related to 𝒮_ MS, evaluated in Ref. <cit.>, and to other FV sums from Ref. <cit.> by 𝒮_3 = 𝒞_F + 𝒞_4 + 𝒞_5 + 8𝒮_ MS.
QED effects will modify Eq. (<ref>), but these modifications can be neglected at the EFT order considered here.
Below, QED effects on three-body forces will be studied by comparing the three-body interaction parameters extracted from LQCD+QED_L results for systems of charged and neutral mesons.



§ RESULTS FOR CHARGED MULTI-HADRON SYSTEMS
 

This section combines the LQCD+QED_L results from Section <ref> with the NRQED_L results from Section <ref> in order to obtain  QCD+QED predictions for scattering lengths and other hadronic interaction parameters at the values of the quark masses and α used here.



 §.§ Charged meson scattering










The LQCD+QED_L results for the FV spectrum results in Table <ref> can be used to constrain the low-energy EFTs for charged and neutral meson interactions.
It is convenient to focus on results for the FV energy shifts

    Δ E_nM(L) = E_nM(L) - n E_M(L), M ∈{K^0,  π^+}.

Results for correlated differences between n-particle ground-state energies and n times the one-particle ground-state energy, as defined in Eq. (<ref>), are more precise than n-particle energies alone. Furthermore, this subtraction nonperturbatively removes single-particle FV effects from n-meson FV energy shift results.
For multi-π^+ systems, LQCD+QED_L results for these FV energy shifts can be identified with the interaction energy shifts Δ E_n computed perturbatively in NRQED_L in Sec. <ref>.
For multi-K^0 systems, LQCD+QED_L results can be identified with the same EFT results after setting α to zero. 

In the numerical LQCD+QED_L calculation, FV energy shifts are computed in a correlated manner using bootstrap resampling as detailed in Appendix <ref>, and the results are shown in Table <ref>.
To access QED-specific effects, the differences of these differences between the n charged pions and n neutral kaon systems are also computed similarly.
The double subtraction suppresses any strong isospin breaking effects arising from mistuning of the quark masses for different charge quarks.
In the numerical LQCD+QED_L calculation, correlated differences of FV energy shifts are computed using bootstrap resampling as detailed in Appendix <ref>, and the results are shown in Table <ref>.






Results for the two-particle FV energy shifts Δ E_π^+π^+ and Δ E_K^0 K^0 are shown in Fig. <ref>.
Both energy shifts are clearly resolved from zero with relative uncertainties in the range of 15 - 25% for both volumes, although Δ E_π^+π^+ and Δ E_K^0 K^0 on a given volume are indistinguishable.
The small magnitude of QED effects on Δ E_π^+π^+ might appear surprising because α m_π L ∼ 0.74 for this volume, but as discussed in Sec. <ref> the appropriate FV analog of the Coulomb expansion parameter is η_L = α m_π L / (4π) ∼ 0.06 for the L/a = 48 volume.
Eqs. (<ref>)-(<ref>) therefore predict that in addition to differences arising from a_C^π^+π^+≠ a^K^0K^0, NLO corrections from Coulomb photon exchange modify the LO FV energy shift by ∼ 20% on the L/a=48 lattice volume, which is not expected to be distinguishable given the statistical uncertainties on Δ E_π^+π^+ and Δ E_K^0 K^0.
This expectation is consistent with LQCD+QED_L results, as shown in Fig. <ref>.





The scattering lengths a_C^π^+π^+ and a^K^0K^0 can be extracted from a combined fit to the results for Δ E_π^+π^+ and Δ E_K^0 K^0 shown in Table <ref> and the results for the precisely determined correlated differences Δ E_π^+π^+ - Δ E_K^0 K^0 shown in Table <ref>.
Fitting to Eq. (<ref>), as is appropriate for the power counting PC1 of Eq. (<ref>), gives the results,


    NNLO, PC1: a^K^0K^0 m_K^0 = 0.335(26),   a^π^+ π^+_C m_K^0 = 0.463(41),

   
where the common scale m_K^0 has been included for both π^+π^+ and K^0K^0 to facilitate comparison of a^K^0K^0 and a^π^+π^+_C.
The lowest-order QED effect on Δ E_π^+π^+ from Coulomb photon exchange in Eq. (<ref>) decreases Δ E_π^+π^+ compared to the FV energy shift for neutral particles.[Coulomb photon exchange leads to a decrease in the energy of a system of π^+ mesons in QED_L because of zero-mode subtraction. Physically, the energy decrease can be understood as arising from attraction between the charged particle system and the uniform background of opposite charge associated with zero-mode subtraction <cit.>. Formally, zero-mode subtraction removes the LO one-photon-exchange diagrams associated with repulsion between charged particles. The dominant QED contribution therefore arises at NLO and necessarily lowers the ground-state energy since it appears at second order in perturbation theory.]
The scattering length results in Eq. (<ref>) show that this effect from Coulomb photon exchange competes with additional QED effects that lead to a_C^π^+π^+ > a^K^0K^0.
Fitting to Eq. (<ref>), as is appropriate for the power counting PC2 of Eq. (<ref>), gives consistent results,


    N^3LO, PC2: a^K^0K^0 m_K^0 = 0.332(27),   a^π^+ π^+_C m_K^0 = 0.465(42),


demonstrating that the fit is not overly sensitive to higher-order terms absent in one or the other power counting.

Both Eq. (<ref>) and Eq. (<ref>) neglect relativistic effects from radiation photon exchange leading to 𝒪( (ML)^-1) FV effects.
   These effects are estimated in Sec. <ref> to lead to a shift in a_C^π^+π^+ of order ∼ 3%, which is smaller than the 6 - 9% statistical uncertainty on a_C^π^+π^+ in Eqs. (<ref>)-(<ref>) and can be consistently neglected.




 §.§ Charged multi-nucleon systems





Two-proton states receive QED_L FV effects from Coulomb photon exchange proportional to α M_p L that are enhanced compared with those in the π^+ π^+ case discussed above.
For both lattice volumes α M_p L > 1, and according to the scaling estimates of Ref. <cit.> Coulomb effects should be nonperturbative.
However, η_L^p = α M_p L / (4π) ∼ 0.15 for the L/a=48 lattice volume and the NRQED_L results in Sec. <ref> expressed as power series in η_L^p show signs of convergence.
   Examining the QED_L contributions proportional to a_C^pp / L in Eq. (<ref>), the NLO mixed strong-Coulomb contribution is suppressed compared to the LO strong contribution by (η_L^p / π^2)2𝒥∼ 0.51, while the corresponding NNLO contribution is suppressed compared to the LO contribution by (η_L^p/π^2)(2ℛ_24 + 𝒥^2 -ℒ) ∼ 0.15.
This suggests that Coulomb effects should be perturbative and subdominant compared to strong-interaction FV effects, which are enhanced by the large size of baryon-baryon scattering lengths <cit.>.
The quantization condition in Eq. (<ref>), which neglects 𝒪( (η_L^p)^2 ) perturbative Coulomb effects and (ML)^-1 relativistic effects but is nonperturbative in strong interaction effects, is therefore needed to relate the pp FV energy shifts to the infinite volume pp phase shift and determine a_C^pp.








   The two-nucleon isospin I=1 systems pp, nn, np(^1S_0), as well as the deuteron, are studied on both lattice volumes.
   Relatively clean signals are seen for each system, and their ground-state energies are determined with total (statistical plus fitting systematic) uncertainties at the 2% level as shown in Table <ref>.
   FV energy-level shifts are determined from combined analyses of the two-nucleon and single-nucleon correlation functions as described in Appendix <ref>, and fit results for all systems are shown in Appendix <ref>.
   As shown in Table. <ref> and Fig. <ref>, the statistical precision of this calculation is insufficient to resolve either the proton-proton or neutron-neutron FV energy shift from zero on either volume studied.
   Resolving non-zero FV shifts of the 𝒪(10  MeV) size expected for two-nucleon systems without QED at these quark masses at a 95% confidence level for the pp, nn, and np systems on the L/a=32 lattice volume would require statistical ensembles approximately ∼ 50 - 100 times larger  than the one used here, as estimated by extrapolating the uncertainties of the four different two-nucleon systems considered assuming 1/√(N) scaling of uncertainties.
Determination of a_C^pp through the QED_L quantization condition of Eq. (<ref>) is therefore left to future work.

The two I=1/2 three-nucleon systems ^3He and ^3H are also investigated, and their ground-state energies are given in Table <ref>.
Correlated differences between three-nucleon ground-state energies and the sums of their constituent nucleon masses are show in Table <ref>.
As with the two-nucleon systems, fit results are shown in Appendix <ref>.
The results for ^3He and ^3H are not precise enough to allow FV effects to be reliably determined.
Precision in the three-nucleon sector is significantly worse than in the two-nucleon sector, as expected. 
   The absolute size of FV energy shifts is also expected to be larger for three-nucleon systems than two-nucleon systems, and for instance resolving an 𝒪(50 MeV) FV energy shift at a 95% confidence level for the ^3He and ^3H systems on the L/a=32 lattice volume would require a statistical ensemble approximately ∼ 100 times larger than the one considered here, based on an extrapolation analogous to that described for the two-nucleon case.
   
Future high-precision LQCD+QED_L calculations of these multi-nucleon systems will provide insight into QED effects on nucleon-nucleon and three-nucleon interactions through a determination of the ^3He - ^3H binding-energy difference and its decomposition into QED and strong isospin breaking effects from LQCD+QED_L.




 §.§ Systems of many charged mesons




Multi-pion correlation functions do not suffer from significant exponential signal-to-noise degradation with increasing particle number and can be used to study QED in the regime where the charge Z ∼ 1/α.
In particular, the correlation functions for systems with n ≤ 12 π^+ mesons described in Sec. <ref> can be used to study systems with Z α≤ 1.2, reaching a charge density of n / L^3 ∼ 1.2  fm^-3.
The dominant strong interactions and QED effects on many-particle FV energy shifts in Eqs. (<ref>)-(<ref>) both scale with n^2, and both Δ E_n π^+ and Δ E_n K^0 can be extracted for larger n with better relative precision than from the n=2 case discussed in Sec. <ref>.


Results for Δ E_n π^+, Δ E_nK^0, and the correlated differences (Δ E_nπ^+ - Δ E_nK^0) for n∈{2,…,12} are shown in Tables <ref>-<ref>.
QED effects leading to non-zero (Δ E_n π^+ - Δ E_n K^0) can be resolved to better than 1 σ on the L/a=32 lattice volume for 3 ≤ n ≤ 7, and on the L/a = 48 lattice volume for n ≥ 9.
These 33 FV energy shifts and correlated differences Δ E_n π^+, Δ E_nK^0, and (Δ E_nπ^+ - Δ E_nK^0) can be used to constrain the low-energy interaction parameters { a_C^π^+π^+, a^K^0K^0, η_3^π^+π^+π^+(m_K^0), η_3^K^0K^0K^0(m_K^0) } appearing in Eqs. (<ref>)-(<ref>).[The renormalization scale used to evaluate η_3(μ) should be chosen close to the “high” energy scale where NRQED_L is matched to QED_L in order to avoid large logarithms that can worsen EFT convergence. For simplicity, μ = m_K^0 is used as the renormalization scale for η_3 throughout this work.]
Fits to the NNLO expression given in Eq. (<ref>) in PC1, which includes 𝒪(η_L^3) Coulomb effects but neglects three-body forces, underpredict LQCD+QED_L energy-shift results for n ≳ 8 meson systems on both lattice volumes and obtain a minimum χ^2 / N_dof∼ 1.3.
The N^3LO expression Eq. (<ref>) in PC2 includes additional free parameters related to three-body forces not present at NNLO.
Including three-body force parameters improves the quality of the fit, and Eq. (<ref>) provides a better description of the LQCD+QED_L results with χ^2 / N_dof∼ 0.8.
Fit results using the N^3LO expression Eq. (<ref>) and uncertainties computed using bootstrap resampling of the global fitting procedure are compared to LQCD+QED_L results for the π^+ and K^0 energy shifts in Figs. <ref>-<ref>.
The results for the meson scattering lengths are consistent with, but more precise than, the results obtained from the two-meson FV energy shifts alone,


    N^3LO, PC2: a^K^0K^0 m_K^0 = 0.337(19),  a^π^+ π^+_C  m_K^0 = 0.464(29).


It is noteworthy that results with Z α≥ 1 can be fit by the NRQED_L formula given in Eq. (<ref>) without additional modifications to account for relativistic QED effects or additional nonperturbative effects.
Some tensions between LQCD+QED_L results and N^3LO NRQED_L fits can be observed for n ≳ 8 meson systems on the L/a = 48 lattice volume in Fig. <ref>; however, since these tensions are more significant for multi-K^0 than multi-π^+ systems they are unlikely to be signals of nonperturbative QED effects and might result from correlations between LQCD+QED_L results with different n not accounted for in the fitting procedure employed here.












The three-pion scattering amplitude was calculated at LO in χPT in Ref. <cit.> and is given by ℳ_ 3,th = 108 m_K^0^2/f_K^0^4 (with conventions for the LO Lagrangian such that f_π^+∼ 130 MeV).[The SU(2) χPT result for ℳ_ 3,th given in Ref. <cit.> is valid for tree-level K^0K^0 scattering after reinterpreting the SU(2) isospin χPT Lagrangian in terms of the SU(2) V-spin doublet (π^+, K^0). More formally, the V-spin analog of G-parity acts as (π^+,K^0) → (K^0, -π^0) and (K^0, -π^0) → (-π^+, K^0) and therefore relates the three-body contact operators of the SU(2) isospin and SU(2) V-spin Lagrangians by (π^+ π^-)^3 → (K^0 K^0)^3. For the SU(3) flavor symmetric quark mass scheme used here, V-spin is an exact symmetry of the leading order χPT Lagrangian broken only by QED corrections at higher orders.]
This can be combined with Eq. (<ref>) to provide a χPT prediction for the non-relativistic 3K^0 contact interaction,


    η_3^K^0K^0K^0(μ)    = 3/4 m_K^0 f_K^0^4 + m_K^0/64 π^3 f_K^0^8(3√(3) - 4π)ln(m_K^0/2πμ) 
        + 3 m_K^0^3/1024 π^6 f_K^0^8𝒮_ 3,th + 3m_K^0/2 f_K^0^6𝒞_3 ,

where 𝒮_ 3,th and 𝒞_3 are constants defined below Eq. (<ref>), and to obtain a result entirely in terms of m_K^0 and f_K^0, the LO χPT relations <cit.>


    a^K^0K^0 m_K^0 = m_K^0^2/8 π f_K^0^2,
             r^K^0K^0 a^K^0K^0 m^2_K^0 = 3  ,


have been used.
The first of these relations also allows f_K^0 and therefore η_3^K^0K^0K^0(μ) to be predicted numerically at LO in χPT for the quark masses used in this work.
Inserting the LQCD+QED_L results for m_K^0 a^K^0K^0 from Eq. (<ref>) into the first relation in Eq. (<ref>) provides a prediction for f_K^0 at the parameters of this LQCD+QED_L calculation that is valid at LO in χPT:


    af_K^0 = 0.0476(35),  f_K^0 = 139(10)(4) MeV.


Here, the first uncertainty is statistical and the second uncertainty is from the uncertainty in the lattice spacing.
In this calculation m_K^0 = 404(1)(12) MeV is between the physical pion and kaon masses; this can be compared with f_π^+∼ 130 MeV and f_K^0∼ 156 MeV extracted from experiments <cit.>.
Inserting this result for f_K^0, and the result for a^K^0K^0 in Eq. (<ref>), into Eq. (<ref>) then gives the numerical result m_K^0 f_K^0^4 η_3^K^0K^0K^0(m_K^0) = 0.62(27).
This result is valid at LO in χPT and can be compared to LQCD+QED_L results for η_3^K^0K^0K^0(m_K^0) combined with the result for f_K^0 given in Eq. (<ref>).








Dimensionless LQCD+QED_L results for the three-body coupling that are expected to be 𝒪(1) in χPT are given by


    η_3^K^0K^0K^0(m_K^0) m_K^0 f_K^0^4 = 2.28(86), η_3^π^+ π^+ π^+(m_K^0) m_K^0 f_K^0^4 = 0.72(73).


The result of this work for η_3^K^0K^0K^0(m_K^0)  is consistent within 1σ uncertainties with the LQCD results of Ref. <cit.>, which were obtained in a calculation using quark masses similar to those used in this work, corresponding to a pion mass of 352 MeV, and extracted η_3 by fitting to the same 𝒪(L^-6) threshold expansion as used here for multi-K^0 systems.
The corresponding result for η_3^π^+ π^+ π^+(m_K^0) is about 2σ smaller than η_3^K^0K^0K^0(m_K^0), although it is also consistent within 1σ uncertainties with the LQCD results of Ref. <cit.>.
Results for η_3^K^0K^0K^0(m_K^0) and η_3^π^+ π^+ π^+(m_K^0) as well as comparisons to LO χPT and to the LQCD results of Ref. <cit.> are shown in Fig. <ref>.
Instead of using a global fit to extract the scattering lengths and three-body interaction parameters as above, one can instead fix the scattering lengths to the results obtained using fits to Δ E_π^+π^+, Δ E_K^0K^0 shown in Eq. (<ref>) and extract the three-body interaction parameters from fits to multi-meson results on each volume separately in order to provide an estimate of higher-order FV corrections. Using this alternative fitting procedure to extract the three-body interaction parameters from fits using only the L/a = 32 results gives   η_3^K^0K^0K^0(m_K^0) m_K^0 f_K^0^4 = 2.1(9) and  η_3^π^+ π^+ π^+(m_K^0) m_K^0 f_K^0^4 = 0.5(9), while using only the L/a = 48 results gives η_3^K^0K^0K^0(m_K^0) m_K^0 f_K^0^4 = 7.0(2.6) and  η_3^π^+ π^+ π^+(m_K^0) m_K^0 f_K^0^4 = 5.6(2.4), indicating that the L/a=32 results are primarily responsible for constraining the three-body interaction parameters and suggesting that higher-order FV corrections to the extracted three-body interaction parameters may be significant. Additional consistency checks on the three-body interaction determinations are discussed in Ref. [53], where combinations of n-meson FV energy shifts with fixed n are presented that isolate the three-body interaction parameter in the threshold expansion. In all cases, differences between global fit results and results for three-body interactions parameters obtained with fixed n and fixed L/a in this way are smaller than the corresponding differences between global fits results and the combined fits to all n and fixed L/a discussed above.

Differences between  η_3^π^+ π^+ π^+(m_K^0) and η_3^K^0K^0K^0(m_K^0) might arise from QED effects on ℳ_ 3,th beyond LO in χPT, or from QED effects on the matching between ℳ_ 3,th^π^+π^+π+ and η_3^π^+ π^+ π^+(m_K^0).
Differences between the extracted η_3^π^+ π^+ π^+(m_K^0) and η_3^K^0K^0K^0(m_K^0) might also spuriously arise from mismodeling of 𝒪((ML)^-1) relativistic effects estimated in Sec. <ref> to modify FV energy shifts by ∼ 3%.
This estimated shift is larger or comparable to the statistical uncertainties on the three-body energy shift for all nπ^+ systems on the L/a = 48 volume and nπ^+ systems with n≲ 6 on the L/a = 32 volume.
Given these systematic uncertainties in conjunction with the statistical uncertainties on η_3^π^+ π^+ π^+(m_K^0) and η_3^K^0K^0K^0(m_K^0), the results of this work do not provide significant evidence for differences in non-relativistic short-range three-meson interactions arising from QED effects.






§ CONCLUSIONS


In this work, lattice QCD+QED_L has been used to study systems of up to 12 charged or neutral mesons as well as systems of one, two, and three nucleons. Calculations were performed in two lattice volumes with charge-dependent quark masses tuned such that strong isospin breaking effects are negligible and energy differences between charged and neutral systems are primarily QED effects.
While the ground-state energies of two- and three-nucleon systems are determined with few-percent-level precision, QED effects leading to differences between two-nucleon and three-nucleon FV energies are not resolved.
Significantly higher precision will be needed in future calculations of QED effects in multi-nucleon systems.
Differences between charged and neutral FV energy shifts are resolved at the level of 1-2σ for systems of 3-12 mesons, demonstrating the presence of QED effects on meson-meson interactions.
Analysis of the FV energy levels for multi-meson systems using non-relativistic EFT has allowed the extraction of the π^+π^+ and K^0K^0 scattering lengths as well as the 3 π^+ and 3 K^0 interaction parameters.
Differences between a_C^π^+π^+ and a^K^0 K^0 are clearly resolved,
demonstrating that additional QED effects on meson-meson interactions can be resolved beyond the Coulomb photon exchange explicitly included in the EFT.
Differences between the three-body interactions for charged and neutral mesons are not well resolved.

The QED effects on multi-meson systems determined from LQCD+QED_L in this work are well described by NRQED_L results that incorporate short-range two- and three-body contact interactions as well as perturbative Coulomb photon exchange.
Although Coulomb photon exchange must be treated nonperturbatively in sufficiently large volumes, the expansion parameter describing the size of FV Coulomb effects is found to be α / v = α M L/(4π) by examining the convergence pattern of the NRQED_L expansion.
This includes a numerically significant factor of 1/(4π) compared to the parameter α M L discussed in Ref. <cit.>.
For systems with unphysically large α and quark masses such as those studied here, α M L / (4π) ≪ 1 is satisfied for volumes satisfying L ≪ 20 fm for nucleons (L ≪ 50 fm for pions), and Coulomb corrections to the LO strong interaction FV energy shift appear perturbative for L ≲ 6 fm for nucleons (L ≲ 15 fm for pions).
For calculations with physical α and quark masses, FV Coulomb effects are reduced by a factor of 20 for nucleons (40 for pions) and are expected to be perturbative for all practically accessible lattice volumes.
The EFT analysis of this work also demonstrates that NRQED_L results for FV energy shifts are unaffected by the complications of photon zero-mode subtraction up to effects suppressed by 𝒪( (ML)^-3) that are consistently neglected along with other relativistic effects.
Future LQCD+QED_L calculations, especially those using lighter quark masses and/or smaller values of the QED fine structure constant, can therefore be interpreted using hadronic EFTs with perturbative Coulomb effects with a similar procedure to the one undertaken here.
Such calculations will give insight into the quark mass dependence of QED effects on meson-meson interactions, and, combined with higher-precision calculations of multi-nucleon FV energy levels, will permit first principles predictions of QED effects on nucleon-nucleon interactions and QED effects in light nuclei.




§ ACKNOWLEDGEMENTS



We thank the other members of the NPLQCD collaboration, in particular Z. Davoudi and M.J. Savage, for discussions during the initial stages of this work. 
Calculations in this project were performed using the Hyak High Performance Computing and Data Ecosystem at the University of Washington (https://itconnect.uw.edu/research/hpc/), supported, in part, by the U.S. National Science Foundation Major Research Instrumentation Award, Grant Number 0922770, and on clusters at MIT with support from the NEC Corporation Fund.  
The numerical configuration generation (using the BQCD lattice  QCD program <cit.>) was carried out on the IBM BlueGene/Q and HP Tesseract using DIRAC 2 resources  (EPCC, Edinburgh, UK), the IBM BlueGene/Q (NIC, Jülich, Germany)  and the Cray XC40 at HLRN (The North-German Supercomputer  Alliance), the NCI National Facility in Canberra, Australia  (supported by the Australian Commonwealth Government) and Phoenix (University of Adelaide).
The Chroma software library <cit.> was used in the data analysis.
SRB is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under grant Contract Number DE-FG02-97ER-41014.
WD, PES, and MLW are supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under grant Contract Number DE-SC0011090.  
WD is also supported within the framework of the TMD Topical Collaboration of the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, and  by the SciDAC4 award DE-SC0018121.
RH is supported by STFC through grant ST/P000630/1.
MI is supported by the Universitat de Barcelona through the scholarship APIF, by the Spanish Ministerio de Economia y Competitividad (MINECO) under Project No. MDM-2014-0369 of ICCUB and with additional European FEDER funds, under the contract FIS2017-87534-P.
HP is supported by DFG Grant No. PE 2792/2-1.
PELR is supported in part by the STFC under contract ST/G00062X/1.
GS is supported by DFG Grant No. SCHI 179/8-1.
PES is supported in part by the National Science Foundation under CAREER Award 1841699.
MLW was supported in part by an MIT Pappalardo Fellowship.
RDY and JMZ are supported by the Australian Research Council Grants FT120100821, FT100100005, DP140103067 and DP190100297.  
This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.









§ MATCHING NRQED_L AND QED_L


Neglecting relativistic effects suppressed by M^-1 (including nonlocal effects from zero-mode subtraction discussed in Refs. <cit.>), the NRQED_L Lagrangian is identical for bosons and fermions. For a particle with charge Q=1, it is given by


    ℒ^NRQED_L   = ψ^†( i D_0 - D_i D^i/2M) ψ -  2π a/M (ψ^†ψ)^2  + ℒ_γ^ξ
       = ψ^†( i ∂_0 - ∂_i ∂^i/2M) ψ - e A_0 (ψ^†ψ) -  2π a/M (ψ^†ψ)^2  + ℒ_γ^ξ,


where in this section we work in Minkowski spacetime with (-+++) signature, A_μ = A_μ^† is the photon field, F_μν = ∂_μ A_μ - ∂_ν A_μ is the field strength tensor, D_μ = ∂_μ + i A_μ, and in generic R_ξ gauge the photon Lagrangian is


    ℒ_γ^ξ   =  -1/4F_μνF^μν + 1/2ξ(∂_μ A^μ)^2.


Results for the Landau gauge QCD+QED_L calculations performed in the main text are obtained by setting ξ=0.
In what follows, a finite spatial volume of extent L^3 with PBCs is considered.
Zero mode subtraction can be implemented in NRQED_L by defining the FV photon field as a Fourier transform of the zero-mode subtracted momentum-space field,


    A_μ(x) = ∫dp^0/2π1/L^3∑_ň∈ℤ^3 ∖{0̌} e^-ip^0 x^0 + 2π i/Lň·x̌Ã_μ(p^0,ň), Ã_μ(p^0,ň) = ∫ d^4x e^ip^0 x^0 - 2π i/Lň·x̌ A_μ(x),


where ň∈ℤ^3 ∖{0̌} excludes the photon zero mode.
Zero-mode subtraction can be defined at the path integral level as a constraint on the photon field <cit.>, but for perturbative matching with QED_L Eq. (<ref>)-(<ref>) define NRQED_L.

The NRQED_L massive particle propagator G^NRQED_L is given by


    G^NRQED_L(p^0, ň)    = ∫ d^4x e^ip^0x^0 - 2π i/Lň·x̌< ψ(x) ψ^†(0) >  = i/p^0 - 1/2M(2πň/L)^2 + i ϵ.


The photon propagator is given by


    G_μν^γ_L(p^0, ň)    = ∫ dx^0 ∑_x̌ e^ip^0x^0 - 2π i/Lp̌·x̌< A_μ(x) A_ν(0) > = i [ g_μν - (1-ξ) p_μ p_ν/p^2]/(p^0)^2 - (2πň/L)^2 + i ϵ.


Introducing Fourier transformed fields,


    ψ̃(p^0,ň) = ∫ d^4x e^ip^0x^0 - 2π i/Lň·x̌ψ(x), ψ(x) = ∫dp^0/2π1/L^3∑_ň∈ℤ^3 e^-ip^0x^0 + 2π i/Lň·x̌ψ̃(p^0,ň).


Matching between QED_L and NRQED_L is performed for four-point correlation functions describing (off-shell) particles with energy M and three-momentum p̌≡2π/Lř with ř∈ℤ_3,


    ℳ^NRQED_L   = G^NRQED_L(0,ř)^-4< ψ̃(0,ř)^2 ψ̃(0,ř)^† 2> .


At tree level this correlation function is given by


    ℳ^NRQED_L_ LO   = -8π a/M .


The choice ř≠ 0, which does not affect the tree-level correlation function, is made in order to regulate IR divergences in the one-loop correlation function.
Although one-photon-exchange contributions lead to an IR divergence (regulated for instance by considering non-zero momentum transfer) in the NRQED analog of Eq. (<ref>), one-photon-exchange contributions to the NRQED_L amplitude in Eq. (<ref>) vanish because of zero-mode subtraction.

Matching is performed by expanding ℳ^NRQED_L and its QED_L analog ℳ^QED_L perturbatively in the small parameters α, 1/(ML), and a/L and defining higher-order terms in the NRQED_L Lagrangian so that QED_L and NRQED_L agree order by order in this expansion.
Matching will be performed to leading order in 1/(ML) and third order in η_L = α M L / (4π) and a/L.
This corresponds to NNLO in the power counting of Eq. (<ref>) and is equivalent to the non-relativistic limit with up to two-loop contact interactions and Coulomb photon exchange.
This matching is also sufficient for N^3LO in the power counting of Eq. (<ref>), which decreases the relative importance of Coulomb photon exchange and only requires one-loop Coulomb photon exchange effects.


This LO NRQED_L amplitude in Eq. (<ref>) can be straightforwardly matched to its QED_L analog.
The scalar QED_L Lagrangian is


    ℒ^QED_L   = -(D_μφ)^† D^μφ - M^2 φ^†φ  - 8π a M (φ^†φ)^2  + ℒ_γ^ξ
       = -(∂_μφ)^†∂^μφ - M^2 φ^†φ + i e A^μ ( φ^†∂_μφ - (∂_μφ)^†φ) 
        - e^2 A_μ A^μφ^†φ    -  8π a M (φ^†φ)^2  + ℒ_γ^ξ,


where ℒ_γ^ξ is defined in Eq. (<ref>).
The fields appearing in the QED_L and NRQED_L Lagrangians are related by


    ψ(x) = √(2 M) e^i M tφ(x).


The scalar QED_L particle propagator is


    G^QED_L(p^0, ň)    = ∫ dx^0 ∑_x̌ e^i p^μ x_μ< φ(x) φ^†(0) > = i/(p^0)^2 - ( 2πň/L)^2 - M^2 + i ϵ.


The QED_L photon propagator is identical to the NRQED_L propagator Eq. (<ref>).
Introducing Fourier transformed fields φ̃ as in Eq. (<ref>), the FV two-particle amplitude in QED_L normalized identically to the NRQED_L amplitude is given by


    ℳ^QED_L_ LO   = 1/4M^2 G^QED_L(M,ř)^-4< φ̃(M,ř)^2 φ̃(M,ř)^† 2> 
       = -8π a/M ,


in agreement with Eq. (<ref>). This demonstrates that the four-fermion contact interactions in Eq. (<ref>) and Eq. (<ref>) are normalized consistently at tree level.

Corrections to the LO amplitude arise from one-loop diagrams shown in Figs. <ref>-<ref>.
The one-loop diagrams shown in Fig. <ref> involve similar sums/integrals over loop momenta and differ only in the number of photon propagators n_γ and contact interactions present.
The contribution to the amplitude from each diagram is denoted ℳ^NRQED_L_ NLO(n_γ) where n_γ∈{ 0, 1, 2} labels the numbers of photon propagators present,


    ℳ^NRQED_L_ NLO(n_γ) =  i 16 π^2 a^2 /M^2 L^3( 2 α M/a)^n_γ∑_ň^(')∫dk^0/2π( 1 / k^0 - 2π^2(ň - ř)^2/ML^2 + i ϵ)  
       ×( 1/ k^0 + 2π^2(ň + ř)^2/ML^2 - i ϵ) (1 - (1-ξ)(k^0)^2/(k^0)^2 - 4π^2ň^2/L^2 + iϵ/ (k^0)^2 - 4π^2ň^2/L^2 + iϵ)^n_γ
       = -16 π^2 a^2 /M^2 L^3( α M/a)^n_γ∑_ň^(')[ Res(- 2π^2(ň - ř)^2/ML^2 + iϵ) +  Res(-2π|ň|/L + iϵ) ],


where  Res(x) indicates the residue of the integrand in the first line at the pole x and ∑_ň^(') corresponds to ∑_ň∈ℤ^3 for n_γ = 0 and to ∑_ň∈ℤ^3∖{0̌} for n_γ≥ 1.
After taking the ϵ→ 0 limit, the residue at the particle pole 2π^2(ň+ř)^2/(ML^2) is given by 


    Res(-2π^2(ň-ř)^2/ML^2)    = -ML^2/4π^2(ň^2 + p̌^2)[  1/4π^4 (ň - ř)^4/M^2 L^4  - 4π^2ň^2/L^2 - (1 - ξ)4π^4(ň-ř)^4/M^2L^4/( 4π^4 (ň - ř)^4/M^2 L^4  - 4π^2 ň^2/L^2)^2]^n_γ
       = - ML^2/4π^2( - L^2/4π^2)^n_γ1/(ň^2 + ř^2)ň^2n_γ[ 1 + 𝒪( (ML)^-2) ],


where the last line includes an expansion in powers of (ML)^-1.
This expansion is legitimate provided that the sum over ň converges, which holds for n_γ≥ 1.
For n_γ = 0 there is a linear UV divergence that can be removed by adding a UV counterterm,


    ∑_ň∈ℤ^3 1/ň^2 + ř^2→lim_Λ→∞∑_ň∈ℤ^3 ^|ň| < Λ1/ň^2 + ř^2- 4πΛ,


and after making the replacement of Eq. (<ref>) the 1/(ML) expansion can be performed in Eq. (<ref>).
The photon pole residue at -2π|ň|/L involves energy denominators of order 1/L and is suppressed by 𝒪( (ML)^-1) compared to the particle-pole contributions in Eq. (<ref>) where particle propagator energy denominators are of order M/L^2.
For all diagrams in Figs. <ref>-<ref>, particle poles only appear in either the upper or lower half of the complex plane and only photon poles contribute.
After performing the energy integrals these diagrams can be straightforwardly verified to be suppressed by 𝒪( (ML)^-1) or 𝒪( (ML)^-2) compared to Eq. (<ref>).
The full NLO amplitude in NRQED_L is therefore a sum of the amplitudes ℳ^NRQED_L_ NLO(n_γ) associated with the diagrams in Fig. <ref> obtained by substituting Eq. (<ref>) into Eq. (<ref>),


    ℳ^NRQED_L_ NLO(n_γ)    = 4π a/M( - α M L/2π^3)^n_γ( a/π L)^1-n_γ∑_ň^(')1/(ň^2 + ř^2)ň^2n_γ[ 1 + 𝒪( (ML)^-1) ].


Contributions with different n_γ differ parametrically and can be matched independently between NRQED_L and QED_L.


The QED_L amplitudes associated with the NLO diagrams in Fig. <ref> are given by


    ℳ^QED_L_ NLO(n_γ) = -i 64 π^2 a^2/L^3( -α/2aM)^n_γ∫dk^0/2π∑_ň^(')1/(k^0 + M)^2 - 4π^2(ř+ň)^2/L^2 - M^2 + iϵ
       ×1/ (k^0 - M)^2 - 4π^2(ř-ň)^2/L^2 - M^2 + iϵ(  N_γ^ξ(k^0, ň)  /(k^0)^2 - 4π^2ň^2/L^2 + i ϵ)^n_γ,


where


    N_γ^ξ(k^0,ň)    = (k^0)^2 - 4M^2 + ( 4π^2/L^2) (4ř^2 - ň)^2 
        - (1-ξ)[ -(k^0)^4 + 4M^2 (k^0)^2 + ( 4π^2/L^2) (4ň·ř - ň^4)/(k^0)^2 - (4π^2/L^2) ň^2 + iϵ].


The k^0 integral can be performed with contour integration.
Closing the contour in the upper half plane, the result includes contributions from a particle pole k^0 = M - √(M^2 + 4π^2(ň+ř)^2/L^2) = -2π^2(ň+ř)^2/(ML^2)[1 + 𝒪((ML)^-2)], a photon pole at k^0 = 2π|ň|/L, and an antiparticle pole at k^0 = -2M[1 + 𝒪((ML)^-2)],


    ℳ^QED_L_ NLO(n_γ) = 64 π^2 a^2/L^3( -α/4aM)^n_γ∑_ň^(')[ Res(M - √(M^2 + 4π^2(ň+ř)^2/L^2 - iϵ)) . 
       . + Res(-2π|ň|/L + iϵ) + Res(-M - √(M^2 + 4π^2/L^2(ň-ř)^2 + iϵ)) ].


Taking the ϵ→ 0 limit and expanding to leading order in (ML)^-1, the residue at the particle pole is given by


    Res(M - √(M^2 + 4π^2(ň+ř)^2/L^2 - iϵ))    = L^2/16π^2 M( M^2L^2/π^2)^n_γ1/(ň^2+ř^2)ň^2n_γ.


As in the NRQED_L case, the residue at the photon pole -2π|ň|/L is suppressed compared to the residue at the particle pole by 𝒪(ML)^-1.

The residue at the antiparticle pole at -2M(1+𝒪(ML)^-1) does not appear in the corresponding NRQED_L expression Eq. (<ref>) and is therefore associated with contributions that do not appear in loop diagrams in NRQED_L.
This residue is given by


    Res(-M - √(M^2 + 4π^2/L^2 + iϵ)) = -1/16M^3( - π^2  /M^2 L^2)^n_γ
       ×[ 2ř^2(-4 + ξ)  + ň^2(-3+2ξ) - 4ř·ň(-2+ξ)  ]^n_γ.


The term involving two contact interactions involves the UV divergent sum ∑_ň^(') 1, which for n_γ = 0 must be consistent with the infinite-volume result ∫ d^Dk   1 = 0 after subtracting UV counterterms,


    ∑_ň∈ℤ^3 1 → 0.


This leads to a vanishing contribution from the n_γ = 0 diagram appearing in the absence of QED_L, which is consistent with the expectation that antiparticle-pole contributions from off-shell intermediate states that are absent non-relativistically do not lead to power-law FV effects in local field theories.
The terms with n_γ≥ 1 differ by zero mode subtraction.
Since the ň = 0 contribution to ∑_ň 1 is unity, it follows from Eq. (<ref>) that


    ∑_ň∈ℤ^3 ∖{0̌}  1 → -1.


Sums with ň^2k with k > 0 similarly require UV counterterms and vanish after including them.
Zero-mode contributions to these sums vanish, and so


    ∑_ň∈ℤ^3 ∖{0̌}ň^2k→ 0,


for k > 0. After subtracting UV counterterms the antiparticle pole contribution becomes


    ∑_ň∈ℤ^3 ∖{0̌} Res(-M - √(M^2 + 4π^2/L^2 + iϵ)) →1/16M^3( - π^2 2ř^2(-4 + ξ) /M^2 L^2)^n_γ ,


for n_γ≥ 1.
This leads to a finite contribution to ℳ^QED_L_ NLO that does not arise in the corresponding loop diagrams contributing to ℳ^NRQED_L_ NLO.
Local counterterms involving powers of L^-1 must be added to the NRQED_L Lagrangian in order to reproduce this relativistic effect arising from zero-mode subtraction.
Inserting Eq. (<ref>) into Eq. (<ref>) shows that the dominant contribution to ℳ^QED_L_ NLO arising from the n_γ = 1 diagram antiparticle pole is


    -4π a/M(απ^3 (-4 + ξ) ř^2/4 M^5 L^5) = -8π a/M(π (-4 + ξ) /16 )  α/(M L)^3(p̌^2/M^2).


This contribution vanishes for FV systems in the two-particle rest frame where p̌ = 0.
For boosted systems with a non-zero center-of-mass velocity (assumed to be non-relativistic p̌^2 ≪ M^2 in the (ML)^-1 expansion above), this contribution is proportional to the LO contact interaction times the velocity squared times a relativistic QED suppression factor of α / (ML)^-3.
For boosted systems, the NRQED_L contact interaction in Eq. (<ref>) must therefore be supplemented with a nonlocal counterterm suppressed by the same factor of α / (ML)^-3.
This is analogous to the self-energy of a single scalar field, which is shown in Ref. <cit.> to require a nonlocal counterterm equal to the scalar mass times the velocity squared times a relativistic QED suppression factor of α / (ML)^-3.
The LQCD+QED_L calculations discussed in the main text are performed in the center-of-mass rest frame, and these nonlocal counterterms vanish.
Non-vanishing nonlocal counterterms appear for fermion masses in NRQED_L even in the two-particle rest frame, but since effects proportional to (ML)^-1 are neglected throughout this work it is consistent to neglect all nonlocal counterterms suppressed by α / (ML)^-3.

The NLO QED_L amplitude is therefore given by inserting the particle-pole residue in Eq (<ref>) into Eq. (<ref>),


    ℳ^QED_L_ NLO(n_γ)    = 4π a/M( -α M L/2π^3)^n_γ( a/π L)^1-n_γ∑_ň^'1/(ň^2 + ř^2) ň^2n_γ[  1 + 𝒪( (ML)^-1) ],


which is identical to Eq. (<ref>).
The (ML)^-1 suppression of diagrams in Figs. <ref> and <ref> arising from the absence of particle pole contributions is identical in QED_L and NRQED_L.
Additional diagrams appearing in QED_L but not NRQED_L associated with the two-particle-two-photon vertex or particle-antiparticle pair creation can be similarly verified to be suppressed by powers of (ML)^-1.
The NRQED_L Lagrangian in Eq. (<ref>) therefore reproduces QED_L at NLO.

Matching at NNLO proceeds similarly.
The NNLO diagrams shown in Fig. <ref> can all be expressed in terms of the amplitude 


    ℳ^NRQED_L_ NNLO(n_1,n_2,n_3) 
       = 64π^3 a^3/M^3 L^6( α M/a)^n_1+n_2+n_3∑_ň,m̌^(')∫dk^0/2π1/(k^0 + 2π^2(ř-ň)^2/ML^2 - iϵ) (k^0 - 2π^2(ř+ň)^2/ML^2 + iϵ) 
       ×(1 - (1-ξ)(k^0)^2/(k^0)^2 - 4π^2ň^2/L^2/ (k^0)^2 - 4π^2ň^2/L^2 + iϵ)^n_1∫dq^0/2π1/(q^0 + 2π^2(ř-m̌)^2/ML^2 - iϵ) (q^0 - 2π^2(ř+m̌)^2/ML^2 + iϵ)  
       ×(1 - (1-ξ)(k^0 - q^0)^2/(k^0 - q^0)^2 - 4π^2(ň-m̌)^2/L^2/ (k^0 - q^0)^2 - 4π^2(ň - m̌)^2/L^2 + iϵ)^n_2(  1 - (1-ξ)(q^0)^2/(q^0)^2 - 4π^2m̌^2/L^2/(q^0)^2 - 4π^2m̌^2/L^2 + iϵ)^n_3,


where n_i ∈{ 0, } for i∈{1,2,3} labels whether each interaction is a four-particle contact interaction or photon exchange and ∑_ň,m̌^(') excludes ň=0 if n_1 = 1 or n_2 =1 and excludes m̌=0 if n_2 = 1 or n_3 = 1.
Both energy integrands include poles where a particle is on-shell as well as poles where a photon is on-shell, and, as above, photon-pole contributions are suppressed by 𝒪( (ML)^-1).
Evaluating the q^0 and k^0 integrals then gives


    ℳ^NRQED_L_ NNLO(n_1,n_2,n_3)    = - 4π a/M( a/π L)^2 - n_1 - n_2 - n_3( - α M L/2π^3)^n_1+n_2+n_3
       ×∑_ň,m̌^(')1/(ň^2 + ř^2) (m̌^2+ř^2) ň^2n_1 (ň - m̌)^2n_2m̌^2n_3[ 1 + 𝒪((ML)^-1) ],


where UV counterterms should be included as in Eq. (<ref>).
Similarly to the NLO case, the NNLO diagrams in Figs. <ref> and <ref> only include photon pole contributions and are suppressed by 𝒪( (ML)^-1) or 𝒪( (ML)^-2) compared to Eq. (<ref>).

The corresponding QED_L NNLO amplitudes are given by


    ℳ^QED_L_ NNLO(n_1,n_2,n_3) = 1024 π^3 a^3 M/L^6( -α/4 a M)^n_1+n_2+n_3∑_ň,m̌^(')∫dk^0/2π(  N_γ^ξ(k^0, ň)  /(k^0)^2 - 4π^2ň^2/L^2 + i ϵ)^n_1
       ×1/( (k^0 - M)^2 -  4π^2(ň-ř)^2/L^2 - M^2 + iϵ) ((k^0 + M)^2 - 4π^2(ň+ř)^2/L^2 - M^2 + iϵ) 
       ×∫dq^0/2π(  N_γ^ξ(k^0-q^0, ň-m̌) /(k^0-q^0)^2 - 4π^2(ň-m̌)^2/L^2 + i ϵ)^n_2(  N_γ^ξ(q^0, m̌)  /(q^0)^2 - 4π^2m̌^2/L^2 + i ϵ)^n_3
       ×1/( (q^0 - M)^2 - 4π^2(m̌-ř)^2/L^2 - M^2 - iϵ)   ( (q^0 + M)^2 - 4π^2(m̌+ř)^2/L^2 - M^2 + iϵ)  .


Contributions from photon and antiparticle poles are again suppressed by powers of (ML)^-1, and evaluating the q^0 and k^0 energy integrals gives 


    ℳ^QED_L_ NNLO(n_1,n_2,n_3)    =  - 4π a /M( a/π L)^2 - n_1 - n_2 - n_3( -α ML/2π^3)^n_1 + n_2 + n_3
       ×∑_ň,m̌^(')1/(ň^2 + p̌^2) (m̌^2 + p̌^2)ň^2n_1(ň-m̌)^2n_2m̌^2n_3[ 1 + 𝒪( (ML)^-1) ].


All other diagrams are again suppressed by powers of (ML)^-1, and agreement between Eq. (<ref>) and Eq. (<ref>) shows that nonlocal counterterms are suppressed by powers of (ML)^-1 and can be neglected to the accuracy considered here.

While the matching in this section has been explicitly performed for scalar QED_L, the (ML)^-1 suppression of photon and antiparticle poles only relies on the structure of QED_L propagator denominators that are identical for scalars and fermions.
The (ML)^-3 suppression of antiparticle pole contributions leading to nonlocal two-body counterterms also arises from the denominator structure of the QED_L propagators and is expected to be generic for bosons and fermions.
The numerator structure of the scalar QED_L propagators above is relevant for (ML)^-1 suppressed relativistic effects, and in particular the scalar QED_L antiparticle pole in Eq. (<ref>) includes vanishing numerator factors for a system at rest that lead to p̌^2/M^2 velocity suppression of nonlocal counterterms in scalar NRQED_L.
This cancellation might be absent for fermions, and nonlocal two-body counterterms might be relevant for QED_L calculations of charged fermions in the center-of-mass rest frame.
This would parallel the situation for one-body nonlocal counterterms, which arise at 𝒪(α/(ML)^3) for fermions with any center-of-mass velocity and for scalars with non-zero center-of-mass velocity.




§ FITTING PROCEDURES




 §.§ Energy level determination


In general the spectral representation in Eq. (<ref>) cannot be inverted to determine the full energy spectrum from finite samples of correlation functions over a finite range of source/sink separations.
Any fitting procedure to extract energies from correlation functions involves making several choices, in particular the range of t to include in the fit, the number of excited states to include in a truncation of Eq. (<ref>) to use as a fitting model, and how to estimate the covariance matrix from a finite statistical ensemble.
In order to assess the systematic uncertainties associated with these fitting choices and provide a reproducible procedure for extracting energy levels from correlation function results, we use an approach detailed here for making fitting choices based on well-defined statistical criteria and random sampling over the space of possible fitting choices.

The first step in this fitting procedure is choosing the maximum source/sink separation t_ max included in the fit.
For (multi-)baryon correlation functions, the signal-to-noise (StN) problem implies that results with larger temporal separation t make exponentially smaller contributions to χ^2 when fitting energy levels from correlation functions.
For the nucleon, the scaling StN(G) = G(t)/√(Var(G(t)))∼ e^-(M_N - 3/2m_π)t predicted by Parisi <cit.> and Lepage <cit.> applies in the limit of a large statistical ensemble size N →∞ and shows that fit results are exponentially insensitive to the choice of t_ max.
Similar results apply for multi-nucleon systems with baryon number A where StN(G) = G(t)/√(Var(G(t)))∼ e^-A(M_N - 3/2m_π)t for large t.
For fixed N (and, since volume-averaging increases StN, for fixed L) it is important to choose a fixed t_ max before the StN has degraded to the point where correlation function estimates are unreliable.[The complex phases of baryon correlation functions are circular random variables, and there is therefore a noise region at large t where ln N is not much larger than the variance of the phase distribution and the sample mean is a systematically unreliable estimator of the average correlation function <cit.>.]

The StN ratio decreases with t at small and intermediate t before saturating at an 𝒪(1) value in the noise region, and to avoid the noise region t_ max should not exceed the smallest t where StN(G) = G(t) / √(VarG(t)) reaches a specified 𝒪(1) cutoff.
For positive-definite meson correlation functions this issue does not arise, and the maximum t that can be reliably included in the fit is only limited by the accuracy by which finite-temperature effects are modeled.
In this work, finite-temperature excited-state effects are neglected, and finite-temperature effects on correlation functions with t∼ T/2, where T is the length of the Euclidean time direction, are not reliably modeled.

In our fitting procedure, t_ max is therefore chosen to be the minimum t satisfying either: the correlation function noise-to-signal ratio at (t+a) is smaller than a specified tolerance (final results use tol_ noise = 1.0), the correlation function sample mean at t + a is negative,[If the average correlation function is not expected to be positive definite, then this condition should not be enforced. In principle, the source and sink interpolating operators used in this work differ from one another, and it is possible for the sign of an average correlation function to fluctuate at small t where excited-state contributions with opposite sign to the ground-state contribution can be significant. In practice, small t fluctuations of the sign of the sample mean correlation function that could be attributed to excited-state effects are not observed in this work, and a negative sample mean correlation function is taken as an indicator of large statistical noise.] or t+a is larger than a finite-temperature cutoff (final results use tol_ temp = 3 T /8).
The values of tol_ noise and tol_ temp are free parameters in our fitting procedure that must be varied to assess the sensitivity of fit results to these choices; for concreteness the parameter choices are presented here that lead to the final results quoted in the main text. 
Results are found to be relatively insensitive to the parameter choices controlling t_ max, and for example varying the noise tolerance in the range tol_ noise∈ [0.1, 1] leads to results with consistent central values at the 1σ level and few-percent variation of the corresponding uncertainties.




The next step in our fitting procedure is to choose t_ min, the minimum t included in the fit.
The choice of t_ min significantly impacts how well excited-state effects at small t can be resolved and how many excited states should be included in fits.
Furthermore, the uncertainties of energy-level determinations are exponentially sensitive to the choice of t_ min for baryons because of the StN problem.
Fit results are therefore more sensitive to the choice of t_ min than to the choice of t_ max.
Rather than choose a single t_ min, it is preferable to sample from many possible choices of t_ min and quantify the sensitivity to this choice as the systematic error.
The minimum permissible t_ min is fixed by the temporal nonlocality in the lattice action, and for the improved action used in this work  the transfer matrix involves fields on two adjacent timeslices <cit.> and t_ min≥ 2 is required.

The largest allowed t_ min is limited by t_ min≤ t_ max - t_ plateau, where t_ plateau is a free parameter that is not found to significantly affect final results when varied over the range 2 ≲ t_ plateau≲ 8 (final results use t_ plateau = 4).
For each type of interpolating operator included in a combined fit, t_ min is sampled randomly within this range until either all possible values of t_ min have been chosen or a maximum of N_ fits = 200 fits have been performed.

With t_ min and t_ max specified, the covariance matrix 𝒞^ij_tt^' must be estimated for t_ min≤ t,t^'≤ t_ max and interpolating operators i,j ∈{1,… N_ op}.
Fits involving a large number N_ pts of time separations and interpolating operator choices may not satisfy the condition N ≫ N_ pts^2 needed to ensure that the N terms contributing to the  N_ pts× N_ pts sample covariance can accurately estimate the true underlying covariance matrix, where N_ pts is the total number of source/sink separations from all interpolating operators included in the fit.
Shrinkage techniques have been developed to provide more accurate estimates of the underlying covariance matrix than the sample covariance matrix when N ≫ N_ pts^2 is not satisfied <cit.>.
Shrinkage estimators 𝒮_tt^'^ij(λ) of the covariance matrix are constructed as mixtures of a well-conditioned target matrix T_tt^'^ij and the covariance matrix 𝒞_tt^'^ij estimated using standard bootstrap techniques from N_ boot = 200 samples of N correlation functions G_i,a^b(t) with a ∈{1,…, N} and b∈{1,…,N_ boot} drawn from the original correlation function ensemble with replacement,

    𝒮_tt^'^ij(λ) = 𝒞_tt^'^ij(1-λ) + T_tt^'^ijλ,

where 0 ≤λ≤ 1 is a shrinkage parameter.
A common choice of well-conditioned target matrix for many problems in statistics is the identity matrix; however, this does not accurately describe the underlying covariance matrix for correlation functions whose diagonal entries decrease exponentially with t.
Following applications of shrinkage to lattice QCD in Ref. <cit.>, we take T_tt^'^ij = diag(𝒞_tt^'^ij).
In this case, shrinkage corresponds to an interpolation between a fully correlated fit with λ = 0 and an uncorrelated fit with λ = 1.
Shrinkage gives an unbiased estimator of the underlying covariance matrix in the infinite-statistics limit provided λ vanishes sufficiently quickly in this limit.
It can be shown <cit.> that for finite N an optimal λ^* ≠ 0 satisfying this restriction can be chosen in order to minimize the average mean-squared difference between 𝒮_tt^'^ij(λ) and the underlying covariance matrix, and that a sample estimate for λ^* is given by

    λ^*    = Max(0, Min(1, ∑_a=1^N ∑_t,t^',i,j[ G_a^i(t) G_a^j(t^') - S^ij(t,t^') ]^2 /N^2 ∑_t,t^',i,j[ S^ij(t,t^') - δ^ijδ_t,t^']^2) ),

where

    G_a^i(t)    = G_a^i(t) - G^i(t)/√(S^ii(t,t)), 
    S^ij(t,t^')    = S^ij(t,t^')/√(S^ii(t,t)S^jj(t^',t^')),

are defined in terms of the sample mean correlation function and sample covariance as in Ref. <cit.>

    G^i(t)    = 1/N∑_a=1^N G_a^i(t), 
    S^ij(t,t^')    =  1/N-1∑_a=1^N [ G_a^i(t) G_a^j(t^') - G^i(t) G^j(t^') ],

such that shrinkage of S^ij(t,t^') with T_tt^'^ij = diag(𝒮_tt^'^ij) corresponds to shrinkage of  S^ij(t,t^') with the identity matrix as a target, and the results of Ref. <cit.> assuming an identity matrix target can be applied.
The covariance matrix estimate with optimal shrinkage is then given by 𝒮_tt^'^ij(λ^*).
Fits to truncations of Eq. (<ref>) including e excited states can then be performed by minimizing the corresponding χ^2 function defined by

    χ^2 = ∑_t,t^' = t_ min^t_ max∑_i,j = 1^N_ op( G_i^B(t) - f(t,𝐄,𝐙) ) [ 𝒮(λ^*)^-1]_tt^'^ij( G_j^B(t^') - f(t^',𝐄,𝐙) ),

where G_i^B = 2G_i - 1/N_ boot∑_b=1^N_ bootG_i^b includes a 1/N bias correction estimated using bootstrap techniques and 𝐄 and 𝐙 denote the energies and overlap factors appearing in Eq. (<ref>), including excited-state and thermal effects.
Since the overlap factors enter f_N linearly, the values of Ž minimizing χ^2 for fixed Ě can be determined by solving a system of linear equations analogous to variable projection techniques <cit.>.
χ^2-minimization can therefore be efficiently performed by using a nonlinear optimization method to determine Ě with the optimal Ž determined by solving a system of linear equations at each step of nonlinear optimization for Ě.
In order to ensure positivity of the spectrum and remove fitting degeneracies, the parameters used for nonlinear optimization are ln E_0 and ln(E_k - E_k-1) for 1≤ k ≤ e.

For each randomly sampled choice of t_ min, the next step in the fitting procedure is to determine the number of excited states to be included in the sum of exponentials used as a fit function.
This is done by first performing a fit including zero excited states and then adding successively more excited states until the addition of the next excited state does not improve the goodness of fit according to an information criterion.
This work employs the Akaike Information Criterion <cit.> (AIC) with a cutoff chosen to penalize overfitting in which a fit with e excited states is only preferred over a fit with e-1 excited states if AIC(e) - AIC(e-1) < -𝒜 N_ dof(e) where N_ dof(e) = N_ pts - N_ params(e) is the number of degrees of freedom of the fit, N_ params(e) is the number of fit parameters for a fit with e excited states, and AIC(e) = 2N_ params(e) + χ^2(e) + k with χ^2(e) the (unreduced) χ^2 of the fit defined in Eq. (<ref>) and k is an irrelevant e-independent constant.
This choice corresponds to a preference for an e state fit only if it improves the χ^2/N_ dof by an 𝒪(1) value 𝒜 compared to the (e-1) state fit.
For baryon correlation functions a value of 𝒜 = -0.5 is used, while for multi-meson correlation functions a value of 𝒜 = -0.1 is used.[Optimal shrinkage values of λ^* ≳ 0.1 appear for n-meson correlation functions with n ≳ 6 and χ^2 values are correspondingly lower than would be expected for fully correlated χ^2 minimization, leading to smaller absolute changes in AIC for multi-meson correlation functions than for multi-baryon correlation function. Increasing -𝒜 from 0.1 over the range -𝒜∈ [0.1, 2] leads to consistent results with larger uncertainties because precise and accurate two-state fits with small t_ min are rejected in favor of one-state fits more frequently.]

Bootstrap resampling techniques are then used to estimate the uncertainty on the ground-state energy extracted from the fit with the preferred number of excited states in each fit region <cit.>.
The same χ^2-minimization procedure and estimated covariance matrix 𝒮(λ^*) are used to determine the spectrum for each of N_ boot bootstrap resampled ensembles.
Results are found to be insensitive to the choice of N_ boot, and final results use  N_ boot = 200.
The 67% confidence interval for the ground-state energy is then obtained from the quantiles of the distribution of differences between the b'th bootstrap sample result E^b,f_0 and the fit result E_0^f obtained for fit range f,

    δ E_0^f = 1/2[ Q_5/6( E^b,f_0 - E_0^f )  - Q_1/6( E^b,f_0 - E_0^f ) ],

where Q_p(x^f) is the p-th quantile of the set of fit results with elements x^f.
Using this definition δ E_0^f minimizes the impact of outlier bootstrap samples compared to a definition based on the standard deviation of E^b,f_0 - E_0^f <cit.>.
An analogous procedure is used to estimate uncertainties for excited-state energies and overlap factors.



Several additional checks are used to ensure the robustness of χ^2-minimization results: two different optimization algorithms, Nelder-Mead (NM) and conjugate gradient (CG), are used[Newton's method is used in place of Nelder-Mead if N_ states = 1, since Nelder-Mead does not work for a single fit parameter.] and are verified to give energies that differ by less than a specified tolerance[Fits resulting in ground-state energies less than tol_ sol, which appeared only for the 11K^0 system on the L/a=32 lattice volume, were also rejected.]  (final results use tol_ sol = 10^-5), median results for ground- and excited-state energies from the bootstrap samples are verified to agree with fit results from the average correlation functions for each energy level within a specified tolerance (final results use tol_ med = 2σ), uncorrelated fit results obtained by repeating the χ^2-minimization procedure with 𝒮(λ=1) are verified to give consistent results for each energy level within a specified tolerance (final results use tol_ corr = 5σ), the χ^2/N_ dof is verified to be less than a specified tolerance (final results use tol_χ^2=2).
This defines a reproducible and automatable procedure for fitting correlation functions, including sampling of possible fit ranges and excited-state model selection, in which fit results are functions of only the tolerances described above and the given correlation functions.
A graphical illustration of the fitting procedure is shown in Fig. <ref>.
This fitting procedure was implemented in the ulia-2017 using the therwise, fits for multi-nucleon energy shifts are performed identically to fits for multi-meson energy shifts not including thermal effects.




 §.§ Fit results


Figs. <ref>-<ref> show fit results for n π^+ systems with L/a = 32.
Results with L/a = 48 are shown in Figs. <ref>-<ref>.
Figs. <ref>-<ref> show analogous fit results for n K^0 systems with L/a∈{32,48}.
Single-nucleon fit results for p and n are shown in Fig. <ref>.
Two-nucleon fit results are shown for pp and nn in Fig. <ref> and for np systems in Fig. <ref>.
Three-nucleon fit results are shown in Fig. <ref>.







































