Charged multi-hadron systems in lattice QCD+QED
S. R. Beane, W. Detmold, R. Horsley, M. Illa, M. Jafry, D. J. Murphy, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, P. E. Shanahan, H. Stüben, M. L. Wagman, F. Winter, R. D. Young, J. M. Zanotti
DDESY 20-28, FERMILAB-PUB-20-123-T, ICCUB-20-007, MIT-CTP/5183, NT@UW-20-03
Charged multi-hadron systems in lattice QCD+QED
S. R. Beane, W. Detmold, R. Horsley, M. Illa, M. Jafry, D. J. Murphy, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, P. E. Shanahan, H. Stüben, M. L. Wagman,
2, 10
F. Winter, R. D. Young, and J. M. Zanotti (NPLQCD and QCDSF collaborations) Department of Physics, University of Washington,Box 351560, Seattle, WA 98195, U.S.A. Center for Theoretical Physics, MassachusettsInstitute of Technology, Cambridge, MA 02139, U.S.A. School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD, UK Departament de Física Quàntica i Astrofísica,Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona,Martí Franquès 1, E08028 Barcelona, Spain RIKEN Center for Computational Science, Kobe, Hyogo 650-0047, Japan Institut für Theoretische Physik, Universität Leipzig, 04109 Leipzig, Germany Theoretical Physics Division, Department of Mathematical Sciences,University of Liverpool, Liverpool L69 3BX, UK Deutsches Elektronen-Synchrotron DESY, 22603 Hamburg, Germany Regionales Rechenzentrum, Universität Hamburg, 20146 Hamburg, Germany Fermi National Accelerator Laboratory, Batavia, IL 60510, U.S.A. Jefferson Laboratory, 12000 Jefferson Avenue, Newport News, VA 23606, USA CSSM, Department of Physics, University of Adelaide, Adelaide SA 5005, Australia
Systems with the quantum numbers of up to twelve charged and neutral pseu-doscalar mesons, as well as one-, two-, and three-nucleon systems, are studiedusing dynamical lattice quantum chromodynamics and quantum electrodynamics(QCD+QED) calculations and effective field theory. QED effects on hadronic inter-actions are determined by comparing systems of charged and neutral hadrons aftertuning the quark masses to remove strong isospin breaking effects. A non-relativisticeffective field theory, which perturbatively includes finite-volume Coulomb effects, isanalyzed for systems of multiple charged hadrons and found to accurately reproducethe lattice QCD+QED results. QED effects on charged multi-hadron systems beyondCoulomb photon exchange are determined by comparing the two- and three-body in-teraction parameters extracted from the lattice QCD+QED results for charged andneutral multi-hadron systems.
I. INTRODUCTION
The interplay between the strong and electromagnetic (EM) interactions of the StandardModel is subtle but central to the complexity of visible matter. While EM interactionsare typically much weaker than the strong interactions, it is their competition with strongisospin-breaking effects that leads to the observed proton–neutron mass difference and de-termines the stability of atomic nuclei. Moreover, the hierarchy between the length scaleswhere strong and EM interactions are important leads to the existence of chemistry andall of the complexity that it entails. Progress towards understanding the combined effects a r X i v : . [ h e p - l a t ] A p r of the strong and EM interactions from first principles has been made in recent years, andthe 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 lat-tice formulation of quantum chromodynamics (QCD) and quantum electrodynamics (QED)[1–3]. Coupled lattice QCD+QED calculations have also been performed for meson andbaryon masses [4–7], leptonic decay rates [8–10], and the anomalous magnetic moment ofthe muon [11–19], as both QCD and QED must be accounted for in precision tests of theStandard Model.Systems of multiple charged hadrons exhibit rich phenomenology in which the interplaybetween QCD and QED effects is less well understood than for single-hadron systems. Thedifferences between pp , np and nn two-nucleon scattering channels are much more poorlyconstrained than their isospin average. Lattice QCD+QED can potentially provide phe-nomenologically useful information on such isospin-breaking effects in nucleon-nucleon scat-tering; furthermore, calculations with a range of quark masses and values of the QED finestructure constant can provide insight into the fine-tuning of QED effects and strong isospinbreaking interactions, providing complementary information to experiment. In large nucleiwith charge Z ∼ /α , where α is the EM fine structure constant, relativistic QED effectsare expected to become nonperturbative and give rise to electron-positron pair-productionand interesting consequences for the structure of elements [20].QCD+QED studies of systems involving two or more protons are of great phenomeno-logical interest, but calculations of large nuclei are difficult using current lattice QCD+QEDtechniques [21]. It is also interesting to consider systems of multiple charged mesons, whichallow an investigation of similar QED effects without of number of the numerical complexi-ties of multi-proton calculations. In particular, charged multi-hadron systems include QEDeffects that become non-perturbative for hadron pairs with sufficiently small relative velocity v (cid:28) α . All calculations of charged hadrons must also reconcile the presence of net chargein 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], inwhich the photon spatial zero-mode is removed, as has been recently studied in detail forsingle-hadron systems in Refs. [1–3, 5, 6, 8, 9, 19, 22–27, 29, 31–37]. Subtleties related tothe nonlocality of zero-mode subtraction in QED L have been understood for single-hadronsystems [2, 29, 31, 35, 37, 38], and similar methods are used to understand nonlocality inmulti-hadron systems in this framework as described below.In lattice QCD calculations, hadronic scattering phase shifts are determined by relat-ing them to the FV energies of two-hadron systems using FV relations first derived byLüscher [39, 40] and generalized in recent years (see Ref. [41] for a review). Understand-ing scattering in lattice QCD+QED is complicated by the lack of a FV formalism thatincludes nonperturbative Coulomb effects, and thereby relates FV energy shifts calculablein lattice QCD+QED to infinite-volume (Coulomb-corrected) scattering phase shifts. In theQED L formulation, Ref. [42] argues that Coulomb effects can be treated perturbatively forsufficiently small volumes where L (cid:28) / ( αM ) and therefore the Bohr radius of the charged-particle system is larger than the infrared cutoff provided by the FV. However, the volumemust also satisfy L (cid:29) /M in order for relativistic FV effects not accounted for in Lüscher’sanalysis to be negligible. For some systems there may be a window of intermediate-sizedvolumes where relativistic FV effects and nonperturbative Coulomb effects can both be ne-glected. Non-relativistic effective field theory (EFT) can be used to investigate this issueand to relate FV energy shifts for systems with more than two particles to two-, three-, andhigher-body EFT low-energy constants (LECs).This work studies QCD+QED effects for systems of multiple hadrons including up totwelve charged or neutral mesons and up to three nucleons in multiple finite volumes. Alarger than physical value of the fine-structure constant, α (cid:39) . , is used in order to increasethe size of QED effects and permit the study of systems with total charge Z satisfying Zα (cid:38) . As described in Ref. [4] and summarized below, the quark masses are tunedto remove strong isospin breaking effects so that splittings between particles that are de-generate in QCD can be identified as pure QED effects. Calculations are performed inboth lattice QCD+QED L (LQCD+QED L ) and in an EFT appropriate for non-relativisticcharged hadrons, referred to as non-relativistic QED L (NRQED L ) below. Two-body scat-tering lengths and three-body interactions of charged and neutral mesons are determinedby tuning the hadronic interaction parameters appearing in the NRQED L Lagrangian toreproduce the LQCD+QED L FV energy levels. QED effects on hadron interactions beyondCoulomb photon exchange are determined by comparing the two-body and three-body in-teraction parameters determined for charged and neutral mesons. Studies of multi-baryonsystems are used to probe volumes where
L > / ( αM ) and Coulomb effects may not beperturbative according to the criteria of Ref. [42]. It is noteworthy that no significant non-perturbative Coulomb effects or relativistic effects are seen at the level of precision of theresults obtained here, even for systems with L > / ( αM ) or Zα > .The practical window of L where relativistic effects can be neglected and Coulomb effectscan be treated perturbatively is explored by comparing LQCD+QED L and NRQED L resultsfor a variety of FV systems. From these studies, it is argued that this window requires thespatial extent L of a cubic FV to satisfy M L (cid:28) , αM L π (cid:28) . (1)The combination αM L/ (4 π ) in Eq. (1) that quantifies the size of FV Coulomb effects isequivalent to the infinite-volume Coulomb expansion parameter α/v for a pair of hadronsmoving back-to-back with one unit of FV momentum, i.e. ± π/L . Eq. (1) indicates thatCoulomb effects can be treated perturbatively in current and future LQCD+QED L calcula-tions over a wide range of volumes and in particular for systems with L ∼ / ( αM ) , wherethe Bohr radius is commensurate with the volume. This scaling differs from that suggestedin Ref. [42] by the factor of π in the denominator of the αM L/ (4 π ) .This manuscript is structured as follows. Section II presents the details and methodologyof the LQCD+QED L calculations that are performed. Section III discusses the constructionof NRQED L for charged multi-hadron systems and provides the formalism needed to extracthadronic scattering lengths and other parameters from the results of these LQCD+QED L calculations. Section IV describes the determination of hadronic scattering parameters fromLQCD+QED L and presents results for multi-meson and multi-nucleon systems. Conclusionsare presented in Section V. Appendix A contains further technical details on the matchingbetween QED L and NRQED L , and Appendix B contains details on the fitting procedure usedto 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 intermediatestages implements an ultraviolet regulator defined by a lattice spacing a (where /a is L × T β β E κ u κ d κ s N cfg N src × × L × T is the Euclidean spacetime volume, β is the SU (3) gauge coupling as defined in Refs. [43–45], β E = 1 /e is related to the QED gauge coupling appearing in Eq. (4), κ q are quark mass parametersfor flavors q ∈ { u, d, s } appearing in the quark action, Eq. (3), N cfg is the number of gauge fieldconfigurations analyzed in this work, and N src is the average number of quark propagator sourcesrandomly distributed on each gauge field configuration. assumed to be much smaller than the QED Landau pole). Calculations are performed inEuclidean spacetime with a cubic spatial volume of extent L × L × L and a finite temporalextent 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-particlecorrelation functions as detailed below, which defines the LQCD+QED L formalism. TheLQCD+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]. Forcompleteness, 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 S = S G + S A + S F , (2)where S G is the tree-level Symanzik-improved SU (3) gauge action as defined in Refs. [43–45].The quark dynamics are encoded by an O ( a ) -improved stout link non-perturbative clover(SLiNC) action: S F = (cid:88) q ∈{ u,d,s } (cid:88) x (cid:40) (cid:88) µ (cid:104) ¯ q ( x )( γ µ − e − ie q A µ ( x ) ˜ U µ ( x ) q ( x + ˆ µ ) − ¯ q ( x )( γ µ + 1) e ie q A µ ( x − ˆ µ ) ˜ U † µ ( x − ˆ µ ) q ( x − ˆ µ ) (cid:105) + 12 κ q ¯ q ( x ) q ( x ) − c SW (cid:88) µν ¯ q ( x ) σ µν F µν q ( x ) (cid:41) , (3)where ˜ U µ is a single-iterated stout-smeared SU (3) link [46], A µ is a non-compact U (1) gaugefield, and the U (1) quark charges are given by e q . The field-strength F µν appearing in theSheikholeslami-Wohlert or “clover term” [47] involves the unsmeared SU (3) gauge-field asin Ref. [48]. The clover coefficient c SW was non-perturbatively determined for pure QCDin Ref. [48]. Electromagnetic gauge fields are not included in the clover term, however withthese simulation parameters, the O ( αa ) effects are no larger than the residual O ( a ) effectsof pure QCD. The photon action is described by the non-compact form: S A = 12 e (cid:88) x,µ<ν (cid:18) A µ ( x ) + A ν ( x + µ ) − A µ ( x + ν ) − A ν ( x ) (cid:19) (4)where e is the U (1) gauge coupling corresponding to α = e / (4 π ) . Gauge fixing and thetreatment 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—seeRef. [45] for further discussion. With dynamical QED, this is complicated by the fact thatthe quark charges explicitly break the SU (3) flavor symmetry. An approximate SU (3) fla-vor symmetry is then realized by tuning the quark mass parameters such that the connected flavour-neutral pseudoscalar mesons are degenerate. As it is inspired by Dashen’s theo-rem [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 pseu-doscalar mesons are protected from receiving an electromagnetic self-energy correction inthe chiral limit. To reach the physical quark masses and charges, the SU (3) flavor symmetricpoint can then be lowered to the physical value before the symmetry is broken to fix theindividual quark masses. In this exploratory work a single set of approximately SU (3) flavorsymmetric 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 aspure QED effects. The explicit breaking of the quark masses from this point then allows oneto independently isolate the effects of strong (or quark mass) flavor symmetry breaking andQED effects. The action parameters used in this work correspond to the SU (3) symmetricpoint 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 ensemble-averaged quantities are only meaningful within some gauge fixing prescription. In this work,Landau gauge is adopted by enforcing the condition (cid:88) µ ( A µ ( x ) − A µ ( x − ˆ µ )) = 0 . (5)This condition can be imposed after generation of the gauge fields. However, the Landaugauge condition leaves the 4-dimensional zero mode unconstrained. These uniform back-ground fields do not contribute to the gauge action, however they do couple to the fermionicaction, Eq. (3). The action remains invariant under discrete shifts of the zero mode in unitsof π/L µ . This redundant gauge degree of freedom can be eliminated by mapping the 4-dimensional zero modes onto the finite interval − π/L µ < (cid:101) A µ ( k = 0) ≤ π/L µ [50], where Since κ s = κ d and e s = e d , this theory exhibits exact U -spin symmetry, and therefore the connected d ¯ d and s ¯ s correlation functions are identical to the d ¯ s correlation functions The connected part of the u ¯ u correction function can be interpreted in a partially quenched theory, see Ref. [4] for further discussion. (cid:101) A µ ( k ) is the Fourier transform of A µ ( x ) defined explicitly in Eq. (A3). The leading effect ofthese zero modes is to induce a charge-dependent twist of the single-hadron momenta. Thissmall energy shift can be corrected in single-hadron states [3], however this would introducean unnecessary complication in the analysis of many-body interactions. Instead, this workadopts the QED L prescription [2, 25] of setting the spatial (3-dimensional) zero modes ofthe gauge potential to zero on every timeslice.With respect to the action used to generate the gauge configurations, the elimination ofthe 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 experiencea different zero mode to the quarks in the sea. The zero modes cannot affect closed fermionloops, and their only contribution to the determinant will be associated with fermion lineswrapping around the boundary of the lattice. Consequently, this partial-quenching effect isexponentially suppressed in m π L and is negligible in comparison to the power-law FV effectsstudied in this work. C. Correlation functions
The particular gauge field ensembles used in this work are described in Table I; alongwith the parameters used in their generation, the number of configurations and the averagenumber of correlation function source locations that are used per configuration are alsoreported. Up and down/strange (equivalent for the masses used here) quark propagatorsare computed from each of the randomly chosen source locations using three-dimensionallyJacobi smeared sources [51] (100 iterations with ρ = 0 . ) using a solver tolerance of − .Local and smeared fields are used in the sink interpolating operator to construct smeared-point (SP) and smeared-smeared (SS) correlation functions with the former providing cleanersignals in all cases.FV energy levels are determined by analyzing two-point correlation functions G h ( t ; x ) = (cid:42)(cid:88) x ˜ χ h ( x , t ) χ † h ( x , (cid:43) , (6)where χ h ( ˜ χ h ) is a source (sink) interpolating operator with the quantum numbers of thestate being considered and x is the source location. Two-point correlation functions are con-structed for n ∈ { , . . . , } charged pions ( u ¯ d ) and neutral kaons ( s ¯ d ) using the techniquesdeveloped in Refs. [52–56]. At the sink, each colour-singlet meson bilinear is separatelyprojected to zero momentum; for example, for a system of n pions ( h = nπ + ), G nπ + ( t ; x ) = (cid:42)(cid:32)(cid:88) x ¯ u ( x , t ) γ d ( x , t ) (cid:33) n (cid:0) u ( x , γ ¯ d ( x , (cid:1) n (cid:43) . (7)The correlation functions for the single proton and neutron make use of standard localinterpolating operators χ p ; α = (cid:15) ijk ( u i Cγ d j ) u kα and χ n ; α = (cid:15) ijk ( d i Cγ u j ) d kα where the paren-theses indicate contraction of spin indices. For two-baryon and three-baryon systems, thecontraction techniques of Ref. [57, 58] are used, again with each baryon separately projectedto zero momentum at the sink. D. Finite-volume energy level determinations
Finite-volume energy levels are extracted from the correlation functions for each systemusing correlated fits to their time dependence. For the multi-meson systems, which factorizeeasily into colour-singlet components, thermal contributions where one component propa-gates forward in time and another propagates backwards in time are particularly relevant.The corresponding functional forms that are used to fit these correlation functions are f nM ( t, E , Z ) = n (cid:88) m =0 Z nM ; m e − E mM t e − E ( n − m ) M ( T − t ) + N states − (cid:88) e =1 Z ( e ) nM e − E ( e ) nM t , (8)where M ∈ { π + , K } labels the type of meson, N states is the total number of non-thermalstates included in the fit, E nM is the ground-state energy of the system with the quantumnumbers of n mesons of type M , Z nM ; m is the overlap factor describing the amplitude ofthermal contributions with m forwards propagating mesons and n − m backwards propagat-ing mesons, and non-thermal excited states are also included with energies E ( e ) nM and overlapfactors Z ( e ) nM . In practice, fits with N states ∈ { , , } are used in this work. The vectors E and Z indicate the energy and overlap factor parameters to be constrained in the fit. In orderto determine the many parameters of these fitting functions, fits are performed sequentiallyfor 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 betweenthe source and the sink, and the contributions of thermal states are negligible with respectto the statistical uncertainties. For these systems a simpler fit function is used: f b ( t, E , Z ) = Z b e − E b t + N states − (cid:88) e =1 Z ( e ) b e − E ( e ) b t , (9)where b labels the type of baryon system and, as for meson systems, the second termcorresponds to a sum over the excited states included in the fitting model. This workinvestigates single-nucleon systems with b ∈ { p, n } as well as two-nucleon systems with b ∈ { nn, np ( S ) , np ( S ) , pp } and three-nucleon systems with b ∈ { H , He } .Best-fit parameters are determined from minimization of the correlated χ function χ ( E , Z ) = t max (cid:88) t,t (cid:48) = t min ( G ( t ) − f ( t, E , Z ))( S ( λ ∗ ) − ) t,t (cid:48) ( G ( t (cid:48) ) − f ( t (cid:48) , E , Z )) (10)for the appropriate correlation functions, G ( t ) , and fit function, f ∈ { f nM , f b } , [ t min , t max ] is the range of times included in the fit, and S ( λ ∗ ) is an estimate of the covariance matrixdescribed below. Finite sample-size fluctuations may make the sample covariance matrixill-conditioned, and shrinkage techniques [59, 60] are used to obtain a numerical stableinverse covariance matrix. Following the application of shrinkage to LQCD in Ref. [61], thecovariance matrix including shrinkage is defined as S ( λ ) = (1 − λ ) C + λ T , where C is thebootstrap covariance matrix and T = diagonal ( C ) , and therefore χ -minimization with S ( λ ) interpolates between correlated χ -minimization for λ = 0 and uncorrelated χ -minimizationfor λ = 1 . The optimal shrinkage parameter λ ∗ appearing in Eq. (10) is chosen to maximize aE nK aE nπ + n L/a = 32 L/a = 48
L/a = 32
L/a = 48 n ∈ { , . . . , } neutral K mesons andcharged π + mesons with lattice volumes L/a shown. Results are determined by fitting Eq. (8) tothe LQCD+QED L Euclidean correlation functions using the methods described in Appendix B.Dimensionful energy results can be obtained using the lattice spacing a = 0 . fm obtained forthis gauge field ensemble in Ref. [3]. the expected closeness to the true covariance matrix and defined in Eq. (B2), see Appendix Band Ref. [60] for further discussion.Systematic uncertainties arise from the dependence of the fits on the functional forms andtime 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 usedto select the appropriate number of excited states to include for each fit range choice. Aweighted average of the results from all successful fits passing various reliability checks isused 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 functionsare performed using common energies but different overlap factors for local and smearedsources. For mesons and three-nucleon systems, combined fits performed in this way areonly marginally more precise than fits using SP correlation functions alone, and fits usingonly SP correlation functions are therefore used in what follows for simplicity.An effective energy plot that removes thermal effects from backwards propagating statesfor n = 1 systems, as well as constant contributions from thermal effects on n ≥ correlationfunctions, is defined as E ( t ) = ArcCosh (cid:18) G ( t ) − G ( t + 4)2 [ G ( t + 1) − G ( t + 3)] (cid:19) . (11)Figs. 1-3 show effective energy plots, including correlation function results and fit resultsfor E ( t ) for each case of n ∈ { , . . . , } π + on the L/a = 48 lattice volume as well asthe ground-state energy results from all successful fit range choices. Analogous results forthe
L/a = 32 volume and for n = { , . . . , } K systems on the L/a ∈ { , } volumesare shown in Appendix B. The resulting ground-state energies for n ∈ { , . . . , } mesonsystems on both lattice volumes are shown in Table II. π + L/a = 48 t / a . . . . . . E a SP f . . . . . . E f a π + L/a = 48 t / a . . . . . . E a SP f . . . . . . . . E f a π + L/a = 48 t / a . . . . E a SP f . . . . . . . E f a π + L/a = 48 t / a . . . . . . E a SP f . . . . . . . E f a FIG. 1: Fit results for n ∈ { , . . . , } π + systems with L/a = 48 . Blue points in the left plots showLQCD+QED L results for E ( t ) defined in Eq. (11). Blue bands show confidence intervals fromfits to Eq. (8) described in Appendix B. Horizontal light (dark) gray bands show the statistical(total) uncertainty of the fitted ground-state energy. The right plot shows ground-state energyresults from each successful fit range with opacity equal to their relative weight in the averagedetermining the total statistical plus systematic uncertainty shown as a pink band. π + L/a = 48 t / a . . . . . . . . E a SP f . . . . . . E f a π + L/a = 48 t / a . . . . . . . . E a SP f . . . . . E f a π + L/a = 48 t / a . . . . E a SP f . . . . . . . E f a π + L/a = 48 t / a . . . . . E a SP f . . . . . . . . E f a FIG. 2: Fit results for systems of n ∈ { , . . . , } π + mesons for the L/a = 48 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. π + L/a = 48 t / a . . . . . . E a SP f . . . . . E f a π + L/a = 48 t / a . . . . . . . E a SP f . . . . . . . . E f a π + L/a = 48 t / a . . . . . . . . E a SP f . . . . . . . E f a π + L/a = 48 t / a . . . . . . . . E a SP f . . . . . . . . . E f a FIG. 3: Fit results for systems of n ∈ { , . . . , } π + mesons for the L/a = 48 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. □ □○ ○ □ ○ FIG. 4: The blue and red points show the LQCD+QED L results from Table II for the K and π + single-particle energies for the L/a ∈ { , } volumes. The blue band shows a constant fit to E K results, and the red band shows a fit to QED L power-law FV effects at O ( α/ ( M L ) ) derivedin Refs. [2, 29] and presented in Eq. (12). The width of the bands corresponds to 67% confidenceintervals estimated using bootstrap resampling. E. Hadron mass results
Single-meson masses computed with the gauge field ensembles used here have already beenpresented in Ref. [3]; the meson mass results obtained using the fitting produce employed hereare shown for completeness in Fig. 4. The K energy is equal within statistical uncertaintiesfor the L/a = 48 and
L/a = 32 lattice volumes. Volume dependence in E π + ( L ) anddifferences in the FV single-particle energies E π + ( L ) − E K ( L ) , however, are clearly visible.Because of the quark mass tuning described in Sec. II, which is designed to remove strongisospin breaking effects, these energy differences are ascribed to QED effects. In order tointerpret these QED effects, results for E π + ( L ) can be compared to the O ( α, L − ) predictionof QED L [2] or equivalently NRQED L [29], E π + ( L ) = m π + α L c (cid:18) m π L (cid:19) + O (cid:0) α , L − (cid:1) , (12)where c = − . is a geometric constant that does not depend on the system under con-sideration. Nonlocal effects from zero mode subtraction and charge radius effects introducinga new parameter both arise at N LO and are neglected here. Fitting the
L/a ∈ { , } results to Eq. (12) gives a m π + = 0 . , m π + = 449(1)(13) MeV , (13)where the lattice spacing a = 0 . fm is determined in Ref. [3]. The first uncertainty ineach term of Eq. (13) includes statistical and systematic fitting uncertainties as detailed inAppendix B, and the second uncertainty in the expression for m π + arises from the uncertainty3 aE b b L/a = 32 L/a = 48 p n pp np ( S ) nn np ( S ) He 1.294(77) 1.268(42) H 1.291(53) 1.294(31)TABLE III: Ground-state energies of systems of n ∈ { , , } protons and neutrons determined byfitting Eq. (9) to the LQCD+QED L Euclidean correlation functions as described in Appendix B. □ □□○ ○○ □ ○
FIG. 5: The blue and red points show the LQCD+QED L results from Table III for the neutronand proton single-particle energies for the L/a ∈ { , } volumes. The blue band shows a fit to aconstant for the L/a ∈ { , } 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 widthscorrespond to 67% bootstrap confidence intervals. A small horizontal offset is applied symmetricallyto proton and neutron results. in a . A constant fit to the L/a ∈ { , } results for E K ( L ) gives a m K = 0 . , m K = 404(1)(12) MeV . (14)For the values of the quark masses and α (cid:39) . used in this work, strong isospin breakingeffects approximately vanish and the difference between the charged and neutral pseudoscalarmeson masses m π + − m K = 45(2) MeV is attributed entirely to QED effects. It is noteworthythat QED effects account for approximately of the π + mass with these parameters.Results for the proton and neutron ground-state energies are shown in Table III andFig. 5. Volume dependence of the neutron mass is expected to be exponentially suppressed4and is found to be consistent with zero within statistical uncertainties. The proton massincludes 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 ) − ) and are thereforesuppressed in E p ( L ) compared to E π + ( L ) by m π + /M p (cid:28) . The proton energy is foundto 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 . , M p = 1171(7)(34) MeV , (15) aM n = 0 . , M n = 1152(7)(34) MeV , (16)where the uncertainties are as defined for the π + after Eq. (13). Combining these resultsgives M p − M n = 20(10) MeV. This result is approximately ten times larger than the QEDcontribution to the proton-neutron mass difference at the physical values of the quark massesand α [1–3], which is consistent with the expected linear dependence of the proton-neutronmass difference on α , given the value α (cid:39) . that is used here. III. FINITE-VOLUME NON-RELATIVISTIC QED
Hadronic EFT results relating FV energy levels to the LECs appearing in EFT La-grangians are useful for interpreting LQCD+QED L results, as evidenced by the use of fitsto 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 neededto extract hadronic scattering information from LQCD+QED L results for multi-hadron FVenergy levels, but EFT for charged multi-hadron systems is complicated by the presence ofCoulomb interactions that are nonperturbative for hadron pairs with sufficiently small rela-tive momentum as discussed below. Further complications arise for non-relativistic EFTs forQCD+QED L because of the nonlocality inherent in zero-mode subtraction. However, non-relativistic EFTs have the advantage that FV energy shifts for systems of arbitrary particlenumber can be computed more simply than in relativistic EFTs, where the relation betweenscattering parameters and FV energy levels is not yet known for systems of four or morehadrons. This work therefore pursues the application of non-relativistic EFT to chargedmulti-hadron FV systems, and this section develops the formalism necessary for studyingmulti-hadron systems in NRQED L including nonlocal effects arising from zero-mode sub-traction. A. Finite-volume formalism for two charged hadrons
Interactions of two electrically neutral particles with mass M and relative momentum p are described at low-energies by a scattering phase shift δ ( p ) , where p = | p | . The phase shiftis an analytic function of the center-of-mass energy E ∗ = 2 (cid:112) p + M for energies below boththe t -channel cut and the s -channel inelastic particle production threshold. For p (cid:28) M , theenergy can be described by the non-relativistic expansion E ∗ = 2 M + p /M + . . . , and belowthe t -channel cut and inelastic threshold the phase shift admits a convergent effective rangeexpansion p cot δ ( p ) = − a + r p + . . . , where a is the scattering length (not to be confusedwith the lattice spacing a appearing elsewhere) and r is the effective range. For neutralparticles, 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.5Interactions of electrically charged particles are complicated by the fact that the t -channelcut associated with one-photon-exchange and the inelastic photon production threshold startat p = 0 . This leads to p = 0 singularities in contributions to the scattering amplitudefrom 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 diagramscorresponds to a perturbative expansion in powers of η ( p ) ≡ αv ( p ) = αM p , (17)where v is the relative velocity of the two charged particles and p (cid:28) M is assumed through-out this section. For η ( p ) (cid:38) , QED becomes nonperturbative and Coulomb ladder diagramsmust be resummed to all orders in α . As shown in nonrelativistic quantum mechanics byBethe [64], and in EFT by Kong and Ravndal [65, 66], the resummed scattering amplitudeis nonanalytic in η and the effective range expansion is modified as (cid:18) πηe πη − (cid:19) p cot δ ( p ) + 2 pη [ Re ψ ( iη ) − ln( η )] = − a C + 12 r C p + . . . , (18)where ψ ( x ) = Γ (cid:48) ( x ) / Γ( x ) is the digamma function, a C is the Coulomb-corrected scatteringlength, and r C is the effective range of the charged particle system. To leading order in α , andat all orders in η and in the four-particle contact interactions as in pionless EFT, the chargedparticle effective range is unaffected by QED interactions and r C = r [66]. Conversely, thecontact interaction associated with the scattering length must be renormalized to absorbdivergences from Coulomb ladder diagrams. In the MS scheme, the running coupling thatwould be identified with the scattering length in the absence of QED is related to theCoulomb-corrected scattering length by [66] a MS ( µ ) = 1 a C + αM (cid:20) ln (cid:18) µ √ παM (cid:19) + 1 − γ E − µr (cid:21) , (19)and can be understood as the physical scattering length with Coulomb effects from infraredlength scales > /µ removed.The FV two-particle energy spectrum can be related to the scattering phase shift innon-relativistic quantum mechanics [67] and in quantum field theory [39, 40]. In a finitespatial volume of size L with PBCs, the system exhibits reduced spatial symmetries char-acterized by covariance under the cubic group, and the momentum carried by a free particleis quantized as p = 2 π n /L with n ∈ Z . In QED L , zero-mode subtraction implies thatthe momentum carried by a photon propagator is restricted to p ≥ π/L . The expansionparameter η ( p ) in Eq. (17) is therefore restricted to η ( p ) ≤ η L ≡ αM L π (20)in QED L Coulomb ladder diagrams. For sufficiently small volumes, η ≤ η L (cid:28) andCoulomb photon exchange can be treated perturbatively.The O ( α ) quantization condition relating the s -wave scattering amplitude to the two-particle spectrum in the A +1 representation of the cubic group in non-relativistic EFT in the6approximation of negligible partial wave mixing was derived in Ref. [42], (cid:18) πηe πη − (cid:19) p cot δ ( p ) + 2 pη [ Re ψ ( iη ) − ln( η )] = 1 πL S C ( p ) + αM (cid:20) ln (cid:18) παM L (cid:19) − γ E (cid:21) , (21)where η ( p ) is defined in Eq. (17) and S C ( x ) = S ( x ) − αM L π S ( x ) + αM a C r C π L I S ( x ) + O ( α ) , (22a) S ( x ) = lim Λ n →∞ | n | < Λ n (cid:88) n ∈ Z | n | − x − π Λ n , (22b) S ( x ) = lim Λ n →∞ | n | < Λ n (cid:88) n ∈ Z (cid:88) m ∈ Z \{ n } | n | − x | m | − x | n − m | − π ln Λ n , (22c)where I ≈ − . is a geometric constant whose evaluation is detailed in Refs. [39] andsums 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 FVenergy levels below three-particle thresholds can be related to charged particle scatteringphase shifts by Eq. (21) and used to constrain parameterizations of the phase shift such asEq. (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 thatmust be resummed for sufficiently large volumes, Eq. (21) is necessarily valid only for anintermediate range of L . An important goal of this work is to test the range of volumesover which Eq. (21) can be reliably used to extract Coulomb-corrected scattering parametersfrom 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 ef-fective range expansion coefficients. Ref. [42] determined this expansion to O ( αM L ) and O ( a C /L ) under the assumption that αM L (cid:28) a C /L (cid:28) . Following Refs. [68, 69], thesame result can be derived in non-relativistic quantum mechanics using a Hamiltonian thatperturbatively includes the effects of relativity and allows straightforward generalization tomany-particle systems. The NRQED L Lagrangian is given by L = ψ † (cid:18) iD − D M (cid:19) ψ − (cid:18) πaM (cid:19) ( ψ † ψ ) − η ( µ )3! ( ψ † ψ ) + L ξγ + L r . (23)In this expression ψ is a non-relativistic hadron field, D µ = ∂ µ + ieQA µ where Q is theelectric charge operator, the gauge-fixed photon Lagrangian L ξγ is given in Eq. (A2), η ( µ ) is a renormalization-scale-dependent coupling associated with short-range three-body inter-actions, and four- and higher-body interactions are neglected. L r includes effective rangecontributions and relativistic corrections involving two derivatives, and is given by L r = ψ † (cid:18) D M (cid:19) ψ − (cid:16) πaM (cid:17) (cid:18) ar − M (cid:19) ( ψ † ψ )( ψ † D ψ + ψ D ψ † ) . (24)7The coefficient of the ( ψ † ψ ) operator in Eq. (24) can be fixed by demanding that thestrong interaction EFT given by replacing D µ with ∂ µ reproduces Eq. (21) with α = 0 whenboth are expanded perturbatively in powers of L − to O ( L − ) . This O ( L − ) thresholdexpansion of Lüscher’s quantization condition is given in Ref. [72] and is verified below to bereproduced by the non-relativistic EFT defined by Eqs. (23)-(24) . Operators in L r lead tocontributions to the threshold expansion suppressed by ar/L as well as relativistic effectssuppressed by O (( M L ) − ) that will be neglected in the power counting schemes discussedbelow. Additional relativistic corrections to Eq. (23) arise from photon loops and operatorswith four and more derivatives, but these give rise to FV effects suppressed by powers of O (( M L ) − ) that will be neglected as discussed below and detailed in Appendix A.Introducing Fourier transformed fields (cid:101) ψ k = 1 L / (cid:90) d x e i k · x ψ ( x ) , ψ ( x ) = 1 L / (cid:88) k e − i k · x (cid:101) ψ k , (25)the FV Hamiltonian for NRQED L is given by H = (cid:88) k (cid:101) ψ † k (cid:18) k M (cid:19) (cid:101) ψ k + H ξγ + H int , (26)where H ξγ is the gauge-fixed photon Hamiltonian, whose explicit form will not be used below,and H int = − (cid:88) k (cid:101) ψ † k (cid:18) k M (cid:19) (cid:101) ψ k + 12 L (cid:88) p (cid:48) , p , Q V ( p , p (cid:48) ) (cid:101) ψ † Q − p (cid:48) (cid:101) ψ † Q + p (cid:48) (cid:101) ψ Q − p (cid:101) ψ Q + p + η ( µ )(3!) L (cid:88) Q , p , p (cid:48) , q , q (cid:48) (cid:101) ψ † Q + p (cid:48) (cid:101) ψ † Q + q (cid:48) (cid:101) ψ † Q − p (cid:48) − q (cid:48) (cid:101) ψ Q + p (cid:101) ψ Q + q (cid:101) ψ Q − p − q , (27)where the two-body potential includes strong interaction and Coulomb terms, V ( p (cid:48) , p ) = 4 πM (cid:18) a + a (cid:18) ar − M (cid:19) ( p + p (cid:48) ) + . . . (cid:19) + 4 πα | p (cid:48) − p | (1 − δ p , p (cid:48) ) . (28)Relativistic corrections to the potential, including photon loop effects as well as nonlocaleffects of zero-mode removal, can in principle be calculated by carrying out the FV analogof higher-order matching between QED, NRQED, and pNRQED [73–75], but these termsgive rise to O ( M L ) − effects that are neglected below.Specializing to the case of identical bosons, operators are associated with the fields inEq. (26) and are defined to satisfy commutation relations [ (cid:101) ψ † p (cid:48) , (cid:101) ψ p ] = δ p , p (cid:48) . FV two-particlestates defined by | p , p (cid:105) = 1 √ (cid:101) ψ † p (cid:101) ψ † p | (cid:105) (29) The Lagrangian in Eq. (24) can also be obtained by studying the nonrelativistic limit of relativistic scalarfield theory as in Refs. [70, 71] . Refs. [68, 69] include the operators in Eq. (24), but a factor of 2 discrepancy in the coefficient of the lastterm leads to a difference in the O ( L − ) threshold expansion result. (cid:104) p , p | p , p (cid:105) = 1 . The ground-state of the two-particle system in its center-of-mass frame is | p = , p = (cid:105) . Splitting the Hamiltonian into kinetic energy terms and aninteraction term H int that will be treated perturbatively, the ground-state FV energy shiftat leading order in Rayleigh-Schrödinger perturbation theory is then given by ∆ E LO = (cid:104) , | H int | , (cid:105) = 1 L V ( , ) = 4 πaM L . (30)Standard perturbation theory techniques can be used to extended this result to higher ordersin a/L and α . At next-to-next-to-leading order (NNLO) in H int , the result is given by ∆ E NNLO,PC1 = 4 πa C M L (cid:26) − (cid:16) a C πL (cid:17) I + (cid:16) a C πL (cid:17) (cid:2) I − J (cid:3)(cid:27) + απL (cid:26) − (cid:16) a C πL (cid:17) J + (cid:16) a C πL (cid:17) (cid:26) IJ + R − K + 4 π (cid:20) ln (cid:18) παM L (cid:19) − γ E (cid:21)(cid:27) − (cid:18) αM L π (cid:19) K + (cid:16) a C πL (cid:17) (cid:18) αM L π (cid:19) (cid:2) R + J − L (cid:3) + (cid:18) αM L π (cid:19) R (cid:41) , (31)where PC1 is a label for the power counting scheme discussed below and higher-order cor-rections in a C /L ∼ r/L , as well as relativistic effects suppressed by O (( M L ) − ) , have beenneglected. All sum-integral differences I , J , K , L , and R nm which appear in this expressionare evaluated in Ref. [42], which includes an evaluation of Eq. (31) in EFT to leading orderin αM L and to higher order in a/L , except for R which is a convergent sum given by R ≈ . .Eq. (31) can be expressed as an expansion in a C /L and the FV Coulomb parameter η L defined in Eq. (20), ∆ E NNLO,PC1 = 4 πa C M L (cid:26) − (cid:16) a C πL (cid:17) I + (cid:16) a C πL (cid:17) (cid:2) I − J (cid:3)(cid:27) + 4 η L M L (cid:26) − (cid:16) a C πL (cid:17) J + (cid:16) a C πL (cid:17) (cid:2) IJ + R − K − π [ln ( η L ) + γ E ] (cid:3) − (cid:16) η L π (cid:17) K + (cid:16) a C πL (cid:17) (cid:16) η L π (cid:17) (cid:2) R + J − L (cid:3) + (cid:16) η L π (cid:17) R (cid:27) . (32)This resembles a double power series expansion in the parameters η L and a C /L , suggestingthat the NRQED L threshold expansion should provide a good approximation of FV energyshifts for PC1: η L ∼ a C L (cid:28) . (33)It is noteworthy that η L (cid:28) only requires αM L (cid:28) π and is less restrictive than thecondition αM L (cid:28) discussed in Ref. [42]. In the matching to the LQCD+QED L simulations See also Refs. [39, 40, 42, 76, 77] for more details on evaluating FV sums. η L ∼ a C /L numerically overestimates the size of QEDeffects, particularly on the smaller volumes studied, and the alternative power counting PC2: η / L ∼ a C L (cid:28) , (34)will also be used in fits to LQCD+QED L results. Higher-order results for the thresholdexpansion for short-range contact interactions without QED [68, 69, 72, 78] can be used toextend Eq. (32) from NNLO in Eq. (33) to next-to-next-to-next-to-leading-order (N LO) inEq. (34), ∆ E N LO,PC2 = 4 πa C M L (cid:40) − (cid:18) a C πL (cid:19) I + (cid:18) a C πL (cid:19) (cid:2) I − J (cid:3) − (cid:18) a C πL (cid:19) (cid:2) I − IJ + K (cid:3)(cid:41) + 4 η L M L (cid:26) − (cid:18) a C πL (cid:19) J − (cid:16) η L π (cid:17) K + (cid:18) a C πL (cid:19) (cid:2) IJ + R − K − π [ln ( η L ) + γ E ] (cid:3)(cid:41) . (35)The parameter a C ( L ) is equal to a C plus a /L suppressed correction arising from theinteraction terms in Eq. (24), a C = a C ( L ) − πa C ( L ) L (cid:18) a C ( L ) r − M (cid:19) , (36)where ( M L ) − suppressed effects are shown for completeness and lead to agreement with the O ( L − ) strong interaction relativistic threshold expansion of Ref. [72] after taking α = 0 .In principle, LQCD+QED L results for π + π + FV energy shifts on multiple lattice volumescould be used to extract both a C and r by constraining the O ( L − ) difference between a C and a C . In the LQCD+QED L calculations discussed below, ( a C ( L ) − a C ) /a C can be estimatedat LO in chiral perturbation theory ( χ PT) to be and for the L/a = 48 and
L/a = 32 lattice volumes respectively. To see whether a C − a C can be reliably determined, this estimatemust be compared with an estimate of relativistic effects neglected in Eq. (35), which asdiscussed in Sec. III C below modify the dominant QED FV effects by O (( M L ) − ) . 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 − J ( η L /π ) a C , and radiation photon effects can be estimatedto lead to a O (( M L ) − ) suppressed shift in ( a C ( L ) − a C ) /a C of order ∼ J ( η L /π ) / ( M L ) = J α/ (2 π ) ∼ for α (cid:39) . . Relativistic effects are therefore comparable to a C − a C for π + π + systems and prevent the effective range contribution to a C ( L ) from being disentangledfrom other FV effects neglected in Eq. (35). Therefore, a C ( L ) − a C will be neglected whenfitting π + π + LQCD+QED L results to Eq. (35). The description of FV effects on charged hadron masses in QED L as a dual expansion in α and ( M L ) − and the possibility of using alternative power countings in LQCD+QED L calculations depending on thevalues of these parameters is explored in Ref. [35]. The operators in Eq. (24) lead to additional effective range and relativistic corrections to the right-hand-side of Eq. (36), namely an additional term of the form − a C rη L I / ( M L ) − η L J / (4 π M L ) , but theseare higher-order than the effective range term in Eq. (36) according to the power counting of Eq. (34). C. Zero-mode effects
The derivation of Eq. (21) in Ref. [42] and the form of the NRQED L Lagrangian inEq. (23) assume that charged particle NRQED L contact interactions in FV are equal totheir infinite-volume counterparts up to exponentially suppressed corrections. More re-cently, however, it has been shown in Refs. [2, 29, 31, 35, 37, 38] that this assumption isviolated in the single-particle sector of NRQED L because of the inherent nonlocality of zero-mode subtraction. NRQED L parameters can be obtained by calculating masses, scatteringamplitudes, or other observables in both NRQED L and QED L and tuning the parametersof the NRQED L Lagrangian to reproduce QED L results. In QED L , on-shell photon ex-change leads to power-law FV effects on charged particle masses suppressed by powers of α and / ( M L ) that have been studied in Refs [2, 29, 31, 35, 37, 38]. The O ( α/ ( M L )) and O ( α/ ( M L ) ) corrections are independent of the structure of the charged particle andare described by one-loop diagrams in both QED L and NRQED L . At order O ( α/ ( M L ) ) ,structure-dependent effects involving magnetic moments and charge radii arise. Nonlocaleffects from zero-mode subtraction also enter at order O ( α/ ( M L ) ) because zero-mode sub-traction leads to power-law FV effects from off-shell antiparticle modes in QED L that arenot reproduced by NRQED L loop diagrams. These effects can be included in NRQED L byadjusting the Lagrangian to include additional particle-antiparticle interactions [31, 38], ormore simply by adjusting the coefficients of mass operators in the NRQED L Lagrangianby factors proportional to α/ ( M L ) [37]. For charged scalars, although not for chargedfermions, these effects vanish in the charged particle rest frame [37]. This non-decouplingof antiparticle modes is a consequence of the nonlocality of NRQED L , and for EFTs withbreakdown scale Λ generically produces O ( α/ (Λ L ) ) corrections to LECs [37].Nonlocal effects from zero-mode subtraction could lead to power-law FV effects that mod-ify four-hadron contact interaction couplings proportional to a C in NRQED L . Consideringthe O ( M L ) enhancement of FV effects associated with Coulomb ladder diagrams, it is neces-sary to analyze nonlocal effects of zero-mode subtraction on a C in order to determine whetherEq. (21) is modified within the order of approximation considered. The effects of zero-modesubtraction on a C can be determined by matching any QED L and NRQED L correlationfunctions sensitive to four-hadron contact interactions. One-particle FV self-energies onlyreceive contributions from four-hadron contact interactions in O ( α ) diagrams containingclosed loops of particle-antiparticle pairs. Two-particle FV Green’s functions receive con-tributions from four-hadron contact interactions at O ( α ) , and it is convenient to calculatenonlocal FV effects on four-hadron contact interactions in NRQED L by directly matchingtwo-particle FV Green’s functions in QED L and NRQED L . This matching is detailed inAppendix A and summarized below.The O ( M L ) enhancement of Coulomb ladder diagrams arises when one intermediateparticle propagator is placed on-shell and another intermediate particle propagator is nearlyon-shell with virtuality k /M , as compared with virtuality k when a photon propagatoris 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 row1 FIG. 6: The strong-interaction and Coulomb scattering diagrams contributing to the two-bodyFV energy shift in NRQED L . The top-left section shows the LO diagram . The top-right sectionshows the NLO diagrams. The bottom section shows the NNLO diagrams in the power countingof Eq. (33). Diagrams that vanish because of zero-mode subtraction, including the tree-level one-photon-exchange diagram, are not shown.FIG. 7: Radiation photon diagrams making power-suppressed contributions to the two-body FVenergy shift in NRQED L . The left and right sections show NLO and NNLO diagrams in the powercounting of Eq. (33), respectively, which lead to FV effects suppressed by O (cid:0) ( M L ) − (cid:1) . in Fig. 7). In NRQED L , the box diagram involves the integral i (cid:90) dk π (cid:18) k − k / M + i(cid:15) (cid:19) (cid:18) − k − k / M + i(cid:15) (cid:19) (cid:18) k ) − k + i(cid:15) (cid:19) = − (cid:18) Mk (cid:19) (cid:18) k / M ) − k (cid:19) + ddk (cid:34)(cid:18) k − | k | (cid:19) (cid:18) k − k / M (cid:19) (cid:18) − k − k / M (cid:19)(cid:35) k →− k = − Mk (cid:2) O ( k /M ) (cid:3) − k (cid:2) O ( k /M ) (cid:3) . (37)In this expression, the first contribution involving the particle pole places the second particlepropagator nearly on-shell with kinetic energy k / M . The second contribution from thephoton double pole gives both particle propagators off-shell kinetic energies k . When FV2 FIG. 8: Radiation photon diagrams, including the jellyfish diagram on the forth row, that appearat NNLO in the power counting of Eq. (33) and lead to FV effects suppressed by O (cid:0) ( M L ) − (cid:1) . effects are computed, k is replaced by the quantized values πn/L with n ∈ Z , and (afteradding all necessary UV counterterms) amplitude suppression by powers of k/M impliessuppression of FV effects by the corresponding power of ( M L ) − . The NRQED L crossed-box diagram involves the integral i (cid:90) dk π (cid:18) k − k / M + i(cid:15) (cid:19) (cid:18) k ) − k + i(cid:15) (cid:19) = ddk (cid:34)(cid:18) k − | k | (cid:19) (cid:18) k − k / M (cid:19) (cid:35) k →− k = 34 k (cid:2) O ( k /M ) (cid:3) , (38)where only the photon pole contributes and leads to particle propagators with off-shell kineticenergies k . FV effects associated with the crossed-box diagram are therefore suppressed by O (( M L ) − ) compared to the dominant contribution of the box diagram.The O ( M L ) power enhancement of the box diagram is only present in the diagramsinvolving repeated s -channel Coulomb interactions shown in Fig. 6 but not in the diagramsin Figs. 7-8 or other diagrams involving particle-antiparticle pair creation that vanish non-relativistically . Furthermore, the O ( M L ) enhancement only occurs when the intermediate-state charged particles are both nearly on-shell. As detailed in Appendix A, power-law FVeffects in QED L that are not reproduced by loop diagrams with the leading order NRQED L Lagrangian arise from antiparticle poles where intermediate states have a large virtuality oforder M . These do not receive the O ( M L ) enhancement of particle pole contributions, andzero-mode effects are found to be suppressed by O ( α/ ( M L ) ) . Matching between QED L and NRQED L is explicitly performed for charged scalars in Appendix A, and zero-mode More details on the power-counting of diagrams involving FV photon exchange are given in Ref. [42] anddiscussions of analogous power counting arguments for radiation pions are given in Ref. [79]. FIG. 9: The strong-interaction and Coulomb diagrams contributing to three- and higher-body FVenergy shifts in NRQED L at leading order in ( M L ) − and NNLO in short-range and Coulomb inter-actions in the power counting of Eq. (33). Diagrams that vanish because of zero-mode subtraction,partially disconnected diagrams involving pairs of two-body interactions among four and more par-ticles, and diagrams that involve on-shell internal propagators and vanish in Rayleigh-Schrödingerperturbation theory, are not shown. effects are found to modify the four-scalar coupling in NRQED L at order O ( α/ ( M L ) ) forboosted systems but not to modify the coupling for scalars at rest at this order. This showsthat Eq. (21) is valid for charged scalars in NRQED L up to O ( α ) effects and O ( M L ) − relativistic effects. The vanishing of these effects in the center-of-mass frame arises fromcancellations due to the specific form of the scalar-photon vertex functions. It is possiblethat for charged fermions nonlocal effects also arise at O ( α/ ( M L ) ) in the charged fermionrest frame. The O ( α/ ( M L ) ) suppression factor is consistent with the physical argumentsof Ref. [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 lawis satisfied for the FV system [25], and high-energy modes that have been integrated outof 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 expan-sion for FV effects on systems of n non-relativistic particles using Rayleigh-Schrödingerperturbation theory. Unit-normalized many-particle states are given by | p , . . . , p n (cid:105) = 1 √ n (cid:101) ψ † p × . . . × (cid:101) ψ † p | (cid:105) , (39)and the leading order FV energy shift for the ground state of n identical bosons in thecenter-of-mass frame is given by (cid:104) , . . . | H int | , . . . (cid:105) = 1 L (cid:18) n (cid:19) V ( , ) + 1 L (cid:18) n (cid:19) η ( µ ) . (40)4 FIG. 10: Radiation photon diagrams with negligible contributions to three- and higher-body FVenergy shifts in NRQED L at NNLO. The left and center diagrams are suppressed by O (cid:0) ( M L ) − (cid:1) ,while the right diagram is suppressed by O (cid:0) ( M L ) − (cid:1) . As in Fig. 9, partially disconnected diagramsand diagrams that vanish because of zero-mode subtraction or kinematical constraints are notshown. Working to NNLO in the power counting of Eq. (33), the energy shift of an n -hadron state isequal to (cid:0) n (cid:1) times the two-body energy shift of Eq. (32), plus additional contributions frominduced three-body and four-body forces shown in Fig. 9. Additional diagrams associatedwith radiation photon exchange shown in Fig. 10 are suppressed by O (( M L ) − ) . Theresulting threshold expansion for the FV energy shift for a system of n like-charged hadronsat rest is given by ∆ E NNLO,PC1 n = 4 πa C M L (cid:18) n (cid:19) (cid:26) − (cid:16) a C πL (cid:17) I + (cid:16) a C πL (cid:17) (cid:2) I + (2 n − J (cid:3)(cid:27) + 4 η L M L (cid:18) n (cid:19) (cid:26) − (cid:16) a C πL (cid:17) J − (cid:16) η L π (cid:17) K + (cid:16) η L π (cid:17) R + (cid:16) a C πL (cid:17) (cid:2) IJ + R + (4 n − K − π [ln ( η L ) + γ E ] (cid:3) + (cid:16) a C πL (cid:17) (cid:16) η L π (cid:17) (cid:2) R + J + (2 n − L (cid:3)(cid:111) , (41)where omitted terms are: quartic or higher in η L ∼ a C /L ; relativistic effects suppressedby O (( M L ) − ) ; or three-body contact interactions where η ∼ a C /M is assumed so that O ( L − ) terms of the strong interaction threshold expansion appear at the same order.At N LO in the power counting of Eq. (34) there are contributions from short-distancethree-body interactions well as the induced few-body interactions discussed above. Theseintroduce a new free parameter η ( µ ) that parameterizes the strength of six-particle opera-tors in the NRQED L Hamiltonian, Eq. (26). Combining Eq. (41) with the N LO results forthe QCD threshold expansion from Ref. [68], the N LO FV energy shift for a system of n ∆ E N LO,PC2 n = 4 πa C M L (cid:18) n (cid:19) (cid:40) − (cid:18) a C πL (cid:19) I + (cid:18) a C πL (cid:19) (cid:2) I + (2 n − J (cid:3) − (cid:18) a C πL (cid:19) (cid:2) I + (2 n − IJ + (5 n − n + 63) K (cid:3)(cid:41) + 4 η L M L (cid:18) n (cid:19) (cid:26) − (cid:18) a C πL (cid:19) J − (cid:16) η L π (cid:17) K + (cid:18) a C πL (cid:19) (cid:2) IJ + R + (4 n − K − π [ln ( η L ) + γ E ] (cid:3)(cid:41) + (cid:18) n (cid:19) L (cid:20) η ( µ ) + 64 πa C M (3 √ − π ) ln( µL ) − a C π M S MS (cid:21) . (42)The non-relativistic EFT three-body coupling η ( µ ) is renormalization scheme- and scale-dependent. The scale-dependence of η ( µ ) cancels the explicit ln( µL ) scale-dependenceshown in Eq. (42), and the scheme-dependence is compensated by scheme-dependence inthe finite term S MS . This scale-dependence arises from the ambiguity in separating short-distance three-body interactions described by contact operators in NRQED L from long-distance two-body rescattering effects. In relativistic theories this ambiguity does not arise,and the particle mass plays the role of the scale µ in relativistic descriptions of the L − ln L term in the three-body threshold expansion derived for generic relativistic field theories inRef. [78]. Equating the O ( L − ) term in the relativistic threshold expansion of the three-particle threshold amplitude M , th in Ref. [78] with the corresponding O ( L − ) nonrelativisticthreshold expansion provides a relation between the non-relativistic coupling η ( µ ) and thescale-independent M , th in the absence of QED, η ( µ ) = − M , th M + 64 πa M (3 √ − π ) ln (cid:18) M πµ (cid:19) + 48 a π M + 48 a π rM + 12 a π M S + 768 π a M C , (43)where a is the scattering length for a neutral two-particle system, C = − . is a FVsum evaluated in Ref. [78], and S = 571 . is related to S MS , evaluated in Ref. [69], andto other FV sums from Ref. [78] by S = C F + C + C + 8 S MS . QED effects will modifyEq. (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 interactionparameters 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 andother hadronic interaction parameters at the values of the quark masses and α used here.6 a ∆ E nK a ∆ E nπ + n L/a = 32 L/a = 48
L/a = 32
L/a = 48 n ∈ { , . . . , } neutral K mesons and charged π + mesons for lattice volumes with L/a ∈ { , } . Results are determined by taking correlateddifferences between LQCD+QED L ground-state energies during the fit range sampling proceduredescribed in Appendix B. A. Charged meson scattering
The LQCD+QED L results for the FV spectrum results in Table II can be used to constrainthe low-energy EFTs for charged and neutral meson interactions. It is convenient to focuson results for the FV energy shifts ∆ E nM ( L ) = E nM ( L ) − nE M ( L ) , M ∈ { K , π + } . (44)Results for correlated differences between n -particle ground-state energies and n times theone-particle ground-state energy, as defined in Eq. (44), are more precise than n -particleenergies alone. Furthermore, this subtraction nonperturbatively removes single-particle FVeffects from n -meson FV energy shift results. For multi- π + systems, LQCD+QED L re-sults 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 systems, LQCD+QED L re-sults can be identified with the same EFT results after setting α to zero. In the numericalLQCD+QED L calculation, FV energy shifts are computed in a correlated manner usingbootstrap resampling as detailed in Appendix B, and the results are shown in Table IV. Toaccess QED-specific effects, the differences of these differences between the n charged pionsand n neutral kaon systems are also computed similarly. The double subtraction suppressesany strong isospin breaking effects arising from mistuning of the quark masses for differentcharge quarks. In the numerical LQCD+QED L calculation, correlated differences of FVenergy shifts are computed using bootstrap resampling as detailed in Appendix B, and theresults are shown in Table V.Results for the two-particle FV energy shifts ∆ E π + π + and ∆ E K K are shown in Fig. 11.Both energy shifts are clearly resolved from zero with relative uncertainties in the rangeof − for both volumes, although ∆ E π + π + and ∆ E K K on a given volume areindistinguishable. The small magnitude of QED effects on ∆ E π + π + might appear sur-7 a ∆ E nπ + − a ∆ E nK n L/a = 32 L/a = 48 n ∈ { , . . . , } charged π + and neutral K mesons for lattice volumes with L/a ∈ { , } . Results are obtained by taking correlateddifferences between fitted energies as in Table IV. prising because αm π L ∼ . for this volume, but as discussed in Sec. III the appropri-ate FV analog of the Coulomb expansion parameter is η L = αm π L/ (4 π ) ∼ . for the L/a = 48 volume. Eqs. (32)-(35) therefore predict that in addition to differences arisingfrom a π + π + C (cid:54) = a K K , NLO corrections from Coulomb photon exchange modify the LO FVenergy shift by ∼ on the L/a = 48 lattice volume, which is not expected to be distin-guishable given the statistical uncertainties on ∆ E π + π + and ∆ E K K . This expectation isconsistent with LQCD+QED L results, as shown in Fig. 11.The scattering lengths a π + π + C and a K K can be extracted from a combined fit to theresults for ∆ E π + π + and ∆ E K K shown in Table IV and the results for the precisely deter-mined correlated differences ∆ E π + π + − ∆ E K K shown in Table V. Fitting to Eq. (32), asis appropriate for the power counting PC1 of Eq. (33), gives the results,NNLO, PC1 : a K K m K = 0 . , a π + π + C m K = 0 . , (45)where the common scale m K has been included for both π + π + and K K to facilitatecomparison of a K K and a π + π + C . The lowest-order QED effect on ∆ E π + π + from Coulombphoton exchange in Eq. (32) decreases ∆ E π + π + compared to the FV energy shift for neutralparticles. The scattering length results in Eq. (45) show that this effect from Coulomb Coulomb photon exchange leads to a decrease in the energy of a system of π + mesons in QED L becauseof zero-mode subtraction. Physically, the energy decrease can be understood as arising from attractionbetween the charged particle system and the uniform background of opposite charge associated with zero-mode subtraction [25]. Formally, zero-mode subtraction removes the LO one-photon-exchange diagramsassociated with repulsion between charged particles. The dominant QED contribution therefore arisesat NLO and necessarily lowers the ground-state energy since it appears at second order in perturbationtheory. □ □○ ○ □ ○ FIG. 11: The blue and red points show the LQCD+QED L results from Table IV for the K K and π + π + FV energy shifts for the
L/a ∈ { , } 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 ∈{ , } results to Eq. (32). The blue band shows the prediction of Eq. (32) with α = 0 usingthe best-fit result for a K K in Eq. (45). As in Fig. 4, the widths of the bands correspond to67% confidence intervals estimated using bootstrap resampling. A small horizontal offset is appliedsymmetrically to π + and K results. photon exchange competes with additional QED effects that lead to a π + π + C > a K K . Fittingto Eq. (35), as is appropriate for the power counting PC2 of Eq. (34), gives consistent results,N LO, PC2 : a K K m K = 0 . , a π + π + C m K = 0 . , (46)demonstrating that the fit is not overly sensitive to higher-order terms absent in one or theother power counting.Both Eq. (32) and Eq. (35) neglect relativistic effects from radiation photon exchangeleading to O (( M L ) − ) FV effects. These effects are estimated in Sec. III B to lead to a shiftin a π + π + C of order ∼ , which is smaller than the − 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 proportionalto αM p L that are enhanced compared with those in the π + π + case discussed above. For bothlattice volumes αM p L > , and according to the scaling estimates of Ref. [42] Coulomb effectsshould be nonperturbative. However, η pL = αM p L/ (4 π ) ∼ . for the L/a = 48 latticevolume and the NRQED L results in Sec. III expressed as power series in η pL show signs ofconvergence. Examining the QED L contributions proportional to a ppC /L in Eq. (32), the NLOmixed strong-Coulomb contribution is suppressed compared to the LO strong contributionby ( η pL /π )2 J ∼ . , while the corresponding NNLO contribution is suppressed compared9 a ∆ E b b L/a = 32 L/a = 48 pp np ( S ) nn np ( S ) He -0.011(96) 0.038(56) H 0.015(75) 0.080(45)TABLE VI: FV energy shift results for systems of n ∈ { , , } protons and neutrons determined byfitting Eq. (9) to LQCD+QED L Euclidean correlation function results as described in Appendix B. □ □□○ ○○△ △△▽ ▽▽ □ ○ △ ▽ - - FIG. 12: Points show the LQCD+QED L results from Table VI for FV energy shifts determinedfrom the correlated differences of two-nucleon and one-nucleon ground-state energies for the neutron-neutron and proton-proton FV energy shifts for the L/a ∈ { , } volumes. A smaller (larger)horizontal offset is applied symmetrically to pp and nn ( np ( S ) and np ( S ) ) results. to the LO contribution by ( η pL /π )(2 R + J − L ) ∼ . . This suggests that Coulombeffects 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–87]. Thequantization condition in Eq. (21), which neglects O (( η pL ) ) perturbative Coulomb effects and ( M L ) − relativistic effects but is nonperturbative in strong interaction effects, is thereforeneeded to relate the pp FV energy shifts to the infinite volume pp phase shift and determine a ppC .The two-nucleon isospin I = 1 systems pp , nn , np ( S ) , as well as the deuteron, arestudied on both lattice volumes. Relatively clean signals are seen for each system, and theirground-state energies are determined with total (statistical plus fitting systematic) uncer-tainties at the level as shown in Table III. FV energy-level shifts are determined fromcombined analyses of the two-nucleon and single-nucleon correlation functions as describedin Appendix B, and fit results for all systems are shown in Appendix B 4. As shown in Ta-0ble. VI and Fig. 12, the statistical precision of this calculation is insufficient to resolve eitherthe 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 with-out QED at these quark masses at a 95% confidence level for the pp , nn , and np systems onthe L/a = 32 lattice volume would require statistical ensembles approximately ∼ − times larger than the one used here, as estimated by extrapolating the uncertainties of thefour different two-nucleon systems considered assuming / √ N scaling of uncertainties. De-termination of a ppC through the QED L quantization condition of Eq. (21) is therefore left tofuture work.The two I = 1 / three-nucleon systems He and H are also investigated, and theirground-state energies are given in Table III. Correlated differences between three-nucleonground-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 Heand H are not precise enough to allow FV effects to be reliably determined. Precision in thethree-nucleon sector is significantly worse than in the two-nucleon sector, as expected. Theabsolute size of FV energy shifts is also expected to be larger for three-nucleon systems thantwo-nucleon systems, and for instance resolving an O (50 MeV ) FV energy shift at a 95%confidence level for the He and H systems on the
L/a = 32 lattice volume would require astatistical ensemble approximately ∼ times larger than the one considered here, basedon an extrapolation analogous to that described for the two-nucleon case.Future high-precision LQCD+QED L calculations of these multi-nucleon systems will pro-vide insight into QED effects on nucleon-nucleon and three-nucleon interactions through adetermination of the He - H binding-energy difference and its decomposition into QEDand 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-noisedegradation with increasing particle number and can be used to study QED in the regimewhere the charge Z ∼ /α . In particular, the correlation functions for systems with n ≤ π + mesons described in Sec. II can be used to study systems with Zα ≤ . , reaching acharge density of n/L ∼ . fm − . The dominant strong interactions and QED effects onmany-particle FV energy shifts in Eqs. (41)-(42) both scale with n , and both ∆ E nπ + and ∆ E nK can be extracted for larger n with better relative precision than from the n = 2 casediscussed in Sec. IV A.Results for ∆ E nπ + , ∆ E nK , and the correlated differences (cid:0) ∆ E nπ + − ∆ E nK (cid:1) for n ∈{ , . . . , } are shown in Tables II-V. QED effects leading to non-zero (cid:0) ∆ E nπ + − ∆ E nK (cid:1) can be resolved to better than σ on the L/a = 32 lattice volume for ≤ n ≤ , and on the L/a = 48 lattice volume for n ≥ . These 48 energy levels andcorrelated differences can be used to constrain the low-energy interaction parameters (cid:110) a π + π + C , a K K , η π + π + π + ( m K ) , η K K K ( m K ) (cid:111) appearing in Eqs. (41)-(42). Fits to the The renormalization scale used to evaluate η ( µ ) should be chosen close to the “high” energy scale whereNRQED L is matched to QED L in order to avoid large logarithms that can worsen EFT convergence. Forsimplicity, µ = m K is used as the renormalization scale for η throughout this work. O ( η L ) Coulomb effects but ne-glects three-body forces, underpredict LQCD+QED L energy-shift results for n (cid:38) mesonsystems on both lattice volumes and obtain a minimum χ /N dof ∼ . . The N LO expres-sion Eq. (42) in PC2 includes additional free parameters related to three-body forces notpresent 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 χ /N dof ∼ . .Fit results using the N LO expression Eq. (42) and uncertainties computed using bootstrapresampling of the global fitting procedure are compared to LQCD+QED L results for the π + and K energy shifts in Figs. 13-14. The results for the meson scattering lengths areconsistent with, but more precise than, the results obtained from the two-meson FV energyshifts alone,N LO, PC2 : a K K m K = 0 . , a π + π + C m K = 0 . . (47)It is noteworthy that results with Zα ≥ can be fit by the NRQED L formula given inEq. (42) without additional modifications to account for relativistic QED effects or additionalnonperturbative effects. Some tensions between LQCD+QED L results and N LO NRQED L fits can be observed for n (cid:38) meson systems on the L/a = 48 lattice volume in Fig. 13;however, since these tensions are more significant for multi- K than multi- π + systems theyare unlikely to be signals of nonperturbative QED effects and might result from correlationsbetween LQCD+QED L results with different n not accounted for in the fitting procedureemployed here.The three-pion scattering amplitude was calculated at LO in χ PT in Ref. [88] and is givenby M , th = 108 m K /f K (with conventions for the LO Lagrangian such that f π + ∼ MeV). This can be combined with Eq. (43) to provide a χ PT prediction for the non-relativistic K contact interaction, η K K K ( µ ) = 34 m K f K + m K π f K (3 √ − π ) ln (cid:18) m K πµ (cid:19) + 3 m K π f K S , th + 3 m K f K C , (48)where S , th and C are constants defined below Eq. (43), and to obtain a result entirely interms of m K and f K , the LO χ PT relations [89] a K K m K = m K πf K , r K K a K K m K = 3 , (49) The SU (2) χ PT result for M , th given in Ref. [88] is valid for tree-level K K scattering after rein-terpreting the SU (2) isospin χ PT Lagrangian in terms of the SU (2) V -spin doublet ( π + , K ) . Moreformally, the V -spin analog of G -parity acts as ( π + , K ) → ( K , − π ) and ( K , − π ) → ( − π + , K ) andtherefore relates the three-body contact operators of the SU (2) isospin and SU (2) V -spin Lagrangians by ( π + π − ) → ( K K ) . For the SU (3) flavor symmetric quark mass scheme used here, V -spin is an exactsymmetry of the leading order χ PT Lagrangian broken only by QED corrections at higher orders. ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○□ □ □ □ □ □ □ □ □ □ □ □ ○ (a) Multi-meson FV energy shifts for the L/a = 32 lattice volume. ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○□ □ □ □ □ □ □ □ □ □ □ □ ○ (b) Multi-meson FV energy shifts for the
L/a = 48 lattice volume.
FIG. 13: Points show the LQCD+QED L results from Table IV for the FV energy shifts ofmulti- K 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 best-fit parameters (cid:110) a π + π + C , a K K , η π + π + π + ( m K ) , η K K K ( m K ) (cid:111) obtained from a global fit to the L/a ∈ { , } results for n ∈ { , . . . , } mesons in Tables IV-V as described in the main text. Asmall horizontal offset is applied symmetrically to π + and K results. have been used. The first of these relations also allows f K and therefore η K K K ( µ ) to bepredicted numerically at LO in χ PT for the quark masses used in this work. Inserting theLQCD+QED L results for m K a K K from Eq. (47) into the first relation in Eq. (49) providesa prediction for f K at the parameters of this LQCD+QED L calculation that is valid at LO3 ○ ○ ○ ○ ○ ○ ○ ○ ○ - - - (a) π + and K energy shift differences for the L/a = 32 lattice volume. ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ - - - - (b) π + and K energy shift differences for the L/a = 48 lattice volume.
FIG. 14: Points show the LQCD+QED L results from Table V for the FV energy-shift differencesbetween multi- π + and multi- K systems as a function of meson number n on both lattice volumes.The shaded band shows the 67% bootstrap confidence interval for the corresponding predictionof Eq. (42) for the best-fit parameters (cid:110) a π + π + C , a K K , η π + π + π + ( m K ) , η K K K ( m K ) (cid:111) obtainedfrom a global fit to the L/a ∈ { , } results for n ∈ { , . . . , } mesons, given in Tables IV-V, asdescribed in the main text. in χ PT: af K = 0 . , f K = 139(10)(4) MeV . (50)Here, the first uncertainty is statistical and the second uncertainty is from the uncertainty inthe lattice spacing. In this calculation m K = 404(1)(12) MeV is between the physical pion4 □□ □□□□□□
FIG. 15: The blue and red points show the best-fit values and 67% bootstrap confidence intervalsfor the three-body interaction parameters for neutral K K K and charged π + π + π + systems re-spectively, including a common normalization factor of f K m K to obtain a dimensionless quantity.The green point shows the LQCD result of Ref. [52], which was obtained in a calculation using apseudoscalar meson mass of 352 MeV similar to m K here. The black point shows the LO χ PTprediction of Eq. (48) multiplied by the same normalization factor. and kaon masses; this can be compared with f π + ∼ MeV and f K ∼ MeV extractedfrom experiments [90]. Inserting this result for f K , and the result for a K K in Eq. (47),into Eq. (48) then gives the numerical result m K f K η K K K ( m K ) = 0 . . This resultis valid at LO in χ PT and can be compared to LQCD+QED L results for η K K K ( m K ) combined with the result for f K 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 η K K K ( m K ) m K f K = 2 . , η π + π + π + ( m K ) m K f K = 0 . . (51)The result of this work for η K K K ( m K ) is consistent within σ uncertainties with theLQCD results of Ref. [52], which were obtained in a calculation using quark masses similarto those used in this work, corresponding to a pion mass of 352 MeV, and extracted η byfitting to the same O ( L − ) threshold expansion as used here for multi- K systems. Thecorresponding result for η π + π + π + ( m K ) is about σ smaller than η K K K ( m K ) , althoughit is also consistent within σ uncertainties with the LQCD results of Ref. [52]. Resultsfor η K K K ( m K ) and η π + π + π + ( m K ) as well as comparisons to LO χ PT and to the LQCDresults of Ref. [52] are shown in Fig. 15.Differences between η π + π + π + ( m K ) and η K K K ( m K ) might arise from QED effectson M , th beyond LO in χ PT, or from QED effects on the matching between M π + π + π +3 , th and η π + π + π + ( m K ) . Differences between the extracted η π + π + π + ( m K ) and η K K K ( m K ) O (( M L ) − ) relativistic effects estimated inSec. III B to modify FV energy shifts by ∼ . This estimated shift is larger or compa-rable to the statistical uncertainties on the three-body energy shift for all nπ + systems onthe L/a = 48 volume and nπ + systems with n (cid:46) on the L/a = 32 volume. Given thesesystematic uncertainties in conjunction with the statistical uncertainties on η π + π + π + ( m K ) and η K K K ( m K ) , the results of this work do not provide significant evidence for differencesin 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 orneutral mesons as well as systems of one, two, and three nucleons. Calculations were per-formed in two lattice volumes with charge-dependent quark masses tuned such that strongisospin breaking effects are negligible and energy differences between charged and neutral sys-tems are primarily QED effects. While the ground-state energies of two- and three-nucleonsystems are determined with few-percent-level precision, QED effects leading to differencesbetween two-nucleon and three-nucleon FV energies are not resolved. Significantly higherprecision will be needed in future calculations of QED effects in multi-nucleon systems. Dif-ferences between charged and neutral FV energy shifts are resolved at the level of − σ for systems of 3-12 mesons, demonstrating the presence of QED effects on meson-mesoninteractions. Analysis of the FV energy levels for multi-meson systems using non-relativisticEFT has allowed the extraction of the π + π + and K K scattering lengths as well as the π + and K interaction parameters. Differences between a π + π + C and a K K are clearly resolved,demonstrating that additional QED effects on meson-meson interactions can be resolvedbeyond the Coulomb photon exchange explicitly included in the EFT. Differences betweenthe 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 workare well-described by NRQED L results that incorporate short-range two- and three-bodycontact interactions as well as perturbative Coulomb photon exchange. Although Coulombphoton exchange must be treated nonperturbatively in sufficiently large volumes, the expan-sion parameter describing the size of FV Coulomb effects is found to be α/v = αM L/ (4 π ) byexamining the convergence pattern of the NRQED L expansion. This includes a numericallysignificant factor of / (4 π ) compared to the parameter αM L discussed in Ref. [42]. For sys-tems with unphysically large α and quark masses such as those studied here, αM L/ (4 π ) (cid:28) is satisfied for volumes satisfying L (cid:28) fm for nucleons ( L (cid:28) fm for pions), andCoulomb corrections to the LO strong interaction FV energy shift appear perturbative for L (cid:46) fm for nucleons ( L (cid:46) fm for pions). For calculations with physical α and quarkmasses, FV Coulomb effects are reduced by a factor of 20 for nucleons (40 for pions) and areexpected to be perturbative for all practically accessible lattice volumes. The EFT analysisof this work also demonstrates that NRQED L results for FV energy shifts are unaffected bythe complications of photon zero-mode subtraction up to effects suppressed by O (( M L ) − ) 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 QEDfine structure constant, can therefore be interpreted using hadronic EFTs with perturbativeCoulomb effects with a similar procedure to the one undertaken here. Such calculations willgive insight into the quark mass dependence of QED effects on meson-meson interactions,6and, combined with higher-precision calculations of multi-nucleon FV energy levels, willpermit first principles predictions of QED effects on nucleon-nucleon interactions and QEDeffects in light nuclei. Acknowledgements
We thank the other members of the NPLQCD collaboration, in particular Z. Davoudi andM.J. Savage, for discussions during the initial stages of this work. Calculations in this projectwere performed using the Hyak High Performance Computing and Data Ecosystem at theUniversity of Washington (https://itconnect.uw.edu/research/hpc/), supported, in part, bythe U.S. National Science Foundation Major Research Instrumentation Award, Grant Num-ber 0922770, and on clusters at MIT with support from the NEC Corporation Fund. Thenumerical configuration generation (using the BQCD lattice QCD program [91]) was carriedout on the IBM BlueGene/Q and HP Tesseract using DIRAC 2 resources (EPCC, Edin-burgh, 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, Aus-tralia (supported by the Australian Commonwealth Government) and Phoenix (Universityof Adelaide). The Chroma software library [92] was used in the data analysis. SRB is sup-ported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physicsunder grant Contract Number DE-FG02-97ER-41014. WD, PES, and MLW are supportedin part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics un-der grant Contract Number DE-SC0011090. WD is also supported within the frameworkof the TMD Topical Collaboration of the U.S. Department of Energy, Office of Science,Office of Nuclear Physics, and by the SciDAC4 award DE-SC0018121. RH is supportedby STFC through grant ST/P000630/1. MI is supported by the Universitat de Barcelonathrough the scholarship APIF, by the Spanish Ministerio de Economia y Competitividad(MINECO) under Project No. MDM-2014-0369 of ICCUB and with additional EuropeanFEDER funds, under the contract FIS2017-87534-P. HP is supported by DFG Grant No. PE2792/2-1. PELR is supported in part by the STFC under contract ST/G00062X/1. GS issupported by DFG Grant No. SCHI 179/8-1. PES is supported in part by the National Sci-ence Foundation under CAREER Award 1841699. MLW was supported in part by an MITPappalardo Fellowship. RDY and JMZ are supported by the Australian Research CouncilGrants FT120100821, FT100100005, DP140103067 and DP190100297. This manuscript hasbeen authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.
Appendix A: Matching NRQED L and QED L Neglecting relativistic effects suppressed by M − (including nonlocal effects from zero-mode subtraction discussed in Refs. [2, 29, 31, 35, 37, 38]), the NRQED L Lagrangian isidentical for bosons and fermions. For a particle with charge Q = 1 , it is given by L NRQED L = ψ † (cid:18) iD − D i D i M (cid:19) ψ − πaM ( ψ † ψ ) + L ξγ = ψ † (cid:18) i∂ − ∂ i ∂ i M (cid:19) ψ − eA ( ψ † ψ ) − πaM ( ψ † ψ ) + L ξγ , (A1)7where 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 ingeneric R ξ gauge the photon Lagrangian is L ξγ = − F µν F µν + 12 ξ ( ∂ µ A µ ) . (A2)Results for the Landau gauge QCD+QED L calculations performed in the main text areobtained by setting ξ = 0 . In what follows, a finite spatial volume of extent L with PBCsis considered. Zero mode subtraction can be implemented in NRQED L by defining the FVphoton field as a Fourier transform of the zero-mode subtracted momentum-space field, A µ ( x ) = (cid:90) dp π L (cid:88) n ∈ Z \{ } e − ip x + πiL n · x ˜ A µ ( p , n ) , ˜ A µ ( p , n ) = (cid:90) d xe ip x − πiL n · x A µ ( x ) , (A3)where n ∈ Z \ { } excludes the photon zero mode. Zero-mode subtraction can be defined atthe path integral level as a constraint on the photon field [37], but for perturbative matchingwith QED L Eq. (A1)-(A3) define NRQED L .The NRQED L massive particle propagator G NRQED L is given by G NRQED L ( p , n ) = (cid:90) d xe ip x − πiL n · x (cid:10) ψ ( x ) ψ † (0) (cid:11) = ip − M (cid:0) π n L (cid:1) + i(cid:15) . (A4)The photon propagator is given by G γ L µν ( p , n ) = (cid:90) dx (cid:88) x e ip x − πiL p · x (cid:104) A µ ( x ) A ν (0) (cid:105) = i (cid:104) g µν − (1 − ξ ) p µ p ν p (cid:105) ( p ) − (cid:0) π n L (cid:1) + i(cid:15) . (A5)Introducing Fourier transformed fields, ˜ ψ ( p , n ) = (cid:90) d xe ip x − πiL n · x ψ ( x ) , ψ ( x ) = (cid:90) dp π L (cid:88) n ∈ Z e − ip x + πiL n · x ˜ ψ ( p , n ) . (A6)Matching between QED L and NRQED L is performed for four-point correlation functionsdescribing (off-shell) particles with energy M and three-momentum p ≡ πL r with r ∈ Z , M NRQED L = G NRQED L (0 , r ) − (cid:68) ˜ ψ (0 , r ) ˜ ψ (0 , r ) † (cid:69) . (A7)At tree-level this correlation function is given by M NRQED L LO = − πaM . (A8)The choice r (cid:54) = 0 , which does not affect the tree-level correlation function, is made in order toregulate IR divergences in the one-loop correlation function. Although one-photon-exchangecontributions lead to an IR divergence (regulated for instance by considering non-zero mo-mentum transfer) in the NRQED analog of Eq. (A8), one-photon-exchange contributions tothe NRQED L amplitude in Eq. (A7) vanish because of zero-mode subtraction.8Matching is performed by expanding M NRQED L and its QED L analog M QED L perturba-tively in the small parameters α , / ( M L ) , and a/L and defining higher-order terms in theNRQED L Lagrangian so that QED L and NRQED L agree order by order in this expansion.Matching will be performed to leading order in / ( 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 equivalentto the non-relativistic limit with up to two-loop contact interactions and Coulomb photonexchange. This matching is also sufficient for N LO in the power counting of Eq. (34), whichdecreases the relative importance of Coulomb photon exchange and only requires one-loopCoulomb 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 L QED L = − ( D µ ϕ ) † D µ ϕ − M ϕ † ϕ − πaM ( ϕ † ϕ ) + L ξγ = − ( ∂ µ ϕ ) † ∂ µ ϕ − M ϕ † ϕ + ieA µ ( ϕ † ∂ µ ϕ − ( ∂ µ ϕ ) † ϕ ) − e A µ A µ ϕ † ϕ − πaM ( ϕ † ϕ ) + L ξγ , (A9)where L ξγ is defined in Eq. (A2). The fields appearing in the QED L and NRQED L La-grangians are related by ψ ( x ) = √ M e iMt ϕ ( x ) . (A10)The scalar QED L particle propagator is G QED L ( p , n ) = (cid:90) dx (cid:88) x e ip µ x µ (cid:10) ϕ ( x ) ϕ † (0) (cid:11) = i ( p ) − (cid:0) π n L (cid:1) − M + i(cid:15) . (A11)The QED L photon propagator is identical to the NRQED L propagator Eq. (A5). Introduc-ing 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 M QED L LO = 14 M G QED L ( M, r ) − (cid:10) ˜ ϕ ( M, r ) ˜ ϕ ( M, r ) † (cid:11) = − πaM , (A12)in agreement with Eq. (A8). This demonstrates that the four-fermion contact interactionsin Eq. (A1) and Eq. (A9) are normalized consistently at tree-level.Corrections to the LO amplitude arise from one-loop diagrams shown in Figs. 6-8. Theone-loop diagrams shown in Fig. 6 involve similar sums/integrals over loop momenta anddiffer only in the number of photon propagators n γ and contact interactions present. Thecontribution to the amplitude from each diagram is denoted M NRQED L NLO ( n γ ) where n γ ∈ { , , } labels the numbers of photon propagators present, M NRQED L NLO ( n γ ) = i π a M L (cid:18) αMa (cid:19) n γ (cid:88) n ( (cid:48) ) (cid:90) dk π (cid:32) k − π ( n − r ) ML + i(cid:15) (cid:33) × (cid:32) k + π ( n + r ) ML − i(cid:15) (cid:33) − (1 − ξ ) ( k ) ( k ) − π n L + i(cid:15) ( k ) − π n L + i(cid:15) n γ = − π a M L (cid:18) αMa (cid:19) n γ (cid:88) n ( (cid:48) ) (cid:20) Res (cid:18) − π ( n − r ) M L + i(cid:15) (cid:19) + Res (cid:18) − π | n | L + i(cid:15) (cid:19)(cid:21) , (A13)where Res ( x ) indicates the residue of the integrand in the first line at the pole x and (cid:80) n ( (cid:48) ) corresponds to (cid:80) n ∈ Z for n γ = 0 and to (cid:80) n ∈ Z \{ } for n γ ≥ . After taking the (cid:15) → limit,the residue at the particle pole π ( n + r ) / ( M L ) is given byRes (cid:18) − π ( n − r ) M L (cid:19) = − M L π ( n + p ) π ( n − r ) M L − π n L − (1 − ξ ) π ( n − r ) M L (cid:16) π ( n − r ) M L − π n L (cid:17) n γ = − M L π (cid:18) − L π (cid:19) n γ n + r ) n n γ (cid:2) O (cid:0) ( M L ) − (cid:1)(cid:3) , (A14)where the last line includes an expansion in powers of ( M L ) − . This expansion is legitimateprovided that the sum over n converges, which holds for n γ ≥ . For n γ = 0 there is a linearUV divergence that can be removed by adding a UV counterterm, (cid:88) n ∈ Z n + r → lim Λ →∞ | n | < Λ (cid:88) n ∈ Z n + r − π Λ , (A15)and after making the replacement of Eq. (A15) the / ( M L ) expansion can be performed inEq. (A14). The photon pole residue at − π | n | /L involves energy denominators of order /L and is suppressed by O (( M L ) − ) compared to the particle-pole contributions in Eq. (A14)where particle propagator energy denominators are of order M/L . For all diagrams inFigs. 7-8, particle poles only appear in either the upper or lower half of the complex planeand only photon poles contribute. After performing the energy integrals these diagramscan be straightforwardly verified to be suppressed by O (( M L ) − ) or O (( M L ) − ) comparedto Eq. (A16). The full NLO amplitude in NRQED L is therefore a sum of the amplitudes M NRQED L NLO ( n γ ) associated with the diagrams in Fig. 6 obtained by substituting Eq. (A14)into Eq. (A13), M NRQED L NLO ( n γ ) = 4 πaM (cid:18) − αM L π (cid:19) n γ (cid:16) aπL (cid:17) − n γ (cid:88) n ( (cid:48) ) n + r ) n n γ (cid:2) O (cid:0) ( M L ) − (cid:1)(cid:3) . (A16)Contributions with different n γ differ parametrically and can be matched independentlybetween NRQED L and QED L .0The QED L amplitudes associated with the NLO diagrams in Fig. 6 are given by M QED L NLO ( n γ ) = − i π a L (cid:18) − α aM (cid:19) n γ (cid:90) dk π (cid:88) n ( (cid:48) ) k + M ) − π ( r + n ) L − M + i(cid:15) × k − M ) − π ( r − n ) L − M + i(cid:15) (cid:32) N ξγ ( k , n )( k ) − π n L + i(cid:15) (cid:33) n γ , (A17)where N ξγ ( k , n ) = ( k ) − M + (cid:18) π L (cid:19) (4 r − n ) − (1 − ξ ) − ( k ) + 4 M ( k ) + (cid:16) π L (cid:17) (4 n · r − n )( k ) − (cid:0) π L (cid:1) n + i(cid:15) . (A18)The k integral can be performed with contour integration. Closing the contour inthe upper half plane, the result includes contributions from a particle pole k = M − (cid:112) M + 4 π ( n + r ) /L = − π ( n + r ) / ( M L ) [1 + O (( M L ) − )] , a photon pole at k =2 π | n | /L , and an antiparticle pole at k = − M [1 + O (( M L ) − )] , M QED L NLO ( n γ ) = 64 π a L (cid:18) − α aM (cid:19) n γ (cid:88) n ( (cid:48) ) (cid:34) Res (cid:32) M − (cid:114) M + 4 π ( n + r ) L − i(cid:15) (cid:33) + Res (cid:18) − π | n | L + i(cid:15) (cid:19) + Res (cid:32) − M − (cid:114) M + 4 π L ( n − r ) + i(cid:15) (cid:33)(cid:35) . (A19)Taking the (cid:15) → limit and expanding to leading order in ( M L ) − , the residue at the particlepole is given byRes (cid:32) M − (cid:114) M + 4 π ( n + r ) L − i(cid:15) (cid:33) = L π M (cid:18) M L π (cid:19) n γ n + r ) n n γ . (A20)As in the NRQED L case, the residue at the photon pole − π | n | /L is suppressed comparedto the residue at the particle pole by O ( M L ) − .The residue at the antiparticle pole at − M (1 + O ( M L ) − ) does not appear in thecorresponding NRQED L expression Eq. (A13) and is therefore associated with contributionsthat do not appear in loop diagrams in NRQED L . This residue is given byRes (cid:32) − M − (cid:114) M + 4 π L + i(cid:15) (cid:33) = − M (cid:18) − π M L (cid:19) n γ × (cid:2) r ( − ξ ) + n ( − ξ ) − r · n ( − ξ ) (cid:3) n γ . (A21)The term involving two contact interactions involves the UV divergent sum (cid:80) n ( (cid:48) ) , whichfor n γ = 0 must be consistent with the infinite-volume result (cid:82) d D k after subtractingUV counterterms, (cid:88) n ∈ Z → . (A22)1This leads to a vanishing contribution from the n γ = 0 diagram appearing in the absenceof QED L , which is consistent with the expectation that antiparticle-pole contributions fromoff-shell intermediate states that are absent non-relativistically do not lead to power-law FVeffects in local field theories. The terms with n γ ≥ differ by zero mode subtraction. Sincethe n = 0 contribution to (cid:80) n is unity, it follows from Eq. (A22) that (cid:88) n ∈ Z \{ } → − . (A23)Sums with n k with k > similarly require UV counterterms and vanish after includingthem. Zero-mode contributions to these sums vanish, and so (cid:88) n ∈ Z \{ } n k → , (A24)for k > . After subtracting UV counterterms the antiparticle pole contribution becomes (cid:88) n ∈ Z \{ } Res (cid:32) − M − (cid:114) M + 4 π L + i(cid:15) (cid:33) → M (cid:18) − π r ( − ξ ) M L (cid:19) n γ , (A25)for n γ ≥ . This leads to a finite contribution to M QED L NLO that does not arise in the cor-responding loop diagrams contributing to M NRQED L NLO . Local counterterms involving powersof L − must be added to the NRQED L Lagrangian in order to reproduce this relativisticeffect arising from zero-mode subtraction. Inserting Eq. (A25) into Eq. (A19) shows thatthe dominant contribution to M QED L NLO arising from the n γ = 1 diagram antiparticle pole is − πaM (cid:18) απ ( − ξ ) r M L (cid:19) = − πaM (cid:18) π ( − ξ )16 (cid:19) α ( M L ) (cid:18) p M (cid:19) . (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 (cid:28) M in the ( M L ) − expansion above), this contribution is proportional to the LOcontact interaction times the velocity squared times a relativistic QED suppression factorof α/ ( M L ) − . For boosted systems, the NRQED L contact interaction in Eq. (A1) musttherefore be supplemented with a nonlocal counterterm suppressed by the same factor of α/ ( M L ) − . This is analogous to the self-energy of a single scalar field, which is shown inRef. [37] to require a nonlocal counterterm equal to the scalar mass times the velocity squaredtimes a relativistic QED suppression factor of α/ ( M L ) − . The LQCD+QED L calculationsdiscussed in the main text are performed in the center-of-mass rest frame, and these nonlocalcounterterms vanish. Non-vanishing nonlocal counterterms appear for fermion masses inNRQED L even in the two-particle rest frame, but since effects proportional to ( M L ) − are neglected throughout this work it is consistent to neglect all nonlocal countertermssuppressed by α/ ( M L ) − .The NLO QED L amplitude is therefore given by inserting the particle-pole residue inEq (A20) into Eq. (A19), M QED L NLO ( n γ ) = 4 πaM (cid:18) − αM L π (cid:19) n γ (cid:16) aπL (cid:17) − n γ (cid:88) n (cid:48) n + r ) n n γ (cid:2) O (cid:0) ( M L ) − (cid:1)(cid:3) , (A27)2which is identical to Eq. (A16). The ( M L ) − suppression of diagrams in Figs. 7 and 8 arisingfrom the absence of particle pole contributions is identical in QED L and NRQED L . Addi-tional diagrams appearing in QED L but not NRQED L associated with the two-particle-two-photon vertex or particle-antiparticle pair creation can be similarly verified to be suppressedby powers of ( M L ) − . 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 beexpressed in terms of the amplitude M NRQED L NNLO ( n , n , n )= 64 π a M L (cid:18) αMa (cid:19) n + n + n (cid:88) n , m ( (cid:48) ) (cid:90) dk π (cid:16) k + π ( r − n ) ML − i(cid:15) (cid:17) (cid:16) k − π ( r + n ) ML + i(cid:15) (cid:17) × − (1 − ξ ) ( k ) ( k ) − π n L ( k ) − π n L + i(cid:15) n (cid:90) dq π (cid:16) q + π ( r − m ) ML − i(cid:15) (cid:17) (cid:16) q − π ( r + m ) ML + i(cid:15) (cid:17) × − (1 − ξ ) ( k − q ) ( k − q ) − π n − m )2 L ( k − q ) − π ( n − m ) L + i(cid:15) n − (1 − ξ ) ( q ) ( q ) − π m L ( q ) − π m L + i(cid:15) n , (A28)where n i ∈ { , } for i ∈ { , , } labels whether each interaction is a four-particle contactinteraction or photon exchange and (cid:80) n , m ( (cid:48) ) excludes n = 0 if n = 1 or n = 1 and excludes m = 0 if n = 1 or n = 1 . Both energy integrands include poles where a particle is on-shellas well as poles where a photon is on-shell, and, as above, photon-pole contributions aresuppressed by O (( M L ) − ) . Evaluating the q and k integrals then gives M NRQED L NNLO ( n , n , n ) = − πaM (cid:16) aπL (cid:17) − n − n − n (cid:18) − αM L π (cid:19) n + n + n × (cid:88) n , m ( (cid:48) ) n + r )( m + r ) n n ( n − m ) n m n (cid:2) O (cid:0) ( M L ) − (cid:1)(cid:3) , (A29)where UV counterterms should be included as in Eq. (A15). Similarly to the NLO case, theNNLO diagrams in Figs. 7 and 8 only include photon pole contributions and are suppressedby O (( M L ) − ) or O (( M L ) − ) compared to Eq. (A29).3The corresponding QED L NNLO amplitudes are given by M QED L NNLO ( n , n , n ) = 1024 π a ML (cid:16) − α aM (cid:17) n + n + n (cid:88) n , m ( (cid:48) ) (cid:90) dk π (cid:32) N ξγ ( k , n )( k ) − π n L + i(cid:15) (cid:33) n × (cid:16) ( k − M ) − π ( n − r ) L − M + i(cid:15) (cid:17) (cid:16) ( k + M ) − π ( n + r ) L − M + i(cid:15) (cid:17) × (cid:90) dq π (cid:32) N ξγ ( k − q , n − m )( k − q ) − π ( n − m ) L + i(cid:15) (cid:33) n (cid:32) N ξγ ( q , m )( q ) − π m L + i(cid:15) (cid:33) n × (cid:16) ( q − M ) − π ( m − r ) L − M − i(cid:15) (cid:17) (cid:16) ( q + M ) − π ( m + r ) L − M + i(cid:15) (cid:17) . (A30)Contributions from photon and antiparticle poles are again suppressed by powers of ( M L ) − ,and evaluating the q and k energy integrals gives M QED L NNLO ( n , n , n ) = − πaM (cid:16) aπL (cid:17) − n − n − n (cid:18) − αM L π (cid:19) n + n + n × (cid:88) n , m ( (cid:48) ) n + p )( m + p ) n n ( n − m ) n m n (cid:2) O (cid:0) ( M L ) − (cid:1)(cid:3) . (A31)All other diagrams are again suppressed by powers of ( M L ) − , and agreement betweenEq. (A29) and Eq. (A31) shows that nonlocal counterterms are suppressed by powers of ( M L ) − 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 ) − 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 ) − suppres-sion of antiparticle pole contributions leading to nonlocal two-body counterterms also arisesfrom the denominator structure of the QED L propagators and is expected to be genericfor bosons and fermions. The numerator structure of the scalar QED L propagators aboveis relevant for ( M L ) − suppressed relativistic effects, and in particular the scalar QED L antiparticle pole in Eq. (A21) includes vanishing numerator factors for a system at restthat lead to p /M velocity suppression of nonlocal counterterms in scalar NRQED L . Thiscancellation might be absent for fermions, and nonlocal two-body counterterms might berelevant for QED L calculations of charged fermions in the center-of-mass rest frame. Thiswould parallel the situation for one-body nonlocal counterterms, which arise at O ( α/ ( M L ) ) for fermions with any center-of-mass velocity and for scalars with non-zero center-of-massvelocity. Appendix B: Fitting procedures1. Energy level determination
In general the spectral representation in Eq. (8) cannot be inverted to determine the fullenergy spectrum from finite samples of correlation functions over a finite range of source/sink4separations. Any fitting procedure to extract energies from correlation functions involvesmaking several choices, in particular the range of t to include in the fit, the number of excitedstates to include in a truncation of Eq. (8) to use as a fitting model, and how to estimatethe covariance matrix from a finite statistical ensemble. In order to assess the systematicuncertainties associated with these fitting choices and provide a reproducible procedure forextracting energy levels from correlation function results, we use an approach detailed herefor making fitting choices based on well-defined statistical criteria and random sampling overthe 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 smallercontributions to χ when fitting energy levels from correlation functions. For the nucleon, thescaling StN ( G ) = G ( t ) / (cid:112) Var ( G ( t )) ∼ e − ( M N − m π ) t predicted by Parisi [93] and Lepage [94]applies in the limit of a large statistical ensemble size N → ∞ and shows that fit resultsare exponentially insensitive to the choice of t max . Similar results apply for multi-nucleonsystems with baryon number A where StN ( G ) = G ( t ) / (cid:112) Var ( G ( t )) ∼ e − A ( M N − m π ) t for large t . For fixed N (and, since volume-averaging increases StN, for fixed L ) it is important tochoose a fixed t max before the StN has degraded to the point where correlation functionestimates are unreliable. The StN ratio decreases with t at small and intermediate t beforesaturating at an O (1) value in the noise region, and to avoid the noise region t max shouldnot exceed the smallest t where StN ( G ) = G ( t ) / (cid:112) Var G ( t ) reaches a specified O (1) cutoff.For positive-definite meson correlation functions this issue does not arise, and the maximum t that can be reliably included in the fit is only limited by the accuracy by which finite-temperature effects are modeled. In this work, finite-temperature excited-state effects areneglected, and finite-temperature effects on correlation functions with t ∼ T / , where T isthe length of the Euclidean time direction, are not reliably modeled. In our fitting procedure, t max is therefore chosen to be the minimum t satisfying either: the correlation function noise-to-signal ratio at ( t + a ) is smaller than a specified tolerance (final results use tol noise = 1 . ),the correlation function sample mean at t + a is negative, or t + a is larger than a finite-temperature cutoff (final results use tol temp = 3 T / ). The values of tol noise and tol temp arefree parameters in our fitting procedure that must be varied to assess the sensitivity of fitresults to these choices; for concreteness the parameter choices are presented here that leadto the final results quoted in the main text. Results are found to be relatively insensitive tothe parameter choices controlling t max , and for example varying the noise tolerance in therange tol noise ∈ [0 . , leads to results with consistent central values at the σ level andfew-percent variation of the corresponding uncertainties. The complex phases of baryon correlation functions are circular random variables, and there is thereforea noise region at large t where ln N is not much larger than the variance of the phase distribution and thesample mean is a systematically unreliable estimator of the average correlation function [95]. If the average correlation function is not expected to be positive definite, then this condition should notbe enforced. In principle, the source and sink interpolating operators used in this work differ from oneanother, and it is possible for the sign of an average correlation function to fluctuate at small t whereexcited-state contributions with opposite sign to the ground-state contribution can be significant. Inpractice, small t fluctuations of the sign of the sample mean correlation function that could be attributedto excited-state effects are not observed in this work, and a negative sample mean correlation function istaken as an indicator of large statistical noise. t min , the minimum t included inthe 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, theuncertainties 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 thechoice of t min than to the choice of t max . Rather than choose a single t min , it is preferableto sample from many possible choices of t min and quantify the sensitivity to this choice asthe systematic error. The minimum permissible t min is fixed by the temporal nonlocalityin the lattice action, and for the improved action used in this work the transfer matrixinvolves fields on two adjacent timeslices [96] and t min ≥ 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 foundto significantly affect final results when varied over the range (cid:46) t plateau (cid:46) (final resultsuse t plateau = 4 ). For each type of interpolating operator included in a combined fit, t min issampled randomly within this range until either all possible values of t min have been chosenor a maximum of N fits = 200 fits have been performed.With t min and t max specified, the covariance matrix C ijtt (cid:48) must be estimated for t min ≤ t, t (cid:48) ≤ t max and interpolating operators i, j ∈ { , . . . N op } . Fits involving a large number N pts of time separations and interpolating operator choices may not satisfy the condition N (cid:29) N needed to ensure that the N terms contributing to the N pts × N pts samplecovariance can accurately estimate the true underlying covariance matrix, where N pts isthe total number of source/sink separations from all interpolating operators included in thefit. Shrinkage techniques have been developed to provide more accurate estimates of theunderlying covariance matrix than the sample covariance matrix when N (cid:29) N is notsatisfied [59, 60]. Shrinkage estimators S ijtt (cid:48) ( λ ) of the covariance matrix are constructed asmixtures of a well-conditioned target matrix T ijtt (cid:48) and the covariance matrix C ijtt (cid:48) estimatedusing standard bootstrap techniques from N boot = 200 samples of N correlation functions G bi,a ( t ) with a ∈ { , . . . , N } and b ∈ { , . . . , N boot } drawn from the original correlationfunction ensemble with replacement, S ijtt (cid:48) ( λ ) = C ijtt (cid:48) (1 − λ ) + T ijtt (cid:48) λ, (B1)where ≤ λ ≤ is a shrinkage parameter. A common choice of well-conditioned targetmatrix for many problems in statistics is the identity matrix; however, this does not accu-rately describe the underlying covariance matrix for correlation functions whose diagonalentries decrease exponentially with t . Following applications of shrinkage to lattice QCD inRef. [61], we take T ijtt (cid:48) = diag ( C ijtt (cid:48) ) . In this case, shrinkage corresponds to an interpolationbetween a fully correlated fit with λ = 0 and an uncorrelated fit with λ = 1 . Shrinkage givesan unbiased estimator of the underlying covariance matrix in the infinite-statistics limitprovided λ vanishes sufficiently quickly in this limit. It can be shown [60] that for finite N an optimal λ ∗ (cid:54) = 0 satisfying this restriction can be chosen in order to minimize the averagemean-squared difference between S ijtt (cid:48) ( λ ) and the underlying covariance matrix, and that asample estimate for λ ∗ is given by λ ∗ = Max , Min , (cid:80) Na =1 (cid:80) t,t (cid:48) ,i,j (cid:20) (cid:101) G ia ( t ) (cid:101) G ja ( t (cid:48) ) − (cid:101) S ij ( t, t (cid:48) ) (cid:21) N (cid:80) t,t (cid:48) ,i,j (cid:20)(cid:101) S ij ( t, t (cid:48) ) − δ ij δ t,t (cid:48) (cid:21) , (B2)6where (cid:101) G ia ( t ) = G ia ( t ) − G i ( t ) (cid:113) S ii ( t, t ) , (cid:101) S ij ( t, t (cid:48) ) = S ij ( t, t (cid:48) ) (cid:113) S ii ( t, t ) S jj ( t (cid:48) , t (cid:48) ) , (B3)are defined in terms of the sample mean correlation function and sample covariance as inRef. [61] G i ( t ) = 1 N N (cid:88) a =1 G ia ( t ) ,S ij ( t, t (cid:48) ) = 1 N − N (cid:88) a =1 (cid:104) G ia ( t ) G ja ( t (cid:48) ) − G i ( t ) G j ( t (cid:48) ) (cid:105) , (B4)such that shrinkage of (cid:101) S ij ( t, t (cid:48) ) with T ijtt (cid:48) = diag ( (cid:101) S ijtt (cid:48) ) corresponds to shrinkage of (cid:101) S ij ( t, t (cid:48) ) with the identity matrix as a target, and the results of Ref. [60] assuming an identity matrixtarget can be applied. The covariance matrix estimate with optimal shrinkage is then givenby S ijtt (cid:48) ( λ ∗ ) . Fits to truncations of Eq. (8) including e excited states can then be performedby minimizing the corresponding χ function defined by χ = t max (cid:88) t,t (cid:48) = t min N op (cid:88) i,j =1 (cid:16) G Bi ( t ) − f ( t, E , Z ) (cid:17) (cid:2) S ( λ ∗ ) − (cid:3) ijtt (cid:48) (cid:16) G Bj ( t (cid:48) ) − f ( t (cid:48) , E , Z ) (cid:17) , (B5)where G Bi = 2 G i − N boot (cid:80) N boot b =1 G bi includes a /N bias correction estimated using bootstraptechniques and E and Z denote the energies and overlap factors appearing in Eq. (8), includ-ing excited-state and thermal effects. Since the overlap factors enter f N linearly, the valuesof Z minimizing χ for fixed E can be determined by solving a system of linear equationsanalogous to variable projection techniques [97, 98]. χ -minimization can therefore be effi-ciently 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 optimizationfor E . In order to ensure positivity of the spectrum and remove fitting degeneracies, theparameters used for nonlinear optimization are ln E and ln( E k − E k − ) for ≤ k ≤ e .For each randomly sampled choice of t min , the next step in the fitting procedure is todetermine the number of excited states to be included in the sum of exponentials used asa fit function. This is done by first performing a fit including zero excited states and thenadding successively more excited states until the addition of the next excited state doesnot improve the goodness of fit according to an information criterion. This work employsthe Akaike Information Criterion [99] (AIC) with a cutoff chosen to penalize overfittingin which a fit with e excited states is only preferred over a fit with e − excited statesif AIC ( e ) − AIC ( e − < −A N dof ( e ) where N dof ( e ) = N pts − N params ( e ) is the numberof 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 ) = 2 N params ( e ) + χ ( e ) + k with χ ( e ) the (unreduced) χ of the fitdefined in Eq. (B5) and k is an irrelevant e -independent constant. This choice corresponds7to a preference for an e state fit only if it improves the χ /N dof by an O (1) value A comparedto the ( e − state fit. For baryon correlation functions a value of A = 0 . is used, whilefor multi-meson correlation functions a value of A = 0 . is used. Bootstrap resampling techniques are then used to estimate the uncertainty on the ground-state energy extracted from the fit with the preferred number of excited-states in each fitregion [100, 101]. The same χ -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 confidence interval for the ground-state energy is then obtained from the quantilesof the distribution of differences between the b ’th bootstrap sample result E b,f and the fitresult E f obtained for fit range f , δE f = 12 (cid:104) Q / (cid:16) E b,f − E f (cid:17) − Q / (cid:16) E b,f − E f (cid:17)(cid:105) , (B6)where Q p ( x f ) is the p -th quantile of the set of fit results with elements x f . Using thisdefinition δE f minimizes the impact of outlier bootstrap samples compared to a definitionbased on the standard deviation of E b,f − E f [101]. An analogous procedure is used toestimate uncertainties for excited-state energies and overlap factors.Several additional checks are used to ensure the robustness of χ -minimization results:two different optimization algorithms, Nelder-Mead (NM) and conjugate gradient (CG),are used and are verified to give energies that differ by less than a specified tolerance (final results use tol sol = 10 − ), median results for ground- and excited-state energies fromthe bootstrap samples are verified to agree with fit results from the average correlationfunctions for each energy level within a specified tolerance (final results use tol med = 2 σ ),uncorrelated fit results obtained by repeating the χ -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 χ /N dof is verified to be less than a specified tolerance(final results use tol χ = 2 ). This defines a reproducible and automatable procedure forfitting correlation functions, including sampling of possible fit ranges and excited-state modelselection, in which fit results are functions of only the tolerances described above and thegiven correlation functions. A graphical illustration of the fitting procedure is shown inFig. 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 energyspectrum, and the final estimate of the ground-state energy E and its uncertainty δE are Optimal shrinkage values of λ ∗ (cid:38) . appear for n -meson correlation functions with n (cid:38) and χ valuesare correspondingly lower than would be expected for fully correlated χ minimization, leading to smallerabsolute changes in AIC for multi-meson correlation functions than for multi-baryon correlation function.Increasing A from . over the range A ∈ [0 . , leads to consistent results with larger uncertaintiesbecause precise and accurate two-state fits with small t min are rejected in favor of one-state fits morefrequently. Newton’s method is used in place of Nelder-Mead if N states = 1 , since Nelder-Mead does not work for asingle fit parameter. Fits resulting in ground-state energies less than tol sol , which appeared only for the K system on the L/a = 32 lattice volume, were also rejected. Correlation functions G ( t ) t min ∈ [2 a, t max − t plateau ] t max = min { t | (cid:2) / StN ( G ( t + a )) > tol noise (cid:3) ∨ (cid:2) G ( t + a ) < (cid:3) ∨ [ t + a > tol temp ] } Excited-state model selection: f ( t, E , Z ) = e (cid:80) n =0 Z n e − E n t , e = 0 χ minimization with Nelder-Mead+VarPro using S ( λ ∗ ) → { E f , Z f } ∆ AIC < −A N dof | E f (cid:48) − E f | > tol sol yes e ← e + 1 e ← e − Reject fit χ N dof > tol χ Accept fitno χ minimization withCG+VarPro using S ( λ ∗ ) → { E f (cid:48) , Z f (cid:48) } Confidence intervals: χ minimization withNM+VarPro using S ( λ ∗ ) over bootstrap ensemble → { E b,f , Z b,f } Reject fitnoyes yes no Reject fitReject fit δ E f = Q / ( E b,f − E f ) − Q / ( E b,f − E f ) | E f (cid:48)(cid:48) − E f | > tol corr χ minimization withNM+VarPro using S (1) → { E f (cid:48)(cid:48) , Z f (cid:48)(cid:48) } yesno Q / (cid:0) E b,f − E f (cid:1) > tol med yesnoCovariance 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. Inputparameters to the fitting procedure are shown in blue. As described in the text, the steps illustratedhere are repeated N fits times with different random choices of t min , and final results are obtainedfrom weighted averages of fit results for the t min choices leading to the “Accept fit” rectangle. N success successful fit results E f , E = N success (cid:88) f =1 w f E f ,δ stat E = N success (cid:88) f =1 w f ( δE f ) ,δ sys E = N success (cid:88) f =1 w f (cid:16) E f − E (cid:17) ,δE = (cid:113) δ stat E + δ sys E , (B7)where f labels the choice of fit range specified by t min for each interpolating operator . Eachfit 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 largestatistics; in practice it is advantageous to choose weights that penalize poor fits with larger χ /N dof and unconstraining fits with larger uncertainties δE f . Following Ref. [61], we usethe weights (cid:101) w f = p f (cid:16) δE f (cid:17) − (cid:80) N success f (cid:48) =1 p f (cid:48) (cid:16) δE f (cid:48) (cid:17) − , (B8)where p f = Γ( N dof / , χ f / / Γ( N dof / is the p -value assuming χ -distributed goodness-of-fit parameters with χ f obtained by inserting E f into Eq. (B5) . Variation to the particularchoices of specified tolerances have been studied, and the subsequent variation in the ensem-ble of successful fits is found to have little impact on the results of this weighted averaging.The results E and δE obtained with this procedure are shown as the central values anduncertainties for single-particle energy results E π + ( L ) , E K ( L ) , E n ( L ) , and E p ( L ) in Ta-bles II and III. Effective mass plots showing the smallest t min fit with weight over / of themaximum weight fit as well as E f and (cid:101) w f / Max ( (cid:101) w f ) are shown in Appendix B 4.
2. Multi-meson correlation functions
To determine results for multi-meson ground-state energies with thermal effects takeninto account, fits are performed iteratively starting with fits for n = 1 mesons and then The total error δE describes the combined statistical uncertainty on E plus systematic uncertaintyarising from the choice of fit range and fit model. The partitioning of this error into δ stat E and δ sys E only partially separates statistical and systematic uncertainties because δ stat E includes statistical errorsplus systematic uncertainties related to fluctuations among the δE f . For large λ ∗ , the χ function being minimized approaches an uncorrelated χ and the values of χ will notbe distributed as χ -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 comparativelyless accurate descriptions of the results being fit and their correlations as estimated by S ijtt (cid:48) ( λ ∗ ) , but theabsolute sizes of the p f should not be interpreted as p -values for each fit. n . The excited-state fit form in Eq. (8) includes thermaleffects describing k forwards-propagating mesons and n − k backwards-propagating mesonsfor k < n/ that are included by using the central values E 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/ , and for simplicity are not incorporated into the description of thermaleffects. The overlap factors for these thermal states are determined using linear algebratechniques [97, 98] during each step of nonlinear optimization for the N -particle energyspectrum analogously to the procedure described above for other overlap factors. This fitfunction is found to provide acceptable fits to multi-meson correlation functions without theneed for additional free parameters describing excited-state thermal effects.Before beginning the fitting procedure described above, correlation function results fromall quark propagator sources on a given configuration are averaged and meson correlationfunctions are further blocked along the Markov chain to form N block = 200 approximatelyindependent samples from the N cfg configurations for each volume shown in Table I. Furtheraveraging is found to give statistically consistent results, suggesting that autocorrelations canbe neglected after this blocking. To determine correlated differences of ground-state energiesfor different hadron type ( π + , K ) and hadron number, fits are performed independentlyto determine E f for each hadron type and number for each fit range sampled. A fit rangeis considered to give a successful fit only if the checks on fit robustness described aboveare passed for each hadron type and number involved in the correlated difference. For eachsuccessful fit range, bootstrap resampling is used to determine the uncertainties on correlateddifferences of the resulting E f . 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 are then formed, and confidence intervalsare computed by applying Eq. (B6) to the these correlated differences. Finally, weightedaverages of the resulting correlated differences and their associated uncertainties are takenusing Eq. (B7)-(B8). The results of this procedure are used to determine the FV energyshifts and differences between FV energy shifts for charged and uncharged hadrons shownin Tables IV-III.
3. Multi-nucleon correlation functions
Differences between multi-nucleon ground-state energies and the corresponding sums oftheir constituent nucleon masses are computed using correlated differences of bootstrapresults E f,b for multi-nucleon and single-nucleon correlation functions analogously to themulti-meson case described above. Correlated differences between one-nucleon and multi-nucleon 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 sys-tems. Otherwise, fits for multi-nucleon energy shifts are performed identically to fits formulti-meson energy shifts not including thermal effects.1
4. 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 systems with L/a ∈{ , } . Single-nucleon fit results for p and n are shown in Fig. 26. Two-nucleon fit resultsare shown for pp and nn in Fig. 27 and for np systems in Fig. 28. Three-nucleon fit resultsare shown in Fig. 29. [1] S. Borsanyi et al. (Budapest-Marseille-Wuppertal), Phys. Rev. Lett. , 252001 (2013),1306.2287.[2] S. Borsanyi et al., Science , 1452 (2015), 1406.4088.[3] R. Horsley et al., J. Phys. G43 , 10LT02 (2016), 1508.06401.[4] R. Horsley et al., JHEP , 093 (2016), 1509.00799.[5] S. Basak et al. (MILC), Phys. Rev. D99 , 034503 (2019), 1807.05556.[6] R. Horsley et al. (CSSM, QCDSF, UKQCD), Journal of Physics G: Nuclear and ParticlePhysics , 115004 (2019), 1904.02304.[7] Z. R. Kordov, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, H. Stüben,R. D. Young, and J. M. Zanotti (CSSM/QCDSF/UKQCD), Phys. Rev. D101 , 034517 (2020),1911.02186.[8] V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanfilippo, S. Simula, and N. Tantalo, Phys.Rev.
D95 , 034504 (2017), 1611.08497.[9] D. Giusti, V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanfilippo, S. Simula, N. Tantalo,and C. Tarantino, Phys. Rev. Lett. , 072001 (2018), 1711.06537.[10] M. Di Carlo, D. Giusti, V. Lubicz, G. Martinelli, C. T. Sachrajda, F. Sanfilippo, S. Simula,and N. Tantalo, Phys. Rev.
D100 , 034514 (2019), 1904.08731.[11] T. Blum, S. Chowdhury, M. Hayakawa, and T. Izubuchi, Phys. Rev. Lett. , 012001 (2015),1407.2923.[12] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Phys. Rev.
D96 , 034515 (2017), 1705.01067.[13] S. Borsanyi et al. (Budapest-Marseille-Wuppertal), Phys. Rev. Lett. , 022002 (2018),1711.04980.[14] P. Boyle, V. Gülpers, J. Harrison, A. Jüttner, C. Lehner, A. Portelli, and C. T. Sachrajda,JHEP , 153 (2017), 1706.05293.[15] D. Giusti, V. Lubicz, G. Martinelli, F. Sanfilippo, and S. Simula, JHEP , 157 (2017),1707.03019.[16] T. Blum, P. A. Boyle, V. Gülpers, T. Izubuchi, L. Jin, C. Jung, A. Jüttner, C. Lehner,A. Portelli, and J. T. Tsang (RBC, UKQCD), Phys. Rev. Lett. , 022003 (2018),1801.07224.[17] J. Bijnens, J. Harrison, N. Hermansson-Truedsson, T. Janowski, A. Jüttner, and A. Portelli,Phys. Rev. D100 , 014508 (2019), 1903.10591.[18] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner (2019),1911.08123.[19] D. Giusti, V. Lubicz, G. Martinelli, F. Sanfilippo, and S. Simula, Phys. Rev.
D99 , 114502(2019), 1901.10462. π + L/a = 32 . . . . . . . . . t / a . . . . . . E a SP f . . . . . E f a π + L/a = 32 . . . . . . . . . t / a . . . . . . E a SP f . . . . . . . . . E f a π + L/a = 32 . . . . . . . . . t / a . . . . E a SP f . . . . . . . E f a π + L/a = 32 . . . . . . . . . t / a . . . . . . . E a SP . . . . . . . f . . . . . . . . E f a FIG. 17: Fit results for systems of n ∈ { , . . . , } π + mesons for the L/a = 32 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. π + L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP . . . . . . . f . . . . . E f a π + L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP f . . . . . . . . E f a π + L/a = 32 . . . . . . . . . t / a . . . . . E a SP f . . . . . E f a π + L/a = 32 . . . . . . . . . t / a . . . . . . E a SP f . . . . . E f a FIG. 18: Fit results for systems of n ∈ { , . . . , } π + mesons for the L/a = 32 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. π + L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP f . . . . . . . . E f a π + L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP f . . . . . E f a π + L/a = 32 . . . . . . . . . t / a E a SP f . . . . . . . E f a π + L/a = 32 . . . . . . . . . t / a E a SP f . . . . . . . . E f a FIG. 19: Fit results for systems of n ∈ { , . . . , } π + mesons for the L/a = 32 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. K L/a = 48 t / a . . . . . E a SP f . . . . . . E f a K L/a = 48 t / a . . . . . . E a SP f . . . . . . . . E f a K L/a = 48 t / a . . . . . . . . E a SP f . . . . . . . E f a K L/a = 48 t / a . . . . . . E a SP f . . . . . . . . E f a FIG. 20: Fit results for systems of n ∈ { , . . . , } K mesons for the L/a = 48 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. K L/a = 48 t / a . . . . . . . . E a SP f . . . . . . E f a K L/a = 48 t / a . . . . . . . . . E a SP f . . . . . . . . E f a K L/a = 48 t / a . . . . . . . . E a SP f . . . . . E f a K L/a = 48 t / a . . . . . E a SP f . . . . . . . . E f a FIG. 21: Fit results for systems of n ∈ { , . . . , } K mesons for the L/a = 48 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. K L/a = 48 t / a . . . . . . E a SP f . . . . . E f a K L/a = 48 t / a . . . . . . E a SP f . . . . . . . . . E f a K L/a = 48 t / a . . . . . . . E a SP f . . . . . E f a K L/a = 48 t / a . . . . . . . . . E a SP f . . . . . . . E f a FIG. 22: Fit results for systems of n ∈ { , . . . , } K mesons for the L/a = 48 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. K L/a = 32 . . . . . . . . . t / a . . . . . E a SP f . . . . E f a K L/a = 32 . . . . . . . . . t / a . . . . . . E a SP f . . . . . . . . E f a K L/a = 32 . . . . . . . . . t / a . . . . . . . . . E a SP f . . . . . . E f a K L/a = 32 . . . . . . . . . t / a . . . . . . E a SP . . . . . . . f . . . . . . . . . E f a FIG. 23: Fit results for systems of n ∈ { , . . . , } K mesons for the L/a = 32 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. K L/a = 32 . . . . . . . . . t / a . . . . . . . . . E a SP . . . . . . . f . . . . . E f a K L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP f . . . . . E f a K L/a = 32 . . . . . . . . . t / a . . . . . E a SP f . . . . . . E f a K L/a = 32 . . . . . . . . . t / a . . . . . . E a SP f . . . . . . E f a FIG. 24: Fit results for systems of n ∈ { , . . . , } K mesons for the L/a = 32 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. K L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP f . . . . . . . . . E f a K L/a = 32 . . . . . . . . . t / a E a SP f . . . . . . E f a K L/a = 32 . . . . . . . . . t / a . . . . E a SP f − − E f a K L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP f − − E f a FIG. 25: Fit results for systems of n ∈ { , . . . , } K mesons for the L/a = 32 lattice volume. Thefigures are analogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. pL/a = 32 t / a . . . . . . . . . E a SPSS f . . . . . . . . . E f a pL/a = 48 t / a . . . . . . . . E a SPSS f . . . . . . . E f a nL/a = 32 t / a . . . . . . . . . E a SPSS f . . . . . . . . . E f a nL/a = 48 t / a . . . . . . . . E a SPSS f . . . . . . E f a FIG. 26: Fit results for proton and neutron systems. The figures are analogous to Fig. 1, seeAppendix B for a definition of the fitting procedure employed. ppL/a = 32 t / a . . . . . E a SPSS f . . . . . . . . E f a ppL/a = 48 t / a . . . . . . . . E a SPSS f . . . . . . E f a nnL/a = 32 t / a . . . . . . . . . E a SPSS f . . . . . . . E f a nnL/a = 48 t / a . . . . . . E a SPSS f . . . . E f a FIG. 27: Fit results for systems of two protons and systems of two neutrons. The figures areanalogous to Fig. 1, see Appendix B for a definition of the fitting procedure employed. pn ( S ) L/a = 32 t / a . . . . . E a SPSS f . . . . . . . . E f a pn ( S ) L/a = 48 t / a . . . . . . E a SPSS f . . . . . . E f a pn ( S ) L/a = 32 t / a . . . . . . . . . E a SPSS f . . . . . . E f a pn ( S ) L/a = 48 t / a . . . . . . E a SPSS f . . . . . E f a FIG. 28: Fit results for systems of one neutron and one proton in the cubic irreps associatedwith S and S systems in the infinite-volume limit. The figures are analogous to Fig. 1, seeAppendix B for a definition of the fitting procedure employed. H L/a = 32 . . . . . . . . . t / a . . . . . . . . E a SP f . . . . . . . . E f a H L/a = 48 t / a . . . . . . . . . E a SP f . . . . . . . . E f a He L/a = 32 . . . . . . . . t / a . . . . . . E a SP f . . . . . . E f a He L/a = 48 t / a . . . . . . . E a SP . . . . . . . f . . . . . . . E f a FIG. 29: Fit results for three-nucleon systems H and He. The figures are analogous to Fig. 1, seeAppendix B for a definition of the fitting procedure employed. [20] P. Schwerdtfeger, L. F. Pa s teka, A. Punnett, and P. O. Bowman, Nuclear Physics A , 551 (2015), ISSN 0375-9474, special Issue on Superheavy Elements, URL .[21] W. Detmold, R. G. Edwards, J. J. Dudek, M. Engelhardt, H.-W. Lin, S. Meinel, K. Orginos,and P. Shanahan (USQCD), Eur. Phys. J. A55 , 193 (2019), 1904.09512.[22] A. Duncan, E. Eichten, and H. Thacker, Phys. Rev. Lett. , 3894 (1996), hep-lat/9602005.[23] A. Duncan, E. Eichten, and R. Sedgewick, Phys. Rev. D71 , 094509 (2005), hep-lat/0405014.[24] T. Blum, T. Doi, M. Hayakawa, T. Izubuchi, and N. Yamada, Phys. Rev.
D76 , 114508 (2007),0708.0484.[25] M. Hayakawa and S. Uno, Prog. Theor. Phys. , 413 (2008), 0804.2044.[26] T. Blum, R. Zhou, T. Doi, M. Hayakawa, T. Izubuchi, S. Uno, and N. Yamada, Phys. Rev.
D82 , 094508 (2010), 1006.1311.[27] S. Aoki et al., Phys. Rev.
D86 , 034507 (2012), 1205.2961.[28] N. Tantalo, PoS
LATTICE2013 , 007 (2014), 1311.2797.[29] Z. Davoudi and M. J. Savage, Phys. Rev.
D90 , 054503 (2014), 1402.6741.[30] M. G. Endres, A. Shindler, B. C. Tiburzi, and A. Walker-Loud, Phys. Rev. Lett. , 072002(2016), 1507.08916.[31] Z. Fodor, C. Hoelbling, S. D. Katz, L. Lellouch, A. Portelli, K. K. Szabo, and B. C. Toth,Phys. Lett.
B755 , 245 (2016), 1502.06921.[32] B. Lucini, A. Patella, A. Ramos, and N. Tantalo, JHEP , 076 (2016), 1509.01636.[33] Z. Fodor, C. Hoelbling, S. Krieg, L. Lellouch, T. Lippert, A. Portelli, A. Sastre, K. K. Szabo,and L. Varnhorst, Phys. Rev. Lett. , 082001 (2016), 1604.07112.[34] M. Hansen, B. Lucini, A. Patella, and N. Tantalo, JHEP , 146 (2018), 1802.05474.[35] M. E. Matzelle and B. C. Tiburzi, Phys. Rev. D95 , 094510 (2017), 1702.01296.[36] A. Patella, PoS
LATTICE2016 , 020 (2017), 1702.03857.[37] Z. Davoudi, J. Harrison, A. Jüttner, A. Portelli, and M. J. Savage, Phys. Rev.
D99 , 034510(2019), 1810.05923.[38] J.-W. Lee and B. C. Tiburzi, Phys. Rev.
D93 , 034012 (2016), 1508.04165.[39] M. Lüscher, Commun. Math. Phys. , 153 (1986).[40] M. Lüscher, Nucl. Phys.
B354 , 531 (1991).[41] R. A. Briceño, J. J. Dudek, and R. D. Young, Rev. Mod. Phys. , 025001 (2018), 1706.06223.[42] S. R. Beane and M. J. Savage, Phys. Rev. D90 , 074511 (2014), 1407.4846.[43] K. Symanzik, Nucl. Phys.
B226 , 187 (1983).[44] M. Lüscher and P. Weisz, Commun. Math. Phys. , 59 (1985), [Erratum: Commun. Math.Phys.98,433(1985)].[45] W. Bietenholz et al., Phys. Rev. D84 , 054509 (2011), 1102.5300.[46] C. Morningstar and M. J. Peardon, Phys. Rev.
D69 , 054501 (2004), hep-lat/0311018.[47] B. Sheikholeslami and R. Wohlert, Nucl. Phys.
B259 , 572 (1985).[48] N. Cundy et al., Phys. Rev.
D79 , 094507 (2009), 0901.3302.[49] R. F. Dashen, Phys. Rev. , 1245 (1969).[50] M. Göckeler, R. Horsley, P. E. L. Rakow, G. Schierholz, and R. Sommer, Nucl. Phys.
B371 ,713 (1992).[51] C. R. Allton et al. (UKQCD), Phys. Rev.
D47 , 5128 (1993), hep-lat/9303009.[52] S. R. Beane, W. Detmold, T. C. Luu, K. Orginos, M. J. Savage, and A. Torok, Phys. Rev.Lett. , 082004 (2008), 0710.1827.[53] W. Detmold, M. J. Savage, A. Torok, S. R. Beane, T. C. Luu, K. Orginos, and A. Parreño, Phys. Rev.
D78 , 014507 (2008), 0803.2728.[54] W. Detmold and M. J. Savage, Phys. Rev.
D82 , 014511 (2010), 1001.2768.[55] W. Detmold and B. Smigielski, Phys. Rev.
D84 , 014508 (2011), 1103.4362.[56] W. Detmold, K. Orginos, and Z. Shi, Phys. Rev.
D86 , 054507 (2012), 1205.4224.[57] W. Detmold and K. Orginos, Phys. Rev.
D87 , 114512 (2013), 1207.1452.[58] S. R. Beane, E. Chang, S. D. Cohen, W. Detmold, H. W. Lin, T. C. Luu, K. Orginos,A. Parreño, M. J. Savage, and A. Walker-Loud (NPLQCD), Phys. Rev.
D87 , 034506 (2013),1206.5219.[59] C. Stein, in
Proceedings of the Third Berkeley Symposium on Mathematical Statistics andProbability, Volume 1: Contributions to the Theory of Statistics (University of Califor-nia Press, Berkeley, Calif., 1956), pp. 197–206, URL https://projecteuclid.org/euclid.bsmsp/1200501656 .[60] O. Ledoit and M. Wolf, Journal of Multivariate Analysis , 365 (2004), ISSN 0047-259X,URL .[61] E. Rinaldi, S. Syritsyn, M. L. Wagman, M. I. Buchoff, C. Schroeder, and J. Wasem, Phys.Rev. D99 , 074510 (2019), 1901.07519.[62] U. van Kolck, Nucl. Phys.
A645 , 273 (1999), nucl-th/9808007.[63] J.-W. Chen, G. Rupak, and M. J. Savage, Nucl. Phys.
A653 , 386 (1999), nucl-th/9902056.[64] H. A. Bethe, Phys. Rev. , 38 (1949).[65] X. Kong and F. Ravndal, Phys. Lett. B450 , 320 (1999), [Erratum: Phys.Lett.B458,565(1999)], nucl-th/9811076.[66] X. Kong and F. Ravndal, Nucl. Phys.
A665 , 137 (2000), hep-ph/9903523.[67] K. Huang and C. N. Yang, Phys. Rev. , 767 (1957).[68] S. R. Beane, W. Detmold, and M. J. Savage, Phys. Rev.
D76 , 074507 (2007), 0707.1670.[69] W. Detmold and M. J. Savage, Phys. Rev.
D77 , 057502 (2008), 0801.0763.[70] M. H. Namjoo, A. H. Guth, and D. I. Kaiser, Phys. Rev.
D98 , 016011 (2018), 1712.00445.[71] E. Braaten, A. Mohapatra, and H. Zhang, Phys. Rev.
D98 , 096012 (2018), 1806.01898.[72] M. T. Hansen and S. R. Sharpe, Phys. Rev.
D93 , 014506 (2016), 1509.07929.[73] W. E. Caswell and G. P. Lepage, Phys. Lett. , 437 (1986).[74] A. Pineda and J. Soto, Nucl. Phys. Proc. Suppl. , 428 (1998), [,428(1997)], hep-ph/9707481.[75] A. Pineda and J. Soto, Phys. Rev. D59 , 016005 (1999), hep-ph/9805424.[76] M. Göckeler, R. Horsley, M. Lage, U.-G. Meissner, P. Rakow, A. Rusetsky, G. Schierholz,and J. Zanotti, Phys.Rev.D , 094513 (2012), 1206.4141.[77] L. Leskovec and S. Prelovsek, Phys. Rev. D85 , 114507 (2012), 1202.2145.[78] M. T. Hansen and S. R. Sharpe, Phys. Rev.
D93 , 096006 (2016), [Erratum: Phys.Rev.D96,no.3,039901(2017)], 1602.00324.[79] T. Mehen and I. W. Stewart, Nucl. Phys.
A665 , 164 (2000), nucl-th/9901064.[80] S. R. Beane, W. Detmold, H.-W. Lin, T. C. Luu, K. Orginos, M. J. Savage, A. Torok, andA. Walker-Loud (NPLQCD), Phys. Rev.
D81 , 054505 (2010), 0912.4243.[81] S. R. Beane, E. Chang, W. Detmold, H. W. Lin, T. C. Luu, K. Orginos, A. Parreño, M. J. Sav-age, A. Torok, and A. Walker-Loud (NPLQCD), Phys. Rev.
D85 , 054511 (2012), 1109.2889.[82] T. Yamazaki, K.-i. Ishikawa, Y. Kuramashi, and A. Ukawa, Phys. Rev.
D86 , 074514 (2012),1207.4277.[83] S. R. Beane et al. (NPLQCD), Phys. Rev.
C88 , 024003 (2013), 1301.5790.[84] K. Orginos, A. Parreño, M. J. Savage, S. R. Beane, E. Chang, and W. Detmold, Phys. Rev.
D92 , 114512 (2015), 1508.07583. [85] E. Berkowitz, T. Kurth, A. Nicholson, B. Joó, E. Rinaldi, M. Strother, P. M. Vranas, andA. Walker-Loud, Phys. Lett. B765 , 285 (2017), 1508.00886.[86] T. Yamazaki, K.-i. Ishikawa, Y. Kuramashi, and A. Ukawa, Phys. Rev.
D92 , 014501 (2015),1502.04182.[87] M. L. Wagman, F. Winter, E. Chang, Z. Davoudi, W. Detmold, K. Orginos, M. J. Savage,and P. E. Shanahan, Phys. Rev.
D96 , 114510 (2017), 1706.06550.[88] T. D. Blanton, F. Romero-López, and S. R. Sharpe (2019), 1909.02973.[89] J. Bijnens, G. Colangelo, G. Ecker, J. Gasser, and M. E. Sainio, Nucl. Phys.
B508 , 263(1997), [Erratum: Nucl. Phys.B517,639(1998)], hep-ph/9707291.[90] M. Tanabashi et al. (Particle Data Group), Phys. Rev.
D98 , 030001 (2018).[91] T. R. Haar, Y. Nakamura, and H. Stüben, EPJ Web Conf. , 14011 (2018), 1711.03836.[92] R. G. Edwards and B. Joó (SciDAC, LHPC, UKQCD), Nucl. Phys. Proc. Suppl. , 832(2005), [,832(2004)], hep-lat/0409003.[93] G. Parisi, Phys. Rept. , 203 (1984).[94] G. P. Lepage, in
Boulder TASI 1989:97-120 (1989), pp. 97–120.[95] M. L. Wagman and M. J. Savage, Phys. Rev.
D96 , 114508 (2017), 1611.07643.[96] M. Lüscher and P. Weisz, Nucl. Phys.
B240 , 349 (1984).[97] G. Golub and V. Pereyra, Inverse Problems , R1 (2003), URL http://stacks.iop.org/0266-5611/19/i=2/a=201 .[98] D. P. O’Leary and B. W. Rust, Computational Optimization and Applications , 579 (2013).[99] H. Akaike, IEEE Transactions on Automatic Control , 716 (1974), ISSN 2334-3303.[100] P. Young, Everything you wanted to know about data analysis and fitting but were afraid toask (2012), 1210.3781.[101] A. C. Davison and D. V. Hinkley,
The Basic Bootstraps (Cambridge University Press, 1997),p. 11–69, Cambridge Series in Statistical and Probabilistic Mathematics.[102] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Review , 65 (2017).[103] P. K. Mogensen and A. N. Riseth, Journal of Open Source Software3