Possible bound states of Heavy Baryonium and Heavy Dibaryon systems
Jing-Juan Qi Zhen-Hua Zhang Xin-Heng Guo Zhen-Yang Wang

 
qijj@mail.bnu.edu.cn
zhangzh@usc.edu.cn
xhguo@bnu.edu.cn
wangzhenyang@nbu.edu.cn






§ ABSTRACT 
In this work, we systematically study the heavy baryonium and heavy dibaryon systems using the Bethe-Salpeter equation in the ladder and instantaneous approximations for the kernel. Our results indicate that all the heavy baryonium systems, specifically Λ_QΛ̅_Q, Ξ_QΞ̅_Q, Σ_QΣ̅_Q, Ξ'_QΞ̅'_Q, and Ω_QΩ̅_Q (Q=c, b), can form bound states. Among the heavy dibaryon systems, only the Ξ_QΞ_Q system with I=0 and the Σ_QΣ_Q systems with I=0 and I=1 can exist as bound states. Additionally, the Σ_QΣ̅_Q system with I=2 and the Σ_QΣ_Q system with I=1 are not deeply bound.

Possible bound states of Heavy Baryonium and Heavy Dibaryon systems
    Zhen-Yang Wang0000-0002-4074-7892
    Received 16 July 2024; accepted 04 September 2024
===================================================================








§ INTRODUCTION


Since the discovery of the X(3872) by the Belle Collaboration in 2003 <cit.>, numerous exotic states have been discovered in the charmed sector by various experiments, including BES3, BaBar, Belle, D0, ATLAS and LHCb, et al. (see, e.g., Refs. <cit.> for recent reviews). A common feature of these exotic states is that their masses are mostly located near the threshold of two hadrons. For example, X(3872) and Z_c(3900) are near the DD̅^∗ threshold, T^+_cc is near the DD^∗ threshold, Z_cs(3985) is near the D^∗D̅_s and DD̅^∗_s thresholds, P_c states are near the D̅^(∗)Σ_c thresholds, and X(6900) is near the χ_c0χ_c1 threshold. Therefore, these exotic hadrons are naturally considered as candidates for hadronic molecular states and believed to have four or five quarks. Their exotic spectra and decay widths have made them popular and intriguing topics in both theoretical and experimental research, deepen our understanding of the nature of QCD. 

Many heavy tetraquark and pentaquark states have already been discovered. Therefore, it is urgent to extend the research to heavy hexaquark states. The existence of the baryon-antibaryon (baryonium) and the baryon-baryon (dibaryon) molecular states has naturally become a significant research topic. In the light hexaquark sector, the deuteron is a well known molecular state composed of a proton and a neutron, with a binding energy of 2.225 MeV <cit.>. Recently, the BES3 experiment group reported the observation of a pp̅ bound state in the 3(π^+π^-) invariant mass spectrum <cit.>, which has been predicted by many theoretical works to favor decays into the final states with these pions <cit.>. In the charm sector, the Belle Collaboration observed Y(4630) (located 61 MeV above the Λ_cΛ̅_c threshold) in the e^+e^-→Λ_cΛ̅_c process in 2008 <cit.>, but no resonance structure was observed around 4.63 GeV by the BES3 Collaboration <cit.>. The nature of Y(4630) as a Λ_cΛ̅_c molecular state is highly debated in theory <cit.>. Compared with light baryon molecules, the larger masses of heavy baryons reduce the system's kinetic energy, facilitating the formation of molecules. Thus, the existence of heavy baryon molecules has attracted significant theoretical interest, and has been studied through various models such as the chiral constituent quark model <cit.>, the color flux-tube model <cit.>, the quark delocalization color screening model <cit.>, lattice QCD <cit.>, chiral effective field theory <cit.>, QCD sum rules <cit.>, the one-boson-exchange model <cit.>, and the quasipotential Bethe-Salpeter (BS) equation <cit.>.


In this work, we systematically investigate the existence of S-wave bound states composed of a heavy baryon and an antiheavy baryon or double heavy baryons in the BS equation approach within the ladder approximation and the instantaneous approximation for the kernel. Our model incorporates one free parameter, the cutoff parameter Λ, which is actually not entirely free as it governs the range of interaction and is directly related to the hadron size. Considering that the system may involve contributions from multiple exchange particles, with different interaction ranges, we reparametrize the cutoff parameter Λ as Λ = m + αΛ_QCD with m being the mass of the exchange particle. This approach allows for different cutoffs for various exchange particles through the varying parameter α which is of order unity.

This work is organized as follows. After the introduction,
we present the formalism in Sec. <ref>, which contains the
Lagrangians and the BS equations for the heavy baryonium and heavy dibaryon systems. In Sec. <ref>, we show the numerical results for the the heavy baryonium and heavy dibaryon systems. Finally, Sec. <ref> provides a brief summary and discussion. The isospin conventions and the wave functions for the charmed baryonium and charmed dibaryon systems are given in Appendix.




§ FORMALISM


To study whether the S-wave bound states of heavy baryonium and dibaryon exist, we fist construct the Lagrangians for heavy baryons and light mesons. Then the interaction kernels for the BS equations will be derived from the four-point Green’s function with the relevant Lagrangians.





 §.§ Effective chiral Lagrangians


A heavy baryon contains a heavy quark and two light quarks, which will be refered to as a diquark in the following. Each light quark is in a triplet representation of the flavor SU(3), thus the diquark can form either an antisymmetric antitriplet or a symmetric sextet. The diquark in the flavor-antisymmetric antitriplet has spin 0, and the diquark in the flavor-symmetric sextet has spin 1. Considering a ground state heavy baryon, the diquark combined with the heavy quark can form an antitriplet baryon with spin-1/2 (B^(Q)_3̅) and two sextet baryons with spin-1/2 (B^(Q)_6) and spin-3/2 (B^(Q)∗_6), respectively. The heavy baryon matrices are

    B^(c)_3̅=(
    [      0  Λ_c^+  Ξ_c^+; -Λ_c^+      0  Ξ_c^0; -Ξ_c^+ -Ξ_c^0      0;        ]),
    B^(c)_6=(
    [       Σ_c^++  1/√(2)Σ_c^+ 1/√(2)Ξ_c^'+;  1/√(2)Σ_c^+        Σ_c^0 1/√(2)Ξ_c^'0; 1/√(2)Ξ_c^'+ 1/√(2)Ξ_c^'0        Ω_c^0;              ]),


    B^(b)_3̅=(
    [      0  Λ_b^0  Ξ_b^0; -Λ_b^0      0  Ξ_b^-; -Ξ_b^0 -Ξ_b^-      0;        ]),
    B^(b)_6=(
    [        Σ_b^+  1/√(2)Σ_b^0 1/√(2)Ξ_b^'0;  1/√(2)Σ_b^0        Σ_b^- 1/√(2)Ξ_b^'-; 1/√(2)Ξ_b^'0 1/√(2)Ξ_b^'-        Ω_b^-;              ]),

and the matrices for B^(Q)∗_6 are similar to those for B^(Q)_6.

For convenience, while performing chiral-loop calculations, the two sextet heavy baryons can be combined to a superfield,

    S^μ=   B_6^∗μ-1/√(3)(γ^μ+v^μ)γ_5B_6,
    
        S̅^μ=   B̅_6^∗μ-1/√(3)(γ^μ+v^μ)γ_5B̅_6,

where v^μ is the velocity of the heavy baryon.

Then the general chiral Lagrangian for heavy baryons is <cit.>

    ℒ_B=ℒ_3̅+ℒ_S+ℒ_int,

with

    ℒ_3̅=   1/2tr[B̅_3̅(i v· D)B_3̅]+iβ_Btr[B̅_3̅v^μ(𝒱_μ-ρ_μ)B_3̅]+ℓ_Btr[B̅_3̅σ B_3̅],


    ℒ_S=   -tr[S̅_α(i v· D-Δ_B)S^α]+3/2g_1(iv_κ)ϵ^μνλκtr[S̅_μ𝒜_ν S_λ]
       +iβ_Str[S̅_μ v_α(𝒱^α-ρ^α)S^μ]+λ_Str[S̅_μ F^μνS_ν]+ℓ_Str[S̅_μσ S^μ],


    ℒ_int =   g_4tr[S̅_μ𝒜^μ B_3̅]+iλ_I ϵ^μνλκtr[S̅_ν F^μνS_ν]+H.c.,

where D_μ B=∂_μ B+𝒱_μ B+B𝒱_μ^T, Δ_B=M_6-M_3̅ is the mass difference between the sextet and the antitriplet, 𝒱_μ=1/2(ξ^†∂_μξ+ξ∂_μξ^†) and 𝒜_μ=1/2(ξ^†∂_μξ-ξ∂_μξ^†) are the vector and axial vector fields, respectively, F_μν=∂_μρ_ν-∂_νρ_μ+[ρ_μ,ρ_ν], ξ=exp[iP/f_π] and ρ=ig_V/√(2)V with 

    P=(
    [  π^0/√(2)+η/√(6)              π^+              K^+;              π^- -π^0/√(2)+η/√(6)              K^0;              K^-             K̅^0         -√(2/3)η;                  ]),

and

    V=(
    [ ω/√(2)+ρ^0/√(2)             ρ^+            K^∗+;             ρ^- ω/√(2)-ρ^0/√(2)            K^∗0;            K^∗-           K̅^∗+               ϕ;                 ]),

being the pseudoscalar and vector matrices, respectively. The phases of the fields B̅_3̅, B̅^(∗)_6, 𝒱_μ^T, and 𝒜_μ^T can be fixed by the following charge conjugation convention:

    B̅_3̅=𝒞B_3̅𝒞^-1,  B̅^(∗)_6=𝒞B^(∗)_6𝒞^-1,  𝒱_μ^T=-𝒞𝒱_μ𝒞^-1,  𝒜_μ^T=𝒞𝒜_μ𝒞^-1.


After expanding the effective Lagrangians in Eqs.(<ref>)-(<ref>) to the leading order of the light meson field, we can obtain the following effective interactions needed for our work:

    ℒ_B_3̅B_3̅V   =iβ_Bg_V/2√(2m_B̅_3̅m_B_3̅)tr[B̅_3̅∂_μ V^μ B_3̅],
    ℒ_B_3̅B_3̅σ   =ℓ_B tr[B̅_3̅σ B_3̅],
    ℒ_B_6B_6P   =-g_1/4f_π√(m_B̅_6m_B_6)ϵ_μνλκtr[B̅_6γ^μγ^λ∂^κ∂^ν P B_6],
    ℒ_B_6B_6V   =-iβ_Sg_V/2√(2m_B̅_6m_B_6)tr[B̅_6∂_ν V^ν B_6]-iλ_Sg_V/3√(2)tr[B̅_6γ^μ(∂_μ V_ν-∂_ν V_μ)γ^ν B_6],
    ℒ_B_6B_6σ   =-ℓ_Str[B̅_6σ B_6],

where v is replaced by i∂/(2√(m_B̅m_B)) and the pion decay constant is f_π=132 MeV. The values of relevant coupling constants are listed in Table <ref> <cit.>.









 §.§ The BS equation for the heavy baryonium system

 In this section, we will discuss the general BS formalism for the heavy baryonium composed of a heavy baryon and an anti-heavy baryon. In this case, the BS wave function is defined as: 

    χ_P(x_1,x_2,P)_αβ   =⟨0|Tψ_α(x_1)ψ̅_β(x_2)|P⟩,
       =e^-iP X∫d^4p/(2π)^4e^-ip xχ_P(p)_αβ,

where α and β are spinor indices, ψ(x_1) and ψ̅(x_2) are
the field operators of heavy baryon and anti-heavy baryon, respectively, P (=Mv) is the total momentum of the heavy baryonium and v represents its velocity, X=λ_1x_1-λ_2x_2 and x=x_1-x_2 are the center-of-mass coordinate and the relative coordinate of the heavy baryonium, respectively, with λ_1(2)=m_1(2)/m_1+m_2, where m_1 and m_2 are the masses of heavy baryon and anti-heavy baryon, respectively, p is the relative momentum of the heavy baryonium. The momenta of constituent particles can be expressed in terms of the relative momentum p and the total momentum P as p_1=λ_1P+p and p_2=λ_2P-p, respectively.

The BS equation for the heavy baryonium can be written as

    χ_P(p)=S(p_1)∫d^4q/(2π)^4K̅(P,p,q)χ_P(q)S(-p_2),

where K̅(P,p,q) is the interaction kernel, can derived from the irreducible Feynman diagrams, S(p_1) and S(-p_2) are the propagators of the heavy baryon and the anti-heavy baryon, respectively. For convenience, we define p_l≡ v· p as the longitudinal projection of p along v, and p_t≡ p-p_l v as the transverse component with respect to v.


In the leading of a 1/m_Q expansion, the propagators of the heavy baryon and the anti-heavy baryon can be expressed as:

    S(p_1) =im_1(1+v)/2w_1(λ_1M+p_l-w_1+iϵ),

and

    S(p_2)=im_2(1+v)/2w_2(λ_2M-p_l-w_2+iϵ),

where the energy w_1(2)=√(m^2_1(2)-p_t^2), and ϵ is an infinitesimal parameter.

Substituting Eqs. (<ref>) and (<ref>) into Eq. (<ref>), we obtain the following two constraint relations for the BS wave function χ_P(p):

    vχ_P(p)=χ_P(p),


    χ_P(p)v=-χ_P(p).


The S-wave heavy baryonium can have J^PC=0^-+ and J^PC=1^– states. With the constraints imposed by parity and Lorentz transformations, the BS wave functions can be expressed as the following: 

    χ_P(p)=γ_5f_1+γ_5vf_2+γ_5p_tf_3+(-i)σ_μνv^μ p_t^ν f_4,

and 

    χ_P^(r)(p)=   [p_t^ρ g_1+γ_μ(v^μ p_t^ρ g_2+p_t^μ p_t^ρ g_3+g^μρg_4)+γ_5γ_μϵ^μραβp_tαv_β g_5
       +σ_μν(p_t^μ v^ν p_t^ρ g_6+g^μρp_t^ν g_7+g^μρv^ν g_8)]ϵ_ρ^(r),

for the J^PC=0^-+ and J^PC=1^– S-wave heavy baryonia, respectively, where f_i (i=1,...,4) and g_j (j=1,...,8) are the Lorentz-scalar functions of p_t^2 and p_l, and ϵ^(r)_μ is the polarization vector of the vector heavy baryonium.


By applying the constraint relations (<ref>) and (<ref>), the BS wave functions for the S-wave pseudoscalar (J^PC=0^-+) and vector (J^PC=1^–) heavy baryonia can be simplified to the following forms, respectively:

    χ_P(p)=(1+v)γ_5f_1,

and

    χ^(r)_P(p)=(1+v)ϵ^(r)g_1.





 §.§ The BS equation for the heavy dibaryon

For the heavy dibaryon bound states composed of double heavy baryons, the general
form of the BS equationin in momentum space is: 

    χ_P(p)=S(p_1)∫d^4q/(2π)^4K̅(P,p,q)χ_P(q)S(p_2),

with the BS wave function for the heavy dibaryon being defined as:

    χ_P(x_1,x_2,P)_αβ=⟨0|Tψ_α(x_1)ψ_β(x_2)|P⟩.

For convenience, we define a deformed BS wave function,

    χ̃_P(p)_αβ=χ_P(p)_αγ𝒞_γβ^-1=(𝒞χ_P(p))_αβ,

where 𝒞 is the charge conjugation matrix. 

With this deformed BS wave function, the BS equation (<ref>) can be written in a more conventional matrix form

    χ̃_P(p)^T=S(p_1)∫d^4q/(2π)^4K̅(P,p,q)χ̃_P(q)^TS(-p_2),

where the superscript “T" represents the transpose of the spinor index.

From the BS equations (<ref>) and (<ref>), we see that the BS wave function χ_P(p) in Eq. (<ref>) for the heavy baryonium and the deformed BS wave function χ̃_P(p) in Eq. (<ref>) for the heavy dibaryon satisfy the same equation. And the deformed BS wave function χ̃_P(p) have the same forms as given in Eqs. (<ref>) and (<ref>). 

To simplify the BS equations (<ref>) and (<ref>), we impose the so-called covariant instantaneous approximation in the kernel: p_l=q_l. In this approximation, the projection of the momentum of each constituent particle along the total momentum P is not changed, i.e., the energy exchanged between the constituent particles of the binding system is neglected. This approximation is appropriate since we consider the binding energy of heavy baryonium and heavy dibaryon bound states to be very small compared to the masses of heavy baryons. Under this approximation, the kernel in the BS equation is reduced to K̅(P,p_t,q_t), which will be used in the following calculations.

After some algebra, we find that the BS scalar wave functions f_1 and g_1 satisfy the same integral equation as follows (in the following we will use f uniformly):

    f(p)=-m_1m_2/2w_1w_2(λ_1M+p_l-w_1+iϵ)(λ_2M-p_l-w_2+iϵ)∫d^4q/(2π)^4K̅(P,p_t,q_t)f(q).

We integrate both sides of the above equation with respect to p_l to obtain:

    f̃(p_t)=-im_1m_2/w_1w_2(M-w_1-w_2)∫d^3q_t/(2π)^3K̅(P,p_t,q_t)f̃(q_t),

where we have defined f̃(p_t)=∫ dp_l f(p).



Based on the effective Lagrangians in Eq.(<ref>), the lowest-order interaction kernel can be derived as follows:

    K̅_B_3̅B̅_3̅^V(P,p,q)=   c_I(g_Vβ_B/2√(2m_B_3̅m_B̅_3̅))^2(p_1+q_1)_μ(p_2+q_2)_νΔ^μν_V(k),
    K̅_B_3̅B̅_3̅^σ(P,p,q)=   -c_Iℓ_B^2Δ_σ(k),
    K̅_B_6B̅_6^P(P,p,q)=   -c_I(g_1/4f_π√(m_B_6m_B̅_6))^2ϵ_μνλκϵ_αβτδγ^μγ^λγ^αγ^τ(p_1+q_1)^κ(p_2+q_2)^δ k^ν k^βΔ_P(k),
    K̅_B_6B̅_6^V(P,p,q)=   c_I[β_Sg_V/2√(2m_B_6m_B̅_6)(p_1+q_1)^μ g_μτ-λ_Sg_V/3√(2)γ^μγ^ν(k_μ g_ντ-k_ν g_μτ)]
       ×[β_Sg_V/2√(2m_B_6m_B̅_6)(p_2+q_2)^α g_ακ-λ_Sg_V/3√(2)γ^αγ^β(k_α g_βκ-k_β g_ακ)]Δ_V^τκ(k),
    K̅_B_6B̅_6^σ(P,p,q)=   -c_Iℓ_S^2Δ_σ(k),

where Δ^μν_V(k), Δ_P(k), and Δ^σ(k) are the propagators of the exchanged vector, pseudoscalar and σ mesons, respectively, k represents the momentum of the exchanged meson, and c_I is the isospin coefficient, given in Table <ref>. In our model, the BS wave function depends only on the isospin I but not on its
component I_3 because we consider only strong interactions that preserve the isospin symmetry.





To account for the structure and finite size effects of the interacting hadrons, it is necessary to introduce the form factor at the vertices. For t-channel vertices, we use the monopole form factor:

    F_M(k^2)=   Λ^2-m^2/Λ^2-k^2,

where m and Λ represent the mass and cutoff parameter of the exchanged meson, respectively. Since the heavy baryonium and heavy dibaryon systems can have interaction by exchanging multiple particles, different masses of exchanged particles correspond to different interaction ranges and, consequently, different cutoff parameters. Thus, we further reparameterize the cutoff Λ as Λ=m+αΛ_QCD with Λ_QCD = 220 MeV, where the parameter α is of order one. The value of α depends on the exchanged and external particles involved in the strong interaction vertex and cannot be obtained from the first principle.



§ NUMERICAL RESULTS


In the numerical calculations, we first present the masses of the relevant mesons and heavy baryons in Table <ref> <cit.>, which are essential for investigating whether heavy baryonium and heavy dibaryon systems can exist as bound states. In our model, we have two parameters, the cutoff Λ and the bounding energy E_b. The cutoff Λ is reparameterized as a variable, α, with Λ=m+αΛ_QCD. The parameter α is not a completely free parameter, since its varying range is related to the sizes of hadrons <cit.>. Based on the experience with deuteron, the parameter α is typically of order unity. The other parameter E_b (defined as E_b= m_1+m_2-M, where we consider the heavy baryonium and dibaryon systems as shallow bound states with E_b ranging from 0 to 50 MeV), is dependent on the value of the parameter α, and is therefore not absolutely determined. In this work, we allow the parameter α to vary over a wide range (0.3–8) to search for possible solutions in the heavy baryonium and dibaryon systems.

To solve the three-dimensional integral BS equation (<ref>), we fist simplify it to a one-dimensional integral equation by completing the azimuthal integration. This one-dimensional integral BS equation is further discretised into a matrix eigenvalue equation by the Gaussian quadrature method. By solving the eigenvalue equation, we can find the possible bound states of the heavy baryonium and heavy dibaryon systems depending on the parameter α.













 §.§ The results of charmed baryonium and charmed dibaryon systems

The results for the possible bound states of the charmed baryonium and charmed dibaryon systems are shown in Figs. <ref> and <ref>, respectively. Our research indicates that all the charmed baryonium systems, specifically Λ_cΛ̅_c, Ξ_cΞ̅_c, Σ_cΣ̅_c, Ξ'_cΞ̅'_c, and Ω_cΩ̅_c, can exist as bound states. Among the charmed dibaryon systems, only the Ξ_cΞ_c system with isospin I=0 and the Σ_cΣ_c system with isospin I=0 and I=1 can exist as bound states. 

For the Λ_cΛ̅_c (Λ_c) system, since Λ_c is an isoscalar state, the interaction kernel arises from the exchanges of ω and σ mesons. Both ω and σ mesons induce attractive interaction in the Λ_cΛ̅_c system, allowing it to form a bound state in our model. The values of the parameter α along with the corresponding binding energy E_b are displayed in Fig. <ref>(a). It is also found that this system can exist as a bound state in various models <cit.>. However, the results of Refs. <cit.> show that the binding energy is sensitive to the cutoff Λ. 

For the Λ_cΛ_c system, the interaction contributed from ω is repulsive and that from σ is attractive. Our result indicates that the Λ_cΛ_c system cannot form a bound state which is consistent with the results in Refs. <cit.> that the Λ_cΛ_c system can not be a bound state by itself. However, the coupling of the Λ_cΛ_c to the strongly attractive Σ^(∗)_cΣ^(∗)_c system may lead to a state below the Λ_cΛ_c threshold <cit.>. On the contrary, in Refs. <cit.>, it is pointed out that the single channel Λ_cΛ_c can form a bound state. This discrepancy arises from the fact that in these models, the attractive contribution from the σ meson is stronger than the repulsive contribution from the ω meosn. However, in our model, we find that even only considering the contribution from the σ meson is not sufficient for the Λ_cΛ_c system to form a bound state. 



The Ξ_c baryon contains a strange quark and has an isospin of 1/2. Therefore, the Ξ_cΞ̅_c (Ξ_c) system can have isospins of both I=0 and I=1 and the interaction kernel can arise from the exchanges of ρ, ω, ϕ, and σ. For the Ξ_cΞ̅_c system with I=0, the exchanges of ρ, ω, ϕ, and σ mesons all induce attractive interaction. The Ξ_cΞ̅_c system with I=0 can form a bound state with a binding energy in the range from 0 to 50 MeV when the parameter α ranges from 1.12 to 2.88, which is presented in Fig. <ref>(b). In the Ξ_cΞ̅_c system with I=1, the interaction magnitudes due to ρ and ω are the same, but ρ provides a repulsive contribution, thus their contributions almost cancel each other considering the similar masses of ρ and ω. Then the Ξ_cΞ̅_c system with I=1 is able to exist as a bound state with a binding energy in the range of 0 to 50 MeV when the parameter α ranges from  1.32 to 3.97. The relevant results are presented in Fig. <ref>(c). That the system Ξ_cΞ̅_c with I=0 and I=1 can exist as a bound state is also supported by Refs. <cit.>. It is worth mentioning that in the lattice QCD <cit.> and the chromomagnetic interaction model <cit.> it is found that the masses of the hidden-charm and hidden-strange hexaquarks are below the Ξ_cΞ̅_c threshold by 700-1000 MeV, which cannot be obtained within a reasonable range of the parameter α in our model.







For the Ξ_cΞ_c system, only the isospin I=0 configuration can exist as a bound state, as depicted in Fig. <ref>(a). In the Ξ_cΞ_c system with I=1, contributions from the vector mesons ρ, ω, and ϕ are repulsive, while that from the  σ meson is attractive but insufficient to form a bound state in our model. However, the Ξ_cΞ_c system with I=1 could be a loosely bound state with a binding energy of only a few hundred keV within the one-boson-exchange model <cit.>, and it could be a deeply bound state within the quasipotential BS equation framework <cit.>. In our model, unlike the Λ_cΛ_c system, the Ξ_cΞ_c system can form a bound state when only considering the σ meson exchange due to the greater mass of the Ξ_c compared with the Λ_c. However, given that the total contributions of the vector mesons in the I=1 Ξ_cΞ_c system are repulsive, it is quite inconceivable that this system could form a bound state.







For the Σ_cΣ̅_c(Σ_c) system, the interaction kernels are induced by the exchange of the pseudoscalar mesons π and η, the vector mesons ρ and ω, and the scalar meson σ. In the Σ_cΣ̅_c system, our results indicate that the isospin states with I = 0, 1, and 2 can all exist as bound states, consistent with Refs. <cit.>. However, in our model, the Σ_cΣ̅_c system with I = 2 is a loosely bound state with the binding energy very sensitive to the parameter α compared with the isospin states I = 0 and 1, as shown in Figs. <ref>(d)-<ref>(f). This sensitivity is due to the contributions from ρ and ω mesons nearly canceling out, while the attraction from the π meson accounts for the long-range interaction. In contrast, Ref. <cit.> reports a binding energy of 149.66 MeV for this state with a cutoff Λ = 0.95.

In the Σ_cΣ_c system, the isospin states with I=0 and I=1 can form bound states, as presented in Figs. <ref>(b) and <ref>(c). Specifically, the Σ_cΣ_c system with I=1 is a loosely bound state, consistent with the one-boson-exchange model <cit.> and the chiral effective theory <cit.>. The binding energy E_b of the Σ_cΣ_c system with I=1 is very sensitive to the parameter α. In our model, contributions from pseudoscalar mesons π and η, and vector mesons ρ and ω are repulsive. Only the σ meson provides an attractive force, which is insufficient to form a bound state for the I=2 Σ_cΣ_c system. In contrast, Ref. <cit.> suggests that the Σ_cΣ_c system with I=2 can form a bound state.


For the Ξ'_cΞ̅'_c(Ξ'_c) system, only the charmed baryonium Ξ'_cΞ̅'_c system can exist as bound state, with results presented in Figs. <ref>(g) and <ref>(h). The charmed dibaryon Ξ'_cΞ'_c system cannot exist as a bound state in our model, which is inconsistent with Refs. <cit.> and <cit.>. According to Ref. <cit.>, as the root mean square (rms) radius increases, the vector mesons that originally provided repulsive contributions become attractive, allowing the Ξ'_cΞ'_c system to exist as a loosely bound state. For the Ξ'_cΞ'_c system with I=0, contributions from η, ρ, and σ are attractive, while contributions from π, ω, and ϕ are repulsive, but no bound state is found. Thus, except for the scalar σ meson, which provides an attractive contribution, all other particles contribute repulsive forces to the Ξ'_cΞ'_c system with I=1, so that it cannot exist as a bound state in our model.

For the Ω_cΩ̅_c(Ω_c) system, only the Ω_cΩ̅_c can exist as a bound state. In the Ω_cΩ_c system, the ρ and η provide repulsive contributions, and the attractive contribution from σ exchange alone is insufficient to form a bound state. However, Ref. <cit.> suggests that the Ω_cΩ_c system can exist as a loosely bound state, because the repulsive contributions of η and ϕ decrease rapidly as the rms radius increases, making the attraction contribution provided by σ greater than the repulsion.





 §.§ The results of bottom baryonium and bottom dibaryon systems



The interaction kernels in the bottom sector are the same as those in the charm sector. Thus, similar to the charmed baryonium and charmed dibaryon systems, all the bottom baryonium systems, such as Λ_bΛ̅_b, Ξ_bΞ̅_b, Σ_bΣ̅_b, and Ω_bΩ̅_b, as well as the bottom dibaryon systems Ξ_bΞ_b with isospin I=0 and Σ_bΣ_b with isospin I=0 and I=1, can exist as bound states. The results for the parameter α and the corresponding binding energy E_b are displayed in Figs. <ref> and <ref>. Because of the much heavier reduced masses of hidden-bottom systems, it is easier to form bound states than in charmed systems. Therefore, for the same binding energy, the bottom region corresponds to a smaller parameter α.

Similar to the Σ_cΣ̅_c system with I=2 and the Σ_cΣ_c system with I=1 in the charm region, the I=2 Σ_bΣ̅_b and I=1 Σ_bΣ_b systems show the binding energies being very sensitive to the parameter α. These two systems are also unable to bind very deeply, with the corresponding maximum binding energies of 70 MeV (α = 5.32) and 62 MeV (α = 4.28), respectively. To reasonably apply the instantaneous approximation in solving the BS equation (<ref>), we choose a maximum binding energy of 50 MeV. Therefore, the results for binding energies larger than 50 MeV are not shown in Figs. <ref>, <ref>, <ref>, and <ref>.






















§ SUMMARY AND DISCUSSION



In this work, we utilized the BS equation to systematically study whether heavy baryonium and heavy dibaryon systems can exist as bound states. Our research indicates that all the heavy baryonium systems, including Λ_QΛ̅_Q, Ξ_QΞ̅_Q, Σ_QΣ̅_Q, Ξ'_QΞ̅'_Q, and Ω_QΩ̅_Q (Q=c, b), can exist as bound states. Among the heavy dibaryon systems, only the Ξ_QΞ_Q system with isospin I=0 and the Σ_QΣ_Q systems with isospin I=0 and I=1 can exist as bound states. Additionally, we found that the Σ_QΣ̅_Q system with I=2 and the Σ_QΣ_Q system with I=1 cannot exist as very deeply bound states. Furthermore, the large mass of heavy baryons reduces the kinetic energy of the system, making it easier to form bound states. Therefore, as shown in Figs. <ref>, <ref>, <ref>, and <ref>, the parameter α required to form bound states in the bottom region is smaller than that in the charm region, implying that the binding in the bottom region is deeper than that in the charm region.

However, there is considerable debate among different models regarding whether heavy baryonium and heavy dibaryon systems can exist as bound states, especially for the heavy dibaryon systems. In our model, the contribution of the ω meson in the Λ_cΛ_c system is repulsive, and the attractive contribution of the σ meson is insufficient to form a bound state in the Λ_cΛ_c system. Nevertheless, in many other models <cit.>, the Λ_cΛ_c system can exist as a bound state. In Refs. <cit.>, the Λ_cΛ_c system cannot be a bound state by itself but it is shown that the coupling to the strongly attractive Σ_c^(∗)Σ_c^(∗) system may lead to a state below the Λ_cΛ_c threshold. In the one-boson-exchange model <cit.>, the Ξ'_cΞ'_c system and the Ω_cΩ_c system can exist as shallow bound states, while the I=2 Σ_cΣ̅_c system can exist as a deeply bound state. Therefore, the existence of these bound states requires further theoretical studies and experimental verification.

The charmed baryonium bound states can be studied via B decays and e^+e^- collisions at LHCb, RHIC, Belle II, and BES3. With the upcoming BEPCII upgrade to 5.6 GeV by the end of 2024, as well as the completion of PANDA and the Super Tau-Charm Factory, the detailed study of charmed baryon-antibaryon bound states will become possible. Compared with the production of charmed baryon-antibaryon bound states, the production of charmed dibaryon bound states is significantly more challenging and faces immense difficulties, although it can still occur at LHC and RHIC. Charmed dibaryon bound states are highly stable because their constituent particles primarily decay through weak interactions, leading to long lifetimes (except for Σ_c which primarily decays via the strong process Σ_c →Λ_c π). However, due to the larger masses, bound states in the bottom region are more difficult to be produced than those in the charm region.




This work was supported by National Natural Science Foundation of China (Project Nos. 12105149, 12405115 12475096 and 12275024). 










99


                               


Belle:2003nnu
S. K. Choi et al. [Belle],

Phys. Rev. Lett. 91, 262001 (2003).







Hosaka:2016pey
A. Hosaka, T. Iijima, K. Miyabayashi, Y. Sakai and S. Yasui,

PTEP 2016, 062C01 (2016).





Meng:2022ozq
L. Meng, B. Wang, G. J. Wang and S. L. Zhu,

Phys. Rept. 1019, 1-149 (2023)





Chen:2022asf
H. X. Chen, W. Chen, X. Liu, Y. R. Liu and S. L. Zhu,

Rept. Prog. Phys. 86, 026201 (2023).





Brambilla:2019esw
N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan,

Phys. Rept. 873, 1-154 (2020).





Liu:2019zoy
Y. R. Liu, H. X. Chen, W. Chen, X. Liu and S. L. Zhu,

Prog. Part. Nucl. Phys. 107, 237-320 (2019).





Ali:2017jda
A. Ali, J. S. Lange and S. Stone,

Prog. Part. Nucl. Phys. 97, 123-198 (2017).





Guo:2017jvc
F. K. Guo, C. Hanhart, U. G. Meißner, Q. Wang, Q. Zhao and B. S. Zou,

Rev. Mod. Phys. 90, 015004 (2018),
[erratum: Rev. Mod. Phys. 94, 029901 (2022)].





Olsen:2017bmm
S. L. Olsen, T. Skwarnicki and D. Zieminska,

Rev. Mod. Phys. 90, 015003 (2018).





Lebed:2016hpi
R. F. Lebed, R. E. Mitchell and E. S. Swanson,

Prog. Part. Nucl. Phys. 93, 143-194 (2017).





Esposito:2016noz
A. Esposito, A. Pilloni and A. D. Polosa,

Phys. Rept. 668, 1-97 (2017).






Weinberg:1962hj
S. Weinberg,

Phys. Rev. 130, 776-783 (1963).




Weinberg:1963zza
S. Weinberg,

Phys. Rev. 131, 440-460 (1963).




Weinberg:1965zz
S. Weinberg,

Phys. Rev. 137, B672-B678 (1965).




BESIII:2023vvr
M. Ablikim et al. [BESIII],

Phys. Rev. Lett. 132, 151901 (2024).





Yan:2004xs
M. L. Yan, S. Li, B. Wu and B. Q. Ma,

Phys. Rev. D 72, 034027 (2005).





Ding:2005ew
G. J. Ding and M. L. Yan,

Phys. Rev. C 72, 015208 (2005).





Yang:2022kpm
Q. H. Yang, D. Guo and L. Y. Dai,

Phys. Rev. D 107, 034030 (2023).







Belle:2008xmh
G. Pakhlova et al. [Belle],

Phys. Rev. Lett. 101, 172001 (2008).





BESIII:2023rwv
M. Ablikim et al. [BESIII],

Phys. Rev. Lett. 131, 191901 (2023).





Cotugno:2009ys
G. Cotugno, R. Faccini, A. D. Polosa and C. Sabelli,

Phys. Rev. Lett. 104, 132005 (2010).





Guo:2010tk
F. K. Guo, J. Haidenbauer, C. Hanhart and U. G. Meissner,

Phys. Rev. D 82, 094008 (2010).





Simonov:2011jc
Y. A. Simonov,

Phys. Rev. D 85, 105025 (2012)





Liu:2021gva
Z. Liu, H. T. An, Z. W. Liu and X. Liu,

Phys. Rev. D 105, 034006 (2022).





Agaev:2022iha
S. S. Agaev, K. Azizi and H. Sundu,

Phys. Rev. D 106, 014025 (2022).





Song:2022yfr
L. Q. Song, D. Song, J. T. Zhu and J. He,

Phys. Lett. B 835, 137586 (2022).





Mei:2022msh
X. H. Mei, Z. Yu, M. Song, J. Y. Guo, G. Li and X. Luo,

Chin. Phys. C 47, 033104 (2023).
doi:10.1088/1674-1137/aca959
[arXiv:2212.02218 [hep-ph]].



Salnikov:2023cug
S. G. Salnikov, A. E. Bondar and A. I. Milstein,

Nucl. Phys. A 1041, 122764 (2024).





Wang:2013exa
Z. G. Wang,

Eur. Phys. J. C 74, 2874 (2014).





Liu:2016sip
X. Liu, H. W. Ke, X. Liu and X. Q. Li,

Eur. Phys. J. C 76, 549 (2016).





Guo:2016iej
X. D. Guo, D. Y. Chen, H. W. Ke, X. Liu and X. Q. Li,

Phys. Rev. D 93, 054009 (2016).





Wang:2016fhj
Y. Y. Wang, Q. F. Lü, E. Wang and D. M. li,

Phys. Rev. D 94, 014025 (2016).





Dai:2017fwx
L. Y. Dai, J. Haidenbauer and U. G. Meißner,

Phys. Rev. D 96, 116001 (2017).





Sundu:2018toi
H. Sundu, S. S. Agaev and K. Azizi,

Phys. Rev. D 98, 054021 (2018).





Anwar:2018sol
M. N. Anwar, J. Ferretti and E. Santopinto,

Phys. Rev. D 98, 094015 (2018).





Cao:2019wwt
Q. F. Cao, H. R. Qi, Y. F. Wang and H. Q. Zheng,

Phys. Rev. D 100, 054040 (2019).







Carames:2015sya
T. F. Carames and A. Valcarce,

Phys. Rev. D 92, 034015 (2015).





Garcilazo:2020acl
H. Garcilazo and A. Valcarce,

Eur. Phys. J. C 80, 720 (2020).





Deng:2013aca
C. Deng, J. Ping, Y. Yang and F. Wang,

Phys. Rev. D 88, 074007 (2013).





Huang:2013rla
H. Huang, J. Ping and F. Wang,

Phys. Rev. C 89, 035201 (2014).






Junnarkar:2022yak
P. M. Junnarkar and N. Mathur,

Phys. Rev. D 106, 054511 (2022).





Chen:2022iil
K. Chen, B. L. Huang, B. Wang and S. L. Zhu,

[arXiv:2204.13316 [hep-ph]].



Lu:2017dvm
J. X. Lu, L. S. Geng and M. P. Valderrama,

Phys. Rev. D 99, 074026 (2019).





Oka:2013xxa
M. Oka,

Nucl. Phys. A 914, 447-453 (2013).




Chen:2013sba
Y. D. Chen, C. F. Qiao, P. N. Shen and Z. Q. Zeng,

Phys. Rev. D 88, 114007 (2013).





Chen:2011cta
Y. D. Chen and C. F. Qiao,

Phys. Rev. D 85, 034034 (2012).





Wang:2021pua
X. W. Wang and Z. G. Wang,

Adv. High Energy Phys. 2022, 6224597 (2022).





Wang:2021qmn
X. W. Wang, Z. G. Wang and G. l. Yu,

Eur. Phys. J. A 57, 275 (2021).





Wan:2019ake
B. D. Wan, L. Tang and C. F. Qiao,

Eur. Phys. J. C 80, 121 (2020).





Cheng:2022vgy
J. B. Cheng, D. x. Zheng, Z. Y. Lin and S. L. Zhu,

Phys. Rev. D 107, 054018 (2023).





Ling:2021asz
X. Z. Ling, M. Z. Liu and L. S. Geng,

Eur. Phys. J. C 81, 1090 (2021).





Chen:2017vai
R. Chen, A. Hosaka and X. Liu,

Phys. Rev. D 96, 116012 (2017).





Meguro:2011nr
W. Meguro, Y. R. Liu and M. Oka,

Phys. Lett. B 704, 547-550 (2011).





Lee:2011rka
N. Lee, Z. G. Luo, X. L. Chen and S. L. Zhu,

Phys. Rev. D 84, 014031 (2011).





Li:2012bt
N. Li and S. L. Zhu,

Phys. Rev. D 86, 014020 (2012).





Song:2022svi
D. Song, L. Q. Song, S. Y. Kong and J. He,

Phys. Rev. D 106, 074030 (2022).





Song:2023vtu
D. Song, S. Chen, S. Y. Kong and J. He,

Chin. Phys. C 47, 113102 (2023).




Kong:2023dwz
S. Y. Kong, J. T. Zhu and J. He,

Eur. Phys. J. C 83, 436 (2023).






Liu:2011xc
Y. R. Liu and M. Oka,

Phys. Rev. D 85, 014015 (2012).





Cheng:1993kp
H. Y. Cheng, C. Y. Cheung, G. L. Lin, Y. C. Lin, T. M. Yan and H. L. Yu,

Phys. Rev. D 49 (1994), 5857-5881
[erratum: Phys. Rev. D 55 (1997), 5851-5852].





ParticleDataGroup:2024prd
S. Navas et al. [Particle Data Group],

Phys. Rev. D 110, 030001 (2024).







Dong:2021juy
X. K. Dong, F. K. Guo and B. S. Zou,

Progr. Phys. 41, 65-93 (2021).






Gerasyuta:2011zx
S. M. Gerasyuta and E. E. Matskevich,

Int. J. Mod. Phys. E 21, 1250058 (2012).





Yu:2021lmb
Z. Yu, M. Song, J. Y. Guo, Y. Zhang and G. Li,

Phys. Rev. C 104 (2021), 035201.





Dong:2021bvy
X. K. Dong, F. K. Guo and B. S. Zou,

Commun. Theor. Phys. 73 (2021), 125201.





Liu:2022gxf
H. Liu, J. He, L. Liu, P. Sun, W. Wang, Y. B. Yang and Q. A. Zhang,

Sci. China Phys. Mech. Astron. 67 (2024), 211011.












§ THE FLAVOUR WAVE FUNCTIONS


For the isospin conventions, we use the following ones:

    | u⟩=|1/2,1/2⟩,  | d⟩=|1/2,-1/2⟩,  |u̅⟩=|1/2,-1/2⟩,  |d̅⟩=-|1/2,1/2⟩.

Then, we have 

    |Λ_c^+⟩   =|0,0⟩,  |Λ_c^-⟩=-|0,0⟩,  |Σ_c^++⟩=|1,1⟩,  |Σ_c^–⟩=|1,-1⟩,
    |Σ_c^+⟩   =|1,0⟩,  |Σ_c^-⟩=-|1,0⟩,  |Σ_c^0⟩=|1,-1⟩,  |Σ̅_c^0⟩=|1,1⟩,
    |Ξ_c^(')+⟩   =|1/2,1/2⟩,  |Ξ_c^(')-⟩=|1/2,-1/2⟩,  |Ξ_c^(')0⟩=|1/2,-1/2⟩,  |Ξ̅_c^(')0⟩=-|1/2,1/2⟩,
    |Ω_c^0⟩   =|0,0⟩,  |Ω̅_c^0⟩=|0,0⟩.

The flavour wave functions of the charmed baryon and anti-charmed baryon (charmed baryon) systems can be construct with the Clebsch-Gordan coefficients and above conventions,

    Ψ_Λ_cΛ̅_c^0,0   =-|Λ_cΛ̅_c⟩,
    Ψ_Σ_cΣ̅_c^0,0   =1/√(3)|Σ_c^++Σ_c^–+Σ_c^+Σ_c^-+Σ_c^0Σ̅_c^0⟩,
    Ψ_Σ_cΣ̅_c^1,0   =1/√(2)|Σ_c^++Σ_c^–-Σ_c^0Σ̅_c^0⟩,Ψ_Σ_cΣ̅_c^1,1=-1/√(2)|Σ_c^++Σ_c^-+Σ_c^+Σ̅_c^0⟩,Ψ_Σ_cΣ̅_c^1,-1=1/√(2)|Σ_c^+Σ_c^–+Σ_c^0Σ_c^-⟩,
    Ψ_Σ_cΣ̅_c^2,0   =1/√(6)|Σ_c^++Σ_c^–-2Σ_c^+Σ_c^-+Σ_c^0Σ̅_c^0⟩,Ψ_Σ_cΣ̅_c^2,1=-1/√(2)|Σ_c^++Σ_c^--Σ_c^+Σ̅_c^0⟩,
    Ψ_Σ_cΣ̅_c^2,2   =|Σ_c^++Σ̅_c^0⟩, Ψ_Σ_cΣ̅_c^2,-1=1/√(2)|Σ_c^+Σ_c^–-Σ_c^0Σ_c^-⟩,Ψ_Σ_cΣ̅_c^2,-2=|Σ_c^0Σ_c^–⟩,
    Ψ_Ξ_c^(')Ξ̅_c^(')^0,0   =1/√(2)|Ξ_c^(')+Ξ̅_c^(')-+Ξ_c^(')0Ξ̅_c^(')0 ⟩,
    Ψ_Ξ_c^(')Ξ̅_c^(')^1,0   =1/√(2)|Ξ_c^(')+Ξ_c^(')--Ξ_c^(')0Ξ̅_c^(')0⟩, Ψ_Ξ_c^(')Ξ̅_c^(')^1,1=-|Ξ_c^(')+Ξ̅_c^(')0⟩, Ψ_Ξ_c^(')Ξ̅_c^(')^1,-1=|Ξ_c^(')0Ξ_c^(')-⟩,
    Ψ_Ω_cΩ̅_c^0,0   =|Ω_c^0Ω̅_c^0⟩,


and

    Ψ_Λ_cΛ_c^0,0   =|Λ_cΛ_c⟩,
    Ψ_Σ_cΣ_c^0,0   =1/√(3)|Σ_c^++Σ_c^0-Σ_c^+Σ_c^++Σ_c^0Σ_c^++⟩,
    Ψ_Σ_cΣ_c^1,0   =1/√(2)|Σ_c^++Σ_c^0-Σ_c^0Σ_c^++⟩,Ψ_Σ_cΣ_c^1,1=1/√(2)|Σ_c^++Σ_c^+-Σ_c^+Σ_c^++⟩,Ψ_Σ_cΣ_c^1,-1=1/√(2)|Σ_c^+Σ_c^0-Σ_c^0Σ_c^+⟩,
    Ψ_Σ_cΣ_c^2,0   =1/√(6)|Σ_c^++Σ_c^0+2Σ_c^+Σ_c^++Σ_c^0Σ_c^++⟩,Ψ_Σ_cΣ_c^2,1=1/√(2)|Σ_c^++Σ_c^++Σ_c^+Σ_c^++⟩,
    Ψ_Σ_cΣ_c^2,2   =|Σ_c^++Σ_c^++⟩, Ψ_Σ_cΣ_c^2,-1=1/√(2)|Σ_c^+Σ_c^0+Σ_c^0Σ_c^+⟩,Ψ_Σ_cΣ_c^2,-2=|Σ_c^0Σ_c^0⟩,
    Ψ_Ξ_c^(')Ξ_c^(')^0,0   =1/√(2)|Ξ_c^(')+Ξ_c^(')0-Ξ_c^(')0Ξ_c^(')+ ⟩,
    Ψ_Ξ_c^(')Ξ_c^(')^1,0   =1/√(2)|Ξ_c^(')+Ξ_c^(')0+Ξ_c^(')0Ξ_c^(')+⟩, Ψ_Ξ_c^(')Ξ_c^(')^1,1=|Ξ_c^(')+Ξ_c^(')+⟩, Ψ_Ξ_c^(')Ξ_c^(')^1,-1=|Ξ_c^(')0Ξ_c^(')0⟩,
    Ψ_Ω_cΩ_c^0,0   =|Ω_c^0Ω_c^0⟩,

The wave functions of the bottom baryonium and bottom dibaryon systems can be obtained analogously.



