The energy-dependent π + π + π + scattering amplitude from QCD
Maxwell T. Hansen, Raul A. Briceño, Robert G. Edwards, Christopher E. Thomas, David J. Wilson
CCERN-TH-2020-147, JLAB-THY-20-3242
The energy-dependent π + π + π + scattering amplitude from QCD Maxwell T. Hansen, ∗ Raul A. Brice˜no,
2, 3, † Robert G. Edwards, ‡ Christopher E. Thomas, § and David J. Wilson ¶ (for the Hadron Spectrum Collaboration) Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland Thomas Jefferson National Accelerator Facility, 12000 Jefferson Avenue, Newport News, VA 23606, USA Department of Physics, Old Dominion University, Norfolk, VA 23529, USA DAMTP, University of Cambridge, Wilberforce Road, Cambridge, CB3 0WA, UK (Dated: September 11, 2020)Focusing on three-pion states with maximal isospin ( π + π + π + ), we present the first non-perturbativedetermination of an energy-dependent three-hadron scattering amplitude from first-principles QCD.The calculation combines finite-volume three-hadron energies, extracted using numerical lattice QCD,with a relativistic finite-volume formalism, required to interpret the results. To fully implementthe latter, we also solve integral equations that relate an intermediate three-body K matrix to thephysical three-hadron scattering amplitude. The resulting amplitude shows rich analytic structuresand a complicated dependence on the two-pion invariant masses, represented here via Dalitz-likeplots of the scattering rate. Introduction — The three-body problem lies at the coreof a broad range of outstanding questions in quantumchromodynamics (QCD), from the determination of three-body nuclear forces to the nature of the Roper resonance,the lowest-lying excitation of the nucleon. In particular,the largest systematic uncertainty in QCD-based structurecalculations of light nuclei is the estimate of the three-nucleon force (see, for example, Ref. [1]). More gener-ally, QCD resonances often couple strongly to asymptoticstates containing three or more hadrons. For example, theRoper has defied simple quark-model descriptions, duein part to its nature as a broad resonance with a ∼ N ππ . A rigorous QCD calculationwould elucidate the role of non-perturbative dynamicsin the Roper’s peculiar properties, e.g. the fact that ithas a lower mass than the negative-parity ground state,which seems unnatural from the perspective of the quarkmodel [2, 3].As a necessary step towards studying a broad class ofhadronic systems with significant three-hadron branchingfractions, in this work we present the first study of anenergy-dependent three-body scattering amplitude fromQCD. This non-perturbative result is achieved by thecoalescence of three novel techniques: a calculation offinite-volume three-hadron energies based in numericallattice QCD using the distillation method [4], a relativis-tic finite-volume formalism to relate the energies to Kmatrices, and a numerical evaluation of corresponding in-tegral equations to convert the latter into the three-hadronscattering amplitude. The theoretical basis required toachieve these final two steps was derived in Refs. [5, 6]. This work considers the scattering of three-pion states ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] A large body of formal work has gone in to understanding with maximal isospin ( I = 3) in QCD with three dynami-cal quarks ( N f = 2+1): two degenerate light quarks, withheavier-than-physical mass corresponding to a pion mass m π ≈
391 MeV, and a strange quark. The calculation isperformed on two lattice gauge ensembles, differing onlyin their spatial and temporal extents.Note that both the maximal-isospin three-pion systemand its two-pion subsystem are expected to be weaklyinteracting and non-resonant. This is in contrast to I = 2three-pion states, for which the two-pion subsystem formsthe ρ resonance, as well as to I = 0 and I = 1 three-pionstates, for which both two- and three-pion resonances areexpected. The formal framework to describe all such sys-tems using lattice QCD has been worked out in Ref. [43].Many numerical studies of three-hadron states havebeen published over the last decade, ranging from earlyperturbative work describing the three-pion ground state[46] to more recent results using quantization conditionsto study ground [47] and excited states [48–50], with thelatter set each analyzing the lattice QCD spectrum pub-lished in Ref. [51]. Independent sets of finite-volume ener-gies have also been calculated and analyzed in Refs. [52]and [53], with the latter making the impressive step ofcalculating at the physical value of the pion mass. Thepresent investigation goes beyond this previous work, byproviding the first complete numerical determination ofphysical scattering amplitudes for three-body systems.In the following, we first discuss our determinationof two- and three-pion finite-volume energies, before de-scribing the fits used to relate these to infinite-volumeK matrices. The latter then serve as inputs to knownintegral equations, which we solve numerically to extractthe 3 π + → π + scattering amplitude. Various details ofthe analysis are presented in the supplementary material. general methods for relating finite-volume energies to scatteringamplitudes for both two- and three-body states. See Refs. [7–25]and Refs. [26–45], respectively. a r X i v : . [ h e p - l a t ] S e p FIG. 1. The π + π + and π + π + π + finite-volume spectra in thecenter-of-momentum frame for the relevant finite-volume irrepswith various overall momenta, as explained in the text. Pointsare computed energy levels on the two volumes with errorbars showing statistical uncertainties. Each rectangular insertshows a vertical zoom of the region indicated by the smallneighboring rectangle. Grey curves are the “non-interacting”finite-volume energies, i.e. the energies in the absence of anyinteractions between pions. Orange curves are predictions fromthe finite-volume formalism based only on the two-particlescattering length, given in Eq. (4) (here with the local three-body interaction set to zero). Spectral Determination — Figure 1 summarizes the two-and three-pion finite-volume spectra calculated in thiswork. The figure also shows the “non-interacting” finite-volume energies in the absence of interactions betweenpions and the predictions from the finite-volume formalismdiscussed later.Computations were performed on anisotropic latticeswhich have a temporal lattice spacing, a t , finer thanthe spatial lattice spacing, a s ( a t = a s /ξ with ξ =3 . Two-pion energies on the larger volume have already appeared inRef. [54]. level Symanzik-improved gauge action. The bare pa-rameters and basic lattice properties are detailed inRefs. [55, 56], and some specific details of the ensemblesused in this work are summarized in Table I. Setting thescale via a − t = m expΩ ( a t m lattΩ ) − , and combining with a t m π = 0 . a t m K = 0 . m π ≈
391 MeV and m K ≈
550 MeV.The spectrum of energies in a finite volume is discreteand each energy level provides a constraint on the scatter-ing amplitudes at the corresponding center-of-momentumenergy. To obtain more constraints, we perform calcu-lations where the the two-pion and three-pion systemshave overall non-zero momentum, P , as well as for P = .Momenta are quantized by the cubic spatial boundaryconditions, P = πL ( n , n , n ), where { n i } are integers,and we write this using a shorthand notation as [ n n n ].In this work we restrict attention to S -wave scattering.The reduced symmetry of a cubic lattice means that totalangular momentum, J , is not a good quantum numberand instead channels are labelled by the irreducible repre-sentation (irrep, Λ) of the octahedral group with parity for P = or the relevant subgroup that leaves P invariantfor P (cid:54) = [60, 61]. We consider the relevant irreps whichcontain J = 0: A − ( A +1 ) for πππ ( ππ ) at rest and A ( A )for πππ ( ππ ) with non-zero P . Isospin, I , and G -parity, G , are good quantum numbers in our lattice formulation;these distinguish the two-pion ( I G = 2 + ) and three-pion( I G = 3 − ) channels. We neglect higher partial waves here,in particular the two-particle D -wave which mixes withthe S -wave in the finite-volume energies. As described inRef. [54], a nonzero D -wave interaction can be extracted,in particular if aided by the consideration of other, non-trivial finite-volume irreps. However, as can also be seenfrom the earlier work, this has a small influence on the A +1 and A two-pion energies considered here. To reliably extract the discrete set of energy levelsin the finite volumes considered here, in particular themany excited energies in each channel, we have com-puted two-point correlation functions using a large basisof appropriate creation and annihilation interpolatingoperators. From these, the finite-volume spectra were de-termined using the variational method [62–64], with ourimplementation described in Refs. [65, 66]. This amountsto calculating a matrix of correlation functions, G ij ( t ) = (cid:104)O i ( t ) O † j (0) (cid:105) , (1)and diagonalizing M ( t, t ) = G ( t ) − / · G ( t ) · G ( t ) − / for a fixed t . One can show that the corresponding eigen-values satisfy λ n ( t, t ) → e − E n ( L )( t − t ) , where E n ( L ) is where a t m lattΩ = 0 . m expΩ is the experimentally determined Ω baryon mass from Ref. [58] There is nonetheless, in principle, a systematic uncertainty as-sociated with neglecting the D -wave contribution. Given ourconsistency with the results of Ref. [54], this appears to be well-below the statistical uncertainty in the present fits. See alsoSecs. VIII A and B of that work for more discussion. volume m π L N cfg N vec n [000]3 π n [001]3 π n [011]3 π n [111]3 π ×
256 4.76 256 128 8 8 11 724 ×
128 5.71 553 162 10 16 20 7TABLE I. Summary of lattice ensemble volumes (presented as(
L/a s ) × T /a t ), number of gauge-field configurations ( N cfg ),number of distillation vectors ( N vec ), and maximum numberof operators used for each total momentum ( n P π ). the n ’th energy level with overlap to some of the operatorsin the basis. This basic methodology has been appliedto a wide range of two-hadron scattering observables forseveral phenomenologically interesting channels [54, 67–78]. See Sec. 1 of the supplementary material for someexample plots of λ n ( t, t ).In order to robustly interpolate the two- and three-pionenergy eigenstates we use operators with two- and three-meson-like structures in the appropriate irrep, constructedfrom products of single-meson-like operators projected todefinite spatial momentum. The latter are built from lin-ear combinations, chosen to optimize overlap to the single-pion states, of fermion bilinears of the form, ¯ ψ Γ D . . . Dψ ,where ψ is a quark field and D is a discretized covariantderivative. Details of these operator constructions aregiven in Sec. 5 of the supplementary material, and thetwo-pion and three-pion operators are listed in respec-tively Table VI and Tables VII - XI. We also summarizethe maximum number of operators considered for thethree-pion spectra in Table I. Using such a wide varietyof optimized operators, and especially multi-hadron opera-tors with momentum-projected single-hadron components,allows one to minimize excited state contamination andextract the energies reliably and precisely from smallvalues of t . This approach is made feasible due to thedistillation method [4] which we employ to efficiently com-pute the numerous quark-field Wick contractions that arerequired.Returning to the two- and three-pion spectra summa-rized in Fig. 1, we observe a one-to-one correspondence be-tween the computed energy levels and the non-interactingenergies in all panels, with the computed values slightlyhigher in energy than the non-interacting levels. This sug-gests that the system is weakly interacting and repulsivein both the two- and three-hadron sectors. Analyzing the finite-volume spectra — We now describeour method for determining two- and three-body K matri-ces from the extracted finite-volume energies, beginningwith an overview of scattering observables:The two-pion scattering amplitude is defined as theconnected part of the overlap between an incoming π + π + In extracting the energies, subsets of the full operator set werealso considered, in order to investigate sensitivity to the detailedchoice of operator basis.
FIG. 2. Example of data and fits for K and K , iso , as describedin the text. The red points are given by substituting finite-volume energies into − /F ( E , P , L ) and − /F , iso ( E , P , L )for the two- and three-particle energies, respectively, withthe volume and P value indicated in the legend. A symbolappearing at the very top or bottom represents a case wherethe central value falls outside the plotted region. The darkcyan bands represent the fit shown in Eq. (4) while the lighterbands show the spread covered by the various fits described inthe supplemental material. For the bottom panel we normalizeto m π K LO3 , iso = 4608 π ( m π a ) , with m π a taken from Eq. (4).This simple relation between K , iso and the two-particle scat-tering length holds at leading order in chiral perturbationtheory, as was first derived in Ref. [48]. asymptotic state (with momenta p , − p ) to an outgoing π + π + state (with p (cid:48) , − p (cid:48) ). Without loss of generality,here we have assumed the center-of-momentum frame.We also define p = | p | = | p (cid:48) | , where we have used that themagnitudes must be equal to satisfy energy conservation.In addition, s ≡ E (cid:63) ≡ p + m π ) defines the squaredcenter-of-momentum frame energy. The only remainingdegree of freedom is the scattering angle between p and p (cid:48) .In this work we focus on the S -wave scattering amplitude,denoted M , in which this angle is integrated to projectonto zero-angular-momentum states. Finally we recallthe simple relation between M and the K matrix in theelastic region, K − = Re M − . In contrast to M , K is real for real s and is meremorphic in a region of thecomplex s plane around s = 4 m π . In this work we alsoconsider an analogous, three-body K matrix, introduced In addition, the imaginary part of M − is completely fixed byunitarity so that K is the only part free to depend on thedynamics of the system. We work with the simple phase spacefactor, proportional to the momentum magnitude. See, e.g.,Ref. [24] for more details. in Ref. [5] and denoted by K df , .In the two-pion sector, in the case that the S -waveinteractions are dominant, the scalar-irrep finite-volumeenergies satisfy the quantization condition [7, 8, 10], K ( E (cid:63) ) + F − ( E , P , L ) = 0 , (2)where E (cid:63) ≡ (cid:112) E − P is the center-of-momentum en-ergy and F ( E , P , L ) is a known geometric function. Forthe three-body sector, we use the isotropic approximationof the general formalism derived in Ref [5], which takesan analogous form, now for pseudoscalar-irrep energies K , iso ( E (cid:63) ) + F − , iso [ K ]( E , P , L ) = 0 , (3)where the notation is meant to stress that F , iso [ K ]( E , P , L ) is a functional of K ( E (cid:63) ). F , iso is defined in Eq. (39) of Ref. [5]. Here K , iso is thecomponent of K df , that only depends on the totalthree-hadron energy, i.e. is “isotropic”. Equation (3)holds only when K df , is well approximated to beisotropic and our fits give evidence that this is a goodapproximation for this system.Combining these two conditions with the energies plot-ted in Fig. 1 allows one to constrain both the two- andthree-hadron K matrices. One possible strategy is to fita parametrization of K and use this to then determinethe energy dependence of K , iso . This is summarized inFig. 2. An alternative approach is to parametrize both Kmatrices and fit these simultaneously to the entire set offinite-volume energies. A detailed discussion with a widerange of fits is given in Sec. 2 of the supplementary ma-terial. Both strategies give consistent results and the keymessage is that the full data set is well described by a con-stant K , iso that is consistent with zero, together with theleading-order effective range expansion: tan δ ( p ) = − a p with K ( E (cid:63) ) = − πE (cid:63) tan δ ( p ) /p . Here the second equa-tion defines the S -wave scattering phase shift, δ ( p ), andthe first defines the scattering length, a . Our best fit, per-formed simultaneously to all spectra shown in Fig. 1 butwith a cutoff in the center-of-momentum frame energiesincluded, yields m π a = 0 . ± . m π K , iso = − ± (cid:20) . . . (cid:21) , (4)with a χ per degree-of-freedom of 64 . / (37 −
2) = 1 . This fit is denoted by B in Sec. 2 of the supplementarymaterial. As explained there, the fitted data includes all two-pion energies below E (cid:63) , cut = 3 . m π and all three-pion energiesbelow E (cid:63) , cut = 4 . m π , with both cutoffs applied to energies inthe center-of-momentum frame. with K , iso = 0). In Fig. 2 we illustrate the same fit usingthe darker cyan curves. In addition, we include the lighterbands as a systematic uncertainty, estimated from thespread of various constant and linear fits, as detailed inSec. 2 of the supplementary material.3 π + scattering amplitude — Following the relativisticintegral equations presented in Ref. [6], we can write the J = 0 and K , iso = 0 amplitude as follows: M ( u,u )3 ( p, k ) = −M ( E (cid:63) ,p ) G s ( p, k ) M ( E (cid:63) ,k ) − M ( E (cid:63) ,p ) (cid:90) k (cid:48) G s ( p, k (cid:48) ) M ( u,u )3 ( k (cid:48) , k ) , (5)where (cid:82) k ≡ (cid:82) dk k / [(2 π ) ω k ] and we have introduced G s ( p, k ) ≡ − H ( p, k )4 pk log (cid:20) α ( p, k ) − pk + i(cid:15)α ( p, k ) + 2 pk + i(cid:15) (cid:21) , (6) α ( p, k ) ≡ ( E − ω k − ω p ) − p − k − m . (7) M is the S -wave two-particle scattering amplitude,introduced above, which depends on the invariant E (cid:63) ,k ≡ ( E − ω k ) − k , with ω k = √ k + m . The ( u, u ) super-script emphasizes that specific spectator momenta, k and p , are singled out in the initial and final states respectively.For example k is the magnitude of the three-momentumfor the incoming pion that is disconnected from the first M insertion. The function G s encodes the spectatorexchange, projected to the S -wave. It inherits a schemedependence through the smooth cutoff function H , de-fined in Eqs. (28) and (29) of Ref. [5] [and in Eq. (30)of Sec. 3 in the supplementary material]. This schemedependence is matched by that inside of K , iso such thatthe resulting scattering amplitude is universal.To use Eq. (5) in practice, one requires a parameteri-zation for M . As described in the previous section, the π + π + system is well described using the leading ordereffective range expansion for M , M ( E (cid:63) ) = 16 πE (cid:63) − /a − i (cid:112) E (cid:63) / − m π . (8)Following the derivation of Ref. [6], the final step is tosymmetrize with respect to the spectators, to reach M ( s , m (cid:48) , m (cid:48) , m , m ) = (cid:88) p i ∈P p (cid:88) k ∈P k M ( u,u )3 ( p, k ) , (9)where P p = { p , a (cid:48) , − p − a (cid:48) } and P k = { k , a , − k − a } . Wehave presented the left-hand side as a function of the fiveLorentz invariants that survive after truncating to J = 0in both the two and three particle sector: the squaredthree-hadron center-of-momentum frame energy, s , aswell as two pion-pair invariant masses for each of the initialand final states. These are defined by introducing thenotation { k , a , − k − a } = { p , p , p } , then for example m = ( p + p ) = ( E (cid:63) − [ m π + p ] / ) − p , (10) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m / m ⇡ m / m ⇡ m /m ⇡ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ ⇥ FIG. 3.
Top:
Dalitz-like plot of m π |M | for √ s = 3 . m withfinal kinematics fixed to { p (cid:48) , p (cid:48) } = { . m π , . m π } = ⇒{ m (cid:48) , m (cid:48) } = { . m π , . m π } . Bottom:
Same total energy,now with incoming and outgoing kinematics set equal, asdiscussed in the text. where the middle expression depends on on-shell four-vectors with p = m π .In the top panel of Fig. 3 we show a Dalitz-like plotof |M | as a function of ( m , m ), with all other kine-matics fixed as indicated in the caption. Note that, ina usual Dalitz description, the incoming energy is fixedby the decaying particle so that only the outgoing kine-matics can vary, whereas here we fix the other kinematicsto some chosen values. The inputs to this plot are thebest-fit scattering length, given in Eq. (4), together with K , iso = 0. The bottom panel of Fig. 3 shows the same √ s but varies incoming and outgoing kinematics accord-ing to m = m (cid:48) and m = m (cid:48) . These kinematicsmay be instructive in future resonant analyses since theylead to a residue at the pole that can be roughly inter-preted as an overlap of the resonance state to a particularthree-hadron state.Additional details concerning the S -wave integral equa- tions, their numerical solutions, and the properties of theresulting amplitude are presented in Secs. 3 and 4 ofthe supplementary material, where we also describe thepropagation of the uncertainties of m π a and K , iso intothe predicted scattering amplitude. Summary — In this work we have presented the first lat-tice QCD determination of the energy-dependent three-to-three scattering amplitude for three pions with maximalisospin. The calculation proceeded in three steps: (i) determining finite-volume energies with π + π + π + quan-tum numbers, (ii) using the framework of Ref. [5] toextract two- and three-body K matrices from these, and (iii) applying the results of Ref. [6] to convert these tothe three-hadron scattering amplitude, by solving knownintegral equations. The three steps are summarized, re-spectively, by Figs. 1, 2 and 3 of the text.Having established this general workflow, it is now wellwithin reach to rigorously extract three-hadron resonanceproperties from lattice QCD calculations. In particularthe formalism has recently been extended to three-pionstates with any value of isospin in Ref. [43]. This shouldenable studies, for example, of the ω , h and a resonances.The main outstanding challenges here include rigorousresonant parametrizations of the intermediate three-bodyK matrix, as well as a better understanding of the ana-lytic continuation required to identify the resonance poleposition. ACKNOWLEDGMENTS
The authors would like to thank M. Bruno, J. Dudek,A. Jackura, L. Leskovec and A. Rodas for useful con-versations, as well as our other colleagues within theHadron Spectrum Collaboration. RAB, RGE and CETacknowledge support from the U.S. Department of En-ergy contract DE-AC05-06OR23177, under which Jef-ferson Science Associates, LLC, manages and operatesJefferson Lab. RAB also acknowledges support from theUSDOE Early Career award, contract de-sc0019229. CETand DJW acknowledge support from the U.K. Scienceand Technology Facilities Council (STFC) [grant num-ber ST/P000681/1]. DJW acknowledges support froma Royal Society University Research Fellowship. MTH,CET and DJW acknowledge the MITP topical workshop“Scattering Amplitudes and Resonance Properties for Lat-tice QCD” for stimulating this project and for hospitalityduring the initial discussions. MTH and CET also ac-knowledge the CERN-TH Institute “Advances in LatticeGauge Theory”, which provided the opportunity to makesignificant progress on this work. CET further acknowl-edges CERN TH for hospitality and support during a See also Ref. [79] for a discussion of related integral equationsand their solutions in a resonant three-hadron channel. visit in January and February of this year. The softwarecodes
Chroma [80] and
QUDA [81–83] were used. The au-thors acknowledge support from the U.S. Department ofEnergy, Office of Science, Office of Advanced ScientificComputing Research and Office of Nuclear Physics, Sci-entific Discovery through Advanced Computing (SciDAC)program. Also acknowledged is support from the U.S.Department of Energy Exascale Computing Project. Thiswork was also performed on clusters at Jefferson Labunder the USQCD Collaboration and the LQCD ARRAProject. This research was supported in part under anALCC award, and used resources of the Oak Ridge Lead-ership Computing Facility at the Oak Ridge NationalLaboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This research used resources of theNational Energy Research Scientific Computing Center(NERSC), a DOE Office of Science User Facility sup-ported by the Office of Science of the U.S. Department ofEnergy under Contract No. DE-AC02-05CH11231. Theauthors acknowledge the Texas Advanced ComputingCenter (TACC) at The University of Texas at Austin forproviding HPC resources. Gauge configurations were gen-erated using resources awarded from the U.S. Departmentof Energy INCITE program at the Oak Ridge LeadershipComputing Facility, the NERSC, the NSF Teragrid atthe TACC and the Pittsburgh Supercomputer Center, aswell as at Jefferson Lab. [1] M. Piarulli et al. , Phys. Rev. Lett. , 052503 (2018),arXiv:1707.02883 [nucl-th].[2] N. Isgur and G. Karl, Phys. Lett. B , 109 (1977).[3] N. Isgur and G. Karl, Phys. Rev. D , 2653 (1979),[Erratum: Phys.Rev.D 23, 817 (1981)].[4] M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek,R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, andK. J. Juge (Hadron Spectrum), Phys. Rev. D80 , 054506(2009), arXiv:0905.2160 [hep-lat].[5] M. T. Hansen and S. R. Sharpe, Phys. Rev.
D90 , 116003(2014), arXiv:1408.5933 [hep-lat].[6] M. T. Hansen and S. R. Sharpe, Phys. Rev.
D92 , 114509(2015), arXiv:1504.04248 [hep-lat].[7] M. L¨uscher, Nucl. Phys.
B354 , 531 (1991).[8] K. Rummukainen and S. A. Gottlieb, Nucl. Phys.
B450 ,397 (1995), arXiv:hep-lat/9503028 [hep-lat].[9] P. F. Bedaque, Phys. Lett.
B593 , 82 (2004), arXiv:nucl-th/0402051 [nucl-th].[10] C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl.Phys.
B727 , 218 (2005), arXiv:hep-lat/0507006 [hep-lat].[11] Z. Fu, Phys. Rev.
D85 , 014506 (2012), arXiv:1110.0319[hep-lat].[12] L. Leskovec and S. Prelovsek, Phys. Rev.
D85 , 114507(2012), arXiv:1202.2145 [hep-lat].[13] M. G¨ockeler, R. Horsley, M. Lage, U. G. Meißner, P. E. L.Rakow, A. Rusetsky, G. Schierholz, and J. M. Zanotti,Phys. Rev.
D86 , 094513 (2012), arXiv:1206.4141 [hep-lat].[14] S. He, X. Feng, and C. Liu, JHEP , 011 (2005),arXiv:hep-lat/0504019 [hep-lat].[15] M. Lage, U.-G. Meißner, and A. Rusetsky, Phys. Lett. B681 , 439 (2009), arXiv:0905.0069 [hep-lat].[16] V. Bernard, M. Lage, U. G. Meißner, and A. Rusetsky,JHEP , 019 (2011), arXiv:1010.6018 [hep-lat].[17] M. D¨oring, U.-G. Meißner, E. Oset, and A. Rusetsky,Eur. Phys. J. A47 , 139 (2011), arXiv:1107.3988 [hep-lat].[18] M. D¨oring and U. G. Meißner, JHEP , 009 (2012),arXiv:1111.0616 [hep-lat].[19] D. Agadjanov, U.-G. Meißner, and A. Rusetsky, JHEP , 103 (2014), arXiv:1310.7183 [hep-lat].[20] M. D¨oring, U. G. Meißner, E. Oset, and A. Rusetsky,Eur. Phys. J. A48 , 114 (2012), arXiv:1205.4838 [hep-lat].[21] M. T. Hansen and S. R. Sharpe, Phys. Rev.
D86 , 016007(2012), arXiv:1204.0826 [hep-lat].[22] R. A. Brice˜no and Z. Davoudi, Phys. Rev.
D88 , 094507 (2013), arXiv:1204.1110 [hep-lat].[23] P. Guo, J. Dudek, R. Edwards, and A. P. Szczepaniak,Phys. Rev.
D88 , 014501 (2013), arXiv:1211.0929 [hep-lat].[24] R. A. Brice˜no, J. J. Dudek, and R. D. Young, Rev. Mod.Phys. , 025001 (2018), arXiv:1706.06223 [hep-lat].[25] R. A. Brice˜no, Phys. Rev. D89 , 074507 (2014),arXiv:1401.3312 [hep-lat].[26] F. Romero-L´opez, S. R. Sharpe, T. D. Blanton, R. A.Brice˜no, and M. T. Hansen, JHEP , 007 (2019),arXiv:1908.02411 [hep-lat].[27] R. A. Brice˜no, M. T. Hansen, and S. R. Sharpe, Phys.Rev. D99 , 014516 (2019), arXiv:1810.01429 [hep-lat].[28] R. A. Brice˜no, M. T. Hansen, and S. R. Sharpe, Phys.Rev.
D95 , 074510 (2017), arXiv:1701.07465 [hep-lat].[29] M. T. Hansen and S. R. Sharpe, Ann. Rev. Nucl. Part.Sci. , 65 (2019), arXiv:1901.00483 [hep-lat].[30] R. A. Brice˜no, M. T. Hansen, and S. R. Sharpe, Phys.Rev. D98 , 014506 (2018), arXiv:1803.04169 [hep-lat].[31] T. D. Blanton, F. Romero-L´opez, and S. R. Sharpe,JHEP , 106 (2019), arXiv:1901.07095 [hep-lat].[32] H.-W. Hammer, J.-Y. Pang, and A. Rusetsky, JHEP ,109 (2017), arXiv:1706.07700 [hep-lat].[33] H. W. Hammer, J. Y. Pang, and A. Rusetsky, JHEP ,115 (2017), arXiv:1707.02176 [hep-lat].[34] R. A. Brice˜no and Z. Davoudi, Phys. Rev. D87 , 094507(2013), arXiv:1212.3398 [hep-lat].[35] K. Polejaeva and A. Rusetsky, Eur. Phys. J.
A48 , 67(2012), arXiv:1203.1241 [hep-lat].[36] M. Mai and M. D¨oring, Eur. Phys. J.
A53 , 240 (2017),arXiv:1709.08222 [hep-lat].[37] P. Guo and V. Gasparian, Phys. Lett.
B774 , 441 (2017),arXiv:1701.00438 [hep-lat].[38] P. Guo and V. Gasparian, Phys. Rev.
D97 , 014504 (2018),arXiv:1709.08255 [hep-lat].[39] P. Guo and T. Morris, Phys. Rev.
D99 , 014501 (2019),arXiv:1808.07397 [hep-lat].[40] F. Romero-L´opez, A. Rusetsky, and C. Urbach, Eur.Phys. J.
C78 , 846 (2018), arXiv:1806.02367 [hep-lat].[41] A. W. Jackura, S. M. Dawid, C. Fern´andez-Ram´ırez,V. Mathieu, M. Mikhasenko, A. Pilloni, S. R. Sharpe,and A. P. Szczepaniak, Phys. Rev.
D100 , 034508 (2019),arXiv:1905.12007 [hep-ph].[42] R. A. Brice˜no, M. T. Hansen, S. R. Sharpe, andA. P. Szczepaniak, Phys. Rev.
D100 , 054508 (2019), arXiv:1905.11188 [hep-lat].[43] M. T. Hansen, F. Romero-L´opez, and S. R. Sharpe, JHEP , 047 (2020), arXiv:2003.10974 [hep-lat].[44] T. D. Blanton and S. R. Sharpe, (2020), arXiv:2007.16190[hep-lat].[45] T. D. Blanton and S. R. Sharpe, (2020), arXiv:2007.16188[hep-lat].[46] W. Detmold, M. J. Savage, A. Torok, S. R. Beane, T. C.Luu, K. Orginos, and A. Parreno, Phys. Rev. D78 , 014507(2008), arXiv:0803.2728 [hep-lat].[47] M. Mai and M. D¨oring, Phys. Rev. Lett. , 062503(2019), arXiv:1807.04746 [hep-lat].[48] T. D. Blanton, F. Romero-L´opez, and S. R. Sharpe,Phys. Rev. Lett. , 032001 (2020), arXiv:1909.02973[hep-lat].[49] M. Mai, M. D¨oring, C. Culver, and A. Alexandru, Phys.Rev. D , 054510 (2020), arXiv:1909.05749 [hep-lat].[50] P. Guo and B. Long, Phys. Rev. D , 094510 (2020),arXiv:2002.09266 [hep-lat].[51] B. H¨orz and A. Hanlon, Phys. Rev. Lett. , 142002(2019), arXiv:1905.04277 [hep-lat].[52] C. Culver, M. Mai, R. Brett, A. Alexandru, andM. D¨oring, Phys. Rev. D , 114507 (2020),arXiv:1911.09047 [hep-lat].[53] M. Fischer, B. Kostrzewa, L. Liu, F. Romero-L´opez,M. Ueding, and C. Urbach, (2020), arXiv:2008.03035[hep-lat].[54] J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys.Rev.
D86 , 034031 (2012), arXiv:1203.6041 [hep-ph].[55] R. G. Edwards, B. Joo, and H.-W. Lin, Phys. Rev.
D78 ,054501 (2008), arXiv:0803.3960 [hep-lat].[56] H.-W. Lin et al. (Hadron Spectrum), Phys. Rev.
D79 ,034502 (2009), arXiv:0810.3588 [hep-lat].[57] R. G. Edwards, J. J. Dudek, D. G. Richards, and S. J.Wallace, Phys. Rev.
D84 , 074508 (2011), arXiv:1104.5152[hep-ph].[58] P. Zyla et al. (Particle Data Group), PTEP , 083C01(2020).[59] D. J. Wilson, J. J. Dudek, R. G. Edwards, and C. E.Thomas, Phys. Rev.
D91 , 054008 (2015), arXiv:1411.2004[hep-ph].[60] R. C. Johnson, Phys. Lett. , 147 (1982).[61] D. C. Moore and G. T. Fleming, Phys. Rev.
D73 ,014504 (2006), [Erratum: Phys. Rev.D74,079905(2006)],arXiv:hep-lat/0507018 [hep-lat].[62] C. Michael, Nucl. Phys.
B259 , 58 (1985).[63] M. L¨uscher and U. Wolff, Nucl. Phys. B , 222 (1990).[64] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes,and R. Sommer, JHEP , 094 (2009), arXiv:0902.1265[hep-lat].[65] J. J. Dudek, R. G. Edwards, N. Mathur, andD. G. Richards, Phys. Rev. D77 , 034501 (2008),arXiv:0707.4162 [hep-lat].[66] J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G.Richards, and C. E. Thomas, Phys. Rev.
D82 , 034508(2010), arXiv:1004.4930 [hep-ph].[67] J. J. Dudek, R. G. Edwards, and C. E. Thomas (HadronSpectrum), Phys. Rev.
D87 , 034505 (2013), [Erratum:Phys. Rev.D90,no.9,099902(2014)], arXiv:1212.0830 [hep- ph].[68] R. A. Brice˜no, J. J. Dudek, R. G. Edwards, C. J. Shultz,C. E. Thomas, and D. J. Wilson, Phys. Rev. Lett. ,242001 (2015), arXiv:1507.06622 [hep-ph].[69] R. A. Brice˜no, J. J. Dudek, R. G. Edwards, C. J. Shultz,C. E. Thomas, and D. J. Wilson, Phys. Rev.
D93 , 114508(2016), arXiv:1604.03530 [hep-ph].[70] R. A. Brice˜no, J. J. Dudek, R. G. Edwards, andD. J. Wilson, Phys. Rev. Lett. , 022002 (2017),arXiv:1607.05900 [hep-ph].[71] J. J. Dudek, R. G. Edwards, and D. J. Wilson(Hadron Spectrum), Phys. Rev.
D93 , 094506 (2016),arXiv:1602.05122 [hep-ph].[72] G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, andD. J. Wilson, JHEP , 011 (2016), arXiv:1607.07093[hep-lat].[73] A. J. Woss and C. E. Thomas, Proceedings, 34th In-ternational Symposium on Lattice Field Theory (Lattice2016): Southampton, UK, July 24-30, 2016 , PoS
LAT-TICE2016 , 134 (2016), arXiv:1612.05437 [hep-lat].[74] R. A. Brice˜no, J. J. Dudek, R. G. Edwards, and D. J.Wilson, Phys. Rev.
D97 , 054513 (2018), arXiv:1708.06667[hep-lat].[75] A. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards,and D. J. Wilson, JHEP , 043 (2018), arXiv:1802.05580[hep-lat].[76] D. J. Wilson, R. A. Brice˜no, J. J. Dudek, R. G. Ed-wards, and C. E. Thomas, Phys. Rev. D92 , 094502(2015), arXiv:1507.02599 [hep-ph].[77] G. K. C. Cheung, C. O’Hara, G. Moir, M. Peardon, S. M.Ryan, C. E. Thomas, and D. Tims (Hadron Spectrum),JHEP , 089 (2016), arXiv:1610.01073 [hep-lat].[78] D. J. Wilson, R. A. Brice˜no, J. J. Dudek, R. G. Edwards,and C. E. Thomas, Phys. Rev. Lett. , 042002 (2019),arXiv:1904.03188 [hep-lat].[79] D. Sadasivan, M. Mai, H. Akdag, and M. D¨oring, Phys.Rev. D , 094018 (2020), arXiv:2002.12431 [nucl-th].[80] R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD), Lattice field theory. Proceedings, 22nd International Sym-posium, Lattice 2004, Batavia, USA, June 21-26, 2004 ,Nucl. Phys. Proc. Suppl. , 832 (2005), arXiv:hep-lat/0409003 [hep-lat].[81] M. A. Clark, R. Babich, K. Barros, R. C. Brower, andC. Rebbi, Comput. Phys. Commun. , 1517 (2010),arXiv:0911.3191 [hep-lat].[82] R. Babich, M. A. Clark, and B. Joo, in
SC 10 (Supercom-puting 2010) New Orleans, Louisiana, November 13-19,2010 (2010) arXiv:1011.0024 [hep-lat].[83] M. A. Clark, B. Jo´o, A. Strelchenko, M. Cheng, A. Gamb-hir, and R. C. Brower, in
SC ’16: Proceedings of the In-ternational Conference for High Performance Computing,Networking, Storage and Analysis (2016) pp. 795–806.[84] A. Jackura, C. Fern´andez-Ram´ırez, V. Mathieu,M. Mikhasenko, J. Nys, A. Pilloni, K. Salda˜na, N. Sherrill,and A. P. Szczepaniak (JPAC), Eur. Phys. J.
C79 , 56(2019), arXiv:1809.10523 [hep-ph].[85] A. J. Woss, C. E. Thomas, J. J. Dudek, R. G. Ed-wards, and D. J. Wilson, Phys. Rev.
D100 , 054506(2019), arXiv:1904.04136 [hep-lat].
FIG. 4. GEVP eigenvalues used to determine the P = [000], Λ = A − , three-pion spectrum on the 20 ×
256 ensemble. Asexplained in the text, the diagonalization is performed with t /a t = 10 and the eigenvalues are rescaled by their expected large- t fall-off. SUPPLEMENTARY MATERIAL1. Spectra
In this section we give details concerning the finite-volume spectra described in the main text. We focus here on tworepresentative examples for the correlators used to extract the three-pion energies. The quality of two-pion correlatorscan be inferred from the earlier work presented in Ref. [54], which includes a partially overlapping data set.As a first example, consider the three-pion spectrum for the P = [000], Λ = A − irrep on the 20 ×
256 ensemble. Inthis case G ij ( t ) = (cid:104)O i ( t ) O † j (0) (cid:105) is an 8 × λ n ( t, t ), of the matrix M ( t, t ) ≡ G ( t ) − / · G ( t ) · G ( t ) − / , (11)entering the generalized eigenvalue problem (GEVP). These are determined for t /a t = 10 and are plotted vs. t/a t fora range of values both before and after the reference time. To display the eigenvalues in a useful manner, we plot thecombination e E n ( t − t ) λ n ( t, t ), where E n has been determined from a two-state fit to λ n ( t, t ) λ n ( t, t ) = (1 − A n ) e − E n ( t − t ) + A n e − E (cid:48) n ( t − t ) . (12)The quality of the fit is indicated by the χ / dof in each panel. The plotted combination behaves as expected for asuccessful GEVP, showing a reasonable plateaux over a range of t/a t .This result also exhibits no evidence for thermal states on this lattice, as expected given the length of the temporalextent, m π T ≈ .
7. We detour slightly, to explain this point in more detail:In general, for multi-pion systems, the leading finite- T effects are given by G ij ( t ) = (cid:104) |O i (0) e − ˆ Ht O † j (0) | (cid:105) + e − m π ( T − t ) (cid:104) π − |O i (0) e − ˆ Ht O † j (0) | π − (cid:105) + · · · , (13)where ˆ H is the Hamiltonian and the ellipsis represents thermal effects falling faster than e − m π T . For concreteness, wehave assumed that O † (0) creates three- π + quantum numbers, so that O † (0) | π − (cid:105) has the quantum numbers of two FIG. 5. GEVP eigenvalues used to determine the P = [000], Λ = A − , three-pion spectrum on the 24 ×
128 ensemble. Thediagonalization is performed with t /a t = 9 and the eigenvalues are rescaled as in Fig. 4. pions with isospin two. A spectral decomposition of Eq. (13) then yields G ij ( t ) = (cid:88) n c ( n ) i c ( n ) ∗ j e − E πππn t + e − m π ( T − t ) (cid:88) n b ( n ) i b ( n ) ∗ j e − E ππn t + · · · , (14)where the sum in the first term (second term) runs over all maximal-isospin three-pion (two-pion) finite-volumestates with specified P . In the case of P = [000], the two- and three-pion ground states take the form N m π + ∆ E N where N = 2 , E N ∼ /L for weakly-interacting systems. Taking the leading ( n = 0) terms of Eq. (14) andsubstituting this scaling then yields G ij ( t ) = c (0) i c (0) ∗ j e − (3 m π +∆ E ) t + b (0) i b (0) ∗ j e − m π ( T + t ) e − ∆ E t + · · · . (15)When this same exercise is performed for a two-pion correlator, again focusing on the case of P = [000], one findsthat the leading thermal contamination is a constant in t in the non-interacting limit. As discussed in Ref. [54] thiscan thus be removed by applying a shift to the correlator G ij ( t ) → G ij ( t ) − G ij ( t + δt ). In the present case, however,the leading contaminations are t -dependent. One option is to reweight and shift, i.e. G ij ( t ) → e − m π t (cid:0) G ij ( t ) e m π t − G ij ( t + δt ) e m π ( t + δt ) (cid:1) . (16)This approach, already used in Ref. [54] for ππ systems with non-zero total momentum, reduces thermal effects at thecost of generally degrading the signal quality. Fortunately, for the 20 ×
256 lattice, this is not required. Comparingthe leading and subleading terms of Eq. (15), and neglecting the interactions, one finds that the relative size of thethree-pion thermal contamination is e − m π ( T − t ) . Thus, assuming that the relevant matrix elements have the sameorder of magnitude, for the range of t considered these effects are ∼ − and are safely negligible, despite the highstatistical precision of the extracted energies. This concludes our general comments on thermal effects.As a second example, in Fig. 5 we consider the same three-pion quantum numbers ( P = [000], Λ = A − ) on the24 ×
128 ensemble. Because the larger spatial volume lowers the value of the n th level, here we include 2 additionaloperators to better absorb the excited states. As with the previous example, the χ / dof and the plotted curves providestrong evidence of a successful GEVP extraction. For this case, e − m π ( T − t ) ∼ − so that finite- T effects potentiallypresent a more significant issue. As a result we have also considered shifting and reweighting, as summarized byEq. (16), in our various fits. However, across all values of P , we find that more stable fits are achieved via the0unmodified correlators, relying on the basis of operators to push the extraction to earlier times and examining theresulting GEVP eigenvalues. This is in contrast to the 24 ×
128 two-pion fits, where the extractions are improved byshifting in certain cases, as described in Ref. [54].
2. K-matrix fits
In this subsection we give additional details concerning the K-matrix fits, summarized in the main text. We presentfour basic types of fits:1. Fitting only the π + π + spectra to various choices of p cot δ ( p ) [Table II].2. Fitting only the π + π + π + spectra to various choices of p cot δ ( p ), with K , iso = 0 [Table III].3. Fitting only the π + π + π + spectra to various choices of K , iso , with p cot δ ( p ) fixed by independent π + π + fits[Table IV].4. Fitting all spectra simultaneously to various choices of p cot δ ( p ) and K , iso [Table V].Here δ ( p ) is the S -wave, π + π + scattering phase shift, related to the scattering amplitude via M ( E (cid:63) ) = 16 πE (cid:63) p cot δ ( p ) − ip , (17)where p = E (cid:63) / − m π . One standard parametrization of the scattering amplitude follows from the effective rangeexpansion p cot δ ( p ) = − a + 12 r p + O ( p ) , (18)and below we present fits to the leading term as well as to the leading two terms.In the case of π + π + scattering, chiral perturbation theory predicts the Adler zero, which leads to a pole in p cot δ ( p ),limiting the range of convergence for the effective range expansion. This motivates the alternative form p cot δ ( p ) = A ( c, p ) (cid:104) − a + c (cid:48) p + O ( p ) (cid:105) , (19)where we have introduced A ( c, p ) ≡ c m π (cid:112) p + m π p + c m π . (20)The function A ( c, p ) is chosen to match the analytic structure predicted by leading-order chiral perturbation theory(both the energy numerator and the pole in the denominator) and is normalized so that A ( c,
0) = 1. The leading-orderprediction for the pole position corresponds to c = 1 and in the following we present fits both with c fixed and allowedto vary.As we explain in detail in Sec. 3 below, the three-particle scattering amplitude, M , is determined from thetwo-particle scattering amplitude together with a local three-particle K matrix, first introduced in Ref. [5] anddenoted by K df , . As already described in the main text, we work here in the isotropic approximation, for which thisthree-particle K matrix reduces to a simple function of the total center-of-momentum frame energy, denoted K , iso ( E (cid:63) ).This quantity admits an expansion similar to the effective range expansion K , iso ( E (cid:63) ) = c /m π + c ∆ /m π + O (∆ ) , (21)where ∆ ≡ E (cid:63) − m π . The fits presented below take either the first or else the first two terms in this expansion.In each case, the fit is performed by minimizing the χ ( { η i } ), where χ ( { η i } ) ≡ [ E d − E ( { η i } )] · C − · [ E d − E ( { η i } )] T . (22)Here E d is a vector built from two and three-particle energies extracted from the lattice calculation, C is the covariancematrix, and E ( { η i } ) is a vector of solutions to the two- and three-particle quantization conditions. The solved energies1 Fit E (cid:63) , cut p cot δ ( p ) fit result χ / dofA . m π − /a m π a = 0 . ± .
007 89 . / (32 −
1) = 2 . . m π − /a m π a = 0 . ± .
010 26 . / (21 −
1) = 1 . . m π − /a + r p / m π a = 0 . ± . m π r = − . ± . (cid:20) . − . . (cid:21) . / (32 −
2) = 2 . . m π − /a + r p / m π a = 0 . ± . m π r = 0 . ± . (cid:20) . − . . (cid:21) . / (21 −
2) = 1 . . m π − A ( c , p ) /a m π a = 0 . ± . m π c = 11 . ± . (cid:20) . − . . (cid:21) . / (32 −
2) = 3 . . m π − A ( c , p ) /a m π a = 0 . ± . m π c = 3 . ± . (cid:20) . − . . (cid:21) . / (21 −
2) = 1 . . m π A (1 , p )( − /a + c p ) m π a = 0 . ± . m π c = − . ± . (cid:20) . . . (cid:21) . / (32 −
2) = 2 . . m π A (1 , p )( − /a + c p ) m π a = 0 . ± . m π c = − . ± . (cid:20) . . . (cid:21) . / (21 −
2) = 1 . [ r = 0 .
01 (9 / . m π − /a m π a = 0 . ± .
008 40 . / (32 −
1) = 1 . [ r = 0 .
01 (9 / . m π − /a + r p / m π a = 0 . ± . m π r = − . ± . (cid:20) . − . . (cid:21) . / (32 −
2) = 1 . π + π + finite-volume energies, for various choices of p cot δ ( p ). All values of total momentum P (from [000] to [002]) and both volumes (20 and 24 ) are used in each fit. The entries below the lower double horizontal line aredetermined using a regularized covariance matrix as explained in the text. The function A ( c , p ) = c m π (cid:112) p + m π / (2 p + c m π )encodes the effect of the Adler zero, with c = 1 corresponding to the pole position from leading-order chiral perturbation theory.The columns are understood as follows: “Fit” gives a label to the fit (and defines the regularized covariance matrix for the finaltwo fits); “ E (cid:63) , cut ” gives the two-particle center-of-momentum frame energy cutoff (i.e. only points with central values below thisthreshold enter the fit); “ p cot δ ( p )” indicates the fit function; “fit result” displays the extracted parameters and their correlation;“ χ / dof” gives the value of χ ( { η i } ) (evaluated at the best fit parameters) divided by the number of degrees of freedom. depend on { η i } , which stands for all two- and three-particle K-matrix parameters, over which the minimization isperformed.In certain cases the low-lying eigenvalues of C cannot be estimated reliably and, if underestimated, can lead toartificially enhanced eigenvalues in C − , and therefore inflated values for χ ( { η i } ). To study this problem we have alsoconsidered an alternative method in which the low lying eigenvalues of C are adjusted. To do so, one first diagonalizes C C = R T · D · R , (23)where R is an orthogonal matrix of eigenvectors and D a diagonal matrix of eigenvalues. We label the eigenvalues(ordered from smallest to largest) by λ , · · · , λ N and note that these are positive. We assume also that the rows andcolumns of R are organized such that D = diag[ λ , λ , · · · , λ N ] . (24) χ ( { η i } ) may be poorly estimated if there is a large hierarchy between the smallest and largest eigenvalues, λ and λ N , respectively. This motivates the definition C r = R T · D r · R , (25)where D r is a diagonal matrix defined as D r = diag (cid:104) max[ λ , rλ N ] , max[ λ , rλ N ] , · · · , max[ λ N , rλ N ] (cid:105) . (26)Note, if r = 0, then D = D r = ⇒ C = C r . As this parameter is increased, the lowest eigenvalues are adjusted to somefixed fraction of the largest value. This approach defines a new test statistic and, in principle, one can sample its2corresponding distribution to define p -values and assess the quality of fits. This goes beyond the scope of this workand we only perform the modified fits as a cross check to show that the extracted fit parameters are robust underthese regularizations of the covariance matrix. Such fits are reported in the tables of this section with labels of theform [ r = 0 .
01 ( n/m )] where r indicates the adjustment parameter and ( n/m ) indicates the number of eigenvaluesthat have been changed versus the total number. We have additionally re-done many of the fits summarized TablesII-V with the pion mass, m π , and the anisotropy, ξ , varied by one standard deviation. We find in all cases that theeffect of this shift is well below the statistical uncertainties on the extracted fit parameters. Fit E (cid:63) , cut p cot δ ( p ) fit result χ / dofA K . m π − /a m π a = 0 . ± .
011 31 . / (16 −
1) = 2 . K [ r = 0 .
01 (5 / . m π − /a m π a = 0 . ± .
014 24 . / (16 −
1) = 1 . π + π + π + finite-volume energies, for p cot δ ( p ) = − /a with K , iso = 0 fixed. All values of totalmomentum P (from [000] to [111]) and both volumes (20 and 24 ) are used in each fit. Columns as in Table II, with “ E (cid:63) , cut ”indicating the three-particle center-of-momentum frame energy cutoff.Fit p cot δ ( p ) E (cid:63) , cut K , iso fit result χ / dofA K − /a (B ) 4 . m π c /m π c = − ±
874 31 . / (16 −
1) = 2 . K − /a (B ) 4 . m π c /m π + c ∆ /m π c = 5039 ± c = − ± (cid:20) . − . . (cid:21) . / (16 −
2) = 1 . K [ r = 0 .
01 (5 / − /a (B ) 4 . m π c /m π c = − ±
892 24 . / (16 −
1) = 1 . π + π + π + finite-volume energies, for various choices of K , iso (with p cot δ ( p ) given by fit B of Table II). Columns as in Tables II and III, with “ K , iso ” indicating the fit function used and ∆ ≡ E (cid:63) − m π encoding alinear-dependence in the squared center-of-momentum frame energy.Fit E (cid:63) , cut E (cid:63) , cut p cot δ ( p ) K , iso fit result χ / dofA . m π . m π − /a m π a = 0 . ± .
007 64 . / (37 −
1) = 1 . . m π . m π − /a c /m π m π a = 0 . ± . c = − ± (cid:20) . . . (cid:21) . / (37 −
2) = 1 . [ r = 0 .
005 (11 / . m π . m π − /a m π a = 0 . ± .
008 50 . / (37 −
1) = 1 . [ r = 0 .
005 (11 / . m π . m π − /a c /m π m π a = 0 . ± . c = − ± (cid:20) . . . (cid:21) . / (37 −
2) = 1 .
3. Details of the three-particle integral equations
In this subsection we prove Eq. (5) and give details on its numerical implementation. We begin by reviewing theintegral equations presented in Ref. [6]. As the results of the fits summarized in the main text (and detailed in theprevious subsection) are consistent with K , iso = 0, we focus here on the case of a weak three-body interaction, keepingonly the linear contribution in this term. We begin with the unsymmetrized three-body scattering amplitude M ( u,u )3 ( p , k ) = D ( u,u ) ( p , k ) + E ( u ) ( p ) K , iso E ( u ) ( k ) + O ( K , iso ) , (27)3where the superscripts indicate the lack of exchange symmetry, and k and p specify the momenta of the so-calledspectator particles in the initial and final state, respectively. In general, the factors appearing in Eq. (27) carry angularmomentum indices. However, as discussed in the main text, the π + π + system at low energies is dominated by the S -wave component. Thus, we restrict attention here to π + π + with zero angular momentum, such that M ( u,u )3 ( p , k )and D ( u,u ) ( p , k ) are simple functions, with no implicit indices.In this limiting case, D ( u,u ) ( p , k ) satisfies the implicit equation D ( u,u ) ( p , k ) = −M ( E (cid:63) ,p ) G ∞ ( p , k ) M ( E (cid:63) ,k ) − M ( E (cid:63) ,p ) (cid:90) d k (cid:48) (2 π ) ω k (cid:48) G ∞ ( p , k (cid:48) ) D ( u,u ) ( k (cid:48) , k ) , (28)where E (cid:63) ,k ≡ ( E − ω k ) − k is the center-of-momentum energy for the non-spectator pair (with k = | k | defined in thethree-particle zero-momentum frame). In words, the unsymmetrized amplitude D ( u,u ) can be evaluated by solving anintegral equation depending only on the two-particle scattering amplitude M ( E (cid:63) ,p ) and the exchange propagator G ∞ ( p , k ) ≡ H ( p, k ) b pk − m + i(cid:15) . (29)Here we have suppressed the π subscript on the mass and have introduced b pk ≡ ( E (cid:63) − ω p − ω k ) − ( p + k ) as well as H ( p, k ), a cut-off function built into the relation between finite-volume energies and K , iso , as well as that between K , iso and the physical scattering amplitude. The definition used here is H ( p, k ) ≡ J (cid:0) E (cid:63) ,k / [2 m ] (cid:1) J (cid:0) E (cid:63) ,p / [2 m ] (cid:1) , J ( x ) ≡ , x ≤ (cid:16) − x exp (cid:104) − − x (cid:105)(cid:17) , < x ≤ , < x . (30)(See also Refs. [5, 6] for more discussion on this technical aspect.) Finally E ( u ) ( p ) is a function closely related to D ( u,u ) ( p , k ) and defined, in a specific limiting case, in Eq. (38) below. In this work it is only relevant to demonstratethat K , iso contributes negligibly to both the central value and uncertainty of the three-hadron scattering amplitude,as we describe in the following section.The next step is to project the remaining directional freedom within D ( u,u ) onto vanishing three-particle angularmomentum, i.e. to J = 0. Defining D ( u,u ) s ( p, k ) ≡ (cid:90) d Ω ˆ k π d Ω ˆ p π D ( u,u ) ( p , k ) , (31)one can show that D s satisfies a one-dimensional integral equation of the form D ( u,u ) s ( p, k ) = −M ( E (cid:63) ,p ) G s ( p, k, (cid:15) ) M ( E (cid:63) ,k ) − M ( E (cid:63) ,p ) (cid:90) k max k (cid:48) dk (cid:48) (2 π ) ω k (cid:48) G s ( p, k (cid:48) , (cid:15) ) D ( u,u ) s ( k (cid:48) , k ) , (32)where G s ( p, k, (cid:15) ) ≡ (cid:90) d Ω ˆ p π d Ω ˆ k π G ∞ ( p , k ) = − H ( p, k )4 pk log (cid:20) pk − ( E − ω k − ω p ) + p + k + m − i(cid:15) − pk − ( E − ω k − ω p ) + p + k + m − i(cid:15) (cid:21) . (33)In Fig. 6 we plot G s for a range of kinematic values. Setting K , iso = 0 and combining Eqs. (27), (32) and (33), wearrive at Eq. (5) of the main text. Note that, in Eq. (32), we have included an explicit cutoff at k max = ( s + m ) / (2 √ s ).This is done without any additional approximation as H ( p, k ) has vanishing support for k > k max .We comment here that, in general, projecting M ( u,u )3 to definite partial waves is not equivalent to doing the same forthe symmetrized three-hadron scattering amplitude. However, as discussed around Eqs. (15) and (16) of Ref. [84] anydefinite-angular-momentum component of the the fully symmetric M can be assembled from known combinationsof the various angular-momentum components of M ( u,u )3 . Thus the latter form a basis for constructing the physicalthree-hadron amplitudes.To solve Eq. (32) numerically we replace the integral (cid:82) k max dk (cid:48) with a discrete sum (cid:80) k (cid:48) ∆ k containing N terms.Then a discretized version of the equation can be written in a matrix form D ( N, (cid:15) ) = − M · G ( (cid:15) ) · M − M · G ( (cid:15) ) · P · D ( N, (cid:15) ) , (34)4 . . . . . E ? ,k /m π − . − . . . . R e [ p k G s ] √ s = 3 . m π √ s = 3 . m π √ s = 3 . m π √ s = 3 . m π . . . . . E ? ,k /m π − . − . . . . I m [ p k G s ] √ s = 3 . m π √ s = 3 . m π √ s = 3 . m π √ s = 3 . m π FIG. 6. Plots of the kernel function G s ( p, k, (cid:15) ) [defined Eq. (33)] vs E (cid:63) ,k , here for four different energies as shown in the legendand with p = 0 . m π and (cid:15) → + . For Im[ p k G s ( p, k, (cid:15) )], in the right panel, we have vertically off-set the curves for improvedreadability. where we have introduced the following N × N matrix representations G pk ( (cid:15) ) = G s ( p, k, (cid:15) ) , M k (cid:48) k = δ k (cid:48) k M ( E (cid:63) ,k ) , P k (cid:48) k = δ k (cid:48) k k ∆ k (2 π ) ω k , (35)as well as D ( N, (cid:15) ), which becomes our target quantity in the ordered double limit D ( u,u ) s ( p, k ) = lim (cid:15) → lim N →∞ D pk ( N, (cid:15) ) . (36)Here, in a slight abuse of notation, the indices p k represent the choices that are closest to physical momenta p and k for a given N value. Eq. (34) can then be solved through a matrix inverse to yield D ( N, (cid:15) ) = − (cid:2) I + M · G ( (cid:15) ) · P (cid:3) − · M · G ( (cid:15) ) · M . (37)Finally we return to the endcap factors appearing on either side of K , iso in Eq. (27). These are defined in Eqs. (105)and (106) of Ref. [5] and in the overall S -wave approximation they take the form E ( u ) s ( p ) = 13 − M ( E (cid:63) ,p ) ρ ( E (cid:63) ,p ) − (cid:90) dp (cid:48) p (cid:48) (2 π ) ω p (cid:48) D ( u,u ) s ( p, p (cid:48) ) ρ ( E (cid:63) ,p (cid:48) ) , (38)where ρ ( E (cid:63) ,p ) ≡ − i J (cid:0) E (cid:63) ,k / [2 m ] (cid:1) πE (cid:63) ,k (cid:113) E (cid:63) ,p / − m , (39)and it is understood that the the i √− x branch is taken for x < K , iso is consistent with zero in all fits we have considered, the main purpose in keeping track of these quantitiesis to estimate the propagation of uncertainties into M . In particular we note∆ M = (cid:88) η,η (cid:48) [ ∂ η M ] C ηη (cid:48) [ ∂ η (cid:48) M ] , (40)where C ηη (cid:48) represents the fit-parameter covariance matrix and the sums run over all inputs to the scattering amplitude,in particular the scattering length a as well as K , iso . In the next section we present numerical solutions for M basedon the parameters extracted from the finite-volume energies.5
4. Solutions for M . . . . . . E ? ,k /m π − R e [ p k M ( u , u ) s ] . . . . . . E ? ,k /m π I m [ p k M ( u , u ) s ] FIG. 7. Real (left) and imaginary (right) parts of p k M ( u,u ) s , determined by solving Eq. (37) using m π a = 0 .
296 and K , iso = 0for the central values. As we explain in the text, the uncertainties here follow from propagating the uncertainties on K , iso and m π a , taken from Fit B in Table V. Figure 7 shows the result of solving Eq. (37) and varying N and (cid:15) in order to ensure that the ordered double limithas been saturated. Symmetrizing this function over the three incoming and three outgoing momenta yields the fullthree-hadron scattering amplitude, plotted in Fig. 3 of the main text.To assign an uncertainty estimate to M ( u,u ) s ( p, k ) we have applied a straightforward adaptation of Eq. (40).Specializing to the case of only m π a and K , iso as input parameters, and neglecting the correlations between them,this becomes (cid:0) ∆ M ( u,u ) s ( p, k ) (cid:1) = (cid:0) ∆[ m π a ] (cid:1) (cid:18) ∂ M ( u,u ) s ( p, k ) ∂ ( m π a ) (cid:19) + (cid:0) ∆ K , iso / (cid:1) + O (∆ K , iso ρ M ) , (41)where ∆( m π a ) and ∆ K , iso denote the uncertainties on the input parameters and ∆ M ( u,u ) s ( p, k ) is the resultingamplitude uncertainty plotted in the figure. In practice one finds that ∆ K , iso contributes negligibly to the overalluncertainty, simply because the series of M and G s insertions dominates the value of the amplitude for weakly-interacting systems. For this reason we have taken the leading part of E ( u ) s , resulting in the ∆ K , iso / m π a . Since M ( u,u ) s ( p, k ) is dominated by thecontribution proportional to M , in this case the overall uncertainty is well approximated by∆ M ( u,u ) s ( p, k ) M ( u,u ) s ( p, k ) = 2 ∆[ m π a ] m π a . (42)To produce Fig. 7 we have used ∆[ m π a ] = 0 .
016 and ∆ K , iso = 770. These are taken from B in Table V with theuncertainty on m π a doubled to account for systematic variations between the various fits that have been performed.
5. Operator construction and lists
Following Ref. [54], to determine the ππ I = 2 finite-volume energies we compute correlation functions featuringoperators constructed to resemble a ππ structure. Schematically, for an operator in lattice irrep Λ and row µ withoverall momentum k ,( ππ ) [ k , k ] † Λ µ ( k ) = (cid:88) k , k k + k = k C ( k , Λ , µ ; k , Λ ; k , Λ ) π † Λ ( k ) π † Λ ( k ) , (43)where the sum is over all momenta related to k and k by allowed lattice rotations, C is an appropriate generalizedClebsch-Gordan coefficient for Λ ⊗ Λ → Λ , and flavor indices and the projection onto I = 2 are not writtenexplicitly. Here π † Λ i ( k i ) is the optimal linear combination of operators to interpolate a π with momentum k i in irrep These are all one dimensional and so we omit the irrep row index. P Λ L/a s = 20 L/a s = 24[000] A +1 π [000] π [000] , π [100] π [100] , π [110] π [110] , π [111] π [111] , π [000] π [000] , π [100] π [100] , π [110] π [110] , π [111] π [111] , π [200] π [200] π [200] π [200] [100] A π [000] π [100] , π [100] π [110] , π [110] π [111] , π [100] π [200] π [000] π [100] , π [100] π [110] , π [110] π [111] , π [100] π [200] , π [110] π [210] , π [200] π [210] , π [111] π [211] [110] A π [000] π [110] , π [100] π [100] , π [100] π [111] , π [110] π [110] π [000] π [110] , π [100] π [100] , π [100] π [111] , π [110] π [110] , π [110] π [200] , π [100] π [210] , π [111] π [210] , π [110] π [211] [111] A π [000] π [111] , π [100] π [110] π [000] π [111] , π [100] π [110] , π [111] π [200] , π [110] π [210] , π [100] π [211] [200] A π [100] π [100] , π [000] π [200] , π [110] π [110] , π [111] π [111] π [100] π [100] , π [000] π [200] , π [110] π [110] , π [100] π [210] , π [111] π [111] , π [110] π [211] , π [210] π [210] TABLE VI. The ππ I = 2 operators, π k π k , used to compute the finite-volume energy levels shown in Fig. 1 (upper plots) inirrep Λ with overall momentum P . These are constructed from optimized π operators with momentum types k and k ; differentmomentum directions are summed over as in Eq. (43). Momenta are displayed using the shorthand notation [ ijk ] = πL ( i, j, k ). Λ i using a basis of fermion-bilinear operators featuring various Dirac γ matrices and gauge-covariant derivatives –see Ref. [54] for details. In this work the basis of fermion-bilinear operators used for a π operator has up to threederivatives for π at rest and up to one derivative for π at non-zero momentum, except we use up to two derivatives for1 ≤ | k i | ≤ volume. The operators used to compute the ππ spectra shown in Fig. 1 are listed in Table VI.In a similar way, operators used to compute πππ I = 3 energies resemble a πππ structure and are formed bycombining a ππ I = 2 operator with a π operator, as detailed in Ref. [85]. Schematically, for an operator in latticeirrep Λ and row µ with overall momentum P ,( πππ ) [ k [ k , k ] , k ] † Λ µ ( P ) = (cid:88) k , k k + k = P C ( P , Λ , µ ; k , Λ , µ ; k , Λ ) ( ππ ) [ k , k ] † Λ µ ( k ) π † Λ ( k ) , (44)where the sum is over all momenta related to k and k by allowed lattice rotations and, again, flavor indices and theprojection onto I = 3 are not written explicitly.From Bose symmetry, a πππ system must be symmetric under the interchange of any pair of pions. The pions haveno intrinsic spin and we are considering I = 3 which means that the flavor structure is symmetric under interchangeof any pair; therefore, the spatial structure must also be symmetric under the interchange of any pair of pions. Theoperator construction in Eq. (44) gives operators with the correct symmetry properties because the three pions areidentical, but it does not make this symmetry manifest and two different sets of ( | k | , | k | , | k | , | k | , Λ ) maylead to equivalent operators, or a number of different sets may give linearly-dependent operators. To ensure wehave an appropriate set of independent operators, for each P and Λ, we construct an operator for every possible set( | k | , | k | , | k | , | k | , Λ ) with | k | + | k | + | k | less than some cutoff. We write Eq. (44) schematically as,( πππ ) Λ µ ( P ) = (cid:88) (cid:101) C ( k , k , k ) π † ( k ) π † ( k ) π † ( k ) = (cid:101) C · V [ πππ ] , (45)where (cid:101) C with no labels represents a row vector of coefficients and V [ πππ ] a column of operators, such that the dot-product reproduces the sum. We then introduce the matrix R ( ijk ) which acts on V [ πππ ] by mapping a given entry π † ( p ) π † ( p ) π † ( p ) into π † ( p i ) π † ( p j ) π † ( p k ). This allows us to define the symmetrized vector of Clebsch-Gordancoefficients (cid:101) C · R (123) + (cid:101) C · R (231) + (cid:101) C · R (312) + (cid:101) C · R (132) + (cid:101) C · R (213) + (cid:101) C · R (321) , (46)where R (123) is just the identity matrix. The final step is to check that the resulting vector is non-zero and linearlyindependent from the analogous expressions for the already considered operators. The resulting sets of independentoperators used in this work are listed in Tables VII, VIII, IX, X and XI. Strictly we mean the types of momenta (i.e. the equivalence classof momenta related by rotations in the octahedral group or littlegroup) rather than the magnitudes, but there is no distinction for the momenta we are considering here. This is not always a unique choice and any independent set ofoperators could be used to achieve the same results. d , d , d , eachgiven by d i = L k i / (2 π ). The vectors are assigned a color (orange and green for the first two pions, and blue for thethird) and the absence of any given color corresponds to a vector of magnitude zero. As summarized by Eq. (44), ouroperator construction is based on combining two-pion operators in a definite irrep with the third pion. An alternativebasis is given by summing a given momentum assignment, represented by a given set d d d , over all rotations in theoctahedral group (in the case of P = [000]) or else a little group thereof (for non-zero total momentum) weighted bythe appropriate Clebsch-Gordan coefficients. The operators reached via this alternative construction are equal to alinear combination of those given by Eq. (44).For example, on the third line of Table VIII, two distinct momentum assignments arise from combining the ππ [000] A +1 with the third π [100] operator. In this case the ππ [000] A +1 is built from individual pions with a unit of back-to-backmomentum. When one sums over the coefficients projecting onto A +1 , contributions arise with the back-to-backaxis both aligned and perpendicular to the total momentum direction, P = [001]. Thus, momentum assignmentscorresponding to both diagrams shown in line 3 of Table VIII contribute to the operator on that line. By contrast, online 4 of Table VIII only a single momentum configuration contributes, as indicated. This implies that the two operatorsare independent, since they are built from independent linear combinations of the two momentum configurations.Operators 8 through 11 of Table VIII give a more complicated example. The first two (8 and 9) correspond to twolinear combinations of two configuration types, labeled with subscripts 1 and 3, and the next two (10 and 11) areequal to linear combinations of the operators labeled 1, 2, and 4.The diagrams in Tables VII, VIII, IX, X and XI provide a cross check on the linear-independence of the operators.Each row corresponds to a linear combination of the displayed momentum configurations and the number of linearlyindependent operators is equal to the number of distinct diagrams.8 d d d π k π k ππ k Λ ( I = 2) π k momentum configurations1 π [000] π [000] ππ [000] A +1 π [000] π [100] π [000] ππ [100] A π [100] π [110] π [000] ππ [110] A π [110] π [100] π [100] ππ [110] A π [110] π [111] π [111] ππ [000] A +1 π [000] π [200] π [100] ππ [100] A π [100] π [111] π [100] ππ [110] A π [110] π [110] π [110] ππ [110] A π [110] π [200] π [000] ππ [200] A π [200] π [210] π [100] ππ [110] A π [110] TABLE VII. The πππ I = 3 operators used to compute the finite-volume energy levels shown in Fig. 1 (lower left plot) inirrep A − with overall momentum P = [000], labeled by d d d where d i = k i ( L/ π ). Different momentum directions forthe momentum types k , k , k and k are summed over as in Eq. (44). Operators 1 to 8 are used on the 20 volume andoperators 1 to 10 are used on the 24 volume. The momentum configuration diagrams in the rightmost column are explainedin the text. Operators separated by a single horizontal line correspond to states that are degenerate in the non-relativistic,non-interacting theory. d d d π k π k ππ k Λ ( I = 2) π k momentum configurations1 π [000] π [000] ππ [000] A +1 π [100] π [100] π [000] ππ [100] A π [110] π [100] π [100] ππ [000] A +1 π [100] 1 2 π [100] π [100] ππ [110] A π [100] 1 π [100] π [000] ππ [100] A π [200] π [110] π [000] ππ [110] A π [111] π [100] π [100] ππ [110] A π [111] π [110] π [100] ππ [111] A π [110] 1 3 π [110] π [100] ππ [111] E π [110] π [110] π [100] ππ [100] A π [110] 1 2 4 π [110] π [100] ππ [100] B π [110] π [110] π [000] ππ [110] A π [210] π [210] π [100] ππ [110] A π [100] π [210] π [100] ππ [110] B π [100] π [200] π [100] ππ [210] A π [110] π [200] π [100] ππ [100] A π [110] TABLE VIII. As Table VII but for the A irrep with overall momentum P = [001]. Operators 1 to 8 are used on the 20 volumeand operators 1 to 16 are used on the 24 volume. Operators separated by a single horizontal line correspond to states thatare degenerate in the non-relativistic, non-interacting theory. Operators with no horizontal line correspond to states that aredegenerate in the relativistic, non-interacting theory but are split by the interactions. d d d π k π k ππ k Λ ( I = 2) π k momentum configurations1 π [000] π [000] ππ [000] A +1 π [110] π [100] π [000] ππ [100] A π [100] π [100] π [000] ππ [100] A π [111] π [110] π [000] ππ [110] A π [110] π [100] π [100] ππ [000] A +1 π [110] 1 2 π [100] π [100] ππ [110] A π [110] π [100] π [100] ππ [200] A π [110] π [100] π [100] ππ [000] E + π [110] 1 2 TABLE IX. As Table VII but for the A irrep with overall momentum P = [011] (continued in Table X). Operators 1 to 11 areused on the 20 volume and operators 1 to 20 are used on the 24 volume. d d d π k π k ππ k Λ ( I = 2) π k momentum configurations9 π [100] π [000] ππ [100] A π [210] π [110] π [000] ππ [110] A π [200] π [200] π [100] ππ [210] A π [100] π [111] π [100] ππ [211] A π [110] π [111] π [100] ππ [110] A π [110] 1 2 π [111] π [100] ππ [110] B π [110] 1 2 π [110] π [110] ππ [000] A +1 π [110] 1 2 3 π [110] π [110] ππ [211] A π [110] 1 π [110] π [110] ππ [110] A π [110] 1 4 π [110] π [110] ππ [200] A π [110] 4 3 π [110] π [000] ππ [110] A π [211] π [111] π [000] ππ [111] A π [210] TABLE X. Continuation of Table IX. d d d π k π k ππ k Λ ( I = 2) π k momentum configurations1 π [000] π [000] ππ [000] A +1 π [111] π [110] π [000] ππ [110] A π [100] π [100] π [100] ππ [110] A π [100] π [100] π [100] ππ [000] A +1 π [111] π [100] π [100] ππ [200] A π [111] π [110] π [100] ππ [100] A π [110] 1 2 π [110] π [110] ππ [110] A π [100] 2 TABLE XI. As Table VII but for the A irrep with overall momentum P = [111]. Operators 1 to 7 are used on both the 20 and 243