Research: Science and Education
Accurate and Efficient Perturbation Theory by Matrix Inversion B. Cameron Reed Department of Physics, Alma College, Alma, MI 48801;
[email protected] Given that perturbation methods persist in their value for securing approximate solutions to the Schrödinger equation in cases where analytic solutions are intractable, discussion of such methods remains a popular topic in educational journals. In a recent paper in this Journal, for example, Besalu and Carbo-Dorca (1) describe a recursive matrix-based algorithm for nondegenerate Rayleigh–Schrödinger (RS) perturbation theory (PT), developing explicit expressions for energy eigenvalues up to sixth order. The physics pedagogical literature is particularly rich in discussions of PT methods. In the present context, for example, Gottdiener’s derivation (2) of first-, second-, and third-order RS energy expressions via examination of the secular equation arising from the determinant of a matrix whose elements involve the perturbation integrals ∫ ψ i(* ψ j )dτ (where the ψ ’s are unperturbed wave functions, * is the perturbing Hamiltonian, and dτ is the relevant volume element) is of interest, as is Lynch and Mavromatis’s (3) demonstration of how a variational form of RS perturbation theory can be carried out using a Mathcad electronic scratch pad. As presented in many papers and numerous standard texts (4 ), the central concept underlying PT is that both the eigenvalues and eigenfunctions corresponding to some potential can be expanded in powers of a perturbation parameter, usually designated as λ. Since λ has no physical meaning, it ultimately cancels out of the calculations, leaving a series of successive approximations to the perturbed eigenstates accurate to a given order of (the vanished) λ if the perturbation is small compared to the spacing between the energy levels of the unperturbed states. For students, this traditional approach to PT has two main disadvantages: the introduction of a parameter that has no physical role in the theory seems arbitrary at best, and the computations attending energy corrections beyond the first order become so involved as to preclude any but simple examples. The situation becomes even more complex on considering degenerate systems. The purpose of this paper is to describe an approach to degenerate time-independent PT that has a number of advantages over the traditional approach: there is no need to introduce a perturbation parameter or various orders of perturbation, the perturbation can be of any strength, any desired number of unperturbed states can be made to contribute to a perturbed state, the convergence of the calculation to an accurately defined perturbed state is readily apparent, and the method is easy and efficient to implement. No extensive summations are involved, and no integrals beyond those normally faced in first-order PT need be computed; the only limitation is the ability to invert a matrix. While the method presented here incorporates elements of all the treatments referenced above, it differs from them in significant ways. Like that of Besalu and Carbo-Dorca, my approach is fundamentally a matrix one; however, one does not have to keep track of orders of perturbation, and 1444
the method is applicable to degenerate and nondegenerate systems alike. Like Gottdiener, I make use of a secular equation, but again without the complication of tracking orders. As do Lynch and Mavromatis, I depend on the computational capabilities of desktop computers, but without invoking their variational approach. Since only standard perturbation integrals are used, the formulation of this method is extremely transparent. Perturbation Theory Up to eq 7 below, the development that follows is a form of standard textbook perturbation theory likely to be familiar to many readers. However, it is worthwhile summarizing this briefly in the interests of defining notation and of achieving a self-contained presentation. Assume as usual that one has a set of solutions to the Schrödinger equation for some potential V (0)(rW): wave functions ψn(0)( rW), and energies E n(0). While these unperturbed states may be infinite in number, a practical computation must be limited to using only a finite number N of them. Therefore, designate each of the ψn(0)( rW )—whether or not some of them be degenerate—with a unique integer subscript n = 1, 2, 3, …, N; this numbering can proceed in any arbitrary way. A perturbation V ′( rW) is then introduced, rendering the potential as V (rW) =V (0)(rW) + V ′(rW)
(1)
The boundary conditions for the perturbed and unperturbed potentials are assumed to be the same. A general perturbed state is then designated as ψm( rW) (m < N; energy Em) and, as usual, we assume that it can be expressed as a linear expansion in terms of the unperturbed wave functions: N
ψm = Σ C n mψn(0) n=1
(2)
Substituting eq 2 into the Schrödinger equation gives
Σ Cn mE n(0)ψn(0) + Σ Cn mV′ψn(0) = E m Σ Cn mψn(0)
(3)
It is worthwhile remarking that ψm( rW) represents a perturbed wave function and should not be thought of as “the perturbed mth state.” With no numbering system assigned to perturbed states, such a concept is meaningless. Multiplying through eq 3 by one of the unperturbed wave functions, say ψj(0) ( j = 1, 2, …, N ), and integrating over the range of the potential gives
Σ Cn m
E n – E m ∫ ψj(0) ψn(0)dτ + ∫ V ′ψj(0) ψn(0)dτ = 0 ( j = 1, N ) (4) (0)
Since no orthogonality properties of the unperturbed wave functions have been invoked, eq 4 is general for either de-
Journal of Chemical Education • Vol. 76 No. 10 October 1999 • JChemEd.chem.wisc.edu
Research: Science and Education
generate or nondegenerate systems. On defining an overlap integral and
ψjn = ∫ ψ j(0) ψn(0)dτ
(5)
Vjn = (E n(0) – Em) ψjn + ∫V ′ψ j(0) ψ n(0)dτ
(6)
eq 4, when written out for all possible values of j = 1, 2, 3, …, N, can be put in a matrix form:
V 11
V 1k
V 1r
V 1N
C 1m
V k1
V kk
V kr
V kN
C km =0
V r1
V rk
V rr
V rN
C rm
V N1
V Nk
V Nr
V NN
C Nm
(7)
The foregoing development has essentially reiterated textbook perturbation theory. Nondegenerate first-order perturbation theory, for example, is equivalent to restricting attention to the diagonal elements of the N by N square matrix (hereafter designated [V]) appearing in eq 7. To proceed further, particularly in the case of degenerate systems, it is customary to isolate that submatrix of eq 7 that incorporates only states of some energy E ; diagonalizing the submatrix then leads to various possible Em for perturbed states that derive from unperturbed states of energy E. However, this procedure unnecessarily restricts the calculations to first-order effects. The approach taken here is that if the matrix elements Vjn can be computed in general, why not use information from all of them (or at least N by N of them) in determining the E m? This represents the point of departure of the present method from the traditional one. Equation 7 represents a set of homogeneous linear equations. The only way of obtaining nonzero solutions for the Cjm is that the determinant of [V] be zero, which suggests a strategy for determining the perturbed energies: by trial and error, search for those values of E m that make det [V] = 0; this is essentially the approach taken by Gottdiener in ref 2. In practice, since det [V] can become large, the easiest way to proceed is to isolate values of E m between which the determinant changes sign. There is no need to isolate submatrices or to keep track of orders of perturbation; indeed, the convergence of a putative Em to a stable value can be tracked by monitoring the effect of increasing N, as is illustrated in example 1 below. To establish the Cjm , it is helpful to split up one of the elements of [V] with the help of eq 6, say element Vrk, and cast eq 7 as V 11
V 1k
V 1r
V 1N
C 1m
V k1
V kk
V kr
V kN
C km
0
Note that only the (r,k) element of [V] has been modified. Retaining the designation of [V] for the (modified) square matrix on the left side of eq 8, a general expansion coefficient Cjm can be obtained in terms of Ckm upon inverting [V]: {1
(0)
C j m = C k m Vj r E m – E k
ψr k
(9)
where Vjr{1 designates the ( j,r) element of the inverse of [V]. Consider the special case j = k in eq 9. Provided that the matrix in eq 8 can be inverted (i.e., that element (r,k) is in fact changed by splitting it up, rendering the determinant of [V] nonzero) and that unperturbed state k makes a contribution to ψm( rW ), then one should find Vkr{1(Em – Ek(0))ψrk = 1. In other words, a “viable” set of coefficients Cjm can always be found by selecting an element (r,k) in eq 8 that yields Vkr{1(Em – Ek(0))ψ rk = 1. Many choices of (r,k) will in fact satisfy this situation, since many unperturbed wave functions will contribute to the perturbed wave function, leading to different possible sets of coefficients Cjm for a given ψm( rW ). However, upon normalization, the same computed ψm(rW) will result, with the possible exception of a physically unimportant multiplicative factor of {1. Once the ratios C jm /C km are on hand from eq 9, normalization of the perturbed wave function proceeds by determining Ckm according as N
N
Σ C n mψn(0) pΣ=1 C p mψp(0) n =1
ψm2 dτ = 1 ⇒
dτ = 1
Using eq 9 to express Cnm and Cpm in terms of Ckm, this constraint reduces to
1
Ck m =
N
N
ΣΣ n =1 p =1
(0)
E m – E k ψr k
{1
{1
V n r V p r ψn p
where Vnr{1 and Vpr{1 again denote the (n,r) and ( p,r) elements of the inverse of the matrix in eq 8. Back-substituting this result into eq 9 gives {1
Cj m =
Vj r N
N
ΣΣ
n =1 p =1
{1 Vn r
(10) {1 Vp r
ψn p
Examples
1. One-Dimensional Harmonic Oscillator As a first test of the above procedure it is desirable to consider an example whose solution is known exactly. Such a case is provided by a harmonic oscillator perturbed by an additive linear potential: V(x) = 1 kx + β 2E 0 k x 2 2
0
(8)
=
V r1
(0) ∫ V′ψ(0) r ψk dτ
V rr
V rN
C rm
C km E m – E k ψrk
V N1
V Nk
V Nr
V NN
C Nm
0
(0)
(11)
where β is a dimensionless strength parameter and E0 is the ground-state energy of the unperturbed oscillator, h⁄ ω/2. β = 1 corresponds to a perturbation of 2E0 at the turning points of the unperturbed ground state. First-order perturbation theory predicts no energy shift for this situation, whereas second-order theory predicts the exactly correct result. By
JChemEd.chem.wisc.edu • Vol. 76 No. 10 October 1999 • Journal of Chemical Education
1445
Research: Science and Education
completing the square in eq 11 it can be shown that the Schrödinger equation for this potential reduces to that of a harmonic oscillator shifted to negative x and depressed in energy; the potential is symmetric about x = { β (2E0 /k)1/ 2 and leads to even- (odd-)parity solutions that are symmetric (antisymmetric) about that position. In terms of a dimensionless coordinate χ given by χ = (x/xo) where xo = (2E0/k)1/ 2, the perturbed wave functions and energies are given by
2. Two-Dimensional Isotropic Harmonic Oscillator As an example of a degenerate system, consider a twodimensional isotropic harmonic oscillator subject to an asymmetric perturbation:
ψnperturbed( χ) = AnHn(χ + β) exp{{( χ2/2 + βχ + β2/2)} (12)
where β is again a dimensionless strength parameter and E0 is the ground state energy h⁄ ω (note difference from previous example). The unperturbed wave functions for this case involve two quantum numbers, nx and ny (= 0, 1, 2, ...), and are given by
and E n = [2n + 1 – β2]E0
(n = 0, 1, 2, ...)
(13)
(π1/ 2 2n n!){1/2
and where Hn(χ + β) designates the where A n = nth order Hermite polynomial of argument (χ + β). The matrix elements in this case are given by
H
k
3/2
E0
x 2y
ψ n(0),n (x, y) = x
j = n – 1, n ≠ 0 j=n j=n+1
β E 0[2n!/(n – 1)!]1/ 2
(En(0) – Em) (14) 1/ 2 β E 0[2(n + 1)!/n!] or by zero otherwise. Figure 1 shows the rapidly convergent behavior of Em as a function of the number of contributing unperturbed states N (= dimension of [V]) for the case of the first excited state of such a linearly perturbed oscillator with β = 0.7; the exact result (eq 13) is E = 2.51 E0. At N = 4, the difference between the computed and exact perturbed energies is less than 1%; at N = 7 (i.e., upon incorporating contributions from up to and including the sixth excited state), Em/E0 = 2.5100173. Figure 2 shows the similarly rapid convergence of the computed normalized wave functions to the exact solution; note the symmetry of the latter about χ = { β = {0.7, as alluded to above (E0 is taken here to be unity). Further experimentation with this example shows that the present method successfully recovers perturbed states even in cases where β is large enough to place a perturbed state at an energy lower than that of unperturbed levels lying below the unperturbed state corresponding to the state being investigated. Vjn =
V (x, y) = 1 k x 2 + y 2 + β 2
α π
y
2 2 2 2 1 H n αx e {α x /2 H n αy e {α y /2 n x+ny x y nx !ny! 2
(15)
(16)
where α4 = (mk/h⁄ 2) and where H j(αx) again denotes the Hermite polynomial of order j and argument αx. The unperturbed energies are integral multiples of E0 according as
E
(0) n x ,ny
= nx + ny + 1 E0
(17)
The degeneracy of energy E ((0)nx,ny) is (nx + ny + 1). The perturbation integrals ∫V ′ ψ n(0)ψ(0) j d τ are double integrals over dτ = dx dy and involve four quantum numbers: (nx, ny) from state n, and (jx, jy) from state j. They evaluate to n x+n y+jx+jy ∫∫ V ′ ψ n(0)ψ(0) nx !ny !jx!jy!}{1/2I x I y (18) j dx dy = β E 0{2
where Ix = {2nx–1(2nx + 1)nx! δnjxx + 2nx(nx + 2)! δ njxx+2 + 2nx–2nx! δ njxx–2} (19) and Iy = {2ny–1ny! δ njyy –1 + 2ny(ny + 1)! δ njyy +1}
(20)
3.6
3.4
E(m) / E(0)
3.2
3.0
2.8
2.6
2.4
1
2
3
4
5
6
7
8
N Figure 1. Computed perturbed energy as a function of the number of contributing unperturbed states N for the first excited state of a linearly perturbed harmonic oscillator with β = 0.7 (see eq 11). The curve is interpolated. The exact solution has E = 2.51 E0.
1446
Figure 2. Computed and exact normalized wave functions for the first excited state of a linearly perturbed harmonic oscillator with β = 0.7. Computed curves with N > 4 are essentially indistinguishable from the exact solution.
Journal of Chemical Education • Vol. 76 No. 10 October 1999 • JChemEd.chem.wisc.edu
Research: Science and Education
The degenerate unperturbed levels are Table 1. Lowest 6 split by the perturbation. For α = 1, Energy States of Perturbed 2-D Isotropic β = 0.5, and N = 10 (that is, inHarmonic Oscillator corporating contributions from up m Em /E0 to and including the unperturbed ninth excited state, a total of 55 indi1 0.9376297 vidual states), the lowest six per2 1.6038944 turbed energies are listed in Table 3 2.3997905 1. The splitting of the first and sec4 2.5775212 ond excited states is nearly symmet5 2.9665967 ric. Figure 3 shows a pseudo-three6 3.4464876 dimensional view of the perturbed NOTE: β = 0.5. ground state for { 6 ≤ (αx, αy) ≤ + 6; consistent with the nature of the perturbing potential, the perturbed wave function is symmetric in x and asymmetric in y. Aside from the transparency of its formulation, the present method is easy to execute on even a modest desktop computer. Any specific problem is best approached by writing a routine to compute the matrix elements Vij and reading them into general routines for handling the computation of det[V] and inverting eq 8 upon inputting trial estimates for E m. Fast, compact determinant and matrix-inversion routines can be found in a number of sources (e.g., Press et al. [5]). Summary A conceptually and computationally straightforward approach to PT has been developed that has significant advantages over traditional methods: no unphysical perturbation parameter need be introduced, the user need not be concerned with orders of perturbation or submatrices, the perturbation can be of any desired strength, the effect of any desired number of unperturbed states can be incorporated in the calculation of a perturbed state, the convergence of the calculation to an accurately defined perturbed state is easily tracked, and no integrals beyond those encountered in traditional first-order PT need be computed. Once working determinant and inversion routines are on hand, the difficult work for any particular problem is reduced to establishing the matrix elements Vij . Finally, a few words on the general approach taken in this work are appropriate. Since the approach here is to treat a full Hamiltonian eigenvalue problem with a finite basis set
Figure 3. Ground state wave function for a perturbed two-dimensional isotropic harmonic oscillator subject to a perturbing potential proportional to x2y; β = 0.5 in eq 15.
assumed to span the function space containing the exact eigenfunctions of the perturbed Hamiltonian, the method is not identical to perturbation theory in the usual sense of picking off various orders of approximation from infinite expansions. For pedagogical purposes, however, this shortcoming is unimportant and arguably is more than compensated for by the intuitively powerful idea of approximating more and more accurate solutions by expanding the matrix in eq 7. Acknowledgments I am grateful to Gene Deci and Alan Jackson for a number of comments. Literature Cited 1. 2. 3. 4.
Besalu, E.; Carbo–Dorca, R. J. Chem. Educ. 1998, 75, 502. Gottdiener, L. Am. J. Phys. 1978, 46, 893. Lynch, R.; Mavromatis, H. A. Am. J. Phys. 1991, 59, 270. Liboff, R. L. Introductory Quantum Mechanics, 2nd ed.; AddisonWesley: Reading, MA, 1992; Chapter 13. 5. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: New York, 1986; Chapter 2.
JChemEd.chem.wisc.edu • Vol. 76 No. 10 October 1999 • Journal of Chemical Education
1447
