The Yukawa potential: ground state energy and critical screening

We study the ground state energy and the critical screening parameter of the Yukawa potential in non-relativistic quantum mechanics. After a short review of the existing literature on these quantities, we apply fifth-order perturbation theory to the calculation of the ground state energy, using the exact solutions of the Coulomb potential together with a cutoff on the principal number summations. We also perform a variational calculation of the ground state energy using a Coulomb-like radial wave function and the exact solution of the corresponding minimization condition. For not too large values of the screening parameter, close agreement is found between the perturbative and variational results. For the critical screening parameter, we devise a novel method that permits us to determine it to ten digits. This is the most precise calculation of this quantity to date, and allows us to resolve some discrepancies between previous results.

We study the ground state energy and the critical screening parameter of the Yukawa potential in non-relativistic quantum mechanics. After a short review of the existing literature on these quantities, we apply fifth-order perturbation theory to the calculation of the ground state energy, using the exact solutions of the Coulomb potential together with a cutoff on the principal number summations. We also perform a variational calculation of the ground state energy using a Coulomb-like radial wave function and the exact solution of the corresponding minimization condition. For not too large values of the screening parameter, close agreement is found between the perturbative and variational results. For the critical screening parameter, we devise a novel method that permits us to determine it to ten digits. This is the most precise calculation of this quantity to date, and allows us to resolve some discrepancies between previous results.

Introduction
The Yukawa potential was proposed by Yukawa in 1935 [1] as an effective non-relativistic potential describing the strong interactions between nucleons. It takes the form V (r) = −α e −µr r , (1.1) and thus can be seen as a screened version of the Coulomb potential, with α describing the strength of the interaction and 1/µ its range. The same potential appears under the name of Debye-Hückel potential in plasma physics, where it represents the potential of a charged particle in a weakly nonideal plasma [2], as well as in electrolytes and colloids. In solid state physics it is known as Thomas-Fermi potential, and describes the effects of a charged particle in a sea of conduction electrons.
In quantum mechanics, the physics of this potential depend strongly on the value of the screening parameter µ. While for the Coulomb case µ = 0 there is an infinite number of bound states, for any positive value of µ the screening is sufficient to reduce this number to a finite one [3][4][5], and for µ larger than a certain critical value µ c , bound states cease to exist altogether. This critical value is proportional to [3,[6][7][8][9][10][11] αm: Despite its superficial closeness to the Coulomb potential, the Yukawa one shares hardly any of the exceptional mathematical properties of the former. To this date, for µ = 0 neither the energy eigenvalues, nor the eigenfunctions, nor the critical screening parameter are known in closed form. This combination of physical importance and mathematical intractability makes the Yukawa potential a natural test case for approximation methods in quantum mechanics.
The purpose of the present paper is fourfold. First, in section 2 we will give an overview of the various approximation methods that have been used to date. Our emphasis here is on the ground state energy E 0 (µ) and on the critical screening µ c , since these are the quantities which we are then studying ourselves in the rest of the paper. We will use these literature data to construct a literature average curve E 0 (µ).
Second, in section 3 we will perform a perturbative calculation of the ground state energy, taking the exactly solvable Coulomb Hamiltonian as the unperturbed one. The special properties of the Coulomb case will allow us to push this calculation to an unusual fifth order.
We further improve on this calculation by including in the unperturbed Hamiltonian the contribution from the Yukawa potential that is linear in µ. We also study the dependence of the perturbative calculation for the energy on the cutoff in the principal quantum number of the Coulomb wave functions that becomes necessary starting from at second order.
Third, in section 4 we then compare our perturbative results with a variational calculation, using a trial function of the type of the Coulomb ground state wave function. Although this simple trial function has been used before, to the best of our knowledge the minimization condition, a third-order algebraic equation, was solved only approximately (e.g. in the textbook of [12], see problem 7.14). Here we give its exact solution, and it turns out that, remarkably, it closely matches our result from fifth-order perturbation theory in the whole range of µ except the region close to the critical end-point.
Fourth, in section 5 we will present a novel, but simple, method to determine the critical screening parameter. We obtain for it the value which is the most precise value to date. We compare with previous results given for the screening parameter.
2 Review of the literature A number of general theorems exist on the existence and number of bound states for a given potential. In 1951 Pais and Jost [3] showed that for a 3-dimensional spherical potential such that I = 2m ∞ 0 dr r|V (r)| is finite, bound states must exist for I > 1. One year later, Bargmann [4] proved that the number of bound states, n l , for a given angular momentum quantum number, l, is bounded by (2l + 1)n l < I . (2.1) For the Yukawa potential, this relation (in our units) is (2l + 1)n l < 2m µ α. In particular, no bound state can exist for µ > 2mα. The inequality (2.1) was rederived and further generalized in 1960 by Schwinger [5].
As was mentioned in the introduction, no exact results exist to date for the wave functions and energies of the Yukawa potential. As to approximative calculations, the most widely used method has been the variational principle. In 1962, Harris [6] used trial wave functions constructed from the 1s, 2s and 3s solutions of the Coulomb potential to obtain very good values for the ground state and the first 45 excited energies of the system. In 1990, Garavelli and Oliveira [8] applied the variational method using the 1s Coulomb solution together with a second wave function involving a screening parameter to be determined. In 1993, Gomes et al. [9] devised a two-step procedure where optimized few-parameter trial functions are obtained from an initial linear combination of atomic orbitals (LCAO) with up to 26 basis functions. This allowed them to obtain very precise values for E 0 (µ) and µ c , and also to demonstrate the delocalization of the ground state wave function in the bound-unbound transition, i.e. ψ 0 (r) → 0 for µ → µ c . They were also able to determine the critical exponents for ψ 2 0 (0) and E 0 for this transition. Somewhat less popular in this context has been perturbation theory. The work by Harris already cited [6] was also the first to treat the Yukawa potential as a perturbation of the Coulomb one, although only in first-order perturbation theory. Gönül et al. [13] in 2006 combined the perturbative treatment with an expansion of the exponential factor e −µr . In 1985, Dutt et al. [14] used a scaled Hulthén potential instead of the Coulomb one as the unperturbed Hamiltonian.
As to numerical approximations, in 1970, Rogers et al. [7] solved the Schrödinger equation numerically. For the same purpose, in 2005 Yongyao et al. [11] used Runge-Kutta and Numerov algorithms, as well as Monte Carlo methods.
There have also been less standard approximations to analyze the Yukawa potential.
Garavelli and Oliveira [8] used an iterative process to solve the Schrödinger equation in momentum space. In 2012, Hamzavi et al. [15] used the generalized parametric Nikiforov-Uvarov method for obtaining approximate analytical solutions of the Schrödinger equation, and showed that this works well for µ 0.15mα.
In Fig. 1 we show a plot of E 0 (µ) obtained by averaging over the results given by various authors, based on TABLE I of [8] (we have not included here the results of [9], since they consider only four different values of µ).
The µ c values given by a number of authors are listed in Table 1 of section 5 below.

Perturbative calculation of the ground state energy
Of the many ways of finding approximate solutions to the Schrödinger equation for a system that cannot be solved exactly, probably the most widely used is perturbation theory, where one builds on the known exact solutions of some other, usually simpler system. However, the perturbative expansion becomes quickly cumbersome at higher order, so that most textbooks of quantum mechanics, e.g. Griffiths [12], give the explicit formulas only to second order. Exceptionally, Landau and Lifshitz [16] give them to third order, the fourth order is worked out in an unpublished article [17], and Wikipedia [18] has the expressions for the energy levels to fifth order (for the non-degenerate case). Those expressions, whose correctness we have verified by an independent calculation, are included in appendix A for easy reference. Here we wish to apply them to the ground state energy of the Yukawa Hamiltonian, taking advantage of the fact that the Coulomb case is exactly solvable: We recall that the eigenvalues of H 0 are where n = 1, 2, 3, . . . ; l = 0, 1, 2, · · · , n − 1; m = −l, · · · , l, and a 0 = 1/mα is the Bohr radius. Y m l (θ, φ) are the spherical harmonics, and L k n the associated Laguerre polynomials (our convention for the latter is given in (A.6) in the appendix and corresponds to that of MATHEMATICA). Since both the unperturbed ground state and the perturbation ∆H are spherically symmetric, it is easily seen that the eigenstates with non-vanishing angular momentum, i.e. with l > 0, will not contribute to any order in the perturbative expansion. This greatly simplifies the expansion, and in particular reduces it to the non-degenerate case, so that we can use the formulas for non-degenerate perturbation theory as given in appendix A, restricting them to the spherically symmetric eigenfunctions ψ n ≡ ψ n00 from the beginning. They involve, apart from the energy differences ∆ nm ≡ E n − E m , only the matrix elements V nm ≡ ψ n |∆H|ψ m , which we obtain in closed form in (A.8). From the second order correction onwards the expressions involve infinite sums over the principal quantum number n, which we were unable to do in closed form. However, all these sums converge very rapidly (at least as 1/n 3 ), so that a cutoff could be used on them; using C++, we were able to sum over the first 19 terms for each infinite sum.
In Fig. 2 we show a plot of the results of this calculation for the ground-state energy as a function of µ at various orders of perturbation theory (choosing m = α = 1), together with the literature average curve obtained in the previous section. We observe that perturbation theory works well up to µ 0.5 and breaks down for µ 0.8. As one would expect, the addition of the higher-order terms delays the onset of this breakdown, but not very significantly. It is also interesting to note that, in the range where perturbation theory works, the perturbation series for fixed µ still shows an apparent convergent behavior to fifth order, even though it is known that the perturbation series in quantum mechanics (as well as in quantum field theory) is generically an asymptotic divergent one (see, e.g. [19]). It would be interesting to push this calculation to even higher orders to see the onset of asymptoticity.
As a check on the cutoff that we have used for the principal number summations, let us also show in Fig. 3 the corresponding plot obtained by summing only over the first 10 terms, rather than 19, for all these sums. The plots are indistinguable in the range of µ where perturbation theory works. The expansion in ∆H is effectively an expansion in µ, which suggests that better results might be obtained by using the same perturbation series with a different break-up of the Yukawa Hamiltonian: instead of (3.1), let us try moving the term linear in µ contained in the Yukawa potential, which corresponds to a constant term in the Hamiltonian, from ∆H to H 0 : The new H 0 has the same eigenfunctions as H 0 , and eigenvalues shifted by the constant µα: Redoing the perturbative calculation with this modification, up to fifth order and with the same cutoff at n = 19, we have obtained the results for E 0 (µ) shown in Fig. 4.
Comparing with Fig. 2 we see that the new break-up does not make any improvement. An important remark is that increasing the cut-off above n = 19 leads to numerical instability in the fifth order perturbation theory due to delicate cancellations between large numbers (generated by the combinatoric factorials present in the equations in the appendix). Various ad-hoc modifications to the methods of calculating the summands become necessary to alleviate these problems.

Variational principle
For the case of the ground state energy, the variational principle provides an approximation method that is more universally applicable than perturbation theory, and often yields accurate results with relatively little effort. It states that to obtain the ground state energy of a system described by the Hamiltonian, H, in the case that it is not possible to solve the Schrödinger equation exactly, then one can pick any normalized wave function ψ whatsoever and, using the spectral representation of the Hamiltonian, it is easily seen that one always has E 0 ≤ ψ|H|ψ . (4.1) That is, unless ψ is the true ground state, the expectation value of H in the state ψ is certain to overestimate the ground state energy [12].
For the Yukawa Hamiltonian (1.1), the simplest possible choice is a trial wave function ψ t (r) that mimics the ground state wave function of the Coulomb potential, that is the ψ 100 of equation (3.3): This wave function is normalized, and b is the variational parameter that needs to be adjusted to minimize the expectation value of the Hamiltonian (3.3). Since there is only radial dependence, the calculations are simple and lead to where we have also used the relation α = 1/(ma 0 ).
The minimization of the expression (4.3) leads to an algebraic third-order equation. Curiously, although the variational method has been applied to the Yukawa potential by a number of authors, and the trial function (4.2) was used in [12] (in problem 7.14), the exact solution of this equation seems to be not in the literature. MATHEMATICA gives it as  can be obtained by replacing the exact b 0 (µ) of (4.4) by its small µ approximation (4.5).
Remarkably, the variational and the fifth-order perturbative curves are in close agreement. 5 The critical screening parameter µ c As we mentioned in the introduction, an important difference between the Coulomb and Yukawa cases is that, for the latter, bound states cease to exist if µ becomes larger than a critical value µ c . Such a transition in quantum mechanics holds important information on the dynamics of the system. For example, in solid state physics the existence of bound states makes it possible to condense electrons around protons and get an insulating system, while  GCM [9] (V) GCM [9] (LCAO) Harris [6] (VC) RGH [7] (N) YXK [11] (N) in their absence the system is a metal [20]. The transition between both regimes is called Mott transition. The value of µ c = (·)mα is not known exactly. In Table 1 we give some values for (·) found in the literature, together with the ones following from our approximate calculations above, determined by the intersection of the E 0 (µ) curve with the E 0 = 0 axis.
Further, we will now present a method for the calculation of µ c that is quite simple, but which we have not been able to find in the literature. For µ = µ c we have the lowest eigenvalue E 0 = 0 with corresponding radial Schrödinger equation In terms of g(r) ≡ rψ 0 (r), this can be written as According to (1.2) the critical value µ c is proportional to αm, so that we can as well set αm = 1. Equation (5.2) describes, for µ just below the critical value µ c , a bound state solution, and for µ just above µ c a scattering solution (here we neglect an infinitesimal contribution to the right hand side of (5.2), −2mE 0 g(r), that is present whilst µ = µ c ). In either case we must have g(0) = 0, for regularity of the wave function at the origin, and in the bound state case also g(∞) = 0, since the wave function must decrease (much) faster than 1/r for large r to be square-integrable. This suggests [21] that one could distinguish between µ < µ c and µ > µ c by checking numerically whether (5.2) can be solved with the boundary conditions However, these boundary conditions are not the most convenient ones for numerical purposes. It is thus important to observe that we can as well replace them by the conditions where g 1 is an arbitrary non-zero number, and check whether the numerical solution of (5.2) leads to g(∞) = 0 or not. These procedures are equivalent, since the boundary conditions (5.3) are homogeneous, so that a solution of (5.2) fulfilling them can be rescaled to make g (0) take any given non-zero value g 1 (non-zero, since g (0) = g(0) = 0 leads to the trivial solution g(0) ≡ 0). Moreover, it is easy to check the behaviour of g(r) at large r, since (5.2) implies that g (r) rapidly converges to a constant for large r. Thus, the critical µ c can be found by starting with a µ somewhat below it, and then hiking it up little by little, each time solving (5.2) numerically with some arbitrary g 1 > 0, and checking whether the slope of g(r) still becomes negative for large r, as it must for g(r) to represent a bound state solution (in the bound state case, the slope of g(r) ultimately must go to zero for very large r; in an exact calculation, this would be ensured by the tiny positive contribution to the right hand side of (5.2) that is present whilst µ < µ c ). From these plots we conclude that µ c /αm lies between 1.190612210 and 1.190612211. This is compatible with the LCAO result of [9], µ = 1.19061227(4) (except for their last digit) and [11], but not with [7] and [8]. We have also checked that this result is stable under a further reduction of the radial cutoff r 0 .

Conclusions
In this paper, we have summarized the existing literature on the ground state energy E 0 and the critical screening parameter µ for the Yukawa potential, and we have performed three calculations that, although straightforward, to the best of our knowledge have not been done before: (i) a fifth-order perturbative calculation of E 0 using the exact solutions of the Coulomb potential together with a cutoff on the principal number summations (ii) a variational calculation of E 0 using a simple Coulomb-like radial wave function and the exact solution of the corresponding minimization condition, a third-order equation (iii) a highprecision determination of µ c by numerical integration of the radial Schrödinger equation with appropriate boundary conditions. Our main findings are the close agreement between the fifth-order perturbative result with the variational one, and a calculation of µ c to ten significant numbers. This is one digit more than was obtained in the LCAO calculation of [9], Here we have used the short-hand notation V nm ≡ ψ n |∆H|ψ m and ∆ nm ≡ E n − E m .