Towards the Theory of the Yukawa Potential

# Towards the Theory of the Yukawa Potential

J. C. del Valle Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, A. Postal 70-543 C. P. 04510, Ciudad de México, México.    D. J. Nader Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, A. Postal 70-543 C. P. 04510, Ciudad de México, México.
###### Abstract

Using three different approaches, Perturbation Theory (PT), the Lagrange Mesh Method (Lag-Mesh) and the Variational Method (VM), we study the low-lying states of the Yukawa potential . First orders in PT in powers of are calculated in the framework of the Non-Linerization Procedure. It is found that the Padé approximants to PT series together with the Lag-Mesh provide highly accurate values of the energy and the positions of the radial nodes of the wave function. The most accurate results, at present, of the critical screening parameters () for some low-lying states and the first coefficients in the expansion of the energy at are presented. A locally-accurate and compact approximation for the eigenfunctions of the low-lying states for any is discovered. This approximation used as a trial function in VM eventually leads to energies as precise as those of PT and Lag-Mesh. Finally, a compact analytical expression for the energy as a function of , that reproduce at least decimal digits in the entire physical range of , is found.

## I Introduction

The Yukawa potential, sometimes called the screened Coulomb potential, has a wide range of applications in many branches of physics. Originally, Yukawa Yukawa () proposed this potential to describe the interaction between a pair of nucleons. However, it is often used as a first approximation of the interaction between two dust particles immersed in a plasma Plasma1 (); Plasma2 (); Plasma3 (). The interaction of dark matter particles through a Yukawa potential could explain the recently observed cores in dwarf galaxies Dark (). The Yukawa potential has also been employed to model the systems of colloidal particles in electrolytes Colloidal ().

In the standard form the three-dimensional Yukawa potential is given by

 V(r) = −λe−αrr , (1)

where and are parameters and is the radial coordinate. If , the Yukawa potential degenerates into the Coulomb potential. When , the Yukawa potential vanishes: it describes a free motion.

In the non-relativistic approximation the radial Schrödinger equation (in atomic units ) with the Yukawa potential (1) reads

 [−12d2dr2 − 1rddr+ l(l+1)2r2 − λe−αrr]Ψn,l = En,l(λ,α)Ψn,l , (2)

where , , is the principal quantum number and , , …, is the angular momentum. It is well known that the Schrödinger equation (2) with is not exactly solvable and only holds a finite number of bound states. Naturally, this number depends on the parameters and . The scaling transformation in the equation (2) allow us to remove the dependence of the energy since it is scaled as

 En,l(λ,α) = λ2En,l(1,αλ−1). (3)

From now on we set . Consequently the number of bound states as well as the energy depends only on the value of .

There is a certain critical value of , denoted by , with the following property: for the state is bound being characterized by a normalized wave function, otherwise if the energy of the state passes to the continuous regime and the corresponding state is no longer bound, i.e. the wave function becomes non-normalizable. Hence, for the state the critical (screening) parameter is such that

 En,l(α(n,l)c) = 0 . (4)

There is a large number of papers in the literature devoted to the calculation of the critical screening parameters . Let us list some of the approaches used for this purpose: pseudo spectral method IJQC (), variational calculations lam (); stubbins (); variational (), perturbation theory Logarithmic (), finite scaling finite (), matrix propagation matrix (), direct solution of the Schrödinger equation Rogers (); Weber () and various numerical methods. At the moment there are certain discrepancies in the critical parameters , in general they differ in the sixth decimal digit. A brief historical discussion is given in Weber () and complemented with katriel () about some estimates of the critical parameter of the ground state .

Following the results SIMON (), the energy can be represented as a series expansion at . For states with this expansion is

 En,0(α) = α2∞∑k=2β(n,0)kα−k(α(n,0)c−α)k, (5)

while for states with it reads

 En,l(α) = α2∞∑k=2β(n,l)k−1α−k/2(α(n,0)c−α)k/2, (6)

where and , , are real coefficients. One of the aims of this study is to calculate the critical screening parameter and to find the first terms in the expansion (5) and (6). We focus on the low-lying states with quantum numbers and .

Contrary to the situation of , the first coefficients in (5) and (6) are not known. The simplest way to estimate them is to construct an interpolating function which has the form of terminated expansions (5) and (6) and find coefficients by describing the energy in a close vicinity of . In this manner the coefficients are free parameters to be adjusted in a fit. As we will discuss in Section III.C, and are related to the Hellmann-Feynman theorem, therefore their calculation turns out to be appropriate to measure the accuracy of an approximate wave function corresponding to via an expectation value.

Our calculations of are carried out using two approaches: PT in powers of and Lag-Mesh. The first one is PT within the framework of the Non-Linearization Procedure Turbiner1984 (). Due to the probably divergent nature of this PT we will use a continued fraction representation for the perturbation series Vrscay (). The second one is Lag-Mesh Baye1 () that, within its applicability domain, is one of the most accurate numerical methods to solve the Schrödinger equation.

To compute and , , we proceed with the interpolation procedure already described. The calculation of the first terms in expansions (5) and (6) plays a fundamental role in the construction of an analytical approximation to via rational functions. We present such approximations for some states. Naturally, highly accurate energy estimates are needed to construct accurate analytical approximations. Therefore we complement the energy estimates that come from PT and Lag-Mesh with variational calculations. We construct a locally accurate and compact trial function for the low-lying states based on the interpolation of the series expansions of the wave function at and . The Non-Linearization Procedure allows us to estimate the accuracy of the variational calculations for the energy and also for the trial function.

It is found that the three approaches (PT, Lag-Mesh and VM) are also appropriate to calculate the real radial nodes {} of the wave function where it vanishes

 Ψn,l(r0) = 0 . (7)

A precise calculation of the real nodes (and more general, the nodal surface) is relevant in quantum mechanics. For instance, in nodes () it is shown that the nodes plays an essential role for the variational estimation of upper bounds of the energy of excited states. A similar situation occurs with fermionic systems studied with the diffusion Monte Carlo method within the fixed-node approximation DMC (). As we will present, the Non-Linearization Procedure is adequate to construct a PT for the factor of the wave function that defines the radial nodes. Again, since PT in powers of is probably divergent, we use the continued fraction representation to estimate . It is found that the results of PT are in excellent agreement with those provided by the Lag-Mesh and VM.

## Ii The methods

### ii.1 Perturbation Theory and its Summation

#### ii.1.1 Non-Linerization Procedure

This perturbative approach, sometimes called Logarithmic PT, has shown to be a useful tool in non-relativistic quantum mechanics, see e.g. Turbiner1984 () and Turbiner1984H (). In particular, it is a very efficient method to calculate perturbation series. Let us give a brief review of this procedure by considering a generic spherically symmetric potential .

For convenience we use the following representation of the radial part of the wave function,

 Ψn,l(r) = rlfn,l(r)e−Φn,l(r) . (8)

Note that if we restrict for finite then the function completely characterizes the nodal surface of any excited state . Similar as in the one-dimensional case Turbiner1984 (), we require to be a polynomial of degree () with real (and non-negative) zeros. Substituting (8) in the radial Schrödinger equation (2), we obtain the non-linear differential equation

 y′n,l − yn,l(yn,l − 2(l+1)r) + 2(yn,l − l+1r)f′n,l − f′′n,lfn,l = 2(En,l−V) , (9)

where . The equation (9) is the starting point to develop PT. First, we assume that the potential can be written as an expansion in powers of some parameter ,

 V(r;γ) = ∞∑k=0V(k)(r)γk , (10)

where are some given functions. The energy is also expanded in series of powers of ,

 En,l = ∞∑k=0E(k)n,lγk , (11)

as well as the functions and ,

 fn,l(r) = ∞∑k=0f(k)n,l(r)γk ,yn,l(r) = ∞∑k=0y(k)n,l(r)γk . (12)

One can show that the linear differential equation that determines the th corrections , and is

 (f(k))′′+2(l+1)r(f(k))′−2k∑i=0(f(i))′y(k−i) +k∑i=0{f(i)(k−i∑j=0y(j)y(k−i−j)+2(E(k−i)−V(k−i))−(y(k−i))′−2(l+1)ry(k−i))}= 0

where . The boundary condition

 y(k)n,le−Φ(0)n,l∣∣r=0,{r0},∞→ 0 (14)

must be imposed Turbiner1984 (). For convenience we have omitted the labels in some expressions. Interestingly, the equation (LABEL:eq:kcorrection) can be further developed to obtain an integral representation for the th energy correction Turbiner1984 (). One should note that this approach does not require the previous knowledge of the entire spectrum of the unperturbed () equation to construct perturbation series.

For the Yukawa potential (1) we set and consequently for the expansion (10) we have

 V(k)(r) = (−r)k−1k! . (15)

With the expansion (10) corresponds to the Coulomb potential. In this case the quantum numbers take the values and . The zeroth order corrections , and are well known,

 E(0)n,l = −12n2 ,f(0)n,l(r) = L2l+1n−l−1(r) ,y(0)n,l = 1n , (16)

where is the generalized Laguerre polynomial of degree . For convenience we will denote

 L2l+1n−l−1(r)=n−l−1∑i=0a(0)n,l[i]ri (17)

and choose as a normalization.

#### ii.1.2 Continued Fractions and Padé Approximants

For the ground state of the Yukawa potential, it is well known that PT (11) in powers of is divergent Logarithmic (). For excited states this PT is probably divergent. This statement is also true for the series (12). Even though perturbation series (11) are not Stieldjes Vrscay (), the continued fraction representation can be used to calculate accurately the energy of a given state Lai (); Vrscay (). The basic idea is to assume that the energy has the representation

 En,l(α) = c0 + c1α1 + c2α1 + c3α1 + ⋯ , (18)

where , with , are coefficients to be determined. For convenience we omit the label () in the coefficients . In practice, the continued fraction is truncated by setting if , where is some positive integer. As a result of this truncation, a Padé approximant of the form emerges,

 P⌈M/2⌉⌊M/2⌋(α) = ∑⌈M/2⌉k=0Ckαk∑⌊M/2⌋k=0Dkαk, (19)

where and denote the ceiling and floor functions, respectively. The coefficients and are determined by demanding

 M∑k=0E(k)n,lαk − P⌈M/2⌉⌊M/2⌋(α) = O(αM+1) . (20)

Once the Padé approximant is completely determined it is used to calculate the energy , the value of the critical parameter and the coefficients and and . Additionally, we will show that if the continued fraction representation is assumed for the perturbation series of , see (12), we can calculate the position of the nodes with high accuracy via Padé approximants.

### ii.2 The Lagrange Mesh Method

The Lag-Mesh has shown to be a simple and very accurate method to solve the Schrödinger equation, see Baye1 (); Baye2 (). Essentially, the wave function is expanded in terms of the Lagrange functions while the Gauss quadrature is used to calculate approximately the matrix elements of the Hamiltonian. Once the matrix elements are known, we proceed to calculate the eigenvalues and the corresponding eigenfunctions. In the next two Subsections we give a brief review of the method nevertheless it is described in full detail in Baye1 (); Baye2 (); Baye3 ().

#### ii.2.1 The Lagrange Functions and the Gauss Quadrature

A mesh of dimension involves real zeros of a particular orthogonal polynomial of degree . Given the values of a function at , the polynomial of minimal grade , denoted by , which interpolates the function is of the form

 LN−1(r) = N∑i=1F(ri)fi(r) , (21)

where the Lagrange functions are defined by

 fi(r) = PN(r)(r−ri)P′N(ri) . (22)

Since are the roots of the Lagrange functions satisfy the property , and therefore . The integral of the function in the domain can be approximated using the Gauss Quadrature as follows

 ∫baF(r)dr ≈ N∑iλiF(ri), (23)

where are the associated weights. The Gauss quadrature provides high accuracy on the integrals except when the function contains singularities or discontinuities Baye1 ().

#### ii.2.2 The Lag-Mesh in Quantum Mechanics

For spherically symmetric potentials it is convenient to transform the radial Schrödinger equation (2) into its one-dimensional counterpart. If we assume a wave function of the form , then the function satisfies

 [−12d2dr2 + U(r)]un,l = En,lun,l , (24)

with the effective potential given by

 U(r) = V(r) + l(l+1)2r2 . (25)

Now, we consider a wave function given as an expansion

 un,l(r) = N∑i=1ci^fi(r) , (26)

where are coefficients and are the regularized Lagrange functions

 ^fi(r) = λ−1/2irriw(r)1/2f(r) , (27)

here is the weight function associated to the orthogonal polynomial . One can immediately notice that in this representation the function always vanishes at the origin. The other boundary condition, as , is satisfied by choosing the appropriate polynomial . Since we are interested in the domain the Laguerre mesh is adequate. Therefore, we set , where is the th Laguerre polynomial and the weight function is .

The coefficients are determined by the secular equation related to (24),

 N∑j=1{Tij + Uij }cj =Eci , (28)

where the matrix elements are

 Tij = −12∫^fi(r)d2dr2^fj(r)dr ,Uij = ∫^fi(r)U(r)^fj(r)dr . (29)

With the Gauss quadrature the matrix elements of the effective potential are found

 Uij = U(ri)δij = (V(ri) + l(l + 1)2r2i)δij . (30)

Therefore, we can see that the matrix representation of is diagonal in the Gauss approximation. For the kinetic matrix elements , the discrete variable representation Szalay () of the operator is useful to obtain the elements in closed form within the Gauss quadrature,

 Tii = 4 + (4N + 2)ri − r2i24r2i ,i=j , (31) Tij = (−1)i−j(ri + rj)2(rirj)1/2(ri − rj)2 ,i≠j . (32)

It is worth mentioning that the use of the regularized Lagrange functions circumvents the error of the Gauss quadrature induced by the singularity of the Yukawa potential at . Once all the matrix elements are known, we proceed to solve the secular equation (28) for the Yukawa potential (1). From the solution of the secular equation we obtain the first approximate wave functions and their corresponding energies for a fixed angular momentum and parameter . Interestingly, the expectation value of any function can be computed (in the Gauss quadrature approximation) as follows

 ⟨Ψn,l|g(r)|Ψn,l⟩ = N∑i=1|ci|2g(ri) . (33)

### ii.3 Trial Functions

In order to design a trial function for the state we will follow the approach presented in TURBINER2005 (). In the latter, a trial function was constructed for the one-dimensional quartic anharmonic and double-well potentials which leaded to the most accurate variational energy estimates of the ground state. The basic idea is to construct a minimal interpolation between the expansions of the wave function at and .

We begin the construction by considering the ground state wave function in the representation (8),

 Ψ1,0(r) = e−Φ1,0(r) . (34)

Using (2), it is straightforward to show that the function satisfies a non-linear differential equation, namely

 Φ′′1,0−Φ′1,0(Φ′1,0−2r)=2(E1,0+e−αrr) . (35)

From this equation one can construct the series expansions of the function at and . These expansions read

 Φ1,0(r) =  r + 16(1−2α+2E1,0)r2 + 136(2−4α+3α2+4E1,0)r3 + … , (36)
 Φ1,0(r) = √−2E1,0r + lnr − 4α(α+2√−2E1,0)e−αrr − 4α(α+2√−2E1,0)2e−αrr2 + … , (37)

for and , respectively. In the limit both expansions are truncated and they coincide

 Φ1,0(r) = r . (38)

Together with (34), this is nothing but the ground state wave function of the hydrogen atom. One of the simplest interpolations between (36) and (37) is given by

 Φ(t)1,0 = r(a1,0r + b1,0e−αr + c1,0e−2αrd1,0r + e−ar) + log(d1,0r + e−αr) (39)

where are variational parameters. In this manner, our trial function for the ground state is

 Ψ(t)1,0(r) = e−Φ(t)1,0(r) . (40)

For excited states we use as a building block to construct the trial function of the state

 Ψ(t)n,l(r) = f(t)n,l(r)e−Φ(t)n,l(r) ,f(t)n,l(r) =rlP(t)n,l(r) , (41)

where is a polynomial of degree . Here has the same functional structure as (39) with a different set of parameters . We impose the constraint that the trial function must be orthogonal to the functions , , , . This constraint fixes the value of some parameters of . The remaining free parameters are taken as variational parameters. This parameters are adjusted such that the expectation value of the radial Hamiltonian

 ^h = −12d2dr2 − 1rddr + l(l+1)2r2  −e−αrr, (42)

is minimal.

Interestingly, there is a connection between PT and VM Turbiner1984 () that we explain briefly. Any potential can always be rewritten as

 V(r) = V0(r) + γ(V(r) − V0(r)) ,γ=1 , (43)

where is formal parameter and is the potential for which the trial function is the exact solution of the Schrödinger equation

 [−12d2dr2 − 1rddr + l(l+1)2r2 + V0(r)]Ψ(t)n,l = E(0)n,lΨ(t)n,l . (44)

In equation (43)

 V(1)(r) = V(r) − V0(r) (45)

plays the role of the perturbation potential. Consequently, the variational energy corresponds to the first two terms of a perturbative series,

 ∫∞0Ψ(t)n,l^hΨ(t)n,lr2dr∫∞0(Ψ(t)n,l)2r2dr = E(0)n,l + γE(1)n,l ,E(1)n,l = ∫∞0Ψ(t)n,lV(1)Ψ(t)n,lr2dr∫∞0(Ψ(t)n,l)2r2dr . (46)

Since the Non-Linearization Procedure only requires as entry the unperturbed wave function, we can take the trial function as an unperturbed wave function and then develop PT in order to construct higher order corrections for , and . In this manner we can estimate the accuracy of our calculations by means of the perturbative corrections. As a consequence we can define (and calculate), for example

 En,l;1 = E(0)n,l + E(1)n,l (47)

and

 En,l;2 = E(0)n,l + E(1)n,l + E(2)n,l, (48)

which correspond to the first and second order approximations to the energy. In general we can define the sum of the corrections

 En,l;i = i∑k=0E(k)n,l (49)

as the -th approximation to the exact energy . In particular, if the trial function is chosen appropriately a convergent PT occurs Turbiner1984 (),

 limi→∞En,l;i = En,l . (50)

Similar approximations can be defined for the functions and . However, if the exact position of the node is known and this information is codified in the trial function, no correction of should be calculated. This situation occurs for the states with quantum numbers () which only have a node of order at .

## Iii Results

### iii.1 Some Remarks on PT

We begin by presenting some results about the realization of PT in the framework of the Non-Linearization Procedure.

The first order corrections always vanish, . We found that the corrections and are both polynomials in of the form

 f(k)n,l(r) = n−l−2∑i=0a(k)n,l[i]ri ,k>1 , (51)

and

 y(k)n,l(r) = k−1∑i=0b(k)n,l[i]ri ,k>1 , (52)

respectively. The coefficients and are always real and rational numbers. Since the corrections are polynomials the realization of PT is an algebraic and iterative procedure. More precisely, using expressions (51) and (52), the differential equation (LABEL:eq:kcorrection) transforms into an algebraic equation for , and the energy correction . The expansion of in powers of is an alternating series and the coefficients , , are also rational numbers. The algebraic nature of the realization of the PT allowed us to calculate high orders in PT by using the software Mathematica. The Padé approximants to the series can be easily constructed using also this software. For some selected states the first 10 corrections are presented in the Table 1. However, for the states considered, we were able to calculate exactly the first 400 perturbative corrections in (11) and (12). For example, for the energy of the ground state, some high-order coefficients (rounded) are: , , and .

From the polynomial structure of the corrections we corroborate that the function is also a polynomial,

 fn,l(r) =rln−l−1∑k=0An,l[k]rk , (53)

where

 An,l[n−l−1] = 1 (54)

and

 An,l[k] = ∞∑i=0a(k)n,l[i]αi% for k

Since both series, for the energy and for coefficients (55) of , are probably divergent, we used a continued fraction representation (and hence Padé approximants) for the summation. Therefore, once the coefficients (55) are known, the radial nodes correspond to the real roots of the polynomial (53). At Subsection D we present numerical results of the position of the nodes.

### iii.2 Critical Screening Parameter α(n,l)c

We carried out calculations of the energy of the low-lying states of the Yukawa potential using PT (+ Padé approximants) and Lag-Mesh. For the Lag-Mesh calculations, we wrote a computational code in Fortran 90. The roots , , that define the mesh were calculated with Mathematica. We diagonalized the matrix representation of the secular equation (28) using the DSYEV routine of LAPACK Diagonal ().

Far enough from , the Padé approximants provide high accuracy estimates of the energy . As regards to the Lag-Mesh, a dimension generates accurate results that agree with those of PT.

As the parameter approximates to a remarkable difference between the energy obtained by the two methods appears. Near the state becomes weakly bound, the corresponding wave function is very flat and extended: an extension of the configuration domain of the variable is required. In PT an extension of the domain is achieved by calculating higher order terms in the expansion of the energy (11), while in the Lag-Mesh we need a larger dimension of the mesh, i.e. a larger value of . It is found that the first 100 coefficients in the expansion (11) provide, via Padé approximants, high accuracy for states with near . Otherwise, when , even in (19) is not large enough to reach the accuracy of the Lag-Mesh. In any case, we only present results with diagonal Padé approximants of the form .

For the results of Lag-Mesh we present only stable digits with respect to variations of the dimension . The largest value of considered was . In the case of we scale with a positive parameter in the form . The parameter is such that the first positive energy is minimal.

The calculations of the energy are presented in Tables 2 - 3. In Table 2 we compare the energy of the ground state obtained by the two methods. One can see that the difference appears in the 15th decimal digit. The same difference in digits appears in the energy of the excited states. In Table 3 we only present the digits that are in agreement using both methods. In general, for sufficiently far from the critical value , both methods provides at least 14 decimal digits.

The critical values of are presented in Table 4. The results are compared with those of matrix (). For states with angular momentum , we find that Padé approximants give more accurate estimates of in comparison with those of Lag-Mesh (see results of matrix ()). Otherwise, when , the Lag-Mesh provides more accurate results in comparison with those of PT. As one can see in Table 4, the critical parameters obtained with Lag-Mesh exhibit more precision as the angular momentum increases.

### iii.3 Expansion of the Energy at α(n,l)c

Near the critical parameter , the energy of the states can be represented as an expansion (5), while the energy of the states () is represented as an expansion of the type (6). Interestingly, both expansions (5) and (6) can be re-expanded at to obtain a Taylor series and a Puiseux series in half-integer powers, respectively,

 (56)

The coefficients of (5) and (6) are related with the coefficients of (56) as follows

 l≠ 0 →~β(n,l)2=β(n,l)2,~β(n,l)3=β(n,l)3α(n,l)c, l≠ 0 →~β(n,l)1=β(n,l)1α(n