Three-body unitarity versus finite-volume $\pi^+\pi^+\pi^+$ spectrum from lattice QCD

Strong three-body interactions above threshold govern the dynamics of many exotics and conventional excited mesons and baryons. Three-body finite-volume energies calculated from lattice QCD promise an ab-initio understanding of these systems. We calculate the three-$\pi^+$ spectrum unraveling the three-body dynamics that is tightly intertwined with the $S$-matrix principle of three-body unitarity and compare it with recent lattice QCD results. For this purpose, we develop a formalism for three-body systems in moving frames and apply it numerically.


INTRODUCTION
The dynamics of three-body systems above threshold play a key role to our understanding of strong forces. Many emblematic resonances exhibit significant three-body decay channels, such as the Roper resonance N (1440)1/2 + which, despite its low mass, couples strongly to the ππN channel leading to a very non-standard line shape and complicated analytic structure [1,2]. The ππN channels play also a significant role for other excited baryons and their description needs a quantitative understanding of three-body dynamics. Similarly, axial mesons like the a 1 (1260) and, supposedly, exotics decay into three particles [3].
The quantitative understanding of three-body systems in terms of QCD represents a long-term goal in hadronic physics. In lattice QCD (LQCD), the Hamiltonian is discretized and its eigenvalues are determined. These numerically demanding calculations are necessarily performed in a finite volume with periodic boundary conditions. This leads to a discrete eigenvalue spectrum in contrast to the continuous spectral density of scattering states in the infinite volume. These finite-volume effects are determined by hadron interactions and they offer a key to understanding these interactions arising from quark-gluon dynamics.
In this study, we compare the results of a recently developed infinite-volume mapping technique [4] with new finite-volume energy eigenvalues [5]. These data are calculated with multi-pion operators allowing for the reliable extraction of energy eigenvalues, above threshold and in different irreducible representation, providing, for the first time, access to three-body dynamics from first principles. Similarly to the case of the 2π + system, that represents the first physical application of the original Lüscher formalism [6][7][8][9][10], the 3π + system permits few partial waves and is an ideal system to study the pertinent finite-volume effects. This is a first step towards more complicated resonant systems that usually exhibit a complex pattern into two and three-body final states.
Recent progress in the three-particle sector is summarized in Ref. [11], see also Ref. [12] for a broader overview. In elastic two-particle scattering, each energy eigenvalue can be mapped to a phase shift [13,14]. However, the 3 → 3 reaction has eight independent kinematic variables (not including spin.) This requires a new formalism to map the discrete energy spectrum to infinite volume quantities.
Scattering amplitudes cannot be directly computed as infinite-volume limits of finite-volume observables. However, even without fully resolving the three-body dynamics explicitly, methods exist that take into account the contribution of three-body states [15][16][17][18]. These methods connect finite-volume data with infinite volume properties using either the optical potential or by extracting the spectral density from a correlator.
The first application of a three-body formalism to an actual physical system above threshold was achieved in Ref. [26]. Eigenvalues for the 3π + system were analyzed as calculated by the NPLQCD collaboration [58,59].
One of the problems in the lattice QCD calculation of energies for channels where three-body states are relevant is the need for many-hadron type operators to reliably determine the spectrum, as demonstrated, e.g., in Ref. [60]. Indeed, meson-baryon operators are often included in the operator basis [60][61][62][63]. Also, results on the Roper reso-nance at almost physical masses [60] suggest the need to map out finite-volume effects in two and three-body coupled channels, namely the πN, f 0 (500)N, π∆, ρN, . . . channels.
In the meson sector, lattice QCD results are available for channels where three-body states should be relevant [64][65][66], albeit only for pion masses and/or volumes at which the ρ meson can approximately be considered as stable. At lower pion masses, the three-pion spectrum requires three pion operators which has only recently been done [5].
In this study, we compute the excited two and threebody spectrum of the multi-pion system at maximal isospin and compare it to the calculation by Hörz and Hanlon [5]. The work is based on recent formal developments [4,29]; we use the Inverse Amplitude Method (IAM) to one loop [67][68][69][70][71][72][73] to predict the I = 2 pion-pion S and D-waves and then use the S-wave two-body input to predict the three-body finite-volume spectrum. Several eigenvalues are calculated in moving frames [5] which requires us to extend our formalism to boosted frames.
While this paper was prepared for publication, an independent study appeared [74] presenting a similar analysis of the π + spectra generated by Hörz and Hanlon [5].

FORMALISM
The three-body amplitude can be organized in the isobar-spectator picture; to describe three-body on-shell states, first, two particles are combined in terms of their quantum numbers and two-body interactions to form an isobar; the third particle, called spectator, is then added. Using this parametrization, a relativistic three-body unitary amplitude was derived in Ref. [75]. This provides a complete proof of three-body unitarity above threshold missing in previous work [76]. The amplitude is derived from dispersion relations, and can be matched to a Feynman diagrammatic approach but is a-priori independent of it. The isobar-spectator interaction itself is dictated by unitarity and develops an imaginary part. It can be represented as particle exchange as shown on the left-hand side in Fig. 1. There, solid lines indicate the spectator π + and double lines represent the isospin I = 2 isobar; note that any two-body amplitude, as for example the repulsive I = 2, = 0 can be mapped to the isobar picture [26,77]. In the present scheme, three-body forces arise naturally as real parts that can be added to the interaction without destroying unitarity.
For the 3π + system, the left-hand side of Fig. 1 indicates that the I = 2 ππ amplitude can only have evenspin isobars (S, D, . . . ) due to Bose symmetry. Also, the isobar-spectator interaction can only be in even partial waves = 0, 2, . . . . For both cases, we truncate the expansion at D-waves which is a good approximation for low energies that is backed by phenomenology [78].  [79][80][81]. For comparison, the predicted D-wave at the pion mass of Ref. [5] is also indicated (blue curves/area). The respective elastic regions are indicated with the horizontal bars. Predictions are shown using the LECs from Ref. [82] (GW), Ref. [83] (GL), and Ref. [68] (DP).
In addition, the (I, ) = (2, 2) ππ interaction is very small as shown in Fig. 2 for different low-energy constants (LECs). Perturbative next-to-leading order (NLO) calculations (red curves and band) predict very small phase shifts not in contradiction with the scattered phase shifts from experiment. See also Ref. [73] for a similar calculation, comparing also the LQCD phase shifts of Ref. [84]. If one chirally extrapolates the calculation to the pion mass of Ref. [5], of m π ≈ 200 MeV (blue lines and band), one can see that the size of the D-wave stays below one degree in the elastic region. We also find that there is no apparent sign of D-wave in the lattice data under consideration [5]. In all irreducible representations ("irreps") in which the D-wave is the lowest participating wave, the finite-volume energies coincide with noninteracting levels within uncertainties. For irreps with S/D-wave mixing, no apparent sign of D-wave is found, either, as discussed in the Results section. For our predictions, we will, therefore, neglect the (I, ) = (2, 2) ππ interaction in the following. How-ever, there is no reason to exclude the relative π + -isobar D-wave which will turn out to be important.
In Ref. [4] the finite-volume version of the three-body amplitude was derived, and, for the first time, the systematic cancellation of unphysical singularities and the practical applicability of a such a formalism was demonstrated, and correctly projected to the A − 1 irrep. In Ref. [26], for the first time, a three-body formalism was compared to LQCD data of an actual physical system, π + π + π + , including a fit of the three-body force. In summary, the only missing ingredient for the prediction of the new LQCD data consists in the development of a finite-volume formalism allowing for three-body systems in moving frames.

Moving Three-body System
For the formulation of the three-body T -matrix in finite volume [4], we took advantage of cubic symmetry which enabled us to arrange allowed lattice momenta on "shells" of equal absolute momenta. For moving threebody systems, cubic symmetry is broken and it is more advantageous to work in a three-dimensional momentum basis, suitably labeling the allowed momentar i = (2π/L)ñ i withñ i ∈ Z 3 . In the following, three-momenta with tilde are defined in the lattice rest frame, threemomenta without overscript are defined in the threebody rest frame, and starred three-momenta are defined in the two-body isobar rest frame. With this, the symmetrized three-body scattering amplitude in the threebody rest frame reads where (n,n,n) denotes a circular permutation of (1, 2, 3) etc., and v denotes the decay vertex of the isobar, which is chosen as specified in the appendix to reproduce exactly the Inverse Amplitude Method for the two-body sub-channel amplitudes [26]. Note that this vertex also contains a smooth cutoff function which regulates all two and three-body integrals. This function is chosen as in Ref. [26], where it is shown that the dependence on the particular choice of the cutoff is very weak. The quantity s represents the square of the total four momentum of the three-body system, such that the isobar-spectator amplitudeT reads [4] T nm (s) = τ n (s)T nm (s)τ m (s) − 2E n L 3 τ n (s)δ nm , (2) where m, n, x label the incoming spectator momentum p m , outgoing spectator momentum q n , and intermediate spectator momentum l x . A graphical representation of the second ("rescattering") term of Eq. (3) is given on the right-hand side of Fig. 1. Furthermore, E n = m 2 π + q 2 n and analogously for the other momenta. The Jacobian for the mapping from the lattice frame to the threebody rest frame is denoted byJ x . The quantities B ("exchange+three-body") and τ ("three-body propagation") are defined in the Appendix and graphically indicated in Fig. 1 to the left (up to a real three-body term, see Appendix). The projection to irreps is detailed in the Appendix, as well. In summary, for incoming and outgoing spectator momentap i andq i , the 3π + system has momentumP = q 1 +q 2 +q 3 =p 1 +p 2 +p 3 . A boost of lattice momenta bỹ P provides the three-momenta entering Eqs. (2, 3) that is solved in the three-body rest frame; another boost to the isobar rest frame is necessary as the pertinent summations are carried out in that frame. Schematically, we can represent this two-step process as follows: Lattice rest frame 3-body rest frame 2-body rest frame As a result of the formalism (see Appendix), the threebody system in moving frames is entirely expressed in terms of lattice momentap m ,q n ∈ (2π/L)Z 3 , and the invariant s. Its poles indicate the energy eigenvalues after projection to irreps.

RESULTS
Taking the two-body input from IAM [67][68][69][70][71][72][73] we predict the energy eigenvalues for the π + π + and π + π + π + systems. The uncertainties for our prediction are estimated by using central values for LECs from three different analyses [73,82,83] instead of LEC uncertainties. The reason is that, usually, no correlations on LECs are quoted in the literature which leads to an uncontrolled overestimation of the prediction error. For the LECs of Ref. [83], results are quoted in Tab. I. For all results, the D-wave isobar is neglected as discussed before, and the three-body term is set to zero, C = 0 (see discussion below). The predictions for the two-body (three-body) spectrum are represented in the upper (lower) part of Fig. 3. For some of the irreps, phases are extracted and shown together with chiral predictions in Fig. 4 for illustration. Overall, the predictions from different LECs vary surprisingly little given the different origins of their determination. Furthermore, the predictions are all quantitatively very good. In all A + 1 irreps with non-vanishing boost, S and Dwaves mix in the π + π + system. At higher energies, one could therefore expect deviations of the predictions from the data as we have neglected the D-wave throughout. However, the quality of our predictions, even beyond the elastic threshold, adds another piece of evidence that the D-wave can be neglected. Of course, only a S/D-coupled partial-wave fit can provide ultimate clarity for this point (see, e.g., Refs. [84,88]). In Ref. [74], some evidence for a non-vanishing D-wave was found by fitting only irreps in which the lowest participating wave is the D-wave.
The quality of our predictions can be assessed by evaluating the correlated χ 2 /n with n being the number of lattice eigenvalues in the respective elastic regions. For LECs set to the values of Ref. [83] we have n ≈ 39.4 22 for the two-, three-, and combined two and three-pion sectors, including the cross-correlations of energy eigenvalues. The χ 2 values for the LECs from Refs. [73,82] are very similar to the quoted ones.
The predictions were obtained with a vanishing threebody force C. In Ref. [26], C was fitted to the groundstate level and found to be negligible. However, the present data [5] are more precise than the NPLQCD data [59]. Any deviation of the prediction in the threebody sector, especially at higher energies, could be a sign for a non-vanishing three-body force at the chosen regularization. This is obviously not the case as χ 2 (3) /n ≈ 1. Moreover, a large part of the overall χ 2 (2&3) ≈ 39.4 arises from the correlations of one point at low energies (σ 1/2 ≈ 2.4 m π , π + π + sector, A + 1 (1)) with the 3π + sector. Without this point, χ 2 (2&3) /n ≈ 28/21 ≈ 1.3 and it is difficult to explain this change with the discussed simplifications of our formulation. To conclude, consider the excited-state 3π + energy shifts in A − 1u (0) and E − u (0) highlighted in Fig. 3. The relative and absolute sizes of these shifts are governed by the structure of the exchange term B shown in Fig. 1 because that term determines the strength of S-wave vs. D-wave interactions. On the other hand, that term arises as a consequence of three-body unitarity [89] and contributes to the powerlaw finite-volume effects [4]. In conclusion, for the first time, three-body unitarity is directly visible in LQCD data. This conclusion would hold similarly for any parametrization of the two-body sector, i.e., it is independent of the IAM model we choose for our predictions.

CONCLUSIONS
Using a two-body unitary amplitude, that matches Chiral Perturbation Theory up to next-to-leading order (IAM), the isospin I = 2 two-body body eigenvalues of a recent lattice QCD calculation [74] were predicted in a restriction to S-wave. With this two-body input, three-body unitarity served as the S-matrix constraint to predict the three-body spectrum with a correlated χ 2 (3) /n ≈ 10/11, i.e., no sign of a substantial three-body force was seen for the given regularization. Yet, if correlations of the two-and three-body sector are combined, a χ 2 (2&3) /n ≈ 1.8 indicates a residual tension. We want to stress that the LECs are not fit to the lattice data; the tension is likely to disappear if we adjust the LECs to minimize χ 2 . Overall, the predictions, depending only on low-energy constants from independent studies (and, very weakly, on the regularization), are in good agreement with the data. Furthermore, the correct prediction of the S-wave and D-wave excited-level energy shifts in A − 1u (0) and E − u (0) depends only on the structure of the spectator-isobar interaction, which, in turn, is dictated by three-body unitarity. For the first time, this fundamental S-matrix principle is directly visible in lattice QCD data. to the fact that we use the LECs extracted in the latter scheme. Further details on this technicality are discussed in the Ref. [26]. Overall, the above choice of the coupling λ leads to the form of the two-body sub-channel amplitudes, which match the Inverse Amplitude Method [67,90]. This type of amplitudes matches the ChPT amplitude up to the next-to-leading order exactly, allowing also for addressing all three isospin channels of the ππ system in a large energy region as recently demonstrated in Ref. [82].
For completeness, we also quote the (unprojected) driving term where E 2 ex = m 2 π + (q n + p m ) 2 . With all parts of Eqs. (2) of the main text defined, theT matrix can be calculated; its poles coincide with the three-body energy eigenvalues in moving frames.

Predicted Energy Eigenvalues
In the two-body sector, the positions of the poles of the two-body scattering amplitude T 22 = vτ v give the two-body energy eigenvalues in the A + 1 irrep. In the three-body case, the projections to the corresponding irreps is performed similarly to the method of Refs. [5,91], see Ref. [29]. In particular for Eq. (1) of the main text, where the indices i and j run over all group elements and the coefficients χ are the characters of the group elements, see, e.g., Refs. [5,91]. Here, Γ denotes the irreps A − 1u , E − u , A − 2 , B − 2 , and E − . The predicted energy eigenvalues are shown in Table I. two-body three-bodỹ