Numerov and phase-integral methods for charmonium

This paper applies the Numerov and phase-integral methods to the stationary Schrödinger equation that studies bound states of charm anti-charm quarks. The former is a numerical method well suited for a matrix form of the second-order ordinary differential equations, and can be applied whenever the stationary states admit a Taylor-series expansion. The latter is an analytic method that provides, in principle, even exact solutions of the stationary Schrödinger equation, and well suited for applying matched asymptotic expansions and higher-order quantization conditions. The Numerov method is found to be always in agreement with the early results of Eichten et al., whereas an original evaluation of the phase-integral quantization condition clarifies under which conditions the previous results in the literature on higher-order terms can be obtained.


Introduction
As Dirac pointed out in his book of quantum mechanics [1], if the exact equations of a physical theory are too difficult, they are of no use, and one has to resort to approximations. In particular, it is extremely difficult to evaluate the effective form of interquark forces from quantum chromodynamics, and hence since the late seventies it became of interest to consider non-relativistic models of charmed quarks. For this purpose, the work in Ref. [2] assumed an instantaneous potential which is a superposition of a linear and a Coulomb term. Since then, many original investigations of bound states for charmonium appeared in the literature, for which we do not even attempt to write a comprehensive list here.
We instead remark that the resulting stationary Schrödinger equation can be studied with the help of advanced numerical and analytic tools, and it has been our aim to test the efficiency of such tools for the non-relativistic Coulomb plus linear potential model. For this purpose, Sect. 2 outlines the Numerov method [3,4], which has found so far a wide range of applications: maximal adaptation to the Schrödinger equation [5], reduction of the number of stages in the algorithm [6], second-order differential equations with oscillating solutions [7], two-point boundary-value problems [8], 7-stage eighth-order methods [9], vibrational eigenstates of linear triatomic molecules [10]. Section 3 studies instead the phase-integral method [11,12], which is an improvement of the JWKB method, better suited for studying the one-directional nature of connection formulas and higher-order quantization conditions. The integrals occurring in the phase-integral quantization condition are studied in Sect. 4, while numerical results are displayed in Sect. 5, concluding remarks are presented in Sect. 6 and relevant details are provided in the "Appendices".

Numerov method
A very simple and powerful numerical approach to solve the one-dimensional (or the radial part of the 3-dimensional) stationary Schrödinger equation is the Numerov method [13][14][15]. The method is useful to integrate second-order differential equations of the general form: with Cauchy data at some point x 0 y(x 0 ) = y 0 , y (x 0 ) = y 0 . (2.2) In the case of the stationary Schrödinger equation as we will discuss in the next section, one has y(x) = ψ(x) s(x) = 0, g(x) = 2m where V (x) is either the potential in one spatial dimension, or an effective potential that includes also the effects of angular momentum in the three-dimensional spatial case. Let us briefly outline the method. If the stationary state ψ(x) admits a Taylor series expansion, one can write that (2.4) which displays the fourth-order nature of the algorithm; hence, we can obtain the second derivative of ψ(x) by considering Eq. (2.5) can be regarded as an equation for the fourth derivative of ψ, giving at If we consider a finite box (i.e., a closed interval) for the variable x, [x 0 , x N ], and a lattice of N points x i evenly spaced by δ, Eq. (2.6) for i ∈ {1, . . . , N − 1} can be written as: As was pointed out in Ref. [16], by defining the N -dimensional column vector ψ = (. . . ψ i−1 , ψ i , ψ i+1 , . . .), Eq. (2.7) reads as: while I p is the matrix of unit entries along the pth diagonal (I 0 is the Identity N × N matrix). The matrixB is invertible and therefore which can be viewed as an eigenvalue problem and easily solved by a simple numerical code, for example, with the help of Mathematica.

The phase-integral method
Both in one-dimensional problems and in the case of central potentials in three-dimensional Euclidean space, the Schrödinger equation for stationary states leads eventually to a second-order ordinary differential equation having the form: where V (z) is either the potential in one spatial dimension, or an effective potential that includes also the effects of angular momentum. The notation z for the independent variable means that one can study Eq. (3.1) in the complex-z plane, restricting attention to real values of z only at a later stage. In the phase-integral method, one looks for two linearly independent, exact solutions of Eq. (3.1) in the form: Since the Wronskian of two linearly independent solutions of Eq. (3.1) is a non-vanishing constant, while the Wronskian of the functions (3.2) is −2i D 2 dw dz , for consistency one finds One can therefore write (up to a multiplicative constant) where w(z) ≡ z q(τ )dτ is said to be the phase-integral, while q(z) is the phase integrand [11,12]. By virtue of Eqs. (3.1)-(3.4), the exact phase integrand q(z) solves the differential equation which is called the q-equation. Suppose now that it is possible to determine a function Q : z → Q(z) that is an approximate solution of the q-Eq. (3.5). This means that χ 0 , defined by and re-expressible in the form: must be much smaller than 1. The work in Ref. [12] proves that the phase integrand q(z) is related to the freely specifiable base function Q(z) by the asymptotic expansion In order to obtain a stationary state that is regular at the origin, all choices of Q fulfilling the condition [11,12] lim z→0 are admissible. In our work, we have exploited precisely this freedom, by writing R(z) and Q 2 (z) in a form compatible with (3.9), i.e., [11,12,17] where, on denoting by a and b the dimensionful parameters in the linear plus Coulomb potential E being the energy levels for bound states. In other words, following Refs. [17,18], we assume that bound states with energy E exist, and that −Q 2 has distinct zeros x 1 and x 2 on the positive half-line, and a negative zero x 0 on the negative half-line. The phase-integral quantization condition is [19] N j=0 L (2 j+1) = s + 1 2 π, s = 0, 1, 2, . . . , (3.14) up to an error term which is studied in detail in Ref. [11], the L (2 j+1) functions being defined by the integral where the contour encloses clockwise the positive roots x 1 and x 2 and has a branch cut from −∞ to x 0 < 0 and another branch cut from x 1 to x 2 . For example, if j = 0, we deal with where P(z) is the polynomial of third degree such that If one defines [17] one can solve for x 0 , x 1 and x 2 according to The coefficients on the second and third line of Eq. (3.17) take therefore the form: Thus, upon changing variable (see "Appendix A" for the notation on Jacobi elliptic functions) according to [17,18] the polynomial P(z) can be re-expressed, after verifying patiently some cancellations, in the form: At this stage, use is made of the identities and one finds eventually because (sn(u)) = cn(u) dn(u). On passing from the z− to the u− integration, the integral (3.16) becomes proportional to the integral from −K (m = k 2 ) to K (m = k 2 ) of an even function of u. Hence, one finds eventually The method developed in Ref. [17], that we here describe in a more detailed way, amounts to imposing the quantization conditions (3.14), starting of course with L (1) : and then solving numerically the resulting transcendental equation for x 2 . On the other hand, since by assumption Q 2 (x 2 ) = 0, one has also from Eq. (3.11) the equation On multiplying both sides of Eq. (3.30) by (x 2 ) 2 , one finds This tells us that the polynomial of third degree has a root at z = x 2 , and we require that x 0 and x 1 should be the remaining two roots. We can therefore exploit Eq. (3.17) and, upon defining while Eq. (3.35) yields and these solutions satisfy identically Eq. (3.34). Moreover, the terms in Eq. (3.18) can be now expressed in the form: In light of Eqs. More precisely, the energy levels obtained from the phase-integral method can be characterized by means of four numbers: the order N of the phase-integral approximation, the dimensionless strength B of the Coulomb potential, the integer s in the quantization condition (3.14) and the angular momentum quantum number l = 0, 1, 2, . . ..

The L (i) integrals
The evaluation of the integral L (1) in Eq. (3.28) yields (see "Appendix A" for the notation on elliptic integrals K , E, , with the understanding that m = k 2 hereafter) where, from Ref. [17], Next, the integral L (3) reads as: where, as is stressed in Ref. [17], u 0 is (see our analysis after Eq. (4.23)) a point in the complex u-plane such that sn(u 0 ) = 0, cn(u 0 ) = 0, dn(u 0 ) = 0. By exploiting the identities (3.25), here re-expressed in the form: the integrand I (3) of Eq. (4.4) can be decomposed in partial fractions with respect to sn 2 (u), i.e., (setting sn 2 (u) = x for simplicity of notation) where F i is a function of m and α 2 , ∀i = 1, 2, 3. By comparison of left-and right-hand side, we find therefore the three equations By insertion of Eq. (4.7) into Eq. (4.8), and subsequent addition of Eqs. (4.8) and (4.9), we find Eventually, we obtain F 3 from Eqs. (4.7)-(4.10) in the form: At this stage, we can exploit the three indefinite integrals [20] where the function E is defined by where while, upon exploiting the identity (we do not write the second argument (i.e., m) of Jacobi elliptic functions for simplicity of notation) we find from Eqs. (4.12)-(4.16) and "Appendix B" that where by repeated application of Eq. (4.5), while having set A crucial remark is now in order. Since u 0 is in general complex, the quantization condition (3.14), written here in the approximate form (cf. Ref. [21]): would lead to complex energy eigenvalues because C(u 0 , m, α 2 ) would be complex-valued, leading in turn to complex values of L (3) . Of course, this is inconsistent. As far as we can see, the only way out lies in looking for the particular values of u 0 such that C(u 0 , m, α 2 ) = 0. where (4.30) We obtain therefore A particular solution can be obtained upon setting m = 1, which implies, from Eqs. (4.28)-(4.30), and hence which is solved by (4.37) Moreover, for generic values of m and α we find from Eqs. (4.31) and (4.32) the two families of solutions which formally can be written as: For these values of u 0 , our expression of L (3) , as already discussed, is the same as the one in Ref. [17], but before performing the numerical calculations we should verify that these points in the complex u-plane are such that sn(u 0 ) = 0, cn(u 0 ) = 0, dn(u 0 ) = 0. If one bears in mind the expressions in Eqs.

Numerical calculations
In this section, we perform a numerical comparison between the solution of the radial part of the three-dimensional Schrödinger equation with a central potential by using the Numerov method (cf. Sect. 2) and the outcomes of the phase-integral method discussed in the previous sections. In order to test our code, in Table 1 we have collected dimensionless energy eigenvalues A (cf. Eq. (3.13)) for the parameters B and occurring in the central potential discussed in Ref. [17] for different mesh values. As one can see, also for a small grid the results are in good agreement with the final results, which are obtained for N = 5000. In all calculations, we use [10 −4 , 50] as the range for the dimensionless z parameter (cf. (3.13)). In Table 2, for the same set of values of B and in Table 1, we write the dimensionless energy eigenvalues A (cf. Eq. (3.13)). The agreement is very good for low values of B but decreases for large B, for which however there is an improvement when the value of increases. It seems to us that this last feature can be understood upon bearing in mind that, at large l, the functions R and Q 2 in Eqs. (3.10) and (3.11) take almost equal values. More precisely, there is not a unique way of fulfilling the limiting condition (3.9). Once R(z) is given as in Eq. (3.10), our choice of Q 2 (z), inspired by the work of Refs. [11,12], is compatible with (3.9), but one might consider as well a choice of Q 2 (z) in the form: where f is any function such that lim z→0 z 2 f (z) = 0, e.g., f (z) = α e −z . This is not a rigorous argument either, but it displays clearly that the arbitrariness in the choice of Q 2 (z) may account for different theoretical estimates of the energy eigenvalues, since the condition of vanishing Q 2 would no longer lead to Eq. (3.30), but rather to the equation if Eq. (5.1) is taken to hold. In other words, different values of B and different choices of f (z) will affect the theoretical estimate of dimensionless eigenvalues A.  We use six significant digits in order to make it easier to perform a comparison between our findings and the results in Ref. [2] A question of crucial importance is whether the Numerov method displays a good convergence rate. For this purpose, following Refs. [22][23][24][25], one can evaluate from the values in Table 1 the quantity where A k is, for each row, the value of A displayed in the k-th column. We find, for example, for showing that the numerical convergence is quite good for the physical cases we are considering. Of course, as shown in Refs. [22][23][24][25], one can also define the quantity where A Ph I is the value of A provided by the phase-integral method. However, in our opinion such a test merits a separate paper, because the optimal estimate of A from the phase-integral method is affected by different values of B and of the function f as we have just discussed after Eq. (5.2). In order to study the effectiveness of the phase-integral method, we should now consider more accurately the effects on A of the perturbative aspect of the method by taking into account the phase-integral quantization condition in Eq. (3.14). As we have already discussed, in the quantization condition we can fix the "perturbative" order by choosing the value of j and varying s in the set of natural numbers {0, 1, 2, . . .}. This analysis was already performed in Ref. [17], where the modulus of the difference between the values of A, | A| = |A N − A Ph I |, obtained at some perturbative order and the result of a numerical approach was plotted as a function of s. We have done the same by considering the fact that our numerical value of A was obtained from the Numerov method -discussed in Sect. 2. The numerical results of the analysis in the previous section can be read in Fig. 1 where we have collected the plots of | A| as a function of s. The plots are for fixed values of B and, for each value of B, the solid, dashed and dashed-dotted curves correspond to = 0, 1, 2, respectively. As we can see, for values of s larger than s = 6 all the curves are smooth; the agreement between phase-integral approach and numerical results is very good for small values of B and it is the best for = 2, and this is true also for B = 5 and B = 10 where the case = 0 and = 1 shows large | A|. In Fig. 2, the plots of | A| as a function of s are gathered together for the same values of B and in Fig. 1, in this case we are adding, in our approximation, the contribution of L (3) . By comparing the figures, it is clear that for j = 1 | A| is almost two orders of magnitude smaller than the case for j = 0 testifying the quality of the approximation.
As discussed at the end of the previous section, we have numerically checked our hypothesis regarding the existence of a value of u 0 such that, for any solution of the quantization relation, C(u 0 , m, α 2 ) = 0 and sn(u 0 ) = 0, cn(u 0 ) = 0, dn(u 0 ) = 0.

Concluding remarks
As far as we can see, the original contributions of our analysis are as follows.
(i) An improved evaluation of the phase-integral quantization condition (3.14) has been obtained, by proving in Eqs. (4.16)- (4.23) that the L (3) integral contains, for the Coulomb plus linear potential, the additional term (4.20) whose occurrence was not discussed in Ref. [17]. However, our original analysis has also the merit of showing that one can explicitly obtain a countable family of complex values of u 0 for which the undesirable additional term in the quantization condition vanishes.
More precisely, the occurrence of C for arbitrary u 0 is unavoidable because the indefinite integrals (4.12)-(4.14) contain the contribution of Jacobi elliptic functions, and hence, the differences in Eq. (4.16) contain the effect of u 0 by virtue of Eqs. (B.1)-(B.10). Our detailed presentation makes it possible to verify all intermediate steps. Once the order at which the quantization condition (3.14) is studied is fixed, the value of u 0 can be obtained, as we have shown in one case at the end of Sect. 4. On going to the next-to-leading order, one obtains a different value of the appropriate u 0 , because one is correcting L (1) by adding L (3) + L (5) , or L (3) + L (5) + L (7) , . . ., and so on. In other words, the function C(u 0 , m, α 2 ) that pertains to L (3) differs from the function of the same arguments that pertains to L (3) + L (5) , and so on. Upon setting it to zero, one obtains an equation of increasing difficulty.
(ii) The Numerov method has been applied successfully to the non-relativistic Coulomb plus linear potential. In particular, the Numerov method has been found to be always in good agreement with the early results in Ref. [2].
We hope that our research will stimulate further work on the ultimate nature of asymptotic expansions [26,27]. It also appears of interest to investigate error terms in the quantization condition (3.14). As is shown in Refs. [11,21], they are of order O(μ), where μ is defined by the integral The numerical evaluation of such an integral is a hard task. Moreover, the experience gained in assessing the phase-integral method might help in studying the QCD potential approximated in Ref. [28]. More precisely, at large values of a parameter N , the static QCD potential V consists of four parts corresponding to 1 r , r 0 , r 1 , r 2 terms, with logarithmic corrections in the 1 r and r 2 terms, and hence one can write, ρ being the dimensionless form of r [28] V (ρ, N ) = V c (ρ) + β(N ) + Cρ + D(ρ, N ) + W, (6.1) where W denotes terms that vanish as N → ∞, while in particular V c (ρ) ∼ π 2ρ log(ρ) as ρ → 0, (6.2) V c (ρ) ∼ − π ρ as ρ → ∞, Funding Open access funding provided by Università degli Studi di Napoli Federico II within the CRUI-CARE Agreement.

Data Availability Statement
This manuscript has associated data in a data repository. [Authors' comment: The datasets generated during the current study are available from the corresponding author on reasonable request.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. dn(u 0 + K ) sn(u 0 + K ) cn(u 0 + K ) = − cn(u 0 ) dn(u 0 ) sn(u 0 ) , (B.9) cn(u 0 + K ) sn(u 0 + K ) dn(u 0 + K ) = − cn(u 0 ) sn(u 0 ) dn(u 0 ) . (B.10)