Lattice QCD studies on decuplet baryons as meson-baryon bound states in the HAL QCD method

We study decuplet baryons from meson-baryon interactions in lattice QCD, in particular, $\Delta$ and $\Omega$ baryons from P-wave $I=3/2$ $N\pi$ and $I=0$ $\Xi\bar{K}$ interactions, respectively. Interaction potentials are calculated in the HAL QCD method using 3-quark-type source operators at $m_{\pi} \approx 410~\textrm{MeV}$ and $m_{K} \approx 635~\textrm{MeV}$, where $\Delta$ as well as $\Omega$ baryons are stable. We use the conventional stochastic estimate of all-to-all propagators combined with the all-mode averaging to reduce statistical fluctuations. We have found that the $\Xi\bar K$ system has a weaker attraction than the $N\pi$ system while the binding energy from the threshold is larger for $\Omega$ than $\Delta$. This suggests that an inequality $m_{N}+m_{\pi}-m_{\Delta}<m_{\Xi}+m_{\bar K}-m_{\Omega}$ comes mainly from a smaller spatial size of a $\Xi \bar K$ bound state due to a larger reduced mass, rather than its interaction. Root-mean-square distances of bound states in both systems are small, indicating that $\Delta$ and $\Omega$ are tightly bound states and thus can be regarded qualitatively as composite states of 3 quarks. Results of binding energies agree with those obtained from temporal 2-point functions within large systematic errors, which arise dominantly from the lattice artifact at short distances.


I. INTRODUCTION
Most of hadrons have been understood well in the quark model, while exceptions, called exotic hadrons, were found in various experiments recently [1]. Exploring properties or internal structures of such exotic hadrons from QCD is one of the biggest issues in hadron physics. Since exotic hadrons typically appear as resonances due to non-perturbative QCD interactions, theoretical studies from first-principles lattice QCD are mandatory.
Properties of hadron resonances such as masses and decay rates are evaluated from hadron scatterings. The conventional method to analyze hadron scatterings in lattice QCD is the finite-volume method [2,3], which relates energies on finite volume(s) to scattering amplitudes on the infinite volume. An alternative method is the HAL QCD method [4][5][6], in which we extract interaction potentials directly in lattice QCD and then obtain scattering amplitudes from potentials by solving the Schorödinger equations in the infinite volume.
This method is particularly advantageous to study systems involving baryons [7].
As a first step toward understandings of exotic hadrons including pentaquarks, we focus our attention on decuplet baryons, spin 3/2 baryons symmetric under quark flavor exchanges, since all decuplet baryons except Ω appear as resonances. There are studies on decay properties of the process ∆ → N π by extracting the transfer matrix elements in lattice QCD [8,9], and are also several studies on the N π scatterings for the ∆ baryon in the finitevolume method at the lighter quark masses [10][11][12][13][14], where signals of the ∆ as a resonance are observed.
In this paper, as a first analysis for decuplet baryons in the HAL QCD method, we investigate a question why Ω appears as a stable particle below the ΞK threshold while others such as ∆ become resonances. Historically, the decuplet baryons have been studied from the 3-quark state picture in flavor SU(3) symmetry. They belong to a 10-dimensional representation with similar masses and the mass splittings can be also explained by the Gell-Mann-Okubo mass formula. In this picture, the difference between Ω and ∆ seems to be explained by the quark mass (m q ) dependence of the meson and baryon thresholds; the pseudoscalar meson masses are proportional to √ m q while the baryon masses linearly depend on m q . On the other hand, the difference between Ω and ∆ become more non-trivial once we investigate these states from the scattering theory with meson-baryon degrees of freedom (dof). This approach, however, could have a broader utility since it can handle a bound state and a resonance on equal footing and applications to exotic hadrons are possible. We thus reexamine the decuplet baryons from the aspects of meson-baryon dof by first-principles lattice QCD calculations in this study. To this end, we extract the ΞK potential for Ω and the N π potential for ∆ in the HAL QCD method, and make a comparison between them. To reduce huge computational costs required for these scattering channels, we employ heavier quark masses, where u, d quark masses are chosen at the value close to the s quark mass with slightly broken SU(3) flavor symmetry. While both Ω and ∆ appear as stable particles in this setup, an inequality m N + m π − m ∆ < m Ξ + mK − m Ω still holds as seen in previous lattice QCD results [15], and we can investigate the rephrased question, "what is the physical origin which brings this hierarchy?". As we will show later, the HAL QCD method is particularly useful for this purpose, since it can directly extract the interaction potentials and distinguish two possible origins of the hierarchy, one from the difference in interactions and another from difference in kinematics.
This paper is organized as follows. In Sec. II, we briefly review the HAL QCD method in meson-baryon systems. We define N π and ΞK 3-point correlation functions and their radial/spherical decompositions in Sec. III, and show simulation details in Sec. IV. In Sec. V, we present potentials for N π and ΞK, explain the fitting procedure of the potentials, and show numerical results of observables such as scattering phase shifts, binding energies, and root-mean-square distances. Sec. VI is devoted to a conclusion of this paper. Some technical details and additional contents related to our study are discussed in appendices.

II. HAL QCD METHOD IN MESON-BARYON SYSTEMS
We define a meson-baryon potential U αβ (r, r ) in the HAL QCD method [4,5] as where W is a total energy, k is a relative momentum in the center of mass frame, α, β are indices for upper spin components, µ is a reduced mass and H 0 is a free Hamiltonian. The total energy W is related to k as W = k 2 + m 2 M + k 2 + m 2 B , where m M and m B are meson and baryon masses, respectively. An equal-time NBS wave function Ψ W α (r) is defined by where |0 is a vacuum state in QCD, M (x, t) and B α (x, t) are meson and baryon operators at spacetime (x, t), respectively, and |M B; W is a meson-baryon state with energy W .
The potential U αβ (r, r ) in Eq.(1) is energy-independent and non-local. Once the potential U αβ (r, r ) is obtained, we can extract an S-matrix for the meson-baryon scattering by solving the corresponding Schrödinger equation.
In our study, we employ the time-dependent HAL QCD method [6], an improved version of the original HAL QCD method, to extract potentials in lattice QCD. In this method, we first define an R-correlator as where C M (t) and C B (t) are 2-point correlation functions for meson and baryon, respectively, and F is a n-point correlation function (n > 2) of the meson-baryon system, which is given whereJ M B (t 0 ) is the source operator at timeslice t 0 , which creates meson-baryon states with the quantum numbers of interest from the vacuum. The R-correlator can be decomposed into contributions from elastic and inelastic states as where W n is the energy of the nth eigen states, A n is a coefficient independent of r and α, and ∆W n = W n − m M − m B is an energy difference from the threshold. Eq. (1) implies that each term in the elastic part, A n Ψ Wn α (r)e −∆Wnt , satisfies where k 2 n /2µ can be expressed in terms of ∆W n as with M = m M + m B and We thus express k 2 n by an expansion in terms of ∆W n as where δ = (m M − m B )/M . Rewriting ∆W n in this series in terms of time derivatives and summing over n in Eq.(6), we obtain for a large enough t to suppress inelastic contributions.
Applying the Okubo-Marshak expansion [16] to meson-baryon systems, the leading order (LO) term of U αβ (r, r ) in the derivative expansion is given by where V LO (r) can be extracted from R α (r, t) for any α as We truncate an infinite summation over k by k ≤ 2 for N π and k ≤ 3 for ΞK, respectively, since remaining higher order contributions are negligibly small 1 .

TORS
In order to investigate interactions of P-wave I = 3/2 N π and I = 0 ΞK systems in the HAL QCD method at low energies, we use following 3-point correlation functions, where sink operators are defined by Since we expect bound decuplet baryons to appear below thresholds, we employ 3-quark-type decuplet baryon operators at the source, where the explicit forms are given by with Γ ± = 1 2 C(γ 2 ± iγ 1 ) and Γ z = −i √ 2 Cγ 3 , and q = (u, s) for D = (∆ ++ , Ω − ), respectively. Indeed each operator in Eq. (18) couples to a spin-3/2 particle with a different z component of the spin, since they belong to an H g irreducible representation of the cubic group O D h . A summation over z selects zero total momentum states.
To obtain NBS wave functions with J P = 3/2 + , we project F α,jz (r, t) onto the same component in the H g representation ofD jz (t 0 ). Using Clebsch-Gordan coefficients, projected 3-point correlation functions can be described as where Y l,m (r) is the spherical harmonics and f jz (r, t) is a factor that depends only on r = |r| and t. We extract f jz (r, t) using a projection to (l = 1, m), defined on a discrete space as with corresponding (m, α) for each j z . For j z = ±1/2, we can derive f jz (r, t) in two ways by setting either (m, α) = (0, ↑) or (1, ↓) for j z = +1/2 and either (m, α) = (0, ↓) or (−1, ↑) for j z = −1/2, respectively. In this paper, we take an average over the factors calculated from the two choices.
Applying Eq. (19) to Eq.(3) and Eq (12), we obtain with the angular momentum l = 1 and In this study, we use this equation to extract the LO potentials. Since f jz (r, t) for any j z gives the same LO potential thanks to a rotation symmetry, we take an average over j z to increase statistics. Furthermore, a charge conjugation symmetry provides a relations among so that an average of f jz (r, t) over j z is guaranteed to be pure imaginary, and we therefore ignore its real part. We employ a smeared quark source with an exponential smearing function [17], which is given by for the strange quark. We also apply the same smearing to quarks at the sink with (A, B) = (1.0, 1/0.7) to reduce singular behaviors of potentials at short distances [18]. Details of singular behaviors and associated issues for its fitting are explained in Appendix A.
We show quark contraction diagrams corresponding to Eq. estimate them approximately, together with dilutions [19] for color/spinor/time components and the s2 (even/odd) dilution [20] for the position z.
For the spatial coordinate x at the sink, the calculation is performed with fixed x,  agators which requires smaller numerical cost. In order to increase the statistics using the translational invariance, we calculate at 64 different spatial points given by x + ∆x with ∆x = (0, 0, 0), (0, 0, 8), · · · , (24, 24, 24). As a method which keeps this computational cost moderate, we employ the truncated solver method [21] combined with the covariantapproximation averaging [22], namely an all-mode averaging (AMA) without low modes.
For the specific value of x, we choose it randomly for each gauge configuration to keep the exact covariance of the AMA, which might be broken by round-off errors. See Appendix C in Ref. [22]. Masses of single hadrons obtained from them in this setup are listed in Table I. Since ∆ and Ω masses lie below N π and ΞK threshold energies, respectively, they appear as bound states in this setup.
In Fig. 2, we present LO potentials for P-wave N π at t = 8-10 and ΞK at t = 8-11, where effective masses of ∆ and Ω reach plateaux, respectively. We do not observe significant t-dependence of potentials, indicating that inelastic contributions as well as effects of higherorder terms in the derivative expansion of the potential are well under control.
Both potentials have very strong attractions at short distances, which can be responsible to produce bound states corresponding to ∆ and Ω. We also observe that the attraction of the N π potential is deeper than that of the ΞK by a few thousands MeV at short distances.
Both potentials have similar shapes at middle and long distances, where meson-baryon interactions may be dominated by one meson exchanges. One possibility is that the N π system exchanges a ρ meson while the ΞK system exchanges a ρ meson or an octet part of Although both potentials go to zero within errors at long distances, interaction ranges are longer than a half of the box size, L/2 ≈ 1.45 fm. We therefore carefully check possible finite-volume effects on observables as will be explained below.

B. Systematic uncertainties of the potentials and fitting results
We estimate systematic uncertainties in our analysis as follows.
Finite-volume effects are estimated by using two types of fit functions, one is a simple 2 Note that this discussion is only qualitative because even the Compton wave length of the lightest pion is 0.5 fm in this setup.
We estimate effects of the leading order truncation for non-localities, by employing potentials at different t, from t = 8 to 10 for N π and from t = 8 to 11 for ΞK, since t dependences of potentials are mainly caused by the truncation of the derivative expansion.
Finite lattice spacing effects are estimated by two approaches. In the first approach, we evaluate the difference between the laplacian term in Eq. (12) calculated with 2nd and 4th order accuracies, whose explicit forms are given by In the second approach, we estimate the uncertainty from the difference between fits with and without data at r = a. For the fit including data at r = a, we find that modification for fit functions is necessary. In fact, if we naively fit with Eqs. (25), (26)   In summary, for potentials at a given t, there are 2 × 2 × 2 = 8 combinations of fitting schemes, which are listed in Table.II.
Central values of physical observables are calculated from Fit 2 at t = 9 for N π and at t = 10 for ΞK. Fit parameters are listed in Table III,   where the latter error is estimated from other fitting schemes and t dependences. Note that systematic uncertainties at short distances, mainly caused by finite lattice spacing effects, are much larger than those at long distances caused by finite-volume effects.

C. Phase shifts, binding energies, and root-mean-square distances
We solve the Schrödinger equation in the radial direction as where k is the relative momentum, V (r) is the fitted potential, and the angular momentum is fixed to l = 1 for the P-wave scattering. The total energy E is related to k 2 as E = k 2 + m 2 M + k 2 + m 2 B . We then extract a scattering phase shift from its solution. Finally, by varying k 2 (or E), we determine an energy dependence of the scattering phase shift.
where first and second errors represent statistical and systematic errors, respectively. As The binding energy in the ΞK system is larger than that in the N π system by more than a hundred MeV. Combined with the previous observation that the ΞK potential has weaker attraction than the N π, it is indicated that a reason for a larger binding energy in the ΞK than in the N π is not the difference of the two interactions.
We also calculate root-mean-square distances of wave functions obtained via the GEM for bound states, which represent the distances between meson and baryon components in ∆ and Ω baryons 3 . We obtain where first and second errors represent statistical and systematic errors, respectively. Results are also shown in Fig. 6 (horizontal axis). Sizes of both bound states estimated by rootmean-square distances are quite small and similar to ranges of attractive pockets in their potentials. These observations suggest that ∆ and Ω are tightly bound states at this quark mass, and can be regarded qualitatively as composite states of 3 quarks rather than mesonbaryon molecule states, which are already known from the phenomenological studies.
More quantitatively, however, their root-mean-square distances are larger than a range of the attractive pocket r ≈ 0.2 fm in both N π and ΞK potentials with the centrifugal term.
(See Fig. 4.) This can be explained by shapes of wave functions, shown in Fig. 7, which exhibits peak structures at r ≈ 0.2 fm as well as long-range tails. Table. IV represents the expectation value of each term in the Schrödinger equation Eq. (28). By comparing ΞK and N π systems, the kinetic term for ΞK is positively larger than that for N π, while the potential term for ΞK is negatively larger. Since the latter effect is larger than the former, the sum of the kinetic and potential terms (= k 2 /2µ ) for ΞK is negatively larger than that for N π. These behaviors can be intuitively understood in the following way: Since the ΞK system has a larger reduced mass, the system tends to be squeezed into a smaller size, which leads to larger effect from the attractive potential at short distances (but with a larger kinetic energy). In fact, in Fig. 7, we observe that the ΞK system has a larger peak of the wave function at short distance, which gives a negatively larger contribution of the potential energy, d 3 r ψ * (r)V (r)ψ(r).
These observations lead to an interesting picture on how these two systems evolve toward smaller quark masses. Assuming that the potentials are not so sensitive to the quark masses, the bound state disappears and ∆ becomes a resonance in the N π system because of its small reduce mass and the broad structure of ∆, while the bound state Ω remains in the ΞK system due to its large reduced mass leading to the compact size of Ω. This scenario may explain spectra observed in Nature.
Results of binding energies and root-mean-square distances suffer from quite large systematic uncertainties compared to statistical errors. Such large uncertainties come dominantly from lattice artifacts at short distances. Fig. 8 and Fig. 9 shows binding energies and root-mean-square distances, respectively, estimated for different fitting schemes and t.
A dependence on an order of the approximation for the Laplacian (2nd/4th for Fit 1/2, 3/4, 5/6, 7/8), or a dependence on a treatment of data at r = a (without/with for Fit   1/3, 2/4, 5/7, 6/8), are much larger than a dependence on t and a dependence on fit without/with mirrors (Fit 1/5, 2/6, 3/7, 4/8). As seen in Fig. 10, the precision of the Laplacian term affects potential data at short distances. (Compare, for example, magenta and cyan points.) Moreover, fit results without and with data at r = a (without/with for Fit 1/3, 2/4) deviate from each other around the origin. Thus, large dependences of binding energies and root-mean-square distances on the accuracy of the Laplacian and the treatment of the data at r = a are indeed caused dominantly by lattice artifacts at short distances, which is associated with rapid changes of potentials around the origin. On the contrary, the results are not so sensitive to finite-volume effects as seen from comparisons between Fit i and Fit i + 4 in Fig. 8 and Fig. 9. potentials.
The ΞK potential has a weaker attraction than the N π while the binding energy of the ΞK is larger than that of the N π, which indicates that the difference between the two interactions is not a reason for the larger binding energy of the ΞK, in other words, the inequality m N + m π − m ∆ < m Ξ + mK − m Ω . Instead, larger difference of the potential term between the two systems compared with the kinetic term suggests that this inequality is mainly explained by a smaller spatial size of the wave function in ΞK due to a larger reduced mass, together with a strong attraction at short distances. This probably holds even at the physical pion mass, so that ∆ exists as a resonance while Ω is a stable particle.
Root-mean-square distances of bound states in both systems are very small, which indicates that in this setup, ∆ and Ω are tightly bound states and can be regarded qualitatively as 3-quark states.
Binding energies in both systems are consistent with those extracted from 2-point functions of single ∆ and Ω baryons, though systematic errors are rather large, mainly due to lattice artifacts at short distances. Large systematic errors are observed for the root-meansquare distances as well. This suggests that further improvement is desirable for the HAL QCD method to analyze states which behave like a single particle rather than a bound state of two hadrons, such as ∆ and Ω in this study. This is because a short distant part of the HAL QCD potentials between two hadrons, which is relevant for such a compact state, suffers severely from lattice artifacts, while a direct extraction of a single hadron mass from a 2-point correlation function of a single hadron operator is rather insensitive to such problems. In a setup where ∆ appears as a resonance, on the other hand, the HAL QCD method is expected to become more efficient.
In this work, we have performed the LO analysis in the derivative expansion of non-local potentials in the HAL QCD method. In the present study, it is implied that the LO analysis is sufficient since binding energies from potentials are consistent with those extracted from 2-point correlation functions within systematic uncertainties. However, when the systematic uncertainties from lattice discretization are improved, the next-leading order (NLO) analysis might be required as well, since the states are deeply bound below the threshold in this lattice setup. NLO analysis may be also necessary to study the ∆ as a resonance in future, since a resonance peak appears much higher energy than the threshold. The study in this direction is presented in Ref. [18].

VII. ACKNOWLEDGEMENTS
We thank the PACS-CS Collaboration [15] and ILDG/JLDG [27] for providing us their gauge configurations. We use lattice QCD code of Bridge++ [28,29] and its optimized version for the Oakforest-PACS by Dr. I. Kanamori [30]. potential. Due to a complicated structure with rapid changes of the 3-point function around the origin, the corresponding potential shows multi-valued behavior at short distances, which is, however, NOT caused by contaminations of higher partial waves due to a cubic box (long distance effects), as discussed below. Due to quark-antiquark annihilations possible in this channel, the operator-product expansion [32][33][34][35] predicts a behavior of the 3-point function at r → 0 as Y l,m (Ω) is the spherical harmonics with the angular momentum of the system, and thus F (r) has non-trivial angular dependence except for l = 0. We find that the behavior of Eq. (A1) for l = 1 at short distances looks very similar to a plot in Fig.11 (Left), and the corresponding potential also becomes multi-valued as in Fig.11 (Right), if the number of discrete data at short distances is too small to reproduce a smooth single-valued behavior.
Thus the problem is caused by short distance discretization effects. Note that a similar singular behavior was already observed in the I = 1 P-wave ππ system [18].
One of possible prescriptions to tame singular behaviors at short distances is the sink smearing at the quark level, introduced in [18]. Fig.12 (Left) represents the 3-point function with smeared sink quark operator, while Fig.12 (Right) is the corresponding potential. As expected, sink smearings make the behavior of the wave function at short distances much smoother, so that the potential becomes almost single-valued even at short distances where discrete data points are sparsely located.

Non-locality in potentials with smeared sink operators
Although smeared sink operators suppress singular behaviors in potentials at short distances, they may enhance non-locality in potentials. This may cause large systematic errors in the leading order analysis, as was observed in the case of the LapH method [36], In order to quantify effects of non-locality caused by our sink smearing, we compare scattering phase shifts between point and smeared sink operators in the same setup as the main text, taking I = 1 S-wave ΞK system. Fig. 13 show that phase shifts from different sink operators agree with each other at ∆E 100 MeV, where ∆E is an energy difference from the threshold. This suggests that non-locality due to the sink smearing has negligible effects on physical observables in the low energy region in this setup. where