Masses of the conjectured H-dibaryon at different temperatures

We present a lattice QCD determination of masses of the conjectured H-dibaryon $m_H$ at nine different temperatures $T/T_c =0.24, 0.63, 0.76, 0.84, 0.95, 1.09, 1.27, 1.52, 1.90$. In the meantime, the masses of baryon $N$, $\Sigma$, $\Xi$ and $\Lambda$ at different temperatures are also computed. The simulation is performed on anisotropic lattice with $N_f=2+1$ flavours of clover fermion at quark mass which corresponds to $m_\pi=384(4) {\rm MeV} $. The thermal ensembles were provided by the FASTSUM collaboration and the zero temperature ensembles by the Hadspec collaboration. We also calculate the spectral density of the correlation function of those particles. The spectral density distributions show rich peak structure at the lowest temperature, while at intermediate temperatures, the mass values of those particles obtained by extrapolation method reflect a two-peak structure. While the spectral density for octet baryon becomes smooth at $T/T_c = 1.27, 1.52, 1.90$, the spectral density for H-dibaryon becoms smooth at $T/T_c = 1.90$. At $T/T_c =0.24 $, the mass difference of H-dibaryon and $\Lambda$ pair $\Delta m = m_H - 2\,m_{\Lambda} $ is estimated to be $\Delta m = -14.6(6.2) {\rm MeV}$ which suggests there exists a bound H-dibaryon state.


I. INTRODUCTION
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 [1].In contrast with the only known stable dibaryon (deuteron) whose binding energy is about 2.2 MeV, Jaffe predicted that the binding energy of H-dibaryon is about 80 MeV below the ΛΛ threshold 2230 MeV which means that Hdibaryon is a deeply bound state.
The observation of double hypernuclei, A ΛΛ Z, is very important in connection with the existence of Hdibaryon [3].If the mass of H-dibaryon m H is much smaller than the mass of double Λ hyperon 2 m Λ , 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 [3].So the observation of the weak decay of double hypernuclei will put limitation on the mass of H-dibaryon m H > 2 m Λ − B ΛΛ with B ΛΛ being the binding energy of ΛΛ hyperons.
The experiments [3][4][5][6][7][8]11] investigated the nuclear capture of Ξ − 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 ΛΛ [3,5,8,11], or the cross section of enhanced production of ΛΛ pair [4,6,7].The results do not confirm the existence of H-dibaryon, and set the lower limit for the mass of H-dibaryon.
The experiments [9,10] were carried out to search for H-dibaryon or deeply bound singlet uuddss sexaquark S (for the explanation of S, see [2]) in Υ → S ΛΛ 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 Λ , such results support that H-dibaryon is a bound state [12][13][14] , while other quenched studies show that H-dibaryon is not a bound state [15][16][17][18].
Aside from quenched studies, simulations with dynamical fermions have been carried out by NPLQCD, HALQCD and other groups.
The HALQCD collaboration investigated the baryonbaryon interaction in terms of the baryon-baryon potential.They extracted the Nambu-Bethe-Salpeter wavefunction 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 [24][25][26], and on N f = 2 + 1 ensembles [27][28][29].Of the results obtained by the two groups, the simulations on N f = 2 + 1 ensembles [27][28][29] by HALQCD claimed that H-dibaryon may be a ΛΛ resonance.Other results by the two groups agreed on the presence of the H-dibaryon, despite disagreement on the binding energy [30].
Ref. [30] 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. [31] carried out simulations by using O(a)improved Wilson fermions at SU(3) symmetric point with m π = m K ≈ 420 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 [36,37].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, [38][39][40][41][42][43]). 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 [45][46][47][48][49]. 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 Hdibaryon 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 Λ 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. [50].Our simulation results are presented in Sec.IV, followed by discussion in Sec.V.

II. LATTICE CALCULATION AND SETUP
In our simulation, we compute the correlation functions of H-dibaryon and Λ, we also calculate the correlation function of N , Σ and Ξ.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 [30]: where a, b, . . ., f denote generic quark flavors, and P + = (1 + γ 0 )/2 projects the quark fields to positive parity.We choose the operator O H as the H-dibaryon interpolating operator which transforms under the singlet irreducible representation of flavor SU(3) [17,[51][52][53]: The H-dibaryon correlation function can be obtained based on the formulae in Ref. [53].
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 τ in the interval 0 ≤ τ < 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 [τ 1 , τ 2 ] where τ 2 is fixed to the whole time extent, and τ 1 runs over several values from τ 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 [τ 1 , τ 2 ] against 1/τ 1 , and fit a linear expression to those mass values, then, extrapolate the linear expression to τ 1 → ∞.

III. SPECTRAL FUNCTION
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 [57][58][59].The second is the Backus-Gilbert method and their variants [60][61][62][63]( reviews on the spectral function in lattice QCD can be found in Ref. [64,65] and references therein ).Recently, based on the Backus-Gilbert method, a new 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 which was designed 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 written as: with ρ L (E) being the spectral function.We choose the basis function as: We can approximate ρ L (E ⋆ ) by ρL (E ⋆ ), where ρL (E ⋆ ) can be evaluated by after those coefficients g τ (E ⋆ ) are determined.
The coefficients g τ (E ⋆ ) are determined by minimizing the linear combination W [λ, g] of the deterministic functional A[g] and error functional B under the unit area constraint where A[g] is defined as: with ∆σ (E, E ⋆ ) and ∆ σ (E, E ⋆ ) being the smearing function and target smearing function, respectively.These two functions are given by and respectively.B[g] is written as: with Cov is the covariance matrix of the correlation function G(τ ).More details are given in Ref. [50].

IV. MC SIMULATION RESULTS
Before presenting the simulation results, we describe the computation details.The simulations are carried out on N f = 2 + 1 Generation2 (Gen2) FASTSUM ensembles [43] of which the ensembles at the lowest temperature are provided by the HadSpec collaboration [66,67], so the computation details are the same as those used in Ref. [43].We recompile the simulation details in the following three tables I, II, and III from Ref. [43].
The ensembles are generated with a Symanzikimproved 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 in Table I.The N f = 2 + 1 Gen2 ensembles correspond to a physical strange quark mass and a bare light quark mass of a τ m l = −0.0840,yielding a pion mass of m π = 384(4) MeV (see Table II).
The ensemble detail is listed in Table III which is recompiled from Ref. [43] with a slight difference on ensemble N 3 s × N τ = 32 3 × 48.The corresponding physical parameters such as the lattice spacing, and the pion mass etc are collected in Table II.
The quark propagators are computed by using the deflation-accelerated algorithm [68,69].When computing the propagator, The spatial links are stout smeared [71] with two steps of smearing, using the weight ρ = 0.14.For the sources and sinks, we use the Gaussian smearing [70] η where H is the spatial hopping part of the Dirac operator and C an appropriate normalisation [49].
The correlators of Λ and H-dibaryon are presented in Fig. 1 and Fig. 2, respectively.For the correlators of Λ and H-dibaryon, we find similar behavior which was displayed in Fig. 1 in Ref. [49] for N .For the correlator of Λ on large N τ and relatively small N s lattice, especially 24 3 × 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 , Ξ, Σ, Λ, 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 τ = 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 τ = 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 III.We can find that the masses decrease when temperature increases.We compare our results of N and Λ below T c with those in Refs.[45,49].The results are consistent within errors.
We also calculate the spectral density ρL (E ⋆ ) of the correlation function of N , Σ, Ξ, Λ and H-dibaryon by using the public computer program [72].
We present the spectral density with different σ = 0.02, 0.04, 0.06, 0.08 for N , Ξ 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 σ 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 σ = 0.08 has just one position where ρL (E ⋆ ) takes locally maximum value.The position is about at E ⋆ = 0.37.At T /T c = 0.24, the time extent N τ = 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 τ m N which is larger than 0.30.The largest value of a τ 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 IV.So we think taking large σ may lead to missing some peak structure.On the other hand, spectral density ρL (E ⋆ ) obtained by using σ = 0.02 in the upper panel of Fig. 5 has a peak position at E ⋆ ≈ 0.05 with small peak value.This peak structure may be due to the lattice artefact.
The middle panel in Fig. 5 for Ξ at T /T c = 0.95 shows that smaller σ value can make peak structure of spectral density in small E ⋆ region more pronounced.The lower panel for H-dibaryon at T /T c = 1.90 suggests that different σ value has little effect on the computation of spectral density at high temperature.So, we just present the spectral density results computed with σ = 0.020 in the following.
The spectral density ρL (E ⋆ ) of Ξ, Λ and H-dibaryon are given in Fig. 6, 7 and 8, respectively.The spectral density ρL (E ⋆ ) of N and Σ has similar behaviour to that of Λ. From Fig. 6 , 7 and 8, we can find that the spectral extrapolation method are in that range of E ⋆ .However, a τ m H = 0.44 at T /T c = 0.24 for H-dibaryon is in the neighbour of peak position E ⋆ = 0.35 where the ρL (E ⋆ ) value is not very large.Obtaining a τ 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 Ξ and Λ 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 ρL (E ⋆ ) has a two-peak structure.If we take the smaller values of E ⋆ at peak positions as the ground state energies of corresponding particle, then these mass values obtained by the peak position of ρL (E ⋆ ) are smaller than those mass values obtained in Ref. [49] and Ref. [45].Mass values of a τ m Ξ , a τ m Λ and a τ m H presented in Table.III are not consistent with the peak positions of cor- responding 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 δ function structure contributing to the correlation function.We can find this observation from Fig. 9 for N .Similar behaviour can be found for Σ, Ξ and Λ.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.

V. DISCUSSIONS
We have made a simulation in an attempt to determine the masses of the conjectured H-dibaryon with 2 + 1 fla- 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) MeV [43,73].
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 III.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 IV 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 ⋆ becomes smooth gradually when σ increases.Considering the quark mass which corresponds to m π = 384(4)MeV, if we take the smaller E ⋆ 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 , Σ, Ξ and Λ which are presented in table III 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 multibaryon state, this conclusion awaits further investigation.
It is appropriate to consider the lowest temperature ensembles N τ = 128 to be the zero temperature ones since N τ > ξN s [43].Using the mass values of H-dibaryon and Λ in Table III at N τ = 128, T /T c = 0.24, an estimation of ∆m = m H − 2m Λ can be made.∆m = m H − 2m Λ = −0.0026(11) which is converted into physical unit to be ∆m = m H − 2m Λ = −14.6(6.2)MeV.
In our simulation, the correlators of proton, Λ and Hdibaryon at N τ = 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 π = 384(4) 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.

FIG. 1 :FIG. 2 :
FIG.1: Euclidean correlator G(τ )/G(0) of Λ as a function of τ T at different temperatures.At the lowest temperature T /Tc = 0.24, the correlators at some points are not displayed due to minus values.

1 FIG. 3 :
FIG. 3: Mass values of nucleon and H-dibaryon obtained by fitting equation (5) to correlators on Nτ = 128 ensembles.The mass values are the fitting parameter m+ in equation 5 extracted by fitting procedure in different interval [τ1, τ2].Horizontal axis label τ1 represents different number of time slices suppressed which corresponds the lower bound of interval [τ1, τ2].

FIG. 5 :
FIG.5: Spectral density computed at different parameter σ for the target smearing function ∆σ(E, E⋆) for H-dibaryon, Ξ and N at different temperatures.

TABLE I :
[43]meters in the lattice action.This table is recompiled from Ref.[43].

TABLE III :
[66,67] and temporal extent, temperature in MeV, number of configurations, mass of N , Σ, Ξ, Λ and H-dibaryon.Masses of baryons and H-dibaryon are obtained by extrapolation method.Estimates of statistical and systematic errors are contained in the first and second brackets, respectively.The errors of fitting parameters of linear extrapolation are taken to be the systemtatic errors of the masses.The ensembles at the lowest temperatures were provided by HadSpec[66,67](Gen2).

TABLE IV :
On different Nτ lattice, peak position E⋆ of spectral density for N .

TABLE V :
On different Nτ lattice, peak position E⋆ of spectral density for Σ.

TABLE VI :
On different Nτ lattice, peak position E⋆ of spectral density for Ξ.

TABLE VII :
On different Nτ lattice, peak position E⋆ of spectral density for Λ.

TABLE VIII :
On different Nτ lattice, peak position E⋆ of spectral density for H-dibaryon.