Strongly Bound Dibaryon with Maximal Beauty Flavor from Lattice QCD

We report the first lattice QCD study of the heavy dibaryons in which all six quarks have the bottom (beauty) flavor. Performing a state-of-the-art lattice QCD calculation we find clear evidence for a deeply bound $\Omega_{bbb}$-$\Omega_{bbb}$ dibaryon in the $^1S_0$ channel, as a pole singularity in the $S$-wave $\Omega_{bbb}$-$\Omega_{bbb}$ scattering amplitude with a binding energy $-81(_{-16}^{+14})$ MeV. With such a deep binding, Coulomb repulsion serves only as a perturbation on the ground state wave function of the parameterized strong potential and may shift the strong binding only by a few percent. Considering the scalar channel to be the most bound for single flavored dibaryons, we conclude this state is the heaviest possible most deeply bound dibaryon in the visible universe.

Understanding baryon-baryon interactions from first principles is of prime interest in nuclear physics, cosmology and astrophysics [1][2][3][4]. Dibaryons are the simplest nuclei with baryon number 2, in which such interactions can be studied transparently. However, the only known stable dibaryon is deuteron and the possible observation of perhaps just one more unstable light dibaryon [d * (2380)] has recently been reported [5,6] Even so, based on the theory of strong interactions, one expects to have more dibaryons in nature, particularly with the strange and heavy quark contents. Ab initio theoretical investigations using lattice QCD are well suited for studying such hadrons and indeed it can play a major role in their future discovery.
Lattice QCD calculations of dibaryon systems are becoming more feasible now particularly in the light and strange quark sectors [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. Even so, such studies involving heavy flavors are limited to only a few calculations [22][23][24][25]. Among the heavy dibaryons, a system of two Ω QQQ baryons (Q ≡ c, b) provides a unique opportunity to investigate baryon-baryon interactions and associated nonperturbative features of QCD in a chiral dynamics free environment. Such a system in the strange sector have been studied using Lüscher's finite-volume formalism [26], which suggested that the Ω-Ω channel is weakly repulsive [10]. Another study using the HALQCD procedure [27] suggested that the system is not attractive enough to form a bound state [12]. A recent high statistics HALQCD study [14] on a very large volume (∼ 8 fm) claimed that such a system is weakly attractive and the strength of potential is enough to form a very shallow bound state. Although the inferences from different procedures differ, they all agree on the fact that the interaction between two Ω baryons is weak. Another recent HALQCD investigation of Ω ccc -Ω ccc dibaryon reported a shallow bound state in the 1 S 0 channel [24]. While all these investigations suggest that the interactions in two Ω QQQ baryon systems are rather weak with quark masses ranging from light to charm, several lattice studies in the recent years on heavy dibaryons [23,25] and heavy tetraquarks [28][29][30][31] have shown that multihadron systems with multiple bottom quarks can have deep binding. Hence, it is very timely to study Ω bbb -Ω bbb interactions using lattice QCD. Note that very little is known about it through other theoretical approaches. [32][33][34].
The motivation for such a study is multifold. Theoretically it can provide an understanding of the strong dynamics of multiple heavy quarks in a hadron. In cohort with results from single- [10,12,24,35], double- [11, 15, 18-21, 23, 25, 36-38], and triple-flavored dibaryons [7-9, 11, 13, 16-18, 39], particularly those with heavier quarks, one would be able to build a broader picture of the baryon-baryon interactions at multiple scales. This can illuminate the physics of heavy quark dynamics in nonmesonic hadrons. A study of the quark-mass dependence of scattering parameters can further shed light into the dominant dynamics in different regimes. Indication of possible promising channels on any bound heavy dibaryon from such studies can also stimulate future experimental searches for them, as in the case of heavier tetraquarks [40][41][42][43].
In this Letter, we report the first lattice QCD investigation of the ground state of the dibaryons with the highest number of bottom (beauty) quarks in the 1 S 0 channel. We name it D 6b ≡ Ω bbb -Ω bbb , a dibaryon formed out of a combination of two Ω bbb baryons. Using various state-of-the-art lattice QCD utilities and methodologies, we extract the mass of D 6b and find clear evidence for a strongly bound state, with a binding energy of −81( +14 −16 )(14) MeV, and a scattering length of 0.18( +0.02 −0.02 )(0.02) fm. Despite its compactness, we find the Coulomb interactions act only as a perturbation to the strong interactions and do not change the binding in any significant way. Upon comparison to the binding energies of other dibaryons, e.g. 2.2 MeV of deuteron, and other strange or heavy dibaryons [23,24], we conclude D 6b to be the most deeply bound heaviest possible dibaryon in our visible universe.
The lattice setup that we use here is similar to the one used in Refs. [30,44] and we discuss it below. Lattice ensembles:− We employ four lattice QCD ensembles with dynamical u/d, s and c quark fields, generated by the MILC Collaboration [45] with highly improved staggered quark (HISQ) fermion action [46], as shown in Fig. 1. Lattice spacings are determined using r 1 parameter [45], which are found to be consistent with the scales obtained through Wilson flow [47]. Bottom quarks on lattice:− Since the bottom quark is very heavy, we use a nonrelativistic QCD (NRQCD) Hamiltonian [48], including improvement coefficients up to O(α s v 4 ) [49]. Quark propagators are calculated from the evolution of NRQCD Hamiltonian with Coulomb gauge fixed wall sources at multiple source time-slices. We tune the bottom quark mass using the Fermilab prescription for heavy quarks [50] in which we equate the lattice-extracted spin-averaged kinetic mass of the 1S bottomonia states with its physical value [51]. Such a tuning was also used in Refs. [23,30,44] and was found to reproduce the physical value of the hyperfine splitting of 1S bottomonia. (Di)baryon interpolators:− For the single Ω bbb baryon, we use the quasilocal nonrelativistic operator with J P = 3/2 + , as was used in Ref. [10]. This operator was constructed by the LHPC Collaboration and is listed in Table VII of Ref. [52] and also detailed in Ref. [53]. For extracting the ground state mass we assume only S-wave interactions in two baryon systems where the overall state is antisymmetric under the exchange of two baryons. Denoting components of the J=3/2 Ω bbb operator (O Ω bbb ) with χ m , m being the azimuthal component of J, we construct the Ω bbb -Ω bbb dibaryon operators as, Here, [CG] mn are the relevant spin-projection matrix constructed out of the appropriate Clebsch-Gordon coefficients. The J = 0 dibaryon operator that we employ in this work, is given by [10,53], Using these baryon and dibaryon operators (O Ω bbb and O D 6b ) we compute two-point correlation functions be-tween the source (t i ) and sink (t f ) time slices, At the sink time slice, we use several different quark field smearing procedures to identify the reliable ground state plateau and quantify possible excited state contamination (see Ref. [53] for more details). Ground state masses for the single and the dibaryon are obtained by fitting the respective average correlation function with a single exponential at large times (τ = t f − t i ). While determining mass in a lattice calculation it is often useful to plot the effective mass, defined as m ef f a = log[ C(τ ) / C(τ + 1) ], to show the signal saturation and justify the time window to be chosen in the exponential fit. In Fig. 2, we present the effective masses for C 2 Ω bbb (green circles) and C D 6b (blue squares) on the finest ensemble (a ∼ 0.06 fm) using wall quark sources and point quark sinks. We make the following observations from this result: (i) The signal in the effective masses saturates well before the noise takes over, and hence one can reliably extract the respective ground state masses. (ii) The signal in the noninteracting 2Ω bbb level survives until large times. This is because 2Ω bbb level is obtained from the single baryon Ω bbb correlator that decays with an exponent of M Ω bbb < M D 6b , and hence can propagate further than the D 6b state. (iii) Most importantly, it is quite evident that there is a clear energy gap between the ground state energy levels of the noninteracting two-baryon and the dibaryon systems at all times. This clearly shows that the ground state mass of dibaryon M D 6b is smaller than that of the non-interacting level 2M Ω bbb . We find similar energy differences for all the ensembles and we discuss the results below. Based on the t min dependence of the fits, which are discussed in Ref. [53], we make our final choices for the fit ranges and uncertainties arising out of such choices.
In order to gauge the extent of excited state contaminations in our estimates, we carry out two additional calculations: one with a wall source and a Gaussian- smeared sink [54,55], and the other with a wall source and spherical-box sink [56]. The results are detailed in the Supplemental Material [53]. We find that results are clearly consistent between different measurement setups and validate our estimates. We pass the results from all these different smearing procedures through the scattering analysis, as discussed below, to determine uncertainties related to the excited state contamination. Moreover, an effective mass analysis using Prony's method [57][58][59] and a lattice setup with displaced baryons [53], further reinforce the findings of two clearly separated energy levels as in Fig. 2 [53]. Next we calculate the energy difference between the ground state of the dibaryon (D 6b ) and the noninteracting two baryons (2 Ω bbb ) In  Scattering analysis:− To establish the existence of a state from these energy levels in terms of pole singularities in the Ω bbb Ω bbb S-wave scattering amplitudes across the complex Mandelstam s-plane, we use the generalized form of finite-volume formalism proposed by M. Lüscher [26]. For the scattering of two spin-3/2 particles in the S-wave leading to a total angular momentum and parity J P = 0 + , the phase shifts δ 0 (k) are related to the finite-volume energy spectrum via Lüscher's relation: Here, k is the momentum of Ω bbb in the center of momentum frame and is given by where ∆E is the energy differences listed in Table I, and M phys Ω bbb is the mass of Ω bbb in the continuum limit. The S-wave scattering amplitude is given by t = (cotδ 0 −i) −1 , and a pole in t related to a bound state happens when k cotδ 0 = − √ −k 2 . We parameterize k cotδ 0 = −1/a 0 , where a 0 is the scattering length. The scattering analysis is performed following the procedure outlined in Appendix B of Ref. [60], such that the best fit parameters are constrained to satisfy Eq. (5). To estimate the systematic uncertainties from the lattice cut-off effects, we perform several different fits involving different subsets of the four levels with k cotδ 0 parameterized either as a constant or as a constant plus a linear term in the lattice spacing. All of the fits indicate the existence of a deeply bound state. We find that the best fit corresponds to the one that considers all energy levels and incorporates the lattice spacing a dependence of the scattering length with a linear parameterization 0 − a/a [1] 0 . We present this as our main result, leading to a χ 2 /d.o.f = 0.7/2, with the following best fit parameters and binding energy In Figure 3, we present details of our main results. On the top pane, the analytically reconstructed finite-volume energy levels (black stars) from best fit parameters in Eq. (7) can be seen to be in agreement with the simulated energy levels (large symbols), indicating quality of fit. In the middle pane, we plot k cotδ 0 versus k 2 in units of the energy of the threshold. The orange dashed curve is the bound state constraint √ −k 2 and the red solid line is the fitted k cotδ 0 in the continuum limit. The crossing between these two curves, highlighted by the magenta symbol, is the bound state pole position in t. In the bottom pane, we present the continuum extrapolation of binding energy leading to the value in Eq. (8) compared with the simulated energy levels at the respective lattice spacings. The magenta symbol represents the binding energy in the continuum limit, with thick error representing the statistical and fit window error. The thin error includes the systematics related to excited state effects added in quadrature. Coulomb repulsion:− With two units of electric charge in the system, the effect of Coulomb repulsion on the binding energy of this dibaryon could be important. To gauge that, we perform an analysis, as in Ref. [24], and detail that in the Supplemental Material [53]. We model the strong interactions between two interacting Ω − bbb baryons with a quantum mechanical multi-Gaussian attractive potential, constrained to match the binding energy −81( +14 −16 ) MeV that we find in this work. In Fig. 4, we present the model potentials for strong and Coulombic interactions and also their combination, together with the radial probabilities of the ground state wave functions in the strong and combined potentials. Evidently, the Coulombic potential hardly affects the strong interac-    Coulomb (Ve), the parameterized strong potentials (Vs) and their sum are shown by the black, blue and red curves, respectively. The shaded region represents the variation of Vs with respect to its parameters. Ve is evaluated at a rms charge radius equal to the rms radius of the Vs ground state. The radial probability densities of the ground state wave-functions of the strong and combined potentials are shown by the dashed-dotted curves.
tion potential in the length scales where the ground state probabilities peak and infer that it serves only as a perturbation. The associated maximum change in binding energy is found to be between 5 and 10 MeV.
After addressing the systematic errors along with excited state contaminations [53] the final value of the dibaryon mass is determined by adding ∆E D 6b [−81( +14 −16 ) (14) MeV] with the two-baryon mass 2M Ω bbb .
Since the Ω bbb baryon mass is unknown we use its lattice extracted value. To this end, we perform continuum extrapolation of the energy splitting M Ω bbb (a) − 3 2 M 1S (a), and then add 3/2M phys 1S , with M phys 1S = 9445 MeV [61], to that. Thus we arrive at M phys Ω bbb = 14366(7)(9) MeV, which is consistent with other lattice results [62]. Using that, we obtain M phys Possible effects of Coulomb repulsion are included in the systematic errors. Error budget:− Finally we address the possible sources of errors in this calculation. We use a lattice setup with 2+1+1 flavored HISQ fermions where the gauge fields are Symanzik-improved at O(α s a 2 ), and the NRQCD Hamiltonian has improvement coefficients up to O(α s v 4 ). Such a lattice setup has shown to reproduce energy splittings in bottomonia with a uncertainty of about 6 MeV [53]. Note that here we are calculating the energy difference in which some of the systematics get reduced. For the dibaryon ground state in the finite volume, statistical, excited-state-contamination, and fit-window errors are the main sources of error. The energy levels are extracted using single exponential fits to the correlation functions from rigorously identified ground state plateau regions [53]. Correlated averages of various fitting intervals are considered to arrive at conservative fitting-window errors. Statistical and fit window errors are added in quadrature, and then convolved through the Lüscher's analysis and continuum extrapolation. The excited state contamination is determined from differences in the continuum limit estimates from the scattering analysis using results from different sink smearing procedures followed. However, it still would be worthwhile to investigate excited state uncertainties more precisely in future variational calculations. Other possible sources of errors are related to the continuum extrapolation fit forms, scale setting, quark mass tuning and electromagnetic corrections that together are found to be 12 MeV in such energy splittings, as detailed in the Supplemental Material [53]. Various errors are finally added in quadrature, yielding a total error of about 20% for the binding energy. Our results and inferences are robust up to the statistical and systematic uncertainties we have determined. Summary and Outlook:− In this Letter, using lattice QCD we present a first investigation of the dibaryons in which all six quarks have bottom flavor and find a deeply bound dibaryon (D 6b ≡ Ω bbb -Ω bbb ) in the 1 S 0 channel. Following Lüscher's formalism, we determine the relevant scattering amplitude, and after considering possible systematic uncertainties [53], we identify a bound state pole with a binding energy −81( +14 −16 )(14) MeV relative to the threshold 2M Ω bbb . The mass of D 6b dibaryon corresponding to this pole is found to be 28651( +16 −17 )(15) MeV. Although this dibaryon is expected to be compact, we find the Coulomb repulsion within this dibaryon acts only as a perturbation to the strong interactions and may shift the mass only by a few percent. The use of complementary measurements and analysis procedures in identifying the real ground state plateau ensure the robustness of our results. Our results provide intriguing evidence for the existence of the bound D 6b state, and it would surely motivate both phenomenological studies of its detection as well as follow-up lattice QCD studies investigating hardto-quantify excited-state uncertainties more precisely.
It is interesting to observe that the interactions between similar baryons using different procedures at the strange and charm quark masses are found to be very weak [10,12,14,24]. Note that a clear consensus on such systems with possible near threshold features requires complementary investigations of the same system with same high statistics ensembles but with different procedures. In comparison with the light and strange sectors, the binding energy of multiquark hadrons involving more than one bottom quark are predicted to be large [23, 28-31, 63, 64]. In this work we also observe the similar pattern in the Ω bbb Ω bbb channel. Taken together a common interesting pattern is emerging that the presence of more than one bottom quark enhances the binding in multihadron systems, which needs to be understood thoroughly including the quark mass dependence of scattering parameters.
Although a direct identification of D 6b dibaryon is a long way to go, our results on this heavy dibaryon, particularly because of its deep binding, will provide a major impetus in experimental searches for heavy quark exotics. Very much like the discovery of Ξ cc leading to predictions of various possible heavy multiquark systems [63], the discovery of doubly bottom baryons would be an important step in filling up the blanks higher up in the hadronic reaction cascade bringing prospects for discovering various bottom quark exotics, including D 6b . Given the recent excitements in the search for new heavy exotics [65][66][67] with multiple theoretical proposals and ideas [40][41][42][43], it is highly anticipated that substantial efforts, both on the theoretical as well as experimental fronts, would be steered and accelerated in this direction in the coming years.
using correlated χ 2 and maximum likelihood estimators to extract E 0 and W 0 . In Figure 5 we present such a result showing the projections of posterior probability distributions of the parameters E 0 and W demonstrating the reliability of the fits for the example of D 6b correlation functions in the finest ensemble. In order to quantify the uncertainties arising from the choice of fitting window (τ min , τ max ), we do the following. First choose a τ max as large as possible with a good signal-to-noise ratio. Then the τ min is varied over a range to determine the stability of E 0 estimate and a τ min value is chosen where a clear plateau is observed. A conservative estimate taking account of an uncertainty on this choice is arrived at using a correlated average over neighboring τ min values in the plateau. In Figure 6 we present the τ min dependence for all the fits along with the 1-σ statistical errors for the chosen fit window (blue bands), and the final estimate considering the uncertainty from the chosen fitting window (magenta bands). In both figures, we present the estimates for single baryons on the left and for dibaryons D 6b on the right. These estimates are then utilized to arrive at the energy differences in Eq. 4 and Table I

ERROR ANALYSIS
The main source of error in a lattice QCD calculation for a multi-hadron system arises from the rapid decrease of the signal-to-noise ratio in the correlation functions [36]. In heavy hadrons, this is somewhat mitigated due to the presence of heavy quarks. In this calculation since all the valence quarks are of bottom flavor and no chiral dynamics is involved, it is expected to have a relatively better signal-to-noise ratio than that of other dibaryons. Nevertheless, different systematics need to be addressed, particularly those arising from the contamination of excited states and from the lattice discretization, to arrive at a reliable estimate for the binding energy of the D 6b dibaryon. We discuss various relevant systematics involved in our calculation below.

Extraction of the ground state masses:
Coulomb gauge fixed wall sources are utilized for the quark fields, which has been widely used over many years for calculations involving NRQCD and is known to produce good ground state plateau at a large source-sink separation. We also average the correlation functions over multiple source time-slices to improve the statistical uncertainties. In addition to this, we follow the procedure outlined in the previous section to include fit-window uncertainties, and arrive at the final energies and energy differences presented in the main text. For J = 3/2 Ω bbb baryons, we utilize the most symmetric operator ( 1 H + ) with rows given in Table II, also expressed in Eq. (7) of Ref. [10]. It was observed from studies which had previously utilized this operator for the studies of ∆ [69,70], Ω [10,70], Ω ccc [71] and Ω bbb [72] baryons, that the ground state has the largest overlap with this operator, and the ground states for these singly-flavored baryons are best determined with this operator. It is also observed that the radial excitations for the decuplet baryons are very high in energy arising from higher partial waves. For Ω bbb , the first radial excitation is observed to be > 400 MeV above the ground state Ω bbb from the lattice calculation in Ref. [72]. Any reminiscent effects from those should be reflected as significantly different approach to the energy plateau in the effective energy plots and the t min dependence plots in different lattice QCD ensembles. We observe in all our ensembles, which vary in lattice spacing and the volume, the signal plateauing commences from approximately the same physical temporal extent ( 2.5 fm ).
For the dibaryon, we utilize the S-wave projected twobaryon interpolating operator expressed in Eq. (2) in the main text, as was also utilized in Ref. [10]. In our case, the action being non-relativistic we are limited to the two-baryon operator built purely out of 1 H + baryon operators, which is also observed to be the best choice for the ground state determination. We also note that the large time approach to the energy plateau in this case also occurs around the same physical temporal extent ( 2 fm ) across all our ensembles, which confides the reliability of our energy estimates. Nevertheless, to ensure it further, we perform a set of additional calculations, both for single and dibaryon systems, in various setups to arrive at a conservative estimate on possible excited state contaminations in our results. We describe that below.

Excited state contamination:
The use of wall smearing at the quark source and no smearing at the quark sink is an asymmetric setup in building two point correlation functions. This results in an unconventional rising-from-below behavior of the effective energies, as a result of competing overlap factors with different signatures for different states having same quantum numbers. Consequently, there could be low lying plateaus at early times that mimics the real ground state plateau. To remove this complication, we perform a set of calculations with different source-sink setup, study the asymptotic behavior in search of a universal estimate for the ground state energy level along with a reliable estimate for possible uncertainty from excited state contaminations. To this end, in addition to the previous wall-source and point sink setup, we perform various different exercises as described below.
A. Wall-source and Gaussian-smeared sink: In the first additional setup, we use a wall source along with a Gaussian-smeared sink which have been extensively used in heavy hadron calculations [54,55]. Choosing a suitable Gaussian-width one can achieve a reliable ground state plateauing in the correlation functions. We choose the Gaussian width such that the effective energies feature a conventional falling-from-above behavior and yet the statistical noise in the correlation functions do not wash the signal away. We observe that the best value of Gaussian width across all the ensembles we study is ∼ 0.2 fm. Corresponding results obtained on our finest lattice ensemble are shown in Fig. 7. It is clear from this figure that results from the wall-source point-sink (w-p) correlators and the wall-source Gaussian-smeared sink (w-gs) correlators are quite consistent with each other. The effective mass obtained from the w-gs correlators for the dibaryon operator always stays below than that of the corresponding non-interacting two-baryon correlators. We fit these w-gs correlators with one exponential and the fitted results with errorbars are shown by the horizontal bands.

B. Wall-source and spherical-box-smeared sink:
We employ a second setup with a spherical-box-smeared sink. This procedure was utilized in Ref. [56] for doubly heavy tetraquark calculations and have been found to be effective in avoiding the rising behavior in the effective energies, and in getting an early plateau at the correct ground state energy. We have varied the radius (r) of the spherical-box and have tuned its value such that the effective energies show a falling behavior and yet retains a sufficiently good signal-to-noise ratio. In Fig. 8, we plot the effective mass obtained on our finest ensemble with r ∼ 0.34 fm. It is evident that the effective mass falls from above and its asymptotic behavior is consistent with that of the wall-point and wall-Gaussian-smeared setups that we have discussed. We fit these wall-box (wb) correlators and the results with 1-σ and fitting-window errors are shown by the horizontal bands in Fig. 8.
C. Effective mass using Prony's method: Further to the above consistency checks, we have performed another complementary analysis known as the "Prony's-method" [57][58][59] to find the reliability of the ground state plateau and the extent of excited state contaminations in the ground state energy estimation. Although this method was found to be unstable with smaller statistics, it was shown to be quite effective to get a reliable ground state effective mass with the high statistics correlation functions [73]. It was also found to produce the energy levels favorably to that obtained through the variational approach [74]. In Fig. 9 we show the effective mass obtained for the wall-point correlators using this method with two exponentials where the solutions are numerically stable (solving Eq. (16) and (17) of Ref. [73] numerically). We also find that the solutions are often unstable with large errors for arbitrary choices of n and q 1 values of Eq. (16) and (17) of Ref. [73]. Using suitable choices of n and q 1 we find stable solu- tions and show that in Fig. 9. The effective masses are shifted towards the right as per the choice of n and q 1 . It is clear from this plot that the effective masses obtained using one exponential from the wall-point correlators are consistent with that obtained using Prony's method using two exponentials. This provides another consistency check of the findings using wall-point correlators. FIG. 9. Effective mass obtained from using Prony's method [57][58][59]. The band shows the fitted results with its error determined using the wall-source point-sink data.

D. Dibaryon operators with displaced baryons at sink:
As an additional exercise, we investigate the effect of using displaced baryons (B) for two-baryon operators O(r = 0; x) = B(x)B(x) at the sink, and also compare the observations with that obtained with similar setup for the well established deeply bound system of doubly heavy tetraquarks. To this end, we displace the two baryons in the two-baryon operator such that O(r; (x 1 + x 2 )/2) = B(x 1 )B(x 2 ), where the displacement r = |x 1 − x 2 | (here r is symmetrized with respect to all three spatial directions). We observe that the ground state energy estimate from such an operator is consistent with that of the local two-baryon operators until below r ∼ 0.25 fm. Note that this also coincides with the chosen width of the Gaussian smearing and closer to the boundary of the spherical-box-smearing that yields a conventional behavior of falling effective energies. In Fig. 10, we present these results obtained on the finest ensemble with different values of r. To understand the significance of this result, we then carried out the same exercise for the case of wellstudied doubly-bottomed tetraquarks (bbud) with the same asymmetric wall-source and point-sink setup. Note that for this case, the energy level corresponding to the four-quark operator of two-meson type lies below than that of the threshold level obtained from the two noninteracting two mesons. As in the case of dibaryon, we vary the two sink points (x 1 and x 2 ) of two-meson operators (b(x 1 )u(x 1 )b(x 2 )d(x 2 )). Interestingly, we find strikingly similar results as above, where up to a certain displacement (r), the effective mass of two-meson operators with one (x 1 = x 2 ) and two sink points (x 1 = x 2 ) are found to be consistent with each other. Thebbud systems have been studied recently with variational methods by multiple lattice-QCD groups with asymmetric [29,30], box-smeared [56] as well as smeared point-to-all correlators [31], and results on the bindings obtained by those different lattice calculations consistently found a deeply bound state. Since the response to the binding with respect to the displacement of two sink points are strikingly similar both for the dibaryon studied here and for b(x 1 )u(x 1 )b(x 2 )d(x 2 ) tetraquarks, andbbud was found to be deeply bound by multiple studies, we believe our result is robust up to the statistical and systematic uncertainties that we have determined. This finding on the existence of a deeply bound heavy dibaryon calls for fur-ther lattice calculations on heavy dibaryons particularly using multi-operator variational approaches to quantify the systematics related to the hard-to-quantify excited state effects more precisely, which have been found to affect the results obtained for light dibaryons that employed asymmetric correlators [16,17].
In summary, various procedures followed above led us to conclude that the results obtained using asymmetric wall-point setup is robust and the effect of the excited state is minimum as long as the fitting window corresponds to the real plateau from the ground state energy. We find that the ground state plateau saturates 2 fm for both the baryon and dibaryon correlators across all the ensembles. We also observe that the dibaryon correlator overlaps with the ground state maximally when the smearing size is about 0.2 fm. This is also in line with the results from two-baryons operators with displaced baryons at the sink. This observation on the smearing width (∼ 0.2 fm) perhaps is indicating that the observed dibaryon could be a compact state. To account for the effect of contaminations from the excited state we include conservative errors determined based on the difference in energy estimates obtained from different calculations, as discussed above.

Continuum extrapolation:
We employ a set of lattice QCD ensembles in which gauge fields are Symanzik-improved at O(α s a 2 ) and include the effect of u, d, s and c quark vacuum polarization generated with the highly improved staggered quark action [45]. Quark propagators are generated with NRQCD action with improvement coefficients up to O(α s v 4 ). The lattice spacing dependence of the energy differences (in Table I) could be nontrivial. Similar to the approach made in Ref. [16], we account for this by parameterizing kcotδ 0 , that enter the scattering analysis in Eq. 5, with different forms and perform fits with different sets of energy levels determined from the simulation. Choosing the linear parameterization k cotδ 0 = −1/a [0] 0 − a/a [1] 0 that best describes the entire data, we find the total uncertainties arising from statistics, fitting window and continuum extrapolation to be ∼18% of the binding energy from the continuum extrapolation. We find that choosing other forms of continuum extrapolation for the scattering length −1/a 0 leads to a change of at most 8 MeV in the binding energy, which we quantify as the uncertainty arising from the discretization error.

Scale setting:
Scale settings through r 1 parameter [45] and Wilson-flow were found to be consistent [45] for these lattice ensembles. Systematics with the scale settings further gets reduced in the estimation of energy differences (Eq. 4), and as in Ref. [23,44] we find it to be maximum of about 3 MeV.

Quark mass tuning:
We tune the bottom quark mass employing the Fermilab method of heavy quarks [50]. Here, we equate the lattice extracted spin average 1S bottomonia kinetic mass, 1 4 [3M Υ + M η b ] kin , with its physical value. We perform this tuning corresponding to the central value of the chosen scale and also at its error values. We calculate E D 6b for each of the tuned masses and include the variation as the estimation of error due to quark mass tuning. We find it to be less than 2 MeV.
With the above mentioned lattice setup we find the hyperfine splitting in 1S bottomonia, a benchmark observable for the evaluation of the goodness of lattice calculations with bottom quarks, is quite consistent with its experimental value, as demonstrated in Figure 11. The continuum value (green star) is obtained taking the average of estimates from all ensembles and the error (green band) is estimated as a weighted average with respect to the lattice spacings. Continuum extrapolation with the linear as well as and quadratic forms in lattice spacing are also shown by the orange and blue stars respectively with the same color bands for their 1-σ errors. Together with possible other systematics, that we are discussing here, we estimate its value to be 62.6(3)(5) MeV.

Electromagnetism:
The dibaryon investigated here has two units of electric charge which can affect its binding. To gauge that, we perform the following analysis as in Ref. [24]. First, we model the strong interactions between two interacting Ω − bbb baryons with a quantum mechanical multi-Gaussian attractive potential V s [24], constrained to match the binding energy −81( +14 −16 ) MeV that we find in this work. Next, we assume the form of the Coulomb potential (V e ) of Ω − bbb to be similar to that of Ω ++ ccc , except the total electric charge is −2. We present a comparison of the strengths of these potentials as a function of the radial distance in Figure 4 of main text, with the rootmean-square (rms) charge radius r d chosen as the rms radius of the ground state of V s . Next, we solve the energy eigenvalue problem with the effective potential (V ef f = V s + V e ) and determine the scattering length a e+s 0 and effective range r e+s , following the procedure discussed in Ref. [24]. The radial probability densities of the ground state wave-functions (dashed-dotted curves) corresponding to V s and V ef f are shown in Figure 4 of the main text. It is evident that the Coulomb repulsion serves only as a perturbation and hence does not change the binding energy of D 6b in any significant way. We also vary r d and find that the effect of Coulomb repulsion is largely perturbative and binding may reduce at most by 10 MeV even when r d is chosen to be unphysically low as 0.01 fm. We present 1/a e+s 0 for V ef f as a function of the Coulomb interaction strength α e in Figure 12. Note that 1/a e+s 0 remains to be very much positive even at α e = α e phys , confirming that D 6b remains to be a deeply bound state even in the presence of Coulomb repulsion, with a total binding energy of about −75 MeV. For heavy baryons, the possible systematics due to other electromagnetic corrections was found to be 3 MeV [75]. Keeping that in mind as the source of other electromagnetic effects beside the Coulomb repulsion, we take a conservative estimate of 8 MeV corrections for the binding energy (by adding the average of Coulomb repulsion with the above mentioned 3 MeV in quadrature).
No chiral extrapolation is necessary for D 6b . For heavier dibaryons the unphysical sea quark mass effects are expected to be within a percent level [76][77][78], and particularly for D 6b , it would be negligibly small. In Table  III we summarize the error-budget estimate where above mentioned systematics are added in quadrature.

Source
Error (MeV) Statistical + Fit-window + +16 −14 Excited states  8  Discretization  8  Scale setting  3  m b tuning  2  Electromagnetism  8  Total systematics  12   TABLE III. Error budget in the calculation of the binding energy ∆ED 6b . The total systematics quoted above includes those from the discretization, scale setting, bottom quark mass tuning and electromagnetic effects.