Dibaryon with highest charm number near unitarity from lattice QCD

A pair of triply charmed baryons, $\Omega_{ccc}\Omega_{ccc}$, is studied as an ideal dibaryon system by (2+1)-flavor lattice QCD with nearly physical light-quark masses and the relativistic heavy quark action with the physical charm quark mass. The spatial baryon-baryon correlation is related to their scattering parameters on the basis of the HAL QCD method. The $\Omega_{ccc}\Omega_{ccc}$ in the ${^1S_0}$ channel taking into account the Coulomb repulsion with the charge form factor of $\Omega_{ccc}$ leads to the scattering length $a^{\rm C}_0\simeq -19~\text{fm}$ and the effective range $r^{\rm C}_{\mathrm{eff}}\simeq 0.45~\text{fm}$. The ratio $r^{\rm C}_{\mathrm{eff}}/a^{\rm C}_0 \simeq -0.024$, whose magnitude is considerably smaller than that of the dineutron ($-0.149$), indicates that $\Omega_{ccc}\Omega_{ccc}$ is located in the unitary regime.

A pair of triply charmed baryons, ΩcccΩccc, is studied as an ideal dibaryon system by (2+1)-flavor lattice QCD with nearly physical light-quark masses and the relativistic heavy quark action with the physical charm quark mass. The spatial baryon-baryon correlation is related to their scattering parameters on the basis of the HAL QCD method. The ΩcccΩccc in the 1 S0 channel taking into account the Coulomb repulsion with the charge form factor of Ωccc leads to the scattering length a C 0 −19 fm and the effective range r C eff 0.45 fm. The ratio r C eff /a C 0 −0.024, whose magnitude is considerably smaller than that of the dineutron (−0.149), indicates that ΩcccΩccc is located in the unitary regime.
Introduction.− Quantum chromodynamics (QCD) is a fundamental theory of strong interaction and governs not only the interaction among quarks and gluons but also the interaction between color-neutral hadrons. In particular, the nucleon-nucleon (N N ) interaction, which shows a characteristic mid-range attraction and a short-range repulsion, as well as the baryon-baryon (BB) interactions are important for describing the nuclear structure and dense matter relevant to nuclear physics and astrophysics [1][2][3][4][5].
As originally pointed out by Bjorken [15], the triply charmed baryon (the charm number C = 3) Ω ccc is stable against the strong interaction and provides an ideal ground to study the perturbative and non-perturbative aspects of QCD in the baryonic sector. Although it has not been observed yet experimentally 1 there have been numerous LQCD studies on its mass and electromagnetic form factor (see [18] and references therein). Accordingly, it is timely to study the Ω ccc Ω ccc as the simplest possible system to study heavy-baryon interactions. Its recent phenomenological study using the constituent quark model can be found in Ref. [19].
The purpose of this Letter is to study a system with the charm number C = 6 system, Ω ccc Ω ccc in the 1 S 0 channel, for the first time from first principle LQCD ap-proach. 2 The reason why we consider the S-wave and total spin s = 0 system is that the Pauli exclusion between charm quarks at short distance does not operate in this channel, so that the maximum attraction is expected in comparison to other channels. It is of critical importance to examine the scattering parameters such as the scattering length and the effective range to unravel the properties of such heavy dibaryons near threshold. The HAL QCD method [11,22,23], which treats the spatial correlation between two baryons on the lattice, provides a powerful tool for such analysis: Indeed, we show below that Ω ++ ccc Ω ++ ccc ( 1 S 0 ) with both strong interaction and Coulomb repulsion is located near unitarity [24,25] just above the threshold with a large negative scattering length.
HAL QCD Method.− The crucial steps in the HAL QCD method [11,22,23] are to obtain the equaltime Nambu-Bethe-Salpeter (NBS) wave function ψ(r) whose asymptotic behavior at a large distance reproduces the phase shifts, along with the corresponding twobaryon irreducible kernel U (r, r ). Since the same kernel U (r, r ) governs all the elastic scattering states, separating the ground state and the excited states on the lattice, which is exponentially difficult for baryon-baryon interactions [26,27], is not required to calculate the physical observables [23]. The normalized four-point function (the R-correlator) related to the NBS wave function is defined as where ∆W n = 2 m 2 Ωccc + k 2 n − 2m Ωccc with the baryon mass m Ωccc and the relative momentum k n . O(e −(∆E * )t ) denotes the contributions from the inelastic scattering states with ∆E * being the inelastic threshold, which are exponentially suppressed when t (∆E * ) −1 ∼ Λ −1 QCD with Λ QCD ∼ 300 MeV. J (0) is a source operator which creates two-baryon states with the charm number C = 6 at Euclidean time t = 0 and A n = n|J (0)|0 with |n representing the QCD eigenstates in a finite volume with ∆W n < ∆E * . In this study, we take a local interpolating operator, Ω ccc (x) ≡ lmn [c T l (x)Cγ k c m (x)]c n,α (x), where l, m, and n stand for color indices, γ k being the Dirac matrix, α being the spinor index, and C ≡ γ 4 γ 2 being the charge conjugation matrix.
When contributions from the inelastic scattering states are negligible (t (∆E * ) −1 ), the R-correlator satisfies [23] 1 4m Ωccc where H 0 = −∇ 2 /m Ωccc . By using the derivative expansion at low energies, U (r, r ) = V (r)δ(r − r ) + n=1 V 2n (r)∇ 2n δ(r − r ), the central potential V (r) in the leading order (LO) is given as The spatial and temporal derivatives of R(r, t) on the lattice are calculated in central difference scheme by using the nearest neighbor points. To extract the total spin s = 0, the following interpolating operators for the Ω ccc Ω ccc system is adopted, [Ω ccc Ω ccc ] 0 = ccc Ω 3/2 ccc ). Here the spin and its z component of the interpolating operator Ω sz ccc are 3/2 and s z = ±3/2, ±1/2, respectively, and Ω sz ccc is constructed by spin projection as shown in Ref. [28]. To obtain the orbital angular momentum L = 0 on the lattice, the projection to A 1 representation of the cubic group SO(3, Z) is employed; in Eq. (3) contains the channel coupling effect such as 1 S 0 -5 D 0 mixing and should be considered as an "effective" potential projected onto the S-wave state [29]. Lattice setup.− (2+1)-flavor gauge configurations are generated on the L 4 = 96 4 lattice with the Iwasaki gauge action at β = 1.82 and nonperturbatively O(a)-improved Wilson quark action combined with stout smearing at nearly physical quark masses (m π 146 MeV and m K 525 MeV) [30]. The lattice cutoff is a −1 2.333 GeV (a 0.0846 fm), corresponding to La 8.1 fm, which is sufficiently large to accommodate two heavy baryons. For the charm quark, we employ the relativistic heavy quark (RHQ) action in order to remove the leading order and the next-to-leading order cutoff errors associated with the charm quark mass [31]. We use two sets (set 1 and set 2) of RHQ parameters determined in Ref. [32] so as to interpolate the physical charm quark mass and reproduce the dispersion relation for the spin-averaged 1S charmonium, i.e. a weighted average of the spin-singlet state η c and the spin-triplet state J/Ψ.
For the source operator J (0), we use the wall type with the Coulomb gauge fixing. We employ the periodic (Drichlet) boundary condition for spatial (temporal) direction. We use 112 gauge configurations which are picked up one per ten trajectories. In order to reduce statistical fluctuations, forward and backward propagations are averaged, and four times measurements are performed by shifting source position along the temporal direction for each configuration. Then, the total measurements amount to 896 for each set. The statistical errors are estimated by the jackknife method with a bin size of 14 configurations. A comparison with a bin size of 7 configurations shows that the bin size dependence is small. The quark propagators are calculated by the Bridge++ code [33], and the unified contraction algorithm is utilized to obtain the correlation functions [34].  [35]. We have checked that our results for hadron masses are unchanged within errors by the fitting interval t/a = 30 − 35.
Numerical results.− The Ω ccc Ω ccc potential V (r) in the 1 S 0 channel from the interpolation between set 1 and set 2 is shown in Fig. 1 for t/a = 25, 26 and 27 (see Supplemental Material [36] for the t-dependence of V (r) in a wide range of t). Since the potentials from set 1 and set 2 are found to be consistent within statistical errors, the uncertainty in the interpolation is negligible. Λ −1 QCD ∼ 0.7 fm characterizing the inelastic states, and is small enough to avoid large statistical errors. We find that the potentials for t/a = 25, 26 and 27 are consistent with each other within statistical errors. This indicates that systematic errors due to inelastic states and higher order terms of the derivative expansion do not largely exceed the size of statistical errors [23] as we show below.
We find that the potential V (r) is repulsive at shortrange and attractive at mid-range, which has the same qualitative behaviors with the N N potential [37] and the ΩΩ potential [10]. The magnitude of the potential in the repulsive region r < 0.25 fm (corresponding to dV (r)/dr < 0) for Ω ccc Ω ccc is an order of magnitude smaller than that of ΩΩ obtained by the same method [10]. This may be qualitatively explained by the phenomenological quark model [38] as the color-magnetic interaction between constituent quarks is proportional to the square of reciprocal constituent quark mass. Qual- cm is the color-magnetic interaction between the quarks with flavor f and f with m * f being the constituent quark mass. On the other hand, the attraction in the region r > 0.25 fm (corresponding to dV (r)/dr > 0) may originate from the exchange of charmed mesons or rather be attributed to the direct exchange of charm quarks and/or multiple gluons. As can be seen in Fig. 1, the range of the potential is much smaller than the size of the lattice volume, indicating that the finite volume artifact is negligible.
In order to convert the potential to physical observables such as the scattering phase shifts and binding energy, we perform the uncorrelated fit for V (r) in Fig. 1 in the range r ≤ 2.5 fm by threerange Gaussians, V fit (r) = In Fig. 2, we show the Ω ccc Ω ccc scattering phase shifts δ 0 in the 1 S 0 channel calculated by solving the Schrödinger equation with the potential V (r) at t/a = 25, 26, and 27. The relativistic kinetic energy is defined as E CM = 2 k 2 + m 2 Ωccc − 2m Ωccc with a momentum k in the center of mass frame. The error bands reflect the statistical uncertainty of V (r). In all three cases, the phase shifts start from 180 • at E CM = 0, which indicates the existence of a bound state in Ω ccc Ω ccc system without Coulomb repulsion.
The central values and the statistical errors in the first parentheses are obtained at t/a = 26, while the systematic errors in the last parentheses are estimated from the values at t/a = 25, 26 and 27, which originates from the inelastic states and the higher order terms of the derivative expansion.
The binding energy B and the root-mean-square distance r 2 of the bound Ω ccc Ω ccc state are obtained These results are consistent with the general formula for loosely bound states [24,39] with scattering parameters a 0 and r eff : B = 1/(m Ωccc r 2 eff )(1 − 1 − (2r eff /a 0 )) 2 5.7 MeV and r 2 = a 0 / √ 2 1.1 fm. Since the binding energy and the size of the bound state from the strong interaction are not large, we need to take into account the Coulomb repulsion V Coulomb (r) between Ω ++ ccc s with finite spatial size. For this purpose, we consider the dipole form factor for Ω ++ ccc according to the LQCD study on the charge distribution of heavy baryons [18]: In the coordinate space, it corresponds to an exponential charge distribution ρ(r) = 12 6r/r d , where the charge radius r d = | r 2 charge | of Ω ++ ccc is taken to be r d = 0.410(6) fm [18]. Then, we have . The effective range expansion with Coulomb repulsion is written as where δ C 0 (k) is the phase shift in the presence of Coulomb repulsion, C 2 η = 2πη e 2πη −1 , η = 2α e m Ωccc /k, h(η) = Re[Ψ(iη)] − ln(η), and Ψ is the digamma function [40].
To see the effect of the Coulomb repulsion, we vary α e from zero to the physical value α phys. e = 1/137.036 below. Note that the systematic errors originated from the uncertainty in r d are found to be much smaller than the statistical errors and are neglected.
In Fig. 3, we show the inverse of scattering length 1/a C 0 under the change of α e /α phys. e from 0 to 1. Due to the large cancellation between the attractive strong interaction and the Coulomb repulsion, the result at α e /α phys. e = 1 is located very close to unitarity with a large scattering length The ratio r C eff /a C 0 = −0.024(0.010)( +0.006 −0.014 ) is considerably smaller in magnitude than that of the dineutron (−0.149).
Finally, we briefly discuss other possible systematic errors in this work: (i) The finite cutoff effect is O(α 2 s aΛ QCD , (aΛ QCD ) 2 ) thanks to the RHQ action for the charm quark and the non-perturbative O(a) improvement for light (u, d, s) quarks, and thus amounts to be O(1)%. (ii) In the vacuum polarization, light quark masses are slightly heavier than the physical ones and charm quark loop is neglected. The former effect is expected to be small since light quarks are rather irrelevant for Ω ccc Ω ccc system. In fact, the range of the Ω ccc Ω ccc potential is found to be shorter than 1 fm. The latter effect is suppressed due to the heavy charm quark mass, and is typically O(1)% [44]. These estimates for (i) and (ii) are also in line with the observation that our value of m Ωccc is consistent with that in the literature or has deviation of ∼ 1% at most, where we refer to LQCD studies by (2+1)-flavor at the physical point with finite a [35], (2+1)-flavor with chiral and continuum extrapolation [45] and (2+1+1)-flavor with chiral and continuum extrapolation [46,47]. In the future, these systematic The red up(down)-pointing triangle and the blue right(left)pointing triangle correspond to ΩcccΩccc system and ΩΩ system in the 1 S0 channel with(without) the Coulomb repulsion respectively. The black circle represents N N system in the 3 S1-3 D1 channel. The green square (nn) and diamond (pp) correspond to N N system in the 1 S0 channel. The error bars for ΩcccΩccc are the quadrature of the statistical and systematic errors in Eqs. (4) and (8).
errors will be evaluated explicitly. Moreover, a finite volume analysis with proper projection of the sink operator [27] for the Ω ccc Ω ccc system indicates that the truncation effect in the derivative expansion of U (r, r ) is small [36]. Further details will be reported elsewhere. Summary and discussions.− In this Letter, we presented a first investigation on the scattering properties of the Ω ccc Ω ccc on the basis of the (2+1)-flavor lattice QCD simulations with physical charm mass and nearly physical light quark masses. The potential for Ω ccc Ω ccc ( 1 S 0 ) obtained by the time-dependent HAL QCD method without the Coulomb interaction shows a weak repulsion at short distance surrounded by a relatively strong attractive well, which leads to a most charming (C = 6) dibaryon with the binding energy B 5.7 MeV and the size r 2 1.1 fm. By taking into account the Coulomb repulsion between Ω ++ ccc s with their charge form factor obtained from LQCD, the Ω ++ ccc Ω ++ ccc ( 1 S 0 ) system turns into the unitary region with r C eff /a C 0 −0.024. This provides good information toward the understanding of the interaction between heavy baryons. It is an interesting future work to study Ω − bbb Ω − bbb ( 1 S 0 ) for revealing the quark mass dependence of the scattering parameters. Finally, our results may further stimulate the future experimental activities to measure pair-momentum correlations of heavy baryons in high energy pp, pA and AA collisions [8,14].

Supplemental Material
We present more details about the systematics related to the inelastic excited states and the derivative expansion of the non-local potential in this supplemental material.
In Fig. 5, we plot the t-dependence of Ω ccc Ω ccc potential V (r) in the 1 S 0 channel at several distances, r = 0.08, 0.25, 0.46 and 1.00 fm. We find that V (r) at given r varies slowly with t. This indicates the contributions from the inelastic excited states are small irrespective of the value of r. Note that the corresponding t ∈ [2.0, 2.9] fm is sufficiently large compared to the scale relevant for the inelastic excited states, Λ −1 QCD ∼ 0.7 fm. The derivative expansion of the non-local potential up to the next-to-next-to-leading order (N 2 LO) was studied in Ref. [1]. It was found that the leading order (LO) result is accurate enough for heavy valence quarks at low energies. We expect the same to hold for the Ω ccc Ω ccc system in this study, and leave an explicit N 2 LO analysis along the line with the above paper for future study. In the main text, we have estimated the systematic error from N 2 LO by studying the t-dependence of the observables since the truncation effect manifests itself through such t-dependence as discussed in the above mentioned paper.
Yet another way to study the truncation effect in the derivative expansion was proposed in Ref. [2]. In the context of the present paper, it goes as follows. First, we construct a Hamiltonian H in a finite box with the LO poten-tial V (r) obtained when t/a = 26 as H = −∇ 2 /m Ωccc + V (r). Eigenvalues and eigenfunctions of the LO Hamiltonian H in the A 1 representation in a finite box are obtained from Hψ i = i ψ i . Then, an improved twobaryon sink operator for a designated eigenfunction can be constructed as r ψ † i (r)Ω ccc (r, t)Ω ccc (0, t). Equivalently, we can define the generalized temporal correlation function as R i (t) ≡ r ψ † i (r)R(r, t), from which the effective energy shift for the i-th eigenstate can be defined as ∆E eff i (t) = a −1 log(R i (t)/R i (t + 1)). In Table. II, we show ∆E eff i fitting from the interval t/a = 24 − 28 for the ground state (i = 0) and the first excited state (i = 1). Shown together are the energy shifts obtained by ∆E i = 2 i × m Ωccc + m 2 Ωccc − m Ωccc . We find that ∆E eff 0,1 from the generalized temporal correlation function R i (t) agree with the ∆E 0,1 obtained from the LO Hamiltonian H within the statistical errors. Although R i (t) utilizes the information of H through eigenfunctions, the agreement of the energy shifts indicates that higher order terms in the derivative expansion of the potential are not significant. In other words, if the effect of higher order terms were large, ψ i would be so different from the i-th eigenstate of the system that the effective energy shifts from R i (t) would be distorted and do not agree with those from H.