Ground state baryons in the flux-tube three-body confinement model using Diffusion Monte Carlo

The quark model is widely used to study the hadron mass spectrum. The quarks inside a hadron are modeled as the constituent quarks interacting via the effective potentials. Various quark potential models have been used for the conventional and exotic hadron spectrum over the years (see Refs. [1, 2] for the quark model reviews and see Refs. [3, 4, 5, 6, 7] for the reviews about the exotic hadrons.). The quark level interactions usually include the color-dependent Coulomb interaction, spin-dependent chromo-magnetic interaction, tensor interaction, and spin-orbit interaction. Basically, the above interactions can be derived from the one-gluon-exchange mechanism [8]. In addition, there is a confinement term to describe the long-range interaction. Quark potential models behave nicely to describe the conventional hadron spectrum by solving the two-body or three-body problems. For instance, the Coulomb-plus-Linear Cornell potential proposed by Eichten et al. can well reproduce the charmonium and bottomonium spectra [9, 10, 11]. A relativized quark model constructed by Isgur and his collaborators works successfully for all mesons and baryons [8, 12]. After that, Semay and Silvestre-Brac found that if the parameters are chosen correctly, the non-relativistic approach could accurately simulate the spectra from the relativistic one [13], and they built a new non-relativistic potential that works equally well in the meson and baryon sectors [14, 15].

However, the confinement mechanism in quantum chromodynamics (QCD) is still elusive and ambiguous, which is usually introduced phenomenologically [16] or inspired by lattice QCD (LQCD) simulations [17]. For the $q\bar{q}$ mesons, the confinement potential is apparently pairwise. But for the $qqq$ baryons, there is no compelling reason to assume the confinement interaction to be a two-body one as shown in the left panel of Fig. 1, which is called the $\Delta$-type potential. Actually, another Y-type interaction was suggested by Artru in the string model for the baryons [18]. In this scheme, three strings are connected at one point and carry quarks at their ends, as shown in the right panel of Fig. 1. This Y-type confinement mechanism has been investigated in different models, such as the QCD bag model [19], non-relativistic quark model [20, 21], semi-relativistic quark model [22, 23, 24, 25], and relativized quark model with chromodynamics [12]. A comparison between the $\Delta$-type and Y-type interactions was made by LQCD, and its fitting results seemed to favor the Y-type one [26, 27].

Nevertheless, this Y-type potential is in fact difficult to calculate with the conventional variational principle. If one expresses the minimum total length $L_{\mathrm{min}}$ of the flux tubes in terms of the three-body Jacobi coordinates, it will turn into a form containing the square-root as well as some complicated angle-dependence, making the evaluation of the integral in the matrix elements a difficult task [28]. And certainly, extending the multi-body confinement to the multiquark states is a tougher work [29]. Besides, there are other shortcomings in the framework of the variational method. For instance, when the number of particles (equivalently the dimensions of the wave function) is large, the number of the basis increases exponentially, so does the matrix dimension. Thus meeting the demand for the storage space and computing resources could become a challenge [30].

A promising alternative of the variational method is the diffusion Monte Carlo (DMC) method. In this formalism, the distribution of the so-called walkers is used to represent the spatial wave function, which will gradually evolve to the ground state over time in principle. In this scheme, the spatial integrals are replaced by summations over walkers. Consequently, there is no extra complexity in handling the Y-type interaction compared with the $\Delta$-type interaction. Meanwhile, the uncertainties of the Monte Carlo methods decrease as $1/\sqrt{N}$, where $N$ is the number of walkers. This behavior is independent on the dimension of the wave function, which is a promising advantage against most deterministic methods that depend exponentially on the dimensions. As for the multi-particle problem, the DMC method is easily parallelized to carry out the high-performance simulations [31, 32]. The DMC has been well used in the simulations of the molecular physics [33], solid physics [34] and nuclear physics [35].

In hadronic physics, the DMC method has been used in quark models in several works. Bai et al. calculated the $bb\bar{b}\bar{b}$ tetraquark ground state energy in a non-relativistic model with the flux-tube confinement potential [36]. A $0^{++}$ bound state with a mass of 18.69 GeV was predicted, which is 100 MeV below the $\eta_{b}\eta_{b}$ threshold. The $0^{++}$ $bb\bar{b}\bar{b}$ system was also calculated by Gordillo et al. using a similar DMC method [37], with the pairwise confinement interaction. They gave a different mass of 19.199 GeV, which is 300-400 MeV above the $\eta_{b}\eta_{b}$ and $\Upsilon(1S)\Upsilon(1S)$ thresholds. In addition to the $bb\bar{b}\bar{b}$ system, they also calculated other fully-heavy tetraquarks in Ref. [37]. The predicted masses are all above the corresponding meson-meson thresholds, and agree with the results of the variational method in Ref. [38]. The difference between the fully-heavy tetraquark states of these two DMC calculations could arise from the different interactions, especially the confinement part, the different color configurations assumed (the diquark-antidiquark or di-meson configurations ) or the details of the DMC algorithm. Recently, the DMC method was also used to investigate other multiquark systems [39, 40, 41]. Experimentally, an analysis using data from the CMS detector reported an 3.6$\sigma$ enhancement at 18.4 GeV in the invariant mass distribution of the $\Upsilon(1S)\ell^{+}\ell^{-}$ final state [42, 43], which might be a candidate for the $bb\bar{b}\bar{b}$ tetraquark state. But this enhancement was not seen by LHCb [44]. Recently, the LHCb [45], CMS [46] and ATLAS [47] collaborations reported the observation of several resonance structures in the di-$J/\psi$ or $J/\psi\psi(2S)$ channels. All of these structures are above the $J/\psi J/\psi$ and $\eta_{c}\eta_{c}$ thresholds.

Theoretically, in the variational method, it was shown that the tetraquark states above the di-meson thresholds in Ref. [38] will become either the scattering states or resonances if more quark-clustering configurations are considered [48]. However, the primitive DMC is a method to calculate the ground bound states without quark-clustering beforehand (In order to calculate the excited states [33] and resonances [49, 50] with DMC, other nontrivial ingredients should be included.). For the tetraquark system, if there exists no bound state solution, one should expect that the wave function in DMC evolves exactly to the corresponding meson-meson threshold, rather than a state above the threshold as in Ref. [37]. Therefore, the doubt arises whether the DMC really avoids the usual quark-clustering and provides an exact estimate of the ground state energy as well as the wave function.

Considering the above issues, a simpler system is needed as a benchmark to test the DMC method and explore the reason why the lowest state cannot be obtained in some cases. Though the DMC method has been used for the nucleons, electrons etc., the multiquark systems are different because of the color confinement. The baryon system is a suitable platform to examine the DMC method, because of the following reasons: (i) the baryon is a bound state, having no ambiguity with the resonance; (ii) the color component is simple. There is only one possible color configuration; (iii) the experiment results are precise; (iv) the flux-tube confinement can be included; (v) although quarks are Fermions, the baryon is equivalent to a boson system, which avoids the notorious sign problem. This can be easily seen from its wave function. The baryon total wave function is the product of the spatial, spin, flavor and color parts:

For the baryons, the only possible color configuration is that any two quarks form a color $\bar{3}_{c}$ representation and then combine with the third quark in the $3_{c}$ representation to become a color singlet. In this way, the color wave function is fully anti-symmetric and the color factor $\langle\chi_{\mathrm{color}}|\frac{\lambda_{i}}{2}\cdot\frac{\lambda_{j}}{2}|\chi_{\mathrm{color}}\rangle$ for any $(ij)$ quark pair is always $-\frac{2}{3}$. Thus if there are identical particles in the system, the exchanging symmetry renders the remaining $|\phi_{\mathrm{spatial}}\rangle|\chi_{\mathrm{spin}}\rangle|\chi_{\mathrm{flavor}}\rangle$ part to behave like a boson system. Therefore, there is no Fermionic sign problem in the DMC calculation [51].

The fully-heavy baryon system has been calculated by Gordillo et al. [37]. However, the baryon systems containing the light quarks are more interesting, especially the nucleon system. Hopefully the discussions about the baryons can also give enlightenment on the understanding of the threshold problem of the tetraquark system and provide some hints for the future explorations.

This paper is arranged as follows. In Sec. II, the diffusion Monte Carlo method is introduced, including the single-channel and coupled-channel formalisms, respectively. In Sec. III, the quark model Hamiltonian with different types of the confinement interaction is presented. In Sec. IV, the numerical results for all ground state baryons are given. The results among different calculation methods and different confinement potentials are compared. In Sec. V, we discuss the tetraquark threshold problem, clustering problem, and give some prospects on further applications of the DMC method in the field of hadron physics. Finally, a brief summary is given in Sec. VI. In Appendix A, we give the formalism of the convection–diffusion equation. In Appendix C, we illustrate the constituent quark size contribution to the charge radius.

The nonrelativistic Hamiltonian of a 3-quark system reads

where $m_{i}$ and $\bm{p}_{i}$ are the mass and momentum of quark $i$. $T_{CM}$ is the center-of-mass kinematic energy, which automatically vanishes in the evolution, because the system will tend to the lowest energy state.

In order to investigate the effects of different confinement scenarios, we modify a nonrelativistic potential proposed in Ref. [15]. The potential is divided into two parts, $V=V_{0}+V_{\text{conf}}$. For different scenarios, we choose the same $V_{0}$ terms,

with

where $\lambda_{i}$ is the $\text{SU}(3)$-color Gell-Mann matrix, $\bm{s}_{i}$ is the spin operator of quark $i$, and $r_{ij}$ is the relative distance between quark $i$ and $j$. In the above interaction, the Coulomb term and hyperfine term come from the one-gluon exchange interaction. The $\Lambda$ term is a pairwise constant interaction to shift the overall mass spectrum. The $C$ term is a three-body phenomenological contribution to obtain a better agreement with the experimental baryon masses. In Ref. [15], the parameters of the potential have been determined through fitting to a large number of mesons, and here we use the $AL$1 model parameters, which are listed in TABLE 1.

In this work, we investigate two different confinement scenarios, as shown in Fig. 1, the pairwise ($\Delta$-type) confinement and the three-body flux-tube (Y-type) confinement. In the first one, the two-body linear confinement term is introduced as

where $\lambda$ is the confinement coupling constant in the $AL$1 model of Ref. [15]. The value of $\lambda$ is presented in TABLE 1. For the baryons, we introduce an alternative coupling constant $\sigma_{\Delta}$ to replace $\lambda$. Apparently, the confinement potential is proportional to the length of the perimeter of the triangle connecting three quarks.

In the second scenario, the confinement interaction is proportional to the minimal total length of the color flux tubes linking three quarks,

where $\sigma_{Y}$ is the string tension. $L_{min}$ denotes the minimal length by tuning the joint point. In this scenario, three quarks are confined by the three-body force. In Refs. [26, 27], the Y-type confinement interaction was favored by the lattice simulations. Another point of view is that the confinement in baryon is roughly one-half $\Delta$ and one-half Y-string [60, 61], while it is not explored in this work.

We will take two $\sigma_{Y}$ values and give their results respectively later in this work. According to the lattice QCD result in Ref. [26], a universal feature of the string tension $\sigma_{Y}\simeq\sigma_{Q\bar{Q}}$ is found. Herein for the first $\sigma_{Y}$ value, we naively take $\sigma_{Y}=\sigma_{Q\bar{Q}}=2\sigma_{\Delta}=\lambda$. However, the above value is only a rough approximation. To get a more accurate $\sigma_{Y}$ value, we adjust it to make the $\Omega^{-}(sss)$ mass coincide exactly with the experimental mass $1.672$ GeV. Since we expect the confinement coupling constants for the light quark systems and heavy quark systems could be different, the system composed of strange quarks (too heavy to be light and too light to be heavy) could be a good compromise. The benchmark gives us $\sigma_{Y}=0.9204\sigma_{Q\bar{Q}}=0.9204\lambda$. This coefficient is consistent with the best fitting parameters in lattice QCD simulation [26] with $\sigma_{Y}/\sigma_{Q\bar{Q}}=0.1524/0.1629=0.9355$.

For the three quark system, the minimal value of the total length $L_{min}$ linking three quarks can be obtained analytically [27]. In this work, we choose a more general numerical algorithm to determine $L_{min}$ in order to extend the framework to the multiquark systems [36] in the future. This operation is essentially a Euclidean Steiner Tree Problem (ESTP) [62]. The ESTP seeks a network of the minimal length spanning a set of points by allowing the insertion of new points (Steiner points). The Smiths algorithm was designed for ESTP, which is an iterative method to optimize the search of coordinates of Steiner points in $d$-dimensional space, given a topology and terminal positions. The code of this program can be found in Ref. [63].

In the simulation, we use $1\times 10^{4}$ walkers to sample the $f(\bm{R},t)$ or $\mathcal{F}(\bm{R},t)$. We let the ensemble evolve $1\times 10^{4}$ steps for the single-channel system and $2\times 10^{4}$ steps for the coupling-channel systems with the $\Delta t=0.01\mathrm{GeV}^{-1}$ for each step to ensure stability. The resulting energy is averaged over the last 5000 steps to reduce fluctuation. To estimate the statistical uncertainty, the correlations among adjacent steps should be considered. To do this, we divide the steps into blocks and calculate the block averages. When the block size is taken large enough, the averages become uncorrelated and the uncertainty becomes independent of the block size. We find the blocks have become uncorrelated when taking the size as 500 steps. So we divide the last 5000 steps into 10 blocks of size 500 steps. Then by using the Jackknife resampling method [64], the uncertainty turns out to be less than 1 MeV. As an example, the uncertainty of $\Omega^{-}(sss)$ mass is 0.3 MeV. More details are given in Appendix B.

Once the system is stabilized, the wave function will not change with time. We calculate the radial distribution $r_{ij}^{2}\rho(r_{ij})$, mean-square radius and charge mean-square radius using forward walking algorithm introduced in Sec. II.3. The $\rho(r_{12})$ is defined as

The square of the wave function $|\Psi(\bm{r}_{1},\bm{r}_{2},\bm{r}_{3})|^{2}$ can be obtained from Eq. (II.3). The definitions of $\rho(r_{13})$ and $\rho(r_{23})$ are similar. The root mean-square radius is defined as

where $\bm{r}_{i}$ indicates the position of the $i$th quark. For the point-like quark, the charge mean-square radius $R_{c}^{2}$ is defined as

where $\bm{R}_{CM}$ refers to the center of mass, and $e_{i}$ is the charge of the $i$th quark. In the practical calculation, the first term can be obtained from Eq. (II.3). The second term represent the effect of the charge radii of the constituent quarks. It is unreasonable to naively regard the constituent quarks as point-like particles. Details are given in Appendix C.

In this work, the charge radii of the constituent quarks are extracted from experiments and LQCD calculations. The contributions from the $u$ and $d$ quark size are extracted from the $p$ and $n$ experimental $\langle R_{c}^{2}\rangle$ values [65]. In principle, one can extract the charge radius of the constituent $s$ quark from the experimental $\langle R_{c}^{2}\rangle$ of $\Sigma^{-}$ [66]. However, the present experimental uncertainty is large. Hence we choose to extract the $s$ and $c$ charge radii from the $\langle R_{c}^{2}\rangle$ of $\Omega(sss)$ and $\Omega(ccc)$ through lattice QCD simulations in Refs. [67, 68]. These contributions are listed in TABLE 2. We can see that the quark charge radii vary with flavors, which is different from the scenario in Refs. [69, 70]. Meanwhile, their signs are consistent with their charges. The electric size of the constituent quark decreases with the quark mass. In view of the small contribution of the $c$ quark size, we neglect the contribution from the $b$ quark.

The distribution plots related to the wave function are averaged over 500 steps to reduce the noise and make it smooth. And the expectation value for the observables are averaged over 2000 steps.

In order to show the internal structure of the quarks more visually, we provide two additional plots, the angle distribution and rotation-irrelevant distribution following Ref. [71]. For each walker, we introduce $\theta_{1}$, $\theta_{2}$ and $\theta_{3}=180^{\circ}-\theta_{1}-\theta_{2}$ to label the three inner angles of the triangle by linking three quarks. With the spatial probability distribution of walkers, we can easily get the angle distributions and the expectation values of the inner angles. To visualize the rotation-irrelevant structure, we define a way to eliminate the rotation degree of freedom. For each walker, we first put the triangle $R_{1}R_{2}R_{3}$ connecting three quarks into a 2D $x-y$ plane as shown in Fig. 2. We put its center of mass at the origin $O$. Now the only remaining degree of freedom for each walker is the rotation around the origin $O$. We then use the triangle $ABC$ formed by three inner angle expectation values to define a reference frame, where $A,B,C$ corresponds to $1,2,3$ quark respectively, and vertex $C$ is fixed on the y-axis. Finally we rotate the triangle $R_{1}R_{2}R_{3}$ to minimize the quantity $\angle R_{1}OA^{2}+\angle R_{2}OB^{2}+\angle R_{3}OC^{2}$. We can draw the new positions of quarks for all walkers in the plot and get a distribution reflecting the inner structure of the quarks. In principle, we can choose other angle-fixing strategies, for example minimizing $\angle R_{1}OA^{4}+\angle R_{2}OB^{4}+\angle R_{3}OC^{4}$. The different strategies will not change the qualitative properties. We will choose some typical baryons to show their angle distributions and rotation-irrelevant distributions. The complete distributions will be given in the Supplement material.

In the following sections, we will classify all of the ground state baryons according to the number of identical quarks inside them. We treat the $u$ and $d$ quarks as identical particles under SU(2)-flavor symmetry. We will use $n$ to label them. The related spin matrix elements are presented in TABLE 3. For the the $AL$1 model, apart from the DMC, we use a variational method with $20^{2}$ Gaussian basis functions [72] as a benchmark, which is denoted as GEM in the following. Our results are also compared with the Faddeev formalism results in Ref. [15]. Furthermore, for the flux-tube confinement model, we use two sets of parameters as mentioned above. Flux-tube I refers to the universal tension parameter with $\sigma_{Y}=\lambda$, and Flux-tube II refers to $\sigma_{Y}=0.9204\lambda$ determined by the experimental $\Omega^{-}(sss)$ mass. Both of them are solved using DMC algorithm.