Deuteron-like heavy dibaryons from Lattice QCD

We report the first lattice quantum chromodynamics (QCD) study of deuteron($np$)-like dibaryons with heavy quark flavours. These include particles with following dibaryon structures and valence quark contents: $\Sigma_c\Xi_{cc} (uucucc)$, $\Omega_c\Omega_{cc} (sscscc)$, $\Sigma_b\Xi_{bb} (uububb)$, $\Omega_b\Omega_{bb} (ssbsbb)$ and $\Omega_{ccb}\Omega_{cbb} (ccbcbb)$, and with spin ($J$)-parity ($P$), $J^{P} \equiv 1^{+}$. Using a state-of-the art lattice QCD calculation, after controlling relevant systematic errors, we unambiguously find that the ground state masses of dibaryons $\Omega_c\Omega_{cc} (sscscc)$, $\Omega_b\Omega_{bb} (ssbsbb)$ and $\Omega_{ccb}\Omega_{cbb} (ccbcbb)$ are below their respective two-baryon thresholds, suggesting the presence of bound states which are stable under strong and electromagnetic interactions. We also predict their masses precisely. For dibaryons $\Sigma_c\Xi_{cc} (uucucc)$, and $\Sigma_b\Xi_{bb} (uububb)$, we could not reach to a definitive conclusion about the presence of any bound state due to large systematics associated with these states. We also find that the binding of these dibaryons becomes stronger as they become heavier in mass. This study also opens up the possibility of the existence of many other exotic nuclei, which can be formed through the fusion of heavy baryons, similar to the formation of nuclei of elements in the Periodic Table.

A deuteron is a bound state of two baryons, a proton and a neutron, and is made of six light valence quarks.In the early Universe, deuterons were created and their stability is responsible for the creation of other elements.Interestingly, the strong interactions between quarks, which bring stability to deuterons, also allow various other six-quark combinations leading to the possible formation of many other dibaryons.However, no such dibaryons, though speculated about many times [1][2][3][4][5][6][7], have been observed yet.Using a state-of-the-art first principles calculation of lattice quantum chromodynamics (QCD), here we report, for the first time, a definite prediction of the existence of other deuteron-like spin-1 dibaryons.We also predict their masses precisely.These new subatomic particles could either be made of six heavy quarks (charm and bottom) or heavy and strange quarks.Such dibaryons are stable against strong and electromagnetic decays, but, unlike the deuteron, they can decay through weak interactions.We also find that the stability of such dibaryons increases as they become heavier.We expect that prediction from this work will aid in discovering these new subatomic particles at experimental facilities, such as the Large Hadron Collider.In fact this study opens up the possibility of the existence of many other exotic nuclei, which can be formed through the fusion of heavy baryons, similar to the formation of nuclei of elements in the Periodic Table .Formation of these hadrons perhaps also enhances the possibility of a quarklevel analogue of nuclear fusion as discussed recently [8].Further study of these exotic states can also provide information on the strong interaction dynamics at multiple scales.
The particular dibaryons (D) 1 that we investigate are heavy quark analogues of deuteron (np).They have the spin-(J)-parity (P ) quantum numbers: J P = 1 + , with following dibaryon configurations: Here, Σ q , Ξ qq , Ω qq , Ω q1q2q2 's are heavy baryons with the usual nomenclature of the Particle Data Group [9], and u, s, c, b inside brackets are various quark flavours.We find that D cs , D bs and D bc are stable against strong and electromagnetic decays, and therefore, it should be possible to discover them in experimental facilities.However, for D cu and D bu , we find the ground state masses are consistent with their respective two-baryon thresholds, and therefore our results are not currently precise enough to reach a definitive conclusion on their stability.Incidentally, only recently tetra- [10][11][12] and pentaquark [13,14] states have been discovered and those are made of heavy quarks.Recent model and lattice QCD studies also suggest the existence of other heavy tetraquark hadrons [15][16][17][18][19].It is thus natural to search for six-quark states made of heavy quarks and our predictions for bound heavy dibaryons provide crucial information for their discovery.
Dibaryons have been studied through various models over the years [1][2][3][4][5][6][7]20], However, it is quite crucial to have first principles lattice QCD based studies on these states to predict their masses and to understand their structures.In fact a few lattice QCD studies have already been carried out [21][22][23].However, those are mainly focused on light quarks with spin-0 states [21][22][23][24], studies of deuteron in Refs.[25,26] as well as studies of and name them as Dq 1 q 2 , which are made of two baryons with valence quark contents (q 1 q 1 q 2 ) and (q 1 q 2 q 2 ).In this notation the deuteron is D ud ≡ np(uududd).Such dibaryons may be called as D-dibaryons.
spin-2 states [27].This work is the first lattice study on heavy dibaryons.A lattice dibaryon calculation with heavy quarks is advantageous over the light counterparts in two ways.The two point correlators are less noisy in comparison with light dibaryons and secondly the signalto-noise ratio of the two point correlators is far better in comparison to the light quark calculations.Both of these provide an added advantage in our calculation and help in reliable extraction of binding energies of these heavy dibaryons.
The lattice set up that we utilize here is similar to that was used in Refs.[18,28,29].Below we elaborate it further for the sake of completeness.A. Lattice ensembles: Three sets of dynamical 2+1+1 flavours (u/d, s, c) lattice ensembles, generated by MILC collaboration [30] with HISQ fermion action [31] are employed for this study.Lattice spacings, using r 1 parameter, for these ensembles are measured to be 0.1207 (11) 0.0888(8) and 0.0582(5) fm, respectively [30], which are also found to be consistent with scales obtained through Wilson flow [32].B. Quark actions: In the valence sector, from light to charm quarks, we utilize the overlap action which has exact chiral symmetry at finite lattice spacings [33][34][35] and no O n (ma), n = 1, 3, 5, • • • errors.A non-relativistic QCD (NRQCD) formulation [36] is adapted for the bottom quark with an O(α s v 4 ) improved Hamiltonian with non-perturbatively tuned improvement coefficients [37].For reliable extraction of the ground states we employ a wall source.C. Quark mass tuning: We tune both the charm and bottom quark masses using the Fermilab prescription for heavy quarks [38] in which we equate the latticeextracted spin-averaged kinetic masses of the 1S quarkonia states with their physical values [9].Following Ref. [39] we tune the strange quark mass to its physical value.D. Dibaryon interpolators: The dibaryon interpolating operator with spin (J = 1) and antisymmetric in flavours (q, Q) ∈ (s, c, b) is constructed as: where Ω qqQ and Ω QQq are spin-1/2 baryons defined as, (Ω qqQ Here Latin letters indicate color, while Greek letters indicate the spinor degree of freedom.The various deuteron analogues with appropriate flavour antisymmetry are listed in Table I. In a lattice QCD formulation [40], the mass (m) of a particle is extracted in two steps: first by calculating the two-point Euclidean time (τ ) correlator ( C(τ ) ) of the interpolating source and sink operators, over many gauge configurations, and then extracting the exponent (m) by fitting the exponential decay ( C(τ ) ∼ e −mτ ) of that correlator at large Euclidean time.Following the dibaryon structure of deuteron (d = 1 √ 2 (pn − np)), the spin-1 flavour-antisymmetric dibaryon interpolating fields are constructed with the appropriate spin projection of two individual spin-1/2 baryons, as shown in Table I.We then calculate the two-point correlators with those interpolating fields and from the exponential decay of these correlators we calculate the lowest energy states, i.e., the ground state masses of each dibaryons.Masses of individual baryons are also calculated similarly.Detail of the calculation method is provided in the supplementary information (SI).

DQq
Interpolating fields The next step is to find the relative energy levels of the ground state dibaryons with respect to their two-baryon thresholds.Since baryon number is a conserved quantity in the Standard Model of particle physics, these spin-1 dibaryons can have only two strong decay channels: either to two spin-1/2 baryons or to two spin-3/2 baryons.Interestingly, we find that the ordering of masses for these two combinations is different with the charm and the bottom quarks.While the sum of masses of two spin-1/2 charmed baryons (Σ c and Ξ cc or Ω c and Ω cc ) are smaller than that of two spin-3/2 charmed baryons (Ω qqq,q=u,s and Ω ccc ), for the bottom quark this trend is opposite, These observations are consistent with known experimental results and lattice determination of baryon masses [28,29,41].After computing the dibaryon masses we then calculate their mass differences from both the spin-1/2 and spin-3/2 two-baryon thresholds.When the mass difference of a dibaryon from its closest threshold is negative and the finite volume effects are convincingly small, that particular dibaryon is likely to be a bound state.
We compute the aforementioned mass differences at multiple lattice spacings within our lattice set-up.Our results are shown in Figure 1.Mass differences from spin-1/2 and spin-3/2 thresholds are shown by red and green colours respectively.The upper plot is for the dibaryon D bc (Ω ccb Ω cbb ), middle one is for D bs (Ω ssb Ω sbb ) and the bottom one is for D cs (Ω ssc Ω scc ).It should be noted  Mass differences between various spin-1 heavy dibaryons (Dq 1 q 2 ) from their two-baryon threshold states.Red points represent mass differences when the threshold states are with spin-1/2 baryons and the green points are the same with spin-3/2 baryons.Results are shown at three lattice spacings and at the continuum limit with shaded bands as one sigma error.that these dibaryon energy levels are already computed at their physical quark masses and therefore only require a continuum extrapolation and finite volume corrections.For the continuum extrapolation, the energy levels are computed at three lattice spacings, indicated by different marker styles, and we use a linear in a 2 ansatz as well as with an a 2 ln(a) term.To be noted that use of overlap action, which has exact chiral symmetry on lattice, ensures no O n (ma), n = 1, 3, 5, • • • errors, which then assures that higher order discretization errors are smaller, particularly at the finest lattice.In addition, since we calculate mass splittings, rather than masses, errors due to finite lattices are in good control.The errors shown include both statistical as well as systematic errors (see SI for details) which are then added in quadrature to get the final errors.
From these figures it is quite apparent that the ground state masses of D bc , D bs and D cs lie below their respective closest two baryon thresholds by about 4, 2 and 3 sigma errors, respectively.Next, it is natural to ask if these mass differences obtained at finite volume lattices (L 3 ) are indeed the physical binding energies that hold these dibaryons together from decaying to individual baryons.
To sort this out one needs to study finite volume correction of these lattice computed energy levels.Fortunately for multi-hadron systems this has already been worked out [42][43][44] and the finite volume corrections, ∆ F V , for such systems was found to be ∼ O(e −k∞L )/L, where k ∞ = (m 1 + m 2 )B ∞ , m 1 , m 2 being the masses of threshold states and B ∞ the infinite volume binding energy.Since here m 1 , m 2 are masses of two baryons with multiple heavy quarks, ∆ F V receives a large suppression even when B ∞ is of a few MeV.This assures that the binding energies for these dibaryons will be very close to the extracted mass differences that we showed in Figure 1.Since we use the strange, charm and bottom quark masses already at their physical values, the extracted mass differences are indeed the physical binding energies of these dibaryons.
For dibaryons, D cu and D bu we first perform chiral extrapolations with a constant plus a term linear in m 2 π .Due to the presence of light quarks, the signal to noise ratios in the correlation functions of these states are rather poor.Moreover for D bu , one decay product is the resonance state ∆ uuu , which needs to be addressed with adequate finite volume study [45].With the current lattice set up, it is therefore difficult for us to make a precise conclusive statement about the stability of these two dibaryons.In future such a study can be carried out with the availability of more computing resources.However, following the example of deuteron, if we assume they are also bound, our results suggest that their binding energies will be much smaller compared to other three dibaryons mentioned above.The final values of dibaryon masses are calculated by adding the known two-baryon threshold masses to the continuum value of binding energies that we extracted.For masses of yet-to-be discovered baryons we use their lattice values as calculated in this work and in Refs.[28,29,41].We also use subtraction method [28,29,41,46] utilizing spin-average value of 1S charmonia and bottomonia and find results consistent with above mass estimates.The final values of dibaryon masses are shown in the third column of Table II.In Fig. 2 we show the relative energy levels of TABLE II.Energy differences between spin-1 heavy dibaryons and their two-baryon thresholds.The third column shows the predictions of masses for the stable dibaryons.The bottom part of the table is for dibaryons has unphysical quark, q, with mass mq varying over a wide range.
It is quite apparent that the dibaryons D bc , D bs and D cs lie below their respective closest thresholds (horizontal zero line) by about 4, 2 and 3 sigma errors, respectively, while for the other two, D bu , and D bu , we could not reach to a definitive conclusion due to large errors.It may be noted here that we have obtained our results from simulations with a finite number of statistical ensembles, and on space-time grids with finite lattice spacings and finite volume, and hence these are associated with both statistical and systematic errors.Our final results are obtained after carefully addressing those errors and we elaborate that in SI.
Binding energy of the spin-1 dibaryon D bq where we kept the bottom quark fixed and vary the other quark, q, from light to bottom quarks.Along the x-axis we show the ratio of the pseudoscalar mass at (mπ q ) to that at the bottom quark (mη b ).As the mass of the quark q increases (mq = mu to m b ), the dibaryon binding energy also increases.
Interestingly, our results point out that the strong binding energies of these heavy dibaryons increase as they become heavier in mass.To confirm this pattern we also calculate similar dibaryons, D bq , with bottom quark at its physical value while varying the other quark mass (m q ) between charm to bottom quark masses.Of course, such quarks are unphysical and so are these dibaryons, but since they obey the same strong dynamics they can clarify the pattern of stability.In Table II, at the bottom part, we tabulate the binding energies of these unphysical dibaryons (these results are obtained only at our finest lattice).To depict it more clearly, in Fig. 3 we plot these binding energies for dibaryons D bq (Ω bbb Ω qqq ) with q varying from light to all the way to the bottom quark.It clearly shows that the spin-1 dibaryons become more stable when they are heavier.We even calculate such dibaryons with much higher quark masses (which are of course unphysical) and observe that their binding energies decrease as q > b and they vanish at very large quark masses.This suggests that the combination D bq (Ω bbb Ω qqq,q=b ) has the maximum binding.
The stability of these dibaryons against the strong and electromagnetic decays opens up the possibility of finding these and even more complicated higher nuclei with many heavy quarks, similar to nuclei of various elements in our periodic table.Similar to the role of deuteron in the nuclear fusion cycle for the creation of elements, these dibaryons can help the fusion processes of heavy baryons to produce nuclei with heavy quarks.Such nuclei can be studied theoretically in future with adequate computational resources and it may well be possible to discover them in future higher energy heavy ion facilities.It will also be interesting to see if the formation of such states has any implication in cosmology.
Formation of such hadrons also enhances the speculation on the possibility of a quark-level analogue of nuclear fusion which was discussed recently in Ref. [8].For example, formation of D bs through fusion of Ω bb and Ω b , as well as through fusion of Ω bbb and Ω sss , are highly exothermic with the release of energy about 300 and 30 MeV, respectively.We also find that reactions such as qqq , q ≡ c, s, u/d, and qqc , q ≡ s, u/d, are highly exothermic (here we represent B J q 1 q 2 q 3 as a baryon with spin J and quark contents q 1 , q 2 and q 3 ).
Since these dibaryons involve quarks with masses over a wide range, studies of their properties will be helpful to understand the strong dynamics at multiple scales.It is believed that the effective tensor interaction provides stability to deuterons [47] and hence these dibaryons will be an ideal laboratory to study the origin of a such force.Due to the presence of multiple heavy quarks they will decay via various possible weak decay processes.For example, D cb can decay through b → c, b → s and c → s to various light baryons and multiple mesons which can interfere among themselves.Detailed analysis of these multiple ways of decay may also be an ideal place to study the hadronic-interference in weak decay processes.of spin average masses to that as We find results from both methods are consistent with each other.

Error analysis:
Results obtained at finite volume lattices have both statistical as well as systematic errors.We address each of those below which are similar to the error analysis in Refs.[28,29,41].Statistical error: We use wall source to obtain dibaryon correlators and it helps to obtain better ground state plateau.We find a statistical uncertainty of 10 MeV while calculating mass difference for D bc .Fitting window error: With stable plateau we find uncertainty due to different fitting windows to be maximum of 4 MeV.Discretization error: To obtain reliable results for these dibaryons with heavy quarks, a crucial issue is the control of discretization errors.The following three methods help us to reduce discretization errors: i) continuum extrapolation, as mentioned above, of the results obtained at three lattice spacings, the finest one being at 0.06 fm, ii) use of overlap action and iii) extraction of mass differences, which has less discretization error than masses.As mentioned previously, the first two terms that enters in the continuum extrapolation are a 2 and a 2 ln(a).With three data points we fit these two forms separately and also together in a constrained fitting.Differences in central values at the continuum limit obtained with these different fittings are included in the discretization errors.After continuum extrapolation, we find the maximum discretization would be less than 4 MeV for the mass difference in D bc .Scale setting error: On these set of ensembles we have also determined scales by measuring the Ω sss baryon mass and those were found to be consistent with the determinations using r 1 parameter [49].Measurement of scale with Wilson flow [32] was also found to be consistent with the scale used here.The scale setting uncertainty in the mass difference (Eq.( 2)) for D cs is found to be ∼ 4 MeV.Quark mass tuning error: We use the Fermilab method of heavy quarks [38] to tune the charm and bottom quark masses.In this method, we calculate the kinetic masses, instead of pole masses, of the spin average 1S quarkonia and equate those with their experimental values.We perform this process corresponding to the central value of the scale and also with central ± error values.For each of the tuned masses, we calculate hadron masses and include the variation as the errors due to quark mass tuning.Our estimate for the mass tuning errors in the energy difference for D bc due to the charm and bottom quarks are about 2-3 MeV.Finite size effects: The finite volume corrections, ∆ F V , for these multi-hadron systems were found to be ∼ O(e −k∞L )/L, where k ∞ = (m 1 + m 2 )B ∞ , m 1 , m 2 being the masses of threshold states and B ∞ infinite volume binding energy.Since here m 1 , m 2 are masses of two heavy baryons, ∆ F V would be very small even if B ∞ is of a few MeV.Following our works on heavy baryons in Ref. [28] and similar works in Ref. [41], we include an uncertainty of 2 MeV from finite volume effects while calculating the mass difference for D bc .However, the fine volume effects for dibaryons with light-quarks would be bigger and one needs to perform a dedicated finite volume study for that [45].
Other sources: For dibaryons D cu and D bu , we had to carry out chiral extrapolations and we perform that with constant plus m 2 π terms.However, since the pion masses are not so light, particularly on the fine lattice, it is difficult to control the systematic associated with the chiral extrapolation.The dibaryon D bu also involves its decay to ∆-resonance which needs to be treated with finite volume study and that is beyond the scope of this work.We have thus included large errors associated with it coming from chiral and possible finite volume corrections which in turn limits the precision of our predictions for these baryons.In future with the availability of more computer resources this limitation could be addressed.No chiral extrapolation is necessary for D cs , D bs and D bc for which we have made conclusive prediction.
The unphysical sea quark mass effects are expected to be within a percent level for heavier dibaryons with no effective valence light quark content [39,46,50].This effect however would be larger for D cu and D bu dibaryons.Errors due to mixed action effects are found to be small within this lattice set up [51] and are expected to vanish in the continuum limit.Errors from electromagnetism are expected to be within 3 MeV for baryons [52] and should be similar for these dibaryons.
Adding all these systematic errors in quadrature, we found that for the energy difference (Eq.2), i.e., for the binding energy, in the case of D bc is less than 8 MeV.For other dibaryons we also estimated those similarly and included those in table II.We find that systematic errors for dibaryons D cu and D bu are too big to reach a definitive conclusion about their bindings.
Below we summarize the error budget for D bc dibaryon.

D
FIG. 1.Mass differences between various spin-1 heavy dibaryons (Dq 1 q 2 ) from their two-baryon threshold states.Red points represent mass differences when the threshold states are with spin-1/2 baryons and the green points are the same with spin-3/2 baryons.Results are shown at three lattice spacings and at the continuum limit with shaded bands as one sigma error.

2 FIG. 2 .
FIG.2.Relative energy levels of the spin-1 dibaryons and their respective two-baryon non-interacting strong-decay threshold states.The horizontal line is for the closest threshold (normalized to zero) while 1/2 and 3/2 signify the spin of the constituent single baryons of the two-baryon threshold.

TABLE I
. Structures of spin-1 heavy dibaryons that we study in this work.

TABLE III .
Error budget in the calculation of energy splittings for D bc dibaryon.