Closed-form solutions for the Schrödinger wave equation with non-solvable potentials: A perturbation approach

1 Introduction

One of the challenging problems in non-relativistic quantum mechanics is to find closed-form solutions to the radial Schrödinger wave equation for non-solvable potentials that are used in different fields of physics. In particular, the perturbed Coulomb potentials are representing the most important potentials in many scientific models, especially in particle physics [1].

As an example, the experimental observations of the exotic states in the charm sector with c ¯ c c ¯ c content at Belle [2] and BESIII [3], as well as the structures of the form Q ¯ q q ¯ Q in ref. [4,5], have opened new interest to re-study the hadronic structures and their spectroscopy with different phenomenological potential models.

However, the lack of closed-form solutions to the radial Schrödinger wave equation with many potential models makes the study of such potentials one of the most popular theoretical laboratories for examining the validity of the various approximation techniques based on perturbative and non-perturbative approaches [6,7,8, 9,10,11].

Many of the existing theoretical work in the literature can obtain the hadron masses by solving numerically the two-body radial Schrödinger wave equation with Cornell-inspired potentials, and then the respective spin-dependent interaction terms are added and treated by the standard perturbation technique [12]. Nevertheless, the computational schedule in those works is very complicated and can be hardly extended to study other potential models.

On this basis, the present work is aimed to solve the radial Schrödinger wave equation for heavy quark–antiquark structures with the vector Coulomb potential plus linear and harmonic scalar potentials parts, implemented in the same form with their counterpart spin-dependent terms.

We will show how the PAIM can be used directly to find analytically the coefficients in the energy perturbation series as a solution to the radial Schrödinger wave equation, without the need of using the base eigenfunctions of the unperturbed problem. Furthermore, PIAM is unlike the other perturbation methods, with no constraints on the coupling constants or the quantum numbers involved in the phenomenological potential, and even it avoids sticking to numerical computations at the early stages of its procedure. Besides, PIAM can yields very accurate and rapidly converging eigenvalues, and it is applicable in the same form to both the ground and excited energy states.

With this in mind, this article is organized as follows. Section 2 is devoted to the theoretical framework of the phenomenological potential model. We presented the static potential, which is composed of scalar and vector parts, implemented with their counterpart spin-dependent correction terms. In Section 3, the formulation of the PAIM for finding the energy eigenvalues is given. The analytical expressions for PAIM are cast in such a way that allows the reader to use them without proceeding into their derivation. In Section 4, as an example, we tested numerically the method to find some of the P-wave mass spectrums of the charmonium c c ¯ , and bottonium b b ¯ mesons, and therein for comparison we presented the available experimental and theoretical results in the literature. Finally, we give a summary and concluding remarks for possible future works on tetraquark masses.

2 Theoretical framework: Implementation of the spin-dependent terms within the phenomenological potential parts

It is generally believed that the phenomenological potential models used in the hadronic spectroscopy with the Coulomb potential perturbed by a term or terms involving various powers of r seem to provide the clearest and most accurate picture for meson states. The Coulomb potential arises from the one-gluon exchange associated with a Lorentz vector structure, and the various powers of r parts in the potential are responsible for the confinement. These parts are known as the static spin-independent potential parts. For our further discussion, we used a fairly flexible mixing type of such potentials. We took the potential interaction parts as a sum of vector and scalar parts, plus the spin-dependent terms:

with

and

where, V V ( r ) , V S ( r ) , and V S D ( r ) are the vector, scalar, and spin-dependent parts, respectively. Λ , K , and σ represent the independent coupling constants of the potential. In principle, the spin-dependent interaction terms are usually added to account for the mass splitting between states with different quantum numbers. Here, β is the mixing coefficient. Additionally, this potential has distinct features of strong interaction quantum chromodynamics, the asymptotic freedom, and the confinement parts.

To implement the spin-dependent terms within the pre-assumed phenomenological potential parts, we based on the Breit-Fermi one-gluon exchange Hamiltonian [13], where the spin-dependent part for quark masses m Q and m Q ¯ have three distinctive types of radial-interaction parts known as the spin–spin V S S ( r ) , the spin–orbit V L S ( r ) , and the tensor part V T ( r ) , respectively:

Performing the derivatives, simplifying the resulting expressions, and collecting terms of the same powers of r , then we immediately obtain a phenomenological potential model implemented with the spin-dependent terms:

where we have used ∇ 2 1 r = − 4 π δ 3 ( r ) , and

and

Here it should be noted that the obtained phenomenological potential in equation (8) is elaborated in such a way to be comfortable with the application of the PAIM. Furthermore, when ℓ = 0 this leads to z ˜ = Λ − 4 β σ 3 m Q m Q ¯ ⟨ S 1 ⋅ S 2 ⟩ , λ = 0 , and ω ˜ 0 = 4 β K m Q m Q ¯ ⟨ S 1 ⋅ S 2 ⟩ . In this case, the phenomenological potential model in equation (8) is reduced to a potential model, which is responsible for generating the S-wave states:

while, when ℓ ≠ 0 then δ 3 ( r ) = 0 , since only S-wave states have a non-zero value of the wave function at the origin. For this case, the phenomenological potential model in equation (8) is reduced to a potential that is responsible for generating the P-wave, D-wave, etc. states:

In principle, from now on, and for our future discussion, we will only solve the radial Schrödinger wave equation with the potential model of equation (13) for ℓ ≠ 0 , which represents the most general case.

The expectation values of the spin–spin interaction operator in equation (4) are calculated in terms of the spin quantum numbers: ⟨ S 1 ⋅ S 2 ⟩ = 1 2 ( S 2 − S 1 2 − S 2 2 ) , where S 1 and S 2 are the spins of the two quarks m Q and m Q ¯ , respectively, and S stands for the total spin of the structure under consideration. Since each quark (anti-quark) has spin 1 2 , thus they can couple into singlets S = 0 with ⟨ S 1 ⋅ S 2 ⟩ = − 3 4 or triplet S = 1 with ⟨ S 1 ⋅ S 2 ⟩ = 1 4 .

On the other hand, the expectation values of the spin–orbit interaction operator are calculated in terms of the total angular momentum quantum number J , defined by the vector sum J = L + S : ⟨ L ⋅ S ⟩ = 1 2 ( J 2 − L 2 − S 2 ) . For S-wave states ( ℓ = 0 ) and/or for the singlet state case S = 0 , the spin–orbit term ⟨ L ⋅ S ⟩ is always zero. However, for the triplet state cases S = 1 with ℓ ≠ 0 , we have J = ℓ + 1 , ℓ , ℓ − 1 , which, in general, can yield a non-zero value. Similarly, the expectation values of the tensor operator are calculated in terms of the total spin S , and the orbital angular momentum quantum number ℓ , and it is given by ref. [13]:

The results of this operator are equal to zero only for S = 0 or ℓ = 0 , and its expectation values are non-zero for the triplet states S = 1 with ℓ ≠ 0 , and J = ℓ + 1 , ℓ , ℓ − 1 :

3 Formalism of the PAIM for the potential V Q Q ¯ ( r ) for ℓ ≠ 0

The radial Schrödinger wave equation that describes the internal structure between any quark–antiquark with the phenomenological potential model in equation (13) for ℓ ≠ 0 (in atomic units, ℏ = c = 1 ) is as follows:

where μ = m Q m Q ¯ m Q + m Q ¯ is the reduced mass of the quark–antiquark system, n r is the radial quantum number, and for states of definite orbital angular momentum l the radial wavefunction is ψ n ℓ ( r ) = r − 1 R n ℓ ( r ) .

The solution of equation (16) within the framework of PAIM is started with the change of variables by letting r = t 2 , and further we removed the first derivative by employing the ansatz Φ n r ℓ ( t ) = R n r ℓ ( t ) exp ( − p ( t ) / 2 ) , provided that p ′ ( t ) = − 1 t , which in turn implies a Schrödinger-like radial equation:

with ℓ ́ ( ℓ ́ + 1 ) = 4 ℓ ( ℓ + 1 ) + 3 4 . For further simplifications to equation (17), we employed the transformation Φ n r ℓ ( t ) = t ( ℓ ́ + 1 ) e − α t 4 / 4 f n r ℓ ( t ) , simplifying out the mathematics, collecting the terms of equal powers in t , and choosing an arbitrary parameter α = 8 μ K in such a way to remove the terms of higher t powers. This procedure leads to another second-order homogenous linear differential equation of the form [11,14]:

where

with ε n r ℓ = 8 μ ( E n r ℓ − ω 0 ) − ( 2 ℓ ́ + 5 ) α , and γ is an artificially introduced perturbation expansion parameter to be equal to one at the end of the calculations.

It should be noted here that the choice of the ansatz in equation (19) is motivated, that is, when we switch γ off, equation (19) is reduced to an exactly solvable eigenvalue problem within the framework of AIM [14,15]. On the other hand, if we switch γ on, a straightforward application of AIM to equation (19) does not give any reasonable results, where it is observed that the sequence appears not to converge, a result that violates the principle behind the AIM.

Accordingly, to have a closed-form solution for equation (19) within the framework of PAIM, we proposed s 0 ( t , γ ) be expressed as a sum of two parts:

with s 0 0 ( t ) = − ε n r ℓ t 2 which has an exact solution within AIM, and s 0 1 ( t ) = 8 μ σ t 4 + λ t 4 − z to be solved using PAIM.

Following the systematic procedure of PAIM, the energy eigenvalues E n r ℓ are solved by imposing the asymptotic aspect of the method for sufficiently large iteration number i :

with

and

As a step further, to obtain the energy eigenvalues for s 0 1 ( t ) , we expand δ i ( t , γ ) around γ = 0 :

and according to the procedure of PAIM, δ i ( t , γ ) must be zero; if this to be true for every γ value, then every term in the series must be zero. That is to say:

At this stage, it is also convenient to expand ε n r ℓ in terms of γ :

A quantitative estimation for ε n r ℓ expansion terms are obtained by comparing the terms of equation (24) with the same order of γ as the terms of equation (26). Therefore, it is obvious from the procedure of AIM that the roots of δ i ( 0 ) ( t , 0 ) = 0 give the leading contribution energy terms ε n r ℓ ( 0 ) . Likewise, the roots of δ i ( 1 ) ( t , γ ) = 0 give the first correction terms ε n r ℓ ( 1 ) and so on.

To obtain the analytical leading energy term ε n r ℓ ( 0 ) , one should simply switch off γ in equation (19), and for further simplifications, we recall the mathematical definitions of equations (21)–(23). Then the roots of δ i ( 0 ) ( t , 0 ) = 0 leads to an exactly solvable eigenvalue in terms of α :

It should be addressed here that for each iteration the expression δ i ( 0 ) ( t , 0 ) = 0 depends on t , and the calculated eigenvalues ε n r ℓ ( j ) using this condition should, in principle, be independent of the choice of t . Nevertheless, the choice of t is critical only to the stability of the AIM process. In this formalism, it is observed that the starting value for t is the value at which λ 0 ( t , 0 ) = 0 , so that t = t 0 = ℓ ́ + 1 2 α 1 / 4 , and the sequence appears to converge when the number of iterations i ≥ 25 , as shown in Table 1.

To find the other analytical perturbed expansion terms, one should simply go to equation (19) switch γ on, therein, replace ε n r ℓ with ε n r ℓ ( 0 ) + γ ε n r ℓ ( 1 ) , and terminate the iterations by imposing the condition δ i ( 1 ) ( t , γ )  = 0. The first root of the resulting equation gives ε n r ℓ ( 1 ) . Proceeding exactly in the same way, one can obtain the second-order closed-form perturbed expansion terms and so on. The analytical expressions of this procedure are displayed in Table 1. Of course, in Table 1, z and λ are changeable according to the total spin and total angular momentum of the structure state.

Consequently, the final closed-form solution for the energy eigenvalues E n r ℓ in conjunction with equation (19) within PAIM framework is as follows:

where ξ = 8 μ γ .

A further interesting quantity in this respect is the transformed eigenfunctions f n r ℓ ( t ) . Knowledge of f n r ℓ ( t ) will help to construct the radial eigenfunctions R n r ℓ ( r ) . Inspired by equation (18), given that ϱ i ( t , γ ) = s i ( t , γ ) λ i ( t , γ ) , and following again exactly the same systematic procedure of PAIM explained above, by expanding this time ϱ ( t , γ ) :

This procedure allows us to find the coefficients of the transformed eigenfunctions:

and finally this will lead to a general solution to the transformed second-order differential equation (18):

where N is an integration constant.

4 Numerical results for heavy mass spectrum of quark–antiquark structures

Having solved the radial Schrödinger wave equation with the proposed phenomenological potential model, the corresponding mass of a particular quark–antiquark Q ¯ Q structure within this model can be written as:

where the energy eigenvalues E n r ℓ are given in equation (28). Typically, to obtain numerically the meson masses, this approach needs various input parameters that one should fix at this stage. For the quark masses, we used m b = ( 4.992 ± 0.001 ) GeV , and m c = ( 1.55 ± 0.001 ) GeV . Furthermore, the PAIM approach contains four auxiliary parameters with their working regions to be determined, namely, the Coulomb gluon strength part Λ , the linear and harmonic confining parts σ and K , and the mixing parameters β . Our analysis shows that Λ , σ , K , and β = 0.25 are varied in the following regions, respectively, Λ = 0.73 ± 0.01 , σ = 0.018 ± 0.001 GeV 2 , K = 0.026 ± 0.001 GeV 3 , and β = 0.25 ± 0.01 .

Few P-wave mass spectrums of the charmonium c c ¯ are calculated and listed in Table 2. The calculated masses are found to be in good agreement with the recently observed ones. Moreover, with the same input parameters, in a unified way, we extend this study and we predicted the P-wave bottonium b b ¯ mass spectrums. In Table 3, we present our final findings on the P-wave mass spectrums of bottonium b b ¯ . In the tables, we reported also on the predictions of the other methods together with the experimental results. A direct comparison of our results with the experimental results reveals that they are in good agreement.

Evidently, our numerical analysis shows that the contribution of the non-perturbative terms is about 90%, while the maximum contribution of the perturbative terms is about 10% of the total sum.

Finally, although in this work we have tested with PIAM the P-wave mass spectrum of the potential model of equation (1), one can easily find the analytical and numerical solutions for other potential models using the same procedure. It is also important to emphasize that the PAIM described in this article and its analytical expressions are much more useful than pure numerical calculations.

5 Conclusion

The present work is focused on the theoretical implementation of the PAIM and to show how it is simple and easy to solve the radial Schrödinger wave equation with non-solvable potential models. Although we aimed not to make the fittings with the available experimental we only emphasized on the method and how it works, and we believe that the interested readers if necessary can do better fitting. Also, we aimed to show how it is an easy task to use the PAIM without having to worry about the ranges of the couplings in the introduced potential model. Moreover, the degree of precision of the results can also be drastically improved by raising the perturbative order in the energy expansion series one step more when it is necessary without any technical difficulties. We believe that this method and its results will provide a valuable future hint for investigating tetraquark configuration.