Bound states and energy eigenvalues of a radial screened Coulomb potential

We analyze bound states and other properties of solutions of a radial Schrödinger equation with a new screened Coulomb potential. In particular, we employ hypervirial relations to obtain eigen-energies for a Hydrogen atom with this potential. Additionally, we appeal to a sharp estimate for a modified Bessel function to estimate the ground state energy of such a system. Finally, when the angular quantum number ℓ ≠ 0, we obtain evidence for a critical screening parameter, above which bound states cease to exist.


Introduction
The Yukawa potential is a non-relativistic potential describing strong interaction between nucleons, and the Woods-Saxon potential is a mean field potential for nucleons. Recently, in [1], another such radial potential was introduced: within the context of Time Dependent Density Functional Theory. Indeed, this potential was introduced in order to obtain a Runge-Gross uniqueness theorem, which says that the external potential in a many-body system is completely determined by the one particle density [2]. However, it is still a debate on the applicability of this theorem to Coulombic systems [3], and for this reason a smoothed out Coulombic potential was sought. See figure 1 for a plot of the potential (1.1) with various values of the screening parameter C. Exponential screening is particularly relevant in atomic, nuclear, and solid state physics, and can model the potential between an ionized impurity and electron in a metal [4] or semiconductor [5]; see also [6]. Basic properties of this potential are discussed in [1], but let us mention that the potential belongs to the class of so-called dilatation analytic potentials [7]. This means, in particular, that the spectrum of the associated Hamiltonian can be completely characterized. This is particularly useful in understanding so-called quasi-static energies [8,9] and potential energy surfaces [10]. The dilatation analytic potentials also contain the classical Coulomb potential as well as the Yukawa potential V Y .
There are some subtle differences between the potential V P given by (1.1) and the Yukawa potential, especially related to their behavior at the origin. In particular, all derivatives of V P have a removable discontinuity at the origin, but the derivatives of the Yukawa potential are not bounded. As further motivation for a different screened potential, it was seen in [1] that the classical Runge-Gross method does not apply with the Yukawa potential, but will apply for V P . Essentially, the applicability of this method reduces to saying that the potential should preserve the domain of a particular Hamiltonian (see [1], theorem 2.2).
The purpose of this paper is to further study properties of solutions, including the number of bound states, of the Schrödinger equation with potential given by (1.1). Such problems continue to be of interest, especially if explicit or almost explicit solutions can be obtained [11][12][13]. More precisely, in section 2, we use the Hypervirial and Hellmann-Feynman theorems to calculate the energy levels for a Hydrogen atom with potential V P . This method is based on expanding the expectation values 〈r j 〉 for any integer j into power series in the parameter C. It turns out to be necessary to calculate coefficients of this series for negative values of j; this is related to the socalled Kramers-Pasternack inverse relation [14].
In section 3, we appeal to the Hypervirial theorem and the variational principle to obtain an upper bound for the ground state energy of a system with potential V P . It does not appear to be possible to minimize the Hamiltonian directly as it involves a modified Bessel function of the second kind; see (3.22). However, using a sharp estimate for the modified Bessel function from [15], we can provide an upper bound for the ground state energy; see (3.24). Additionally, by using a small argument approximation for the modified Bessel function of the second kind, we obtain an estimate for the ground state energy assuming certain parameters are sufficiently small; see (3.30).
Next, in section 4, we move to the analysis of bound states. For a given radial potential, there are quite a few upper and lower bounds on the number of bound states; see [16] for a review of some such bounds. These bounds can be translated into conditions on the screening parameter; this is indeed the case for the Yukawa potential [17]. However, these bounds (including the famous Bargmann bound [18]) do not apply with V P as these require the finiteness of a certain integral, which for V P does not converge. We discuss the ℓ = 0 and ℓ ≠ 0 bound states, where ℓ denotes the azimuthal quantum number. In particular, for certain values of ℓ and C we see that an effective potential which contains V P can become repulsive. As such, there appears to be a critical screening parameter for which bound states cease to exist at some point; this is what is seen with the Yukawa potential [19].
Finally, in section 5 we apply an extended Nikiforov-Uvarov method to calculate the eigen-energies E n,ℓ for the radial Schrödinger equation (5.39) with potential V P . Using this method we obtain a relationship between the energies E and the perturbation parameter C; see equations (5.49)-(5.50) for the ℓ = 0 case, and equation (5.57) for the ℓ ≠ 0 case.

Energy levels via hypervirial relations
By using the Hypervirial theorem and Hellmann-Feynman Theorem, we can calculate the energy levels for a Hydrogen atom with screened Coulomb potential V P via power series expansions in the parameter C from the potential V P . This method does not require the calculation of any perturbed wavefunctions [20]. In what follows, 〈 · 〉 denotes expectation value.

Consider the Hamiltonian
From [21], equation (4), the hypervirial relations between the energy E and the expectation values 〈r j 〉 for j 0 are given by In fact, we consider the hypervirial relations with j 0 as well by considering a general hypervirial operator of the form ( )  [20]. Notice that (2.3) is obtained by taking f (r) = r j for j 0; here we take f (r) = r j for any integer j.
We assume that the potential V P can be expanded into a series in terms of the parameter C: has an essential singularity at the origin, we interpret (2.5) as a Laurent series in complex r; in this case, the V k (r) can be found to be Then, using (2.2) and (2.5), we compute Now, we find the right hand side of (2.3) to be Writing the right hand side of (2.7) as one sum yields That is, Dividing by 2( j + 1) yields We assume that the energy and expectation values of r j can be expanded into series in the parameter C as well: where we take the energy of the unperturbed series to be known: From the normalization condition that 〈r 0 〉 = 1 we get that where δ 0,k is the Kronecker delta. Combining equations (2.8), (2.9) and (2.10) provides a relationship between the coefficients E ( k) and the ( ) a j k . To obtain a clearer relation, we use the Hellmann-Feynman Theorem, which says that for a Hamiltonian depending on a parameter C: From (2.9), (2.10), and (2.11), we see that Equating the like terms in the first equality above with those terms in the last equality yields the following expression for the perturbed energies E ( k) :  (15) in [21], we see that the potential V P , as opposed to the Yukawa potential, leads to the ( ) E k in terms of the coefficients ( ) a j k with j negative integers. This is unsurprising since we have expanded V P in terms of a Laurent series with infinitely many negative terms, as opposed to the Yukawa potential which can be expanded as a standard power series. Now, to obtain the coefficients ( ) a j k , we equate the coefficients of C 0 in (2.8): In particular: It is straight forward to continue in this way and find the ( ) a j 0 for j > 0. For example, we find and this can easily be continued. However, in order to find the E ( k) ʼs from (2.12), we need to calculate the ( ) a j k ʼs for j < 0.
As is common in textbooks (see e.g. [22][23][24]), it is assumed that 〈r −2 〉 can be determined from an independent calculation (e.g. by integrating Laguerre polynomials). However, this is not always straight forward. Fortunately, the recent result [14] has produced such an independent calculation, namely that ( ) ℓ á ñ = + r a n where a 0 denotes the Bohr radius. From (2.13) we then obtain Note that all of these coefficients can be written in terms of ( ) a -2 0 , which is given here by (2.16). From (2.12) we obtain ( ) E n 1 immediately: Now we equate the coefficients of C 1 on both sides of (2.8): This allows the calculation of ( ) E n 2 from (2.12) and (2.16): Continuing in this way we can calculate the remaining ( ) a j k for j < 0 as well as the ( ) E n k .

Variational Principle
For the ground state, the variational principle provides a more general approximation method than perturbation theory. The goal of this section is to understand what the variational principle implies for a Hamiltonian with potential V P . We approach this problem in two ways: first, we derive an upper bound for the ground state energy using a sharp estimate for a modified Bessel function. The second approach involves approximating the modified Bessel functions with small arguments in order to obtain an approximate ground state energy estimate.

Sharp estimate from K 2
Here, we work in units such that e = ò 0 = ÿ = m = 1. As in [17] with the Yukawa potential, we take the trial function where b is the variation parameter to be adjusted in order to minimize the expectation value of the Hamiltonian. Notice that (3.21) mimics the ground state wavefunction for the classical Coulomb potential. With this trial wavefunction we compute which is the same as for Hydrogen. This follows directly from the Virial Theorem. We then compute where K 2 ( · ) denotes the modified Bessel function of the second kind [25], section 9.6. Bringing this together we find that Now we need to minimize (3.22) with respect to b. It is not immediate how to do this, since where K 1 ( · ) and K 2 ( · ) denote the modified Bessel functions of the first and second kind, respectively. To calculate ( ) ¢ F b we have used the well-known derivative formulas for the modified Bessel functions, see e.g. [25], page 376. So instead, we seek an upper bound for (3.22). This essentially reduces to finding a lower bound for K 2 ( · ). To this end, we let = y C b. Then consider the function which comes from the previous substitution for the second term in F(b). We change variables by letting = w y 2 2 so then we consider To estimate K 2 (w), we appeal to corollary 3.4 from [15], which gives the following sharp estimate: and p 2 is the best constant. Then, (3.23) implies This means that This provides an upper bound for the ground state energy E 0 since E 0 〈ψ 0 |H|ψ 0 〉.

Using the Bohr radius
In this section we instead use the exact Bohr radius p =   a me Taking derivative with respect to b yields Setting the numerator equal to zero yields Now, for small arguments, it is known that where Γ( · ) denotes the Gamma function. Thus and thus Using the value of b from (3.29) in (3.25) yields

Bound States
The number of bound states for a given system with potential V is a deep rooted question in quantum mechanics. For the classical Coulomb potential, it is well known that there are infinitely many bound states. However, for other potentials, it is much less clear. For instance, consider the 1/r 2 potential. In [26] it was shown that no discrete bound states exist for systems at low energies described by such a potential. This is deeply rooted to the problem of constructing self-adjoint extensions [27,28].
For the Yukawa potential, the issue of bound states is an extremely interesting question, given the existence of the critical screening parameter [17]. The problem of counting the number of bound states for this potential has been considered by many, and depends on the dimension. We begin our analysis in dimensions 1 and 2.

Dimensions 1 and 2
We show that there are infinitely many bound states in 1D. The starting point is the zero-energy Schrödinger equation for some potential U(x). In general, this can be transformed into the following 2D zero angular momentum Schrödinger equation through a change of variables (a Liouville transformation to be precise, see e.g. [29], equation (13)) which is known to have infinitely many bound states, see [29], appendix B where this is shown in detail. Since (see figure 2) ( ) -< e r W r r for sufficiently large C r we must have that in 1D, V P has infinitely many bound states. This is because for any two potentials V 1 , V 2 0, if V 1 (r) > V 2 (r) for r > R, then the number of nodes n 2 for V 2 in r > R satisfies n 2 n 1 − c for some constant c, where n 1 denotes the number of nodes for V 1 in r > R (see e.g. [29], equation (4)).
In 2D we have the same result: there are infinitely many bound states; this follows directly from the previous substitutions and equation (15) in [29]; namely, that the potential ( ) = This result is for a general 3D spherically symmetric potential V(r) such that the integral I converges. Since then many attempts have been made at getting sharper bounds for the number of bound states; see e.g. [16] for a review of such results. Unfortunately, for V P (r), Bargmann's result cannot be applied since the integral ò ¥e dr C r 0 diverges. Thus the question is: do bound states exist for a system with potential V P ? If so, how many are there? Since V P is radial, the number of negative energy bound states is equal to the number of nodes of the wavefunction solution of (4.34). We further see numerical evidence for the existence of a bound state; see figure 3 below for the numerical solution of (4.34) with C = 0.1 (a Runge-Kutta solver is used). In fact, this can be made more rigorous by the results of [30]; namely, a bound state can exist provided that one can construct a trial wavefunction ψ t such that 〈H〉 < 0. Thus, we construct a function ψ t such that satisfies 〈H〉 < 0, i.e. 〈T〉 < − 〈V〉. Define the function for a parameter α to be determined. With this wavefunction, we compute (in units where ÿ = m = 1): where Shi(z) and Chi(z) denote the hyperbolic sine and cosine integrals [25], section 5.2, respectively: where γ is Euler's constant. It turns out that for α sufficiently large, one can find C > 0 such that 〈T〉 < − 〈V〉. Indeed, taking α = 100 and C = 0.1, one can evaluate explicitly that 〈T〉 + 〈V〉 = − 0.000730032 < 0. Thus, there exists at least one bound state. It is not clear, though, if there are infinitely many bound states in 3D. Another natural question is the following: does there exist a critical screening parameter C 0 such that bound states cease to exist for C > C 0 ? Now, it is well known that for the classical Coulomb potential in dimension three there are infinitely many bound states. However, for the Yukawa potential, there does exist a critical screening parameter, see e.g. [17,19] and the references therein. We expect something similar in this case: the number of bound states should be limited for C ≠ 0 since interactions are screened.

Case ℓ = 0
Under a suitable change of variables, i.e. u(r) = rψ(r), we can change the 3D equation (4.34) into the following problem for the reduced wavefunction u: The first eigenfunction of (4.37) is plotted in figure 4. Figure 4 was obtained using the Matlab package MATSLISE [31], which uses piecewise constant perturbation methods. The first four eigenfunctions can be seen below in figure 5: all appear to be square integrable.
The eigenvalues of (4.37) can also be computed numerically; with an imposed tolerance of 1 × 10 −8 , we obtain in table 1 the first 10 eigenvalues. Note in particular that the computed eigenvalues are well approximated by the well-known formula for the Hydrogen atom E n = E 0 /n 2 . It is also interesting to note that choosing higher values of C than that used in table 1 show a similar dependence of the eigenvalues on n, which is quickly checked numerically.     When ℓ = 0, we expect to see infinitely many bound states. Indeed, let us further analyze the zero-energy wavefunction. As mentioned before, the number of negative energy bound states is equal to the number of nodes of the zero energy wavefunction, solution of Numerical simulations indicate that, at a maximum integration interval of ( ) r 0, max with = r 1000 max , there are already 19 nodes of the wavefunction; see in particular figure 6 below.
In table 2 we provide a table of nodes of the wavefunction with varying values of the maximum integration interval (the numerical solver uses a Runge-Kutta method).
We see that the number of nodes keeps increasing as r → ∞. Note also that (4.31) is the same as the ℓ = 0 reduced Schrödinger equation in dimension three (except restricted to r > 0), see equation (4.38). Thus we indeed see infinitely many bound states when ℓ = 0.

Case ℓ ≠ 0.
When ℓ ≠ 0, for certain values of C, there is a well, and bound states could exist; see figure 7. See also [19] for a similar analysis for the Yukawa potential. In that case, there is a critical parameter in which bound states cease to exist due to the effective potential becoming repulsive at some point. The question is, do we see the effective potential ℓ ℓ ( ) -+ + e r r 1 C r 2 becoming repulsive? Figure 8 suggests this is indeed the case, i.e. for fixed C, for ℓ large enough, no bound state can exist here. This suggests there does exist a critical screening parameter for which bound states cease to exist when ℓ ≠ 0. This certainly warrants a much deeper investigation; in particular, for fixed ℓ, does there exist C * such that for C > C * no bound states exist? This would require an in-depth analysis of the relationship between the screening parameter C and quantum number ℓ.

Extended Nikiforov-Uvarov method
In [1], the radial equation for the Hydrogen atom was shown to be quasi-exactly solvable. Such systems can be solved using an extended Nikiforov-Uvarov method [32,33]. Following [32], we can calculate the eigen-energies E = E n,ℓ in the radial Schrodinger equation are polynomials of at most second degree, and ( )  t s is a first degree polynomial. We rewrite (5.39) as (dropping indices for convenience) To obtain an equation of the form (5.40), we make the substitution The derivatives transform as   Remark 5.1. Note that the assumption  C r 1 is equivalent to  r C, i.e. r is assumed large enough. Now, as discussed in [34], large distance behavior is not enough to determine the eigenvalues for a Hydrogen Hamiltonian; as such, we also need to impose a boundary condition at the origin. This is accounted for through the general method developed in [32], and comes through in choosing the physical acceptable solution for π in (5.43). Dividing through by ( ) s Dividing by s 2 yields Multiply by C 2 gives We can rewrite this as and so in the notation of (5.40), we have: Notice that σ(s) is actually a polynomial of degree three now. Following [33] (which is necessary since σ is a degree three polynomial rather than degree two), we define the polynomial

The radial equation then becomes
where g(s) is a polynomial chosen such that the quantity under the square root is a square of a polynomial of degree at most two. We find that π(s) in this case is given explicitly by Suppose first that ℓ = 0 for simplicity. We need to find a polynomial g so that the radicand of π in (5.43) is the square of a polynomial of second degree. So, suppose that for some coefficients a, b, f. Since g will also need to be a polynomial, suppose that g(s) = ms + n; it is clear that g cannot be a quadratic polynomial. Thus we wish to solve for the unknown coefficients a, b, f, m, n. We do this by comparing like terms on both sides of (5.44); in this way, f is easily found by comparing the constant terms: The rest of the coefficients can be found by solving the following system:  with f given by (5.45) and Then from (5.43), we see that where C n are integration constants. The method in [33] implies that h n (s) = h(s), so we are able to obtain the eigenvalue equalities: The system (5.49)-(5.50) gives an implicit relationship between E and C.