Finite-volume energy spectrum of the $K^-K^-K^-$ system

The dynamics of multi-kaon systems are of relevance for several areas of nuclear physics. However, even the simplest systems, two and three kaons, are hard to prepare and study experimentally. Here we show how to extract this information using first-principle lattice QCD results. We (1) extend the relativistic three-body quantization condition to the strangeness sector, predicting for the first time the excited level finite-volume spectrum of three kaon systems at maximal isospin, and (2) present a first lattice QCD calculation of the excited levels of this system in a finite box. We compare our predictions with the lattice results reported here and with previous ground state calculations and find very good agreement.


INTRODUCTION
In recent years hadron-hadron scattering information from first-principles lattice QCD (LQCD) has become significantly more accessible. Fostered by advances in theoretical and computational tools a large number of high-precision studies have been performed in the meson sector, see for example Refs. , and Ref. [31] for a review. Several research areas of nuclear physics benefit from these studies. For instance, the study of pion, kaon, and proton correlations in heavy ion collisions by the ALICE@CERN collaboration [32] relies on the value of the K − K − scattering length determined in a lattice calculation [33]. There are however only few results focusing on the strange sector [23,33,34], in contrast to the pion systems explored extensively by a number of collaborations [6,15,20,28,[34][35][36][37][38][39][40]. Furthermore, information about many-K − systems is relevant for the understanding of strange nuclear matter and its implications to the equation of state of neutron stars. In particular, it is well known that ultra dense environments (such as those in the core of neutron stars) allow for an appearance of kaon condensates [41][42][43][44], that can soften the equation of state of neutron stars [42,43,45,46]. Further details on the antikaon interaction with baryonic matter can be found in reviews [47,48].
Today, the frontier of hadronic scattering in LQCD is in the scattering of three mesons. Pioneering lattice calculations have moved from the extraction of the ground states of such systems [49][50][51] to the high-precision determination of multiple excited three-hadron states [52][53][54][55][56]. Significant progress has also been made in the development of formalisms relating the finite-and infinite-volume three-hadron spectrum . Applications of such approaches to LQCD data have thus far been for three pion systems in maximal isospin [54-56, 58, 62, 63].
In this paper, we extend these methods to explore a new area: we present both the first determination of the excited three-kaon finite-volume spectrum from LQCD, along with the first connection to infinite-volume scattering using the formalism of Refs. [29,58,63,68,72]. The latter is extended to the three-flavor sector allowing for chiral extrapolations along arbitrary M K (M π ) trajectories using constraints from chiral symmetry. Such implementations are standard in the two-body sector [98][99][100][101][102][103][104][105], but not yet explored for the three-body systems. The present study closes this gap, using relativistic threebody formalism implementing two-body input from the inverse amplitude approach [106,107].

FINITE-VOLUME SPECTRUM FROM LATTICE QCD
The finite-volume spectrum of hadronic states can be directly accessed by studying correlation functions in the framework of LQCD. Here we review the procedure for extracting the finite-volume spectrum of K − K − K − . The energy levels of hadrons in a finite volume can be extracted from the large time behavior of correlation functions consisting of interpolating operators, O i , which create/annihilate the hadrons of interest, (1) If the operators are constructed to overlap with the states n of interest, we can extract the finite-volume energies E n . An important tool to allow the extraction of multiple finite-volume energies is to perform a variational analysis on a matrix of correlation functions constructed from several operators. This is equivalent to solving a generalized eigenvalue problem [108][109][110], and extracting the finite-volume spectrum from the eigenvalues of the correlation matrix. Due to the precision with which we can measure the correlation functions, thermal effects due to the finite temporal extent must be accounted for as in Ref. [54]. The overlap factor n |O j | 0 is non-zero only if our operators and states n have the same quantum numbers. In a finite cubic volume, the rotational symmetry group is reduced from SO(3) to O h . We have to therefore construct our operators with definite quantum numbers according to the irreducible representations (irreps) of O h . An important consequence is that the irreps of O h mix different angular momentum from the infinite volume. The symmetry is further reduced if the system is studied with non-zero total momentum.
To create operators which overlap with the three-kaon spectrum, we begin by constructing a single kaon interpolator according to where s, u are the quark fields, and the momentum matrix Γ(p) = e ip·x γ 5 projects the operator to definite momentum. Our three-kaon operators are now just a product of three single kaon operators. We project the three kaon operators to irreps of the cubic group. To project to row λ of irrep Λ of group G, we evaluate where p 1 , p 2 , p 3 are the three-momenta of each kaon, R is the three-dimensional rotation matrix associated with g, and U is the representation matrix of g in irrep Λ. The GWUQCD ensembles are generated using two mass-degenerate light quarks (N f = 2 QCD), using the nHYP-smeared clover action. Lattice parameters of the ensembles used here are listed in Table I [33] with statistical (red) and systematic (gray) error bars. The chiral prediction at the physical point is indicated (green star), as well as predictions for different pion masses using fπ, fK from NPLQCD [115] (solid blue line), or NLO extrapolating fπ (blue dash-dotted line).
of the ensemble generation, including some discussion on tuning the bare strange quark mass, can be found in Refs. [18,111]. The pion and kaon decay constants, f π , f K are determined using the procedure outlined in Ref. [112]. The strange quark mass is tuned by setting the ratio R = (M K /M π ) 2 to its physical value. For the valence quarks appearing in (kaon) interpolating operators, both light and strange (all-to-all) quark propagation is treated using the LapH method [113]. The all-to-all LapH propagators were computed using our optimized inverters [114]. The lattice results and predictions are tabulated in the Supplementary Material. Jackknife samples are provided as ancillary files with the arXiv submission.

FINITE-VOLUME SPECTRUM FROM INFINITE-VOLUME PHYSICS
In the present work we utilize the three-body relativistic quantization condition (3bQC) derived in Ref. [72] extended later to higher irreps [68], elongations [54], and boosts [58]. For the I = 3/2, S = −3 three-meson channel the 3bQC reads where E 3 and P denote the center of mass energy and total three-momentum of the three-body system, respectively. Note that the implicit dependence on the latter is suppressed. The determinant is taken with respect to the in/outgoing discrete lattice spectator momenta p/q after projecting the elements in parenthesis to an irrep Γ. The non-diagonal matrix B denotes the one-particle exchange term, while the diagonal matrix ρ L (E 3 , P ) represents the two-body self-energy term, see, e.g., appendix of Ref. [58] for explicit expressions. The propagation of the spectator yields the factor [E L ] pq = δ pq 2L 3 M 2 π + p 2 . The only unknown pieces of the quantization condition are matrices K −1 2 (E 3 ) and C(E 3 ), encoding dynamics of two-(via the usual K-matrix) and three-body interactions, respectively. Since not many data is available yet for the 3K − -system and in analogy to the similar 3π + system [63], we set the latter to zero. The two-body Kmatrix is restricted to the dominant S-wave, noting that due to the nature of the 3bQC all relative partial waves between the spectator and the two-body subsystem are included automatically by the one-particle exchange term B. Specifically, the K-matrix is chosen to match the inverse amplitude method [102,106,107,116,117] -a very successful description of two-meson scattering across wide energy and meson mass ranges and all two-pseudoscalar meson interaction channels [29], Here, T (N)LO refers to the (next-to-)leading chiral order scattering amplitudes [118], and ρ(s) denotes the usual finite part of the two-body self-energy evaluated in dimensional regularization. The K − K − amplitude to one loop is obtained by using crossing symmetry for results of Ref. [106]. A summary of the relevant formulas is included in the Supplementary Material. In particular, the corresponding K-matrix depends on {M π , M K , f π , f K } as well as renormalized low-energy constants (LECs) {L r i }. The effect of the first set of parameters is more important than the latter for not too large meson masses, because the chiral series is ordered in powers of M 2 /f 2 . Thus, we fix the LECs to the results of the most recent global fits to the lattice results [100] (discussion of older LECs is moved to the Supplementary Material), but explore various scenarios for the remaining inputs below.
As a check we evaluate the scattering length M K a 0 = T 2 (4M 2 K ) at different meson masses and compare with the NPLQCD collaboration results [33] along their (M π , M K ) trajectory (M π , M K ) ∈ {(293, 583), (355, 601), (493, 643), (592, 680)} MeV. For the decay constants we compare two scenarios: (S1) by extrapolating the pion decay constant using input at the physical point and NLO chiral expressions [118] with LECs from Ref. [100] and (S2) by using the meson decay The MK (Mπ) trajectory extends linearly from the physical point through the two shown lattice points. The insert in the A1u plot shows the ground state data, predictions, and a prediction for the physical point. Note that the excited states in A1u and Eu are close to or beyond the πKKK threshold, but below the lowest relevant lattice threshold, for which two kaons necessarily have finite back-to-back momenta due to parity conservation.
constants determined on the lattice [115]. These two scenarios differ only by higher chiral orders and are employed as representatives of the systematic uncertainty of our predictions. The results are depicted in Fig. 1. They show that the three-flavor formulation of the inverse amplitude approach (5) is a perfectly suitable parametrization of the two-body dynamics at unphysical meson masses. Higher order terms yield sizable corrections above M π ≈ 4M phys π as expected.
Before coming to the results on three-body spectra we point out the major difference between the 3bQC and its two-body equivalent. The 3bQC remains a determinant equation even for the simplest one-channel case. For a fixed energy and momentum of the three-body system, the two-body input is required for a large kinematic range (s in Eq. (5)) due to the variable spectator momentum. Therefore, the two-body amplitude is often evaluated for subthreshold values of s. Various approaches to this issue have been studied in the past [58,72] and it was found that the obtained finite-volume spectra depend little on the subthreshold region. For the present case, we confirm this observation explicitly by varying the cut in the spectator momentum space in Eq. (4). As we change this from the value used throughout this study, L|p max | = 2π √ 5, to L|p max | = 2π √ 11, the largest change (∼ 3×10 −3 %) among all levels in the GWUQCD setup for M π = 315 MeV happens for the first excited level in A 1u . Similarly, we study the dependence on the subthreshold KK amplitude by replacing K −1 2 (s) with a real-valued constant at s = 3M 2 K and then at s = 3.95M 2 K leading to a maximal change of any energy eigenvalue of 0.02% which is a fraction of the smallest statistical uncertainty in the GWUQCD lattice data. The dependence of the results on the use of modified IAM (mIAM) [119] instead of IAM leads to 0.1% change of the scattering length and is, therefore, of similar size for the three-body ground state energy shift. In summary, these sources of systematic uncertainty are very small.

COMPARISON AND DISCUSSION
First, we turn to previous LQCD results, namely the ground state A 1u (0) levels determined by the NPLQCD collaboration [51,115] in a cubic box of L = 2.5 fm at four pairs of pseudo-scalar meson masses. These results are depicted in the top left panel of Fig. 2 and overlayed by the predictions of the 3bQC. As before, differences between scenarios S1 and S2 are visualized by the light blue band. We observe encouraging agreement between our predictions and the NPLQCD results. We also find a similar increase in the size of NNLO effects at higher pion masses, as observed in the two-body results. Going beyond the ground state level, we extend our predictions to excited states, other irreps, and boosts (lower panel of Fig. 2).
New LQCD results are obtained in this work, including for the first time, excited K − K − K − energies, in multiple irreps. This provides an opportunity for a series of unique tests of the predicted finite-volume spectra. Following both scenarios S1 and S2 along the chiral trajectory (see Table I) the predictions for the GWUQCD setup are shown in Fig. 3. The ground state is in excellent agreement with the predictions as was the case for the heavier pion mass results from NPLQCD. For the excited A 1u (0) levels the slight tension with the prediction could be some hint of the need for a non-zero three-body force. Of course, other possible sources for the discrepancy could be (i) the chiral prediction itself is not perfect, or (ii) that the partial quenching of the strange quark plays a role. This will be investigated in a future study.
The E u (0) levels agree with the predictions well. Note that this irrep is dominated by D-wave. Since the twobody interaction is typically smaller for higher partial waves, the major contribution seems to come from the one-particle exchange term B, with no obvious need for contact terms beyond that. In fact, this is very similar to the observed E − u (0)/A − 1u (0) pattern for the three-pion system noted in Ref. [58]. In both cases, the pattern confirms the dominance of the exchange contribution, which is a direct consequence of the S-matrix principle of three-body unitarity.
In summary, we have traced a pathway for studying multi-kaon systems using lattice QCD. We presented the first LQCD calculation of excited three kaon states, in multiple irreps, and at multiple pion masses. We have also extended the relativistic three-body quantization condition to the strange sector, allowing for chiral extrapolations along arbitrary trajectories. We find that this extension consistently describes the data from two independent lattice calculations of multi-kaon systems. In the long run, this provides an avenue for extracting information relevant for strange resonances, kaon condensates, and heavy-ion collisions. In the following we provide the explicit formulas for the scattering matrices of the two-body input. They rely on the three-flavor chiral Lagrangian of the leading and next-to-leading chiral order [118]. We chose the formulation with all decay constants replaced by the "physical" pion one, which is done consistently at the forth chiral order, i.e., differences are of the order O(p 6 ). In practice, we obtain the K − K − → K − K − amplitude from the K + K − → K + K − amplitude quoted in Ref. [106] (see also Ref. [116] for the original calculation) by crossing symmetry which amounts to exchanging s ↔ u in the latter, Here s, t, u denote the Mandelstam variables and ellipses denote the higher chiral orders not taken into account. The tadpole integrals arising from, e.g., the wave function re-normalization procedure are denoted by µ i and read where µ is the renormalization scale and f 0 is the meson decay constant in the chiral limit. The latter is determined dynamically, from employed sets of LECs and {M π , M K , f π , f K } (see discussion below). The finite parts of the meson-meson loop integralsJ..(s) are given bȳ and resembles the self energy part ρ(s) = 16πJ KK (s) used in the main part of the manuscript.
For the purpose of the present work, these formulas are projected to the S-wave. This implies a factor N = 2 for the two identical kaons [106]. The amplitude is subsequently back-transformed to the (two-body) plane-wave basis as required by the form of the quantization condition in Eq. (4).
There are three mass relations and three relations for the decay constants (see, e.g., Refs. [106,118]) setting these physical quantities in relation to the tree level quantities f 0 , M 0π , M 0K , and M 0η at NLO (we also replace M 0η using the Gell-Mann-Okubo relation). For the chiral predictions, we follow two strategies, S1: The tree level quantities are determined at the physical point from M π , M K , M η , f π , f K , f η . Then, f π is obtained by replacing tree level masses by the meson masses on the lattice in the relation For strategy S2: The tree level quantity f 0 is directly inferred from the available lattice information, in the present cases (GWUQCD and NPLQCD) M π , M K , f π , and f K . We consider this method more reliable as one source of extrapolation uncertainty is removed. We replace also the η-mass in Eq. (6) using the Gell-Mann-Okubo relation. All discussed replacements only lead to O(p 6 ) effects and are, thus, consistent to the chosen chiral order. Note also that the scattering amplitudes are regularization (µ) independent, which implies that the LECs are only defined at the given scale for which we choose the same value of µ = 770 MeV as in Ref. [100] to be able to use their LECs for predictions. Note that we cannot use the LECs of Ref. [102] because they correspond to chiral amplitudes formulated in terms of f π , f K , f η as explained in Ref. [117]. We emphasize that effects due to change of the renormalization scale have been studied thoroughly in the two-flavour case in Ref. [29], where they have been found subdominant to the other effects, e.g., due to pseudo-scalar masses, decay constants and LECs.
While the sets of available SU(2) LECs produce very consistent predictions of the three-pion spectrum [54,58], the eight SU(3) LECs are less well determined; we chose the LECs of Ref. [100] (fit 4) for this study because the pertinent fit includes lattice data; if one uses the older LECs of Ref. [106] obtained from only fits to experimental data, the scattering length turns out to be about 25% smaller than the one shown in Fig. 1, in contradiction with the NPLQCD results [115]; likewise, the three-body ground state energy shift gets about 25% smaller, in contradiction with the lattice data.

Lattice energies
Here we tabulate the lattice energy levels extracted from the ensembles in Table I in the main text. For the 315 MeV ensemble (E 1 ), triple exponential fits were performed as in Ref. [54]. Due to an increase in the noise on the 220 MeV ensemble (E 4 ), only single exponential fits were performed. The energies are tabulated in Table II. The final column shows the energies predicted from SU(3) IAM using lattice decay constants and masses from Table I