1.Introduction

The study of heavy quarkonia in quantum chromodynamics (QCD) is a field theoretical non-perturbative problem which has promoted the development of several techniques appropriate to address this issue. Under certain considerations, it can be simplified to a non-relativistic (NR) quantum mechanical problem by exploring the interaction strength between quarks inside hadrons in the static limit of the interaction potential (SLIP) between quarks. Lattice QCD simulations ^1,2 suggest for this potential a Coulombic term at short distances plus a linearly rising part at large inter quark separations, which to-date has motivated a large number of phenomenological models (3)-(19) that capture the quantum chromodynamic traits of strong interactions. Whilst the SLIP picture for light-quarks must be taken with a grain of salt, for the case of heavy mesons, cc¯ and bb¯, this picture is accurate enough to describe the mass spectra of heavy quarkonia. Assuming spherical symmetry, the SLIP can be cast in the general form V(r)=-Ar-α+Brβ+V0, where r is the inter quark distance, A, B and V₀ are constants and α and β are free parameters assessed either from Lattice QCD, see for example , or fixed by fitting the known experimental masses of hadronic states. The choice α = β =1 and V₀ = 0 corresponds to the Cornell potential. The choice V₀ = 0, α = 2 =β/2 corresponds to the anharmonic potential ²¹ and other choices give rise to commonly used potentials.

The relevant Schrödinger wave equation (SWE) with any of these SLIPs often has to be solved numerically. To that end, a vast number of numerical strategies have been implemented in literature among which, to count a few, we find the shooting method ²², the Asymptotic Iteration Method (AIM) ²³ and many others such as different forms of Runge-Kutta methods. Most numerical strategies implemented for solving the SWE are as precise as going up to O(d²) of grid spacing d. The Matrix Numerov Method (MNM) (see, for instance ²⁴) which in the context of meson physics has been put forward by the Qena group ²⁵ and extended by our group ²⁶, departs from a discretization of the kinetic term in the SWE in such a manner that the problem of solving the radial SWE is cast in the form of a matrix eigenvalue problem. In this form, its accuracy is of O(d⁶). On the other hand, many proposals of quark model (QM) potentials have used quark masses far higher than those reported in the Particle Data Booklet in various years ^3-19. The quark model community has exploited the freedom in choosing quark masses (the so called constituent quark masses) and quark interaction energy (QIE) parameters in mq1+mq2+Σ=Mmeson. This hides the factual and physical values of quark masses mq1,mq2 and those of the quark interaction energy Σ.

In our present work we study the spectrum of heavy cc¯, bb¯ states. For that purpose, we use bare quark masses in SWE extracted from experimental evidence got by various Collaborations. Regarding the SLIP, we are interested in the choice of parameters α = β= 1/2 and V₀ = 0, known in literature as the Song-Lin (SL) potential ²⁷. We select this effective potential beyond the paradigmatic Cornell potential because there are a number of appealing features ^28,29 that the SL potential displays better than the Cornell model. Let us recall that the SL potential has been theoretically calculated by coupling the non-Abelian character of QCD to a heavy dilaton field . From phenomenological perspective ²⁷, the first term in this potential b/r is inspired from leptonic widths of vector mesons p, ω, Φ, J/ψ and Υ measured in experiments ³¹ and perturbative QCD at short distances, while the term -ar) has been motivated by considering a chromo-electric flux tube ²⁷ between a quark Q and anti-quark Q¯ in a quarkonium system where perturbative QCD is hard to manage, and the possible colour-screening effects on strong force which are considered of considerable importance ^33-36. In fact, the colour-screening in our parametrization of SL is far more simply added than in ^33-36 using Cornell potential. The SL potential has been considered in ^23,37,38 and applied only to low lying meson states with no or insignificant application to decay rates and other quarkonia observables. We calculate spin averaged masses of cc¯,bb¯ spectroscopic states, which are relevant for the prediction of hyperfine splittings between |1LL⟩ and |3LL⟩ states defined as ⟨3LL⟩ - 1LL where ⟨3LL⟩ is spin averaged mass of the triplet {3LL-1,3LL+1,3LL+3}. By considering the spin averaged wave functions at contact, we readily obtain the leptonic, diphoton, digluon, triphoton, trigluon and one photon plus 2 gluons decay rates of S-wave states by demanding agreement of 1S state for each one of these decays with the experiment. We do this to find the value of the QCD coupling constant (α_s) and then use it for decay modes of 2S to 5S states. A comparison with known values of these decay widths is presented, rendering our full dynamical picture of the NR framework in fair agreement with the dynamics of these heavy meson states. We also include momentum width and velocity corrections to these decays and get their excellent agreement with experiment. We obtain wavefunctions, masses, radii, quark velocities, quark momentum widths, annihilation decay rates and strong running coupling constant for charmonium and bottomonium states in the following manner: we fit the parameters of the potential by demanding exactness of the theoretical mass of only one 1S spin averaged state of cc¯ with experimental value and then calculate all remaining spin averaged masses for heavy mesons with the parameters α = 0.7011 GeV^3/2, b = 0.8912 GeV^1/2 and this parametrization of SL potential makes the SWE a robust dynamical equation of heavy quarkonia even beyond the original expectations of the model ²⁷ and extends the achievements of other similar calculations ^23,37,38. We compare our findings against experimental results and other theoretical calculations and systematically extend our results up to 8S, 7P, 7D, 7F, 7G. This extension to higher states is allowed and reliable because suitable colour-screening effects are already present in confining term of our SL potential parametrization and this colour screening of the confining term in inter quark potential is seen to produce higher orbital and angular momentum excited states efficiently ³³.

After preparing the firm ground from MNM we obtain the charmonium and bottomonium spin averaged masses and radii, in order to circumvent the struggle for numerical accuracy in the numerical solution, we introduce a variatonal approach combining the Heisenberg uncertainty principle and the variational principle in the Hamiltonian of the SWE. Such an approach, which we refer to simply as the variational principle (VP), imposes an algebraic transcendental relation between the masses and the radii of these states and is partially based in our previous work in ²⁶. The results from our VP are less uncertain than those from MNM. In the remaining part of this article we present the details as follows: Section 2 introduces MNM to solve the SWE. Section 3 presents our model of the quark velocities, and the decay widths for different decay modes of spin averaged S-states of charmonia, and bottomonia. In Sec. 4 we develop a strategy of combining the variation in mesonic Hamiltonian with Heisenberg UP through our momentum width-model to establish analytical relation between spin averaged masses and radii of heavy meson states. Conclusion and further possible developments are discussed in Sec. 5. In the Appendix A to this paper we have calculated spin averaged masses by including hyperfine, spin-orbit, and tensor interactions to SL potential and compare these against the spin averaged masses in Tables I and II obtained without these spin effects while in Appendix B we have estimated the bare quark masses. The color screening is explored in Appendix C.

2.Numerov method

The non-relativistic kinetic energy operator for a two quark system is

where m1,2 and p⃗1,2 are the corresponding quark masses and momenta. Then the Hamiltonian of this system is

If these two quarks are moving in a spherically symmetric potential that depends only on their mutual separation r=|r⃗1-r⃗2|, namely V(r⃗1,r⃗2)=V(r), the problem can be decoupled into the center of mass motion of the two quark system and the relative motion of the quarks about their center of mass. Factorizing the former, the SWE for the relative motion reduces to the one-body eigenvalue problem

where μ=m1m2/(m1+m2) is the reduced mass of the system, Δ is the mass of the meson state and p⃗ is the momentum of quark about the center of mass of the meson. Working in spherical coordinates and separating the angular part in Eq. (3) according to Ψ(r⃗)=Rnl(r)Ylm(θ,ϕ), where Ylm(θ,ϕ) are the spherical harmonic functions normalised to unity (with n, l and m the principal, orbital angular momentum and magnetic quantum numbers), the masses of heavy meson states are obtained by solving the radial one-body, one dimensional-like SWE ψ(r)=rRnl,

where we define

as an operator for QIE. In what follows, we obtain the charmonium and bottomonium spectra by solving the SWE with the SL potential ²⁷

having a and b as free parameters. We fix these parameters to the lightest cc¯ mass, and then derive the rest of the heavy quarkonium full spectrum through the Matrix Numeroved Eq. (10). Explicitly, these values are

whereas the bare quark masses are taken to be

In this form we parametrize charmonium mass spectrum as Δ2S+1LJ=2.4+Σ while for bottomonia we take Δ2S+1LJ = 9.334+ Σ. Our parametrization in Eq. (7) is unique and the same for cc¯ and bb¯ mesons. The masses in Eq. (8) are not arbitrary, see Appendix B for this point.

In order to solve the corresponding SWE in Eq. (4) we discretize its radial coordinate 𝑟 into N equidistant points r_i separated a distance d within a preselected characteristic length r_max suitable for the description of heavy meson spectra and represent the kinetic operator in Eq. (4) in terms of the two trigonal matrices

such that Eq. (4) is expressed in the form as in Refs. ^24-26

where ψi=ψ(ri) is a column matrix (ψ1,ψ2,...,ψN)T, the term in square brackets is a diagonal matrix of order N x N and I0,I-1,I+1 are respectively the square matrices of order N x N containing zeros every where except in the principle diagonal (PD), in the diagonal one step below PD, and in the diagonal one step above PD which are filled with entries as 1. Concretely, we use rmax=4 fm and fix the number of grid points with N = 350. These two numbers are obtained from the stability of the inverse of the masses of S, P and D states against variation of the number of points and the size of the domain of integration. For the sake of illustration, we show the variation of the inverse of the masses for S-states for charmonium and bottomonium with the number of grid points in Fig. 1 and with the interval of integration in Fig. 2.

Figure 1 Stability test of the number of grid points for the prediction of the (inverse) mass of S-states for cc- and bb- at fixed r_max = 4fm. 

Figure 2 Stability test of the size of integration interval for the prediction of the (inverse) mass of S-states for cc- and bb- at fixed N = 350. 

2.1.Spectra

Numerical solution of Eq. (10) gives wavefunctions corresponding to spin averaged masses of charmonium and bottomonium. These are our well behaved spin averaged orthonormal wavefunctions plotted for S, P and D states in Figs. 3 and4. The spin averaged masses are shown in Tables I and II and compared against the experimental values ³⁹ and other theoretical findings. It is important to note that finding the correct interaction energy between quarks (QIE) is a key feature for the success or failure of a quark model and is a real essence of understanding the matter in visible universe. The authors in Ref.²³ and Ref.²⁷ used charm quark mass 1.8 GeV which is 600 MeV more than our bare charm quark mass. It means that these authors estimated interaction energy 600 MeV smaller in charmonium system than ours. As for the bottom quark is concerned, these authors estimate the interaction energy to be 1066 MeV smaller when compared to ours. The correct input mass to SWE produces correct predictions and bare quark masses available so far are not less than being the correct masses of charm in Eq. (B.1) and bottom in Eq. (B.2). So our calculated interaction energies for charmonia and bottomonia must be correct in our bare quark model (BQM). A quick analysis of QIE in Tables I and II shows interaction energy in nL spin averaged states for bottomonia is smaller than that for charmonia, and QIE for (n + 1)L state is larger than for nL. As a complementary note, in the Appendix A we have calculated these mass spectra by including spin-spin, spin-orbit, and tensor interactions between quark Q and antiquark Q¯ as part of the dynamical Eq. (10) instead of using the leading order perturbation method. These are the so called quarkonia mass spectra with spin effects. We have compared these spectra with those in Tables I and II obtained without incorporating the spin effects in the form of mass-plots Fig. 7.

Figure 3 Wave functions for first five S-states (upper panel), P- states (mid pannel) and D-states (lower panel) of cc- mesons in the SL potential 

Figure 4 Wave functions for first five S-states (upper panel), P- states (mid panel) and D-states (lower panel) of bb- meson in the SL potential. 

As for the size of the bound states, the root-mean-square radii rrms can be directly obtained from the numerical solution of the SWE as,

(rrms)2=∫0∞ dr r3|Rnl(r)|2, (11)

where the symbols n and l stand for the principal and orbital angular momentum quantum number of the meson, respectively. Furthermore, the heavy quarkonia states exhibit momentum width β as function of the quarkonium size determined by their quantum numbers in our proposed functional form ⁱ

β=δrrms. (12)

Radii from Eq. (11) and momentum width from Eq. (12) for different states are reported in Table III. A variation of momentum width with the size of the S, P and D states of charmonia and bottomonia are depicted on right panels in Figs. 5 and6, along with the analytical expressions for momentum width as functions of instantaneous quark separation r. From Table III we observe that bottomonia are more compact objects than corresponding charmonia. Furthermore the charmonia sizes are roughly twice as large as those of the similar states of the bottomonia while momentum width for c-quark mesons is approximately half that of the b quark.

Figure 5 The velocities and momenta of charm quark within cc- versus radii of cc- obtained from SWE and our VP. 

Figure 6 The velocities and momenta of bottom quark within bb- versus radii of bb- obtained from SWE and our VP. 

Table III Radii rrmscc- of cc- and rrmsbb- of b¹b states in [fm] and corresponding momentum width βcc- and βbb- in [GeV] in spin averaged S -, P- and D quarkonia states. 

nL	rrmscc- (SWE)/ rrmscc-(VP)/Re. (63)			βSW ECC-/βVPcc-		rrmsbb- (SW E/ rrmsbb-(VP)/Ref.(64)			βSW Ebb-	βVPbb-
2S	0.915/	0.922/	0.91	0.765/	0.760	0.488/	0.501/	0.545	1.434/	1.397
3S	1.352/	1.363/	1.38	0.813/	0.810	0.737/	0.755/	0.805	1.493/	1.457
4S	1.762/	1.773/	1.87	0.851/	0.846	0.972/	0.991/	1.030	1.544/	1.514
5S	2.151/	2.161/	2.39	0.882/	0.879	1.200/	1.216/	1.232	1.588/	1.562
1P	0.697/	0.775/	0.71	0.717/	0.645	0.370/	0.416/	0.435	1.350/	1.202
2P	1.155/	1.195/	1.19	0.779/	0.753	0.628/	0.658/	0.711	1.432/	1.368
3P	1.577/	1.601/	1.67	0.824/	0.812	0.869/	0.893/	0.945	1.496/	1.456
4P	1.976/	1.992/	-	0.860/	0.853	1.097/	1.118/	1.154	1.549/	1.521
5P	2.343/	2.367/	-	0.891/	0.887	1.317/	1.336/	1.346	1.594/	1.572
1D	0.936/	1.097/	0.96	0.748/	0.638	0.507/	0.601/	0.593	1.379/	1.165
2D	1.380/	1.472/	1.44	0.797/	0.747	0.760/	0.818/	-	1.448/	1.345
3D	1.791/	1.850/	1.94	0.838/	0.811	0.994/	1.036/	-	1.508/	1.448
4D	2.178/	2.220/	-	0.871/	0.856	1.218/	1.251/	-	1.559/	1.560
5D	2.497	/2.581/	-	0.900/	0.891	1.434/	1.460/	-	1.604/	1.575

3.Other quarkonia observables

From the discussion of the previous section, we obtain a few parameters that allow us to characterize the dynamics of these heavy quark systems.

4.A Variational Approach

A subtle combination of the Heisenberg uncertainty Principle and the variation of mesonic Hamiltonian can be used to estimate the masses of charmonium and bottomonium as follows. The radial Eq. (4) can be written as,

This means radial Hamiltonian becomes (because radial component of momentum operator is p^r=-i[1r][∂∂r]r),

We have modeled the quark momentum of this equation in a peculiar way which is described next.

We define symmetric point as the one for which all three cartesian components of any 3-vector are equal in magnitude. This means in cartesian space,

This definition turns Eq. (30) into the following effective Hamiltonian:

Here we invoke the Heisenberg uncertainty Principle,

as follows: A little cerebration reveals that the equality sign in Eq. (34) would hold for ideal systems (for which experimental as well as theoretical uncertainties are frozen at the minimum possible values of Δx and Δpx), while the inequality sign in Eq. (34) would come when the system experiences all sorts of “dissipative and interactive" effects. From Tables II andIII, we observe that even in the lowest 1S state the QIE is not zero, which means QQ¯ system always has medium in which quarks move. Therefore, the equality sign in Eq. (34) cannot hold. With these key ideas on board, we parametrize the uncertainty relation in the form (that will become clear below)

where δ is real, positive number more than one and we define it as,

where δn,1-l is the Kronecker delta. This definition is in accordance with the momentum width of known botommonium and charmonium states (see, for instance, Refs. ^25,26,67) and is convenient for the discussion below, as we shorty explain. By keeping in mind the quantum indeterminacy in position and momentum of quark in a meson, let us define the notation in which Δx is x¯ -the mean of the position co-ordinate of quark relative to center of mass of the meson, and similarly Δpx by the mean quark momentum width β¯x about center of mass. With these definitions, whatever physical quantity A we calculate in our VP corresponds to the average ⟨A⟩ of many measurements of a and in essence the same as expectation value of a Hermitian operator A^ from some linear vector space acting on a Hilbert space of states the meson could be in.

Thus, from Eq. (35),

The use of Eq. (31) and Eq. (32) gives,

and, in our notation, Eq. (33) becomes

Now we use this expression for computing meson radii and masses in the follwing way. Equation (37) makes Eq. (39) altogether a function of one variable x¯ only,

This is the Hamiltonian of our meson system QQ¯ that depends on x¯, which we regard as a variational parameter. Thus, minimizing H¯R with respect to x¯, we reach to the constraint

Its solution gives x¯=xmin, and when this value is substituted in Eq. (40) then we identify H¯R(xmin) as the mass of the meson denoted by ΔVP written as,

Let us apply this general procedure to a concrete example in which SLIP is the Song-Lin potential,

Upon inserting this potential into Eq. (40), we reach at the transcendental relation

where y=x¯. Equation (44) has only one real root for which H¯R(x¯) is minimum and which we denote as ymin. All other four roots are complex and hence discarded being extraneous. Thus solution of quintic equation gives,

Equation (45) used in Eq.(31) with our SL parametrization Eqs. (7) and (8) give the radii of cc¯ and bb¯ mesons reported in Table III along with root mean square radii obtained from SWE using Eq. (11) and compared with Refs.^44,45. The agreement of our VP with SWE for quarkonia radii is excellent.

Substitution of the root Eq. (45) in Eq. (52) produces the meson mass,

The ground state mass formula is

and the higher states have the mass formula,

The masses of cc¯ and bb¯ mesons from Eqs. (47) and (48) are reported in Table XXIII and compared with experimental values as well as those from SWE. Again the agreement of masses from our VP with experiment and SWE is remarkable. It is further emphasized that mass values from our VP are closer to experimental findings than those predicted by SWE through MNM. It implies our VP is more accurate than MNM. We also compare our VP masses of S-, D-, F- and G-waves with sophisticated calculations of Refs. in Tables XIX toXXII. The SWE is found not accurate enough in predicting higher quarkonia masses, but Tables XIX toXXIII show our VP works better not only for low lying states but it is very much suitable for spectroscopy of higher quarkonia states as well.

There are some more add ons from our VP. We use Eq. (45) in Eq. (31) and then result in Eq. (37) to get βVPcc¯ and βVPbb¯, which are reported in Table III. We use Eqs. (38), (47) and (48) in Eq. (14) to get vS,VP for different states of heavy quarkonia reported in Table IV. The quantity within parentheses in Eqs. (47) and (48) is the QIE which we report in Tables I andII. So every thing done by SWE numerically is more accurately and efficiently done analytically by our quantum mechanical VP.

5.Conclusion

In this article, we have explored the mass spectra of charmonia and bottomonia in a non-relativistic framework invoking the Song-Lin potential as the effective interaction that binds heavy quarks in these meson systems with a single set of parameters. By numerically solving the resulting SWE through the Numerov strategy, we were able to calculate spin averaged masses of S, P, D, F and G states. For the sake of illustration, we have not included in the discussion the techniques and ideas of spin or tensor interactions. These findings are straightforwardly extended in the Appendix A.

We have calculated S-wave annihilation decay rates without and with momentum width and velocity corrections by using wave functions at contact obtained from MNM. The agreement of theory with experiments after applying the corrections indicate that expressions for these decay rates should be revisited by properly incorporating the quark velocities and momentum widths so that the factor |ψ(r=0)|2 of wave function at contact gets replaced with |ψ(r=0,v,β)|2. This factor may be calculated by identifying the string fragmentation in the Lund Model with decay of a meson and using our VP based momentum Eq. (16) as transverse momentum distribution in the Lund Area Law. The linear rise in scalar part of the Cornell potential is inappropriate for calculation of higher excited states of mesons unless other contributions like the readjustments of the constituent quark mass, the free parameters in SLIPS and the complicated relativistic corrections in heavy quarkonia are added. The remedy to this issue has been sought in the color screening of SLIPS (see Appendix C).

The quark kinetic energies found from Table III using μv2/2 and β2/(2μ) where μ is the reduced mass of two bare quarks, are in good agreement with each other. This establishes the validity of our models of quark- velocities and momentum widths. The EMC effect observed in atomic nuclei for first generation quarks is correctly predicted by our quark velocity model in the mesons for second generation c quark and third generation b quark respectively in cc¯ and bb¯. It would be interesting to calculate these velocities in the B_c meson which should be such that velocity of c quark in B_c is less than that in cc¯ and the velocity of b quark in B_c is more than that in bb¯. This would complete the theoretical testing of EMC effect in heavy quarkonia. From Tables IV andV it is observed that v_F is about (1.3 - 2)v_S for c quark and v_S is approximately v_F for b quark. So we conclude our velocity model is more general than the usual velocity Eq. (13) which becomes a special case of our model, Eq. (14). From the analytical expressions of quark velocity in Eq. (15) and momentum width (16) and the spatial rates of changes (the gradients) their of, the result is: c quark velocity decreases with increasing quark separation but at the same time quark momentum width increases. This is an indication that QIE causes an increase in the constituent mass of quark in the form of so-called quark dressing and, as for the b quark velocity is concerned, it increases with increasing interquark separations but the increase in the b quark momentum width is seen to be more than this, which leads to increase in constituent mass of b quark like the c quark. This is one of the main reasons for which the eigenvalues of SWE with any one fixed value of the constituent quark mass as in Refs. ^3-18 spoils the accuracy of the calculation for mesons in higher excited states. This unpleasant feature is -to some extent- also visible in Table XIX with the MNM, but it is absent in our VP mass spectra. It is the landmark of our VP, and we find it more convenient the use of constituent quark masses in quark model calculations. Coming back to gradients, all of them are observed to decrease such that change in velocity of c quark occurs opposite to the direction in which c quark momentum changes and this happens in the same direction for the case of b quark. This is one of the reasons that QIE is more in cc¯ than in bb¯. Finally, our VP loses its accuracy nowhere in the mass spectrum when we compare its predictions against other techniques as observed from Tables XIX toXXIII. As far as meson wave functions are concerned, these can be obtained by using Virial Theorem as in Ref. ⁷³. This makes our VP a robust dynamical method applicable to any two-body non relativistic bound state held together by a suitable interaction potential V(r).