A spectral Phase–Amplitude method for propagating a wave function to large distances

Computer Physics Communications 191 (2015) 33–42

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

A spectral Phase–Amplitude method for propagating a wave function to large distances George Rawitscher Physics Department, University of Connecticut Storrs, CT 08269, United States

article

info

Article history: Received 30 October 2014 Received in revised form 12 December 2014 Accepted 26 January 2015 Available online 4 February 2015 Keywords: Phase–Amplitude representation Spectral methods Computational techniques Atomic scattering

abstract The phase and amplitude (Ph–A) of a wave function vary slowly with distance, in contrast to the wave function that can be highly oscillatory. Hence the Ph–A representation of a wave function requires far fewer computational mesh points than the wave function itself. In 1930 Milne presented an equation for the phase and the amplitude functions (which is different from the one developed by Calogero), and in 1962 Seaton and Peach solved these equations iteratively. The objective of the present study is to implement Seaton and Peach’s iteration procedure with a spectral Chebyshev expansion method, and at the same time present a non-iterative analytic solution to an approximate version of the iterative equations. The iterations converge rapidly for the case of attractive potentials. Two numerical examples are given: (1) for a potential that decreases with distance as 1/r 3 , and (2) a Coulomb potential ∝1/r. In both cases the whole radial range of [0–2000] requires only between 25 and 100 mesh points and the corresponding accuracy is between 10−3 and 10−6 . The 0th iteration (which is the WKB approximation) gives an accuracy of 10−2 . This spectral method permits one to calculate a wave function out to large distances reliably and economically. © 2015 Elsevier B.V. All rights reserved.

1. Introduction When the Phase–Amplitude (Ph–A) method was first introduced by Milne in 1930 [1], and then taken up by many authors, see Ref. [2], the main motivation was the paucity of numerical mesh points required, compared to the calculation of the wave function itself. This is because both phase and amplitude functions are slowly varying, as opposed to the wave function that can be highly oscillatory. This point was verified by many authors, in particular by Calogero and Ravenhall [3] who state that the solution for the phase is more stable than the solution of the wave function. An additional argument in favor of the Ph–A representation is that it lends itself to analytic expressions to address particular problems. For example, the Ph–A representation facilitates the incorporation of the effect of long range potentials [4,5] or the calculation of resonances [2]. It is also helpful in the quantum defect calculation of atomic wave functions [6], the calculation of Gaunt Factors [7], as well as the description of an electron with an ion embedded in a plasma [8], among others. The emission of an electron from an atom, either by incident photons, or by other processes such as beta decay of the nucleus also makes extensive

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.cpc.2015.01.014 0010-4655/© 2015 Elsevier B.V. All rights reserved.

use of Calogero’s Ph–A description [9]. However that method studied by Calogero [10] is different from the one described here, as is further discussed in Appendix A. The Ph–A description of a carrier wave in radio or television also plays a significant rôle in the compactification of the signal transmission in the field of Information Technology [11]. An additional advantage of the present Ph–A representation is that it provides a method to improve the Wentzel, Kramers Brillouin (WKB) [12,13] approximation of a wave function, an important point since the WKB approximation has led, over the years, to a much improved understanding of the solution of the Schrödinger equation. A further advantage for the case of long range potentials is that Milne’s Ph–A method lends itself to provide the normalization of a conventionally obtained wave functions calculated only out to short distances, such that the amplitude of that wave function would asymptotically approach unity, if carried out to large distances. This normalization method does not require a Wronskian-type matching procedure to two analytically known basis functions. The Ph–A representation consists in writing a wave function ψ(r ) in the form

ψ(r ) = y(r ) sin[φ(r )],

(1)

where y is the amplitude and φ is the phase, and r the distance from the origin. Here ψ is a partial wave function that depends only on

34

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

2. Iterative solution of Milne’s phase–amplitude equation

the single variable r. If an overlap matrix element ∞



ψ1 (r )U (r )ψ2 (r )dr

M =

(2)

0

between two wave functions is required, then in the finite difference method of obtaining integrals, both ψ1 and ψ2 have to be calculated on a sufficiently fine mesh, which can be time consuming and prone to errors. However, the Ph–A representation can provide an estimate of M by decomposing the integrand of the overlap matrix element into a slowly oscillating (S) and a fast oscillating (F ) part M = M (S ) − M (F ) .

(3)

The decomposition makes use of a trigonometric identity for the product of two sine functions with the result M

(F )

=

M (S ) =

1 1 2

d2 ψ/dr 2 + k2 ψ = VT ψ.

(6)

The total potential VT is VT (r ) = L(L + 1)/r 2 + V (r ),

(7)

where V (r ) is the atomic or nuclear potential (including the Coulomb potential) in units of inverse length squared, and L is the orbital angular momentum quantum number. Milne’s non-linear equation for the amplitude y(r ) is given by [1] d2 y/dr 2 + k2 y = VT y + k2 /y3 ,

(8)

∞



2

The Schrödinger equation to be solved for a partial wave function ψ is

y1 (r )U (r )y2 (r )[cos(φ1 + φ2 )]dr .

(4)

y1 (r )U (r )y2 (r )[cos(φ1 − φ2 )]dr .

(5)

0

∞

 0

The matrix element M (S ) (S stands for slow variation, F for fast variation of the integrand) can be calculated on a small set of radial mesh points since the integrand oscillates slowly. Further, since |M (F ) | < |M (F ) |, a rough estimate for M is provided by M (S ) alone. Here U (r ) is an overlap function that depends on the physics application envisaged. In 1962 Seaton and Peach [14] presented an iterative scheme to solve Milne’s non-linear differential equation [1] for the amplitude and phase. It is the purpose of the present work to implement this iterative method by means of a spectral [15,16] expansion of the amplitude in terms of Chebyshev polynomials. A further purpose is to examine the accuracy of the resulting Ph–A wave function by comparison with the direct solution of the Schrödinger equation for the wave function, the latter is also obtained by an accurate spectral integral equation method [17,18], denoted as IEM in what follows. The combination of both objectives has not been presented previously. The great advantage of a spectral expansion is that the calculations utilize all the support points located in a given partition simultaneously, with the result that the errors are shared uniformly across the partition in the case of Chebyshev expansions [19,20]. For the present numerical examples the calculation is done in one great radial partition, extending from r = 0 to r = 2000, containing between 25 and 201 Chebyshev support points, depending on the accuracy required. By contrast, other algorithms (such as finite elements, finite differences, or the IEM method described below) divide such a large radial interval into a number of partitions, with the result that the accumulated error from all previous partitions is propagated into the next one, the last partition having the largest error [21]. An analysis of the accumulation of the Finite Element errors is presented in Appendix B. In addition, for calculations that require the storage of many wave functions with high precision [22–25] the use of the Ph–A representation can be very advantageous because the amount of storage required can be substantially smaller than what is needed for other algorithms. In Section 2 the iterative method is explained, Section 3 contains details of the computational spectral method, Section 4 describes a non-iterative analytic solution to an approximate set of iterative equations, Section 5 establishes a connection between a wave function and the corresponding values of the phase and amplitude, Section 6 contains numerical results, the calculation of overlap integrals is described in Section 7, the Summary and Conclusions are given in Section 8. In Appendix A it is shown that the Ph–A combination defined by Calogero is not the same as the one defined by Milne, even though both methods lead to the same wave function ψ , and in Appendix B it is argued that a Gauss–Lobatto finite element method cannot reach large distances with the same accuracy as the Ph–A method.

where the nonlinearity is given by the last term in Eq. (8). In Eqs. (6) to (8) the factor }2/2m has already been divided into the potential and into the energy, so that both are given in units of inverse length squared, and the wave number k is given in units of inverse length. The unit of length can be either fm for nuclear physics applications, or the Bohr radius a0 for atomic physics applications, but will not be explicitly indicated. The phase φ(r ) is obtained from the amplitude y by [1]

φ(r ) = φ(r0 ) + k



r

[y(r ′ )]−2 dr ′ ,

(9)

r0

but it can also be obtained without the knowledge of y [7]. Eqs. (8) and (9) can be obtained by inserting Eq. (1) into Eq. (6), noting that the terms involving the phase φ(r ) can be separated from the terms involving the amplitude y(r ), and each collection of terms can be set to zero independently. An overall normalization is still arbitrary, but it can be fixed by demanding close agreement with the WKB approximation [12,13]. Eqs. (8) and (9) lead to a phase–amplitude representation of the wave function that is different from the one described by Calogero [10], in that Eqs. (8) and (9) do not require the definition of auxiliary basis functions, as is the case for the Calogero method, and the phase is obtained from the amplitude, while the reverse is the case for the Calogero method. Details are given in Appendix A. Eq. (8) has been solved non-iteratively in the past by using some form of a finite difference computational method, such as one of Milne’s predictor–corrector methods [26], or [8] by a Bulirsch–Stoer limit method [27], none of which will be used in the present study. The iterative method of Seaton and Peach [14] consists in rewriting Eq. (8) in the form k2 y4

=w+

1 d2 y

(10)

y dr 2

where

w(r ) = k2 − VT ,

(11)

and calculating the solution of Eq. (10) by means of the iteration [14]



k y2n+1

= w+

1 d2 yn yn dr 2

1/2

,

n = 0, 1, 2, . . . .

(12)

Here n denotes the order of the iteration, and the initial value of y is given by the WKB approximation [12,13] k y20

= w1/2 .

(13)

The advantage of formulating the iteration according to Eq. (12) is that y varies slowly with r, automatically approaching unity at

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

large distances, and hence (1/yn )d2 yn /dr 2 is small compared to w. Near the origin of r this term may become large. In that case the solution of Eq. (12) should be started at a point sufficiently far from the origin, where [(1/yn )d2 yn /dr 2 ]/w is sufficiently small compared to unity, depending on the accuracy desired, as will be described further below. A feature of the present Ph–A method is that no initial conditions are required to be imposed, and since the potential VT in Eq. (11) approaches zero asymptotically, the amplitude y automatically approaches unity, asymptotically. This last property is very important, since the wave function derived from the present Ph–A method automatically also has unit amplitude asymptotically, even though the calculation is not required to be performed to asymptotic distances. This property does not hold for finite element or finite difference calculations. Eq. (9) combined with the first order result is equivalent to the WKB approximation. The iteration scheme (12) provides a method to iteratively improve the WKB approximation, since after convergence, the resulting wave function is in much better agreement with benchmark wave functions than the WKB result for the numerical cases studied below.

35

Because Eq. (16) lacks the derivatives of Dn , it is to be used in the large distance radial regions where Dn ≪ w . Numerical examples will be given below. The iterations described in Eq. (16) can be replaced by an analytical method based on the solution of a cubic equation. The method is as follows. By adding on to both sides of Eq. (16) the function w , by defining zn = w + Dn ,

(17)

one finds zn+1 =

(5/16)(V ′ )2 zn2

+

(1/4)V ′′ zn

+ w,

n = 0, 1, 2, . . . .

(18)

4. A non-iterative analytical solution of Eq. (18) If the iterations converge, then for a very large value n∞ of n the values of the z’s will no longer change with n, i.e. zn∞ +1 ≃ zn∞ = z .

(19)

3. Computational method

Replacing in Eq. (18) all the zn ’s by z, and after some algebraic rearrangement one obtains the cubic equation

The spectral computational method consists in expanding a function y in the whole radial interval into a series of N + 1 Chebyshev polynomials Ts (x), s = 0, 1, 2, . . . , N,

z 3 = A + Bz + w z 2

y(x) =

N 

as Ts (x),

−1 ≤ x ≤ 1.

(14)

s=0

The corresponding N + 1 support points ξ1 , ξ2 , . . . , ξN +1 are the zeros of the polynomial TN +1 . Since the Chebyshev polynomials are defined in the interval −1 ≤ x ≤ 1 the quantities defined in the radial interval 0 ≤ r ≤ rmax are mapped into the x-variable by a linear transformation. The expansion cutoff value N is set arbitrarily, but once chosen, the location and number of support points on the x-axis (and correspondingly on the r-axis) is determined. In the iterative algorithm that follows, only the values of the function y evaluated at the N + 1 support points are required, hence y, and all the terms in Eq. (12) become vectors of dimension N + 1. Extensive use is made of the Clenshaw–Curtis matrix method (CC) [28] that relates the values of a function evaluated at the N + 1 mesh-points to the expansion coefficients as of that function, and vice-versa, by a simple known matrix [17] relation. The quantity (d2 yn /dr 2 )/yn occurring in Eq. (12) is denoted as Dn 2

Dn =

1 d yn yn dr 2

.

(15)

If the second order derivative were obtained by replacing the Ts in N Eq. (14) by their respective second derivatives, d2 y/dr 2 = s=0 as (d2 Ts (x)/dx2 )(dx/dr )2 , one would run into numerical difficulties [29] because the values of d2 Ts (x)/dx2 increase rapidly with the index s, and they would overwhelm the decrease of as with s. For example, for s = 16 and x = −1, d2 Ts (x)/dx2 = 2 × 104 . Thus the value of N would have to be kept very small, like 10, and accuracy would be lost in calculating the phase, Eq. (9). Instead, the derivatives of yn are calculated directly from their iterative expressions as follows, however other methods have also been developed [30]. By using yn+1 = k1/2 (w + Dn )(−1/4) , one obtains after some algebra Dn+1 = (5/16)(ω + Dn )−2 (ω′ + D′n ) − (1/4)(ω + Dn )−1 (ω′′ + D′′n ), where the primes mean derivatives with respect to r. By dropping all the derivatives of Dn , one obtains Dn+1 =

1 (V ′ )2 V ′′ + , 16 (w + Dn )2 4 (w + Dn ) 5

n = 0, 1, 2, . . . .

(16)

(20)

with A(r ) = (5/16)(V ′ )2

and B(r ) = (1/4)V ′′ .

(21)

This equation can be solved by standard methods as follows. By defining p = −(B + w 2 /3), q = −(2w 3 + 9w B + 27A)/27, and J = −4p3 − 27q2 ,

(22)

the three roots of Eq. (20) can be expressed in terms of p, q, w , and (−3J )1/2 . If J ≤ 0, so that (−3J )1/2 is real, then one of the roots z (1) is real, and the other two are complex conjugates of each other. Here z (1) = (A + B + w)/3

(23)

and where

A = −(27/2)q + (3/2)(−3J )1/2

(24)

and

B = −3p/A.

(25)

In the numerical examples given below, it is found that after 5 iterations, the value of ω + D5 differs from z (1) by less than 10−11 , showing that the iterations of Eq. (16) converge. In Eqs. (16) and (18) the derivatives of the potential with respect to r are required. If the potential V is given by analytic expressions, as is the case for the numerical examples described below, the derivatives can be calculated analytically. Otherwise the derivatives have to be calculated by some numerical interpolation procedure. 4.1. The steps of the algorithm 1. For the first iteration, D0 is set to zero in Eq. (12), and the values of the square root in Eq. (12) are evaluated at the support points {ξ }, and the values of y0 are then also obtained at the support points. That causes no non-linearity problems because these quantities can be manipulated numerically (square roots, inverses, etc.) support point by support point, without the Chebyshev expansion being invoked. 2. Subsequently the Chebyshev expansion coefficients an of y0 are obtained by means of the (CC) method, and in terms of

36

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

those the integral in Eq. (9) can be carried out by means of the Gauss–Chebyshev method [17,31]. The result is transformed back into the support points via the (CC) method, and the wave function defined via Eq. (1) is obtained at the support points. The result is the WKB approximation. 3. In order to obtain the next order function y1 the required second derivative of y0 is obtained via Eq. (16), with Dn replaced by D0 . 4. If iterations of Eq. (12) are performed, then Dn+1 is obtained in terms or Dn from Eq. (16), with n = 1, 2, . . . . Else the cubic equation (20) is solved analytically. The Chebyshev expansion of 1/y2n , Eq. (14) is invoked only in order to calculate the phase, Eq. (9), at the support points, and to subsequently interpolate the results to an equidistant mesh, as required for graphical purposes. 5. If a comparison is required of the Ph–A wave function with a wave function ψC calculated by some other means (‘‘C’’ stands for ‘‘numerically calculated’’), for instance in a short radial interval [0, r0 ], then in addition to the steps 1through 4, the procedure described in Section 5 is required in order to re-normalize ψC and to obtain the Ph–A starting phase φ0 at r0 . In the present numerical examples 5 iterations sufficed in order that the difference between yn+1 and yn was less than 10−8 for all support points in the whole radial interval. The difference between y5 and the non-iterative solution of Eq. (20) was also found to be less than 10−8 , which is an indication that the iterations converged for these cases. The integral in Eq. (9) required to calculate the phase φ is performed by a Gauss–Chebyshev method [17,31] that is well suited to this type of spectral expansion since it only requires the values of the expansion coefficients as . Situations that involve imaginary local wave numbers and the respective turning points, as is the case in the presence of repulsive barriers, or where the potential is everywhere repulsive, are postponed to a future study. 5. Connecting a conventional wave function to the Ph–A representation It is common that at small distances the potentials have a repulsive core, and the centripetal potential also introduces a strong repulsive part. For these cases, and also in the presence of repulsive barriers, it is preferable to calculate the wave function by conventional numerical means from the origin to a point r0 , and connect this wave function to the Ph–A wave function at r0 , so that the wave function can be propagated out to large distances by the phase–amplitude method. A procedure for implementing this connection will be described below. This procedure does not use the conventional matching procedure based on the Wronskian of the wave functions. If the wave function obtained by some numerical algorithm for 0 ≤ r ≤ r0 is denoted as ψC (C for ‘‘computational’’), and if that wave function is to be extended beyond r0 by means of the Ph–A representation for which φ0 and y0 are the initial phase and amplitude at point r0 , then two quantities need to be obtained from ψC , namely φ0 and the normalization factor ~ of ψC so that

~ψC = ψ = y sin(φ).

(26)

The wave function defined in Eq. (26) has the property that at large distances the amplitude goes to unity. But for long range potentials it is not possible to assure that ~ψC satisfies this requirement, unless ~ is obtained by a Wronskian matching procedure involving basis functions g and f that are known to approach unity at large distances. Such basis functions are known analytically for the Coulomb case, and for short ranged potentials, but are difficult to obtain for general long range potentials. The present procedure avoids the need for such basis functions. In the expansion of a wave function in the interval [r0 , rmax ] in terms of Chebyshev polynomials of the first kind, [20] the support

points {ξ } do not include r0 but have values r0 < ξ1 < ξ2 · · · < ξN +1 < rmax that are described below. The corresponding values of the amplitude at the support points are denoted as y1 , y2 , . . . and the phases are given by

φ1 = φ0 + ϕ1 ,

φ2 = φ0 + ϕ2 · · ·

(27)

where

ϕi = k



ri

[y(r ′ )]−2 dr ′ ,

i = 1, 2, . . . , N + 1.

(28)

r0

The values of yi and ϕi , i = 1, 2, . . . are known from the numerical procedure, Eq. (12) or Eq. (20). By defining

∆ = ϕ2 − ϕ1

(29)

and

O=

ψC (r2 )/y2 ψC (r1 )/y1

(30)

and after making use of some trigonometric identities, one obtains [30]

O − cos(∆) sin(∆)

= cot(φ0 + ϕ1 ).

(31)

The value of φ0 can be obtained from Eq. (31) and ~ is given by

~=

y1 sin(φ0 + ϕ1 )

ψC (r1 )

.

(32)

6. Results The feasibility of the present approach will be demonstrated by means of two examples, for which the potential VT is everywhere attractive. In these examples the potential has a long range tail proportional to either r −3 or r −1 respectively. In order to mitigate the singularities at the origin of r, a ‘‘rounding’’ procedure is introduced that replaces r by a constant in the limit r → 0 by means of an analytic transformation from r to R

R(r ) = r /[1 − exp(−r /t )],

(33)

but does not alter the potential at distances much larger than the rounding parameter t. 6.1. The 1/r 3 case In this case the potential is given by V3 (r ) = −1.6224 × 104 /R3 ,

(34)

where the ‘‘rounded’’ radial distance R(r ) is given by Eq. (33). The value of this potential, illustrated in Fig. 1, is appropriate for atomic physics applications [32], when the two interacting molecules have dipole moments each. The reason this 1/r 3 long range nature was chosen, is because this case did not get addressed successfully by means of a Born-approximation method [32], while it is well described in the present study. At r = 2500 the value of V is ≃10−6 , while near r = 0 its value is 1900. The latter is unphysically large, but the example serves to illustrate the numerical method described here. The corresponding wave function is highly oscillatory at small distances, with an amplitude and wave length that vary substantially with distance, as illustrated in Fig. 2. The corresponding phase shift θ (r ), defined as

θ (r ) = φ(r ) − (kr − Lπ /2),

(35)

is illustrated in Fig. 3. Because of the long range nature of the potential, the phase shift still changes for 1000 < r < 2000 by ≃0.3 rad.

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

Fig. 1. The ‘‘rounded’’ 1/r 3 potential, given by Eq. (34) with the rounding parameter t = 2. The units are in inverse length squared, since the potential, in energy units, has been multiplied by the factor 2m/ h¯ 2 . The dashed line indicates the non-rounded potential, with R replaced by r.

Fig. 2. The wave function calculated by the IEM method, corresponding to potential (34), with t = 10 in Eq. (33) for k = 0.01. Asymptotically its amplitude is unity. The dashed lines at top and bottom represent the amplitude y and −y, respectively, calculated from Eq. (12) for the same potential.

37

Fig. 4. Error of Ph–A wave function for the ‘‘rounded’’ 1/r 3 potential (34) with the rounding parameter t = 2, and k = 0.01. The Ph–A calculation uses 201 Chebyshev expansion functions in the radial domain [0, 2000]. The order n of the iteration is shown in parentheses in the legend. It should be noted that the error for the third iteration is larger than that for the first, even though beyond the third iteration the results do not change very much.

Fig. 5. Comparison between the iterative value of (ω + D5 ) obtained after the 5th iteration of Eq. (16) and the corresponding value z obtained from the analytic solution of Eq. (20). The figure shows the absolute value of the difference between the two for the potential used for Fig. 4.

his extensive book [10], as is elucidated in Appendix A. This is because the equations for each method are different from each other, although asymptotically the phase shifts must agree. The agreement of the Ph–A wave function with a benchmark wave function denoted as IEM, calculated by a spectral method [17] and normalized according to Section 5, is shown below. Noteworthy is the fact that in this example only 201 expansion terms in Eq. (14) have been used to calculate the amplitude and phase for the whole radial interval [0, 2000].

Fig. 3. The phase shift, defined in (35) for potential (34) with t = 2 for k = 0.01. Beyond r = 500 the phase shift still changes because the potential is not yet constant, as further explained in the text.

Their values are 207.45057 rad and 207.74413 rad for r = 1000 and 2000, respectively. The amplitude and phase functions obtained here are different from the ones examined by Calogero in

An evaluation of the error of the Ph–A wave function is obtained by plotting the absolute value of the difference of the Ph–A and the IEM wave functions. The result for the cases k = 0.01, calculated in the radial interval [0, 2000] with 201 Chebyshev expansion functions is illustrated in Fig. 4. The comparison between the iterative and non-iterative solutions is illustrated in Fig. 5. A calculation in a smaller radial interval, starting at r0 = 30, and ending at r = 200, using only 51 Chebyshev functions, is illustrated in Fig. 6. The initial value of the phase φ0 was obtained according to Section 5.

38

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

Fig. 6. Same as Fig. 4 in the radial interval [30, 200], using 51 Chebyshev functions. The equations in Section 4 are used to start the Ph–A wave function at r = 30. Please note the change in scale from Fig. 4.

Fig. 7. The phase θ as a function of distance for the Coulomb case described in the text.

6.2. The 1/r Coulomb case The conventional equation for the wave function χ (ρ) describing the Coulomb case is d2 χ dρ 2

 + 1−

2η

ρ

−

L( L + 1 )

ρ2



χ = 0,

(36)

where ρ describes the radial distance, and η describes the strength of the Coulomb interaction. The corresponding potential in Eq. (6) is given by V (r ) =

Ze2 2m r

}

2

=

Z¯ r

,

(37)

where Z¯ has units of inverse length and is related to 2η according to Z¯ = 2ηk

(38)

where 2η = Z

e2

}c



mc 2 E

1/2

.

(39)

In the above m is the mass and E the energy of the incident particle, Z the number of protons in the nucleus, e the charge of a proton, } is Planck’s constant, and ρ = kr. For the attractive case, Z and η are both negative. Thus, by setting k = 1 and Z¯ = 2η in Eq. (6), one obtains a solution to Eq. (36). Asymptotically the Coulomb wave function approaches

χ(ρ) → sin(ρ − η ln(2ρ) − Lπ /2 + θ∞ )

(40)

where θ∞ is the phase shift, whose value depends on how the potential differs from the point Coulomb potential. The position dependent phase is defined as

θ(ρ) = φ(ρ) − [ρ − η ln(2ρ) − Lπ /2] ,

(41)

which asymptotically approaches the phase shift θ∞ . An example for the Ph–A wave function for a ‘‘rounded’’ Coulomb potential with the rounding parameter t = 2, and η = −2, and started at point ρ = 30 has been calculated, and good agreement with the IEM wave function was obtained. The benchmark IEM wave function was needed only in the interval [0, 32], so as to determine the starting phase φ0 , according to the procedure described in Section 5. By means of this procedure the IEM wave function is also normalized, so that asymptotically its amplitude tends toward unity. The phase, given by Eq. (41), and illustrated in

Fig. 8. The accuracy of the Ph–A wave function for the rounded Coulomb potential, with rounding parameter t = 2 and for η = −1. The radial interval is [30, 1500], the number of Chebyshev polynomials used in this interval is 51. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7, still changes in the second significant figure for ρ ≃ 100, and it stabilizes in the fourth significant figure to −2.6399 rad beyond ρ ≃ 2500. That level of stabilization is compatible with the fact that at these distances the rounded Coulomb potential differs from the point Coulomb potential by 4 × 10−2 and 2 × 10−3 , respectively, and it demonstrates that the long range 1/ρ character is taken into account reliably. The accuracy of the Ph–A wave function for η = −1 is illustrated in Fig. 8. The first iteration gives an accuracy that is 104 times better than the WKB accuracy, and subsequent iterations do make only an additional slight improvement, and the iterations do converge. It is worth noting that the accuracy is practically independent of position for ρ > 500. This is a property of the Chebyshev expansion, for which the error distributes itself uniformly across the whole radial domain. For the case of a point Coulomb potential the iterative method does not converge for small values of ρ . Although the solution of Eq. (20) was obtained non-iteratively, it gave an error for the point Coulomb wave function of more than 10%. This can be understood by examining the value of D/ω, illustrated in Fig. 9, where ω is defined in Eq. (11) and D is defined in Eq. (15) with n → ∞. Since the amplitude y changes slowly with distance, the derivatives of D should be smaller than D itself, hence the plot of D/w gives an indication of the upper limit of the error of Eqs. (16) or (20), due to neglecting the derivatives of D.

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

39

Table 1 Computation times and Ph–A wave function accuracy.

Fig. 9. The functions D/w for point and rounded Coulomb potentials for η = −2.

N

Cubic time (s)

5 iterations time (s)

Accuracy

16 25 50 100

0.045 0.048 0.060 0.061

0.89 0.90 0.90 0.91

10−3 10−4 10−5 10−5

times given in Table 1 do not include the time to interpolate the y and φ to a fine equidistant radial mesh. Interpolating to an equispaced radial mesh size of step length h = 0.1 depends on the size of the radial interval. For the radial interval [0, 40] the fine mesh interpolation requires 0.8 s, and for the radial interval [40, 2000] the interpolation takes between 170 s and 180 s. However, the calculation of the slowly oscillating part M (S ) of an overlap matrix element (5) can be done by using the Gauss–Chebyshev integration method [31], which does not require the interpolation to an equispaced radial mesh, and is expected to take approximately 0.2 to 0.4 s for obtaining both the two wave functions and also M (S ) . 7. Overlap integrals The need for overlap matrix elements in quantum mechanics calculations is very common. An example is the Franck–Condon overlap matrix elements that describe the intensity of electronic transitions in atoms or molecules. [33] How to perform such a calculation, given by Eq. (2), by means of the Ph–A method will be presented below. The two wave functions ψ1 and ψ2 are solutions of the one-dimensional radial Schrödinger equation with the potential VM defined in Eqs. (42), (43), and (34), (33)

Fig. 10. Accuracy of the Ph–A wave function for the rounded Coulomb potential with η = −2 in the radial interval 0 < ρ < 200, using 101 Chebyshev polynomials.

Fig. 9 shows that the value of D/w for ρ < 1 is close to unity for the point Coulomb case while for the rounded Coulomb case it is close to 10−3 . Hence for the point Coulomb case the solution of Eq. (20) should not be expected to be accurate, while for the rounded Coulomb case, it should be expected to be accurate to at least 10−3 . The numerical calculations show that the accuracy is of the order of 10−5 , as illustrated in Fig. 10, which is much higher than expected. 6.3. Numerical details The calculations are done with MATLAB on a desk PC using an Intel TM2 Quad, with a CPU Q 9950, a frequency of 2.83 GHz, and a RAM of 8 GB. The calculation uses typically between N + 1 = 51 and 201 Chebyshev expansion polynomials. The computing time in the radial interval [30, 1500] and the corresponding accuracies for a Coulomb wave function calculation are given in Table 1. The first column lists the number minus one of Chebyshev polynomials used in the calculations. Column 2 gives the computing time for the non-iterative solution of Eq. (20). Column 3 gives the computing time for carrying out five iterations of Eq. (16), and column 4 gives the resulting accuracy. The computing time for the Ph–A iterations depends weakly on the number of Chebyshev functions N + 1, regardless of the size of the radial interval, and depends weakly on the value of k. By contrast, the IEM calculation in the radial interval [0, 1500] for the whole wave function depends on the value of k, requiring 0.20 s for k = 0.01 and 0.29 s for k = 0.1. The

VWS (r ) = −3.36 / [1 + exp{(r − 3.5)/0.6}]

(42)

VM = VWS + V3 .

(43)

However the rounding parameter t = 10 used in Eq. (33) is larger than the one used in previous examples, with the effect of extending to larger distances the rounding effect, and reducing the value of the potential V to ≃−20 near the origin. The wave numbers for each of the wave functions ψ1 and ψ2 are k = 0.01 and 0.005, respectively (in units of inverse length). The two wave functions have different amplitudes, but nearly the same phases at distances where |V | > k2 , as can be seen by the near coincidence between the locations of the maxima and minima in Fig. 11. The overlap potential U is taken from Eq. (4) of Ref. [8], and represents the screened interaction of an electron with an ion embedded in a plasma. It is composed of a sum of exponentials divided by the radial distance r, and is illustrated in Fig. 12. It has a 1/r singularity at r → 0. Using the Ph–A representation of ψ1 and ψ2 , the integrand of the overlap integral separates into a fast oscillating and slowly oscillating parts, Eq. (5), as described above. These integrands are illustrated in Fig. 13. The approximate values of M (F ) and M (S ) are −0.073 and 0.258. As expected, the integrand of M (F ) is more oscillatory than the integrand of M (S ) , and hence |M (F ) | < |M (S ) |. Accordingly a crude estimate of M is given by M (S ) , which can be calculated directly within the Ph–A spectral representation, without the necessity to interpolate to small radial meshes. 8. Summary and conclusions A spectral method of solving the Phase–Amplitude nonlinear equations of Milne [1] for the case of attractive potentials has been implemented. The method is based on the iterative solutions of an approximation to the equations described by Seaton and Peach [14], and carried out by means of a spectral expansion in terms of Chebyshev polynomials [15–17]. At large distances this

40

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

Fig. 11. The two wave functions used in the calculation of the matrix M, defined in Eq. (2). Both have unit amplitude at r = 2500.

Fig. 12. The overlap function U that occurs in the matrix element M, defined in Eq. (2). Because of the factor 1/r, it becomes ∞ at r = 0. The units of U are (1/length)2 .

directly. For example, with 100 Chebyshev basis functions it produces wave functions that, after the third iteration, are accurate to between 0.1% and 0.001% out to distances as large a 2000 units of distance. The new results obtained in the present study are as follows: (a) A wave function calculated only for small distances by some conventional method can be normalized by comparison with a Ph–A wave function, such that asymptotically it has unit amplitude. Such normalizations are usually carried out by matching the computed wave function to analytical basis functions whose asymptotic normalization is known. But for potentials of the form 1/r 3 , for example, there are no analytical functions available for matching at a nearby point. This difficulty is now removed. (b) The accuracy of the Phase–Amplitude method is studied. To the knowledge of the author, no other accuracy investigation of this Ph–A method has been presented in the literature before. (c) A comparison between Milne’s and Calogero’s Ph–A methods is made, as described in Appendix A. The two are found to be quite different. (d) The method can calculate Coulomb wave functions for the attractive case without requiring to calculate complicated Hypergeometric confluent functions. (Codes exist, but it is convenient not to have to depend on them). (e) It is shown in Appendix B that the numerical error accumulation of the Ph–A method is uniform in the whole radial interval, while for a finite element method the accumulation is exponential with the number of elements (or partitions). Hence for large distances the Ph–A wave function is more accurate than the wave functions obtained by conventional means, and further, the calculation is substantially faster. In summary, the Ph–A method of Refs. [1,14], in combination with a spectral expansion method [15–18] provides a robust and economical method for propagating a wave function out to large distances for situations where the conventional solutions of the Schrödinger equation may be inadequate. This method is also expected to be useful for the calculation of overlap matrix elements that involve highly oscillatory wave functions, and to provide a very economical method to store wave functions. The present results open the way to generalize the Ph–A method to scattering cases where barriers are present, to bound states, or to the situation of coupled channel equations for which the potentials involved are of long range. Acknowledgments The author is indebted to Dr. Ionel Simbotin for calling attention to the Ph–A representation, and for stimulating conversations, and also to Professors Richard Pratt and Daniel H. Gebremedhin for useful correspondence. Also, for the record, the author is indebted to Professor G. Breit for having introduced the author to the paper of Seaton and Peach in the middle 1960s. Appendix A. Comparison with the Calogero’s Ph–A method

Fig. 13. The integrands of the matrix elements M (F ) and M (S ) . Due to the presence of oscillations of the integrand of M (F ) , it is clear that M (F ) < M (S ) .

approximation, given by Eq. (16), tends toward the original equation (12) [14], and at short distances the whole procedure can be replaced by a conventional solution of the Schrödinger equation. The spectral method substantially improves the accuracy of the solution for large distances in that it is very accurate and requires far fewer mesh points than methods that calculate the wave function

It will be shown that the Phase–Amplitude equations, denoted as PAM , described by Calogero [10], and used by several authors, including Fano [9], is not the same as the one described here, based on the equations of Milne [1]. The start of the PAM equations is [9] u(r ) = α(r )[f (r ) cos δ(r ) − g (r ) sin δ(r )],

(A.1)

with the constraint du(r ) dr

= α(r )



df dr

cos δ(r ) −

dg dr



sin δ(r ) ,

(A.2)

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

41

where u is the solution of a radial Schrödinger equation for potential V (r ) = U0 (r ) − U (r ), U0 is an ad hoc ‘‘comparison potential’’ f and g are regular and irregular solutions of the Schrödinger equation for the potential U0 , with Wronskian W . The resulting nonlinear equation for the phase δ(r ) is

δ(r ) = W

−1



2m

r



h¯ 2

U (r ′ )[f (r ′ ) cos δ(r ′ ) 0

− g (r ′ ) sin δ(r ′ )]2 dr ′ .

(A.3)

Once Eq. (A.3) is solved for δ,the amplitude α(r ) can be obtained by an integral from r to ∞ involving f , g , δ , and U. By comparing Eq. (1) with Eq. (A.1), it is clear that φ(r ) ̸= δ(r ), and y(r ) ̸= α(r ), since no comparison functions f and g are required for Eq. (1). A phase–amplitude combination will now be shown that is totally different from the Milne-based combination, and which is similar to that of Calogero. The latter can be obtained from equations similar to Eqs. (A.1) and (A.2), for which the functions f and g are given by sin(kr ) and cos(kr ), respectively. These are solutions of a Green’s function Fredholm integral equation [17] that is equivalent to the solutions of Eq. (6) for potential VT .

ψ(r ) = f (r )A(r ) + g (r )B(r )

Fig. 14. The phase and amplitude calculated from Eqs. (A.9) and (A.11) for a numerical example described in the text. These results simulate the ones of the Calogero method. The lowest of the three curves represents the phase. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(A.4)

with A(r ) = 1 −

1

∞



k

g (r ′ )VT (r ′ )ψ(r ′ )dr ′

(A.5)

r

and B(r ) = −

1

r



k

f (r ′ )VT (r ′ )ψ(r ′ )dr ′ ,

(A.6)

0

as described in Ref. [17]. From Eq. (A.4) one finds dψ(r ) dr

=

df (r ) dr

A(r ) +

dg (r ) dr

B(r ).

(A.7)

Here ψ is the same as u above, and VT is the same as U above. By defining a phase shift Θ (r ) cos[Θ (r )] =

A(r ) Y (r )

,

sin[Θ (r )] =

B(r ) Y (r )

(A.8)

with Y (r ) = A2 (r ) + B2 (r )



1/2

,

(A.9)

and remembering that f and g are given by sin(kr ) and cos(kr ), respectively, one finds

ψ(r ) = Y (r ) sin[Φ (r )],

Fig. 15. Comparison between the phase shifts Θ (r ) for the Calogero and θ(r ) for the Milne methods, respectively, labeled ‘‘PAM’’ and ‘‘PhA’’, and given by Eq. (A.11) and Eq. (35), respectively. The potential is V3 described in Section 6, with t = 2 and k = 0.01. In the region where Calogero’s phase is nearly constant, the corresponding value of sin(Φ ) is near zero so as to compensate for the value of the amplitude y larger than unity.

(A.10)

with

Φ (r ) = kr + Θ (r ).

(A.11)

Although formally Eq. (A.10) is identical to Eq. (1), the values of Y and Φ are not the same as y and φ of Milne’s equation because the equations these quantities obey are not the same, even though the resulting wave functions ψ are the same. An example is displayed in Fig. 14. The potential VM is the same as described in Section 7, and the wave number is k = 0.01. Another example is given for the potential V3 described in Section 6.1, for k = 0.01. The results for the phases, obtained either by the Milne or the Calogero methods are compared in Fig. 15 and the respective amplitudes are compared in Fig. 16. For the Calogero method the amplitude Y is not monotonic, the maxima are larger than unity, and the phase shift Θ increases with distance in a ‘‘stairway’’ fashion. In the region where Y > 1, sin(Φ ) becomes sufficiently small so that the product yields the correct value for the wave function.

Fig. 16. Comparison between Calogero’s and Milne’s methods for the amplitudes Y (solid line) and y (dashed line), respectively. The potential and wave number are the same as in Fig. 15.

42

G. Rawitscher / Computer Physics Communications 191 (2015) 33–42

a factor of 10 or 20. If instead of 20 partitions, 200 partitions are needed to propagate the wave function to a distance larger than r = 200, the total error factor would be E (200) = [E (20) ]10 ≃ 1010 .

(B.3) −11

Starting with an error of order 10 , the resulting error after 200 partitions would become unacceptably large. References [1] [2] [3] [4] [5] [6] [7] Fig. 17. The accuracy of the FE–DVR wave function for the potential VM as obtained by comparison with the IEM result. The latter is accurate to 10−11 . The wave number is k = 0.5 fm−1 , the number of Lobatto points per partition is 20, the size of each partition is 1 fm.

[8] [9] [10] [11]

Appendix B. Estimate of errors of the FE–DVR method This estimate is based on Ref. [21] entitled ‘‘Accuracy of a new hybrid finite element method for solving a scattering Schrödinger equation’’. The acronym FE denotes ‘‘finite element’’ and DVR denotes ‘‘Discrete Variable Representation’’. The latter uses Lagrange polynomials as basis functions in each element (or, partition), and is considered to be very accurate. Fig. 4 of that paper, reproduced below, shows that the accuracy of a particular FE–DVR calculation decreases exponentially with distance. An explanation for the exponential increase of the error is as follows. The transition of the wave function ψ I from partition I to the adjoining partition I + 1, which provides continuity in the value and derivative of the wave function, can be expressed as

ψ

(I +1)

(I )

= O ψ (1 + ε),

[13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

(B.1)

where O is an operator acting on the expansion coefficients of the wave function in ψ I +1 , based on the value and derivative of the wave function ψ I , evaluated at the end of partition I. This operator is described in Eq. (23) of Ref. [21]. The quantity ε describes the error that is incurred when numerically performing the operation O in the transition from partition I to I + 1. This error is composed of round-off errors, truncation errors in the expansion of the wave function in terms of basis functions in each partition, [34–37], and also depends on the size of the condition number associated with the two matrices involved in Eq. (23). For each new partition an additional factor (1 + ε) applies, so the overall accumulated error E (N ) after N partitions is E (N ) = (1 + ε)N ≃ eN ε .

[12]

(B.2)

In Fig. 17 the loss of accuracy after 20 partitions is approximately

[23] [24] [25] [26] [27]

[28] [29] [30] [31] [32] [33] [34] [35] [36] [37]

W.E. Milne, Phys. Rev. 35 (1930) 863. M.J. Korsch, H. Laurent, J. Phys. B: At. Mol. Phys. 14 (1981) 4213. F. Calogero, D.G. Ravenhall, Nuovo Cimento 32 (1964) 1755. F. Robicheaux, U. Fano, M. Cavagnero, D.A. Harmin, Phys. Rev. A 35 (1987) 3619. J.L. Dehmer, U. Fano, Phys. Rev. A 2 (1970) 304. C.H. Greene, A.R.P. Rau, U. Fano, Phys. Rev. A 26 (1982) 2441. B. Wilson, C. Iglesias, Mau Chen, J. Quant. Spectrosc. Radiat. Transf. 81 (2003) 499. A.B. Ritchie, A.K. Bhatia, Phys. Rev. E 69 (2004) 035402(R). U. Fano, C.E. Theodosiou, J.L. Dehmer, Rev. Modern Phys. 49 (1976) 49–68. F. Calogero, Variable Phase Approach to Potential Scattering, Academic Press, NY, 1967. Lv Wentao, Shen Chen, Fan Gui, Tian Zhongshan, Jiang Daozhuo, Real-time spectrum analyzer based on all phase FFT spectrum analysis, in: 2013 Fourth International Conference on Digital Manufacturing & Automation, ICDMA, IEEE, Qingdao, China, Piscataway, NJ, USA, 2013, pp. 29–30. H. Jeffreys, B.S. Jeffreys, Methods of Mathematical Physics, Cambridge University Press, NY, 1966. H.A. Kramers, Z., Physik 39 (1926) 828. M.J. Seaton, G. Peach, Proc. Phys. Soc. 79 (1962) 1296. http://dx.doi.org/10. 1088/0370-1328/79/6/127. A. Deloff, Ann. Phys. (NY) 322 (2007) 1373–1419. L.N. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia, PA, 2000. R.A. Gonzales, J. Eisert, I. Koltracht, M. Neumann, G. Rawitscher, J. of Comput. Phys. 134 (1997) 134. R.A. Gonzales, S.-Y. Kang, I. Koltracht, G. Rawitscher, J. Comput. Phys. 153 (1999) 160–202. Y.L. Luke, Mathematical Functions and their Approximations, Academic Press, NY, 1975. John P. Boyd, Chebyshev and Fourier Spectral Methods, second revised ed., Dover Publications, Mineola, NY, 2001. J. Power, G. Rawitscher, Phys. Rev. E 86 (2012) 066707. V.A. Dzuba, A. Derevianko, J. Phys. B: Atomic, Molecular and Optical Physics 43 (2010) 074011. V.A. Dzuba, A. Derevianko, V.V. Flambaum, Phys. Rev. A 86 (2012) 054501. S.G. Porsev, A. Derevianko, Phys. Rev. A 74 (2006) 020502. M.S. Safronova, S.G. Porsev, U.I. Safronova, M.G. Kozlov, C.W. Clark, Phys. Rev. A 87 (2013) 012509. M. Abramowitz, I. Stegun (Eds.), Handbook of Mathematical Functions, Dover, 1972, Eq. 25.5.13. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Section 17.3. Richardson extrapolation and the Bulirsch-Stoer method, in: Numerical Recipes: The Art of Scientific Computing, third ed., Cambridge University Press, New York, ISBN: 978-0-521-88068-8, 2007. C.C. Clenshaw, A.R. Curtis, Numer. Math. 2 (1960) 197. G. Rawitscher, 2014. arXiv:1403.5186. A. Bar-Shalom, M. Klapisch, J. Oreg, Comput. Phys. Commun. 93 (1996) 21. G. Rawitscher, I. Koltracht, Comput. Sci. Eng. 7 (2005) 58. G. Rawitscher, Phys. Rev. A 87 (2013) 032708. E. Condon, Phys. Rev. 28 (1926) 1182–1201. D. Baye, Phys. Status Solidi 243 (2006) 1095–1109. M.J. Rayson, Phys. Rev. E 76 (2007) 026704. Hua Wei, J. Chem. Phys. 106 (1997) 6885. B.I. Schneider, N. Nygaard, Phys. Rev. E 70 (2004) 056706.