The Diatomic Anharmonic Oscillator according to Matrix Mechanics
by M. B. Monagan and J. F. Ogilvie, Simon Fraser University, Centre for Experimental and Constructive Mathematics, Vancouver, Canada
NOTE: The calculation undertaken here employs perturbation theory of multiple orders; the extent of the calculation, according to the order, is adjustable by a user within limitations of machine capacity and acceptable duration of calculation.
Introduction
Matrix mechanics is the original form of quantum mechanics that Werner Heisenberg initiated in 1925 and elaborated in association with Max Born and Pascual Jordan. By this means one can solve physical problems, such as that of an anharmonic oscillator that usefully models vibrational properties of a diatomic molecule as we demonstrate here. The full calculation is practicable with Maple statements and procedures that follow, and an augmented description and explanation appears in a paper to be published during 1999. Alternative procedures in Maple to generate energies of vibration-rotational states of a diatomic molecule, based on BKW quasi-classical theory and hypervirial perturbation theory in wave mechanics, are available in separate worksheets.
The calculation undertaken here employs perturbation theory of multiple orders; the extent of the calculation, according to the order, is adjustable by a user within limitations of machine capacity and acceptable duration of calculation.
> restart:
Computations
To calculate terms of perturbative order K in series for matrices E and U , we proceed as follows. First we construct a matrix to a given rank. Next we create a matrix equation for the coefficient of by expanding an equation symbolically as a power series in . We insert all known matrices truncated to a specified rank and solve this equation for and off-diagonal entries of . Next, to solve for diagonal entries of , we make use of an equation . We expand it as a power series in symbolically and select the coefficient of . We insert all known quantities calculated up to this point and solve for diagonal entries of . After solving for and , one can proceed to calculate terms of order .
Setting up Maple
The following statement tells Maple that is a scalar quantity.
> constants := constants, lambda;
The following Maple procedure is used to simplify matrices of formulae; it factors out any common factor present in each entry of a matrix.
> factormatrix := proc(A) local multigcd, B, L, N, D, G; multigcd := proc(a) local i,g; g := a[1]; for i from 2 to nops(a) while g <> 1 do g := frontend(gcd,[g,a[i]]); od; g end; B := map(factor,A); L := {op(map(op,convert(B,listlist)))} minus {0}; if nops(L) = 0 then RETURN( eval(B) ) fi; N := multigcd(map(numer,L)); if N <> 0 then N := N / icontent(N) fi; if N = 0 then RETURN(eval(B)) fi; D := multigcd(map(denom,L)); if D <> 0 then D := D / icontent(D) fi; if D = 0 then RETURN(eval(B)) fi; G := N/D; if G = 1 or G = -1 then RETURN(eval(B)) fi; G * map( unapply( x/G, x ), B ) end:
With the following command we enable to extract a submatrix from a matrix.
> with(linalg,submatrix);
Construction of matrix equations
In construction of these two matrix equations, in the following procedure we assign coefficients of to variables and to a specified order . This part of the computation is truly ``symbolic'' in that equations appear in terms of unknown matrices. We subsequently fix the rank of these matrices, fill in known entries and solve for unknown entries.
> makeEquations := proc(D) local i,H,U,E,G,Ut,UC; H := add(H.i*lambda^i,i=0..D); U := 1+add(U.i*lambda^i,i=1..D); E := H0 + add(E.i*lambda^i,i=1..D); G := expand(U &* H - E &* U); G := algsubs(lambda^(D+1)=0, G); Ut := 1 + add( transpose(U.i)*lambda^i, i=1..D ); UC := expand( U &* Ut - 1 ); UC := algsubs(lambda^(D+1)=0, UC); for i from 1 to D do G.i := coeff(G,lambda,i); UC.i := coeff(UC,lambda,i); od; RETURN(); end:
> makeEquations(3);
For instance here are linear coefficients as determined previously,
> G1; UC1;
and here are quadratic coefficients.
> G2; UC2;
Matrix is diagonal; this condition makes the algorithm work so simply. In the linear coefficient above, a quantity has diagonal entries that are zero; diagonal entries of are thus equal to those in . In the quadratic coefficient, the unknown quantities are and ; diagonal entries of the first expression are again zero. Diagonal entries of are hence given by diagonal entries of a matrix . The same observation holds for terms of greater order that have a general form . Diagonal entries of are simply determined and off-diagonal entries of give rise to simple linear equations. Here is the cubic coefficient.
> G3; UC3;
Construction of matrices H
In the following procedure we construct matrix for K>0 , bearing in mind that all matrices are formally of infinite rank. A parameter N specifies the rank of a truncated matrix. To obtain correct to N rows and N columns we compute ; we begin with N+K rows and columns for X .
> makeH := proc(K,N) local m,X,H; X := matrix(N+K,N+K,0); for m from 0 to N+K-2 do X[m+1,m+2] := A*sqrt(m+1); X[m+2,m+1] := X[m+1,m+2]; od; H := evalm( k.(K+2) * X^(K+2) ); submatrix(H,1..N,1..N); end:
Here is matrix truncated to rank 6.
> H1 := makeH(1,6): factormatrix(H1);
Formule for bands we obtain by interpolation:
for
and
for .
Likewise, here is matrix truncated to rank 6.
> H2 := makeH(2,6): factormatrix(H2);
Diagonal entries are given by
and the next band is given by
In general, matrices H are banded and symmetric with K+3 non-zero bands and with bands with zero entries between non-zero bands. Formulae for the bands can be determined by inspection and interpolation.
Solving for and .
The following programme in the form of a Maple procedure calculates matrices E and U to order K (entered as D ) truncated to rank N .
> AnHarmonic := proc(D,N) global U0,H0,E0,Order; local i,j,m,x,K,T,Z,U,E,H,eqns,vars,sols; makeEquations(D); U0 := array(1..N,1..N,identity); H0 := array(1..N,1..N,diagonal); for m from 0 to N-1 do H0[m+1,m+1] := (m+1/2)*omega od: E0 := H0; for K from 1 to D do print(`Solving for the coefficient of`, lambda^K); H.K := makeH(K,N); E.K := array(1..N,1..N,diagonal); U.K := array(1..N,1..N); U := evaln(U.K); E := evaln(E.K); H := evaln(H.K); printf(`matrix %s\n`,H); print(factormatrix(H.K)); printf(`matrix equation to be solved\n`); print(G.K); Z := evalm(G.K); if K=1 then printf(`expanded matrix equation to be solved\n`); print(Z) fi; eqns := normal(convert(Z,set)); vars := {seq(seq(U[i,j],j=1..N),i=1..N), seq(E[i,i],i=1..N)} minus {seq(U[i,i],i=1..N)}; sols := solve(eqns,vars); assign(sols); printf(`matrix %s\n`, E); print(factormatrix(E)); T := evalm( UC.K ); printf(`condition for the next U matrix`); print(UC.K = 0); eqns := normal(convert(T,set)); vars := {seq(U[i,i],i=1..N)}; sols := solve(eqns,vars); assign(sols); printf(`matrix %s\n`,U); print(factormatrix(U)); od; end:
We calculate and display then ; we display also the symbolic matrix equations. Like , all matrices E of odd order are 0. Before executing this procedure it is advisable to activate the following solver that greatly accelerates these calculations.
Fast solve
> unprotect(solve):
> solve := proc(eqns,unks) if type(eqns[1],`=`) then RETURN( solve( map(proc(x) lhs(x)-rhs(x) end,eqns) minus {0}, unks ) ) fi; map(proc(e,u) local x; if e=0 then RETURN() fi; x := indets(e) intersect(u); if nops(x)=1 then x := x[1]; x = normal( -coeff(e,x,0)/coeff(e,x,1) ) else ERROR(e,u) fi end, eqns, unks ) end;
The following test ensures that this procedure operates correctly.
> solve( {3*x=2,4*y=a}, {x,y} );
> AnHarmonic(2,6);
To verify that all computations performed with programme AnHarmonic are correct, that is, that we have correctly diagonalized H to , we test that equation holds to .
> eqn := (1 + add(lambda^i*U.i,i=1..2)) &* add(lambda^i*H.i,i=0..2) = add(lambda^i*E.i,i=0..2) &* (1 + add(lambda^i*U.i,i=1..2));
> map(simplify@series, evalm(lhs(eqn)-rhs(eqn)), lambda, 3 );
We proceed to examine the contribution to energy of second order; it is the matrix
> map(expand,E2);
Because we truncated matrices H , not all entries in matrices E and U are correct: actually only the first three diagonal entries in are correct. We seek to compare our results with those of Wu [Quantum Mechanics, World Scientific, 1985], who uses the following definition for H .
.
We convert our results for to conform to Wu's notation.
> W2 := subs(k3=k2/3!, k4=k3/4!, omega=omega[0]*h, submatrix(%,1..3,1..3));
These diagonal entries are consistent with equation (II-140) in Wu's book . The correction to energy of the harmonic oscillator in Wu's equation we obtain by interpolation from three values in the above matrix. Using Maple, we obtain
> interp([0,1,2],[W2[1,1],W2[2,2],W2[3,3]],n);
Two alternative ways to write this expression are as a polynomial in A or in . In Maple
> collect(%,A,factor);
> map(collect,series(%,n=-1/2),A);
Corrections to greater order are obtained by running programme AnHarmonic with matrices of sufficient rank. Five statements, as given below starting with "AnHarmonic(4,10):" for instance, suffice to generate corrections of a given order; the duration of calculation increases exponentially with order (according to this, or any alternative, approach).
We end this section with comments on computations. Firstly, with regard to construction of for K>0 , when we first computed these matrices, we used and etc. This procedure yields an algorithm to compute of which the running time is exponential in . An algorithm based on computing the matrix power is superior, and it is sufficiently quick that there is little point in using formulae for bands of .
The linear system of equations that is constructed and solved is a diagonal system, but, because formulae involve several square roots of small integers, Maple runs much more slowly than if coefficients were simple integers: the reason is that Maple is careful to determine algebraic relationships between square roots, e.g. . In our case, if one does not mind seeing and in the answer, one can speed calculations by a factor 10 or more if one solves equations ``blindly''; we devised for this purpose a solver presented in the paragraph above marked Fast solve.
Conclusion
We illustrate here how matrix mechanics can be effectively used to solve an important problem in quantum mechanics using perturbation theory. Because matrix is diagonal for a diatomic anharmonic oscillator, the requisite matrix algebra is simple; a consequence is that computations can be carried out exactly to high order using a system such as Maple for symbolic algebra.
We include in the following appendix statements that generate formulae for terms in energies of fourth, sixth and eighth orders.
Appendix Correction terms of Fourth, Sixth and Eighth Orders
> AnHarmonic(4,10):
> factormatrix(submatrix(E4,1..4,1..4));
> interp([0,1,2,3],[seq(E4[i,i],i=1..4)],n):
> collect(%,[A,k3,k4,k5,k6,omega],factor);
> evaln( E4[n,n] ) = map(collect, taylor(%,n=-1/2), A, factor);
Succeeding calculations require progressively increasing duration and memory. The first argument in each call to AnHarmonic indicates the perturbative order of a calculation, and the second argument specifies the rank of matrices involved in the calculation; this rank must increase exponentially with order to achieve accurate results.
> AnHarmonic(6,18):
> interp([0,1,2,3,4],[seq(E6[i,i],i=1..5)],n):
> evaln( E6[n,n] ) = map(collect, taylor(%,n=-1/2), A, factor);
E8 eighth order
> AnHarmonic(8,28):
> interp([0,1,2,3,4,5],[seq(E8[i,i],i=1..6)],n):
> evaln( E8[n,n] ) = map(collect, taylor(%,n=-1/2), A, factor);
