Analytic Approach To The Harmonik Oscillator

Copyright 2000 Arne Bergström, B&E Scientific Ltd, and Claudia Eberlein, University of Sussex.

Worksheet 8 - Analytic approach to the harmonic oscillator

Instructions for use of this worksheet

This Maple worksheet has been written for Maple 6.01 and has been tested both under Windows98 and under Linux. To run the worksheet you just need to place the cursor anywhere in the first command group below (printed in red and starting with restart) and press <ENTER> to progress from one command group to the next. If there are animations you will be given instructions on how to play them in the text directly above them.

If you are not familiar with Maple then read the text (printed in black, like this line), ignore the Maple commands (printed in red), and just look at the output of Maple's calculations (displayed plots or formulae printed in blue). If you know Maple then you are encouraged to read also the Maple commands and try out what happens when you change them. However, you are still advised first to go through the whole worksheet without changing anything, so that you understand the aim of the calculation.

Start by putting the cursor anywhere in the command group listed below, and then press <ENTER>.

> restart: with(PDEtools,dchange):

with(orthopoly):

interface(showassumed=0,imaginaryunit=i):

assume(h_>0,m>0,omega>0):

Proceed by pressing <ENTER>.

Notation

Planck's constant h/2 which is commonly called "h bar" cannot be printed in Maple. Therefore it is replaced by the symbol h_ throughout the worksheet.

A. The Schrödinger equation for the harmonic oscillator

In Worksheet 7 we wrote down the stationary Schrödinger equation for an harmonic oscillator of mass m and natural frequency . We derived the potential by integrating the equation for the restoring force exerted by a spring, F=-K x, to get V(x) = ½ K , and then we expressed the spring constant K in terms of the natural frequency and the mass m, in order to be able to apply the equation to a general oscillator and not just a mass suspended on a spring. So, we got for the stationary Schrödinger equation:

> -h_^2/2/m*diff(diff(psi(x),x),x)+V(x)*psi(x)=E*psi(x);

V(x)=m/2*omega^2*x^2; Eq(A1):=subs(%,%%);

As before, we want to simplify the equation because we don't want to carry around all sorts of parameters. So, we scale the variable x to a new variable u which is the displacement measured in units of

.

> Eq(A2):=x=sqrt(h_/m/omega)*u; Eq(A3):=expand(dchange(Eq(A2),Eq(A1),[u],params=[h_,m,omega])/h_/omeg a*2);

We shift and scale the energy E and express it in terms of the dimensionless variable N, which simplifies the equation further.

> Eq(A4):= E=h_*omega*(N+1/2); Eq(A5):=subs(%,Eq(A3));

Finally, we can bring the equation into the standard form for differential equations by shuffling the term on the right-hand side over to the left.

> Eq(A6):= collect(lhs(Eq(A5))-rhs(Eq(A5)),psi(u))=0;

B. Asymptotic solution for large arguments

The best approach to solving a second-order differential equation like Eq(A6) analytically, is to deal with the problem step by step and start by figuring out how the solutions behave asymptotically, that is to say, for large arguments.

For large arguments u, the term in Eq(A6) is much bigger than ( 2 N +1 ) . No matter what N is, for sufficiently large arguments u the wavefunction (u) will behave the same for all finite values of N. Let us look at the equation for N=0.

> Eq(B1):=subs(N=0,Eq(A6));

A second-order differential equation has always two linearly independent solutions. The general solution is a linear superposition of these two, with the integration constants _C1 and _C2 as coefficients.

> Eq(B2):=dsolve(Eq(B1),psi(u));

Let's plot the solutions one by one to see what they are like. The Gaussian looks like a good wave function.

> exp(-u^2/2); plot(%,u=-5...5);

The other solution, the one that contains the error function erf, diverges for large arguments.

> subs(_C1=1,_C2=1,rhs(Eq(B2))-exp(-u^2/2)); plot(%,u=-5...5,10..10);

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

Hence, we must rule out the latter solution as a wave function because if we tried to normalise it the normalisation integral would diverge.

> Int(abs(psi(x))^2,x=-infinity..infinity)=limit(x,x=infinity);

So, we gather that asymptotically, for large u, all wave functions must behave like the first solution, like exp(- ½ ). How can we use this knowledge for finding the exact solution ? Well, we can factor out the asymptotic behaviour by expressing in terms of a new function h(u).

> Eq(B3):=psi(u)=h(u)*exp(-u^2/2);

Substituting this into Eq(A6), we get an equation for the function h(u).

> Eq(B4):=combine(expand(subs(Eq(B3),exp(u^2/2)*Eq(A6))));

Provided that h(u) is reasonably well-behaved at infinity, i e not exponentially rising, the substitution in Eq(B3) ensures that the wave function has the correct asymptotic behaviour and is normalisable. This is because an exponential falls off faster than any power, which means that h(u) could be any combination of powers and the asymptotic behaviour of the wave function would still be alright.

> u^n * exp(-u^2/2): Eq(B5):=Limit(%,u=infinity)=limit(%,u=infinity);

C. Solution for small arguments and exact solution

The problem is now to solve Eq(B4) and find h(u); once we have that, we just need to multiply it with the asymptotic behaviour exp(- / 2), as shown in Eq(B3), and we've got the exact wave function . The trouble is that Eq(B4) is still a second-order differential equation, and as such it is just as difficult to solve as the original equation for , Eq(A6). So, what if anything have we gained by factoring out the asymptotic behaviour in Eq(B3) and now facing Eq(B4) instead of Eq(A6)? At first sight - not much. However, if one cannot solve a differential equation then the best approach is usually to try and solve it asymptotically for large arguments, then solve it for small arguments, and finally try and "knit together" the two approximate solutions in the middle between small and large arguments. In quite a few case this strategy even leads to the exact solution of the differential equation, and this will turn out to happen here.

We know that the asymptotic behaviour for large arguments is taken care of by the factorisation in Eq(B3). So, let's try and work out a solution for h(u) for small arguments. For small u we must be able to Taylor expand h(u) about u=0.

> h(u)=h(0)+Diff(h(0),u)*u+Diff(h(0),u$2)*u^2/2+`...`;

Since we don't know the derivatives of h(u) at u=0 (but these are the unknowns we want to find), we can as well call the expansion coefficients something else: the coefficient in front of we could call , say. We have to work out what these coefficients are, but one thing we know for sure - we cannot have infinitely many of them because that would destroy our asymptotic behaviour for large u. (Eq(B5) is not valid for .) In other words, the power series has to break off at some finite number, say, at L. So, h(u) is a ploynomial of degree L, and what L is doesn't matter but it must be finite.

> h(u)=a[0]+a[1]*u+a[2]*u^2+a[3]*u^3+`...`; Eq(C1):=h(u)=sum(a[k]*u^k,k=0..L);

To find the unknown coefficients we simply insert the polynomial Eq(C1) into Eq(B4) for h(u). Differentiating powers is simple.

> combine(simplify(subs(Eq(C1),Eq(B4)),power));

Next we split the sum into two: one with and the other with . In the latter the terms for k=0 and k=1 vanish (because the coefficients work out to zero - constants and linear terms disappear when one differentiates them twice).

> sum(2*(k-N)*a[k]*u^k,k=0..L)-sum(k*(k-1)*a[k]*u^(k-2),k=2..L)=0;

Instead of summing k from 2 to L in the second sum we can as well substitute k=j+2 and sum j from 0 to L-2.

> sum(2*(k-N)*a[k]*u^k,k=0..L)-sum(subs(k=j+2,k*(k-1)*a[k]*u^(k2)),j=0..L-2)=0;

If we replace both k and j by a new index, say, q then we can amalgamate the two sums into one except for the (L-1)st and the Lth term in the first sum, which we have to write separately.

> Eq(C2):=sum((2*(q-N)*a[q]-(q+1)*(q+2)*a[q+2])*u^q,q=0..L-2) + 2*(L1-N)*a[L-1]*u^(L-1) + 2*(L-N)*a[L]*u^L=0;

Now we have got a sum over powers of u, a polynomial in u, on the left-hand side. The right-hand side tells us that it must equal zero, and that must be so for all u - because, after all, we want to find a function h(u) that satisfies the differential equation Eq(B4) for all u and not just for a few special values of u. The only way to make the left-hand side zero for all u is to have zero as coefficients in front of all powers of u. So, the coefficient in front of must be zero:

> Eq(C3):=2*(q-N)*a[q]-(q+1)*(q+2)*a[q+2]=0;

This equation does not let us evaluate the unknown coefficients but it tells us how one depends on the previous but one. So, the 7th depends on the 5th, the 5th on the 3rd, and the third on the 1st. And the 8th depends on the 6th, the 6th on the 4th, the 4th on the 2nd, the 2nd on the 0th.

> Eq(C4):=a[q+2]=solve(Eq(C3),a[q+2]);

Let's look at the remaining two terms. The coefficient of in Eq(C2) must be zero. We know that cannot be zero (because then h(u) wouldn't be a polynomial of degree L but of degree L-1 or lower); hence L-N must be zero. Wow, that means that N must be a positive integer. Exactly what we found numerically in Worksheet 7!

> 2*(L-N)*a[L]=0; Eq(C5):=L-N=0;

And the coefficient of have .

in Eq(C2) must be zero. Since we have just concluded that L=N, we must

> 2*(L-1-N)*a[L-1]=0; Eq(C6):=a[L-1]=0;

Now we have solved the problem. The solution for h(u) is a polynomial of degree L. The coefficient of the next lower power is zero. Applying Eq(C4) we find that this means that , ... must also be zero. If L is even then the coefficients of all odd powers of u are zero, and if L is odd then the coefficients of all even powers of u are zero. So, the solution is either an even or an odd polynomial in u. That is something we actually also already know because we've figured out in Worksheet 7 that the solution must have either even or odd parity. Even polynomials give even-parity solutions, and odd polynomials odd-parity solutions.

Let's look at a few examples.

Ground state

The lowest possible value of L and therefore N is zero. We insert this into Eq(C1).

> N=0;simplify(subs(L=0,Eq(C1)));

Substituting in Eq(A4) and Eq(B3) and reversing the substituion of u for x according to Eq(A2) gives the solution:

> Eq(C7):=subs(N=0,Eq(A4)); subs(h(u)=a[0],Eq(B3));

Eq(C8):=subs(psi(u)=psi(x),u=solve(Eq(A2),u),%);

The coefficient is determined from the normalisation of

.

> Int(abs(psi(x))^2,x=-infinity..infinity)=int(rhs(Eq(C8))^2,x=infinity..infinity); a[0]=sqrt(solve(simplify(rhs(%))=1,a[0]^2));

Warning, solving for expressions other than names or functions is not recommended.

The plot of the wave function is familiar from Worksheet 7.

> plot(subs(m=1,h_=1,omega=1,a[0]=Pi^(-1/4),rhs(Eq(C8))),x=-4..4);

First excited state

The next possible value of L is 1. Inserting this into Eq(C5) and Eq(C1) gives:

> N=1;simplify(subs(L=1,Eq(C1)));

Eq(C6) tells us that must vanish, which simplifies h(u).

> a[0]=0; h(u)=a[1]*u;

As before we substitute the above into Eq(A4) and Eq(B3) and use Eq(A2) to find the solution:

> Eq(C9):=subs(N=1,Eq(A4)); subs(h(u)=a[1]*u,Eq(B3));

Eq(C10):=subs(psi(u)=psi(x),u=solve(Eq(A2),u),%);

Again, the coefficient is determined from the normalisation of

.

> Int(abs(psi(x))^2,x=-infinity..infinity)=int(rhs(Eq(C10))^2,x=infinity..infinity); a[1]=sqrt(solve(simplify(rhs(%))=1,a[1]^2));

Warning, solving for expressions other than names or functions is not recommended.

And the plot of the wave function:

> plot(subs(m=1,h_=1,omega=1,a[1]=sqrt(2)*Pi^(-1/4),rhs(Eq(C10))),x=4..4);

Second excited state

Next is L=2. Eq(C5) and Eq(C1) give:

> N=2;simplify(subs(L=2,Eq(C1)));

Eq(C6) demands that vanishes.

> a[1]=0; h(u)=a[0]+a[2]*u^2;

We use Eq(C4) to determine how depends on .

> subs(q=0,N=2,Eq(C4)); h(u)=a[0]*(1-2*u^2);

Then we get from Eq(A4), Eq(B3), and Eq(A2):

> Eq(C11):=subs(N=2,Eq(A4)); subs(h(u)=a[0]*(1-2*u^2),Eq(B3));

Eq(C12):=subs(psi(u)=psi(x),u=solve(Eq(A2),u),%);

As before, the coefficient is determined from the normalisation of

.

> Int(abs(psi(x))^2,x=-infinity..infinity)=int(rhs(Eq(C12))^2,x=infinity..infinity); a[0]=sqrt(solve(simplify(rhs(%))=1,a[0]^2));

Warning, solving for expressions other than names or functions is not recommended.

And again the plot of the wave function.

> plot(subs(m=1,h_=1,omega=1,a[0]=Pi^(-1/4)/sqrt(2),rhs(Eq(C12))),x=4..4);

Hermite polynomials

We could go on and play this game for all integer L, and thus work out the wave functions for the Lth excited state. But some mathematicians have done this already, and the polynomials h(u) that we have been calculating are called the Hermite polynomials. Here is a listing of the first ten of them. (note that they have got alternately even and odd parity, as we mentioned earlier.)

> for j from 0 to 9 do h[j](u)=H(j,u) od;

So, we can write down the general solution of the stationary Schrödinger equation for the harmonic oscillator. The Nth excited state has the energy

> Eq(A4);

The wave function of this state is the product of the Nth order Hermite polynomial and a Gaussian. Correctly normalised it reads:

> subs(h(u)=(m*omega/Pi/h_)^(1/4)/sqrt(2^N*N!)*h[N](u),Eq(B3)); `with `*u=solve(Eq(A2),u);

>