# References

SINGULARITY IN POTENTIAL, PERTURBATION AND VARIATIONAL METHODS

Yi-Bing Ding<sup>1,2,5</sup>, Xue-Qian Li<sup>1,3</sup> and Peng-Nian Shen <sup>1,4</sup>

1. CCAST (World Laboratory), P.O. Box 8730, Beijing 100080, China

2. Physics Department, Graduate School of Academia Sinica, Beijing, 100039, China

3. Department of Physics, Nankai University, Tianjin, 300071, China

4. Institute of High Energy Physics, Academia Sinica, Beijing, 100039, China

5. Department of Physics, University of Milan , 20133 Milan, Italy

Abstract

In this work, we carefully study the energy eigen-values and splitting of heavy quarkonia as there exist $`1/r^3`$ and $`\delta ^3(\stackrel{}{r})`$ singular terms in the potential which make a direct numerical solution of the Schrödinger equation impossible. We compare the results obtained in terms of perturbation and variational methods with various treatments.

It is well known that the interactions between heavy quarks (or heavy quark-antiquark) can be reasonably described by a combination of the Coulomb-type potential which is induced by the single-gluon exchange and the confinement potential. That form corresponds to the zero-th order potential.

As the first order corrections, the Breit-Fermi potential can be derived directly from the elastic scattering amplitude of the quark and antiquark, where the non-relativistic approximation is employed .

In the derivation, a non-relativistic expansion with respect to $`(|\stackrel{}{p}|/M_Q)^n(n1)`$ powers is carried out and usually, one only keeps terms up to $`|\stackrel{}{p}|^2/M_Q^2`$ and ignore the higher orders. A Fourier transformation would produce $`\delta ^3(\stackrel{}{r})`$ and $`1/r^3`$ terms in the configuration space. The extra potential term which is related to quark spins can be written in a general form as 

$`V_{SD}(r)`$ $`=`$ $`({\displaystyle \frac{\stackrel{}{L}\stackrel{}{s}_1}{2m_1^2}}+{\displaystyle \frac{\stackrel{}{L}\stackrel{}{s}_2}{2m_2^2}})({\displaystyle \frac{1}{r}}{\displaystyle \frac{dS}{dr}}+{\displaystyle \frac{4\alpha _s}{3r^3}})+{\displaystyle \frac{4\alpha _s}{3}}{\displaystyle \frac{1}{m_1m_2}}{\displaystyle \frac{\stackrel{}{L}\stackrel{}{S}}{r^3}}`$ (1)   
$`+`$ $`{\displaystyle \frac{4\alpha _s}{3}}{\displaystyle \frac{2}{3m_1m_2}}\stackrel{}{s}_1\stackrel{}{s}_24\pi \delta (\stackrel{}{r})`$   
$`+`$ $`{\displaystyle \frac{4\alpha _s}{3}}{\displaystyle \frac{1}{m_1m_2}}[3(\stackrel{}{s}_1\widehat{r})(\stackrel{}{s}_2\widehat{r})\stackrel{}{s}_1\stackrel{}{s}_2]{\displaystyle \frac{1}{r^3}},`$   
where $`\stackrel{}{S}\stackrel{}{s}_1+\stackrel{}{s}_2`$, $`\widehat{r}`$ is a unit spatial vector, $`\stackrel{}{L}`$ is the relative orbital angular momentum of the quark-antiquark system.

As noted, even though the Coulomb potential $`1/r`$ is also singular at $`r0`$, it is benign for solving the second-order differential equation. In fact, the Coulomb potential corresponds to the zero-th order of the non-relativistic expansion of the one-gluon-exchange induced potential, and its singular behavior is fully compensated by the small measure at vicinity of $`r0`$. By contraries, the extra $`\delta ^3(\stackrel{}{r})`$ and $`1/r^3`$ come from higher order (at least the first order) of the expansion. They are more singular and make the equation unsolvable. Here the unsolvability means that when one tries to numerically solve the Schrödinger equation with the potential including $`1/r^3`$ and $`\delta ^3(\stackrel{}{r})`$ in expression (1), he would confront infinity, namely the solution blows up due to the singular behavior of the potential at vicinity of $`r0`$, i.e., no physical solution can be obtained. Landau pointed out that as the potential is more singular than $`1/r^2`$, the equation is unsolvable , obviously $`1/r^3`$ and $`\delta ^3(\stackrel{}{r})`$ are more singular than $`1/r^2`$.

Analysis tells us that such high order singularities are not physical, but brought up by artificial truncation of the non-relativistic expansion, namely one only keeps to $`|\stackrel{}{p}|^2/M_Q^2`$ (or even higher, but finite orders).

Gupta et al. suggested to take an alternative expansion with respect to $`|\stackrel{}{p}|/E`$ instead of $`|\stackrel{}{p}|/M_Q`$ . Obviously the divergence for $`|\stackrel{}{p}|\mathrm{}`$ (or $`r0`$) is avoided, because $`|\stackrel{}{p}|/E1`$ as $`|\stackrel{}{p}|\mathrm{}`$. However, $`E`$ is an operator itself, while it exists at denominator of the expression, the equation is not solvable either unless one takes some extra approximations .

As discussed above, the singularity is artificial because of the non-relativistic approximation, so that in principle, if one deals with the problem in a reasonable scheme, for example, the relativistic B-S equation, such singularity would not appear at all. However, the potential form which is non-relativistic, is most useful and convenient. Many good zero-th order results have been achieved , we have all reasons to believe that the non-relativistic Schrödinger equation properly describes the physics. When the first order correction related to quark spins are considered, which determines the fine-structure and hyperfine structure of the quarkonium, the singularity problem emerges, and it turns the whole equation unsolvable. Therefore, as long as we are working with the non-relativistic Schrödinger system, we are forced to face this unphysical singularity problem. Similar singularity problem was discussed in previous literatures. For example, Frank, Land and Spector analyzed the singular potential , however, they mainly concerned the scattering states instead of the bound states. In this work, we only focus on the bound state problems. Our goal is to correctly solve the problem and obtain the physical quantities, even though there is the singularity.

So now the question is can we obtain the eigen energies and wavefunctions which are sufficiently close to physical reality? Of course, the closeness or validity degree can be tested by comparing the obtained results about the splitting, for example, with experimental data.

The effect of the $`\delta ^3(\stackrel{}{r})`$ term on the s-wave has been carefully studied by some authors , so readers are recommended to the original work for details.

In this work, we only concentrate ourselves at $`1/r^3`$ singular term which does not affect the s-wave (see below).

In this paper, we will obtain the eigen energies in both perturbative and variational methods and then compare the results, we also discuss some alternative ways to modify the potential form for avoiding the singularity which have been employed in literatures. Our numerical results show that there exist certain regions for those possible parameters in the approaches where all the results are consistent to each other with sufficient accuracy.   

(i). Perturbative method.

The zero-th order potential does not include the singular parts of (1), so that we can easily solve the corresponding Schrödinger equation and obtain the zero-th order eigen energies and wavefunctions. Then we are going on studying the first order equation which involves the singular terms and obtain corresponding energies in terms of perturbative method.

Let us investigate the singular terms related to $`1/r^3`$ in eq.(1), we find that $`\stackrel{}{L}\stackrel{}{S}=\frac{1}{2}[(\stackrel{}{L}+\stackrel{}{S})^2\stackrel{}{L}^2\stackrel{}{S}^2]`$, so for the S-state ($`l=0`$), its eigenvalue is zero, moreover, the tensor-coupling term in (1) is proportional to an operator

$$V_{tensor}=[3(\stackrel{}{s}_1\widehat{r})(\stackrel{}{s}_2\widehat{r})(\stackrel{}{s}_1\stackrel{}{s}_2)],$$

which is sandwiched between two states which are solutions of the zero-th order equation. If both sides are S-states, then

$$𝑑\mathrm{\Omega }<Y_{00}|V_{tensor}|Y_{00}>=0,$$

where $`\mathrm{\Omega }`$ is the solid angle. The lowest non-zero value comes from P-state, $`<l=1|V_{tensor}|l=1>0`$, and the operator can also result in a small mixing between S-state and D-state, i.e. $`<l=0|V_{tensor}|l=2>0`$. We are not going to discuss the details of the mixing in this work.

Therefore, the singular $`1/r^3`$ terms do not bother the S-states, but only affect the higher orbital-angular momentum states ($`l1`$), this enables us to avoid singularities in practical calculations (see below). With the first order correction the Schrödinger equation for $`l=1,S=1`$ and $`J=2,1,0`$ ($`\stackrel{}{S}=\stackrel{}{s}_1+\stackrel{}{s}_2`$) is recast as

$$[\frac{1}{2\mu }\frac{1}{r^2}\frac{d}{dr}(r^2\frac{d}{dr})+\frac{2}{2\mu r^2}\frac{4\alpha _s}{3r}+\frac{3}{2m^2}(\frac{4\alpha _s}{3})\frac{1}{r^3}\left(\begin{array}{c}1\\ 1\\ 2\end{array}\right)]R(r)=ER(r),$$
(2)   
where we only keep the $`\stackrel{}{L}\stackrel{}{S}`$ term (dropped out the tensor coupling for our purpose of this study) and take $`m_1=m_2=mM_Q`$, the reduced mass $`\mu =m/2`$.

In this study, our purpose is to investigate the effects of the singular potential, so that in the Schrödinger equation we deliberately omit the confinements piece, such $`\kappa r`$ which is well behaved at $`r0`$ for convenience. This treatment does not change the results of the study, but greatly simplifies the calculation.   

(ii) The practical calculation in terms of perturbation.

For heavy quarkonia $`c\overline{c}`$, $`b\overline{b}`$, the coefficients of the $`1/r^3`$ is much smaller than that of $`1/r`$ which has a milder behavior at $`r0`$, so one can employ the perturbation method to calculate the splitting. Then the zero-th order equation is

$$(\frac{1}{2\mu }\frac{1}{r^2}\frac{d}{dr}(r^2\frac{d}{dr})+\frac{1}{\mu r^2}\frac{4\alpha _s}{3r})R^{(0)}(r)=E^{(0)}R^{(0)}(r).$$
(3)   
For the 1p state, there is an analytical solution as

$`R^{(0)}(r)`$ $`=`$ $`{\displaystyle \frac{r}{2\sqrt{6}a^{5/2}}}e^{r/a}`$ (4)   
$`E^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{2a}}({\displaystyle \frac{4\alpha _s}{3}}){\displaystyle \frac{1}{4}},`$ (5)   
and the Bohr radius is

$$a=\frac{1}{\mu (\frac{4\alpha _s}{3})}.$$
(6)   
If we take the conventional parameters $`m_c=1.84`$ GeV, $`\alpha _s=0.39`$, then

$$a=2.09030\mathrm{GeV}^1,E^{(0)}=0.031096\mathrm{GeV}.$$

We can evaluate the contribution of the $`1/r^3`$ term by means of perturbation. The perturbative hamiltonian is

$$\delta H=\frac{3}{2m^2}(\frac{4\alpha _s}{3})\frac{1}{r^3}.$$
(7)   
Using the zero-th order wavefunction (4), we have

$`I_P`$ $`=`$ $`<P|\delta H|P>{\displaystyle R_P^{(0)}(r)\frac{3}{2m^2}(\frac{4\alpha _s}{3})\frac{1}{r^3}R_P^{(0)}(r)r^2𝑑r}`$ (8)   
$`=`$ $`{\displaystyle \frac{\alpha _s}{12a^3m^2}}.`$   
It is worth of noting, $`r`$ exists at numerator of the $`R^{(0)}(r)`$ expression in (4), so that the integration converges.

Numerically, we obtain

$$I_P=0.00105105\mathrm{GeV}.$$
(9)   
For $`2^3P_2,\mathrm{\hspace{0.33em}2}^3P_1`$ and $`2^3P_0`$ states, the energy shifts are $`I_P,I_P`$ and $`2I_P`$ respectively. These values are only 3$``$ 6% of $`E^{(0)}`$, therefore the perturbation is legitimate and one can be convinced that the results are reliable.   

(iii) Variational method.

The Schrödinger equation is not solvable if there is a potential term $`r^s`$ with $`s>2`$  or $`\delta ^3(\stackrel{}{r})`$, it is because the corresponding eigen-energy tends to $`\mathrm{}`$ and average radius turns to 0. ”Particle falls to the point $`r=0`$”. Let us focus on the $`1/r^3`$ singularity. Definitely, the $`\mathrm{}`$ solution is not physical, but a virtual solution of the equation with the singular potential. The key point is how to systematically get rid of the virtual unphysical solution and achieve the physical one.

The hamiltonian is

$$h_P=\frac{1}{2\mu }\frac{1}{r^2}\frac{d}{dr}(r^2\frac{d}{dr})+\frac{1}{\mu r^2}\frac{4\alpha _s}{3r}\frac{3}{2m^2}(\frac{4\alpha _s}{3})\frac{1}{r^3},$$
(10)   
where we take $`J=1`$ as an example. If we employ the P-state wavefunction (4) as a trial function, the average value is

$$<h_P>=\frac{\alpha _s}{3a}+\frac{1}{8\mu a^2}\frac{\alpha _s}{12m^2a^3},$$
(11)   
where $`a`$ is a variational parameter in the expression. Later we need to differentiate the expression with respect to $`a`$ and obtain the minimum.

For a $`c\overline{c}`$ bound state, the relation between $`<h_P>`$ and $`a`$ shows that as $`a0`$, $`<h_P>\mathrm{}`$, there is no lower bound. The corresponding average radius is

$$<r>=5a,$$
(12)   
so that as $`a0`$, $`<r>0`$. It indicates that if one took the minimal value $`\mathrm{}`$ as the binding energy, the radius turns out to be infinitesimal, it cannot exist in principle. It is exactly the Landau’s conclusion. We can also prove that the conclusion is independent of the adopted form of the trial functions. This result indicates that the singular behavior of the potential near $`r0`$ is transferred to $`<h_P>`$ at $`a0`$.

Strictly speaking, the equation is of no solution in this case, does it mean that the Schrödinger equation cannot describe our physics picture? As pointed out in previous section, this singularity is caused by the ill behavior of the non-relativistic approximation near $`r=0`$. One can conjecture that the given Schrödinger equation can still correctly describe the physics even with the unphysical singular terms, but we have to solve it in a proper way. Namely since the singularity at $`r0`$ is artificial and unphysical, the boundless energy corresponding to $`a=0`$ does not reflect physical reality, i.e. the infinite $`<h_P>`$ at $`a=0`$ is not a physical solution and must be eliminated.

From another angle, as $`a0`$, the trial function (4) has a limit as

$$\underset{a0}{lim}\frac{r}{a^{5/2}}e^{r/a}=0.$$
(13)   
The trial function tends to zero, the solution becomes trivial and does not make any sense at all.

In fact, if one differentiate $`<h_p>`$ with respect to $`a`$, besides the minimal value $`\mathrm{}`$ at the boundary $`a=0`$, $`<h_P>`$ has a local minimum at $`a_1`$ and it is

$$a_1=\frac{12+\sqrt{144192\alpha _s^2}}{32\alpha _s}\stackrel{\alpha _s=0.39}{=}1.97832\mathrm{GeV}^1.$$
(14)   
The corresponding average energy value is

$$E_{var}=0.0322362\mathrm{GeV}.$$
(15)   
When we compare this value $`E_{var}`$ with the result obtained in terms of perturbation

$$E_{per}=E^{(0)}I_P=0.0321470\mathrm{GeV},$$
(16)   
the deviation between $`E_{var}`$ and $`E_{per}`$ is only about $`10^3`$.

For a real meson, $`<r>`$ cannot be zero, so the variational parameter $`a`$ must be non-zero. With this constraint, the boundless energy which corresponds to $`a0`$, is ruled out and we can enjoy the reasonable solution determined by the local minimum. Actually, this local minimum should correspond to the real physical value.

Even though we know the negative infinite minimal value of $`E`$ is not physical and should be abandoned, in practical calculations on computer, it is still bothering, therefore some plausible ways which may eliminate the singularity from the concerned problem would be welcome.   

(iv) Variational method and a ”smear” form factor.

One may introduce a form factor to ”smear” the singularity . The most widely adopted form factor is the error function $`erf(\lambda r)`$. It turns the spatial part of the spin-orbit coupling to

$$V_P^{(SM)}=\frac{3}{2m^2}\frac{4\alpha _s}{3}erf(\lambda r)[\frac{1}{r^3}\frac{2}{\sqrt{\pi }}\frac{1}{r}].$$
(17)   
For any finite $`\lambda `$value, $`V_P^{(SM)}`$ is no longer singular as

$$\underset{r0}{lim}V_P^{(SM)}=\frac{8}{3m^2}\frac{\alpha _s\lambda ^3}{\sqrt{\pi }},$$
(18)   
which is finite as long as $`\lambda `$ is finite. In the approach, the hamiltonian turns into a form

$$h_P^{(SM)}=\frac{1}{2\mu }\frac{1}{r^2}\frac{d}{dr}(r^2\frac{d}{dr})+\frac{1}{\mu r^2}\frac{4\alpha _s}{3r}V_P^{(SM)}.$$
(19)   
Because it has a nice behavior at $`r0`$, in $`<h_p>`$ there is no more the troublesome negative infinity at $`a0`$, and one can safely handle it with the variational method. For $`\lambda =5,10,15\mathrm{GeV}`$, we obtain the minimum of $`<h_P^{(SM)}>`$ which is also the minimal value for the whole range of $`0<a<\mathrm{}`$, and tabulate the results in table. 1.   

Table 1.

| $`\lambda (\mathrm{GeV})`$ | 5 | 10 | 15 |
| --- | --- | --- | --- |
| $`a(\mathrm{GeV}^1)`$ | 1.97967 | 1.97868 | 1.97849 |
| $`E`$ (GeV) | -0.0322276 | -0.322339 | -0.0322352 |

For $`\lambda =15\mathrm{GeV}`$, the result is exactly the same as $`E_{var}`$ achieved in last section. But as $`\lambda >15\mathrm{GeV}`$, $`<h_P^{(SM)}>`$ has a boundless negative minimal value for $`a0`$ again. It means that as $`\lambda `$ in the form factor is greater than $`15\mathrm{GeV}`$, the function cannot smear out the singularity behavior at $`r0`$. The meaning of introducing the ”smear” form factor is that once we introduce a ”smear form factor whose parameter resides within a certain range”, we do not need to worry about the boundless minimal value at all, namely then the local minimum is the minimal value and finite, so that we can safely use the variational method to obtain energies. But beyond the parameter range, the boundless minimal value appears again, then we need to look for the local minimum as the real solution, exactly in analog to the treatment of previous section.   

(v) The variational method and the ”regularization”.

(a) Another approach to handle the singularity is the so-called ”regularization” . In this approach, one can ”regularize” $`\frac{1}{r^n}`$ to $`\frac{1}{(r+c)^n}`$ ($`c>0`$), which obviously makes the potential converge at $`r0`$. Thus $`V_P`$ is replaced by

$$V_P^{(r)}=\frac{3}{2m^2}(\frac{4\alpha _s}{3})\frac{1}{(r+c)^3},$$
(20)   
and $`V_P^{(r)}V_P`$ as $`c0`$. We tabulate the results corresponding to various $`c`$values in Table.2.   

Table 2

| $`c(\mathrm{GeV}^1)`$ | 0.1 | 0.05 | 0.035 |
| --- | --- | --- | --- |
| $`a(\mathrm{GeV}^1)`$ | 1.99094 | 1.98525 | 1.98333 |
| $`E(\mathrm{GeV})`$ | -0.0321330 | -0.0321808 | -0.0321963 |

It is noted that the smaller $`c`$ is, the closer to $`E_{var}`$ and values in Table.1 the $`E`$values obtained by the ”regularization” scheme are. But as $`c<0.03\mathrm{GeV}^1`$, the singular behavior of the potential shows up again, namely then $`<h_P^{(r)}>`$ would have a boundless negative minimal value at $`a0`$.

In fact, in sections (iv) and (v), two different approaches are introduced to remedy the singular behavior of the potential at $`r=0`$. Besides the variational parameter $`a`$, there is an extra parameter $`\lambda `$ or $`c`$ being introduced. We notice that if the concerned parameter falls in a certain range, the negative infinite minimal value of $`<h_P>`$ which corresponds to an extreme variational parameter (usually $`a0`$) does not exist at all and a local minimum would give the physical solution. However, once the parameter is beyond the allowed range, such negative infinite $`<h_P>`$ would appear again, in this case, one has to enforce a constraint condition as in section (iii) for the pure variational method, i.e. prior rule out the negative unphysical infinity and keep the local minimum as reasonable solution.

(b) As the simplest ”regularization” scheme, one can replace $`V_P^{(r)}`$ in eq.(20) by

$$V_P^{(r1)}=\theta (rr_0)\frac{3}{2m^2}\frac{4\alpha _s}{3}\frac{1}{r^3}.$$
(21)   
Obviously, it is an equivalent way to eq.(20) to get rid of the singularity at $`r0`$ and $`r_0`$ can be seen as an arbitrarily introduced ”cut-off” in the configuration space or a corresponding $`\mathrm{\Lambda }_01/r_0`$ in momentum space. For the zeroth order $`1^3P_1`$ state the correction of $`V_P^{(r1)}`$ can be calculated and the numerical results are shown in Table.3.   

table 3.

| $`r_0`$ (GeV<sup>-1</sup>) | 0.20 | 0.15 | 0.10 | 0.08 | 0.068 |
| --- | --- | --- | --- | --- | --- |
| $`a`$ (GeV<sup>-1</sup>) | 1.97926 | 1.97886 | 1.97857 | 1.97848 | 1.97844 |
| $`E`$ (GeV) | $`0.0322303`$ | $`0.0322328`$ | $`0.0322347`$ | $`0.0322352`$ | $`0.0322355`$ |

To make more sense, it would be interesting to re-calculate the energy corrections owing to the potential $`V_P^{r1}`$ for the $`1^3P_1`$ state in perturbation and compare

$$E^{r1}<V_P^{r1}>=_{r_0}^{\mathrm{}}𝑑r[\frac{3}{2m^2}\frac{4\alpha _s}{3}\frac{1}{r^3}(R^{(0)})^2r^2],$$
(22)   
with

$$\delta <V_P^{r1}>=_0^{r_0}𝑑r[\frac{3}{2m^2}\frac{4\alpha _s}{3}\frac{1}{r^3}(R^{(0)})^2r^2].$$
(23)   
Now let us tabulate the results in Table 4 as   

Table 4.

| $`r_0`$ (GeV<sup>-1</sup>) | 0.2 | 0.15 | 0.1 | 0.08 | 0.068 |
| --- | --- | --- | --- | --- | --- |
| $`E^{r1}`$ (GeV) | $`0.0010465`$ | $`0.0010485`$ | $`0.0010499`$ | $`0.0010503`$ | $`0.0010505`$ |
| $`\delta `$ (GeV) | $`4.51\times 10^6`$ | $`2.58\times 10^6`$ | $`1.17\times 10^6`$ | $`7.5\times 10^7`$ | $`5.4\times 10^7`$ |
| R | 0.0043 | 0.0024 | 0.0011 | 0.0007 | 0.0005 |

In the table, the ratio R is defined as

$$R=\frac{\delta }{E^{r1}}.$$

The integration over r from 0 to $`r_0`$ denotes the part we ignore when taking the ”regularization” (21) with the step function $`\theta (rr_0)`$. The ratio R reflects the relative errors, and from Table 4, we can see the errors brought by the ”regularization” are about a few thousandths at most. Therefore, according to the experimental accuracy, the scheme is very satisfactory. Indeed, our results restrict the probable error ranges for using the ”regularization” or ”smear” scheme to remedy the singularity.   

(vi) With the Gupta’s approach.

In Gupta’s approach , where an approximation $`\stackrel{}{s}\stackrel{}{p}^{}+\stackrel{}{p}`$ being small is taken, $`\stackrel{}{p},\stackrel{}{p}^{}`$ are the three-momenta of the two quarks. $`V_P`$ is replaced by

$$V_P^g=\frac{4\alpha _s}{3}\frac{3}{2}\frac{f_1(r)}{r},$$
(24)   
where

$$f_1(r)=[1(1+2mr)e^{2mr}]/m^2r^2,$$
(25)   
obviously

$$\underset{r0}{lim}f_1(r)=2.$$
(26)   
Thus this modified $`V_P^{(g)}`$ only possesses the singular behavior of $`1/r`$, so that would not bring up any difficulties for solving the Schrödinger equation by variational method. With this approach, we employ the same trial function (4) and variational parameter $`a`$ to obtain energy as

$$a_{var}^{(g)}=1.98570\mathrm{GeV}^1,\mathrm{and}E^{(g)}=0.0321867\mathrm{GeV}.$$
(27)   
Instead, if we use the perturbation method, we obtain

$$E_{pert}^{(g)}=0.0321085\mathrm{GeV}.$$
(28)   
This is consistent with the values given in Table.1 and 2, as the errors are negligibly small.   

(vii) Conclusion and discussion.

Based on the above calculations, we can draw our conclusion as

(a) As Landau pointed out that an ill-extrapolation of the non-relativistic expansion brings up an artificial singular $`r^3`$ which makes the Schrödinger equation unsolvable. Here the ”unsolvable” means that it is impossible to obtain an exact solution of the equation, so that one is forced to look for a reasonable way to find the physical solution from the given hamiltonian with such singular terms. Our purpose is to seek for a solution which would describe the concerned physics (eigen-energies and eigen-wavefunctions) to a satisfactory degree. The success of the potential model for charmonia encourages this approach and stimulates further studies.

(b) In particular, the variational method without any special treatments for the singularity would result in an infinite negative energy which corresponds to $`<h_P>=\mathrm{},a=0,<r>=0`$. Obviously it is not physical, because it manifests itself as an infinitesimal radius meson, i.e. a point-like particle which would cause singularity as well known. Moreover, as the variational parameter turns to zero, the trial function (4) is zero, so it becomes trivial and unphysical, so must be ruled out.

The local minimum which corresponds to non-zero variational parameter and average radius indeed gives rise to reasonable solution and should be physical. Our results show a satisfactory consistency with that by other approaches, as expected.

(c) The ”smear” and ”regularization” approaches all introduce new parameters whose physical meaning is to decrease the measure near vicinity of $`r0`$ and/or simply eliminate a small neighborhood of $`r=0`$. Such schemes (”smear” or ”regularization”) with certain parameters partly compensate the artificial truncation of the non-relativistic expansion. Our results show that for certain parameter ranges, the obtained results are consistent with the perturbation. For future application, one needs to adopt reasonable values for the parameters, concretely, $`\lambda =15\mathrm{GeV}`$ for the ”smear” approach and $`c=0.035\mathrm{GeV}^1`$ and $`r_0=0.068`$ for the ”regularization”.

Beyond the parameter ranges, the singular behavior appears again.

It may only concern the automatic calculation procedure on computer for obtaining the energies by variational method. Without carefully programming, the computation may overflow once it gets to the artificial unphysical $`a0`$ and refuses to look for the physical local minimum. Therefore when one writes his program, he has to be careful about the possible negative infinity solution.

But if we take the suggested schemes with the concerned parameters residing in the allowed regions, we do not need to worry about the unphysical infinity at all while programming.

(d) The wavefunction problem is not solved yet. One can use the perturbation method to obtain the first order wavefunction, but it needs to sum over contributions from all possible zeroth order states, so practical applications are restricted. With the variational method one can obtain wavefunctions easily, but as we studied in another work , it is not a simple job. Even without the troublesome $`1/r^3`$ terms, it is difficult to obtain wavefunctions with satisfactory accuracy. Corresponding to the well-known potential forms, the Cornell, logarithmic etc. we have found that it is easy to obtain very accurate energy levels, but the errors for wavefunctions, in particular, for their values at origin are rather large unless we carefully choose the forms of the trial functions and use the multi-variational parameters. Moreover, the form of trial functions and number of variational parameters also depend on the form of potential. In  we only succeeded to select proper trial functions and number of parameters for each well-known potential forms. However, that is only for the zeroth order situation, so when we consider the first order correction in the potential, the situation is even worsened.

So how to determine the wavefunctions, especially the values at origin, with existence of $`1/r^3`$ terms, is an open question and we will investigate it in our future works.   

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (NNSFC).