Masses of the conjectured H-dibaryon at different temperatures

Quantum chromodynamics (QCD) describes the dynamics of quarks and gluons. It underlies all of nuclear physics from hadronic mass spectrum to the phase transition of hadronic matter to quark-gluon plasma (QGP). Because of the nature of the strong interaction of QCD at low energy scale, the perturbative method cannot be applied to explain those low energy phenomena of nuclear physics. Fortunately, lattice QCD which is based on first principles can be employed to make precise predictions for hadronic quantities, especially, for those phenomena which are difficult to explore in the laboratory, for example, mass spectrum of baryon at different high temperatures.

In 1976, by using the bag model, Jaffe predicted a flavour-singlet state ($uuddss$) with quantum number $I(J^{P})=0(0^{+})$ which is called H-dibaryon Jaffe:1976yi . In contrast with the only known stable dibaryon (deuteron) whose binding energy is about $2.2\ {\rm MeV}$, Jaffe predicted that the binding energy of H-dibaryon is about $80\ {\rm MeV}$ below the $\Lambda\Lambda$ threshold $2230\ {\rm MeV}$ which means that H-dibaryon is a deeply bound state.

Unlike mesons and baryons, this exotic hadron may be relevant to the hypernuclei, and to the strange matter that could exist in the core of neutron stars. Moreover, it is a potential candidate for dark matter Farrar:2017eqq . As a consequence, this prediction triggered a vigorous search for such a state, both experimentally Aoki:1991ip ; KEK-PSE224:1998trj ; Takahashi:2001nm ; Yoon:2007aq ; KEKE176:2009jzw ; Nakazawa:2010zza ; Belle:2013sba ; BaBar:2018hpv ; Ekawa:2018oqt and theoretically Iwasaki:1987db ; Luo:2011ar ; Luo:2007zzb ; Mackenzie:1985vv ; Pochinsky:1998zi ; Wetzorke:1999rt ; Wetzorke:2002mx ; Beane:2009py ; NPLQCD:2010ocs ; Beane:2011zpa ; NPLQCD:2011naw ; NPLQCD:2012mex ; Inoue:2010hs ; Inoue:2010es ; Inoue:2011ai ; HALQCD:2019wsz ; Sasaki:2016gpc ; HALQCD:2018lur ; Francis:2018qch ; Green:2021qol .

The observation of double hypernuclei, ${}^{A}_{\Lambda\Lambda}Z$, is very important in connection with the existence of H-dibaryon Aoki:1991ip . If the mass of H-dibaryon $m_{H}$ is much smaller than the mass of double $\Lambda$ hyperon $2\,m_{\Lambda}$, a double hypernuclei may decay into H-dibaryon and a residual nucleus by strong interation. In such case, the branching ratio for the decay of double hypernuclei through weak interaction is very small, practically, cannot be observed Aoki:1991ip . So the observation of the weak decay of double hypernuclei will put limitation on the mass of H-dibaryon $m_{H}>2\,m_{\Lambda}-B_{\Lambda\Lambda}$ with $B_{\Lambda\Lambda}$ being the binding energy of $\Lambda\Lambda$ hyperons.

The experiments Aoki:1991ip ; KEK-PSE224:1998trj ; Takahashi:2001nm ; Yoon:2007aq ; KEKE176:2009jzw ; Nakazawa:2010zza ; Ekawa:2018oqt investigated the nuclear capture of $\Xi^{-}$ at rest produced in $(K^{-},K^{+})$ reaction, and observed the sequential weak decay of double hypernuclei, then measured the binding energy and interaction energy of $\Lambda\Lambda$ Aoki:1991ip ; Takahashi:2001nm ; Nakazawa:2010zza ; Ekawa:2018oqt , or the cross section of enhanced production of $\Lambda\Lambda$ pair KEK-PSE224:1998trj ; Yoon:2007aq ; KEKE176:2009jzw . The results do not confirm the existence of H-dibaryon, and set the lower limit for the mass of H-dibaryon.

The experiments Belle:2013sba ; BaBar:2018hpv were carried out to search for H-dibaryon or deeply bound singlet $uuddss$ sexaquark $S$ (for the explanation of $S$, see Farrar:2017eqq ) in $\Upsilon\rightarrow S\bar{\Lambda}\bar{\Lambda}$ decay. Their results show no signal for the existence of H-dibaryon or $S$ particle.

As a theoretically tool, Lattice QCD is also used to investigate H-dibaryon. Some quenched studies show that $m_{H}<2\ m_{\Lambda}$, such results support that H-dibaryon is a bound state Iwasaki:1987db ; Luo:2011ar ; Luo:2007zzb , while other quenched studies show that H-dibaryon is not a bound state Mackenzie:1985vv ; Pochinsky:1998zi ; Wetzorke:1999rt ; Wetzorke:2002mx .

Aside from quenched studies, simulations with dynamical fermions have been carried out by NPLQCD, HALQCD and other groups. The NPLQCD collaboration investigated baryon-baryon scattering, and then, extracted the phase shift by employing Lüscher’s method Luscher:1986pf ; Luscher:1990ux to distinguish scattering states from binding states Beane:2003da ; Beane:2006mx ; Beane:2009py ; NPLQCD:2010ocs ; Beane:2011zpa ; NPLQCD:2011naw ; NPLQCD:2012mex . They used this scenario to determine if there exists H-dibaryon, by making simulation with $N_{f}=2+1$ dyamical fermions on anisotropic ensembles Beane:2011zpa ; NPLQCD:2010ocs ; NPLQCD:2011naw , and with $N_{f}=3$ dynamical fermions on isotropic ensembles NPLQCD:2012mex .

The HALQCD collaboration investigated the baryon-baryon interaction in terms of the baryon-baryon potential. They extracted the Nambu-Bethe-Salpeter wave-function by computing the four-point green function on lattice, and then determined the baryon-baryon potential from the Nambu-Bethe-Salpeter wave-function. They used this method to address the existence of H-dibaryon on $N_{f}=3$ ensembles Inoue:2010hs ; Inoue:2010es ; Inoue:2011ai , and on $N_{f}=2+1$ ensembles HALQCD:2019wsz ; Sasaki:2016gpc ; HALQCD:2018lur . Of the results obtained by the two groups, the simulations on $N_{f}=2+1$ ensembles HALQCD:2019wsz ; Sasaki:2016gpc ; HALQCD:2018lur by HALQCD claimed that H-dibaryon may be a $\Lambda\Lambda$ resonance. Other results by the two groups agreed on the presence of the H-dibaryon, despite disagreement on the binding energy Francis:2018qch .

Ref. Francis:2018qch made simulations on ensembles of two dynamical quarks and one quenched strange quark. They applied Lüscher’s method to determine S-wave scattering phase shift with local and bilocal interpolators, and found that for pion mass of 960 MeV, there exists a bound H-dibaryon.

Ref. Green:2021qol carried out simulations by using O(a)-improved Wilson fermions at SU(3) symmetric point with $m_{\pi}=m_{K}\approx 420\,{\rm MeV}$, and their results show that there exists a weakly bound H-dibaryon.

Apart from the search for H-dibaryon, there are lattice QCD calculations for three-flavored heavy dibaryons Junnarkar:2022yak ; Mathur:2022ovu . Three-flavored heavy dibaryons are states with possible quark flavour combinations with at least one of them as charm ($c$) or bottom ($b$) quark.

Besides at zero temperature, the properties of hadrons at finite temperature are also one of the central goals of lattice QCD simulation (see, for example, Aarts:2010ek ; Aarts:2011sm ; Aarts:2013kaa ; Aarts:2014cda ; Kelly:2018hsi ; Aarts:2020vyb ). In the past decades, mesons at finite temperature have been studied extensively. This is not the case for baryons. Baryons at finite temperature are hardly investigated on the lattice. In fact, there are a few lattice studies of baryonic screening and temporal masses DeTar:1987ar ; Pushkina:2004wa ; Aarts:2018glk ; Datta:2012fz ; Aarts:2015mma . Nevertheless, the behaviour of baryons in a hadronic medium is relevant to heavy-ion collisions. Therefore, there is a need to unambiguously understand the property of baryons at finite temperature.

Current research work on H-dibaryon focuses on the problem of its existence at zero temperature from different aspects. In this paper, we make lattice QCD simulations to investigate the masses of the conjectured H-dibaryon and octet baryons at different temperatures. The change of mass of H-dibaryon with temperature is worth studying in its own right theoretically, moreover, the comparison of its mass with $m_{\Lambda}$ can provide some information on its existence.

The paper is organized as follows: In Sec. II, we present the technique details of the simulation which include the definition of correlation functions and the interpolating operators. Sec. III introduces the method of extracting spectral density from correlation function designed in Ref. Hansen:2019idp . Our simulation results are presented in Sec. IV, followed by discussion in Sec. V.

In our simulation, we compute the correlation functions of H-dibaryon and $\Lambda$, we also calculate the correlation function of $N$, $\Sigma$ and $\Xi$. The generic form of correlation function is:

For the H-dibaryon interpolating operator, we choose the local operator. The starting point is the following operator notation for the different six quark combination Francis:2018qch :

where $a,b,\ldots,f$ denote generic quark flavors, and $P_{+}=(1+\gamma_{0})/2$ projects the quark fields to positive parity. We choose the operator $O_{\mathbf{H}}$ as the H-dibaryon interpolating operator which transforms under the singlet irreducible representation of flavor SU(3) Donoghue:1986zd ; Golowich:1992zw ; Wetzorke:1999rt ; Wetzorke:2001tgi :

The H-dibaryon correlation function can be obtained based on the formulae in Ref. Wetzorke:2001tgi .

For the baryon $\Lambda$ interpolating operator, we choose the standard definition which is given by (see, for example, Refs. Leinweber:2004it ; Montvay:1994cy ; Gattringer:2010zz ):

For the baryon $\Lambda$ correlation function, we choose the standard definition (see, for example, Refs. Leinweber:2004it ; Montvay:1994cy ; Gattringer:2010zz ).

For the baryon $N$, $\Sigma$ and $\Xi$, we take the standard definition of interpolating operator, and the corresponding definition of correlator (also, see, Refs. Leinweber:2004it ; Montvay:1994cy ; Gattringer:2010zz ).

After we get the correlation function, the mass can be obtained by fitting the exponential ansatz:

with $m_{+}$ being the mass of the particle of interest, and $\tau$ in the interval $0\leq\tau<1/T$ on a lattice at finite temperature $T$.

In order to get the ground state energy of the particle concerned, it is best to choose large time extent lattice. However, at finite temperature, if we make simulation on large time extent lattice, the lattice spacing must be chosen to be small. Therefore, it is a dilemma for us to make lattice simulation at finite temperature presently. In order to get the ground state energy as possible as we can on a relatively small time extent lattice, it is expedient for us to take an extrapolation method in our procedure to get the ground state mass. We fit equation (5) to correlators in a series of time range $[\tau_{1},\tau_{2}]$ where $\tau_{2}$ is fixed to the whole time extent, and $\tau_{1}$ runs over several values from $\tau_{1}=1,2,3,4,...$, then we can get a series of mass values which correspond to different early Euclidean time slices suppression. After that, we plot the mass values obtained in different time interval $[\tau_{1},\tau_{2}]$ against $1/\tau_{1}$, and fit a linear expression to those mass values, then, extrapolate the linear expression to $\tau_{1}\to\infty$.

Hadron properties are encoded in spectral functions which can provide us important information on hadrons. Two approaches and their variants are adopted to reconstruct spectral function. The first is the maximum entropy method and their variants Asakawa:2000tr ; Burnier:2013nla ; Aarts:2007wj . The second is the Backus-Gilbert method and their variants Brandt:2015sxa ; Brandt:2015aqk ; Lupo:2023qna ; Barone:2023iat ( reviews on the spectral function in lattice QCD can be found in Ref. Rothkopf:2022fyo ; Meyer:2017ydp and references therein ). Recently, based on the Backus-Gilbert method, a new method was presented in Ref. Hansen:2019idp . This method allows for choosing a smearing function at the beginning of the reconstruction procedure. To render this paper self-contained, we briefly present the method which was designed in Ref. Hansen:2019idp in this section. In the following, the notations and symbols are almost the same as those used in Ref. Hansen:2019idp .

The correlation function can be written as:

with $\rho_{L}(E)$ being the spectral function. We choose the basis function as:

We can approximate $\rho_{L}(E_{\star})$ by $\bar{\rho}_{L}(E_{\star})$, where $\bar{\rho}_{L}(E_{\star})$ can be evaluated by

after those coefficients $g_{\tau}(E_{\star})$ are determined.

The coefficients $g_{\tau}(E_{\star})$ are determined by minimizing the linear combination $W[\lambda,g]$ of the deterministic functional $A[g]$ and error functional $B[g]$

under the unit area constraint

where $A[g]$ is defined as:

with $\bar{\Delta}_{\sigma}(E,E_{\star})$ and $\Delta_{\sigma}(E,E_{\star})$ being the smearing function and target smearing function, respectively. These two functions are given by

and

respectively. $B[g]$ is written as:

with ${\rm Cov}$ is the covariance matrix of the correlation function $G(\tau)$. More details are given in Ref. Hansen:2019idp .

Before presenting the simulation results, we describe the computation details. The simulations are carried out on $N_{f}=2+1$ Generation2 (Gen2) FASTSUM ensembles Aarts:2020vyb of which the ensembles at the lowest temperature are provided by the HadSpec collaboration Edwards:2008ja ; HadronSpectrum:2008xlg , so the computation details are the same as those used in Ref. Aarts:2020vyb . We recompile the simulation details in the following three tables 1,  2, and 3 from Ref. Aarts:2020vyb .

The ensembles are generated with a Symanzik-improved gauge action and a tadpole-improved clover fermion action, with stout-smeared links. The details of the action are given in Ref. Aarts:2020vyb . The parameters in the lattice action are recompiled in Table 1. The $N_{f}=2+1$ Gen2 ensembles correspond to a physical strange quark mass and a bare light quark mass of $a_{\tau}m_{l}=-0.0840$, yielding a pion mass of $m_{\pi}=384(4)$ MeV (see Table 2).

The ensemble detail is listed in Table 3 which is recompiled from Ref. Aarts:2020vyb with a slight difference on ensemble $N_{s}^{3}\times N_{\tau}=32^{3}\times 48$. The corresponding physical parameters such as the lattice spacing, and the pion mass etc are collected in Table  2.

The quark propagators are computed by using the deflation-accelerated algorithm Luscher:2007es ; Luscher:2007se . When computing the propagator, The spatial links are stout smeared Morningstar:2003gk with two steps of smearing, using the weight $\rho=0.14$. For the sources and sinks, we use the Gaussian smearing Gusken:1989ad

where $H$ is the spatial hopping part of the Dirac operator and $C$ an appropriate normalisation Aarts:2015mma .

The correlators of $\Lambda$ and H-dibaryon are presented in Fig. 1 and Fig. 2, respectively. For the correlators of $\Lambda$ and H-dibaryon, we find similar behavior which was displayed in Fig. 1 in Ref. Aarts:2015mma for $N$. For the correlator of $\Lambda$ on large $N_{\tau}$ and relatively small $N_{s}$ lattice, especially $24^{3}\times 128$ lattice, some correlator data points are negative, and these points are not displayed on the plot, because the vertical axis is rescaled logarithmically. At some points, the error bar looks strange, it is because at these points, the errors are the magnitude of the correlator value, and the vertical axis is rescaled. For the plot of H-dibaryon correlator, the same observation can be observed.

We use the extrapolation method to extract the ground state masses for $N$, $\Xi$, $\Sigma$, $\Lambda$, and H-dibaryon. We first fit equation (5) to correlator by suppressing different early time slices to get a series of mass values. We present the results of nucleon and H-dibaryon on lattice $N_{\tau}=128$ in Fig. 3. After we get a series of mass values with different early time slices suppressed, we extrapolate the mass values linearly with the scenario described in the last paragraph in Sec. II. We present the results of linear extrapolation for nucleon and H-dibaryon on lattice $N_{\tau}=128$ in Fig. 4. In the extrapolation procedure, we use one portion of the data presented in Fig. 3.

The results are listed in table 3. We can find that the masses decrease when temperature increases. We compare our results of $N$ and $\Lambda$ below $T_{c}$ with those in Refs. Aarts:2015mma ; Aarts:2018glk . The results are consistent within errors.

We also calculate the spectral density $\bar{\rho}_{L}(E_{\star})$ of the correlation function of $N$, $\Sigma$, $\Xi$, $\Lambda$ and H-dibaryon by using the public computer program Hansen:code .

We present the spectral density with different $\sigma=0.02,\ 0.04,\ 0.06,\ 0.08$ for $N$ , $\Xi$ and H-dibaryon at three temperatures in Fig. 5. The upper panel in Fig. 5 for $N$ at $T/T_{c}=0.24$ indicates that too large $\sigma$ value may skip peak structure of spectral density. From the upper panel in Fig. 5, we can find that the spectral density distribution obtained by using $\sigma=0.08$ has just one position where $\bar{\rho}_{L}(E_{\star})$ takes locally maximum value. The position is about at $E_{\star}=0.37$. At $T/T_{c}=0.24$, the time extent $N_{\tau}=128$ is large enough to extract the ground state energy.

However, even if we do not suppress any early Euclidean time slices in the fitting procedure with equation (5), we cannot get a mass value $a_{\tau}m_{N}$ which is larger than $0.30$. The largest value of $a_{\tau}m_{N}$ we get by suppressing different number of early time slices is about $0.25$ which is smaller than 0.30. It can be seen clearly from Fig. 3. The mass value of 0.30 is somewhat an arbitrary value between the two peak positions of 0.17 and 0.38 for the spectral density from table 4. So we think taking large $\sigma$ may lead to missing some peak structure. On the other hand, spectral density $\bar{\rho}_{L}(E_{\star})$ obtained by using $\sigma=0.02$ in the upper panel of Fig. 5 has a peak position at $E_{\star}\approx 0.05$ with small peak value. This peak structure may be due to the lattice artefact.

The middle panel in Fig. 5 for $\Xi$ at $T/T_{c}=0.95$ shows that smaller $\sigma$ value can make peak structure of spectral density in small $E_{\star}$ region more pronounced. The lower panel for H-dibaryon at $T/T_{c}=1.90$ suggests that different $\sigma$ value has little effect on the computation of spectral density at high temperature. So, we just present the spectral density results computed with $\sigma=0.020$ in the following.

The spectral density $\bar{\rho}_{L}(E_{\star})$ of $\Xi$, $\Lambda$ and H-dibaryon are given in Fig. 6,  7 and  8, respectively. The spectral density $\bar{\rho}_{L}(E_{\star})$ of $N$ and $\Sigma$ has similar behaviour to that of $\Lambda$. From Fig. 6 , 7 and  8, we can find that the spectral density $\bar{\rho}_{L}(E_{\star})$ of $\Xi$ and $\Lambda$ has similar behaviour, while $\bar{\rho}_{L}(E_{\star})$ of H-dibaryon is slightly different. All the peak positions of $\bar{\rho}_{L}(E_{\star})$ are collected in Table 4,5,6,8 and 7.

From Fig. 6,  7 and  8, we can find that at the lowest temperature $T/T_{c}=0.24$, the spectral density for $\Xi$, $\Lambda$ and H-dibaryon has rich peak structure. Despite there are two peaks approximately between $E_{\star}=0.20$ and $E_{\star}=0.40$, the spectral density $\bar{\rho}_{L}(E_{\star})$ of $\Xi$ and $\Lambda$ in the range of $E_{\star}$ from $0.20$ to $0.40$ are almost the same. The mass values $a_{\tau}m_{\Xi}=0.2459$, $a_{\tau}m_{\Lambda}=0.2299$ obtained by extrapolation method are in that range of $E_{\star}$.

However, $a_{\tau}m_{H}=0.44$ at $T/T_{c}=0.24$ for H-dibaryon is in the neighbour of peak position $E_{\star}=0.35$ where the $\bar{\rho}_{L}(E_{\star})$ value is not very large. Obtaining $a_{\tau}m_{H}=0.44$ at $T/T_{c}=0.24$ is just by suppressing more early Euclidean time slices. From the upper panel of Fig. 8, more high frequency components of the spectral density should be suppressed in the extrapolation procedure.

When temperature increases, the multi-peak structure of spectral density distribution turns into two-peak structure for $\Xi$ and $\Lambda$ until at high temperature $T/T_{c}=1.27,1.52,1.90$, the spectral density distribution has one peak.

At the intermediate temperatures, the spectral density $\bar{\rho}_{L}(E_{\star})$ has a two-peak structure. If we take the smaller values of $E_{\star}$ at peak positions as the ground state energies of corresponding particle, then these mass values obtained by the peak position of $\bar{\rho}_{L}(E_{\star})$ are smaller than those mass values obtained in Ref. Aarts:2015mma and Ref. Aarts:2018glk . Mass values of $a_{\tau}m_{\Xi}$, $a_{\tau}m_{\Lambda}$ and $a_{\tau}m_{H}$ presented in Table. 3 are not consistent with the peak positions of corresponding spectral density. This observation shows that the mass values obtained by extrapolation method are affected by the two-peak structure of spectral density.

At high temperature $T/T_{c}=1.27,1.52,1.90$, for the spectral density distribution for $N$, the spectral density exhibits one peak structure, and the peak position shifts towards large value with increasing temperature. In the meantime, we can find the peak broadens and becomes smooth. It means that in the mass spectrum structure of nucleon, there is no $\delta$ function structure contributing to the correlation function. We can find this observation from Fig. 9 for $N$. Similar behaviour can be found for $\Sigma$, $\Xi$ and $\Lambda$. We think the smooth distribution of spectral density implies that there does not exist one-particle state at high temperature.

This is not the case for H-dibaryon. The spectral density distribution for H-dibaryon at $T/T_{c}=1.27,1.52,1.90$ is presented in Fig. 10 from which we can find that at $T/T_{c}=1.27,1.52$, the spectral density distribution still exhibits one peak structure until at $T/T_{c}=1.90$, the spectral density distribution broadens and becomes smooth. This observation may imply that at temperature $T/T_{c}=1.27,1.52$, H-dibaryon still remains as a one-particle state.

We have made a simulation in an attempt to determine the masses of the conjectured H-dibaryon with $2+1$ flavor QCD with clover fermion at nine different temperatures. In the meantime, we also calculate the masses of $N$, $\Sigma$, $\Xi$ and $\Lambda$. The results are collected in table 3. The spectral density distribution of those particle’s correlation function are computed to understand the mass spectrum obtained by extrapolation method.

In our simulation, the change of temperature is represented by the change of $T/T_{c}$. $T_{c}$ is pseudocritical temperature determined via renormalized Polyakov loop and estimated to be $T_{c}=185(4)\,{\rm MeV}$ Aarts:2014nba ; Aarts:2020vyb .

We have compared two scenarios to obtain the ground state mass as possible as we can. One is to suppress more early Euclidean time slices in the fitting procedure with equation (5). The second method is the extrapolation method. We extrapolate some fitting results which are obtained with different Euclidean time slices suppression to time approaching infinity. The results by the two methods are consistent with each other within errors. We just present the results by the extrapolation method in table 3. In fact, among the series of mass values obtained by different early time suppression, it is difficult to choose which mass value is the proper one. However, the extrapolation method can alleviate this difficulty to some extent.

The analysis of spectral density can provide insights on mass spectrum. However, in our simulation, we can find that the mass values obtained by extrapolation method are not consistent with the peak position of spectral density in some situations, especially at high temperature. Under such situations, we think the results of extrapolation method are more reliable. We can take the peak positions for nucleon in Table 4 as an example to give an explanation. The smaller values of peak position are increasing with temperature. However, with increasing temperature, the mass values of particle concerned are supposed to decrease.

At the lowest temperature $T/T_{c}=0.24$, the spectral density distribution has rich peak structure. The mass spectrum of particles approximately reflects the peak position of spectral density distribution.

At the intermediate temperatures, the spectral density distribution exhibits a two-peak structure. The peak structure at the smaller $E_{\star}$ becomes smooth gradually when $\sigma$ increases. Considering the quark mass which corresponds to $m_{\pi}=384(4){\rm MeV}$, if we take the smaller $E_{\star}$ at the peak position as the ground state energy, the mass value is too small. So we think the mass values obtained by using extrapolation method are affected by the two states.

At high temperature, despite we get the mass values for $N$, $\Sigma$, $\Xi$ and $\Lambda$ which are presented in table 3 by extrapolation method, the spectral density distribution appears to become smooth which implies there does not exist one state.

H-dibaryon is a multi-baryon state. From the spectral density distribution, it is found that at the lowest temperature $T/T_{c}=0.24$, the multi-state structure manifests. When temperature increases, the number of peaks decreases until at $T/T_{c}=1.90$, the spectral density distribution becomes almost smooth. It means it is likely that H-dibaryon survives beyond $T_{c}$ until it melts down at $T/T_{c}=1.90$. Considering that H-dibaryon is a multi-baryon state, this conclusion awaits further investigation.

It is appropriate to consider the lowest temperature ensembles $N_{\tau}=128$ to be the zero temperature ones since $N_{\tau}>\xi N_{s}$ Aarts:2020vyb . Using the mass values of H-dibaryon and $\Lambda$ in Table 3 at $N_{\tau}=128,T/T_{c}=0.24$, an estimation of $\Delta m=m_{H}-2m_{\Lambda}$ can be made. $\Delta m=m_{H}-2m_{\Lambda}=-0.0026(11)$ which is converted into physical unit to be $\Delta m=m_{H}-2m_{\Lambda}=-14.6(6.2){\rm MeV}$. Ref. Beane:2006mx ; NPLQCD:2012mex ; Inoue:2010es ; Inoue:2011ai ; Francis:2018qch ; Green:2021qol showed the presence of a binding state of H-dibaryon, despite the disagreement on the binding energy values. Ref. Beane:2006mx ; NPLQCD:2012mex reported a binding energy of $16.6\pm 2.1\ {\rm MeV}$ at $m_{\pi}=389\ {\rm MeV}$, and $74.6\pm 4.7\ {\rm MeV}$ at $m_{\pi}=800\ {\rm MeV}$. Ref. Inoue:2010es obtained a binding energy of H-dibaryon $30-40\ {\rm MeV}$ for pion mass $673-1015\ {\rm MeV}$. Ref. Inoue:2011ai had the similar results. Ref. Francis:2018qch published the binding energy $19\pm 10\ {\rm MeV}$ of H-dibaryon for $m_{\pi}=960\ {\rm MeV}$. Ref. Green:2021qol presented an estimation of the binding energy $4.56\pm 1.13\ {\rm MeV}$ in the continuum limit at the SU(3)-symmetric point with $m_{\pi}=m_{K}\approx 420\ {\rm MeV}$ ( for the binding energy versus pion mass, see also, Fig. 5 in Ref. Green:2021qol ).

In our simulation, the correlators of proton, $\Lambda$ and H-dibaryon at $N_{\tau}=128$ have negative values, and in the fitting process for H-dibaryon, we drop the negative values. We guess the emergence of negative values of the correlators is due to deterioration of the signal-to-noise ratio.

Our simualtions are at $m_{\pi}=384(4)\ {\rm MeV}$ which is far from physical pion mass, so simulations with lower pion mass are expected to give us more information about the properties of H-dibaryon.