Critical properties of bound states with one-boson-exchange potential

Hadronic molecular states have attracted considerable theoretical interest in hadronic physics since the experimental discovery of exotic hadron states, such as $ X(3872) $, $ Z_c(3900) $, $ Z_b(10610) $, $ Z_b(10650) $, and the $ P_c $ and $ P_{cs} $ states [1−6]. Theoretically, these states are widely interpreted as hadronic molecular states. This picture has been investigated using various methods, most notably the Schrödinger equation with a one-boson-exchange (OBE) potential [7−12], the unitary approach [13−18], and the Bethe-Salpeter equation [19−26], in which the interactions are based on effective Lagrangians. In particular, because of its intuitive physical interpretation, the potential model has been progressively refined and is frequently used to describe hadronic molecular states, including the deuteron, on the basis of the Yukawa potential. The Yukawa potential was proposed to explain the strong interaction between nucleons through the exchange of a massive particle [27]. It has also been widely applied in various branches of physics, such as plasma physics, nuclear physics, astrophysics, and solid-state physics [28−37].

In the context of potentials for hadronic molecular states, form factors are introduced to describe the interaction between two non-point-like particles, leading to potentials beyond the Yukawa form. These potentials are commonly used to study the binding energies of such systems. Comparisons between theoretical results and experimental data are used to determine the parameters in the form factors. The inclusion of these parameters makes the description more physical, but it also increases the complexity and obscures certain properties.

On the other hand, studies of the Schrödinger equation with the Yukawa potential reveal several interesting and elegant global properties. For a system of two nonrelativistic particles, each with mass $ 2\mu $, whose interaction is described by the Yukawa potential:

where α is the coupling constant and $ m_{ex} $ represents the mass of the exchanged particle. The binding energies of this system can be described by the Schrödinger equation. When $ m_{ex}=0 $, the Yukawa potential reduces to the Coulomb potential, for which an infinite number of bound states exist, as is well established. When $ m_{ex}\neq 0 $, the situation is fundamentally different. For any finite, positive $ m_{ex} $, the number of bound states becomes finite. Moreover, there exists a critical mass beyond which no bound state can form. This phenomenon has been extensively investigated in previous studies [38−44]. These studies found that the existence of an S-wave ground-state bound state depends on the relationship among the coupling constant α, the reduced mass μ, and the exchange mass $ m_{ex} $, as given by the critical condition:

Furthermore, Refs. [39, 44, 45] systematically calculated the critical values of $ m_{ex} $ for various angular momentum states. These studies also provide the critical masses associated with the existence of the ground state, first excited state, second excited state, and higher excited states.

Physically, the number of bound states is a global property of the system and is not determined by any single binding energy. For physical hadronic molecular states, once the structure of the hadrons is taken into account, a natural question is what happens to the critical behavior in these systems. In this study, we focus on the critical properties of hadronic molecular states rather than solely on their binding energies. These critical properties reflect the global nature of the system and can help us understand its internal properties. To simplify the discussion, we consider a single channel as an example.

The remainder of the article is organized as follows. In Sec. II, we present the theoretical framework, including the derivation of the Schrödinger equation and the potentials from the Bethe-Salpeter equation, and describe the numerical method employed. In Sec. III, we present our numerical results, including the critical mass $ m_{ex} $ required for the existence of bound states, the critical relationship between $ m_{ex} $ and α with and without form factors, and the dependence of relevant quantities on the cutoff Λ. In addition, we give the computed binding energies as functions of α and Λ. Finally, Sec. IV summarizes our results.

For simplicity, we first consider a system of two identical scalar particles, each of mass m, whose interaction is mediated by the exchange of a third scalar particle of mass $ m_{ex} $. In quantum field theory (QFT), the two-body bound state of this system is described by the Bethe-Salpeter (BS) equation:

where $ \chi_{\text{BS}} $ is the BS vertex function, $ K_{\text{BS}} $ is the BS irreducible interaction kernel, and $ G(P,k) $ is the product of the propagators of the two scalar particles, expressed as

where

where $ k_1 $ and $ k_2 $ are the momenta of the corresponding particles.

In general, the BS equation has no analytic solution, even for $ m_{ex} = 0 $. When the binding energy is small compared with m, the instantaneous approximation is appropriate, meaning that:

In the center-of-mass frame, where $ P\equiv(M,\boldsymbol{0}) $, the instantaneous approximation and the form of $ G(P,k) $ imply that $ \chi_{\text{BS}}(P,p) $ is independent of $ p_0 $, yielding the following equation:

where

where $ \omega_k=\sqrt{m^2+\boldsymbol{k}^2} $. By defining

Eq. (7) can be equivalently written as

This is the usual form of the Salpeter equation. Since $ E_b\equiv M-2m $ and $ |\boldsymbol{p}| $ are much smaller than m, we can expand $ \bar{G}(M,\boldsymbol{k}) $ to leading order in these two small quantities, yielding the Schrödinger equation in momentum space:

The corresponding Schrödinger equation in coordinate space is therefore given by:

where $ \mu=\dfrac{1}{2}m $ and

with

In the present work, the Schrödinger equation is derived from the BS equation by adopting the instantaneous approximation, which effectively reduces the relativistic two-body dynamics to a nonrelativistic framework. This approximation neglects certain relativistic effects and assumes that the interaction kernel is independent of the relative time. It may become inadequate when relativistic corrections and energy-dependent effects are sizable and are not fully accounted for in the present treatment. Nevertheless, the Schrödinger equation provides a well-established and widely used framework for investigating bound-state systems. In practical applications, some relativistic effects can be effectively absorbed into model parameters, such as an effective coupling constant or the effective mass of the exchanged particle. Therefore, the present approach is expected to provide a reasonable description at the qualitative and semi-quantitative levels.

For point-like scalar particles, taking the interaction vertex for the three scalar fields as $ \Gamma=-2igm $ gives the interaction kernel:

and the corresponding potential in coordinate space is expressed as

where $ r\equiv |\boldsymbol{r}|,\alpha\equiv \dfrac{g^2}{4\pi} $. This corresponds to a pure Yukawa potential. When the particles are not point-like, form factors are typically introduced to characterize their structure, and the corresponding correction in momentum space modifies the kernel as:

where $ i=2 $ and $ i=4 $ correspond to introducing monopole-type and dipole-type form factors at each vertex, respectively. Because the kernel contains two vertices, the choices $ i=2 $ and $ i=4 $ correspond to monopole and dipole form-factor schemes, respectively. The corresponding coordinate-space potentials with monopole and dipole form factors are given, respectively, by

where

The preceding discussion can be extended to systems involving particles with other quantum numbers, for which the situation remains similar. After separating the angular part, as in the hydrogen case, the radial part of the equation in coordinate space is given, as in standard textbooks, by:

where l denotes the angular momentum quantum number, $ \chi(r)\equiv r R(r) $, and $ \psi(\boldsymbol{r})\equiv R(r)Y_{lm}(\Omega_{\boldsymbol{r}}) $.

Since the Schrödinger equation with the general potential described above cannot be solved exactly [42], we employ a numerical shooting method to solve the system with high precision. In the numerical calculation, we adopt the approximation that the wave function $ \chi(r) $ satisfies the following boundary conditions:

where δ and Δ are parameters corresponding to the physical limits approaching 0 and $ \infty $, respectively.

In the shooting method, we first fix $ E_b $ and $ m_{ex} $ and use the conditions $ \chi(\delta)=0 $ and $ \chi'(\delta) = 1 $ to solve Eq. (20) in the region $ r\in[\delta,\Delta] $. This yields the value of $ \chi(E_b,m_{ex},\Delta) $. To determine the binding energy for a known $ m_{ex} $, we solve the equation $ \chi(E_b,m_{ex},\Delta)=0 $ for the root $ E_b $. To determine the critical value of $ m_{ex} $, we solve the equation $ \chi(0, m_{ex},\Delta)=0 $ for the root $ m_{ex} $. We use a Mathematica script to implement this procedure, using the NDSolve function to solve the differential equation and FindRoot to locate the desired root. High accuracy is ensured by setting WorkingPrecision=200 and PrecisionGoal=40.

In the numerical calculation of the binding energy for fixed $ m_{ex} $ and α with $ \mu=1 $, we use the following parameter range:

where units are omitted.

In the practical calculation of the critical value of $ m_{ex} $ with $ \mu=1 $ and $ \alpha=1 $, we set the parameters within the following range:

It is important to note that a significantly larger Δ is required to determine the critical points of the state, because the binding energy is zero. Consequently, the final results are considered acceptable only when they are stable with respect to the selected ranges of δ and Δ. For example, Table 1 lists the critical ground-state $ m_{ex} $ as a function of the integration parameters δ and Δ (cf. the exponents in columns two and three), demonstrating the stability of the results.

Using the numerical method outlined above, we compare our results with the analytical solutions for the binding energies in the Coulomb potential case, with parameters $ m_{ex}=0 $, $ \alpha=1 $, and $ \mu=1 $. The binding energies obtained using the numerical shooting method are consistent with the analytical solutions for the Coulomb potential. Specifically, the comparison shows agreement to 35 significant digits. This high level of precision confirms the validity of our approximation and numerical method. For the case where $ m_{ex}\neq 0 $, we expect this method to achieve similar precision.

For the Yukawa potential $ V_0(r) $, the critical value of $ m_{ex} $ has been studied in Refs. [39, 42−47]. As an initial test and comparison, we present our numerical results for the critical value of $ m_{ex} $ for the ground state in Table 2, with the coupling constant fixed at $ \alpha = 1 $ and the reduced mass fixed at $ \mu = 1 $. The comparison shows that our results agree with those reported in Refs. [45−47] to approximately 30 significant digits, indicating that the shooting method is effective for obtaining high-precision results.

The critical values of $ m_{ex} $ for states with other quantum numbers n and l are presented in Tables 3,4,5,6. Most of our results are consistent with those reported in Refs. [39, 44, 45]. However, significant discrepancies are observed in some cases. Given that the approach employed in this work demonstrates sufficient accuracy for subsequent calculations, we conclude that the results obtained using the shooting method are reliable. These results also reveal a global property: for fixed n, the critical values of $ m_{ex} $ for $ l=1 $ are always smaller than those for $ l=0 $.

In Tables 3–6, we present our numerical results to 35 significant digits. In practical calculations, parameters such as $ \text{WorkingPrecision} $ and $ \text{PrecisionGoal} $ can be increased to achieve higher accuracy. Furthermore, results for other states, such as those with $ l=2 $, can also be readily obtained but are not listed here. The complete results demonstrate that the shooting method can yield highly accurate critical values of $ m_{ex} $ for states with high angular momentum (l). These high-precision values have not previously been reported in the literature.

In Fig. 1, we show how the number of bound states depends on the parameters $ m_{ex} $ and α. The panels, ordered from left to right, correspond to $ l=0 $, $ l=1 $, and $ l=2 $. The red (I), blue (II), and green (III) regions indicate the presence of one, two, and three bound states, respectively, for the given l. The dotted region denotes an area outside the scope of our discussion. Our calculations specifically consider the region where $ m_{ex}\geq0.1 $. The dashed lines show the critical-boundary curves obtained from our discrete numerical points and extrapolated toward smaller $ m_{ex} $ in the $ (\alpha,m_{ex}) $ plane, as justified below. We extrapolate the critical boundaries toward the origin for two reasons. First, when $ m_{ex}\rightarrow 0 $, bound states exist for any sufficiently small positive α, whereas none exist for $ \alpha=0 $. Therefore, both α and $ m_{ex} $ must approach zero along the critical line. Second, within the computed parameter range, we observe an approximately linear relation between the critical coupling α and $ m_{ex} $. Following this trend, the boundary curves are extended to smaller $ m_{ex} $ and visually appear to converge to the origin as $ \alpha\rightarrow 0 $. For the specific case of $ l=0 $ and $ \alpha=1 $, the critical value of $ m_{ex} $ is approximately $ 1.19 $, consistent with the data presented in Table 2. We observe that the slope of the boundary curve increases with the orbital angular momentum quantum number l. This indicates that a larger coupling constant is required to form a bound state as l increases. Furthermore, the linear behavior of the region boundaries directly indicates that the critical values satisfy relations similar to Eq. (2). Such relations can, in fact, be determined through dimensional analysis.

Figure 1.  (color online) Critical boundaries for the number of bound states as functions of α and $ m_{ex} $ for specified angular momentum l in the $ V_0(r) $ case, with μ fixed at $ 1 $. The panels, ordered from left to right, show the cases $ l=0 $, $ l=1 $, and $ l=2 $. The red (I), blue (II), and green (III) regions indicate the existence of one, two, and three bound states, respectively.

The potential $ V_1(r) $ is widely used to study molecular states within OBE models. Its critical parameters reflect the global properties of the system and are crucial for understanding the internal properties of the bound states. In this subsection, we discuss the properties associated with this potential. For simplicity, we continue to use the scalar system as an example.

Table 7 provides the critical values of $ m_{ex} $ for $ \mu=1 $, $ \alpha=1 $, and $ \Lambda=1 $, classified according to the quantum numbers n and l. Because no comparable results for the critical mass in the $ V_1 $ case are available in the literature, Table 7 reports only our values. Fig. 2 illustrates how the number of bound states evolves with the parameters $ m_{ex} $ and α, with Λ and μ fixed, using the same notation as in Fig. 1. Similarly, Fig. 3 shows how the number of bound states evolves with the parameters α and Λ, with $ m_{ex} $ and μ fixed.

Figure 2.  (color online) Critical boundaries of regions with different numbers of bound states as functions of α and $ m_{ex} $ for fixed angular momentum l and cutoff Λ for the $ V_1(r) $ potential, with μ fixed at $ 1 $. The panels, ordered from left to right, show the cases $ l=0 $, $ l=1 $, and $ l=2 $. The red (I), blue (II), and green (III) regions indicate the existence of one, two, and three bound states, respectively.

Figure 3.  (color online) Critical boundaries for the number of bound states as functions of α and Λ for specified angular momenta l in the $ V_1(r) $ case, with μ fixed at $ 1 $ and $ m_{ex} $ fixed at $ 0.1 $. The panels, from left to right, show the cases $ l=0 $, $ l=1 $, and $ l=2 $. The red (I), blue (II), and green (III) regions indicate the existence of one, two, and three bound states, respectively.

The results in Table 7 show an interesting feature that differs markedly from the $ V_0 $ case: for some quantum numbers n, the critical value of $ m_{ex} $ in the $ l=1 $ case is larger than that in the $ l=0 $ case. The results in Fig. 2 also exhibit behavior distinct from that in the $ V_0 $ case: the boundaries of the regions are no longer straight. This indicates that, for finite Λ, the critical values of $ m_{ex} $ and α do not satisfy a linear relation analogous to Eq. (2). In addition, when $ m_{ex} $ is fixed, the corresponding critical value of α decreases as the cutoff Λ increases. This suggests that a smaller cutoff Λ has a more pronounced effect on the existence of bound states.

The results in Fig. 3, with $ m_{ex}=0.1 $ and $ \mu=1 $, show that the critical values of α for the ground state depend only weakly on Λ. This suggests that varying Λ does not substantially alter the number of bound states in the large-Λ regime, particularly for the ground state. To illustrate this explicitly, we provide a short numerical table, Table 8, from which the trend can be read directly.

As in the previous subsection, we use the dipole form factor with the same cutoff, $ \Lambda = 1 $. Table 9 shows that the critical values of $ m_{ex} $ are smaller than those in Table 7. This is expected because, for a given nominal cutoff, the strength and range of the effective interaction depend on whether a monopole or dipole form factor is adopted.

Fig. 4 shows that the dependence of α on $ m_{ex} $ is similar to that in Fig. 2, while Fig. 5 closely parallels Fig. 3 for the dependence of α on Λ. In addition, for the dipole form factor with the same Λ, a smaller $ m_{ex} $ and a larger coupling α are required for bound states to appear. By readjusting Λ, similar critical values can be reproduced within our framework.

Figure 4.  (color online) Critical boundaries for the number of bound states as functions of α and $ m_{ex} $ for a given angular momentum l and cutoff Λ in the $ V_2(r) $ case, with μ fixed at $ 1 $. The panels, from left to right, show the cases $ l=0 $, $ l=1 $, and $ l=2 $. The red (I), blue (II), and green (III) regions indicate the existence of one, two, and three bound states, respectively.

Figure 5.  (color online) Critical boundaries for the number of bound states as functions of α and Λ for specific angular momentum l in the $ V_2(r) $ case, with μ fixed at $ 1 $ and $ m_{ex} $ fixed at $ 0.1 $. The panels, ordered from left to right, correspond to $ l=0 $, $ l=1 $, and $ l=2 $. The red (I), blue (II), and green (III) regions indicate the existence of one, two, and three bound states, respectively.

In this subsection, we extend the preceding discussion to the $ D\bar{D} $ system. The kernels $ K_{D\bar{D}} $ and the potential $ V_{D\bar{D}} $ introduced below follow Ref. [48] for the $ J=0 $ $ D\bar{D} $ channel; they should not be confused with the dipole form-factor potential denoted by $ V_2(r) $ in Sec. 3.3. The energy spectrum of this system was studied in Ref. [48], where the corresponding potential in the Schrödinger equation for the $ J=0 $ state is given as follows:

where $ \alpha \equiv \dfrac{g^2}{4\pi} $, $ \alpha_{\sigma} \equiv \dfrac{g_{\sigma}^2}{4\pi} $, and

We note that the form factor used here is not the dipole form but rather the specific form adopted in Ref. [48].

In Ref. [48], the coupling constants are taken to be $ g_{\sigma}=-0.76 $ and $ g=5.247 $ (or $ \alpha=2.19 $ and $ \alpha_{\sigma}=0.046 $). Since $ \alpha_{\sigma} $ is much smaller than α, we fix its value in the following discussion to study the effects of α and Λ. The masses of the exchanged bosons are taken to be $ m_\rho=0.776 $ GeV, $ m_\omega=0.783 $ GeV, and $ m_\sigma=0.6 $ GeV. The masses of D and $ \bar{D} $ are taken to be $ m_D=m_{\bar{D}}=1.867 $ GeV.

In Fig. 6, we show the critical boundaries for the number of bound states as functions of α and Λ, using the same notation as in Fig. 1. For comparison, the parameters used in Ref. [48] are also indicated by the two black triangles, which correspond to $ \alpha=2.191 $ with $ \Lambda=1.46 $ GeV and $ \Lambda=1.76 $ GeV.

Figure 6.  (color online) Critical boundaries indicating the number of bound states as functions of α and Λ for specified angular momentum l in the $ D\bar{D} $ system. The panels, from left to right, show the cases for $ l=0 $, $ l=1 $, and $ l=2 $. The red (I), blue (II), and green (III) regions indicate the existence of one, two, and three bound states, respectively.

The results in Fig. 6 suggest that the critical value of α for the ground state depends only weakly on Λ when $ \Lambda>2 $ GeV. To illustrate this behavior in detail, Table 10 lists the corresponding values. These results show that the variation of α with Λ is much milder for $ \Lambda\in[2, 4] $ GeV than for $ \Lambda\in[1, 2] $ GeV.

Figure 7 shows the number of bound states N as a function of Λ for $ \alpha=2.191 $, 5, 10, and 20. For coupling strengths typical of the present phenomenological context, N remains equal to 1 for $ \Lambda>1.4 $ GeV. For $ \alpha=5 $, N also remains equal to 1 for $ \Lambda>1.1 $ GeV. For larger values of α, N increases with Λ.

Figure 7.  (color online) Number of bound states N as a function of Λ for selected values of α at $ l=0 $ in the $ D\bar{D} $ system.

In Fig. 8, we show the binding energy $ E_b $ as a function of the coupling constant α for fixed Λ in the $ l=0 $ case. The parameters used in Ref. [48] are also indicated by the two black triangles. The notations $ 1^{\text{st}} $, $ 2^{\text{nd}} $, and $ 3^{\text{rd}} $ refer to the first, second, and third bound states with $ l=0 $, respectively. The results show that $ E_b $ is nearly linear in α. As α increases, additional bound states appear, demonstrating that the magnitude of the coupling constant plays a critical role in determining the number of bound states at a fixed cutoff Λ.

Figure 8.  (color online) The binding energy $ E_b $ as a function of α for fixed Λ in the $ l=0 $ case of the $ D\bar{D} $ system. The panels, ordered from left to right, show the cases for $ \Lambda=1.46 $ GeV, $ \Lambda=1.76 $ GeV, and $ \Lambda=5 $ GeV. The red, blue, and green lines denote the first, second, and third bound states, respectively. The two black triangles correspond to the results obtained using the parameters given in Ref. [48].

In Fig. 9, we show the binding energy $ E_b $ as a function of the cutoff Λ for fixed coupling constant α and angular momentum l. The results indicate that $ E_b $ is nearly linear with respect to the cutoff Λ at large Λ. When α is fixed at $ 2.19 $ (the value used in Ref. [48]), only one bound state is obtained.

Figure 9.  (color online) Binding energy $ E_b $ as a function of Λ for fixed angular momentum l and coupling constant α in the $ D\bar{D} $ system. The panels, ordered from left to right, show the cases for $ \alpha=2.19 $, $ \alpha=5 $, and $ \alpha=10 $. The red, blue, and green lines denote the first bound state with $ l=0 $, the second bound state with $ l=0 $, and the first bound state with $ l=1 $, respectively. The black triangle indicates the results obtained using the parameters given in Ref. [48].

For comparison, we apply the same interaction parameters directly to the $ B\bar{B} $ system. The corresponding results for the number of bound states N are shown in Fig. 10 and exhibit trends similar to those in the $ D\bar{D} $ case.

Figure 10.  (color online) The number of bound states N as a function of Λ for selected values of α at $ l=0 $ in the $ B\bar{B} $ system, with the interaction parameters taken from the $ D\bar{D} $ analysis as a test.

Taken together, these results indicate that the existence and multiplicity of bound states are strongly correlated with the coupling constant α, particularly when α is small, but depend only weakly on the cutoff Λ over a broad region.

In this study, we investigate the critical behavior of bound states in the Yukawa potential and a modified version of it, both of which are widely used in the OBE model for hadronic molecular states. The analysis is based on a highly accurate numerical method. Our numerical result for the critical mass $ m_{ex} $ of the ground state is consistent with values reported in the literature. Furthermore, our results for the $ l=1 $ case are significantly more precise than previously reported values. For the $ D\bar{D} $ hadronic molecular system, our results indicate that the number of hadronic molecular states is only weakly sensitive to the cutoff Λ over a broad range (for example, $ \alpha\le5 $ and $ \Lambda>1.5\; \mathrm{GeV} $), although the binding energies may still depend on Λ, as discussed above. This behavior suggests an interesting global regularity in hadronic molecular states. Finally, this highly accurate numerical method can be straightforwardly applied to higher-l cases or other systems.