Harmonic Oscilatorjava Printing

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

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>. O 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/2p 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 w. We derived the potential by integrating the equation for the restoring force exerted by a spring, F=-K x, to get V(x) = ½ K x2, and then we expressed the spring constant K in terms of the natural frequency w 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: O -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(%,%%); 2

h_~ 1 K 2

d2 y x dx2 m~

CV x y x = E y x

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

V x =

2 1 m~ w~ x2 2

d2 y x 2 1 dx2 1 Eq A1 := K C m~ w~ x2 y x = E y x (1.2.1) 2 m~ 2 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 h_ units of . mw O Eq(A2):=x=sqrt(h_/m/omega)*u; Eq(A3):=expand(dchange(Eq (A2),Eq(A1),[u],params=[h_,m,omega])/h_/omega*2); h_~2

Eq A2 := x =

h_~ u m~ w~

2Ey u d2 y u Cu2 y u = (1.2.2) 2 h_~ w~ du We shift and scale the energy E and express it in terms of the dimensionless variable N, which simplifies the equation further. O Eq(A4):= E=h_*omega*(N+1/2); Eq(A5):=subs(%,Eq(A3)); 1 Eq A4 := E = h_~ w~ N C 2 2 d 1 Eq A5 := K y u Cu2 y u = 2 N C y u (1.2.3) 2 2 du 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. O Eq(A6):= collect(lhs(Eq(A5))-rhs(Eq(A5)),psi(u))=0; d2 2 Eq A6 := u K2 N K1 y u K y u =0 (1.2.4) du2 Eq A3 := K

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 u2y u in Eq(A6) is much bigger than ( 2 N +1 ) y u . No matter what N is, for sufficiently large arguments u the wavefunction y(u) will behave the same for all finite values of N. Let us look at the equation for N=0. O Eq(B1):=subs(N=0,Eq(A6)); d2 Eq B1 := u2 K1 y u K y u =0 (1.3.1) du2 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. O Eq(B2):=dsolve(Eq(B1),psi(u));

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

1 2 u 2

K

1 2 u 2

K

Eq B2 := y u = _C1 e C_C2 e erf i u (1.3.2) Let's plot the solutions one by one to see what they are like. The Gaussian looks like a good wave function. O exp(-u^2/2); plot(%,u=-5...5); 1 2 u 2

K

e

1

0.8

0.6

0.4

0.2

K4

K2

0

2

4 u

The other solution, the one that contains the error function erf, diverges for large arguments. O subs(_C1=1,_C2=1,rhs(Eq(B2))-exp(-u^2/2)); plot(%,u=-5...5, -10..10); 1 2 u 2

K

e

erf i u

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

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

10

5

K4

K2

0

2

4 u

K5

K10 Hence, we must rule out the latter solution as a wave function because if we tried to normalise it the normalisation integral would diverge. O Int(abs(psi(x))^2,x=-infinity..infinity)=limit(x,x= infinity); N KN

y x

2

dx = N

(1.3.3)

So, we gather that asymptotically, for large u, all wave functions must behave like the first solution, like exp(- ½ u2). How can we use this knowledge for finding the exact solution y u ? Well, we can factor out the asymptotic behaviour by expressing y u in terms of a new function h(u). O Eq(B3):=psi(u)=h(u)*exp(-u^2/2); 1 2 u 2

K

Eq B3 := y u = h u e Substituting this into Eq(A6), we get an equation for the function h(u). O Eq(B4):=combine(expand(subs(Eq(B3),exp(u^2/2)*Eq(A6))));

(1.3.4)

d2 d h u C2 h u u =0 (1.3.5) 2 du du 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 y u would still be alright. Eq B4 := K2 h u N K

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

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

1 2 u 2

K

Eq B5 := u/N lim u e

=0

(1.3.6)

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(-u2/ 2), as shown in Eq(B3), and we've got the exact wave function y u . 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 y u , 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. O h(u)=h(0)+Diff(h(0),u)*u+Diff(h(0),u$2)*u^2/2+`...`; d 1 d2 h u =h 0 C h 0 uC h 0 u2 C... (1.4.1) 2 du 2 du 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 uk we could call ak, 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 n = N.) 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. O 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); h u = a0 Ca1 u Ca2 u2 Ca3 u3 C... L

Eq C1 := h u =

>a u

k

k= 0

(1.4.2)

k

To find the unknown coefficients we simply insert the polynomial Eq(C1) into Eq(B4) for h(u). Differentiating powers is simple. O combine(simplify(subs(Eq(C1),Eq(B4)),power)); L

L

L

v v2 k K a u C2 ak uk u C K2 N ak uk = 0 (1.4.3) k 2 vu k = 0 vu k= 0 k= 0 Next we split the sum into two: one with uk and the other with uk K 2. 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). O sum(2*(k-N)*a[k]*u^k,k=0..L)-sum(k*(k-1)*a[k]*u^(k-2),k=2..

>

>

>

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

L)=0; L

>2

L

k

k= 0

k KN ak u K

>k

k= 2

k K1 ak uk K 2 = 0

(1.4.4)

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. O sum(2*(k-N)*a[k]*u^k,k=0..L)-sum(subs(k=j+2,k*(k-1)*a[k]*u^ (k-2)),j=0..L-2)=0; L K2

L

>2

k= 0

k

k KN ak u K

>

j C2

j= 0

j C1 aj C 2 u j = 0

(1.4.5)

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. O Eq(C2):=sum((2*(q-N)*a[q]-(q+1)*(q+2)*a[q+2])*u^q,q=0..L-2) + 2*(L-1-N)*a[L-1]*u^(L-1) + 2*(L-N)*a[L]*u^L=0; aL K 1 uL K 1 L K1 2 aL K 1 uL K 1 L K1 aL uL K 1 L K1 2 Eq C2 := K C K (1.4.6) u u2 u2 K

aL uL K 1 L K1

C

u K

> Ku

> Ku

q K2 2

q

q aq Cuq K 2 q aq C2 uq aq q K2 uq aq N

q K2 2

q

q aq Cuq K 2 q aq C2 uq aq q K2 uq aq N

q = L K1 K1 KN aL K 1 uL K 1 C2 L KN aL uL = 0

C2 L q =0

Now we have got a sum over powers of u, a polynomial in u, on the left-hand side. The righthand 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 uq must be zero: O Eq(C3):=2*(q-N)*a[q]-(q+1)*(q+2)*a[q+2]=0; Eq C3 := 2 q KN aq K q C1 q C2 aq C 2 = 0 (1.4.7) This equation does not let us evaluate the unknown coefficients ak 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. O Eq(C4):=a[q+2]=solve(Eq(C3),a[q+2]); 2 Kq CN aq Eq C4 := aq C 2 = K (1.4.8) q C1 q C2 Let's look at the remaining two terms. The coefficient of uL in Eq(C2) must be zero. We know that aL 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! O 2*(L-N)*a[L]=0; Eq(C5):=L-N=0; 2 L KN aL = 0 Eq C5 := L KN = 0 (1.4.9)

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

And the coefficient of uL K 1 in Eq(C2) must be zero. Since we have just concluded that L=N, we must have aL K 1 = 0. O 2*(L-1-N)*a[L-1]=0; Eq(C6):=a[L-1]=0; 2 L K1 KN aL K 1 = 0 Eq C6 := aL K 1 = 0

(1.4.10)

Now we have solved the problem. The solution for h(u) is a polynomial of degree L. The coefficient aL K 1 of the next lower power uL K 1is zero. Applying Eq(C4) we find that this means that aL K 3, aL K 5, ... 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). O N=0;simplify(subs(L=0,Eq(C1))); N =0 0

>a u

k

h u =

k= 0

(1.4.1.1)

k

Substituting in Eq(A4) and Eq(B3) and reversing the substituion of u for x according to Eq (A2) gives the solution: O 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),%); 1 Eq C7 := E = h_~ w~ 2 1 2 u 2

K

y u = a0 e

1 x2 m~ w~ 2 h_~

K

Eq C8 := y x = a0 e

(1.4.1.2)

The coefficient a0 is determined from the normalisation of y x . O 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)); N KN

y x

2

dx =

a20 h_~

p

h_~ m~ w~

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

a0 =

m~

w~

(1.4.1.3) h_~ p The plot of the wave function is familiar from Worksheet 7. O plot(subs(m=1,h_=1,omega=1,a[0]=Pi^(-1/4),rhs(Eq(C8))),x=

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

-4..4);

0.7 0.6 0.5 0.4 0.3 0.2 0.1

K4

K3

K2

K1

0 x

1

2

3

4

First excited state The next possible value of L is 1. Inserting this into Eq(C5) and Eq(C1) gives: O N=1;simplify(subs(L=1,Eq(C1))); N =1 1

h u =

>a u

k

(1.4.2.1)

k

k= 0

Eq(C6) tells us that a0 must vanish, which simplifies h(u). O a[0]=0; h(u)=a[1]*u; a0 = 0 h u = a1 u

(1.4.2.2)

As before we substitute the above into Eq(A4) and Eq(B3) and use Eq(A2) to find the solution: O 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),%); 3 Eq C9 := E = h_~ w~ 2 1 2 u 2

K

y u = a1 u e

1 x2 m~ w~ 2 h_~

K

Eq C10 := y x =

a1 x e

(1.4.2.3) h_~ m~ w~

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

Again, the coefficient a1 is determined from the normalisation of y x . O 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)); N KN

y x

2

a21 h_~

1 dx = 2

p

h_~ m~ w~

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

a1 =

w~

m~

2

(1.4.2.4)

p

h_~

And the plot of the wave function: O plot(subs(m=1,h_=1,omega=1,a[1]=sqrt(2)*Pi^(-1/4),rhs(Eq (C10))),x=-4..4);

0.6

0.4 0.2

K4

K3

K2

K1

0

1

2 x

3

4

K0.2 K0.4 K0.6 Second excited state Next is L=2. Eq(C5) and Eq(C1) give: O N=2;simplify(subs(L=2,Eq(C1))); N =2 2

h u =

>a u

k

k= 0

Eq(C6) demands that a1vanishes. O a[1]=0; h(u)=a[0]+a[2]*u^2;

k

(1.4.3.1)

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

a1 = 0 h u = a0 Ca2 u2

(1.4.3.2)

We use Eq(C4) to determine how a2 depends on a0. O subs(q=0,N=2,Eq(C4)); h(u)=a[0]*(1-2*u^2); a2 = K2 a0 h u = a0 1 K2 u2

(1.4.3.3)

Then we get from Eq(A4), Eq(B3), and Eq(A2): O 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),%); 5 Eq C11 := E = h_~ w~ 2 2

y u = a0 1 K2 u

1 2 u 2

K

e

1 x2 m~ w~

K 2 x2 m~ w~ Eq C12 := y x = a0 1 K e 2 h_~ h_~ As before, the coefficient a0 is determined from the normalisation of y x .

(1.4.3.4)

O 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)); N KN

y x

2

dx =

2 a20 h_~

p

h_~ m~ w~

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

a0 =

1 2

2

m~ h_~

w~

(1.4.3.5)

p

And again the plot of the wave function. O plot(subs(m=1,h_=1,omega=1,a[0]=Pi^(-1/4)/sqrt(2),rhs(Eq (C12))),x=-4..4);

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

0.4

0.2

K4

K3

K2

K1

0

1

2 x

3

4

K0.2 K0.4

K0.6 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.) O for j from 0 to 9 do h[j](u)=H(j,u) od; h0 u = 1 h1 u = 2 u h2 u = K2 C4 u2 h3 u = 8 u3 K12 u h4 u = 12 C16 u4 K48 u2 h5 u = 32 u5 K160 u3 C120 u h6 u = K120 C64 u6 K480 u4 C720 u2 h7 u = 128 u7 K1344 u5 C3360 u3 K1680 u h8 u = 1680 C256 u8 K3584 u6 C13440 u4 K13440 u2 h9 u = 512 u9 K9216 u7 C48384 u5 K80640 u3 C30240 u

(1.4.4.1)

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 O Eq(A4); 1 E = h_~ w~ N C (1.4.4.2) 2

Please purchase PDFcamp Printer on http://www.verypdf.com/ to remove this watermark.

The wave function of this state is the product of the Nth order Hermite polynomial and a Gaussian. Correctly normalised it reads: O subs(h(u)=(m*omega/Pi/h_)^(1/4)/sqrt(2^N*N!)*h[N](u),Eq (B3)); `with `*u=solve(Eq(A2),u);

y u =

m~ w~ p h_~

with u =

1/4

1 2 u 2

K

hN u e 2N N! x h_~ m~ w~

(1.4.4.3)