2.4 Numerical solutions

There are hundreds of possible methods for solving the radial equation (2.5) for confining potentials, and many new suggestions can be found by reading journals of mathematical physics or computational physics as well as the American Journal of Physics and similar publications [24]. Let us describe briefly two methods: the Numerov method and the Multhopp method.

In the simplest of the direct numerical solutions of the radial equation (2.5), one separates the numerical integration in r at a given energy E, and the search for the eigenvalues. Hartree [25], for instance, gave an efficient method based on the Numerov algorithm, which nowadays can easily be implemented on computers or even on programmable calculators. Let us rewrite Eq. (2.5) as

(2.18)

and search for an outgoing u_out(r) solution starting from u_out(r) ∝ r^l+1 near the origin and an ingoing solution u_in(r) starting from u_in(r) ∝ exp(−∫ dr) at large r. Using the Taylor expansion and the above differential equation, one easily shows that the values u_n of u_out (or of u_in) at equally spaced points r = nh fulfil

(2.19)

leading to a very accurate (and stable) iterative algorithm for computing the u_n’s.

Now the (n + 1)^th energy E _n is determined as follows. A trial energy E is above E_n if the corresponding outgoing solution has at least (n + 1) nodes in ]0,r_c[, where r_c is the classical turning point, for which E = V (r_c). One ends quickly with an interval a_n < E_n < b_n in which E_n is the only possible eigenvalue. In this interval, E_n is the only value for which one can match the incoming and the outgoing solutions as well as their first derivative. One can always impose

(2.20)

If u′_in(r₀) ⁄= u′_out(r₀) for some trial energy , then a better estimate of E_n is given by

(2.21)

Iterative applications of this recipe converge very rapidly towards E_n. Equation (2.21) may be derived from the Wronskian theorem, whose many uses in quantum mechanics are well described in Ref. [21]. Alternatively, one may consider that the mismatching solution {u(r) = u_out for 0 < r < r_o and u(r) = u_in(r) for r_o < r < ∞} corresponds exactly to the eigenvalue for the modified potential

(2.22)

When first-order perturbation theory is applied to V starting from the solution of , one gets exactly the prescription (2.21).

The application of matrix methods to quarkonium calculations was proposed independently by Basdevant et al.[26] and Olsson et al.[27] (and, probably, many others!). One of the many possible variants is the following. The transformation

(2.23)

leads to the new equation

(2.24)

where v(x) is submitted to v(0) = v(1) = 0. The value of r₀ should be chosen to be of the order of the typical radius of the wave function. One now computes v from its Fourier expansion

(2.25)

For finite N, all information is contained in the a_i’s and thus in the values v_j of v(x) at the points x_j = j∕(N + 1), since

(2.26)

One ends with the matrix eigenvalue equation

(2.27)

(2.28)

which is a matrix with

hidden symmetry. Indeed, the simple change of function

(2.29)

which, not surprisingly, corresponds to the simple normalization

(2.30)

leads to the new matrix equation ∑ B_ijw_j = Ew_i where, again, w_j = w(j∕(N + 1)), and the matrix

(2.31)

is now symmetric. The diagonalization provides an approximation of the N first levels. A level of given radial number n is well described as soon as N ≫ n, and the convergence as N →∞ is reasonably achieved, provided the parameter r₀ has been chosen to be of the right order of magnitude. The wave functions are mutually orthogonal, to the extent that the scalar product is replaced by its discretized approximation

(2.32)

The tail of the wave function as r →∞ is, however, not too well reproduced, since the exponential vanishing of u(r) is replaced by a zero of finite order for w(x) at x = 1. The main advantage of this method, for our purpose, is its immediate generalization to coupled equations. This will be very useful for Faddeev, hyperspherical, or Born–Oppenheimer approaches to the three-body problem.

For illustration, we display in Table 2.1 the lowest eigenvalues of the linear potential V (r) = r in the case μ = 1, as a function of the dimension N of the matrix (2.31). The exact eigenvalues are given by the negative of the zeros of the Airy function [28]. The parameter r₀ is chosen as being α₀^−1∕2 where exp(−α ₀r²∕2) is the best variational approximation of Gaussian type to the ground state.