Quantum Mechanics (comments to: moloney@nextwork.rose-hulman.edu)        

In the Bohr theory of the hydrogen atom the electron is in a circular orbit of radius rb, 4rb, 9rb,.................., where the radius of the first orbit is 52.92 pm. In the quantum theory the electron is not confined to particular circular orbits; instead there is a certain probability that the electron will be found between distance r and r+dr from the nucleus. For the ground state (lowest possible electron energy) this probability is given by

P(r)dr=4rb^(-3)*r^2*exp(-2r/rb)*dr.

(a) For a ground state electron, make a plot of the probability function P versus the distance r of the electron from the nucleus.

(b) In order to find the most probable electron-nucleus distance r, differentiate P with respect to r and set equal to zero. How does the most probable value of r compare with the radius rb of the first Bohr orbit? Is your graph in agreement?

(c) Calculate the probability of finding the electron between r=0 and 2rb .

(d) Find the radius of a sphere which a ground state electron would spend 90% of its time inside

Solution:

Part A

restart; with(plots):
P:=4/(rb^3)*r^2*exp(-2*r/rb); # Probability function.
rb=5292*10^(-2); # Substitute rb into P;
P:=subs(rb=1323/25,P);

Plot graph of P vs r

plot(P,r=0..250,title=`Probability Function P vs Distance r`);

Part B

der:=diff(P,r); # Find the derivative of P with respect to r.
sol:=solve(der=0,r); # Set derivative equal to zero and solve for r.
evalf(sol);

Part C

Find area under P(r) curve between r=0 and 2*rb.

rb:=1323/25;
Prb:=int(P,r=0..2*rb);
evalf(Prb);

Part D

Find radius of sphere for which there is a 90% probability of finding the electron inside .

ans:=int(P,r=0..a);
fsolve(ans=.9,a); # It turns out to be around 141.

Download hatom.ms [Use 'Save File']

levelpar.ms

PHYSICS RESOURCE PACKETS PROJECT            File: LEVELPAR.MS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Dr. Moloney
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 3

Search engine for energy levels in a potential well.
Solves for even and odd solutions.
Other potential functions can be used in the program, but the program is not able to normalize every potential function.

h-bar^2 /2m is coefficient of first term. ( x' = ax, where x in m, x' in nm. a=10^9)

d/dx = dx'/dx d/dx'; x'=1nm

*** Some commands may return errors when you execute the worksheet. If so, run the same command again until it works. ***

restart;with(plots):
hbar:=6.6260755E-34;
Me:=9.1093987E-31;
ec:=1.60217733E-19;
ck:=evalf((hbar*1E9/(2*Pi))^2 /(2*Me * ec));

Enter in maximum width(nm) and depth(eV) of well

depth:=4; width:=2;

Enter in a function for the well shape

V := depth*((x/(width/2))^2-1)*(1-Heaviside(abs(x)-width/2));
g1:=plot(V,x=-width/2-1..width/2+1,color=blue,axes=framed,title=`V vs x`,labels=[nm,eV]): #see if V looks ok
g1;

Schrodinger's steady-state equation.

eq1:=-ck*diff(y(x),x,x)+V*y(x) = eV*y(x);

Solve Schrodinger's equation numerically as a function of eV.

s_odd:= k -> dsolve({subs(eV=k,eq1),D(y)(0)=1,y(0)=0},y(x),numeric);
s_even:= k -> dsolve({subs(eV=k,eq1),D(y)(0)=0,y(0)=1},y(x),numeric);

Look at what the equation does outside the barrier. The goal is to find the energy levels at which the wave function decays to 0 outside the barrier.

test_point:=width/2+4;
r_odd:= p -> (rhs(s_odd(p)(test_point)[2]));
r_even:= p -> (rhs(s_even(p)(test_point)[2]));

Draw vertical lines where the function outside the barrier changes signs, i.e. goes from negative to positive, or positive to negative. This tells us what the allowable values for eV are.

r_oddl:= p -> 3*sign(rhs(s_odd(p)(test_point)[2]));
r_evenl:= p -> 3*sign(rhs(s_even(p)(test_point)[2]));

Plot an energy band diagram of the allowable energy levels.

plot(r_oddl,-depth..0,y=-2..2,color=blue,labels=[eV,` `],ytickmarks=0,title=`odd energy levels`);
plot(r_evenl,-depth..0,y=-2..2,color=blue,labels=[eV,` `],ytickmarks=0,title=`even energy levels`);

Enter in approximations of where the energy levels fall on the graphs produced above. You need only enter the levels of the wave functions you wish to plot.

levels:=[even,-3.6],[odd,-2.8],[even,-2.05],[odd,-1.3],[even,-0.5];

This calulates the number of levels you have chosen to graph.

numlevels:=nops([levels]);

Definition of the derivative, used later to do Newton's method. Newton's method is used to a more exact value of eV for the approximations above.

dr:=x->(r(x+0.000001)-r(x))/0.000001;

Plot the wave functions specified in "levels"

for j from 1 to numlevels do

Pick an energy level out of "levels"

level:=levels[j][1];
mid:=levels[j][2];

level is used to determine which DE to solve.

if level='odd'
then r:=r_odd; s:=s_odd;
else r:=r_even; s:=s_even;
fi;

Use Newton's method to find an energy level that decays to 0 outside the barrier. "mid" is used as an initial guess. The loop exits after 10 iterations, or after a certain precision has been obtained.

mid0:=0: counter:=0:
while (abs(mid-mid0)>0.000000002) and (counter < 10) do
mid0:=mid:
mid:=mid0-r(mid0)/dr(mid0):
print(mid);
counter:=counter + 1:
od:

Plot the well, the wave function, and its energy level on the same graph

energy:=convert(mid,string);
g2:=odeplot(s(mid),[x,y(x)],-width/2-1..width/2+1,axes=framed,title=`Ev = `.energy,labels=[nm,eV]):
print(display([g2,g1]));

Print the energy levels to the screen.

print(s(mid.j));
mid.j:=mid;
s.j:=s(mid.j);
print(Ev=mid.j);
od:

Plot |psi|^2 and find the area under |psi|^2

for j from 1 to numlevels do
s_abs.j:=(x) -> rhs(s.j(x)[2])^2;
energy:=convert(mid.j,string);
plot(s_abs.j,-width/2-1/2..width/2+1/2,title=`original |psi|^2 Ev=`.energy);
area.j:=evalf(Int(s_abs.j,x=-width/2-1/4..width/2+1/4));
od;

Plot a the normalized wave functions

for j from 1 to numlevels do
s_norm.j:= (x) -> rhs(s.j(x)[2])/sqrt(area.j);
energy:=convert(mid.j,string):
g2:=plot(s_norm.j,evalf(-width/2-1..width/2+1),labels=[nm,eV],title=`normalized psi Ev=`.energy,axes=framed):
print(display([g2,g1]));
od:

Animate several wave functions in the well.

Value of hbar in this coordinate system.

hbar1:=hbar/ec;

Make the wave functions time dependent and find |Psi|^2 of the sum of the wave functions. Add or delete equations as desired.

y1:= (x,t) -> rhs(s1(x)[2])*exp(-I*mid1*t/hbar1)/sqrt(area1);
y2:= (x,t) -> rhs(s2(x)[2])*exp(-I*mid2*t/hbar1)/sqrt(area2);
y3:= (x,t) -> rhs(s3(x)[2])*exp(-I*mid3*t/hbar1)/sqrt(area3);
y4:= (x,t) -> rhs(s4(x)[2])*exp(-I*mid4*t/hbar1)/sqrt(area4);
y5:= (x,t) -> rhs(s5(x)[2])*exp(-I*mid5*t/hbar1)/sqrt(area5);
Y:=(x,t)->abs(y1(x,t)+y2(x,t)+y3(x,t)+y4(x,t)+y5(x,t))^2;

View value of function at any point and time. Y(x,t)

Y(1.5,0);
V1:= X -> subs(x=X,V);
animate({Y,V1},-2..2,0..5.1E-14,frames=50,axes=framed,labels=[nm,`|Psi|^2`],title=`Wave function sloshing in potential well`,color=red,numpoints=100);

Download levelpar.ms [Use 'Save File']

levelsqr.ms

PHYSICS RESOURCE PACKETS PROJECT            File: LEVELSQR.MS
Rose-Hulman Institute of Technology         Authors: Perry Peters & Dr. Moloney
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 3

Search engine for energy levels in a potential well.
Solves for even and odd solutions.
Other potential functions can be used in the program, but the program is not able to normalize every potential function.

h-bar^2 /2m is coefficient of first term. ( x' = ax, where x in m, x' in nm. a=10^9)

d/dx = dx'/dx d/dx'; x'=1nm

*** Some commands may return errors when you execute the worksheet. If so, run the same command again until it works. ***

restart;with(plots):
hbar:=6.6260755E-34;
Me:=9.1093987E-31;
ec:=1.60217733E-19;
ck:=evalf((hbar*1E9/(2*Pi))^2 /(2*Me * ec));

Enter in maximum width(nm) and depth(eV) of well

depth:=4; width:=2;

Enter in a function for the well shape

V := -depth*(1-Heaviside(abs(x)-width/2));
g1:=plot(V,x=-width/2-1..width/2+1,color=blue,axes=framed,title=`V vs x`,labels=[nm,eV]): #see if V looks ok
g1;

Schrodinger's steady-state equation.

eq1:=-ck*diff(y(x),x,x)+V*y(x) = eV*y(x);

Solve Schrodinger's equation numerically as a function of eV.

s_odd:= k -> dsolve({subs(eV=k,eq1),D(y)(0)=1,y(0)=0},y(x),numeric);
s_even:= k -> dsolve({subs(eV=k,eq1),D(y)(0)=0,y(0)=1},y(x),numeric);

Look at what the equation does outside the barrier. The goal is to find the energy levels at which the wave function decays to 0 outside the barrier.

test_point:=width/2+4;
r_odd:= p -> (rhs(s_odd(p)(test_point)[2]));
r_even:= p -> (rhs(s_even(p)(test_point)[2]));

Draw vertical lines where the function outside the barrier changes signs, i.e. goes from negative to positive, or positive to negative. This tells us what the allowable values for eV are.

r_oddl:= p -> 3*sign(rhs(s_odd(p)(test_point)[2]));
r_evenl:= p -> 3*sign(rhs(s_even(p)(test_point)[2]));

Plot an energy band diagram of the allowable energy levels.

plot(r_oddl,-depth..0,y=-2..2,color=blue,labels=[eV,` `],ytickmarks=0,title=`odd energy levels`);
plot(r_evenl,-depth..0,y=-2..2,color=blue,labels=[eV,` `],ytickmarks=0,title=`even energy levels`);

Enter in approximations of where the energy levels fall on the graphs produced above. You need only enter the levels of the wave functions you wish to plot.

levels:=[even,-3.9],[odd,-3.7],[even,-3.3],[odd,-2.8],[even,-2.1],[odd,-1.3],[even,-0.4];

This calulates the number of levels you have chosen to graph.

numlevels:=nops([levels]);

Definition of the derivative, used later to do Newton's method. Newton's method is used to a more exact value of eV for the approximations above.

dr:=x->(r(x+0.000001)-r(x))/0.000001;

Plot the wave functions specified in "levels"

for j from 1 to numlevels do

Pick an energy level out of "levels"

level:=levels[j][1];
mid:=levels[j][2];

level is used to determine which DE to solve.

if level='odd'
then r:=r_odd; s:=s_odd;
else r:=r_even; s:=s_even;
fi;

Use Newton's method to find an energy level that decays to 0 outside the barrier. "mid" is used as an initial guess. The loop exits after 10 iterations, or after a certain precision has been obtained.

mid0:=0: counter:=0:
while (abs(mid-mid0)>0.000000002) and (counter < 10) do
mid0:=mid:
mid:=mid0-r(mid0)/dr(mid0):
print(mid);
counter:=counter + 1:
od:

Plot the well, the wave function, and its energy level on the same graph

energy:=convert(mid,string);
g2:=odeplot(s(mid),[x,y(x)],-width/2-1..width/2+1,axes=framed,title=`Ev = `.energy,labels=[nm,eV]):
print(display([g2,g1]));

Print the energy levels to the screen.

print(s(mid.j));
mid.j:=mid;
s.j:=s(mid.j);
print(Ev=mid.j);
od:

Plot |psi|^2 and find the area under |psi|^2

for j from 1 to numlevels do
s_abs.j:=(x) -> rhs(s.j(x)[2])^2;
energy:=convert(mid.j,string);
plot(s_abs.j,-width/2-1/2..width/2+1/2,title=`original |psi|^2 Ev=`.energy);
area.j:=evalf(Int(s_abs.j,x=-width/2-1/2..width/2+1/2));
od;

Plot a the normalized wave functions

for j from 1 to numlevels do
s_norm.j:= (x) -> rhs(s.j(x)[2])/sqrt(area.j);
energy:=convert(mid.j,string):
g2:=plot(s_norm.j,evalf(-width/2-1..width/2+1),labels=[nm,eV],title=`normalized psi Ev=`.energy,axe
s=framed):
print(display([g2,g1]));
od:

Animate several wave functions in the well.

Value of hbar in this coordinate system.

hbar1:=hbar/ec;

Make the wave functions time dependent and find |Psi|^2 of the sum of the wave functions. Add or delete equations as desired.

y1:= (x,t) -> rhs(s1(x)[2])*exp(-I*mid1*t/hbar1)/sqrt(area1);
y2:= (x,t) -> rhs(s2(x)[2])*exp(-I*mid2*t/hbar1)/sqrt(area2);
y3:= (x,t) -> rhs(s3(x)[2])*exp(-I*mid3*t/hbar1)/sqrt(area3);
y4:= (x,t) -> rhs(s4(x)[2])*exp(-I*mid4*t/hbar1)/sqrt(area4);
y5:= (x,t) -> rhs(s5(x)[2])*exp(-I*mid5*t/hbar1)/sqrt(area5);
y6:= (x,t) -> rhs(s6(x)[2])*exp(-I*mid6*t/hbar1)/sqrt(area6);
y7:= (x,t) -> rhs(s7(x)[2])*exp(-I*mid7*t/hbar1)/sqrt(area7);
Y:=(x,t)->abs(y1(x,t)+y2(x,t)+y3(x,t)+y4(x,t)+y5(x,t)+y6(x,t)+y7(x,t))^2;

View value of function at any point and time. Y(x,t)

Y(1.5,0);
V1:= X -> subs(x=X,V);
g1:=animate(V1,-2..2,0..2E-14,frames=20,color=black,numpoints=200):
g2:=animate(Y,-2..2,0..2E-14,axes=framed,labels=[nm,eV],title=`Wave function sloshing in potential well`,color=red,frames=20,numpoints=100,color=blue):
display([g2,g1]);

Download levelsqr.ms [Use 'Save File']

slosh.ms (runs in either release 3 or release 4; some graphs in release 3 have glitches in them)

7/15/96 solve for a particle sloshing around in a double potential well

This is a program to find the lowest solutions of a double well problem so that an animation can be done of the electron tunneling back and forth between the sides of the barrier.

The 'radius' is the value where we run into the infinite well wall. The width is the barrier width and the height is that of the barrier.

restart;with(plots):
const:=(6.626E-34*1E9/(2*Pi))^2 /(2*9.1E-31 * 1.602E-19);
ck:=37645/(100000*Pi^2); # constant for electron mass, length in nm, energy in eV
V := a*Heaviside(x-b);
Decimals:=15;
width:=4; height:=1; radius:= 6;
plot(subs(a=height,b=width,V),x=0..radius); # half of the double well

eq1:=-ck*diff(y(x),x,x) +subs(a=height,b=width,V)*y(x) = eV*y(x);

Trial and error resulted in the values below. The origin of coordinates for the following calcs is in the corner of an infinite well.
The wave function starts off as sin (kx).

for j from 0.0139445*height to 0.013945*height by 0.0000001*height do
k:=sqrt(j/ck): j:
Y_even:=dsolve({subs(eV=j,eq1),D(y)(width)=k*cos(k*width),y(width)=sin(k*width)},y(x),numeric): rhs(Y_even(radius)[2]);
rhs(Y_even(radius)[2]); od;

The following lines of code were revised repeatedly to obtain a solution close to y=0 at the center of the barrier.

Energy:=13944793*height/1000000000; k:=sqrt(Energy/ck); # solve for y=0
Y_even:=dsolve({subs(eV=Energy,eq1),D(y)(width)=k*cos(k*width),y(width)=sin(k*width)},y(x),numeric): odeplot(Y_even,[x,y(x)],0..radius);Y_even(radius);

The following lines of code were revised repeatedly to obtain a solution close to dy/dx=0 at the center of the barrier.

Energy1:=139446326*height/10000000000; k1:=sqrt(Energy1/ck); # solve for dydx=0
Y_even:=dsolve({subs(eV=Energy1,eq1),D(y)(width)=k1*cos(k1*width),y(width)=sin(k1*width)},y(x),numeric):
odeplot(Y_even,[x,y(x)],0..radius);Y_even(radius);

x:=Energy-Energy1;
freq:=x*1.602E-19/6.636E-34;# difference between frequencies

get wave function for symmetric mode: cos (kx + c), cosh kappa x, kR+c = Pi/2

kappa1:=sqrt((height-Energy1)/ck);
eq1:= cos(Pi/2 - k1*width) = B*cosh(kappa1*(radius-width));
B1:=solve(eq1,B);
psi1:=-cos(Pi/2-k1*(y-radius))*Heaviside(y-(radius-width))+B1*cosh(kappa1*y)*Heaviside(radius-width-y);

plot(psi1,y=0..radius); # symmetric wave function
plot(subs(y=abs(z),psi1),z=-radius..radius,title=`symmetric wave fxn`);

get wave function for antisymmetric mode: cos (kx + c), sinh kappa x, kR+c = Pi/2

kappa:=sqrt((height-Energy)/ck);
eq1:= cos(Pi/2 - k*width) = D*sinh(kappa*(radius-width));
D1:=solve(eq1,D);
psi2:=-sin(k*(y-radius))*Heaviside(y-(radius-width))+D1*sinh(kappa*y)*Heaviside(radius-width-y);

plot(psi2,y=0..radius); # antisymmetric wave function
psiasm:=(Heaviside(z)-Heaviside(-z))*subs(y=abs(z),psi2);
plot(psiasm,z=-radius..radius,title=`anti-symm. wave fxn`);

Add symmetric and antisymm functions and animate

psitot:=psiasm+subs(y=abs(z),psi1);
plot(psitot,z=-radius..radius);
f1:=Energy1*1.602E-19/6.626E-34;
f:=Energy*1.602E-19/6.626E-34;
psitime:=psiasm*cos(f*t)+subs(y=abs(z),psi1)*cos(f1*t);

plot(subs(t=0,psitime),z=-radius..radius);
plot(subs(t=5E-8,psitime^2),z=-radius..radius); animate(psitime^2,z=-radius..radius,t=0..1.466666667E-7,frames=21,numpoints=150,color=blue);

well:=Heaviside(1-abs(z));
plot(well,z=-radius..radius);

animate({psitime^2,well},z=-radius..radius,t=0..1.466666667E-7,frames=21,numpoints=150,color=blue);

Download slosh.ms [Use 'Save File']

uncert.ms

HEISENBERG UNCERTAINTY RELATION

According to Heisenberg the minimum uncertainties in a particle's position and its momentum are related. The smallest value that is allowed the product of these two uncertainties is (roughly) Planck's constant h divided by 2pi (called "h bar"). This means that the smaller the uncertainty in its momentum the larger the uncertainty in its position and vice versa.

To represent a moving particle of uncertain momentum, we can use De Broglie's idea that there is a wave of wavelength h/p associated with a particle whose momentum is p. A particle of uncertain momentum can be represented by a superposition of waves of various wavelengths and wave numbers k, over a wavelength range which covers the uncertainty range of p, making up what is called a "wave packet". The purpose of this exercise is to illustrate how this superposition results in a wave function with a built-in position uncertainty, and to show that decreasing the range of momentum uncertainty broadens the position uncertainty. At the same time we would like to have a look at the "phase" and "group" velocities associated with the wave packet. [The latter is the speed at which the envelope moves, the former is the speed at which the wave inside the envelope moves; the group velocity is the most probable speed of the particle.] The group velocity is twice the phase velocity for nonrelativistic particles, so one should be able to see the envelope overtaking individual peaks as the packet moves along.

In order to build up a "Gaussian" wave packet, we start with a wave whose wave number ko represents the most probable particle momentum and add, over a range of wave numbers, smaller amplitude waves whose numbers represent momenta of lesser probability. The results in a summation of waves whose forms are: y:=exp(-c2*(k[i]-ko)^2)*sin(k[i]*x-c1*k[i]^2*t)

Note that the amplitude of a wave of wave number k[i] is smaller the farther it is fom ko, the wave number associated with the most probable particle momentum. The constant c1 is (hbar)/2m and constant c2 is 1/[2*delta_k)], where (delta_k) is the spacing between the k values).

The following assumes a wave packet for an electron with a speed of 5 x 10^5 m/s.

restart; with(plots):
v:=5*10^5;
m:= 9.1 * 10^(-31);
hbar:=1.055*10^(-34);
ko := m*v/hbar;
delta_k := ko/50;
c1:=hbar/(2*m):
c2:=1/(2*delta_k)^2;
N:=40;


kRange:=4*delta_k;


for i from 1 to N do
k[i]:=ko-kRange/2 + (i-1)*(kRange/N);
od:

s:=0;
for i from 1 to N do
s:=s+exp(-c2*(k[i]-ko)^2)*sin(k[i]*x-c1*k[i]^2*t);
od:


y:=subs(t=0,s):
plot(y,x=0..4*10^(-8),title=`Range = 4 ko / 50`);
animate(s,x=0..3*10^(-8),t=0..4*10^(-14),frames=20,numpoints=200);

kRange2:=2*delta_k;


for i from 1 to N do
k[i]:=ko-kRange2/2+(i-1)*(kRange2/N); od:

s2:=0; for i from 1 to N do
s2:=s2+exp(-c2*(k[i]-ko)^2)*sin(k[i]*x-c1*k[i]^2*t); od:
y:=subs(t=0,s2):
plot(y,x=0..4*10^(-8),title=`Range = 2 ko / 50`);
animate(s2,x=0..3*10^(-8),t=0..4*10^(-14),frames=20,numpoints=200);

Download uncert.ms [Use 'Save File']

Back to the Quantum Mechanics Resource Page