Two- and three-pion finite-volume spectra at maximal isospin from lattice QCD

We present the three-pion spectrum with maximum isospin in a finite volume determined from lattice QCD, including, for the first time, excited states across various irreducible representations at zero and nonzero total momentum, in addition to the ground states in these channels. The required correlation functions, from which the spectrum is extracted, are computed using a newly implemented algorithm which reduces the number of operations, and hence speeds up the computation by more than an order of magnitude. The results for the $I = 3$ three-pion and the $I = 2$ two-pion spectrum are publicly available, including all correlations, and can be used to test the available three-particle finite-volume approaches to extracting three-pion interactions.


INTRODUCTION
Lattice QCD calculations of scattering amplitudes have matured significantly over the last decade owing to marked increases in available computational capacity and improved algorithms. A widely used approach for constraining scattering observables from simulations relies on precise measurements of the interacting energy levels of QCD in a finite volume, which encode hadron interactions via the shifts from their noninteracting values [1][2][3][4][5] (see [6] for a survey of extensions of the formalism and numerical results).
So far, practical calculations in lattice QCD have been mostly confined to the two-hadron sector. Though a large abundance of lattice data is currently available for meson-meson scattering (e.g. ππ scattering in all three isospin channels [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], see also [25,26] for results using a potential-based approach), these calculations are formally restricted to energies below thresholds involving three or more hadrons due to the use of a formalism for relating finite-volume spectra to scattering amplitudes that is limited to two-hadron scattering. This limitation has precluded a proper lattice QCD study of systems involving three or more stable hadrons at light pion masses, e.g. the Roper resonance which decays to both two-and three-particle channels, the ω(782) decaying to three pions, many of the X, Y and Z resonances, and three-nucleon interactions relevant for nuclear physics.
However, significant progress has been made recently in developing the necessary formalism to interpret the three-particle finite-volume spectrum (for a review see [27]), both by extending the two-particle derivation to include three-hadron states [28][29][30][31], as well as through alternative approaches [32][33][34][35][36] [72]. Thus, although the three-particle formalism is quite mature-including numerical explorations of the corresponding quantization conditions [36,37][73]-data for three-particle finitevolume QCD spectra is lacking since previous lattice QCD calculations have been restricted to the extraction of multi-meson ground states at rest [38][39][40].
We fill this gap by providing the two-pion and threepion spectra with maximum isospin in various irreducible representations at zero and nonzero total momentum, including not only the ground states but the excited states in the elastic region as well, i.e. for center-of-mass energies E cm /m π below 4 and 5 for isospin I = 2 and I = 3 respectively. This data, which is made public, including all correlations, will allow for an investigation of the various three-particle interaction parameters as well as the effect of higher partial waves, for which the quantization condition has been worked out recently [41].
A technical challenge concerns the growing number of Wick contractions required to compute correlation functions of suitable interpolating operators-from which the spectrum is extracted-as the number of valence quark fields increases. The continued need for improved algorithms to perform these contractions was pointed out recently [42] and indeed was a limiting factor in a recent study of meson-baryon scattering in the ∆ channel [43]. While Refs. [44][45][46][47][48][49] investigated efficient contraction algorithms at the quark level, we employ the stochastic variant [50] of distillation [51] to treat quark propagation. In this framework, it is useful to view the correlation function construction in terms of contractions of tensors associated with the involved hadrons. Then, to reduce the operation count required to evaluate all contractions, we use a method which is well-known in quantum chemistry [52][53][54] and has attracted renewed interest recently in the context of tensor networks [55]. The proposed optimization achieves an operation-count reduction by more than an order of magnitude, and its implementation is made publicly available.
This letter is organized as follows: We first describe the interpolating operators employed in this work and the method used to speed up the construction of their correlation functions. This is followed by a presentation of the analysis and results.

LATTICE QCD METHODS
Interpolating operators: In lattice QCD the spectrum in finite volume is customarily extracted employing interpolating operators which transform irreducibly under the symmetries of a cubic spatial lattice, i.e. the octahedral group O h for zero total momentum P = 0 and the corresponding little groups for P = 0 [11,56]. Correlation functions of such interpolators access only the subblock of the finite-volume Hamiltonian corresponding to the same irreducible representation (irrep), thus greatly facilitating the determination of the spectrum and the subsequent scattering-amplitude analysis.
We employ the simplest single-pion operator destroying a three-momentum p given by This operator transforms in the A − 1u and A − 2 irrep for zero and nonzero momentum respectively, where the superscript specifies the G-parity.
Two-pion interpolators which transform according to the irrep Λ of the little group of total momentum P are obtained by forming appropriate linear combinations of two single-pion interpolators with momenta p 1 +p 2 = P, The relevant Clebsch-Gordan coefficients c (P,Λ) p 1 ,p 2 were worked out in [56] (see also [11]) and used previously to study ππ scattering [15,19].
Three-pion interpolators are obtained by iterating this process, i.e. by first coupling two of the pions into an intermediate irrep, then using the Clebsch-Gordan coefficients again to obtain operators transforming according to one of the total irreps of interest, πππ (P,Λ) (t) = c (P,Λ) p 1 ,p 2 ,p 3 π p 1 (t)π p 2 (t)π p 3 (t).
The interpolators we use in this work are listed in Tables III to VI in the appendix.
Correlation function construction: Quark propagation is treated using the stochastic LapH method [50] by first obtaining smeared solutions of the Dirac equation, with summed color index a and spin index α, and two open noise and dilution indices. In terms of these meson functions, the single-pion correlation function on a single gauge configuration is obtained by the average over noise combinations with proper normalization given by the number of noise combinations used to perform the average. For a given momentum, pair of source and sink time t s and t f , and noise combination, Eq. (7) is a tensor contraction over dilution indices of two rank-2 tensors with index range N dil . Two-and three-pion correlation functions with maximal isospin can be computed using the same building blocks [50] and involve tensor contractions governed by all possible Wick contractions of four and six rank-2 tensors respectively. Linear combinations of products of these meson functions with various momenta, but the same total momentum, are required to compute correlation functions in a definite irrep according to the Clebsch-Gordan coefficients discussed in the previous subsection. In the following we refer to an expression in terms of tensor contractions for a given pair of source and sink times, and each meson function with a single noise combination and momentum, as a diagram. The different topologies of diagrams for three-pion correlation functions in the sector of maximal isospin are shown in Figure 1.
The number of Wick contractions grows factorially as more pions are included [38]. However, across different diagrams there is a lot of redundancy which can be exploited systematically to reduce the number of arithmetic operations required to evaluate all the necessary diagrams. In particular, all diagrams required for the I: Number of elemental contractions with given complexity in N dil required to compute all correlation functions used in this work for a fixed time separation on a single gauge configuration with and without common subexpression elimination (CSE) and diagram consolidation (DC). Diagram consolidation does not reduce the number of required contractions when CSE is used, but speeds up the optimization process. computation of I = 2 two-pion correlation functions appear as subdiagrams of I = 3 three-pion correlation functions. The method we use to reduce the operation count required for the evaluation of tensor contractions was proposed in the context of quantum chemistry [52][53][54] and consists of two parts.
In the first part, for a given diagram, using pairwise contraction of two tensors over all joint indices as the elemental computational kernel, the list of possible locally optimal next computational steps is determined by requiring that the proposed step remove as many contractions as possible at the lowest cost. This step is irrelevant for multi-meson correlation functions as all contractions have complexity of either N 2 dil if both indices are contracted over or N 3 dil if only one index is contracted with two spectator indices; this step is however useful in meson-baryon and baryon-baryon correlation function construction where the arithmetic operation count is reduced by powers of N dil through judicious choice of the order of elemental contractions [74].
In the second part, each of the locally optimal proposed steps is ranked against all diagrams to identify the step which appears most frequently as a subexpression globally. The best-ranking contraction is performed and the resulting intermediary object substituted in all relevant diagrams, thereby reducing the number of contractions required to compute the whole set of correlation functions compared to computing diagrams individually. This procedure is referred to as common subexpression elimination (CSE) for instance in compiler design.
A simple example is discussed in the appendix. In order to speed up this optimization process, duplicate diagrams are eliminated in a preprocessing step, which we refer to as diagram consolidation (DC). Duplicate diagrams are produced when a given momentum combination appears in the Clebsch-Gordan coefficients for more than one irrep (e.g. for the at-rest π(0)π(1)π(1) operators in irreps A − 1u and E − u ), or in the course of noise averaging (e.g. for the ground-state interpolator in A − 1u ). The number of elemental contractions required to com-  [50] for unexplained notation). Other parameters, e.g. for stout smearing [59], are the same as in [19]. . Given a list of diagrams to compute, it performs the global optimization described above and returns a serial list of elemental steps achieving operation count minimization [76].
Ensemble details: The results in this work employ the D200 ensemble generated through the CLS effort [60,61] with N f = 2 + 1 quark flavors and pion mass m π ≈ 200 MeV. The simulation belongs to a chiral trajectory where 2m l + m s is kept fixed as the light quark mass m l = m u = m d is lowered towards its physical value, implying m π > m phys π and m K < m phys K . The ensemble and measurement setup is detailed in Table  II. In order to ensure a Hermitian matrix of correlation functions despite the use of open boundary conditions in the temporal direction [62], our interpolating operators are always separated from the temporal boundaries by at least t s , where m π t s = 2.2.
Pion-pion scattering in the isovector channel has been investigated on this ensemble previously [19], and statistics subsequently improved considerably to provide spectroscopic information for the determination of the hadronic vacuum polarization on the same ensemble [63]. The pion functions can be re-used from that work, and hence no additional meson functions or solutions of the Dirac equation have to be computed.

Analysis strategy:
The procedure to extract the finitevolume spectrum from a matrix of correlation functions C ij (t) in a given irrep is discussed in detail in [15] and we use the analysis suite [77] developed in [19].
We  5a, t * = 10a), corresponding to roughly 0.32 fm and 0.64 fm in physical units [61], in order to extract not only the ground state but also excited states in most irreps. Results from different (t 0 , t * ) are indistinguishable, presumably due to the weak interaction in I = 2 and I = 3 pion scattering which results in little mixing of our interpolating operators, in which each hadron has been projected to definite momentum and is hence expected to overlap predominantly with a single state.
For two-pion states the difference ∆E between interacting and noninteracting energies is determined from single-exponential fits at sufficiently large time separations to the ratios (8) of diagonal elements of the 'optimized' correlation ma-trixĈ (i.e. the matrix formed from rotations by the eigenvectors of the generalized eigenvalue problem) and two single-pion correlation functions, and similarly for the three-pion states [38]. Absolute energies are reconstructed from those energy differences using the singlepion dispersion relation. Two-pion and three-pion spectra: The two-and threepion spectra with maximum isospin are extracted across a number of irreps with zero and nonzero total momentum. The attainable precision is generally at the few-permille level for the energies measured in units of the singlepion mass am π = 0.06504 (33). Figures 2 and 3 show the extracted two-and three-pion spectra together with the noninteracting energies, displaying significant energy shifts in all considered three-pion irreps. In particular, interacting energy levels from different irreps that contain some degeneracy of the noninteracting spectra (e.g. 1u and E − u at zero total momentum) differ substantially, which may suggest sensitivity to different combinations of low-energy scattering parameters.
FIG. 3: Same as Fig. 2 but for the I = 3 three-pion spectrum.
Only a dedicated investigation of those spectra in the framework of one of the available three-particle finitevolume formalisms can disentangle the effects of twoparticle scattering from genuine three-particle scattering effects, which necessitates further work along the lines of [35,41] to apply to the energies in all irreps presented here. In order to facilitate further investigation along these lines, the two-pion and three-pion spectra presented here are made publicly available, including all correlations. The values and covariance matrix of all extracted energies, as well as the single-pion mass, are given in Table VII, and the original bootstrap samples from this analysis are available as ancillary files with the arXiv submission.
The two-and three-pion excited state spectrum was previously predicted in [36], with input from the ground state energies at rest determined in a lattice calculation [38,39]. However, comparison with our results is difficult due to their use of a much smaller volume making their results subject to more significant finite-volume effects, especially at pion masses near m π ≈ 200 MeV where the exponential volume effects may become non-negligible.

CONCLUSIONS AND OUTLOOK
We have presented the I = 3 three-pion spectrum in finite volume from lattice QCD in which, for the first time, the excited states in various irreps at zero and nonzero total momentum, in addition to the ground states, have been extracted. These spectra need to be interpreted in the framework of one of the available three-particle finitevolume formalisms in order to extract infinite-volume information on three-pion interactions. In order to facilitate those investigations, which will require generalizations of the formulae currently available in the literature, all spectra are made public, including their correlations.
We also described a method, applied for the first time in lattice QCD, to reduce the computational re-sources required to construct correlation functions from the building blocks commonly encountered in the distillation method or its stochastic variant. These methods have been at the heart of much of the progress in lattice QCD studies of hadron interactions over the last years, and the optimization code discussed here can be used for arbitrary correlation functions. Our method is based on an algorithm devised in the context of quantum chemistry to speed up large numbers of tensor contractions and leads to a reduction in operation count of more than an order of magnitude for the set of observables considered here.
This algorithmic improvement paves the way to study more complicated systems requiring even more Wick contractions, such as the Roper resonance which has a sizable branching ratio for decays to N ππ as well as N π. It also reduces the computational cost associated with correlation function construction for baryon-baryon systems and will enable the lattice QCD investigation of nucleon as well as nucleon-hyperon interactions relevant for nuclear physics.
We Calculations for this project were performed on the HPC clusters "HIMster II" at the Helmholtz-Institut Mainz and "Mogon II" at JGU Mainz. Our programs use the deflated SAP+GCR solver from the openQCD package [62], as well as the QDP++/Chroma libraries [66]. We are grateful to our colleagues in the CLS initiative for sharing ensembles.
The interpolators used in this work are given in Tables III to VI for zero and nonzero total momenta. In practice, the normalization of an interpolator along each row of the tables is arbitrary. Correlation functions for equivalent total momenta are averaged, and the corresponding interpolators can be obtained using the reference rotations given in [56].
Simple example of the contraction optimization t f ts [2] [4] [3] (a) t f ts [2] [5] [3] Consider the list of two diagrams to compute shown above. They involve five baryon functions, i.e. rank-3 tensors, each with a definite momentum, irrep and noise combination labeled one to five. The most straightforward way to evaluate those two diagrams is to combine tensors on the source and sink times into outer products and perform the resulting contractions in both diagrams at a cost of 2N 6 dil . According to the first step of the algorithm discussed in the main text, a better way to evaluate the single diagram (a) proceeds as follows. Viewing pairwise contraction of two tensors over all common indices as the computational The step complexity is given by the number of summed indices plus the number of spectator indices. The first and second line are the locally good choices in this greedy algorithm, as they reduce the number of remaining contractions the most and at the smallest cost. Labeling the resulting rank-2 tensor from either of those as a new intermediary [6] transforms diagram (a) into a diagram with contractions between two rank-3 tensors and one rank-2 tensor remaining. Iterating this process shows that each of the two diagrams can be evaluated with complexity 2N 4 dil + N 2 dil . The second step of the algorithm discussed in the main text exploits the freedom to choose between the two locally good steps in order to reduce the global operation count. In this example, the contraction [2] - [3] can be re-used in diagram (b), whereas [1] - [4] cannot. Therefore the first step is more beneficial, and the corresponding subexpression, which only has to be computed once, is replaced with the intermediary [6] in both diagrams.
The two diagrams can thus be computed with complexity 3N 4 dil + 2N 2 dil , saving one of the computationally dominant contractions.

Covariance matrix
All energies and their uncertainties, as well as their normalized covariances, for the single-pion mass, the I = 2 two-pion spectrum and the I = 3 three-pion spectrum are given in Table VII. The single-pion mass in the first row is given in lattice units. The scale setting for the D200 ensemble employed here is discussed in [61]. All other energies are quoted in units of the single-pion mass. VII: Complete spectrum and normalized covariances. Different isospin sectors are separated by large row spacing into I = 1 (a single pion at rest), I = 2 and I = 3. The central values and uncertainties of center-of-mass energies Ecm/mπ in the second column are measured in units of the pion mass except for the single-pion energy in the first row, which is given in lattice units. The irrep labels are formatted as in Fig. 3.