Weakly bound $H$ dibaryon from SU(3)-flavor-symmetric QCD

We present the first study of baryon-baryon interactions in the continuum limit of lattice QCD, finding unexpectedly large lattice artifacts. Specifically, we determine the binding energy of the $H$ dibaryon at a single quark-mass point. The calculation is performed at six values of the lattice spacing $a$, using O($a$)-improved Wilson fermions at the SU(3)-symmetric point with $m_\pi=m_K\approx 420$ MeV. Energy levels are extracted by applying a variational method to correlation matrices of bilocal two-baryon interpolating operators computed using the distillation technique. Our analysis employs L\"uscher's finite-volume quantization condition to determine the scattering phase shifts from the spectrum and vice versa, both above and below the two-baryon threshold. We perform global fits to the lattice spectra using parametrizations of the phase shift, supplemented by terms describing discretization effects, then extrapolate the lattice spacing to zero. The phase shift and the binding energy determined from it are found to be strongly affected by lattice artifacts. Our estimate of the binding energy in the continuum limit of three-flavor QCD is $B_H^{\text{SU(3)}_{\rm f}}=4.56\pm1.13_{\rm stat}\pm0.63_{\rm syst}$ MeV.

The H dibaryon is a scalar six-quark state with flavor content uuddss, originally proposed in 1977 by Jaffe [1]. Despite years of effort, experimental searches have not produced any hard evidence for its existence [2-4]. However, an upper bound on its binding energy has been derived from the observed production and decay pattern of a doubly strange 6 ΛΛ He hypernucleus [2, 3]. Studying the properties of a potential Λ-Λ bound state will help our understanding of the hadronic (Λ-Λ) interaction, which is relevant for the physics of double hypernuclei, neutron-rich matter and neutron stars. Recently, experimental data for two-particle correlations in p-p, p-Pb and Au-Au collisions [5][6][7] have been analyzed to constrain the Λ-Λ interaction and provide model estimates for the binding energy of the H dibaryon. In addition, a dedicated experiment is planned to search for it at J-PARC [8]. Other approaches to study the H dibaryon include chiral effective field theory [9][10][11][12] and lattice QCD.
Lattice QCD studies of dibaryons and baryon-baryon scattering are very challenging because of the signal-tonoise problem [13,14] and the complexity of contractions. In response to an inconsistency between results in the nucleon-nucleon sector [15,16], there has been a recent focus on improved baryon-baryon spectroscopy methods [17][18][19]. This work goes beyond that to achieve control over all systematic effects for the H-dibaryon channel at one unphysical quark mass point.
In our previous work [17], using gauge fields with dynamical u and d quarks and a quenched s quark, we found that the distillation method [35] produced a better determination of the two-baryon spectrum than previously used methods. At a heavy SU(3)-symmetric point with a pion mass of 960 MeV, we obtained B H = 19 ± 10 MeV.
In this letter we extend our calculations to lattice QCD with dynamical u, d, and s quarks with degenerate masses set to their physical average value, corresponding to m π = m K ≈ 420 MeV [36]. We present the first sys- Points are from fits to individual ensembles and curves are from the combined fits to the spectra of different subsets of the ensembles; they are not fitted to these points. Gray diamonds show results from the small-volume ensembles and the black cross shows our final estimate.
tematic study of discretization effects in a multibaryon system, by computing finite-volume spectra at several lattice spacings, extrapolating the corresponding scattering phase shift to the continuum limit, and determining the binding energy. As shown in Fig. 1, at vanishing lattice spacing, we find B SU(3) f H = 4.56 ± 1.30 MeV, which is smaller than the result at the coarsest lattice spacing by a factor of about 7.5. We conclude that a thorough investigation of lattice artifacts is indispensable for answering the question whether a bound H dibaryon exists in nature.
Our calculations are based on a set of eight gauge ensembles generated by CLS [37], with a nonperturbatively O(a)-improved Wilson-clover fermion action. These ensembles have six different values of the lattice spacing and multiple box sizes L (all satisfying m π L ≥ 4.4) as shown in the inset of Fig. 3 [38].
For each ensemble, we determine the energy levels in the rest frame and in four moving frames. To this end, in each frame we compute a Hermitian matrix of two-point correlation functions from a basis of interpolating operators: C ij (t) ≡ O i (t)O † j (0) . The finite-volume spectrum {E n } determines the exponential fall-off of C ij (t).
The building blocks of our operator basis are products of two single-baryon operators projected to momenta p 1 and p 2 with total spin zero or one. For each frame momentum P = p 1 + p 2 , we take linear combinations that transform under the trivial irreducible representation of the little group of P , which contains the 1 S 0 scattering channel [38]. Following Refs. [31,39], the flavor content of our interpolating operators is a linear combination of isospin-zero ΛΛ, ΣΣ, and symmetric N Ξ that corresponds to the singlet irreducible representation of SU(3)-flavor.
Calculating the correlation functions of bilocal operators requires the ability to compute "timeslice-to-all" quark propagators. As in our previous study [17], we have used the distillation technique [35,38].
The finite-volume energy levels in each frame are determined by solving a generalized eigenvalue problem (GEVP) [38,40,41], C(τ D )v n = λ n C(τ 0 )v n , for fixed τ 0 and τ D satisfying τ D > τ 0 ≥ τ D /2. We then use the eigenvectors v n to constructC nm (t) ≡ v † n C(t)v m , an approximately diagonalized correlator matrix. We have verified that different combinations of (τ 0 , τ D ) yield consistent results across a wide range of values [38].
Before fitting to the data, we divide the rotated twobaryon correlators by a product of two single-baryon correlators that form the corresponding two-baryon nonin- is a single-Λ correlator with momentum p i , and the total frame momentum is p 1 + p 2 . The leading term in this ratio falls off exponentially with the shift ∆E of the interacting two-baryon energy away from the noninteracting level. In the ratio, we observe a partial cancellation of correlated statistical fluctuations and residual contributions from excited states, which helps in the reliable determination of ∆E. Our finite-volume energies are determined from singleexponential fits to R n (t). For all levels, we choose t min , i.e. the smallest time separation included in the fits, to lie in the plateau region of R n (t). We also aim to have t min lie in the plateau region of the single-baryon correlators, and in the majority of cases we set it to be the first time separation in this plateau region. Since the single-baryon correlators take longer than the two-baryon correlators to reach their asymptotic behavior, this ensures that all correlators entering the ratio have little to no excitedstate contamination. In some cases, however, the signal of R n (t) is already significantly degraded at the start of the single-baryon plateau region, and we are led to choose a slightly lower t min that still lies within the plateau region of the correlator ratio. For all levels, we estimate the sensitivity to t min by extracting an alternative spectrum with t min further lowered, and use it in subsequent analyses to estimate the systematic uncertainty of our energy determination.
The fits also yield the couplings between each energy eigenstate and our operators. For each frame that includes a spin-one operator, we find one eigenstate that has strong overlap with only that operator, allowing for a simple identification of the spin-one dominated states. Figure 2 shows the effective energy difference ∆E eff (t) ≡ − d dt log R(t) and the extracted ∆E for the ground state in frame (0, 0, 1) on four ensembles that differ primarily in their lattice spacing. This level is particularly important because it is the closest to the bound-state pole determined in the phase shift analysis. An overview of the finite-volume spectrum is shown in Fig. 3 ing is reduced, indicating that discretization effects are significant. Given the two-particle scattering amplitude, Lüscher's finite-volume quantization condition [42] and its generalizations [43][44][45] determine the finite-volume spectrum, up to exponentially suppressed corrections, between the tchannel cut (p 2 > −m 2 π /4) and the three-particle threshold (E cm < 2m B + m π ). Since the quantization condition is diagonal in spin [44,45], the spin-one part of the scattering amplitude does not affect the spin-zero finitevolume spectrum, and we choose to ignore the spin-one states. In addition, we neglect higher partial waves starting from 1 D 2 . In this case, the quantization condition yields the 1 S 0 phase shift δ(p) at the momentum corresponding to each finite-volume energy level: where γ = E/E cm and Z D 00 is a generalized zeta function. In addition to excluding levels with too-low or too-high p 2 from our analysis, we must also exclude the first excited levels in frames (0, 1, 1) and (1, 1, 1), as the 1 D 2 partial wave is necessary to describe their position below the lowest noninteracting level [38].
The quantization conditions do not take discretization effects into account; strictly speaking, they are only valid in the continuum. There is no general formalism for finite-volume quantization at nonzero lattice spacing, except for a simple model studied in Ref. [46]. In principle, discretization effects would affect both the scattering amplitude and the finite-volume quantization condition. Effects on the former could include a-dependence and frame-dependence of the scattering amplitude, as well as couplings between J P that are forbidden in the continuum. Effects on the latter could include a modification of the zeta functions [46]. Either way, discretization effects might spoil the factorization that separates spinzero from spin-one. Lacking a rigorous understanding, we have elected to model discretization effects in a simple way, by allowing the parameters of the phase shift to depend on a.
Our primary analysis is based on combined fits of the dependence of the phase shift on both p 2 and a. Specifically, our model is (2) Concerning the dependence on p 2 , we fit in two ways. The first uses the near-threshold region, |p 2 | < ∼ m 2 π /4 (where the effective range expansion converges), with N = 2 terms for the dependence on p 2 . The second uses the full p 2 range, starting from the t-channel cut and stopping just below the three-particle threshold, with N = 3. Given {c ij }, solving Eq. (1) yields a discrete spectrum of p 2 for each volume and frame; we fit these to the lattice spectra. For comparison, we also performed fits to individual ensembles, neglecting discretization effects. Given δ(p), a solution below threshold to p cot δ(p) = − −p 2 corresponds to a bound state pole. All of the fits yielded a bound H dibaryon.
Our preferred fit is to all ensembles using the full p 2 range; the corresponding continuum interacting energy levels are shown as blue curves in Fig. 3. In addition to the alternative spectrum fit range, we estimate the systematic uncertainty using the root-mean-square difference of alternative combined fits that cover all combinations of cuts on p 2 (full range or near threshold), a (all six or the finest four), and L (excluding L ≈ 2.1 fm or not). All of these fits have acceptable fit quality, with p-values between 0.2 and 0.9. We explored adding an a 3 term in Eq. (2) but found that this reduces χ 2 by at most 1.1 for each additional fit parameter, a sign of overfitting.
The phase shifts from the preferred fit, in the continuum and at nonzero lattice spacing corresponding to the four ensembles with L ≈ 2.4 fm (J500, N300, B450, A653), are shown in Fig. 4. Since these ensembles have similar values of L, they allow us to perform a cross check, shown in the lower panel. We select the volume of ensemble B450 as our target and call this box size L * . For the three other lattice spacings, we estimate each energy level at L * by shifting from L using the quantization condition: p 2 (L * ) ≈ p 2 (L) + p 2 q.c. (L * ) − p 2 q.c. (L). For each energy level, we then study the dependence of p 2 (L * ) on the lattice spacing and compare it with the value obtained from applying the quantization condition to the continuum limit of the preferred fit. The cross check shows that a level-by-level continuum extrapolation at L * is consistent with the latter. However, some levels show curvature in the dependence on a 2 and the fixed-L * extrapolation is less precise, making it less useful than the combined fits. Near threshold, we can write p cot δ = −1/a 0 +r 0 p 2 /2+ O(p 4 ), where a 0 is the scattering length and r 0 is the effective range. We obtain where the first error is statistical and the second is systematic. The dependence of the H dibaryon binding energy on a is shown in Fig. 1; in the continuum, we obtain which is substantially lower than the binding energies determined at nonzero lattice spacing, except on the finest two of our ensembles. We have reported the first lattice study of a baryonbaryon system in the continuum limit. The crucial elements of our methodology are the finite-volume quantization condition, supplemented by terms describing discretization effects and applied over a wide range of lattice spacings, as well as the subsequent extrapolation to the continuum limit. We conclude that cutoff effects are large and cannot be ignored in an investigation of the H dibaryon using lattice QCD; it will be essential to study their importance in other multibaryon systems such as the deuteron, where calculations disagree [15,16,18].
Our final result for the binding energy, given in Eq. (5), suggests the existence of a weakly bound H dibaryon, which is not only at variance with Jaffe's original bag model prediction [1] of a deeply bound uuddss state, but is also substantially lower than the binding energies determined in previous lattice calculations [17,[27][28][29][30][31][32][33] at nonzero lattice spacing (see Fig. 5). This adds to the evidence against deeply bound hexaquark dark matter [47][48][49][50][51][52][53][54]. An obvious caveat is that our calculation was performed for one set of degenerate quark masses. The issue of SU(3) symmetry breaking -which is crucial, since the splitting between physical ΛΛ and N Ξ thresholds is larger than B [55]. Previous estimates based on extrapolations of lattice data found a bound state at the physical point unlikely [9,11,33,56,57]; our smaller binding energy should make it even less likely.
We thank Maxwell T. Hansen, Ben Hörz, and Daniel Mohler for many helpful conversations. Calculations for this project used resources on the supercomputers JUQUEEN [58], JURECA [59], and JUWELS [60] at Jülich Supercomputing Centre (JSC). The authors gratefully acknowledge the support of the John von Neumann Institute for Computing and Gauss Centre for Supercomputing e.V. (http://www.gauss-centre.eu) for project HMZ21. The raw distillation data were computed using QDP++ [61], PRIMME [62], and the de- Pale points (displaced vertically) show the levels before adjustment. The spectrum obtained from the continuum phase shift is indicated using blue crosses. Curves show continuum extrapolations of the form b0 + b1a 2 excluding the coarsest lattice spacing (solid magenta) and b0 + b1a 2 + b2a 3 using all four lattice spacings (dashed cyan).
flated SAP+GCR solver from openQCD [63]. Contractions were performed with a high-performance BLAS library using the Python package opt einsum [64]. The correlator analysis was done using SigMonD [65]. Much of the data handling and the subsequent phase shift analysis was done using NumPy [66] and SciPy [67]. The plots were prepared using Matplotlib [68]. This

Supplemental material
In this supplement, we provide additional details for our calculation. Section I specifies the lattice action and ensembles. The precise definitions of our interpolating operators are given in Section II. Our implementation of the distillation approach is described in Section III. We provide further details about our determination of the spectrum in Section IV and our combined fits to the spectra at different lattice spacings in Section V. The analysis of two N f = 2 ensembles is provided in Section VI. Finally, Section VII describes the spectrum data being made available with this article.

I. LATTICE ENSEMBLES
Our calculations are based on a set of gauge ensembles with N f = 2 + 1 flavors of dynamical quarks, generated by CLS using the openQCD code suite [70] and listed in Table SI. The fields are described by the tree-level O(a 2 )improved Lüscher-Weisz action and the O(a)-improved Wilson-Clover action in the quark sector, with the improvement coefficient c sw tuned to the nonperturbative determination of Ref. [71]. Open or periodic boundary conditions in the time direction are employed. All ensembles realize SU(3) symmetry, with m π = m K ≈ 420 MeV, at six different values of the lattice spacing, covering a range between 0.04 and 0.1 fm. Here we also take the opportunity to extend our earlier calculations with N f = 2 flavors of dynamical quarks [17]. The respective simulation parameters are listed in Table SI, and a detailed description can be found in Section VI.
As discussed in Ref. [72], the quark masses are not exactly matched among the different lattice spacings. Given our choice of scale setting, this corresponds to a 3% variation in the pion mass, from 411 to 424 MeV. This is expected to produce a shift in the octet baryon mass of order 10 MeV, preventing a simple study of discretization effects in the octet baryon mass. However, the latter also varies by just 3% among our ensembles, which puts a likely upper bound on the size of discretization effects. For our main study of baryon-baryon interactions, we always determine energy differences from noninteracting levels and convert them to p 2 using the baryon mass determined on the same ensemble, cancelling the leading effect due to slightly varying baryon masses. Our expectation is that the mistuning of the pion mass will affect the energy differences at the few-percent level, which is much smaller than our statistical uncertainty.

II. INTERPOLATING OPERATORS
In our previous study [17], we found that bilocal twobaryon operators are more effective than local hexaquark operators at identifying the low-lying spectrum; therefore, in this work we use only the former. To begin, we define the single-octet-baryon operators, which make use of the three-quark combination Here r, s, and t denote smeared quark fields of generic flavor at the same point and P (S2) The spin-zero and spin-one two-baryon operators are defined as follows: In these operators, the baryon B j is projected to momentum p j , and the total momentum is P = p 1 + p 2 . Each operator constructed in this way can be identified with a noninteracting finite-volume energy level of energy E = j m 2 Bj + p 2 j . These operators satisfy the exchange symmetry relations This work is focused on flavor-symmetric channels, which implies that the spin-zero operators are even under exchange of momenta and are thus associated with even partial waves, and the opposite is true for the spin-one operators. For each total momentum P , we construct operators that transform under the trivial (A + 1 or A 1 ) irreducible representation of the little group of P , which contains the 1 S 0 scattering channel. Generically, these have the form for some coefficients c j or c ij . For each operator, we choose { p j } such that they lie in the group orbit of a Overview of lattice ensembles. Each ensemble is characterized by the gauge coupling parameter β, the quark hopping parameter κ, the lattice size, and the temporal boundary condition. For N f = 3, the lattice spacing a was determined for the second-finest lattice spacing from the result in Ref. [72] and scaled to the other lattice spacings using the gradient flow scale t0 [73] determined at the symmetric point. For the ensembles with N f = 2 we use the lattice spacing determined in Ref. [74]. The masses of the light octets of pseudoscalar mesons and spin-1/2 baryons are given by mπ and mB, respectively. On each of the N conf gauge configurations analyzed, Ntsrc source timeslices were used. Including both forward and backward-propagating states, the total number of measurements used is Nmeas = NtN conf , where Nt = 2(Ntsrc − N skip ). To avoid boundary effects, we omit some potential measurements, such as the backward-propagating states from the first source timeslice; thus, N skip is 0 for the ensembles with periodic boundary conditions and between 1 and Ntsrc/2 for the ensembles with open boundary conditions. Finally, NLapH is the number of low modes of the Laplacian used in the Laplacian-Heaviside smearing.   (0, 0, 1) + (0, 0, 0) (0, 1, 1) + (0, −1, 0) (0, 1, 1) + (0, −1, 0) (0, 1, 1) A1 (0, 1, 1) + (0, 0, 0) (0, 0, 1) + (0, 1, 0) (0, 0, 1) + (0, 1, 0) (1, 1, 1) A1 (1, 1, 1) + (0, 0, 0) (0, 1, 1) + (1, 0, 0) (0, 1, 1) + (1, 0, 0) (0, 0, 2) A1 (0, 0, 1) + (0, 0, 1) reference momentum p under the little group of P . Representative momenta p 1 and p 2 for each of our operators are listed in Table SII, and these operators are given explicitly in the following subsections. In each frame, we make use of one operator for each noninteracting level below a certain threshold. In the noninteracting and nonrelativistic limit, in all cases the energy gap to the first uncontrolled state, i.e. from the highest level for which an operator is included to the lowest level for which an operator is not included, is (2π/L) 2 /m B , except in frame P = (2π/L)(1, 1, 1), where this gap is doubled.
The flavor content of our chosen operators belongs to the strangeness −2, isospin zero sector: These are transformed to the singlet irreducible representation of flavor SU(3) following Refs. [31,39]: In the following subsections we list the spin-zero and spin-one flavor-symmetric interpolators in the trivial irrep in each frame. Each moving frame has several equivalent copies, related by lattice rotations; the listed operators will be given in a generic way for all equivalent frames, such that all operators in each irrep transform in the same way between equivalent frames. (We have performed a cross-check using computer algebra to verify these transformation properties.) For each term [B 1 B 2 ]( p 1 , p 2 ), only p 1 will be given, since p 2 = P − p 1 . The operators will be labeled [BB] s(n1,n2) Λ, P L/(2π) , where Λ is the irrep, s is the spin, and p 2 i = n i (2π/L) 2 .
A. (0,0,0) A + 1 Here we make use of the standard basis vectors e i . [BB] [BB] [BB] (S15) Here the frame momentum is P = ± 2π L e k for some k.
[BB] 0(0,2) [BB] 0(1,1) [BB] 1(1,1) Note that because we only consider flavor symmetric operators, these are insensitive to the exchange of d 1 and d 2 . [BB] [BB] 0(1,2) As in our previous study [17], we evaluate correlator matrices involving two-baryon operators using the method called distillation [35]. In this approach, the interpolating operators are defined using Laplacian-Heaviside (LapH)-smeared quark fields. LapH smearing uses the N LapH lowest-lying eigenmodes {v (n,t) i ( x) : 1 ≤ n ≤ N LapH } of the spatial gauge-covariant Laplacian (constructed using spatially stout-smeared [75] gauge links) on each timeslice t. The smeared quark fields are obtained by projecting onto the space spanned by these eigenmodes: LapH smearing is a projector onto a much smaller subspace [in practice N LapH N c (L/a) 3 ], making it feasible to compute the full timeslice-to-all quark propagator within this subspace, which is called the perambulator : (S27) The other key object required for evaluating correlation functions involving baryons is the mode triplet, (S28) All of our single-and two-baryon correlation functions can be evaluated by performing tensor contractions of perambulators, mode triplets, and spin matrices. For a fixed choice of timeslices and momentum, the perambulator has size 4N 2 LapH and the mode triplet has size N 3 LapH . (Because of the projector P + in our interpolating operators, there are only two independent spin components.) To keep the smearing width fixed, N LapH should be scaled proportional to the spatial lattice volume, and therefore the scaling of the tensor contraction cost with N LapH should be kept small.
The Wick contractions of quark fields yield two topologically distinct classes of diagrams, shown in Fig. S1. One possible strategy would be, in an intermediate step, to construct two-baryon "source" and "sink" tensors, where the former is the outer product of two mode triplets and the latter additionally includes the six perambulators. This would fully factorize the choice of source and sink operators in the correlator matrix. However, the computational cost would scale with N 6 LapH . Instead, we form partially-contracted source-sink "blocks" (Fig. S2) at a cost proportional to N 4 LapH . Computationally, this is the most costly step in the contractions, and therefore we avoid recomputing blocks that are used in multiple correlators. The cost of combining two blocks to complete a two-baryon contraction is proportional to N 2 LapH and is relatively inexpensive. A similar strategy for two-baryon correlators was described recently in Ref. [76].
In larger volumes, the N 4 LapH cost scaling will eventually become prohibitively expensive. One possible solution is to use stochastic distillation [18,77], which would replace N LapH in the cost scaling with the (much smaller) size of the dilution space.

A. Choosing NLapH
Due to the rise in inversion and contraction costs as N LapH is increased, it is computationally advantageous to use as few LapH eigenvectors as possible. However, making N LapH too small will increase the statistical uncertainty. Hence, for comparison, we computed an octetbaryon correlation function using three values of N LapH on a subset of ensemble U103. The effective energies are shown in the left panel of Fig. S3. It is clearly seen that the error on the effective energy increases as the number of LapH eigenvectors is reduced. At the same time, retaining fewer LapH eigenmodes has resulted in less contamination from the excited states; therefore, a more fair comparison between the three is one in which the onset of the plateau for each effective energy has been shifted to the same point. This is shown in the right panel of Fig. S3, indicating that N LapH = 20 is an acceptable choice. For the other N f = 3 ensembles, N LapH is scaled with the physical three-volume to ensure that the smearing radius remain roughly constant. For the N f = 2 ensembles, we have a single volume and we choose to use a slightly larger N LapH , corresponding to a smaller smearing radius.

IV. ANALYSIS OF CORRELATION FUNCTIONS
The correlation functions computed are of the form where {O i } denotes a set of interpolating operators that all transform irreducibly in the same way, and {t 0 } is the set of N t sources shown in Fig. S4 for all ensembles. The sources and time separations that we include assume t T for periodic boundary conditions, and both 0 t 0 and t + t 0 T for open boundary conditions, such that the effects of the finite temporal extent may be ignored. Under these assumptions, the spectral decomposition of the correlators is given by where |Ω is the vaccuum state, |n are the eigenstates of the system, and E n are the eigenenergies.

A. Octet-baryon mass
In order to calculate p 2 , which is needed for the phaseshift analysis, we must obtain an estimate for the octetbaryon mass. To this end, we perform single-exponential fits to correlators constructed from a single-octet-baryon operator projected to zero momentum. We show the resulting fits and effective energies on four ensembles with similar volumes in Fig. S5.

B. Generalized eigenvalue problem
For all momentum frames that include more than one two-baryon operator, we use the variational approach described in Refs. [40,41], in which a generalized eigenvalue problem (GEVP) is solved from the matrix of correlation functions in Eq. (S29): Provided that τ 0 satisfies τ 0 ≥ t/2, the asymptotic behavior of the generalized eigenvalues is given by [41] where N is the size of the correlator matrix, and the argument τ 0 has been dropped. By contrast, the leading corrections to the eigenvalues of C(t) only fall off as e −t∆n , where ∆ n ≡ min m =n |E n − E m |. Thus, by solving the GEVP rather than the simple eigenvalue problem for C(t), one benefits from a stronger suppression of the contamination from higher excitations.
To simplify the analysis, we turn the GEVP into a normal eigenvalue problem, resulting in the following matrix to be diagonalized and only solve for the eigenvectors and eigenvalues at a single time separation τ D > τ 0 . The resulting eigenvectors can be used to rotateĈ(t) for all other time separations where the columns of V (τ D ) contain the orthonormal eigenvectors ofĈ(τ D ). Then the diagonal elements of C(t) approximate the generalized eigenvalues λ n (t). It can be seen in Fig. S6 that the scattering momenta derived from the spectrum show very little dependence on the chosen GEVP parameters τ 0 and τ D . The rotated correlators are inspected by eye to ensure they remain statistically diagonal for all time separations. Finally, extraction of the leading exponential terms for the diagonal elements ofC(t) gives the lowest N levels that overlap with the states created by the operators used in the correlation matrix, and the overlaps themselves are given by These overlaps are used to identify states as being predominantly spin-zero or spin-one.

C. Ratio fits
In a final step before fitting the correlators, we form a ratio of each diagonal element of the rotated correlator matrix to the product of two single-baryon correlators, . (S36) The momenta p 1,2 are chosen to correspond to the constituent momenta of the individual baryons appearing in the operator that has dominant overlap with state n.
The advantage of forming this ratio is the possibility for partial cancellation of correlations and residual contributions from excited states. One drawback, however, is the loss of the monotonic behavior of the effective energy, making an identification of the plateau less reliable. To avoid this issue, in most cases we fix the lower end of the fit range, t min , on each ensemble to the first time separation in the plateau region of the single-baryon correlators, which were observed to take longer to reach their asymptotic behavior than the two-baryon correlators. However, as mentioned in the main text, in some cases this choice of t min corresponds to a poor signal quality in R n (t) and we instead chose a slightly lower t min . For all of these levels that are also used in the phase shift analysis, the decrease of t min below the start of the single-baryon plateau was by less than 0.12 fm, except on ensemble N300 where the decrease was by 0.25 fm. These choices still lie in the plateau region of R n (t).
To estimate the systematic error corresponding to the chosen fit range, we extracted an alternative spectrum, based on a second value of t min that is below our preferred value by somewhere between 0.0865-0.173 fm, and propagated it through to the subsequent analysis. Finally, the upper end of the fit range, t max , is chosen for each correlator ratio to be one time separation smaller than the first time separation in which |R n (t)| < 3 error(R n (t)). Effective energy differences for two additional ground-state levels are shown in Fig. S7.
By fitting the ratio R n (t), one obtains the shift ∆E of the nth interacting energy eigenstate relative to the corresponding noninteracting level. The interacting energy E is then reconstructed from ∆E by adding the noninteracting energy level, p 2 1 + m 2 B + p 2 2 + m 2 B , determined from the continuum dispersion relation using the singleoctet-baryon energy at rest.

V. COMBINED FITS
We begin by describing the selection of levels that are included in the combined fits. In general, an energy level is excluded for one of four reasons: 1. Energy levels with dominant coupling to spin-one interpolating operators are excluded. The corresponding partial waves such as 3 P 1 factorize in the quantization condition. 2. The spin-zero excited state in frames (0, 1, 1) and (1, 1, 1) cannot be described using the simplest form of the quantization condition. Assuming the phase shift does not pass through zero, Eq.
(1) has a solution between the lowest pair of noninteracting levels, whereas the data lie below this range. Examining these levels in the nonrelativistic limit, one sees that they belong to the same degenerate shell of states, which can contain just one 1 S 0 level. Therefore, higher partial waves are relevant. These levels can be described if 1 D 2 is included in the quantization condition, which we leave to future work [78]; they are excluded from all analyses here.
3. Energy levels with too-high p 2 are susceptible to the influence of the three-particle inelastic threshold, which is described by neither our fit ansatz nor the quantization condition. Therefore we exclude the second excited state in frame (0, 0, 0) on all ensembles except for the two largest volumes, H101 and N202.
4. Energy levels with too-low p 2 are susceptible to the influence of the t-channel cut (arising from the exchange of a pseudoscalar meson), which is described by neither our fit ansatz nor the quantization condition. (We note that the method recently proposed in Ref. [79] might be applicable.) However, on our coarser ensembles the bound-state pole also lies close to the t-channel cut. The ground state in frame (0, 0, 1) is essential for constraining the pole position, and therefore we always include it, even though on our coarsest lattice spacing this level lies below the cut.
On the other hand, for almost all ensembles the ground states in frames (0, 0, 0) and (0, 0, 2) lie below the cut and we exclude these levels. The exception is the largest volume, N202. However, these two levels still lie well below the bound-state pole and are very close to the cut; furthermore, we obtain significantly worse fit quality when either of these two levels is included. (For instance, the single-ensemble fit to the low-p 2 region of N202 has χ 2 /dof = 1.7/1. Including the ground state in the rest frame increases this to 12.0/2. For the fit to the full-p 2 range, χ 2 /dof increases from 4.1/3 to 8.2/4 when including this level.) Therefore, we also exclude these two levels on N202.
Our final choice of levels for the full p 2 range is the following: one or two excited-state levels in frame (0, 0, 0), both the ground and excited spin-zero levels in frame (0, 0, 1), and the ground state in frames (0, 1, 1) and (1, 1, 1). For the near-threshold region, we take the ground state in frame (0, 0, 1) and possibly the ground state in frames (0, 1, 1) and (1, 1, 1). The fits are performed by minimizing with respect to the model parameters, where i indexes all of the levels among all ensembles included in the fit and p 2 q.c.,i is obtained by solving Eq. (1) given the model for p cot δ(p). Here Σ = Σ stat + Σ syst is an estimate of the covariance matrix. Bootstrap resampling is used to obtain Σ stat,ij , which is set to zero when i and j correspond to levels from different ensembles. The alternative spectrum fit range is used to estimate a correlated systematic uncertainty: we set Σ syst,ij = (δp 2 ) i (δp 2 ) j , where δp 2 is the the difference between p 2 obtained using the preferred and alternative spectra.
The statistical uncertainty of our fit results is estimated using bootstrap. When fitting to the nearthreshold region, for a small number of bootstrap resamples (up to 4 out of 1000) the minimum of χ 2 is not a point where its gradient vanishes, but instead lies at a discontinuity. In these rare cases, there exists a level (typically the lowest-lying level in the smallest volume) where the left-hand and right-hand sides of Eq. (1) are tangent and a small adjustment of the model parameters causes the solution to disappear. Although this represents a breakdown of the quantization condition and/or unphysical model parameters, we still keep these solutions in our statistical analysis as their effect is negligible. In addition to the bootstrap resamples, we also perform an additional fit using the alternative spectrum and take the difference in fit results as an estimate of systematic uncertainty. A similar problem occurs for the bootstrap estimate of the uncertainty of the interacting spectrum in the continuum obtained using Eq. (1) and shown in Fig. 3. When L is small, for some of the samples the ground state solution in frames (0, 0, 0) and (0, 0, 2) disappears. Because of this, we do not show an error band for these cases, which correspond roughly to energies below the t-channel cut. To estimate additional systematic uncertainty due to the continuum extrapolation and residual finite-volume effects, we apply various cuts to the selection of ensembles. In addition, to probe the ansatz for p cot δ, we use both a quadratic polynomial in p 2 with the full p 2 range and a linear polynomial with the near-threshold region. These fits are summarized in the first three groupings of Fig. S8. When fitting to all six lattice spacings, the resulting binding energy is very stable with respect to the inclusion of the small volumes and the choice of p 2 range. As the coarser lattice spacings are excluded, the variations increase, with the choice of p 2 range becoming more important than the cut on L. Our preferred fit, which provides our central value and statistical uncertainty, is the one that includes the most data, i.e. the first in the figure. We estimate the systematic uncertainty as the root-mean-square difference from the preferred fit of the central values of the seven other fits in the first and third groupings. Figure S8 also shows two additional variations. Including the ground state in the rest frame from N202 has a negligible impact on the binding energy but can substantially increase χ 2 . Including a 3 terms in the dependence on the lattice spacing significantly increases the uncertainty, without improving the fit quality.
The curves showing the dependence of B H on a 2 in Fig. 1 are based on the fits whose results are shown as filled blue squares and orange circles in the first three groupings of Fig. S8. The same is shown for a −1 0 and r 0 in Fig. S9. The inverse scattering length shows a strong dependence on the lattice spacing (varying by a factor of two) and is fairly insensitive to the choice of p 2 range. The effective range has a weaker dependence on the lattice spacing but shows larger variation with the choice of p 2 range; this contributes to its relatively larger systematic uncertainty.

VI. TWO-FLAVOR ENSEMBLES
In addition to our main analysis of N f = 3 lattice ensembles, we have generated new data for two N f = 2 ensembles (i.e. with dynamical u and d quarks and a quenched s quark) used in our previous study of the H dibaryon [17] and listed in the lower part of Table SI. Based on the analysis in Ref. [80], we expect that the quenched s quark is not an obstacle to using finite-volume quantization conditions. On both ensembles we elected to set the strange quark mass equal to that of the light quarks; this means that both ensembles have SU(3) flavor symmetry in the valence sector. For ensemble E5 with a pion mass of 437 MeV, this is a change from Ref. [17] where we tuned the strange quark mass to be near its physical value; as a result, the main difference between E5 and the N f = 3 ensembles is that the strange quark is quenched.
Our analysis on the N f = 2 ensembles is the same as what was done in the N f = 3 case, except that we cannot study the continuum limit. The finite-volume spectra obtained from ensemble E5 are shown in the left panel of Fig. S10; they have the same qualitative features as observed for the N f = 3 ensembles. Performing fits of the phase shift, we obtain a binding energy B H = 12.0 ± 2.7 ± 0.5 MeV (E5), which is consistent with the binding energies in the N f = 3 case at similar nonzero lattice spacing. The right panel of Fig. S10 shows the finite-volume spectra for ensemble E1. As the pion mass is much larger, the t-channel cut and three-particle threshold are further away from the threshold and all of the obtained levels lie in the region where the two-particle quantization condition is applicable. On this ensemble, the uncertainty of both the spectrum and the fitted quantities are dominated by systematics. The phase shift fits yield B H = 17.3 ± 4.0 ± 5.4 MeV (E1), which is consistent with the value 19 ± 10 MeV reported in our previous work [17] but has a smaller error. In Fig. 5 we compare our results for the binding energy B H for N f = 3 and N f = 2 with the estimates from HAL QCD [32,33] and NPLQCD [27,28,30]. Our N f = 2 calculations at nonzero lattice spacing show a dependence on the pion mass compatible with that observed by HAL QCD, although they lack the precision necessary to make an unambiguous statement. Moreover, this plot underscores our observation that discretization effects in this quantity are sizeable.
The first entry of the first dimension contains the average over the ensemble and the next 1000 are the bootstrap samples. For the two-baryon spectrum, the second dimension indexes the preferred and alternative values in the first and second entries, and the third dimension indexes the energy levels in ascending order. As should be evident from the names, the first of these three datasets contains the octet baryon energy with momentum P , the second contains the two-baryon level identified as spin one, and the last contains both two-baryon-levels identified as spin zero. FIG. S10. Two-baryon spectrum in five different reference frames on ensembles E5 (left) and E1 (right). Green points are the spin-zero levels and gray points are the spin-one levels. Horizontal lines indicate two-and three-particle thresholds and the t-channel cut. Horizontal line segments show finite-volume energies in the noninteracting case (dashed red) and from the fit to the wider p 2 range (solid blue).