Lattice QCD study of the $H$ dibaryon using hexaquark and two-baryon interpolators

We present a lattice QCD spectroscopy study in the isospin singlet, strangeness $-2$ sectors relevant for the conjectured $H$ dibaryon. We employ both hexaquark and two-baryon interpolating operators to isolate the ground state in the rest frame and in moving frames. Calculations are performed using two flavors of O($a$)-improved Wilson fermions and a quenched strange quark. Our initial point-source method for constructing correlators does not allow for two-baryon operators at the source; nevertheless, results from using these operators at the sink indicate that they provide an improved overlap onto the ground state in comparison with the hexaquark operators. We also present results, in the rest frame, using a second method based on distillation to compute a hermitian matrix of correlators with two-baryon operators at both the source and the sink. This method yields a much more precise and reliable determination of the ground-state energy. In the flavor-SU(3) symmetric case, we apply L\"uscher's finite-volume quantization condition to the rest-frame and moving-frame energy levels to determine the $S$-wave scattering phase shift, near and below the two-particle threshold. For a pion mass of 960 MeV, we find that there exists a bound $H$ dibaryon with binding energy ${\Delta}E=(19\pm10)$ MeV. In the 27-plet (dineutron) sector, the finite-volume analysis suggests that the existence of a bound state is unlikely.


I. INTRODUCTION
The strong force between quarks and gluons produces a rich spectrum of bound states and resonances, the colorneutral hadrons. Most of these can be described by constituent quark models as either quark-antiquark mesons or three-quark baryons. The existence of exotic hadrons, which cannot be described as such, is an active field of inquiry. Over the past several years, the so-called "X, Y, Z" mesons have been studied intensively, both theoretically and experimentally [1], and in recent years pentaquark baryons have also gained attention [2].
Nearly four decades ago, using the MIT bag model, Jaffe predicted a deeply bound dibaryon with quark content uuddss that is a scalar and a flavor singlet, the H dibaryon [3]. In contrast with the only known stable dibaryon, the deuteron, which can be well described as a loosely-bound proton-neutron state and is bound by just 2.2 MeV, the bag model predicted the H dibaryon as an exotic hexaquark state where all six quarks are in S-wave in the same hadronic bag, bound by about 80 MeV below the ΛΛ threshold.
Experimental evidence disfavors such a large binding energy. The strongest constraint is the "Nagara" event provided by the E373 experiment at KEK [4], which found a 6 ΛΛ He double-hypernucleus with ΛΛ binding energy B ΛΛ = 6.91 ± 0.16 MeV [5] that decayed weakly. A deeply bound H dibaryon would enable the strong decay 6 ΛΛ He → 4 He + H; its absence implies m H > 2m Λ − B ΛΛ . There was also no indication of an H dibaryon from a high-statistics study of upsilon decays at Belle [6].
The first lattice QCD study of the H dibaryon was per-formed more than thirty years ago [7], using a quenched ensemble with lattice size 6 2 × 12 × 18. Quenched studies -which all used hexaquark interpolating operators in keeping with the bag model picture, together with standard lattice spectroscopy techniques -produced inconclusive results: while some found a bound state [8][9][10], others did not [7,[11][12][13]. Early studies of the H dibaryon using lattice QCD are summarized in Ref. [14]. Aside from the present work 1 , calculations with dynamical fermions have been performed by two collaborations, both of which reported a bound H dibaryon at heavier-than-physical quark masses. The NPLQCD collaboration performed lattice spectroscopy calculations using a setup based on clover fermions with local hexaquark operators at the source and bilocal two-baryon operators at the sink. First results were obtained on anisotropic ensembles with N f = 2 + 1 dynamical fermions [14,18,19], followed by isotropic ensembles with three mass-degenerate (N f = 3) quarks [20]. An alternative approach, employed by the HAL QCD collaboration, is based on determining baryon-baryon potentials from Nambu-Bethe-Salpeter wave functions computed on the lattice, followed by solving the Schrödinger equation to study baryon-baryon scattering and bound states. This was done on ensembles with N f = 3 clover fermions for a range of quark masses [21][22][23]. Although these two sets of calculations agreed on the presence of a bound state, they disagreed significantly on the binding energy: in the N f = 3 case with pseudoscalar meson mass near 800 MeV, the value reported by NPLQCD was 74.6 ± 4.7 MeV, whereas HAL QCD reported 37.8 ± 5.1 MeV. Recently, HAL QCD have published a N f = 2 + 1 study of coupled channel (ΛΛ and N Ξ) baryon-baryon interactions with near-physical quark masses, which claims that the H dibaryon may be a ΛΛ resonance just below the N Ξ threshold [24,25].
Given that there are conflicting results for the binding energy of the H dibaryon, we have started a new initiative which may help to resolve the issue. As a first step we present results from a study in two-flavor QCD, i.e. with a mass-degenerate doublet of dynamical u and d quarks. The mass of the (quenched) strange quark is either tuned such that m s = m d = m u or set to a heavier value, implying that the SU(3) flavor symmetry is broken. Clearly, SU(3) symmetry is significantly broken at the physical point, which allows the three flavor multiplets, i.e. the singlet, octet and 27-plet to couple. Therefore, it is advantageous to study the octet and 27-plet even in the case of exact SU(3) symmetry. Furthermore, the 27-plet contains the two-nucleon I = 1 sector which has a possible dineutron bound state. The nucleon-nucleon sector has been studied extensively in experiment and may serve as a benchmark for lattice calculations.
Our work is mainly focused on the methodology of determining the spectrum and the binding energy via the computation of correlation matrices and their diagonalization [26][27][28]. In order to allow for a direct comparison with the results from older quenched studies and from NPLQCD we have chosen a similar setup. As we will describe in more detail in the following sections, we have used point sources to compute correlator matrices with local (hexaquark) interpolating operators at the source and both local and bilocal (two-baryon) interpolators at the sink. In addition, we report initial results from a follow-up study in which we, for the first time, applied the distillation method [29] to the two-baryon sector, which allowed us to compute a correlator matrix with two-baryon operators at both the source and the sink. We shall see that the hermitian setup leads to a more robust and precise identification of the spectrum.
Since the strange quark is quenched in our calculation, one may think that any observed deviation from the findings of Refs. [14,[18][19][20][21][22][23] should be attributed to the different treatment of the quark sea. However, the FLAG report [30] provides ample evidence that observables computed with N f = 2 or N f = 2 + 1 dynamical quarks differ at the percent level at most. This paper is organized as follows. Our methodology is described in Section II: this includes the interpolating operators and our approach for analyzing correlator matrices. We show our determination of the energy levels using point-source methods in Section III A and using distillation in Section III B. In Section IV, we apply Lüscher's finite-volume quantization condition to determine scattering phase shifts at the SU(3)-symmetric point, and identify the presence of a bound H dibaryon. Finally, our conclusions are presented in Section V.

A. Simulation details
Our study has been performed on a set of ensembles with two mass-degenerate dynamical flavors of O(a)improved Wilson quarks [31] and the Wilson plaquette action, which were generated as part of the CLS (Coordinated Lattice Simulations) initiative, using the deflationaccelerated DD-HMC [32,33] and MP-HMC [34] algorithms. The improvement coefficient c sw multiplying the Sheikholeslami-Wohlert term was tuned according to the non-perturbative determination of Ref. [35]. An overview of the ensembles can be found in Table I. All our calculations were performed in the SU(3)-flavor symmetric limit, with the exception of ensemble E5 for which the valence strange quark mass was tuned so that the combination (2m 2 K − m 2 π )/m 2 Ω takes its physical value. The values of the lattice spacing in physical units were determined using the kaon decay constant [36].
Quark propagators were computed using the SAP domain-decomposed, deflated GCR-solver of the DD-HMC package [32] with smeared point sources on a grid of source positions that was randomly displaced on each gauge configuration. For the distillation calculation, we used a similar solver in OpenQCD [37], computed low modes of the spatial Laplacian using PRIMME [38], and contracted them to form "perambulators" and mode triplets using QDP++ [39]. Baryon and multi-baryon correlators computed in lattice QCD suffer from a severe signal-to-noise problem, since the noise grows with a rate proportional to exp{(m B − 3/2m π )t} per baryon, where m B denotes the baryon mass. This makes it difficult to identify a "window" in which the asymptotic behavior has been reached while the signal is not yet lost in the statistical noise. In order to allow for a significant increase in statistics while keeping the numerical effort at a manageable level, we have employed the method of allmode-averaging (AMA) [40]. This entails computing a high number of samples with lower-precision propagator solves, followed by applying a bias correction using a relatively small number of high-precision solves. We used this for the calculation with point-source propagators, but obtained only modest cost savings due to our use of a highly efficient solver.

B. Interpolating operators
Accurate determinations of the spectrum in the H dibaryon channel require a set of efficient interpolating operators whose projection properties onto the ground state may also give qualitative insights into the nature of the H dibaryon. For instance, the local hexaquark operators defined below resemble more closely Jaffe's original interpretation of the H dibaryon as a deeply bound state of six quarks forming a color singlet. By contrast, bilocal two-baryon operators may be more appropriate to describe loosely-bound states such as the deuteron. While the two types of operators are defined according to a qualitative physical picture that is suggestive of the nature of a given state, this is not a rigorous way to study its properties. The generic form of a lattice QCD correlation function is given by where the interpolating operator O i carries the quantum numbers of the continuum state under study, and it is understood that O i has been projected onto spatial momentum P . When constructing interpolators in the H dibaryon channel, one can think of two generic configurations. Jaffe's original analysis was based on a compact color-singlet comprising six quarks, which gives rise to a hexaquark operator composed of flavors uuddss. The alternative possibility is the product of two individual color-singlets at different positions, i.e. a two-baryon operator.
The starting point for the construction of hexaquark operators is the object where r, s, . . . , w denote generic quark flavors, and P + = (1 + γ 0 )/2 projects the quark fields to positive parity. One can form two operators that transform under the singlet [12,41,42] and 27-plet irreducible representations of flavor SU(3): The continuum quantum numbers of these operators are S = −2, I(J P ) = 0(0 + ), as required by the original bagmodel proposal. We note that an octet state cannot be represented in terms of this simplest class of hexaquark operators. Finally, we project onto the momentum of each lattice frame: An alternative configuration in the H dibaryon channel is described by a bilocal two-baryon operator [21]: where the individual baryons have been projected onto spatial momenta p 1 and p 2 , and P = p 1 + p 2 denotes the total momentum. The index m labels a particular configuration of momenta p 1 and p 2 . In the rest frame, the momentum combinations have a particularly simple construction: We form dibaryon operators from combinations of octet baryons with S = −2, I = 0, i.e.
where the subscript S on (N Ξ) denotes the flavorsymmetric combination. 2 Using the rotation matrices listed in Appendix B of Ref. [21] we form the appropriate linear combinations of (ΛΛ), (N Ξ) S and (ΣΣ) that correspond to different flavor multiplets, i.e. the singlet, octet and 27-plet. The projected operators are then called BB {1,8,27},m .
In the interpolating operators defined above, we use smeared quark fields; the smearing helps to increase the coupling of an operator to the low-lying states. For the point-source calculation, we used Wuppertal smearing [43], with the hopping term constructed using spatially APE-smeared [44] gauge links. To increase the size of our operator basis, we used two different smearing widths: "narrow" (70 steps, denoted N ) and "medium" (140 steps, denoted M ). The distillation approach makes use of Laplacian-Heaviside (LapH) smearing [29], in which the quark fields are smeared by projecting them onto the low-lying modes of the spatial gauge-covariant Laplacian, itself constructed using stout-smeared [45] gauge links. We used 56 modes in all cases, and label this type of smearing as L. We include the smearing as part of the label for each interpolating operator, yielding names such as H 1,N or BB 27,L,0 .

C. Operator basis and correlation matrices
The determination of hadronic energy levels in lattice QCD usually proceeds by computing a correlation matrix C ij (t) for the chosen basis of interpolating operators [see Eq. (1)] and solving a generalized eigenvalue problem (GEVP) [26][27][28].
In the following the correlation matrices C ij (t) = C ij (t, P ) are evaluated at rest, i.e. with total momentum P = 0, and in moving frames, i.e. P 2 > 0. We included two-baryon operators with individual momenta p 1,2 = 2πd 1,2 /L and d 2 = 0, 1, 2, 3 in the calculation. The components were chosen such as to realize a total momentum P = p 1 + p 2 ≡ 2πD/L with D 2 = 0, 1, 2, 3. For each D 2 , we average over all equivalent frames that are related by a lattice rotation.
Performing the Wick contractions of correlators involving bilocal two-baryon operators of the type in Eq. results in diagrams that contain quark lines that start or end at two distinct spatial points within a timeslice, This has important consequences for computing correlation matrices, since point-to-all propagators do not allow for the computation of such diagrams whenever a twobaryon operator is placed at the source. We therefore opted for an asymmetric setup in which only hexaquark operators were put at the source, while at the sink both hexaquark and two-baryon operators were used; this is illustrated in Figure 1. In this setting, the familiar GEVP or, more generally, the diagonalization procedure must be modified in order to allow for a non-hermitian correlator matrix. We select subsets of N op source operators and N op sink operators and form the corresponding square correlator matrix C(t). We then perform the following steps, starting from this N op × N op matrix: 1. Determine the right and left eigenvectors of C(t) by solving for n = 1, . . . , N op , where t 0 and t 1 denote fixed timeslices in the region where the lowest N op states are expected to dominate.
2. Compute the (approximately) diagonal matrix Λ(t) whose elements are given by 3 3. The effective nth energy level is then obtained from the diagonal element Λ nn (t) via the well-known formula We have used a timestep of ∆t = 3a in our analysis. While the choice of interpolators at the source is restricted to hexaquark operators, we can probe a number of different operators at the sink and study their relevance for determining the ground state.
Relying on non-hermitian correlator matrices makes it more difficult to identify the ground state reliably, owing to the fact that the diagonalization does not project exactly on the correlator corresponding to the nth energy eigenstate. It is then not guaranteed that the effective energies computed via Eq. (15) approach the asymptotic value monotonically from above, since the statistical weights of different states are not strictly positive.
In order to overcome this difficulty we have implemented distillation and LapH smearing [29]. Since this is a timeslice-to-all rather than point-to-all method (see Figure 2), it allows us to compute a hermitian correlator matrix using the two-baryon operators listed in Eqs. (9)- (11), both in the center-of-mass frame and for non-vanishing total momentum P . The main objects that we compute are the perambulator, which is the projection of the quark propagator onto the low modes of the spatial Laplacian, and the mode triplets, where u (l,t) is the lth low mode on timeslice t. The correlator matrix is then computed by contracting these two objects. Further details of our implementation will appear in future work; here we present initial results in the rest frame using operators with both baryons at rest, i.e. p 1 = p 2 = 0. For a hermitian correlator matrix, the right and left eigenvectors are obviously identical.
In the case of broken SU(3) symmetry, the diagonalization method also allows us to associate a state with a particular multiplet, by identifying which operators couple most strongly to it. Furthermore, an essential feature of this method -as will be shown in Section III A 2 -is that it can identify the presence of more than one state, even when the statistical signal is too poor to distinguish the energies of those states. 3 Note that Λ(t) is exactly diagonal for t = t 0 and t = t 1 .

III. ENERGY LEVELS
We now present our results for the ground-state energies in the dibaryon channel, as determined using either point-to-all or timeslice-to-all propagators on the ensembles listed in Table I.

A. Point-to-all propagators
Correlators based on point-to-all propagators were computed for ensembles A1, E1, E5 and N1. On all four ensembles we have computed correlators in the rest frame and in three moving frames. Our main findings are most easily explained for the data extracted from ensembles E1 and E5 in the rest frame, which show the highest level of statistical precision.
While A1, E1 and N1 realize an SU(3)-symmetric situation, SU(3) symmetry is broken for ensemble E5. Note that the analysis for the SU(3)-symmetric case is simplified, due to the fact that different flavor multiplets (singlet, octet and 27-plet) cannot mix. Therefore, we discuss the two cases separately.

Analysis of the SU(3) symmetric case
As outlined in Section II C and illustrated in Figure 1, our setup based on point-to-all operators does not allow us to put two-baryon operators at the source. Hence, we have constructed 2 × 2 correlator matrices using hexaquark operators of the two different smearing types (H 1,N and H 1,M for the singlet case and similarly for the 27-plet) at the source, while at the sink we have used either hexaquark or two-baryon operators. In the latter situation the resulting correlator matrix is non-hermitian, as described in Section II C. After applying the appropriate diagonalization procedure we have studied the behavior of the effective energies defined according Eq. (15), as a function of the Euclidean time separation t.
Results for the ground-state energy of the singlet channel on ensemble E1 are shown in the top panel of Figure 3. The effective energy for the first excited state is too noisy to be displayed in the plot. We find that the "narrow" smearing (N ) is the most effective in obtaining clean signals and therefore use only narrow-smeared operators at the sink to determine our final results.
The green, blue and red filled symbols in Figure 3 denote the energy levels extracted from different correlator matrices, whose operator choice at the sink is described in the legend. The open red circles denote the effective energy determined from a single correlator, comprising a hexaquark operator at the source and a two-baryon operator at the sink, with the two momenta each set to zero. The horizontal line and grey band represent the energy level corresponding to 2m Λ .
Our main observations are as follows: The effective energy determined via the diagonalization of the hermitian  2 × 2 hexaquark correlator matrix approaches its plateau from above. In the plateau region the data are noisy and compatible with the ΛΛ threshold. By contrast, correlator matrices including at least one two-baryon operator at the sink yield effective energies that are below the threshold and which have smaller uncertainties. However, care must be taken when deciding at which value of t the asymptotic behavior has been isolated. Owing to the non-hermitian setup, it is possible that residual excitedstate contributions enter the projected ground-state correlator with negative weights, so that the plateau is approached from below. This can also cause local minima to appear, which could be difficult to distinguish from a true plateau. Indeed, we see evidence for this behavior, with the energies showing a dip for t ≈ 0.4 fm before moving closer to the threshold for t 0.9 fm. The effective energy determined from the mixed single correlator approaches a plateau below the threshold, but without showing any dip. We conclude that the ground-state effective energy shows a consistent plateau for t 0.9 fm, which is interpreted as the ground-state energy. The 2×2 hexaquark correlator matrix yields consistent results for t 1.0 fm, given the large statistical noise.
Our estimate of the ground-state energy in the singlet channel on ensemble E1 is obtained from a fit to the effective energy determined from the diagonalization of the 2 × 2 correlator matrix with one hexaquark and one twobaryon operator at the sink. The panel on the top right of Figure 3 shows a blow-up of the plateau region, with the fitted value of the energy and error displayed as a band across the fitting interval. We find that the ground state in the singlet channel lies below the energy of two non-interacting Λ hyperons by 2.5 standard deviations. Numerical values for the fitted energies are listed in Tables II and III, while the results for the energy difference 2m Λ − E are shown in Table IV. The qualitative features observed for the 27-plet are very similar to the singlet channel (see bottom of Figure 3): Whereas the effective energy extracted from a correlator matrix constructed from only hexaquark operators shows no sign of lying below the ΛΛ threshold, we find that using two-baryon operators results in significantly lower values. Fitting the latter in the region where t > 0.9 fm produces an estimate that lies within one standard deviation below the threshold (see Tables V  and VI).

Analysis of the SU(3)-broken case
The breaking of SU(3) symmetry induces mixing among the relevant multiplets (singlet, octet and 27plet), which makes their identification a more involved task. The appropriate strategy is to construct correlator matrices from operators projected onto all relevant flavor multiplets. We recall that the octet can only be represented in terms of two-baryon operators, which, however, cannot be placed at the source when using point-to-all propagators. Therefore, the current analysis is focused only on the singlet and 27-plet states. This deficiency will be rectified by the calculation of hermitian correlator matrices using distillation, described in the next subsection.
In the presence of mixing among different SU(3) multiplets, we have constructed 4 × 4 correlator matrices from the same combinations of operator types as in the SU(3)symmetric case, while including both projections onto the singlet and 27-plet.
The matrix correlator and its diagonalization provides us with important information for the interpretation of our results. First, the amount of SU(3) breaking is small, as evidenced by the fact that the geometric mean of the flavor off-diagonal correlators is less than 0.5% of that of the corresponding flavor-diagonal ones. Second, the diagonalization provides clear evidence for the existence of two different low-lying states, although it is not possible to resolve their energy difference with the current statistics. Furthermore, states corresponding to different flavor multiplets can be identified with the help of the eigenvectors computed via the diagonalization procedure. We find that the eigenvectors corresponding to the two low-lying states are dominated by operators in a single multiplet. In this way we can unambiguously assign the energy levels to a particular flavor multiplet.
Results for the effective energies for the SU(3) singlet and 27-plet are shown in the top and bottom panels of Figure 4, respectively. Data points shown in a given color correspond to a particular choice of sink operators, as indicated in the legend. For both the singlet and the 27plet, it is evident that the asymptotic behavior has not set in below t ≈ 1.0 fm. Furthermore, the rapid growth of statistical noise for t > 1.0 fm precludes the identification of a clear plateau below the ΛΛ threshold, irrespective of the composition of the operator basis at the sink. Statis-tically, the most precise signal is obtained by considering a 2 × 2 correlator matrix, composed of "narrow" smeared hexaquark interpolating operators for the singlet and 27plet at the source and a corresponding set of two-baryon operators with both momenta set to zero at the sink. The corresponding effective energy is represented by the open red circles in Figure 4. After fitting the effective energy to a constant in the interval t = 1.0 − 1.3 fm we obtain the estimates represented by the red bands in the right panels. We find that the fitted ground-state energies of the singlet and 27-plet are of the order of one standard deviation below the ΛΛ threshold (see Tables II -VI).
As discussed in Section II C we have computed correlation matrices in frames with total momentum P = 2πD/L with D 2 = 0, 1, 2, 3. Estimates for the corresponding ground-state energies are listed in Tables III -VI. The combined results in all frames have been used to obtain information on the scattering phase shifts and binding energies in infinite volume. A detailed discussion is deferred to Section IV.

B. Dibaryon analysis with distillation
The distillation technique enables the use of twobaryon operators at both the source and the sink. This results in a hermitian correlator matrix and allows us to properly account for the mixing between the singlet, octet and 27-plet states under SU(3) breaking. Within this setup we do not use hexaquark operators. For our initial study, we have computed correlation functions on ensembles E1 and E5 in the rest frame only, which allows for a direct comparison with results from the previous subsection.
As described in Section II B, we start from the twobaryon interpolating operators of Eqs. (9)-(11), projected on p 1 = p 2 = 0. Using the transformation of Appendix B in Ref. [21] we then perform the rotation to the singlet, octet and 27-plet, which yields the operators BB 1,L,0 , BB 8,L,0 and BB 27,L,0 , respectively. In each channel we compute the corresponding correlation functions and determine the effective energies. Results for the effective energy on E1 are shown in the top left panel of Figure 5. In all three channels, the effective energy approaches a plateau monotonically from above and do not show any of the irregularities observed in the nonhermitian setup based on point-to-all propagators. The top right panel of Figure 5 shows the effective energy difference E eff − 2m Λ,eff . The plot demonstrates that three separated levels can be resolved, corresponding to the singlet, the 27-plet and octet (in ascending order). The fact that the effective energy differences are nearly constant over time indicates that excited-state contributions in the dibaryon correlation functions are very similar to those in the single Λ correlator. Still, it is important to determine the energy difference in the regime where the individual correlators have reached their respective asymptotic behavior. In the singlet channel (color-coded Due to the breaking of SU(3) symmetry, the different multiplets cannot be treated independently anymore and it is here that the inclusion of the octet operator at the source and sink allows for a consistent treatment of the relevant flavor multiplets. In order to take account of the mixing of different multiplets, we compute a 3 × 3 correlator matrix using the operators BB 1,L,0 , BB 8,L,0 and BB 27,L,0 , which is then subjected to the diagonalization procedure. The resulting energy levels are identified with a particular multiplet according to the overlaps with the interpolating operator. The effective energy of the singlet (blue circles) is below the ΛΛ threshold and, in contrast to the case of point-to-all propagators, can be clearly distinguished from the energy corresponding to two noninteracting Λs. The distillation technique allows us to distinguish the singlet and 27-plet as two distinct states with our statistics. The octet, by contrast, lies far above the threshold. For our final estimates of the energies of the different states in the SU(3)-symmetric situation, we have fitted the effective energy to a constant for t starting from 0.8 fm. In the SU(3)-broken situation, the plateau starts later and we have fitted from 1.0 fm.
It is instructive to compare the energy levels computed using point-to-all propagators to the results extracted using the distillation technique. In Figure 6 we show the effective energy gap for the singlet ground state computed on ensemble E1. The plot demonstrates the dramatic improvement in the overall precision provided by the distillation technique in conjunction with LapH smearing: The effective energy gap computed using distillation can be clearly distinguished from zero owing to its tiny statistical errors. Furthermore, it exhibits a flat behavior across a wide range of Euclidean times. While the results obtained using point-to-all propagators are consistent with distillation for t 0.9 fm, the larger statistical noise makes it more difficult to determine a non-vanishing energy gap.
We end this section with a discussion of the relative cost of calculations based on point-to-all propagators and distillation. A simple and fairly accurate cost estimate is provided by the number of propagator solves per configuration. For point sources we have to perform N pt = N flav · N src · N smear · 12 (17) inversions, where N flav denotes the number of different quark flavors, N src is the number of sources, and the number of different smearings is given by N smear . For distillation the number of propagator solves is given by where N tsrc denotes the number of timeslices used to compute timeslice-to-all propagators, and N LapH is the number of low-lying modes of the spatial Laplacian. Using the information from Table I and the fact that we have used N LapH = 56 modes, we find N pt = 3072 and N dist = 1792 for ensemble E1. For E5 the above expressions evaluate to N pt = 1536 and N dist = 1792. We conclude that the better data quality of the distillation technique is achieved for comparable or even lower cost.
For a more precise cost comparison one should also take into account that the cost of contractions in the distillation approach is typically larger. The associated computational overhead may become significant on large volumes. We will present a more thorough discussion of this issue in a future publication.

IV. FINITE VOLUME ANALYSIS
The finite-volume rest frame energy level E, when below the two-particle threshold, provides a naïve estimate of the mass of the bound state. However, this can suffer from significant finite-volume effects that are asymptotically only suppressed as e −κL , where κ is the binding momentum defined via E = 2 m 2 Λ − κ 2 . In addition, there can be an energy level below threshold without a bound state, in which case it has a power-law dependence on L that depends on the scattering length.
For a more careful study of the presence of a bound state and its mass, we turn to Lüscher's finite volume quantization condition [46] and its extension to moving  6. Comparison of the effective energy difference between the singlet ground state and two non-interacting Λs, computed with point-to-all propagators and distillation, respectively, for ensemble E1. Open blue circles derive from the same non-hermitian 2 × 2 correlator matrix that was used to determine the ground state energy in the singlet channel, represented by the blue points in the top panel of Figure 3.
frames [47]. Below the three-particle threshold, this is a relation between the two-particle scattering amplitude and the finite-volume energy levels. In the case of one pair of identical particles scattering in the S wave, if we ignore the influence of higher partial waves due to the breaking of rotational symmetry, the condition takes the form where p is the scattering momentum satisfying E cm ≡ √ E 2 − P 2 = 2 m 2 Λ + p 2 , δ(p) is the scattering phase shift, γ = E/E cm is the moving-frame boost factor, and Z D 00 is a generalized zeta function defined in Ref. [47]. Our numerical implementation of the zeta function is based on Ref. [48].
As recently reviewed in Ref. [49], a bound state, corresponding to a pole in the scattering amplitude on the real p 2 axis below zero, is determined by the condition p cot δ(p) = − −p 2 . That reference also provides a check that can be applied to lattice data: at the pole, the slope of p cot δ(p) (versus p 2 ) must be smaller than that of − −p 2 .
For small p 2 , the phase shift can be described by the effective range expansion. We use the first two terms, where a 0 is the scattering length and r 0 is the effective range. Fitting this to lattice data is nontrivial, since p cot δ(p) and p 2 are not independent variables, being related by the quantization condition. Therefore, we choose to fit to the squared scattering momentum p 2 determined from each lattice energy, similarly to what was done for fitting to the lattice energies in Ref. [50]. Fitting to the momentum rather than energy benefits from cancellations of statistical uncertainties that are correlated between E and m Λ . Given a set of fit parameters (a 0 , r 0 ), in each frame the fit momentum is determined by finding the solution to the quantization condition that is nearest to the corresponding momentum determined from the lattice calculation, i.e. by numerically solving We only consider the ensembles with SU(3) symmetry, since otherwise this is a much more complicated coupledchannel (ΛΛ, N Ξ, and ΣΣ) system. We are only able to obtain a reliable fit in the flavor singlet sector on ensemble E1, where the precise energy level in the rest frame from the distillation method provides a stronger constraint than the other data. This is shown in Figure 7; we find the scattering length to be 1.3 (5) fm and the effective range 0.4(3) fm. There is a clear intersection with the bound-state curve, which has a slope with the correct behavior, indicating the presence of a bound H dibaryon. This intersection yields a binding energy of ∆E = 19(10) MeV, somewhat smaller in magnitude than the naïve value obtained from the rest-frame energy difference 2m Λ − E.
For the ensembles A1 and N1, we are unable to obtain reliable fits. However, in the flavor singlet sector, the data in the rest frame and in the higher moving frames sit on opposite sides of the bound-state curve, which indicates that there will be an intersection between it and the phase shift; see Figure 8. For the frame D 2 = 1, the point -together with its error along the quantization curve -is almost on top of the bound-state curve, and therefore we make a conservative estimate of the binding energy, based on the spread of E cm in this frame and in the rest frame. For N1, we get a∆E = 0.016 (13), or 65(53) MeV. For A1, we get a∆E = 0.028 (25), or 73(66) MeV.
In the 27-plet sector, our data are generally not precise enough to distinguish the energy levels from the nonin-teracting ones; this means that the phase shift is consistent with zero. On E1, where distillation provides a relatively precise value, the rest-frame energy is slightly below threshold. Using the quantization condition, we find that just below threshold, p cot δ is probably positive; see Figure 9. This means that a bound state (which would be a dineutron) is unlikely, however an intersection with p cot δ = + −p 2 , which corresponds to a virtual bound state (i.e., a pole on the unphysical sheet), is possible. If we neglect the dependence on p 2 , the distillation energy level implies a scattering length a 0 = −0.2 +0.1 −0.2 fm. However, note that the quantization condition is very nonlinear and at 3σ the statistical uncertainty on a 0 becomes as large as +0.3 −1.2 , assuming a symmetric Gaussian uncertainty on the p 2 determined from the energy level.

V. DISCUSSION AND CONCLUSION
We have presented a detailed study of the spectrum of the H dibaryon in the SU(3) flavor-symmetric and broken cases, using lattice QCD simulations with dynamical up and down quarks and a quenched strange quark. We have analyzed the efficiency of different interpolating operators in the determination of hadronic finite volume energy eigenvalues via variational analysis.
Our findings indicate that point-like hexaquark operators provide a poor overlap onto the SU(3) singlet ground state. This becomes evident through the slow convergence to the ground state plateau seen in the time dependence of the lowest energy eigenvalue. The inclusion of bilocal two-baryon operators at the sink improves the overlap considerably, as indicated by an earlier onset of the plateau and a stronger statistical signal. However, in a setup based on point-to-all propagators, the improvement comes at the expense of having to deal with non-hermitian correlator matrices. As a consequence, effective energies determined from the eigenvalues are not guaranteed to approach their asymptotic behavior from above. In this way an additional systematic is introduced, as the onset of the plateau is obscured or can be misidentified.
A fully hermitian setup with bilocal two-baryon operators at both the source and sink can be realized with the help of the distillation technique which allows for the calculation of timeslice-to-all propagators. In our study of the SU(3) symmetric and broken cases we have been able to determine the energy levels reliably and with good statistical precision for the flavor singlet, 27-plet and octet.
In particular, we could identify a clear gap between the ground-state singlet energy level and the 2m Λ threshold, which suggests that there is a bound H dibaryon. The energy difference provides a naïve estimate of the binding energy, ∆E ≡ 2m Λ − E. For ensembles E1 and E5, which realize the SU(3) symmetric and broken situations respectively, the estimates are E1: ∆E = 39.0 ± 2.2 MeV, m π = 960 MeV, (22)  In addition, we have applied Lüscher's finite-volume quantization condition to determine the scattering phase shift, which we use to obtain a more reliable estimate of the binding energy. Including data from the rest frame as well as several moving frames in the finite-volume analysis of the singlet case results in a shallower binding energy compared to the naïve estimate, viz.
Repeating the analysis for the SU(3) 27-plet is made more difficult through the relative closeness of the energy levels to the non-interacting levels, and the same is true in the case of broken SU(3)-flavor symmetry. We expect the situation to improve in our ongoing work, since we will obtain more precise results by also using distillation in the moving frames. In Figure 10 we show a compilation of our results for the binding energy on all our ensembles, plotted against the pion mass. The comparison with the estimates of the NPLQCD [14,18,20] and HAL QCD [22,23] collaborations is made in Figure 11. In the SU(3)-symmetric case we find that our estimate is considerably smaller than the result quoted by NPLQCD at a similar value of the pion mass [18]. Potentially, uncontrolled systematics such as the incorrect identification of the plateau, quenching of the strange quark or finite-volume effects could be the source of this discrepancy. In particular, recalling that finite-volume effects are asymptotically suppressed as e −κL , where κ is the binding momentum, the naïve estimates of the binding energy in Eqs.  (24) to the estimates quoted by NPLQCD [14,18,20] and HAL QCD [22,23]. Green and blue symbols refer to the SU(3)-symmetric and broken cases, respectively. and (23) may not be very reliable, given that κL evaluates to 2.99(8) for E1 and 1.74(25) for E5. We will address these issues in a future publication based on ensembles with N f = 2 + 1 flavors of dynamical quarks. Our findings suggest that the combination of distillation and Lüscher's finite-volume formalism will allow for a considerably improved calculation of the binding energy. Thus, there are good prospects for a reliable determination of this quantity at the physical point.