Three-body interactions from the finite-volume QCD spectrum

We perform a fit of the finite-volume QCD spectrum of three pions at maximal isospin to constrain the three-body force. We use the unitarity-based relativistic three-particle quantization condition, with the GWUQCD spectrum obtained at 315 MeV and 220 MeV pion mass in two-flavor QCD. For the heavier pion mass we find that the data is consistent with a constant contact term close to zero, whereas for the lighter mass we see a statistically significant energy dependence in tension with the prediction of leading order ChPT. Our results also suggest that with enough three-body energy levels, the two-body amplitude could be constrained.


I. INTRODUCTION
It is a long-term quest of nuclear physics to understand hadron interactions as they emerge from quark and gluon dynamics. The main challenge lies in the fact that perturbation theory fails at low energies, because the interactions are strong. Lattice QCD (LQCD) offers a non-perturbative method which has quarks and gluons as fundamental degrees of freedom while keeping all systematics under control. LQCD calculations are performed in a finite volume and in Euclidean time, leaving only indirect methods to study real-time infinite-volume scattering. The relation between finite-volume spectrum and infinite-volume scattering amplitudes is called quantization condition, which has been known for two-hadron systems since the pioneering work of Lüscher [1]. The last decade has witnessed significant progress using this approach for a variety of interacting two-particle systems. Only recently has the quantization condition been extended to the three-hadron sector.
Recent theoretical advances  of the three-body formalism as well as related numerical studies [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] (see Refs. [50,51] for reviews) open a possibility for studying hadronic processes that involve three-body interactions from first principles. Many scattering channels of interest receive contributions from three-particle states. For example in the meson sector the a 1 (1260) resonance, seen experimentally in τ decays, couples first to ρπ and σπ, and then due to the instability of the ρ and σ mesons to a final state with three pions. In the baryonic sector an example where three-body states are relevant is the Roper resonance N (1440)1/2 + which has both two-and three-particle final states, as it decays to πN and ππN .
Finite-volume formalisms for excited states were developed later. In Ref. [62] different approaches are reviewed in more detail, including a timeline and a comparison of the methods. In the following we provide a very brief overview. Among other recent developments [28,47], there are two relativistic formalisms available for threeparticle scattering usually referred to as relativistic finitevolume unitarity (FVU) [16] and relativistic effective field theory (RFT) [9] approaches. The former uses unitarity of the S-matrix as a guiding principle, while the latter relies on the re-summation of diagrams. A key ingredient to these formalisms is the parameterization of two and three-body scattering amplitudes. The two-body input can be determined from experiment in combination with chiral extrapolations, but lattice QCD data, together with the quantization conditions, can also be used.
The first determination of the three-body force from the NPLQCD data [52] using such a formalism was achieved in Ref. [40] with the FVU framework. There, within the uncertainties of the lattice data, the three-body force was found to be zero. This study contains also the first predictions for excited levels at different pion masses. Later, in Ref. [54] excited levels for different boosts and irreducible representations (irreps) were calculated that led to a refined determination of the three-body force with the RFT formalism [43]. In particular, the three-body force was found to be non-zero, and even attempts to determine its energy dependence could be made. Within the FVU formalism [44] the data of Ref. [54] were predicted using only chiral extrapolations of two-body input up to next-to-leading order (NLO) and assuming a vanishing three-body force.
The GWUQCD collaboration calculated the three-π + spectrum for different quark masses, box geometries, and boosts, mapping out a plethora of states and comparing also to FVU predictions that were made under the assumption of vanishing three-body forces [55]. The agreement found was fair. Subsequently, the ETMC collaboration calculated the three-π + spectrum at three different pion masses, including the physical point for the first time [56]. The extraction of the three-body force with the RFT formalism was compared with the leading order (LO) Chiral Perturbation Theory (ChPT) prediction for the three-tothree process. Similar to Ref.
[43], the three-body term was found to be non-zero, its energy-independent part being in qualitative agreement with LO ChPT, but its energy-dependent part not.
The evaluation of the infinite-volume three-body amplitude is a challenge by itself. For example, the FVU framework was recently extended to the infinite volume (albeit without lattice input), in order to study the decay a 1 (1260) → πρ in coupled S-and D-waves [63]. For three positive pions, the infinite-volume scattering amplitude was solved in Ref. [45] in the RFT framework, for the first time with actual lattice input. The input was calculated by the Hadron Spectrum collaboration and the analysis included an extraction of the three-body force. Within uncertainties and different fit strategies/parametrizations tried, the three-body term K 3,iso was found to be compatible with zero.
Obviously, the role of three-body forces is not settled. In this paper, we continue the investigation of the three-body force with the FVU formalism. In particular, the mentioned discrepancy of FVU predictions with GWUQCD data (χ 2 dof ≈ 2.68) [55] could be the result of a finite three-body term. Therefore, we perform here a fit to the lattice data of Ref. [55] using the FVU formalism with special focus on the three-body contact term. In Section II we review the extraction of the three-pion finite-volume spectrum from the GWUQCD collaboration. Then in Section III, we review the three-body quantization condition. We include technical implementation details such as parameterizations of the two-and three-body scattering amplitudes, and establish a connection to a three-to-three contact term. In Section IV we present the results of various fit scenarios. Finally in Section V we discuss the impact of the results.

II. FINITE-VOLUME QCD SPECTRUM
Here we review briefly the details of the calculation of the finite-volume spectrum performed by our collab- oration, referring the reader to Ref. [55] for additional material. The ensembles employed use N f = 2 dynamical fermions, with two sets of quark masses tuned so that the pions have masses of 315 MeV or 220 MeV. These data sets were used to compute the two-pion spectrum in all three isospin channels [64][65][66][67], the three-pion I = 3 channel [55] and the K − K − K − channel [57]. The quark propagators are estimated using LapH smearing [68] and an optimized code is employed to compute the required matrix inversions [69]. The parameters describing the ensembles in this study are listed in Table I. The lattice spacing was tuned using the Wilson flow parameter t 0 . For details about tuning the lattice spacing and other parameters we refer the reader to Ref. [70].
For each pion mass, we calculated the I = 3 threepion finite-volume spectrum on three different ensembles, which feature two different geometries. There is one cubic volume and two volumes which are elongated in the zdirection, which is also the direction of the boost when calculating the energies of moving states. The advantage of using elongated ensembles is an increase in energy resolution, since the momenta are quantized in units of 2π/L. For ensemble E 3 , we are able to extract 16 energy levels below the inelastic threshold, compared to 3 in the cubic volume at the same pion mass. The numerical cost of generating these ensembles is reduced since the volume increases linearly with the elongation, as opposed to a cubic increase for symmetric boxes. In total we have 30 energy levels below 5m π across the six ensembles listed in Table I. They form the basis for our quantitative analysis of the three-body interaction in this work. Precise values for the levels can be found in Appendix B of Ref. [55].
Regardless of geometry, the finite volume introduces a critical change to the angular momentum quantum number of the system of interest. In the infinite volume, the symmetry group of rotations is SO (3), and the quantum number of the states are labelled by the angular momentum l, the irreducible representations (irreps) of the group. For a finite volume the symmetry group is reduced from SO (3) to O h in the case of cubic volumes, and D 4h in the case of elongated volumes. If the system is at non-zero total momentum the relevant symmetry group is C 4v , provided the boost is along the direction of any spatial elongation. In this study we only consider a boost of P = [001], which is aligned with the elongation. The rotation quantum number of states in a finite volume are thus labelled by the irreps of the respective symmetry groups. Here we focus only on states that overlap with l = 0, which will be labelled by A 1u for symmetry groups O h and D 4h , or labelled by A 2 for the C 4v symmetry group. For a full discussion on the relation between angular momentum l and irreps of the finite-volume symmetry groups in elongated boxes we refer the reader to Ref. [71].

III. QUANTIZATION CONDITION
The dynamics of a (multi-)hadron system is accessed in lattice QCD by calculating correlation functions on a discretized Euclidean space-time in a finite volume. The so-called discretization effects are related to the finite lattice spacing a. In principle, a continuum limit a → 0 needs to be taken to relate the lattice QCD results to the physical (real-world) quantities. Since only a finite number of lattice sites can be considered in any practical calculation, the system of interest is necessarily evaluated in a finite volume. Imposing boundary conditions restricts the allowed momenta in this system. For example, in a cubic volume of side length L, the frequently applied (also here) periodic boundary conditions yield an infinitely large discrete set of allowed momenta S L = {(2π/L)n|n ∈ Z 3 }. Unavoidably, using threemomenta from such a set changes the infinite-volume spectrum fundamentally, making it a discrete set of energy eigenvalues. This holds for any finite L, such that a simple extrapolation (lim L→∞ ) is not useful, calling for a non-trivial mapping between infinite-and finite-volume. Such mappings are referred to as quantization conditions. An alternative possibility to use ordered double limit [72] techniques to extract (complex-valued) amplitudes was explored in Ref.
As one of the currently available relativistic formulations of the three-body quantization condition (see Refs. [50, 51] for reviews) the FVU approach derives from the relativistic unitary three-body formalism [76] in infinite volume. The formalism differs from the diagrammatic approach followed in Ref. [9], but yields an equivalent finite-volume spectrum given the same input, as shown using time ordered perturbation theory in Ref. [34]. Also in the infinite volume both formalisms yield a unitary three-body amplitude [77,78]. Both formulations are currently being applied to a variety of calculations of simpler three-hadron systems such as 3π + [40, 43-45, 55,79] or 3K − [57]. However, the practical implementation of Schematic representation of a two-plus-one decomposition of the three-body system. P , p/q denote three-momenta of the three-body system and that of the in/outgoing spectators, respectively. Total three-body energy squared and the invariant mass of the two-body systems are denoted by s and σ p/q , respectively. each differs. See also Ref.
[47] for a calculation based on combination of a variational approach and the Faddeev method, including relativistic effects for the pions and lattice spacing effects. In the following, we will demonstrate and discuss an example of contrasting representations of the three-body contact term.

A. Design and implementation
We avoid discussing the derivation of the FVU formalism here [16, 40, 44, 55], but review the results and unify the notation. At its core, the condition derives from a relativistic unitary three-body scattering amplitude, which resolves the three-body dynamics as a cluster of a twobody (also related to 'isobar' or 'dimer' notation) state and a third particle -a 'spectator' [76] (see also Refs. [80][81][82] for recent progress in this direction). The kinematic notation of such a configuration is depicted schematically in Fig. 1. Unitarity constrains the correct interplay between the two-and three-body interactions accounting for on-shell configurations, which are the source of singularities in the finite-volume correlator. The net result is that in a finite volume √ s is the energy corresponding to an interacting three-body state when The determinant is taken over the S Lη × S Lη space, referring to the momentum of the in-and outgoing spectator (third particle).
Here, S Lη = {2π/L(n 1 , n 2 , n 3 /η)|(n 1 , n 2 , n 3 ) ∈ Z 3 } refers to the momentum configuration from elongated boxes used by GWUQCD andP is the total three-momentum of the system. The remaining building blocks of this condition (∈ Mat S Lη ×S Lη ) are ?
This singular and non-diagonal matrix originates from the one-particle exchange contribution, and is a direct consequence of three-body unitarity in the infinite volume. Here and in the following E x := x 2 + m 2 π denotes the on-shell energy of a pion with a momentum x.
. This singular, diagonal matrix accounts for the on-shell configurations in the two-body subsystem with total energy squared σ p (s) = ( √ s − E p ) 2 − p 2 . The equation above is derived in the Appendix with the final form given in Eq. (A29). Note that the summation is performed overk in the lattice frame, while the summands are expressed in terms of the k * momenta in the two-body rest frame. The corresponding two-step boost reads The explicit formulas for the boosts and corresponding Jacobians J andJ are collected in Appendix A. ?
This regular and in general non-diagonal matrix parametrizes the isobar-spectator interaction and, thus, also the three-body contact interaction via a non-trivial mapping discussed in the Sec. III C. In both cases, this contribution is not fixed and needs to be obtained from a fit to the lattice results as will be discussed in Sec. III C and IV. Note that the normalization used here differs slightly from previous work [16,76].
. This regular, diagonal matrix parametrizes the dynamics of the two-body sub-system. In that, it corresponds to the usual K-matrix as explained in detail in Appendix A 2. Similarly to the three-body term C 0 , this contribution needs to be determined from a fit to the lattice eigenvalues. Such a fit can take two-body, three-body, or both data types into account as will be discussed in Sec. IV.
Before proceeding with the discussion of the two-and three-body interaction terms K −1 and C 0 , we note a generic feature of the three-body quantization condition. In contrast to the well-established Lüscher's method in the two-body case, the three-body quantization condition in Eq. (1) emerges as a full-fledged determinant equation, even in the simplest case of identical particles in S-wave. Thus, it prevents any one-to-one mapping between infinite-and finite-volume quantities. Instead, even for the simplest three-body systems, one must fix the volume-independent quantities (C 0 and K −1 ) from a fit to lattice data and then use those to evaluate the Leading contributions to the K-matrix utilized for the two-body input of the quantization condition in Eq. (1). Diagrams represent scattering amplitudes obtained from ChPT with dots and squares denoting the leading and next-to-leading chiral order vertices, respectively (tadpole contributions not shown).
infinite-volume scattering amplitude. Note that this threeto-three amplitude provides only the three-body unitary final state interaction of production processes; to actually obtain mass-spectra or Dalitz plots as, e.g., in Ref. [63], one needs additional information.

B. Parametrization of the two-body input
Depending on the system at hand and the research objectives, various techniques in parametrizing the two-body dynamics (K −1 ) may be more or less advantageous. For example, a model based on ChPT was used in Ref.
[40] to bridge the lattice results obtained at different pion mass. Another approach, related to an effective range expansion was employed in Ref.
[43] allowing for a systematic extraction of threshold quantities.
With the lattice data spanning over large energy ranges and two different pion masses we proceed here with the path traced out in Ref. [40], and implement the modified Inverse Amplitude Method [83][84][85][86] (mIAM) into the three-body formalism. This is also motivated by our previous applications to the isoscalar channel [66,87] and more recently cross-channel (I = 0, 1, 2 ππ scattering) analysis [88] of GWUQCD lattice results [65][66][67] which is based on this approach.
A practical implementation of this requires equating the two-body scattering amplitude used in the infinite volume analog of Eq. (1) to the mIAM amplitude [86]. This procedure has been developed in Ref.
[40], but is improved in the current study. In particular, the matching procedure is now exact and without the need of a redundant regulator. Details can be found in Appendices A 2 and A 3 leading to (see Eq. (A17) there). Here, T n denotes the chiral ππ amplitude of order n projected to the isospin I = 2, angular momentum = 0 partial wave. The barred symbol denotes a reduced amplitude, where the s-channel one-loop diagram is subtracted. This algebraic re-shuffling ensures that the above expression indeed contains only quantities which do not go on-shell in the physical region. Similarly, A(σ) is a function of the chiral amplitudes, introduced in Ref. [86] to improve the analytic behavior of the (infinite volume) scattering amplitude below threshold 0 < √ σ < 2m π . In particular, adding A(σ) prevents the two-body scattering amplitude (∼ 1/(T 2 − T 4 )) from diverging around the Adler zero (σ ⇔ T 22,∞ (σ) = 0, where T 22,∞ is the S-wave scattering amplitude). Recall that the presence of the Adler zero itself is demanded by chiral symmetry. The diagrammatic representation of the leading contributions is given in Fig. 2, while formulas are provided for convenience in Appendix A 2.
Note that the two-body subthreshold amplitude enters the three-body quantization condition Eq. (1) since for somep, σ p(p,P ,s) < 4m 2 π forp ∈ S Lη . The mIAM approximates the left-hand cut by the next-to leading chiral order. Thus, it is expected that at too low invariant masses the latter model is a poor approximation of reality. To account for this, we have explored various possibilities, such as employing smooth form-factors, etc.. We found that fixing is sufficient. Variating the matching point in a reasonable range leads to uncertainties orders of magnitude smaller than other effects, see Ref. [57]. Also, changes in sub-threshold behavior are supposed to be absorbed in three-body contact terms.
Overall, the expression in Eq. (2) is a regular, volumeindependent function of two-body energy σ, ensuring in infinite volume an exact matching of the two-body dynamics of the mIAM approach. In that, Eq. (2) depends on four renormalized low-energy constants {l r 1 , l r 2 , l r 3 , l r 4 }, see Ref. [89] for explicit formulas relating those to Lagrangian constants. As it has been shown in previous studies of twobody lattice results [87,90], l r 3 and l r 4 contribute weakly to the two-body ππ dynamics. This is simply due to the fact that their appearance is solely rooted in the procedure of replacing the Lagrangian quantities -pion decay constant in the chiral limit and leading order pion mass -by their physical/lattice values. Thus, we simply fix the symmetry breaking LECs {l r 3 = 8.94 · 10 −6 , l r 4 = 9.05 · 10 −3 } to their typical values [88,91] at the regularization scale of µ = 770 MeV. Another possibility is to directly use both quantities as fit parameters as introduced in Ref. [92]. This has an advantage of avoiding the scale setting discussion, which is beyond the current stage of three-hadron spectroscopy from QCD.

C. Parametrization of the three-body input
The implementation of the three-body input into the quantization condition Eq. (1) is much less explored than that of the two-body. One of the goals of this paper is to improve on this, in exploring various ways of parametriz-ing three-body dynamics theoretically and in a practical application to the lattice data.
In the three-body finite-volume formalism (FVU) [16], the three-body dynamics is included via a regular (realvalued, infinite volume) function -the isobar-spectator contact term C ∼ C 0 . This term is not fixed by unitarity and, other than demanded by data, there are no constrains on its functional dependence with respect to energies and masses. The RFT formalism [9] on the other hand, relies on the diagrammatic counting of the on-shell configurations in derivation of the three-body quantization condition. Thus, it includes by construction a regular three-body contact term (K 3,df ), which again needs to be determined from a fit to the data.
Interestingly, in the pioneering application to the lattice 3π + finite-volume spectrum [40], a simplistic fit 1 to the NPLQCD data [52] Later, an analysis of more recent data at one pion mass [54] was performed within both RFT [43] and FVU [44]. In that, the FVU formalism with vanishing three-body term C 0 led to a χ 2 dof ≈ 0.9 prediction of threebody energies, while a non-zero m 2 π K 3,df = 270(160) was obtained in the RFT approach by fitting two-and threebody lattice energies. The apparent discrepancy between these results was raised in the community, but, so far, not resolved quantitatively.
To find an algebraic connection between the three-tothree contact term (T 3 ) and the isobar-spectator coupling C 0 we utilize the language of Ref. [16,76] in matrix notation for convenience. There, a fully connected threebody amplitude takes the form for T IS , τ and v denoting isobar-spectator amplitude, isobar propagator and its coupling to asymptotically stable states, respectively. Explicit definitions of the latter do not matter for the derivation but can be found in Refs. [16,76] and Appendix A 3. For simplicity, momentum and energy labels are omitted. Symmetrization over external momenta is taken into account by ... . For three degenerate scalars in (relative) S-wave this leads to Taking now the limit m π → ∞ for all intermediate pions 1 Note that projection to shells (sets of momenta related by octahedral symmetry transformation), was employed in this work, which is not part of the more general form of Eq. (1). There the employed fit form was chosen as C 0 = 1 · c 0 (mπ/m phys π ) 2 . provides an expression equivalent to a contact term (T 3 ) or equivalently whereT 3 denotes a real three-body contact term. This equation is schematically illustrated in Fig. 3 to the top. Note again that appearance of factors (32π) is caused by the fact that we are working in the plane wave basis, see Appendix A 2 for more details. This relation is of the same form as the relation between the K 3,df and C 0 term, which can be obtained from a matching of FVU and RFT formalisms, see for example Refs. [50,51]. Also, as noted there, it implies that in general an isotropicT 3 leads to anisotropic C 0 and vice versa.
To expand on this further, we consider the following example. Chiral perturbation theory at leading chiral order (this was used for the RFT formalism in Ref. [43]) yields for our formalism a three-to-three contact term of the formT Using now Eq. (6) with the K-matrix from Eq. (2) we obtain a prediction for the isobar-spectator interaction C 0 ( √ s, p, q), for the momenta belonging to the in/outgoing spectators. The resulting symmetric (in p, q) isobar-spectator contact interaction is depicted for several values of total three-body energy 3m π < √ s < 5m π in Fig. 3. We observe that several orders of magnitude difference between the overall scales of theT 3 and C 0 (e.g. C 0 ≈ 10 −6 MeV −2 ↔ m πT3 ≈ 10 2 ) occurs naturally, connecting the results of Refs. [43, 44].

D. Three-body state energies
Given parameterizations of the two-and three-body interactions, the finite-volume spectrum can be determined by searching for energies at which Eq. (1) is satisfied. To find the energies associated with a particular irrep Λ of the symmetry group G, we first block-diagonalize the matrix Q Q = diag(Q Λ1 , Q Λ2 , . . . ), The determinant of Q Λ block can then be scanned separately in s for the zeros which determine the energies finite-volume states for that irrep.
To block-diaganolize the matrix Q we need to change to an appropriate basis. The change of basis from the plane-wave basis in Eq. (1) can be done by constructing projectors for the row λ of the irrep Λ of the group G: where d Λ is the dimension of the irrep Λ, Γ Λ λλ (R) is the representation matrix of the group element g in row λ of the irrep Λ, R(g) is the rotation corresponding to g, and U (R(g)) is the unitary transformation of R on the plane-wave space.
We must truncate the plane-wave space S Lη to include only momenta below a sufficiently large threshold, p max ≈ 1 GeV. In general the values of the parameters would need to be adjusted based on the choice of the cutoff p max . In our case we found that there is a very mild dependence on the choice of p max , for the energy levels predicted from the quantization condition, well below the stochastic errors in the lattice data. By construction this truncated space is invariant under the symmetry transformations, since the momenta magnitudes are invariant under rotations. For each irrep Λ we project the entire plane-wave basis using the projector above and the dimensionality of this space represents the number of independent multiplets associated with Λ that appear in the selected plane-wave basis. The restriction of Q to this subspace is Q Λ which has zeros corresponding to energies in the irrep Λ.

IV. RESULTS
In this section we present the results for the three-body terms as extracted from fitting the finite-volume spectrum. For sake of clarity we will first discuss the extraction of the three-body terms when the two-body dynamics is fixed by the mIAM parameterization with LECs determined from fitting our lattice two-pion spectrum in all isospin channels [44]. This model provides a good description of our two-pion spectrum across all channels for the two quark masses we studied. For the three-body terms we parameterizeT 3 , using no energy dependence:T 3 (s) = t 0 , orT 3 (s) = t 0 + t 1 s/m 2 π , a linear function in s. The second parametrization is consistent with the leading order prediction from ChPT in Eq. (7).
The data points included in these fits are all three-body energies below the inelastic threshold in the irreducible representations sensitive to the s-wave three-body terms. The relevant representations are A 1u for the states at rest and A 2 for the moving states. The total number of points is 12 (2, 3, and 7 for ensembles E 1 , E 2 , and E 3  Dark bands indicate 1σ uncertainties for fits of the three-body data using fixed mIAM LECs from Ref. [88]. Light bands indicate the 1σ uncertainties from fits to the three-body data including l r 1 and l r 2 as fit parameters and priors from the cross-channel two-pion fits as described in Sec. IV. respectively) for the heavy quark mass and 7 (2, 2, and 3 for ensembles E 4 , E 5 , and E 6 respectively) for the light quark mass. These energy levels are plotted with errors in Fig. 4.
All fits to the three-body energy levels are performed separately for the two quark masses since our parametrization for the three-body term does not constrain the quark mass dependence. We perform various fits, keeping some of the fit parameters fixed and varying others as shown in the Table III. Below we discuss these results.
To get a sense whether these parametrizations are reasonable, we first extractedT 3 (s) by fitting each individual lattice energy level alone. This offers a rough profile of T 3 as a function of energy. To stabilize these fits we fix the LECs to the values extracted in our cross-channel two-pion study [44]. In Fig. 5 these results are represented by the data points, for the heavy and light quark masses. We note that while the errors are large the results sup-portT 3 (s) being (weakly) dependent on energy. For the heavier quark mass the data is consistent with a constant parameterization forT 3 , close to zero, whereas for the lighter mass we see a statistically significant fall in the data with energy. Note that the magnitude ofT 3 in absolute terms is not that different, but the energy dependence is more pronounced for the lighter mass. Additionally, for the light massT 3 is well constrained to be non-zero away from threshold. We also present in these plots the expectation from the leading order ChPT, plotted with a dashed line. At the light mass there is some tension with the leading order ChPT prediction, in particular at higher energies. Also, the energy dependent term has opposite sign to the ChPT result. This is similar to the tension in the energy dependence of the RFT three-body term in Refs. [43,56].
The simultaneous fit to all the energy levels are indi-cated with color bands in Fig. 5. For these fits we allowed the LECs to vary but constrained their variation using priors based on the probability distribution that was determined in our cross-channel two-pion data fits. We prefer this strategy over the simultaneous fit of the twoand three-body energies, since it makes the results for the three-body fit easier to analyze by isolating the two-body contribution. These priors are included by augmenting the χ 2 -function that is minimized by the fitter: where δl i := l r i −l i , withl being the values of the LECs as determined in Ref.
[44] and Σ ij their covariance matrix. Note that since the LECs variations tend to be strongly correlated, in order to include these priors properly we need to consider the full correlation matrix Σ, not just the diagonal terms associated with individual LECs error estimates. In Ref.
[44] we only reported the errors on the LECs. We include the estimate for the covariance matrix in Table II. The constant fit works reasonably well for the heavy quark mass ensembles and is consistent with the linear fit in s. For the lighter quark mass the linear and constant fit are quite different and the linear term is needed to describe the data well. The values extracted for the fit constants and the χ 2 value for the fit are listed in Table III.  [88], others are left as free parameters of the fit. The final row in each table is for a fit using relaxed priors, as indicated by an asterisk and described in the text.
We also indicate in the figures with narrower bands the fit results for {t 0 , t 1 } when the LECs are not allowed to vary, corresponding to the second and third rows in Table III. We see that while the central values are almost the same, the error bands are almost doubled when we allow the LECs to vary. This shows that even with small error bars (few percent level) in the two-body parametrization, there is a large impact on the error of the three-body terms in this channel. This is partly due to the smallness of three-body terms. The smallness is unsurprising since the three-body effects are suppressed by an additional volume factor relative to two-body interactions.
In terms of fit quality, we note that for both masses, the fits for the linear form with varying LECs produce a χ 2 per degree of freedom around 2. This indicates that there is a slight tension between the data and the parametrization used here. At this point it is not easy to determine whether this is the result of the quality of the lattice data or due to the lack of flexibility in the fit form used for the three-body terms. We note that since the three-body predictions are sensitive to the two-body inputs, some of this tension might have as a source small discrepancies in the two-body amplitudes used in the quantization conditions. We note that when analyzing two-meson energy levels using mIAM framework, a similar level of agreement between data and predictions was found [88].
An interesting question is whether we can extract the LECs parametrizing the two-body interactions directly from three-body energy levels. The three-body energy data set does not provide enough constraints to pin down both the LECs and the three-body terms. We are however able to fit the LECs when the three-body terms are set to zero. Settingt 0,1 to zero makes sense, since their effect is rather small. To stabilize the root finding routines used to predict the three-body energy levels as a function of the LECs, we constrained the region scanned for the LECs to a reasonable window, within one order of magnitude of the values determined from the two-body fits. Procedurally this was accomplished using a set of relaxed priors. We used a correlation matrix Σ relaxed = 30 2 × Σ, so that the equivalent error bands on the LECs were at the level of 100%, in effect constraining only the order of magnitude of the LECs.
The results for these fits are included in the last rows of Table III. We find that the values of the LECs are close to the ones generated from the two-body fits, albeit with larger error bars. This provides a good cross-check for the formalism and suggests that with enough three-body energy levels, we should be able to also constrain the two-body amplitudes.
To put the results on the three-body force in perspective, we compare our determination of the three-body term with those obtained in the literature [43,56] in Fig. 6. In doing so, the matching of corresponding three-body terms can be made on the level of scattering amplitudes applying the procedure discussed in Sec. III C. We note that this yields an approximate identification K iso,0 df,3 6(t 0 + 9t 1 ) and K iso,1 df,3 54t 1 . We see reasonable agreement between different collaborations, not too different from the leading order ChPT prediction. This indicates the rapid progress made in the community in mapping out the three-body force.

V. CONCLUSIONS
The field of three-body physics is rapidly advancing, fueled by progress on two fronts. On the one hand, precise energy levels are being produced in LQCD for interacting systems such as three pions or kaons. On the other hand, formalisms that connect the finite-volume QCD spectrum and infinite-volume three-body scattering amplitude, called quantization conditions, are reaching maturity. Such progress has allowed the possibility of extracting quantitative information on the three-body force from first principles.
In this work we apply the FVU formalism to analyze the spectrum obtained previously in Ref. [55]. We used a minimal parametrization for the three-body contact term and constrain the parameters from fits to the spectrum extracted using lattice QCD. We find that the heavy quark mass results are compatible with expectations from leading order ChPT, but our lower mass results are in tension with the predictions. Note that this is similar to other LQCD determinations of this term in the RFT framework [43,56]. The effects of the three-body force terms are small, in broad agreement with other lattice QCD extractions. We also perform a fit to the three-body spectrum to constrain the parameters of the two-body amplitude. We find that the results are in agreement with the values extracted from the two-body spectrum, indicating that the two-body amplitude can also be determined consistently from the three-body data. While this is expected theoretically, it is an important feasibility check.
This study serves as a benchmark for the fitting strategy and the tools developed to generate predictions for the finite-volume three-body spectrum at maximal isospin. In this channel the effects of the three-body force are small, so to constrain it better we need more energy levels and/or better precision for the lattice data. This is likely to be done in the near future. Another direction that is being explored is to study other three-body channels where resonant amplitudes contribute [93]. To this end, both lattice QCD data needs to be generated and the FVU quantization condition must be extended.   where k * ≡ k * (k(k), l(l)), with the dependence on spectator momentum l explicitly denoted. The boost momenta areP from lattice to three-body rest frame, leading to the JacobianJ, and −l from three-body rest frame to twobody rest frame, leading to the Jacobian J in Eq. (A3).
To further discuss the kinematics, we consider the incoming and outgoing pion momentaq j andp j , j = 1, 2, 3. The 3π + system has then momentumP =q 1 +q 2 +q 3 = p 1 +p 2 +p 3 . Momenta in the three-body rest frame are [94] l =l + EP √ s − 1 l ·P and analogously for the other momenta k, p and q. In Eq. (A4), EP = s +P 2 . The boost of Eq. (A4) leads to the JacobianJ The isobar is not at rest in the three-body rest frame. Thus, an additional boost (by −l) has to be performed for the pertinent summation of momenta k * in the self energy of the isobar of Eq. (A3). This is detailed in Eqs. (11, 12) of Ref. [16] and reads in the current notation leading to the Jacobian i.e. the symmetry normalization is not contained in the definition of T 22,∞ . Furthermore, note that v = 2ṽ where the factor of 2 accounts for the possibilities to connect a decay vertexṽ to two external pions [76]. In other words, v contains the same multiplicity as a Feynman vertex would carry from the two possible contractions with external fields.
This expression may be used instead of Eq. (A3) as it converges rapidly and provides exact matching to the ChPT loopJ ππ .
The difference between Σ (2) P Lη and ΣP Lη from Eq. (A3) is exponentially suppressed for all σ. Indeed, whenever a finite-volume pole arises at σ − σ = 0 (remember that σ is discrete in finite volume), we have (σ/σ ) 2 = 1, i.e., the finite-volume poles in both expressions have the same residues. Furthermore, the difference is of the form as an inspection of the arguments of the summations shows: The exponentially suppressed difference can be sizable and becomes, in general, larger for more subtractions in the dispersion relation; if one aims at regaining a particular form of the self energy, one can reduce the difference.
To do this we solve explicitly for A 1 , with the crucial advantage that it has to be calculated only once for given L, η, three-body boostP , and two-body boost by −l. The quantity A 1 can then be recycled, e.g., in fitting, where the entire σ-dependence of Σ needs to be known. In other words, one can still take advantage of the speedup provided by Σ (2) P Lη at the cost of having to determine A 1 , or A 0 and A 1 , once for every L, η,P , −l. The quantity