Charged multi-hadron systems in lattice QCD+QED

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.


I. 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) [1][2][3]. Coupled lattice QCD+QED calculations have also been performed for meson and baryon masses [4][5][6][7], leptonic decay rates [8][9][10], and the anomalous magnetic moment of the muon [11][12][13][14][15][16][17][18][19], 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 [20].
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 [21]. It is also interesting to consider systems of multiple charged mesons, which allow an investigation of similar QED effects without of 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 [1-3, 5, 6, 8, 9, 19, 22-37]. The current study makes use of the QED L formulation [25], in which the photon spatial zero-mode is removed, as has been recently studied in detail for single-hadron systems in Refs. [1-3, 5, 6, 8, 9, 19, 22-27, 29, 31-37]. Subtleties related to the nonlocality of zero-mode subtraction in QED L have been understood for single-hadron systems [2,29,31,35,37,38], 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 [39,40] and generalized in recent years (see Ref. [41] 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. [42] argues that Coulomb effects can be treated perturbatively for sufficiently small volumes where L 1/(αM ) and therefore the Bohr radius of the chargedparticle 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. [4] 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. [42]. 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 The combination αM L/(4π) in Eq. (1) 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. (1) 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. [42] by the factor of 4π in the denominator of the αM L/(4π). This manuscript is structured as follows. Section II presents the details and methodology of the LQCD+QED L calculations that are performed. Section III 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 IV 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 V. Appendix A contains further technical details on the matching between QED L and NRQED L , and Appendix B contains details on the fitting procedure used to extract energy levels from Euclidean correlation functions computed in LQCD+QED L .

II. 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 TABLE I: Parameters of the gauge field ensembles and calculations performed in this work. L 3 × T is the Euclidean spacetime volume, β is the SU (3) gauge coupling as defined in Refs. [43][44][45], β E = 1/e 2 is related to the QED gauge coupling appearing in Eq. (4), κ q are quark mass parameters for flavors q ∈ {u, d, s} appearing in the quark action, Eq. (3), N cfg is the number of gauge field configurations analyzed in this work, and N src is the average number of quark propagator sources randomly distributed on each gauge field configuration.
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 [25] 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. [3,4]. For completeness, relevant aspects of the ensemble generation are described below.

A. Lattice action and parameters
The LQCD+QED L action used in this study is given by where S G is the tree-level Symanzik-improved SU (3) gauge action as defined in Refs. [43][44][45].
The quark dynamics are encoded by an O(a)-improved stout link non-perturbative clover (SLiNC) action: whereŨ µ is a single-iterated stout-smeared SU (3) link [46], 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" [47] involves the unsmeared SU (3) gauge-field as in Ref. [48]. The clover coefficient c SW was non-perturbatively determined for pure QCD in Ref. [48]. 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: 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. II B. 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. [45] 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 1 are degenerate. As it is inspired by Dashen's theorem [49], this prescription for separating strong and electromagnetic effects is known as the "Dashen scheme" [4]. 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 I. B. Gauge fixing and the U (1) zero mode The correlation functions of charged particles are not gauge invariant and hence ensembleaveraged quantities are only meaningful within some gauge fixing prescription. In this work, Landau gauge is adopted by enforcing the condition 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. (3). 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 4dimensional zero modes onto the finite interval −π/L µ < A µ (k = 0) ≤ π/L µ [50], where 1 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 uū correction function can be interpreted in a partially quenched theory, see Ref. [4] for further discussion.
A µ (k) is the Fourier transform of A µ (x) defined explicitly in Eq. (A3). 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 [3], however this would introduce an unnecessary complication in the analysis of many-body interactions. Instead, this work adopts the QED L prescription [2,25] 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.

C. Correlation functions
The particular gauge field ensembles used in this work are described in Table I; 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 [51] (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 smearedpoint (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 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. [52][53][54][55][56]. At the sink, each colour-singlet meson bilinear is separately projected to zero momentum; for example, for a system of n pions (h = 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 Ref. [57,58] are used, again with each baryon separately projected to zero momentum at the sink.

D. 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 colour-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 where M ∈ {π + , K 0 } labels the type of meson, N states is the total number of non-thermal states included in the fit, E nM is the ground-state energy of the system with the quantum numbers of n mesons of type M , Z nM ;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 nM . In practice, fits with N states ∈ {1, 2, 3} are used in this work. The vectors E and Z 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 mM for m < n used as input for f nM as in Ref. [55]. 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: 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( 1 S 0 ), np( 3 S 1 ), pp} and three-nucleon systems with b ∈ { 3 H, 3 He}. Best-fit parameters are determined from minimization of the correlated χ 2 function for the appropriate correlation functions, G(t), and fit function, is the range of times included in the fit, and S(λ * ) is an estimate of the covariance matrix described below. Finite sample-size fluctuations may make the sample covariance matrix ill-conditioned, and shrinkage techniques [59,60] are used to obtain a numerical stable inverse covariance matrix. Following the application of shrinkage to LQCD in Ref. [61], the covariance matrix including shrinkage is defined as S(λ) = (1 − λ)C + λT , where C is the bootstrap covariance matrix and T = diagonal(C), and therefore χ 2 -minimization with S(λ) interpolates between correlated χ 2 -minimization for λ = 0 and uncorrelated χ 2 -minimization for λ = 1. The optimal shrinkage parameter λ * appearing in Eq. (10)   the expected closeness to the true covariance matrix and defined in Eq. (B2), see Appendix B and Ref. [60] 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 B.
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    Table II for the K 0 and π + single-particle energies for the L/a ∈ {32, 48} volumes. The blue band shows a constant fit to E K 0 results, and the red band shows a fit to QED L power-law FV effects at O(α/(M L) 2 ) derived in Refs. [2,29] and presented in Eq. (12). The width of the bands corresponds to 67% confidence intervals estimated using bootstrap resampling.

E. Hadron mass results
Single-meson masses computed with the gauge field ensembles used here have already been presented in Ref. [3]; the meson mass results obtained using the fitting produce employed here are shown for completeness in Fig. 4. 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. II, 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 O(α, L −2 ) prediction of QED L [2] or equivalently NRQED L [29], 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 N 3 LO and are neglected here. Fitting the L/a ∈ {32, 48} results to Eq. (12) gives a m π + = 0.15419 (29), where the lattice spacing a = 0.068(2) fm is determined in Ref. [3]. The first uncertainty in each term of Eq. (13) includes statistical and systematic fitting uncertainties as detailed in Appendix B, and the second uncertainty in the expression for m π + arises from the uncertainty    Table III for the neutron and proton single-particle energies for the L/a ∈ {32, 48} volumes. The blue band shows a fit to a constant for the L/a ∈ {32, 48} results for the neutron, and the red band shows a fit to Eq. (12) (which is valid for arbitrary charged hadrons as well as the pion) for the proton. The band widths correspond to 67% bootstrap confidence intervals. A small horizontal offset is applied symmetrically to proton and neutron results.
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), 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 III and Fig. 5. 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. (12). These relativistic FV effects are suppressed by O ((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. (12). Fits to the NLO QED L prediction for E p (L) and to a constant for E n (L) give the results aM p = 0.4037 (23), aM n = 0.3971 (25), where the uncertainties are as defined for the π + after Eq. (13). 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 α [1][2][3], which is consistent with the expected linear dependence of the proton-neutron mass difference on α, given the value α 0.1 that is used here.

III. 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. (12) 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, nonrelativistic 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.

A. Finite-volume formalism for two charged hadrons
Interactions of two electrically neutral particles with mass M and relative momentum 2p 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 cot δ(p) = − 1 a + r 2 p 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 [62,63], 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. 6. 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 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 [64], and in EFT by Kong and Ravndal [65,66], the resummed scattering amplitude is nonanalytic in η and the effective range expansion is modified as 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 [66]. Conversely, the contact interaction associated with the scattering length must be renormalized to absorb divergences from Coulomb ladder diagrams. In the MS 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 [66] 1 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 [67] and in quantum field theory [39,40]. 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πn/L with n ∈ Z 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. (17) is therefore restricted to in QED L Coulomb ladder diagrams. For sufficiently small volumes, η ≤ η L 1 and Coulomb photon exchange can be treated perturbatively.
The O(α) quantization condition relating the s-wave scattering amplitude to the twoparticle 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. [42], where η(p) is defined in Eq. (17) and where I ≈ −8.9136 is a geometric constant whose evaluation is detailed in Refs. [39] and sums over integer triplets n are restricted to |n| ≤ Λ 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. (21) and used to constrain parameterizations of the phase shift such as Eq. (18). Since Eq. (21) 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. (21) 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. (21) can be reliably used to extract Coulomb-corrected scattering parameters from LQCD+QED L results.

B. Charged two-hadron systems in NRQED L
Eq. (21) can be perturbatively expanded in powers of a C /L and other higher-order effective range expansion coefficients. Ref. [42] determined this expansion to O(αM L) and O(a C /L) 3 under the assumption that αM L a C /L 1. Following Refs. [68,69], 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 In this expression ψ is a non-relativistic hadron field, D µ = ∂ µ + ieQA µ where Q is the electric charge operator, the gauge-fixed photon Lagrangian L ξ γ is given in Eq. (A2), η 3 (µ) is a renormalization-scale-dependent coupling associated with short-range three-body interactions, and four-and higher-body interactions are neglected. L r includes effective range contributions and relativistic corrections involving two derivatives, and is given by The coefficient of the (ψ † ψ) 2 operator in Eq. (24) can be fixed by demanding that the strong interaction EFT given by replacing D µ with ∂ µ reproduces Eq. (21) with α = 0 when both are expanded perturbatively in powers of L −1 to 2 O(L −6 ). This O(L −6 ) threshold expansion of Lüscher's quantization condition is given in Ref. [72] and is verified below to be reproduced by the non-relativistic EFT defined by Eqs. (23)-(24) 3 . Operators in L r lead to contributions to the threshold expansion suppressed by ar/L 2 as well as relativistic effects suppressed by O ((M L) −2 ) that will be neglected in the power counting schemes discussed below. Additional relativistic corrections to Eq. (23) arise from photon loops and operators with four and more derivatives, but these give rise to FV effects suppressed by powers of O ((M L) −1 ) that will be neglected as discussed below and detailed in Appendix A. Introducing Fourier transformed fields the FV Hamiltonian for NRQED L is given by where H ξ γ is the gauge-fixed photon Hamiltonian, whose explicit form will not be used below, and where the two-body potential includes strong interaction and Coulomb terms, 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 pNRQED [73][74][75], but these terms give rise to O(M L) −1 effects that are neglected below. Specializing to the case of identical bosons, operators are associated with the fields in Eq. (26) and are defined to satisfy commutation relations [ ψ † p , ψ p ] = δ p,p . FV two-particle states defined by 2 The Lagrangian in Eq. (24) can also be obtained by studying the nonrelativistic limit of relativistic scalar field theory as in Refs. [70,71] . 3 Refs. [68,69] include the operators in Eq. (24), but a factor of 2 discrepancy in the coefficient of the last term leads to a difference in the O(L −6 ) threshold expansion result. satisfy p 1 , p 2 | p 1 , p 2 = 1. The ground-state of the two-particle system in its center-ofmass 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 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 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 O ((M L) −1 ), have been neglected. All sum-integral differences I, J , K, L, and R nm which appear in this expression are evaluated 4 in Ref. [42], which includes an evaluation of Eq. (31) in EFT to leading order in αM L and to higher order in a/L, except for R 44 which is a convergent sum given by R 44 ≈ 55.47.
Eq. (31) can be expressed as an expansion in a C /L and the FV Coulomb parameter η L defined in Eq. (20), 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 It is noteworthy that η L 1 only requires αM L 4π and is less restrictive than the condition αM L 1 discussed in Ref. [42]. 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 5 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 [68,69,72,78] can be used to extend Eq. (32) from NNLO in Eq. (33) to next-to-next-to-next-to-leading-order (N 3 LO) in Eq. (34), The parameter a C (L) is equal to a C plus a 1/L 3 suppressed correction arising from the interaction terms in Eq. (24), where (M L) −2 suppressed effects are shown for completeness and lead to agreement with the O(L −6 ) strong interaction relativistic threshold expansion of Ref. [72] after taking 6 α = 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 O(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. (35), which as discussed in Sec. III C below modify the dominant QED FV effects by O ((M L) −1 ). For π + π + systems on the L/a = 48 lattice volume, the dominant (NLO) QED effect in Eq. (35) amounts to a shift a C → a C − 2J (η L /π 2 )a C , and radiation photon effects can be estimated 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. (35). Therefore, a C (L) − a C will be neglected when fitting π + π + LQCD+QED L results to Eq. (35). 5 The description of FV effects on charged hadron masses in QED L as a dual expansion in α and (M L) −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. [35]. 6 The operators in Eq. (24) lead to additional effective range and relativistic corrections to the right-hand-

C. Zero-mode effects
The derivation of Eq. (21) in Ref. [42] and the form of the NRQED L Lagrangian in Eq. (23) 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. [2,29,31,35,37,38] that this assumption is violated in the single-particle sector of NRQED L because of the inherent nonlocality of zeromode 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/(M L) that have been studied in Refs [2,29,31,35,37,38]. The O (α/(M L)) and O (α/(M L) 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 O (α/(M L) 3 ), structure-dependent effects involving magnetic moments and charge radii arise. Nonlocal effects from zero-mode subtraction also enter at order O (α/(M L) 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 [31,38], or more simply by adjusting the coefficients of mass operators in the NRQED L Lagrangian by factors proportional to α/(M L) 3 [37]. For charged scalars, although not for charged fermions, these effects vanish in the charged particle rest frame [37]. This non-decoupling of antiparticle modes is a consequence of the nonlocality of NRQED L , and for EFTs with breakdown scale Λ generically produces O (α/(ΛL) 3 ) corrections to LECs [37].
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 O(M L) 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. (21) 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 O(α 2 ) diagrams containing closed loops of particle-antiparticle pairs. Two-particle FV Green's functions receive contributions from four-hadron contact interactions at O(α 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 A and summarized below.
The O(M L) 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. 6) with the "crossed-box" diagram (leftmost bottom row in Fig. 7). In NRQED L , the box diagram involves the integral 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 ∈ Z 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 (M L) −1 . The NRQED L crossedbox diagram involves the integral where only the photon pole contributes and leads to particle propagators with off-shell kinetic energies k.  [37] 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 [25], 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.
D. Charged many-hadron systems in NRQED L The two-particle energy shifts in Eqs. (32)- (35) 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 and the leading order FV energy shift for the ground state of n identical bosons in the center-of-mass frame is given by  Fig. 9, partially disconnected diagrams and diagrams that vanish because of zero-mode subtraction or kinematical constraints are not shown.
Working to NNLO in the power counting of Eq. (33), the energy shift of an n-hadron state is equal to n 2 times the two-body energy shift of Eq. (32), plus additional contributions from induced three-body and four-body forces shown in Fig. 9. Additional diagrams associated with radiation photon exchange shown in Fig. 10 are suppressed by O ((M L) −1 ). The resulting threshold expansion for the FV energy shift for a system of n like-charged hadrons at rest is given by where omitted terms are: quartic or higher in η L ∼ a C /L; relativistic effects suppressed by O ((M L) −1 ); or three-body contact interactions where η 3 ∼ a 4 C /M is assumed so that O(L −6 ) terms of the strong interaction threshold expansion appear at the same order.
At N 3 LO in the power counting of Eq. (34) 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 parameterizes the strength of six-particle operators in the NRQED L Hamiltonian, Eq. (26). Combining Eq. (41) with the N 3 LO results for the QCD threshold expansion from Ref. [68], the N 3 LO FV energy shift for a system of n like-charged hadrons in the power counting of Eq. (34) is given by The non-relativistic EFT three-body coupling η 3 (µ) is renormalization scheme-and scaledependent. The scale-dependence of η 3 (µ) cancels the explicit ln(µL) scale-dependence shown in Eq. (42), and the scheme-dependence is compensated by scheme-dependence in the finite term S MS . This scale-dependence arises from the ambiguity in separating shortdistance three-body interactions described by contact operators in NRQED L from longdistance 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 −6 ln L term in the three-body threshold expansion derived for generic relativistic field theories in Ref. [78]. Equating the O(L −6 ) term in the relativistic threshold expansion of the threeparticle threshold amplitude M 3,th in Ref. [78] with the corresponding O(L −6 ) nonrelativistic threshold expansion provides a relation between the non-relativistic coupling η 3 (µ) and the scale-independent M 3,th in the absence of QED, where a is the scattering length for a neutral two-particle system, C 3 = −0.05806 is a FV sum evaluated in Ref. [78], and S 3 = 571.398 is related to S MS , evaluated in Ref. [69], and to other FV sums from Ref. [78] by S 3 = C F + C 4 + C 5 + 8S MS . QED effects will modify Eq. (43), 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.

IV. RESULTS FOR CHARGED MULTI-HADRON SYSTEMS
This section combines the LQCD+QED L results from Section II with the NRQED L results from Section III 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.

A. Charged meson scattering
The LQCD+QED L results for the FV spectrum results in Table II 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 Results for correlated differences between n-particle ground-state energies and n times the one-particle ground-state energy, as defined in Eq. (44), 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. III. 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 B, and the results are shown in Table IV. 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 B, and the results are shown in Table V.
Results for the two-particle FV energy shifts ∆E π + π + and ∆E K 0 K 0 are shown in Fig. 11. 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 sur-  prising because αm π L ∼ 0.74 for this volume, but as discussed in Sec. III the appropriate FV analog of the Coulomb expansion parameter is η L = αm π L/(4π) ∼ 0.06 for the L/a = 48 volume. Eqs. (32)-(35) therefore predict that in addition to differences arising from a π + π + C = a K 0 K 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. 11.
The scattering lengths a π + π + C and a K 0 K 0 can be extracted from a combined fit to the results for ∆E π + π + and ∆E K 0 K 0 shown in Table IV and the results for the precisely determined correlated differences ∆E π + π + − ∆E K 0 K 0 shown in Table V. Fitting to Eq. (32), as is appropriate for the power counting PC1 of Eq. (33), gives the results, NNLO, PC1 : where the common scale m K 0 has been included for both π + π + and K 0 K 0 to facilitate comparison of a K 0 K 0 and a π + π + C . The lowest-order QED effect on ∆E π + π + from Coulomb photon exchange in Eq. (32) decreases ∆E π + π + compared to the FV energy shift for neutral particles. 8 The scattering length results in Eq. (45) show that this effect from Coulomb 8 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 zeromode subtraction [25]. 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.  Table IV for the K 0 K 0 and π + π + FV energy shifts for the L/a ∈ {32, 48} volumes. The red band shows the NRQED L predictions of Eq. (32) using the best-result for a π + π + C in Eq. (45) obtained by fitting the L/a ∈ {32, 48} results to Eq. (32). The blue band shows the prediction of Eq. (32) with α = 0 using the best-fit result for a K 0 K 0 in Eq. (45). As in Fig. 4, the widths of the bands correspond to 67% confidence intervals estimated using bootstrap resampling. A small horizontal offset is applied symmetrically to π + and K 0 results.
photon exchange competes with additional QED effects that lead to a π + π + C > a K 0 K 0 . Fitting to Eq. (35), as is appropriate for the power counting PC2 of Eq. (34), gives consistent results, demonstrating that the fit is not overly sensitive to higher-order terms absent in one or the other power counting. Both Eq. (32) and Eq. (35) neglect relativistic effects from radiation photon exchange leading to O ((M L) −1 ) FV effects. These effects are estimated in Sec. III B 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. (45)- (46) and can be consistently neglected.

B. 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.    Table VI for FV energy shifts determined from the correlated differences of two-nucleon and one-nucleon ground-state energies for the neutronneutron and proton-proton FV energy shifts for the L/a ∈ {32, 48} volumes. A smaller (larger) horizontal offset is applied symmetrically to pp and nn (np( 1 S 0 ) and np( 3 S 1 )) results.
to the LO contribution by (η p L /π 2 )(2R 24 + J 2 − L) ∼ 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 [58,[80][81][82][83][84][85][86][87]. The quantization condition in Eq. (21), which neglects O((η p L ) 2 ) perturbative Coulomb effects and (M L) −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 pp C . The two-nucleon isospin I = 1 systems pp, nn, np( 1 S 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 III. FV energy-level shifts are determined from combined analyses of the two-nucleon and single-nucleon correlation functions as described in Appendix B, and fit results for all systems are shown in Appendix B 4. As shown in Ta-ble. VI and Fig. 12, 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 O(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 pp C through the QED L quantization condition of Eq. (21) is therefore left to future work.
The two I = 1/2 three-nucleon systems 3 He and 3 H are also investigated, and their ground-state energies are given in Table III. Correlated differences between three-nucleon ground-state energies and the sums of their constituent nucleon masses are show in Table VI. As with the two-nucleon systems, fit results are shown in Appendix B 4. The results for 3 He and 3 H 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 O(50 MeV) FV energy shift at a 95% confidence level for the 3 He and 3 H 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 3 He -3 H binding-energy difference and its decomposition into QED and strong isospin breaking effects from LQCD+QED L .

C. 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. II 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. (41)-(42) both scale with n 2 , and both ∆E nπ + and ∆E nK 0 can be extracted for larger n with better relative precision than from the n = 2 case discussed in Sec. IV A.
Results for ∆E nπ + , ∆E nK 0 , and the correlated differences ∆E nπ + − ∆E nK 0 for n ∈ {2, . . . , 12} are shown in Tables II-V. QED effects leading to non-zero ∆E nπ + − ∆E nK 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 48 energy levels and correlated differences can be used to constrain the low-energy interaction parameters (42). 9 Fits to the 9 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.
NNLO expression given in Eq. (41) in PC1, which includes O(η 3 L ) 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 3 LO expression Eq. (42) 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. (42) provides a better description of the LQCD+QED L results with χ 2 /N dof ∼ 0.8. Fit results using the N 3 LO expression Eq. (42) 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. 13-14. 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, It is noteworthy that results with Zα ≥ 1 can be fit by the NRQED L formula given in Eq. (42) without additional modifications to account for relativistic QED effects or additional nonperturbative effects. Some tensions between LQCD+QED L results and N 3 LO NRQED L fits can be observed for n 8 meson systems on the L/a = 48 lattice volume in Fig. 13; 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. [88] and is given by M 3,th = 108m 2 K 0 /f 4 K 0 (with conventions for the LO Lagrangian such that f π + ∼ 130 MeV). 10 This can be combined with Eq. (43) to provide a χPT prediction for the nonrelativistic 3K 0 contact interaction, where S 3,th and C 3 are constants defined below Eq. (43), and to obtain a result entirely in terms of m K 0 and f K 0 , the LO χPT relations [89] 10 The SU (2) χPT result for M 3,th given in Ref. [88] is valid for tree-level K 0 K 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  Table IV for the FV energy shifts of multi-K 0 and multi-π + meson systems as a function of meson number n on both lattice volumes.

Shaded bands show 67% bootstrap confidence intervals for the predictions of Eq. (42) for the bestfit parameters
obtained from a global fit to the L/a ∈ {32, 48} results for n ∈ {2, . . . , 12} mesons in Tables IV-V as described in the main text. A small horizontal offset is applied symmetrically to π + and K 0 results.
have been used. The first of these relations also allows f K 0 and therefore η K 0 K 0 K 0 3 (µ) 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 0 K 0 from Eq. (47) into the first relation in Eq. (49) provides a prediction for f K 0 at the parameters of this LQCD+QED L calculation that is valid at LO  Table V for the FV energy-shift differences between multi-π + and multi-K 0 systems as a function of meson number n on both lattice volumes.
The shaded band shows the 67% bootstrap confidence interval for the corresponding prediction of Eq. (42) for the best-fit parameters a π + π + C , 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 The blue and red points show the best-fit values and 67% bootstrap confidence intervals for the three-body interaction parameters for neutral K 0 K 0 K 0 and charged π + π + π + systems respectively, including a common normalization factor of f 4 K 0 m K 0 to obtain a dimensionless quantity. The green point shows the LQCD result of Ref. [52], which was obtained in a calculation using a pseudoscalar meson mass of 352 MeV similar to m K 0 here. The black point shows the LO χPT prediction of Eq. (48) multiplied by the same normalization factor. and kaon masses; this can be compared with f π + ∼ 130 MeV and f K 0 ∼ 156 MeV extracted from experiments [90]. Inserting this result for f K 0 , and the result for a K 0 K 0 in Eq. (47), into Eq. (48) then gives the numerical result m K 0 f 4 (27). This result is valid at LO in χPT and can be compared to LQCD+QED L results for η K 0 K 0 K 0 3 (m K 0 ) combined with the result for f K 0 given in Eq. (50).
Dimensionless LQCD+QED L results for the three-body coupling that are expected to be O(1) in χPT are given by The result of this work for η K 0 K 0 K 0 3 (m K 0 ) is consistent within 1σ uncertainties with the LQCD results of Ref. [52], 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 O (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 η K 0 K 0 K 0 3 (m K 0 ), although it is also consistent within 1σ uncertainties with the LQCD results of Ref. [52]. Results for η K 0 K 0 K 0 3 (m K 0 ) and η π + π + π + 3 (m K 0 ) as well as comparisons to LO χPT and to the LQCD results of Ref. [52] are shown in Fig. 15.
Differences between η π + π + π + might also spuriously arise from mismodeling of O ((M L) −1 ) relativistic effects estimated in Sec. III B 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 η K 0 K 0 K 0 3 (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.

V. 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 0 K 0 scattering lengths as well as the 3π + and 3K 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. [42]. 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 O ((M L) −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. 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 [91]) 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 [92] [2,29,31,35,37,38]), the NRQED L Lagrangian is identical for bosons and fermions. For a particle with charge Q = 1, it is given by 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 µ = ∂ µ + iA µ , and in generic R ξ gauge the photon Lagrangian is 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, where n ∈ Z 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 [37], but for perturbative matching with QED L Eq. (A1)-(A3) define NRQED L . The NRQED L massive particle propagator G NRQED L is given by The photon propagator is given by Introducing Fourier transformed fields, 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 r with r ∈ Z 3 , At tree-level this correlation function is given by The choice r = 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. (A8), one-photon-exchange contributions to the NRQED L amplitude in Eq. (A7) vanish because of zero-mode subtraction.
Matching is performed by expanding M NRQED L and its QED L analog M QED L perturbatively in the small parameters α, 1/(M L), 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/(M L) and third order in η L = αM L/(4π) and a/L. This corresponds to NNLO in the power counting of Eq. (33) 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 3 LO in the power counting of Eq. (34), 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. (A8) can be straightforwardly matched to its QED L analog. The scalar QED L Lagrangian is where L ξ γ is defined in Eq. (A2). The fields appearing in the QED L and NRQED L Lagrangians are related by The scalar QED L particle propagator is The QED L photon propagator is identical to the NRQED L propagator Eq. (A5). Introducing Fourier transformed fieldsφ as in Eq. (A6), the FV two-particle amplitude in QED L normalized identically to the NRQED L amplitude is given by where Res(x) indicates the residue of the integrand in the first line at the pole x and n ( ) corresponds to n∈Z 3 for n γ = 0 and to n∈Z 3 \{0} for n γ ≥ 1. After taking the → 0 limit, the residue at the particle pole 2π 2 (n + r) 2 /(M L 2 ) is given by where the last line includes an expansion in powers of (M L) −1 . This expansion is legitimate provided that the sum over n converges, which holds for n γ ≥ 1. For n γ = 0 there is a linear UV divergence that can be removed by adding a UV counterterm,

(A16)
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. 6 are given by 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 ], a photon pole at k 0 = 2π|n|/L, and an antiparticle pole at Taking the → 0 limit and expanding to leading order in (M L) −1 , the residue at the particle pole is given by As in the NRQED L case, the residue at the photon pole −2π|n|/L is suppressed compared to the residue at the particle pole by O(M L) −1 . The residue at the antiparticle pole at −2M (1 + O(M L) −1 ) does not appear in the corresponding NRQED L expression Eq. (A13) and is therefore associated with contributions that do not appear in loop diagrams in NRQED L . This residue is given by The term involving two contact interactions involves the UV divergent sum n ( ) 1, which for n γ = 0 must be consistent with the infinite-volume result d D k 1 = 0 after subtracting UV counterterms, 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 n = 0 contribution to n 1 is unity, it follows from Eq. (A22) that Sums with n 2k with k > 0 similarly require UV counterterms and vanish after including them. Zero-mode contributions to these sums vanish, and so for k > 0. After subtracting UV counterterms the antiparticle pole contribution becomes (A26) 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 (M L) −1 expansion above), this contribution is proportional to the LO contact interaction times the velocity squared times a relativistic QED suppression factor of α/(M L) −3 . For boosted systems, the NRQED L contact interaction in Eq. (A1) must therefore be supplemented with a nonlocal counterterm suppressed by the same factor of α/(M L) −3 . This is analogous to the self-energy of a single scalar field, which is shown in Ref. [37] to require a nonlocal counterterm equal to the scalar mass times the velocity squared times a relativistic QED suppression factor of α/(M L) −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 (M L) −1 are neglected throughout this work it is consistent to neglect all nonlocal counterterms suppressed by α/(M L) −3 .
The NLO QED L amplitude is therefore given by inserting the particle-pole residue in Eq (A20) into Eq. (A19),  Figs. 7 and 8 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-twophoton vertex or particle-antiparticle pair creation can be similarly verified to be suppressed by powers of (M L) −1 . The NRQED L Lagrangian in Eq. (A1) therefore reproduces QED L at NLO.
Matching at NNLO proceeds similarly. The NNLO diagrams shown in Fig. 6 can all be expressed in terms of the amplitude αM a where n i ∈ {0, } for i ∈ {1, 2, 3} labels whether each interaction is a four-particle contact interaction or photon exchange and n,m ( ) excludes n = 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 Contributions from photon and antiparticle poles are again suppressed by powers of (M L) −1 , and evaluating the q 0 and k 0 energy integrals gives All other diagrams are again suppressed by powers of (M L) −1 , and agreement between Eq. (A29) and Eq. (A31) shows that nonlocal counterterms are suppressed by powers of (M L) −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 (M L) −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 (M L) −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 (M L) −1 suppressed relativistic effects, and in particular the scalar QED L antiparticle pole in Eq. (A21) 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 O(α/(M L) 3 ) for fermions with any center-of-mass velocity and for scalars with non-zero center-of-mass velocity.
Appendix B: Fitting procedures

Energy level determination
In general the spectral representation in Eq. (8) 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. (8) 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 2 mπ)t predicted by Parisi [93] and Lepage [94]  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 finitetemperature 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 noiseto-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, 12 or t + a is larger than a finitetemperature cutoff (final results use tol temp = 3T /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. 11 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 [95]. 12 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. 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 [96] 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 C 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 2 pts 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 2 pts is not satisfied [59,60]. Shrinkage estimators S ij tt (λ) of the covariance matrix are constructed as mixtures of a well-conditioned target matrix T ij tt and the covariance matrix C ij tt estimated using standard bootstrap techniques from N boot = 200 samples of N correlation functions G b i,a (t) with a ∈ {1, . . . , N } and b ∈ {1, . . . , N boot } drawn from the original correlation function ensemble with replacement, 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. [61], we take T ij tt = diag(C ij tt ). 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 [60] that for finite N an optimal λ * = 0 satisfying this restriction can be chosen in order to minimize the average mean-squared difference between S ij tt (λ) and the underlying covariance matrix, and that a sample estimate for λ * is given by where are defined in terms of the sample mean correlation function and sample covariance as in Ref. [61] ) with the identity matrix as a target, and the results of Ref. [60] assuming an identity matrix target can be applied. The covariance matrix estimate with optimal shrinkage is then given by S ij tt (λ * ). Fits to truncations of Eq. (8) including e excited states can then be performed by minimizing the corresponding χ 2 function defined by where i includes a 1/N bias correction estimated using bootstrap techniques and E and Z denote the energies and overlap factors appearing in Eq. (8), including excited-state and thermal effects. Since the overlap factors enter f N linearly, the values of Z minimizing χ 2 for fixed E can be determined by solving a system of linear equations analogous to variable projection techniques [97,98]. χ 2 -minimization can therefore be efficiently performed by using a nonlinear optimization method to determine E with the optimal Z determined by solving a system of linear equations at each step of nonlinear optimization for E. In order to ensure positivity of the spectrum and remove fitting degeneracies, the parameters used for nonlinear optimization are ln E 0 and ln( 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 [99] (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) < −AN 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. (B5) 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 O(1) value A compared to the (e − 1) state fit. For baryon correlation functions a value of A = 0.5 is used, while for multi-meson correlation functions a value of A = 0.1 is used. 13 Bootstrap resampling techniques are then used to estimate the uncertainty on the groundstate energy extracted from the fit with the preferred number of excited-states in each fit region [100,101]. The same χ 2 -minimization procedure and estimated covariance matrix S(λ * ) 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 f 0 obtained for fit range 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 f 0 minimizes the impact of outlier bootstrap samples compared to a definition based on the standard deviation of E b,f 0 − E f 0 [101]. 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 14 and are verified to give energies that differ by less than a specified tolerance 15 (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 S(λ = 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. 16. This fitting procedure was implemented in the Julia language [102] using the Optim optimization package [103] to obtain the results of this work.
Fits that pass all of the checks above are considered reliable estimates of the energy spectrum, and the final estimate of the ground-state energy E 0 and its uncertainty δE 0 are 13 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 A from 0.1 over the range A ∈ [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. 14 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. 15 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.
Excited-state model selection: Z n e −Ent , e = 0

Confidence intervals:
Covariance matrix S(λ) with optimal shrinkage parameter λ * FIG. 16: Flowchart representing the steps of the fitting procedure for one specific fitting range. Rectangular shapes represent process steps, while diamond shapes represent decision steps. Input parameters to the fitting procedure are shown in blue. As described in the text, the steps illustrated here are repeated N fits times with different random choices of t min , and final results are obtained from weighted averages of fit results for the t min choices leading to the "Accept fit" rectangle.
obtained by taking a weighted average of the N success successful fit results E f 0 , where f labels the choice of fit range specified by t min for each interpolating operator 16 . Each fit result provides an unbiased estimate of the ground-state energy. The relative weights w f of each fit in the weighted average can therefore be chosen arbitrarily in the limit of large statistics; in practice it is advantageous to choose weights that penalize poor fits with larger χ 2 /N dof and unconstraining fits with larger uncertainties δE f 0 . Following Ref. [61], we use the weights where p f = Γ(N dof /2, χ 2 f /2)/Γ(N dof /2) is the p-value assuming χ 2 -distributed goodness-offit parameters with χ 2 f obtained by inserting E f 0 into Eq. (B5) 17 . Variation to the particular choices of specified tolerances have been studied, and the subsequent variation in the ensemble of successful fits is found to have little impact on the results of this weighted averaging. The results E 0 and δE 0 obtained with this procedure are shown as the central values and uncertainties for single-particle energy results E π + (L), E K 0 (L), E n (L), and E p (L) in Tables II and III. Effective mass plots showing the smallest t min fit with weight over 1/2 of the maximum weight fit as well as E f 0 and w f /Max( w f ) are shown in Appendix B 4.

Multi-meson correlation functions
To determine results for multi-meson ground-state energies with thermal effects taken into account, fits are performed iteratively starting with fits for n = 1 mesons and then 16 The total error δE 0 describes the combined statistical uncertainty on E 0 plus systematic uncertainty arising from the choice of fit range and fit model. The partitioning of this error into δ stat E 0 and δ sys E 0 only partially separates statistical and systematic uncertainties because δ stat E 0 includes statistical errors plus systematic uncertainties related to fluctuations among the δE f 0 . 17 For large λ * , the χ 2 function being minimized approaches an uncorrelated χ 2 and the values of χ 2 will not be distributed as χ 2 -distributed random variables with N dof degrees-of-freedom. In this regime where finite N artifacts are not negligible, the weights in Eq. (B8) still serve the purpose of penalizing comparatively less accurate descriptions of the results being fit and their correlations as estimated by S ij tt (λ * ), but the absolute sizes of the p f should not be interpreted as p-values for each fit. moving to fits with increasing n. The excited-state fit form in Eq. (8) includes thermal effects describing k forwards-propagating mesons and n − k backwards-propagating mesons for k < n/2 that are included by using the central values E 0 calculated for E k and E n−k in Eq. (8). Uncertainties in E k are found to be significantly smaller than uncertainties in E n for k < n/2, and for simplicity are not incorporated into the description of thermal effects. The overlap factors for these thermal states are determined using linear algebra techniques [97,98] during each step of nonlinear optimization for the N -particle energy spectrum analogously to the procedure described above for other overlap factors. This fit function is found to provide acceptable fits to multi-meson correlation functions without the need for additional free parameters describing excited-state thermal effects.
Before beginning the fitting procedure described above, correlation function results from all quark propagator sources on a given configuration are averaged and meson correlation functions are further blocked along the Markov chain to form N block = 200 approximately independent samples from the N cfg configurations for each volume shown in Table I. Further averaging is found to give statistically consistent results, suggesting that autocorrelations can be neglected after this blocking. To determine correlated differences of ground-state energies for different hadron type (π + , K 0 ) and hadron number, fits are performed independently to determine E f 0 for each hadron type and number for each fit range sampled. A fit range is considered to give a successful fit only if the checks on fit robustness described above are passed for each hadron type and number involved in the correlated difference. For each successful fit range, bootstrap resampling is used to determine the uncertainties on correlated differences of the resulting E f 0 . During bootstrap resampling, the same elements of these N block samples are used to construct bootstrap ensembles for each hadron type and number. Correlated differences of the bootstrap results E f,b 0 are then formed, and confidence intervals are computed by applying Eq. (B6) to the these correlated differences. Finally, weighted averages of the resulting correlated differences and their associated uncertainties are taken using Eq. (B7)-(B8). The results of this procedure are used to determine the FV energy shifts and differences between FV energy shifts for charged and uncharged hadrons shown in Tables IV-III.

Multi-nucleon correlation functions
Differences between multi-nucleon ground-state energies and the corresponding sums of their constituent nucleon masses are computed using correlated differences of bootstrap results E f,b 0 for multi-nucleon and single-nucleon correlation functions analogously to the multi-meson case described above. Correlated differences between one-nucleon and multinucleon energies can be determined much more precisely than multi-nucleon energies alone, and differences of one-nucleon and multi-nucleon fit results with different values of N states are found to describe correlated differences of LQCD+QED L effective energy results poorly. N states is therefore restricted to be identical between single-nucleon and multi-nucleon systems. Otherwise, fits for multi-nucleon energy shifts are performed identically to fits for multi-meson energy shifts not including thermal effects.

Fit results
Figs. 17-19 show fit results for nπ + systems with L/a = 32. Results with L/a = 48 are shown in Figs. 1-3. Figs. 20-25 show analogous fit results for nK 0 systems with L/a ∈ {32, 48}. Single-nucleon fit results for p and n are shown in Fig. 26. Two-nucleon fit results are shown for pp and nn in Fig. 27 and for np systems in Fig. 28. Three-nucleon fit results are shown in Fig. 29.