Masses of the conjecturedH-dibaryon at different temperatures

Quantum chromodynamics (QCD) describes the dynamics of quarks and gluons, underpinning nuclear physics from the 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 a low energy scale, the perturbative method cannot be applied to explain low energy phenomena of nuclear physics. Fortunately, lattice QCD, which is based on first principles, can be employed to make precise predictions of hadronic quantities, especially phenomena that are difficult to explore in the laboratory, for example, the mass spectrum of baryons at different high temperatures.

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

Unlike mesons and baryons, this exotic hadron may be relevant for the hypernuclei and the strange matter that could exist in the core of neutron stars. Moreover, it is a potential candidate for dark matter [2]. Thus, this prediction triggered an intense search for such a state, both experimentally [3–11] and theoretically [12–31].

The observation of double hypernuclei, $ ^A_{\Lambda\Lambda}Z $ , is critical in connection with the existence of theH-dibaryon [3]. If the mass of theH-dibaryon, $ m_H $ , is much smaller than the mass of the double Λ hyperon, $ 2\,m_{\Lambda} $ , a double hypernuclei may decay into anH-dibaryon and a residual nucleus by strong interaction. In this case, the branching ratio for the decay of double hypernuclei through weak interaction is very small, which means that it cannot be observed in practice [3]. Therefore, the observation of the weak decay of double hypernuclei will set a limitation on the mass of theH-dibaryon, $m_H > 2\,m_{\Lambda} - B_{\Lambda\Lambda}$ , with $ B_{\Lambda\Lambda} $ being the binding energy of $ \Lambda\Lambda $ hyperons.

Experiments [3–8,11] investigated the nuclear capture of $ \Xi^- $ at rest produced in the $ (K^-,K^+) $ reaction, analyzed the sequential weak decay of double hypernuclei, and measured the binding and interaction energies of $ \Lambda\Lambda $ [3,5,8,11], or the cross section of the enhanced production of the $ \Lambda\Lambda $ pair [4,6,7]. The results do not confirm the existence of theH-dibaryon, but set a lower limit for its mass.

Experiments [9,10] were also conducted to search for theH-dibaryon or deeply bound singlet $ uuddss $ sexaquarkS(an explanation ofSis provided in [2]) in the $ \Upsilon \rightarrow S\bar\Lambda\bar\Lambda $ decay. The results show no evidence for the existence of theH-dibaryon orSparticle.

Lattice QCD is used as a theoretical tool to investigate theH-dibaryon. Some quenched studies show that $ m_H < 2\ m_{\Lambda} $ , suggesting that theH-dibaryon is a bound state [12–14]. However, other quenched studies show that theH-dibaryon is not a bound state [15–18].

In addition to quenched studies, simulations with dynamical fermions have been carried out by NPLQCD, HALQCD, and other groups. The NPLQCD collaboration investigated baryon-baryon scattering and extracted the phase shift by employing Lüscher's method [32,33] to distinguish scattering states from binding states [19–23,34,35]. They applied this methodology to determine whether theH-dibaryon exists by simulating with $ N_f=2+1 $ dyamical fermions on anisotropic ensembles [20–22], and with $ N_f=3 $ dynamical fermions on isotropic ensembles [23].

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 applied this method to investigate the existence of theH-dibaryon on $ N_f=3 $ [24–26] and $ N_f=2+1 $ [27–29] ensembles. The simulations on $ N_f=2+1 $ ensembles [27–29] conducted by HALQCD suggest that theH-dibaryon may be a $ \Lambda\Lambda $ resonance. Other results obtained by the two groups agreed on the presence of theH-dibaryon, despite disagreement on the binding energy [30].

Ref. [30] reports on the simulation of ensembles of two dynamical quarks and one quenched strange quark. The authors applied Lüscher's method to determine theS-wave scattering phase shift with local and bilocal interpolators, and found a boundH-dibaryon for a pion mass of 960 MeV.

Ref. [31] reports on simulations based onO(a)-improved Wilson fermions at aSU(3) symmetric point with $m_\pi=m_K\approx 420\;{ {\rm{MeV}}}$ . The results show that there exists a weakly boundH-dibaryon.

In addition to the search for theH-dibaryon, lattice QCD calculations for three-flavored heavy dibaryons have been reported [36,37]. These dibaryons are states with possible quark flavour combinations with at least one of them being charm (c) or bottom (b) quark.

Investigating the properties of hadrons at finite temperature, besides those at zero temperature, is among the central goals of lattice QCD simulations (see, for example, [38–44]). In the past decades, mesons at finite temperature have been studied extensively. This is not the case for baryons, which have been barely investigated at finite temperature on the lattice. There are a few lattice studies on baryonic screening and temporal masses [45–49]. Nevertheless, the behaviour of baryons in a hadronic medium is relevant to heavy-ion collisions. Therefore, there is a need for unambiguously understanding the property of baryons at finite temperature.

The present study on theH-dibaryon focused on the problem of its existence at zero temperature from different viewpoints. We conducted lattice QCD simulations to investigate the masses of the conjecturedH-dibaryon and octet baryons at different temperatures. The change in the mass of theH-dibaryon with temperature is worth studying 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 details of the simulation technique, including the definition of correlation functions and interpolating operators. Sec. III introduces a method for extracting the spectral density from the correlation function designed in Ref. [50]. Our simulation results are presented in Sec. IV, followed by a discussion in Sec. V.

In our simulations, we computed the correlation functions of theH-dibaryon and Λ. We also calculated the correlation function ofN, Σ, and Ξ. The generic form of a correlation function is

and for theH-dibaryon interpolating operator, we chose the local operator. The starting point is the following operator notation for the different combinations of six quarks [30]:

where $ a, b,\ldots,f $ denote generic quark flavors, and $P_+= (1+\gamma_0)/2$ projects the quark fields to positive parity. We chose the operator $O_{H}$ as theH-dibaryon interpolating operator which transforms under the singlet irreducible representation of flavorSU(3) [17,51–53],

and theH-dibaryon correlation function can be obtained from the formulae in Ref. [53].

For the baryon Λ interpolating operator, we chose the standard definition (see, for example, Refs. [54–56]),

and for the correlation function of the baryon Λ, we also chose the standard definition (see, for example, Refs. [54–56]).

For the baryonsN, Σ, and Ξ, we selected the corresponding standard definitions of the interpolating operator and correlator (see Refs. [54–56]).

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

where $ m_{+} $ is the mass of the particle of interest, andτtakes values in the interval $ 0 \leq \tau<1/T $ on a lattice at finite temperatureT.

To obtain the ground state energy of the particle concerned, the best approach is to choose a large time extent lattice. However, at finite temperature, if the simulation is conducted on such a lattice, the lattice spacing must be small. Therefore, performing lattice simulations at finite temperature constitutes a dilemma presently. To obtain the ground state energy as accurately as possible on a relatively small time extent lattice, it is suitable for our procedure to incorporate an extrapolation method to obtain the ground state mass. We fit Eq. (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 $ is swept across several values: $\tau_1=1, 2, 3, ...$ Thus, we obtained a series of mass values corresponding to the suppression of different early Euclidean time slices. Subsequently, we plotted the mass values obtained in different time intervals $ [\tau_1,\tau_2] $ against $ 1/\tau_1 $ and fit a linear expression to those mass values. Finally, we extrapolated the linear expression to $ \tau_1 \to \infty $ .

Hadron properties are encoded in spectral functions, which can provide important information on hadrons. Two approaches and their variants are usually adopted to reconstruct the spectral function. The first is the maximum entropy method and its variants [57–59]. The second is the Backus-Gilbert method and its variants [60–63] (reviews on the spectral function in lattice QCD can be found in Refs. [64,65] and references therein). Recently, a new method based on the Backus-Gilbert method was presented in Ref. [50]. 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 proposed in Ref. [50] in this section. In the following, the notations and symbols are almost the same as those used in Ref. [50].

The correlation function can be expressed as

where $ \rho_L(E) $ is the spectral function. We choose the basis function as

and $ \rho_L(E_\star) $ can be approximated by $ \bar{\rho}_L(E_\star) $ , where $ \bar{\rho}_L(E_\star) $ is evaluated by

once the 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

where $ \bar\Delta_\sigma (E,E_\star) $ and $ \Delta_\sigma (E,E_\star) $ are the smearing and target smearing functions, respectively. These two functions are expressed as

and

respectively. $ B[g] $ is expressed as

where $ { {\rm{Cov}}} $ is the covariance matrix of the correlation function $ G(\tau) $ . Further details are given in Ref. [50].

Before presenting the simulation results, we describe the computation details. The simulations were performed on $ N_f=2+1 $ Generation2 (Gen2) FASTSUM ensembles [43]; the ensembles at the lowest temperatures were those provided by the HadSpec collaboration [66,67]. The computation setup was the same as that used in Ref. [43]. We summarize the simulation details inTables 1,2, and3from Ref. [43].

Gauge coupling (fixed-scale approach)	$ \beta = 1.5 $
tree-level coefficients	$ c_0=5/3,\,c_1=-1/12 $
bare gauge, fermion anisotropy	$ \gamma_g = 4.3 $ , $ \gamma_f = 3.399 $
ratio of bare anisotropies	$ \nu = \gamma_g / \gamma_f = 1.265 $
spatial tadpole (without, with smeared links)	$ u_s = 0.733566 $ , $ \tilde{u}_s = 0.92674 $
temporal tadpole (without, with smeared links)	$ u_\tau = 1 $ , $ \tilde u_\tau = 1 $
spatial, temporal clover coefficient	$ c_s = 1.5893 $ , $ c_\tau = 0.90278 $
stout smearing for spatial links	$ \rho = 0.14 $ , isotropic, 2 steps
bare light quark mass for Gen2	$ \hat m_{0, {\rm{light}}} = -0.0840 $
bare strange quark mass	$ \hat m_{0, {\rm{strange}}} = -0.0743 $
light quark hopping parameter for Gen2	$ \kappa_{{\rm{light}}} = 0.2780 $
strange quark hopping parameter	$ \kappa_{{\rm{strange}}} = 0.2765 $

Table 1.Parameters in the lattice action. This table is recompiled from Ref. [43].

$ a_\tau $ /fm	0.0350(2)
$ a_\tau^{-1} $ /GeV	5.63(4)
$ \xi=a_s/a_\tau $	3.444(6)
$ a_s $ /fm	0.1205(8)
$ N_s $	24
$ m_\pi $ /MeV	384(4)
$ m_\pi L $	5.63

Table 2.Parameters such as lattice spacing and pion mass from Ref. [43] for Generation 2 ensemble.

$ N_s $	$ N_\tau $	$T /{\rm MeV}$	$ T/T_c $	$ N_{{\rm{cfg}}} $	$ a_\tau m_N $	$ a_\tau m_\Sigma $	$ a_\tau m_\Xi $	$ a_\tau m_{\Lambda} $	$ a_\tau m_H $
24	128	44	0.24	304	0.2133(24)(6)	0.2349(21)(5)	0.2459(18)(5)	0.2299(23)(6)	0.457(27)(5)
32	48	117	0.63	601	0.208(2)(3)	0.231(1)(2)	0.243(1)(2)	0.226(2)(3)	0.448(17)(4)
24	40	141	0.76	502	0.203(2)(5)	0.228(2)(4)	0.239(2)(4)	0.221(2)(4)	0.437(14)(9)
24	36	156	0.84	501	0.196(2)(6)	0.221(2)(5)	0.231(2)(5)	0.214(2)(5)	0.42(1)(2)
24	32	176	0.95	1000	0.181(2)(9)	0.204(2)(8)	0.215(2)(7)	0.199(2)(8)	0.393(7)(21)
24	28	201	1.09	1001	0.179(2)(12)	0.191(2)(12)	0.201(2)(11)	0.190(2)(11)	0.38(1)(2)
24	24	235	1.27	1001	0.172(3)(15)	0.179(3)(15)	0.191(3)(14)	0.182(3)(14)	0.36(1)(3)
24	20	281	1.52	1000	0.159(4)(18)	0.164(4)(18)	0.176(4)(17)	0.169(4)(17)	0.33(1)(4)
24	16	352	1.90	1000	0.154(6)(24)	0.158(6)(24)	0.171(6)(23)	0.164(6)(23)	0.31(2)(4)

Table 3.Spatial and temporal extent, temperature in MeV, number of configurations, and masses ofN, Σ, Ξ, Λ, andH-dibaryon. Masses of baryons andH-dibaryon are obtained by the extrapolation method. Estimates of statistical and systematic errors are contained in the first and second brackets, respectively. The errors of the fitting parameters of linear extrapolation are the systematic errors of the masses. The ensembles at the lowest temperatures are those provided by HadSpec [66,67] (Gen2).

The ensembles were 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. [43]. The parameters in the lattice action are recompiled inTable 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 (seeTable 2).

The ensemble details are listed inTable 3, which is recompiled from Ref. [43] with a slight difference on the ensemble $ N_s^3\times N_\tau = 32^3\times 48 $ . The corresponding physical parameters, such as lattice spacing and pion mass, are listed inTable 2.

The quark propagators were computed using the deflation-accelerated algorithm [68,69]. When computing the propagator, the spatial links were stout smeared [70] with two steps of smearing, using the weight $ \rho = 0.14 $ . For the sources and sinks, we used Gaussian smearing [71],

whereHis the spatial hopping part of the Dirac operator, andCis an appropriate normalization factor [49].

The correlators of the Λ andH-dibaryon are presented inFigs. 1and2, respectively. For the correlators of the Λ andH-dibaryon, we found a similar behavior to that displayed inFig. 1in Ref. [49] forN. For the correlator of Λ on large $ N_\tau $ and relatively small $ N_s $ lattice, especially the $ 24^3\times 128 $ lattice, some correlator data points take negative values. These points are not depicted in 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 reach the magnitude of the correlator value, and the vertical axis is rescaled. For the plot of theH-dibaryon correlator, the same observation can be made.

Figure 1.(color online) Euclidean correlator $ G(\tau)/G(0) $ of Λ as a function of $ \tau T $ at different temperatures. At the lowest temperature $ T/T_c =0.24 $ , the correlators at some points are not displayed because they take negative values.

Figure 2.(color online) Euclidean correlator $ G(\tau)/G(0) $ ofH-dibaryon as a function of $ \tau T $ at different temperatures. At the lowest temperature $ T/T_c =0.24 $ , the correlators at some points are not displayed because they take negative values.

We used the extrapolation method to extract the ground state masses forN, Ξ, Σ, Λ, andH-dibaryon. We first fit Eq. (5) to the correlator by suppressing different early time slices to obtain a series of mass values. We present the results of nucleon andH-dibaryon on lattice $ N_\tau=128 $ inFig. 3. After obtaining a series of mass values with different early time slices suppressed, we extrapolated the mass values linearly according to the scenario described in the last paragraph in Sec. II. We present the results of linear extrapolation for the nucleon andH-dibaryon on lattice $ N_\tau=128 $ inFig. 4. In the extrapolation procedure, we used one portion of the data presented inFig. 3.

Figure 3.(color online) Mass values of nucleon andH-dibaryon obtained by fitting Eq. (5) to correlators on $ N_\tau=128 $ ensembles. The mass values are the fitting parameter $ m_+ $ in Eq. (5) extracted by the fitting procedure in different intervals $ [\tau_1, \tau_2] $ . The horizontal axis label $ \tau_1 $ represents different number of time slices suppressed corresponding to the lower bound of the interval $ [\tau_1, \tau_2] $ .

Figure 4.(color online) Linear extrapolation of mass values for nucleon andH-dibaryon on $ N_\tau=128 $ ensembles. The horizontal axis represents inverse values of time slices suppressed.

The results are listed inTable 3. Note that the masses decrease when the temperature increases. We compare our results ofNand Λ below $ T_c $ with those in Refs. [45,49]. The results are consistent within errors.

We also calculated the spectral density $ \bar\rho_L(E_\star) $ of the correlation function ofN, Σ, Ξ, Λ, andH-dibaryon using a public computer program [72].

We present the spectral density for $ \sigma = $ 0.02, 0.04, 0.06, and 0.08 andN, Ξ, andH-dibaryon at three temperatures inFig. 5. The upper panel inFig. 5forNat $ T/T_c = 0.24 $ indicates that too largeσvalues may skip the peak structure of spectral density. Note from the upper panel inFig. 5that the spectral density distribution obtained for $ \sigma = 0.08 $ has only one position where $ \bar\rho_L(E_\star) $ takes a local maximum value. This position is approximately 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.

Figure 5.(color online) Spectral density computed with different values of the parameterσfor the target smearing function $ \Delta_\sigma (E,E_\star) $ forH-dibaryon, Ξ, andNat different temperatures.

However, even if we do not suppress any early Euclidean time slices in the fitting procedure based on Eq. (5), we cannot obtain a mass value $ a_\tau m_N $ larger than $ 0.30 $ . The largest value of $ a_\tau m_N $ that we obtained by suppressing different numbers of early time slices was approximately 0.25, which is smaller than 0.30, as can be clearly seen fromFig. 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 presented inTable 4. Thus, we conclude that setting a largeσvalue may lead to missing some peak structures. By contrast, the spectral density $ \bar\rho_L(E_\star) $ obtained for $ \sigma = 0.02 $ in the upper panel ofFig. 5has a peak position at $ E_\star \approx 0.05 $ with a small peak value. This peak structure may be due to lattice artifacts.

	$ N_\tau $	$ E_\star $	$ E_\star $	$ E_\star $
N	128	0.05	0.17	0.38
	48	0.06	0.35	–
	40	0.09	0.42	–
	36	0.10	0.52	–
	32	0.11	0.56	–
	28	0.18	0.57	–
	24	0.26	–	–
	20	0.24	–	–
	16	0.43	–	–

Table 4.For different $ N_\tau $ lattices, peak position $ E_\star $ of the spectral density forN.

	$ N_\tau $	$ E_\star $	$ E_\star $	$ E_\star $
Σ	128	0.05	0.17	0.37
	48	0.06	0.35	–
	40	0.08	0.41	–
	36	0.09	0.49	–
	32	0.11	0.56	–
	28	0.17	0.57	–
	24	0.20	–	–
	20	0.22	–	–
	16	0.39	–	–

Table 5.For different $ N_\tau $ lattices, peak position $ E_\star $ of the spectral density for Σ.

	$ N_\tau $	$ E_\star $	$ E_\star $
Ξ	128	0.05	0.33
	48	0.06	0.34
	40	0.08	0.41
	36	0.09	0.48
	32	0.10	0.55
	28	0.13	0.63
	24	0.17	0.71
	20	0.20	–
	16	0.32	–

Table 6.For different $ N_\tau $ lattices, peak position $ E_\star $ of the spectral density for Ξ.

	$ N_\tau $	$ E_\star $	$ E_\star $	$ E_\star $
Λ	128	0.05	0.18	0.38
	48	0.06	0.35	–
	40	0.08	0.42	–
	36	0.09	0.49	–
	32	0.11	0.56	–
	28	0.14	0.61	–
	24	0.18	0.64	–
	20	0.22	–	–
	16	0.35	–	–

Table 7.For different $ N_\tau $ lattices, peak position $ E_\star $ of the spectral density for Λ.

	$ N_\tau $	$ E_\star $	$ E_\star $	$ E_\star $	$ E_\star $
H-dibaryon	128	0.10	0.19	0.35	0.69
	48	0.05	0.24	0.64	–
	40	0.06	0.32	0.78	–
	36	0.07	0.35	–	–
	32	0.08	0.45	–	–
	28	0.10	0.52	–	–
	24	0.11	–	–	–
	20	0.14	–	–	–
	16	0.17	–	–	–

Table 8.For different $ N_\tau $ lattices, peak position $ E_\star $ of the spectral density forH-dibaryon.

The middle panel inFig. 5for Ξ at $ T/T_c = 0.95 $ shows that smallerσvalues can lead to a more pronounced peak structure of the spectral density in small $ E_\star $ regions. The lower panel for theH-dibaryon at $ T/T_c = 1.90 $ suggests that setting differentσvalues has little effect on the computation of spectral density at high temperature. Therefore, we only present the spectral density results computed for $ \sigma =0.020 $ in the following.

The spectral density $ \bar\rho_L(E_\star) $ of Ξ, Λ, andH-dibaryon is presented inFigs. 6,7, and8, respectively. The spectral density $ \bar\rho_L(E_\star) $ ofNand Σ has a similar behaviour to that of Λ. Note fromFigs. 6,7,and8that the spectral density $ \bar\rho_L(E_\star) $ of Ξ and Λ has a similar behaviour, while $ \bar\rho_L(E_\star) $ of theH-dibaryon is slightly different. All the peak positions of $ \bar\rho_L(E_\star) $ are provided inTables 4−8.

Figure 6.(color online) Spectral density of Ξ at different temperatures.

Figure 7.(color online) Spectral density of Λ at different temperatures.

Figure 8.(color online) Spectral density ofH-dibaryon at different temperatures.

Note fromFigs. 6,7,and8that, at the lowest temperature $ T/T_c = 0.24 $ , the spectral density for Ξ, Λ, andH-dibaryon has a rich peak structure. Despite the two peaks approximately located between $ E_\star=0.20 $ and $ E_\star=0.40 $ , the spectral density $ \bar\rho_L(E_\star) $ of Ξ and Λ in the range of $ E_\star $ from $ 0.20 $ to $ 0.40 $ are approximately the same. The mass values $ a_\tau m_\Xi = 0.2459 $ and $ a_\tau m_\Lambda = 0.2299 $ obtained by the extrapolation method lie in that range of $ E_\star $ .

However, $ a_\tau m_H = 0.44 $ at $ T/T_c = 0.24 $ for theH-dibaryon is in the neighborhood of the peak position $ E_\star = 0.35 $ , where the $ \bar\rho_L(E_\star) $ value is not very large. Note that $ a_\tau m_H = 0.44 $ at $ T/T_c = 0.24 $ is obtained by suppressing more early Euclidean time slices. The upper panel ofFig. 8shows that more high frequency components of the spectral density should be suppressed in the extrapolation procedure.

When temperature increases, the multi-peak structure of the spectral density distribution turns into a two-peak structure for Ξ and Λ, and at high temperatures, i.e., $ T/T_c = $ 1.27, 1.52, and 1.90, the spectral density distribution presents one peak.

At intermediate temperatures, the spectral density $ \bar\rho_L(E_\star) $ exhibits a two-peak structure. If we take the smaller values of $ E_\star $ at peak positions as the ground state energies of the 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 Refs. [49] and [45]. The mass values of $ a_\tau m_\Xi $ , $ a_\tau m_\Lambda $ , and $ a_\tau m_H $ presented inTable 3are not consistent with the peak positions of the corresponding spectral density. Note that the mass values obtained by the extrapolation method are affected by the two-peak structure of the spectral density.

At high temperatures, i.e., $ T/T_c = $ 1.27, 1.52, and 1.90, the spectral density distribution forNexhibits one peak structure, and the peak position shifts towards large values with increasing temperature. Note that the peak broadens and becomes smooth. It means that in the mass spectrum structure of the nucleon, there is noδfunction structure contributing to the correlation function, as shown inFig. 9forN. A similar behaviour can be found for Σ, Ξ, and Λ. The smooth distribution of the spectral density implies that a one-particle state does not exist at high temperature.

Figure 9.(color online) Spectral density distribution ofNat different temperatures: $ T/T_c = $ 0.84, 1.27, 1.52, and 1.90. At $T/T_c = 0.84,~ N_\tau = 36$ , the spectral density ofNexhibits two peaks. At $ T/T_c = $ 1.27, 1.52, and 1.90, the spectral density distribution becomes approximately smooth.

This is not the case for theH-dibaryon. Its spectral density distribution at $ T/T_c = 0 $ , 1.27, 1.52, and 1.9 is presented inFig. 10. Note that, at $ T/T_c = $ 1.27 and 1.52, the spectral density distribution still exhibits a one peak structure; at $ T/T_c = 1.90 $ , the spectral density distribution broadens and becomes smooth. This observation may imply that, at temperatures $ T/T_c = $ 1.27 and 1.52, theH-dibaryon still exhibits a one-particle state.

Figure 10.(color online) Spectral density distribution of theH-dibaryon at different temperatures: $T/T_c = 0.84, 1.27, 1.52, $ $ {\rm and}\;1.90$ . At $ T/T_c = 0.84, N_\tau = 36 $ , the spectral density of theH-dibaryon exhibits two peaks. At $ T/T_c = 1.27, 1.52 $ , the spectral density distribution has one peak. At $ T/T_c = 1.90 $ , the spectral density distribution becomes approximately smooth.

We conducted simulations aiming to determine the masses of the conjecturedH-dibaryon using the $ 2+1 $ flavor QCD with clover fermion at nine different temperatures. We also calculated the masses ofN, Σ, Ξ, and Λ. The results are presented inTable 3. The spectral density distribution of the correlation function of these particles was computed to understand the mass spectrum obtained by the extrapolation method.

In our simulations, the change in temperature was represented by the change in $ T/T_c $ , where $ T_c $ is the pseudocritical temperature determined via the renormalized Polyakov loop and estimated to be $T_c = 185(4)~ {\rm MeV}$ [43,73].

We compared two scenarios to obtain the ground state mass as accurately as possible. One is based on suppressing more early Euclidean time slices in the fitting procedure with Eq. (5). The second is based on the extrapolation method. We extrapolated some fitting results obtained for different suppressions of Euclidean time slices and time approaching infinity. The results from both methods are consistent within errors. We only present the results from the extrapolation method inTable 3. Among the series of mass values obtained by different early time suppressions, 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 simulations, we found that the mass values obtained by the extrapolation method are not consistent with the peak position of the spectral density in some cases, especially at high temperature. In this case, we consider that the results provided by the extrapolation method are more reliable. The peak positions for the nucleon inTable 4can be selected as an example to give an explanation. The smaller values of the peak positions increase with temperature. However, the mass values of the corresponding particle are supposed to decrease with increasing temperature.

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

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

At high temperature, we obtained the mass values forN, Σ, Ξ, and Λ, listed inTable 3, using the extrapolation method. Note that the spectral density distribution appears to become smooth, which implies that a one-particle state does not exist.

TheH-dibaryon is a multi-baryon state. Note from the spectral density distribution that, at the lowest temperature $ T/T_c = 0.24 $ , the multi-state structure manifests. When temperature increases, the number of peaks decreases until reaching $ T/T_c = 1.90 $ , at which the spectral density distribution becomes approximately smooth. This means that it is likely that theH-dibaryon survives beyond $ T_c $ until it melts down at $ T/T_c = 1.90 $ . Considering that theH-dibaryon is a multi-baryon state, this conclusion calls for further investigation.

It is appropriate to consider the lowest temperature ensembles $ N_\tau=128 $ to be those of the zero temperature given that $ N_\tau>\xi N_s $ [43]. Using the mass values of theH-dibaryon and Λ inTable 3at $ 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 units to be $ \Delta m = m_H - 2m_\Lambda = -14.6(6.2) $ MeV. Refs. [23,25,26,30,31,35] reported the presence of a binding state of theH-dibaryon, despite disagreement on the binding energy values. Refs. [23,35] reported a binding energy of $ 16.6 \pm 2.1 $ MeV at $ m_\pi= 389 $ MeV, and $ 74.6 \pm 4.7 $ MeV at $ m_\pi= 800 $ MeV. Ref. [25] reported a binding energy for theH-dibaryon of $ 30-40 $ MeV for a pion mass of $ 673-1015 $ MeV. Ref. [26] reported similar results. Ref. [30] published a binding energy of $ 19 \pm 10 $ MeV for theH-dibaryon at $ m_\pi = 960 $ MeV. Ref. [31] presented an estimation of the binding energy, $ 4.56 \pm 1.13 $ MeV, in the continuum limit at theSU(3)-symmetric point with $ m_\pi = m_K \approx 420 $ MeV (for the binding energy versus pion mass, seeFig. 5in Ref. [31]).

In our simulations, the correlators of proton, Λ, andH-dibaryon at $ N_\tau = 128 $ took negative values; in the fitting process for theH-dibaryon, we dismissed these negative values. We hypothesize that the emergence of negative values in the correlators is due to deterioration of the signal-to-noise ratio.

Our simulations were conducted at $ m_\pi = 384(4) \ {{\rm{MeV}}} $ , which is far from the physical pion mass. Therefore, simulations with a lower pion mass are expected to provide more information about the properties of theH-dibaryon.

We thank Gert Aarts, Simon Hands, Chris Allton, and Jonas Glesaaen for their valuable help, and Chris Allton for fruitful discussions about the extrapolation method. We modified the adapted version of OpenQCD code[74]to carry out the simulations and used a computer program[72]to calculate the spectral density of the correlation function. The adaptation of OpenQCD code is publicly available[75].The simulations were carried out on $ N_f=2+1 $ Generation2 (Gen2) FASTSUM ensembles[43];the ensembles at the lowest temperature were provided by the HadSpec collaboration[66,67].