UNIVERSIDAD DE PUERTO RICO EN HUMACAO
www.uprh.edu

 

 

 Departamento de Física

 

LECTURES ON QUANTUM MECHANICS

by Reinaldo Baretti Machín

 reibaretti2004@yahoo.com

 

  Free counter and web stats

Chapter 5.  Potential Wells

 

 

www1.uprh.edu/rbaretti/LQMIntro.htm

www1.uprh.edu/rbaretti/LQMch1.htm

www1.uprh.edu/rbaretti/LQMch2.htm

www1.uprh.edu/rbaretti/LQMIch3.htm

www1.uprh.edu/rbaretti/LQMch4.htm 

 

 

www1.uprh.edu/rbaretti

www1.uprh.edu/rbaretti/methodsoftheoreticalphysics.htm

www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart2.htm

www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart3.htm

 

 

 

 

Para preguntas y sugerencias escriba a :

reibaretti2004@yahoo.com

 

 

 

 

References:

 

 1. Principles of quantum mechanics; nonrelativistic wave mechanics with illustrative applications.  

by W. V. Houston

 

2. Introduction to Quantum Mechanics with Applications to Chemistry by Linus Pauling and E. Bright Wilson Jr.

 

3. Quantum Mechanics: Non-Relativistic Theory, Volume 3, Third Edition (Quantum Mechanics) by L. D. Landau and L. M. Lifshitz

 

4.Wikipedia   http://en.wikipedia.org/wiki/Old_quantum_theory

 

5. The Rise of the New Physics ** 2 Volumes Complete** by A. D'Abro

 

6. Foundations of Physics by Robert B. Lindsay , Henry Margenau

 

7. LECTURES ON QUANTUM MECHANICS by Gordon Baym

 

 8. Wave Mechanics (Vol. 5 of Pauli Lectures on Physics) (Pauli Lectures on Physics Volume 5)

 

 

FREE DOWNLOAD FORTRAN-FORCE 2.0.8

 

 

In a broad sense every particle in a bound state is caught in a certain well. By bound state it is meant

 that both  one or both Ψ (+∞) →0, Ψ (-∞) → 0.Or the function is finite at one boundary and zero at the other. The particle in the box, though simple, exhibits the features of all eigenvalue problems.

Using Schrodinger's equation solutions ( plural) are sought that satisfy the boundary conditions . In this case they become zero at the boundaries. From this condition a set of discrete energy eigenvalues E1 , E 2... En  results accompanied by a set of corresponding eigenfunctions  Ψ1 ,  Ψ2  ....Ψn  .   

 

 One starts with the general statement     H Ψ = E Ψ    and ends with a set  { E i } ,  { Ψn } such that ,

 

                                                            H Ψn = En Ψn     .                   

 

 

The infinite one  dimensional  well

   

 

Fig 1. Square well with infinte potential at the walls.

 

Let   V= ∞ when x < 0 and x > L . V(x)=0   ,   0<x<L  and the boundary condition be

 Ψ (0) = Ψ (L) =0. The walls are impenetrable because they have an infinite height.

 

 

 Schrodinger equation is

 

    d2 Ψ /dx2 = -(2mE/h'2 ) Ψ = -kΨ                           .  (1)

 

The two independent solutions are

Ψ (x) = A sin(kx) + B cos(kx)                                 .      (2)

Applying the boundary condition at the origin ,gives , Ψ (0) =0 = B  . The cosine part drops out.

Applying the boundary condition at x= L , gives Ψ (L) ==0 = A sin(kL)  which implies kL= nπ  , which is  the quantization condition.

 

The constant A will be determined from the normalization condition. It is not necessary to find the energy since

                                En  =  (h'2 /(2m) )  (π/L)2 n2   , n=1,2,3.....                        (3)

 

The normalization condition requires that   ∫ Ψ2 dx = 1 = A2 ∫ sin (kx)2 dx    ,     0 < x < L

                                                                         = A2 ∫ ( 1-cos (2kx) )/2  dx = A2 L/2   .

 

Thus the normalization constant is simply   A = (2/L)1/2 .The normalized eigenfunctions are

                      

                                  Ψ n (x) =  (2/L)1/2 sin ( n π x/L)                              .           (4)

 

The set  Ψ n (x) is orthonormalized   ∫ Ψ n (x) Ψ m (x)  dx = δ n,m  .  ( δ n,m =1 ,if n=m ;δ n,m =0 for n≠m ).

 

 

 

The delta function potential

Some properties of the Delta function in one dimension  are

 

                             ∫ f(x) δ(x- x0 ) dx = f (x 0 )                    ,        -∞ < x < +∞       (5)

 

      

 

                                              ∫ δ(x- x0 ) dx =1                    ,                              (6)

 

 

and                                          ∫ δ( a x ) dx =1/a                .                               (7)

 

The dimensions of   δ( x ) are 1/length  .

A dirac sequence is established when certain functions approach the delta function  with the help of an adjustable parameter. Two example are given. First consider the gaussian function

 

                              P(x) = (1/σ) ( 1/(2π)1/2 ) exp[ -(x-μ)2 /(2σ2 )]                 ,        (8)

 

μ is the mean and σ2 the variance. In other words

                                           ∫ x P(x) dx =  μ       ,      -∞ < x < +∞                      (9)

 

and                                     ∫ x2 P(x) dx =  σ2  .                                             . (10)                                           

                               

 

 

 

 

Fig 2. The gaussian function P(x) will resemble Dirac delta function as σ→ 0 .

 

Figure 2 shows that P(x) becomes sharply peaked as  the variance ( σ ) goes to zero.

 

 

 

A second  sequence that converges to the Dirac delta function is 

                                                                                        (11)          

It resembles the dirac delta function as n becomes large .

 

Its integral   ∫ δn (x) dx = 1    ,  - ∞ < x < + ∞ .

 

 

 

Let the potential be V(x) =- η δ(x)  where η is the strength of the delta function. This peculiar potential produces a single energy eigenvalue  given by E = - {h'2 /(2m)} η2  . To show this start with

 

                      d2 Ψ /dx2  - (2m/h'2 ) V(x) Ψ = - (2m/h'2 ) E Ψ                         .   (12)

Integrate the whole expression around the origin  from  -ε to  +ε . We get 

 

 

 ∫ (d2 Ψ /dx2 )  dx =  ( d Ψ( +ε ) /dx   - d Ψ( -ε ) /dx ) = 2 ( d Ψ( +ε ) /dx )           ,  (13)

see next figure 3.

Fig 3. Shape of  the eigenfunction of the delta function.

 

The potential term integration gives

 

- (2m/h'2 )∫ V(x) Ψ dx = - (2m/h'2 )∫  -η δ(x)  Ψ dx =  (2m/h'2 ) η  Ψ(0)    .      (14)

 

The last term is zero , i.e.

 

- (2m/h'2 ) E ∫ Ψ  dx =  - (2m/h'2 ) E  Ψ( 0) (2ε ) → 0                           .     (15)

 

From (13) and (140 it follows that

 

                     d Ψ( +ε ) /dx )  = (-1/2)(2m/h'2 ) η  Ψ(0)  =  -(m/h'2 ) η  Ψ(0)  . (16)

The function near the origin ( x > 0) is thus of the form   

 

                 Ψ   ~     exp{-(m η /h'2 ) x }  .                                                  (17)

To get the energy go back to (12) at x near the origin (x>0) and use (17)

 

      E Ψ  = - {h'2/(2m) } d2 Ψ /dx2 = - (1/2) (m/h'2 ) η2 Ψ                ,        (18) 

 

                   E = -m η2 /(2 h'2 )  ,

or                E = - η2 /2     in units where m=1,h'=1.  

 

This is the only allowed energy eigenvalue by an attractive delta function of strength η.

Figure 4 shows a delta like potential. Its eigenvalue is E = - η2 /2 = -50.

 

Sage code for the plot

 

var('x')
mu=0 ; sigma=1.e-3 ;eta=10;
gauss= (1/sigma)*(1/(2*pi)^(1/2))*exp(-(x-mu)^2/(2*sigma^2));
v=-eta*gauss;
#v =x;
y=plot(v,x,-3*sigma,3*sigma)
show(y)

eta,-eta**2/2.= 0.1000E+02 -0.5000E+02

 

Fig 4. The potential   V(x)= - η (1/σ) ( 1/(2π)1/2 ) exp[ -(x-μ)2 /(2σ2 )]  , η=10,  μ=0 ,σ = 1.E-3.

 

We integrated Schrodinger  equation  numerically ,starting at x=0, with initial conditions Ψ(0)=1 , dΨ(0) /dx =-1. The integration is carried out up to  xfinal = 250 σ  for different values of the energy. The eigenvalue becomes apparent when Ψ ( xfinal) changes sign.

e,psif= -0.9667E+02 0.4139E+01
e,psif= -0.9333E+02 0.3697E+01
e,psif= -0.9000E+02 0.3279E+01
e,psif= -0.8667E+02 0.2885E+01
e,psif= -0.8333E+02 0.2513E+01
e,psif= -0.8000E+02 0.2162E+01
e,psif= -0.7667E+02 0.1832E+01
e,psif= -0.7333E+02 0.1522E+01
e,psif= -0.7000E+02 0.1230E+01
e,psif= -0.6667E+02 0.9578E+00
e,psif= -0.6333E+02 0.7027E+00
e,psif= -0.6000E+02 0.4646E+00
e,psif= -0.5667E+02 0.2428E+00
e,psif= -0.5333E+02 0.3659E-01
e,psif= -0.5000E+02 -0.1546E+00

e,psif= -0.4667E+02 -0.3315E+00
e,psif= -0.4333E+02 -0.4946E+00
e,psif= -0.4000E+02 -0.6446E+00
e,psif= -0.3667E+02 -0.7820E+00
e,psif= -0.3333E+02 -0.9073E+00
e,psif= -0.3000E+02 -0.1021E+01
e,psif= -0.2667E+02 -0.1124E+01
e,psif= -0.2333E+02 -0.1216E+01
e,psif= -0.2000E+02 -0.1299E+01
e,psif= -0.1667E+02 -0.1371E+01
e,psif= -0.1333E+02 -0.1435E+01
e,psif= -0.1000E+02 -0.1490E+01
e,psif= -0.6667E+01 -0.1537E+01
e,psif= -0.3333E+01 -0.1576E+01
e,psif= 0.0000E+00 -0.1608E+01

 

 


 

                

c FORTRAN code for the numericalintegration of delta function potential
implicit real*8(a-h,o-z)
real*8 k
dimension psif(200),energy(200),steps(200)
data ne,epsi,nx/30,1.d-2,2000/
gauss(x) =(1.d0/(sigma*dsqrt(2.d0*pi)))*
$ dexp(-(x-amu)**2/(2.d0*sigma**2))
V(x)=-eta*gauss(x)
c v(u)=-eta*delta(u,epsi)
pi=2.d0*asin(1.d0)
amu=0.d0
sigma=1.d-3
eta=10.d0
c ei= -eta/(2.d0*epsi)
ei=-eta**2
ef=0.d0
de=(ef-ei)/dfloat(ne)
do 20 ie=1,ne
e=ei + de*dfloat(ie)
k=dsqrt(-2.d0*e)
xi=0.d0
c xf= -(1.d0/k)*log(1.d-2)
xf=250.d0*sigma
c dx=epsi/10.d0
dx=(xf-xi)/dfloat(nx)
c epsi=dx
kp=int(float(nx)/40)
kount=kp
psi0=1.d0
deriv=-1.d0
psi1=psi0+dx*deriv
do 10 i=2,nx
x=xi+dx*dfloat(i)
psi2=2.d0*psi1-psi0 -dx**2*(2.d0*(E-V(x-dx))*psi1)
c psi2=2.d0*psi1-psi0 -dx**2*(2.d0*(E-V(x-dx,eta,epsi))*psi1)
c if(i.eq.kount)then
c print 130,x,psi2
c kount=kount+kp
c endif
psi0=psi1
psi1=psi2
10 continue
energy(ie)=e
psif(ie)=psi2
20 continue
print 120 ,eta,-eta**2/2.d0
print*,' '
do 30 i=1,ne
print 100,energy(i),psif(i)
30 continue
100 format('e,psif=',3(3x,d11.4))
120 format('eta,-eta**2/2.=',2(3x,d11.4))
130 format('x,ps=',2(3x,d10.3))
stop
end

 

 

Fig 5. Eigenfunction (not normalized) for the delta function potential.

 

 

 

 

 

The Linear Oscillator

The problem of the harmonic oscillator is one of paramount importance in quantum mechanics.

 

For the one dimensional oscillator the  hamiltonian is

 

                              H = -  {h'2/(2m)} d2/dx2  + (1/2) k x2                                    ,       (19)

accordingly     

                                 -  {h'2/(2m)} d2 Ψ/dx2  + (1/2) k x2  Ψ = E Ψ             .         (20)   

 

 

scale units             h'2/(mL2) ~ k L2      , Lscale = ( h'2 /(mk) )1/4                          (21)

 

                             Escale ~ k  (Lscale )2  = h' (k/m)1/2 = h' ω  ,                           (22)

where ω is the classical frequency of the oscillator. 

It was shown in chapter two that  the energy dependence on the quantum number n has the form  En ~  n2α/(2+α)  if V(x) ~ xα .

For the harmonic oscillator  α =2   and  En  ~ n , or more explicitly    En ~ n h' ω + constant term. 

We cannot physically have simply En ~ n h' ω , because   n=0 would be satisfied only by Ψ(x)=0. Thus a constant term must be present in the formula for En .                              

( See Ref  http://www1.uprh.edu/rbaretti/LQMch2.htm for a discussion of the dependency of E upon the quantum number n.)

Replacing  h'=1,m=1,k=1 in (20)   transforms the equation in

 

                                                  - (1/2) d2 Ψ/dx2  + (1/2)  x2  Ψ = E Ψ        .  (23) 

 

 

The set of eigenfunctions Ψis an orthogonal set. Let    H Ψn = En  Ψn   and  H Ψm = EmΨm   . Multiply the first by

Ψm and the second by Ψn .

 

 

Substract   ∫ ( Ψm H Ψnn H Ψm   ) dx =  0   =  ( En - Em)∫  Ψm Ψn  dx . Therefore

 

                                                        ∫  Ψm Ψn  dx =0 .                

 

 We will seek the first eigenvalues and eigenfunctions using what is called the creation operator.

 

 

When (1/2)x2 >>E,    the asymptotic  form of (23)  is      - (1/2) d2 Ψ/dx2  + (1/2)  x2  Ψ ≈ 0    and the approximate solution for x large (x>>1) is

                                                   

                                                          Ψ asymptotic = exp(-x2 / 2)     .               (24)

 

- (1/2) d2 Ψasymptotic /dx2  + (1/2)  x2  Ψasymptotic = (1/2) exp(-x2 /2)→ 0  as x → +∞ or  -∞ .

Consider now the substitution of the normalized function  Ψ = (1/π)1/4 exp(-x2 /2) ≡ N exp(-x2 /2)  in eq (23).

 

The left side is  

                   - (1/2) d2 Ψ/dx2  + (1/2)  x2  Ψ = (1/2) (1/π)1/4 exp(-x2 /2) =(1/2) Ψ  .  (25)

Equated to E Ψ    gives    E0 =1/2 the energy of the ground state. Hence the complete expression for the eigenvalues is 

                                        E n = (n +1/2) h'ω , or  

                                        En   = (n+1/2)                                                   (26)

 

Figure 6 shows the normalized ground state wave function Ψ0 = (1/π)1/4 exp(-x2 /2).

 

Fig 6 . Ground state wave function of the harmonic oscillator.

 

Part a: Creation and annihilation operators

 

The first excited state has n=1 , E= 3/2 and can be obtained with the rising or creation operator defined by

   

                       a+ = (1/21/2 ) ( -i p +x) = (1/21/2 ) ( - d/dx + x )                                        . (27)

 

Operating with with the rising operator a+ upon the ground state ,

 

a+ Ψ0 = (1/21/2 ) ( - d/dx + x ) (1/π)1/4 exp(-x2 /2) = (  21/2 / π1/4 )  x exp(-x2/2) = Ψplus  ,  (28)

 

the subscript  label (plus)  attached to psi  is provisional.

 

SAGE CODE

psi0= (1/pi^(1/4))*exp(-x^2/2);
psiplus=(1/2^(1/2) )*(-diff(psi0,x)+x*psi0 );

psiplus 

sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4)

 

We integrate next  psi squared square,  (Ψplus )2 , finding   ∫  (Ψplus )2 dx  =  1.       ,     -∞ < x < + ∞ .

 

 

 

Sage code

psiplus =sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4)

integral(psiplus^2,x,-oo,+oo)

1

 

Evaluate now   the average energy    E = ∫ Ψplus H Ψ plusdx  /  ∫ (Ψplus )2 dx

 

                                                        = ∫ Ψplus {- (1/2) d2 /dx2  + (1/2)  x2 }  Ψ plus dx 

 

 

                                                        = 3/2                                                           .  (30)

 

When (30) is equated to   En    we have (n+1/2)   = 3/2   ; n=1.

Function Ψplus is renamed    Ψ1 and corresponds to the first excited state of the quantum oscillator,

                                  Ψ1  =   ( 21/2 / π1/4 )  x exp(-x2/2)    ,  E1 =3/2   .                (31)

In conclusion, the effect of the rising operator upon Ψ0 creates a function proportional ,or equal , to Ψ1 .

Fig 7. Normalized Ψ1 .

 

Evidently    ∫ Ψ1Ψ0 dx = 0 since the integrand is an odd function in the interval  -∞ < x < + ∞ .

 

 

 

Applying      a+ on Ψ1 , produces    a+ Ψ1 = Ψplus = (1/π1/4 )*(2x2 -1)*exp(-x2/2)  . To normalize Ψplus , we find

 

                           ∫  (Ψplus )2 dx  = ∫ {(1/π1/4 )*(2x2 -1)*exp(-x2/2) } 2 dx= 2 .

 

It turns out that the  function  obtained  from  Ψplus / 21/2 = (1/21/2) (1/π1/4 )*(2x2 -1)*exp(-x2/2) , is  Ψ2 as we verify next. Its energy is

 

                                 E = ∫ Ψ2  {- (1/2) d2 /dx2  + (1/2)  x2 }  Ψ2 dx

 

                                     =  5/2 = ( 2 + 1/2)    ,                                                        (32)

 

i.e.    n=2 .

 

So                      Ψ2   =   (1/21/2) (1/π1/4 )*(2x2 -1)*exp(-x2/2)  ,  E2  = (2 +1/2)     . (33)            

 

Sage code

f=sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4);
(1/2^(1/2))*(-diff(f,x) +x*f)

 

 

psiplus= (1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
integral(psiplus^2,x,-oo,+oo)

2

 

Sage code energy eigenvalue (n=2)

psi2=(1/2^(1/2))*(1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
integral(psi2*( (-1/2)*diff(psi2,x,2) +(1/2)*x^2*psi2),x,-oo,+oo)

 

 

psi2=(1/2^(1/2))*(1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
psi3=(1/3^(1/2))*(1/2^(1/2))*(-diff(psi2,x) +x*psi2);
#integral(psi3^2,x,-oo,+oo)
# 1
integral(psi3*( (-1/2)*diff(psi3,x,2) +(1/2)*x^2*psi3),x,-oo,+oo)
# 7/2

Sage code

f=sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4);
(1/2^(1/2))*(-diff(f,x) +x*f)

 

 

psiplus= (1/pi^(1/4))*(2*x^2 -1)*exp(-x^2/2);
integral(psiplus^2,x,-oo,+oo)

2

 

 

 

Sage Code

 

psiplus=sqrt(2)*x*e^(-1/2*x^2)/pi^(1/4);
hpsiplus= (-1/2)*diff(psiplus,x,2) +(1/2)*x^2*psiplus;
integral(psiplus*hpsiplus,x,-oo,+oo)

3/2

 

 

 

 

Do it one more time. Apply  a+ upon   Ψ2  . Guided by the previous example, we suppose that

 

                       a+  Ψ2 = (3)1/2  Ψ3 = Ψplus                                                 .              (34)

 

  Hence                                     Ψ3 = (1/31/2 )(1/π1/4 ) *exp( -x2/2)*( 2 x3 - 3x )  .   (35)

We verify that ∫ Ψ2 dx =1 ; it is normalized.

psi3= (1/3^(1/2) )*(1/pi^(1/4))*exp( -x^2/2)*( 2*x^3 - 3*x );
integral(psi3^2,x,-oo,oo)

 

The energy eigenvalue is  

                                   E3 = ∫ Ψ3  {- (1/2) d2 /dx2  + (1/2)  x2 }  Ψ3 dx

                                     = 7/2 = ( 3+1/2) .                                                       (36)

d2=-(1/2)*diff(psi3,x,2)
integral(psi3*(d2+(1/2)*x^2*psi3),x,-oo,oo)

 

Thus applying the  creation operator on a normalized eigenfunction    gives,

                                                       a+ Ψn =  (n)1/2 Ψn+1                     ,                (37)

or in  a more abstract form  

                                                    a+ │ n>  = (n+1)1/2 │ n+1 >      ,               (38) 

The vectors │ n> are an abstraction whose projection in x space is written as

                                      

          Ψn (x)  = < x│n> . Alternatively the function in p -space (momentum ) is

 

          φn (p)  = < p│n> .

 

From the Fourier transform defintion     φn (p) = ∫ dx exp(-px) Ψn (x) ,  this can be rewritten as

 

                             φn (p) = ∫ dx <p│ x> <x │n  >  .

The term    <p│ x> =exp(-ipx)    and  <x│ p> = exp( +ipx)    ( remember we use  h'=1)  .                                  

 

The annihilation operator , a , is the conjugate of a+ ,

                                                               

                                a ≡ (1/21/2 ) ( i p +x) = (1/21/2 ) (  d/dx + x )              . (39) 

 

One readily verifies that

                                 a Ψ0 =(1/21/2 ) (  d/dx + x )(1/π)1/4 exp(-x2 /2)=0       , (40)

 

 psi0=(1/pi^(1/4))*exp(-x^2/2) ;
(1/2^(1/2))*(diff(psi0,x)+x*psi0)

0                          

                               

 and                                a Ψ1 = (1/21/2 ) (d/dx + x ) ( 21/2 / π1/4 )  x exp(-x2/2)

                                         = (1/π)1/4 exp(-x2 /2)   =  Ψ0                            .     (41)

 

psi1=(2^(1/2)/pi^(1/4))*x*exp(-x^2/2);
(1/2^(1/2))*(diff(psi1,x)+x*psi1)

 

 
e^(-1/2*x^2)/pi^(1/4)

Also                              a Ψ2 = (2 /π1/4 ) x  exp(-x2/2) = 21/2 Ψ1             (42)  

and in general                a  Ψn = n1/2  Ψn-1                                     ,       (43)

or                                    a│ n>  = (n)1/2 │ n-1 >                         .   (44)

2*x*e^(-1/2*x^2)/pi^(1/4)

Consider   the average    < n │a+  a│n> = < n │a+  n1/2 │n-1 > = n1/2 < n │a+  │n-1 >

                                                                         =  n1/2 < n │ (n-1+1)1/2  │n-1+1 >

                                                           = n  < n │n > = n                (45)

   Thus    

                                          En  = n + 1/2 =   <n  │ a+  a    +1/2│ n>      , (46)

                                               =      <n  │ H│ n>    .                                 (47)

                                                =     Hn,n  , all non diagonal values are zero.

 

The physical dimensions of these particular operators , a+ and a, are  energy1/2 ~ (h' ω)1/2 = h'1/2(k/m)1/4 .

In many fields of physics  the use of creation and annihilation operators is preferred  over the use of wave functions. The technique is also called second quantization.

Example: Find the matrix elements of  a , a + , x ,  p , x2 and p2.

 

                                  <n │a│ m > =   <n│(m)1/2 │ m-1 >  

                                 a n,m              = m1/2 δ n,m-1           . (48)

The only non zero elements are those when the column index is one more than the row index.  

Thus               

                               m =    0      1      2       3     4

                           an,m  =      0      1      0      0     0 .......

                                         0      0     21/2  0     0

                                         0      0     0      3/2   0        

      

    and                          <n │a+│ m > =   <n│(m+1)1/2 │ m+1 >

                                            a+ n,m      =  (m+1)1/2  δ n,m+1               .       (49)

 

                                m =     0      1      2       3     4

                           a+ n,m  =      0      0      0      0     0 .......

                                           1       0       0    0     0

                                           0      21/2     0    0   0

                                           0      0       31/2  0   0

 

Evidently        <n │a+ a │ m> =         0    0    0   0     0 ......

                                                         0    1    0   0    0

                                                         0    0    2    0   0

                                                         0    0     0   3   0

The matrix elements of x are

 

<n │x │ m> = (1/2)1/2 ( a+  +  a ) n,m 

                   

                                                    0      1       0      0          0 .......

                        =  (1/2)1/2             1       0       21/2    0       0

                                                    0      21/2     0      31/2    0

                                                    0      0        31/2    0       41/2 

The matrix elements of p are

 

<n │p │ m> = (i/2)1/2 ( a+  -  a ) n,m 

                                      0      -1        0        0       0

                   =  (i/2)1/2    1       0       -21/2    0       0

                                       0      21/2     0     - 31/2    0

                                       0      0        31/2    0      - 41/2 

 

 

<n │x2 │ m> = <n │p2 │ m> = (1/2)   1   0   0    0.......

                                                          0   3    0    0....

                                                          0   0     5    0......

                                                          0   0     0    7....

For more details on creation and annihilation operators see ,

http://en.wikipedia.org/wiki/Creation_and_annihilation_operators

 

 

 Part b: Hermite polynomials

We saw above that   exp(-x2/2) is both the asymptotic solution and the ground state wave function ( except for a

normalization constant). Now  write

 

                                        Ψ(x) = f(x) exp(-x2/2)                      ,       (50)

and seek the differential equation for the function f(x) . The solution of f (x) is sought by series method with the proviso that it is a finite polynomial.  Susbstitution of (50) in  (23) gives ,

 

 

                                d2 f/dx2 -2x df/dx -( 1-2E)f = 0                  .         (51)                        

Suppose  first  f (x) = xs  ∑n=0 an xn . Near the origin , xs is the leading power. Let f ~ xs in (51)   ,

 

                              s(s-1) xs-2 - 2 sxs -( 1-2E)xs ~ 0      .

The leading power near the origin is  xs-2 thus  either s=0  or s= +1 . As will be shown consequence this gives rise to even or odd power series.Let s=0 , then  f (x) =  ∑n=0 an xn  and

 

 d2 f/dx2 = ∑n an  n(n-1) xn-2                                      (52-a)

 

-2x df/dx  =∑n -2 n an  xn                                           (52-b)

 

-( 1-2E)f=  -(1-2E)∑ an xn                            .              (52-c)

 

Alternatively setting s=1 and f(x) =  ∑n=0 an xn+1 leads essentially to the same algebraic relations as above.

 

The factors of the xn power must ad up to zero , i.e. letting n→ n+2 in term (52-a) we get

 

 

an+2 (n+2)(n+1)  + ( - 2n -(1-2E) ) an = 0    , producing the recursion formula

 

                       

                     an+2  = {2n +(1-2E)} an /((n+2)(n+1))                  .    (53)

If a0 ≠ 0   the recursion relation (53 ) , yields the even power series a ,  a4 ,   a6 ..... .

If a1≠ 0 the recursion relation , yields an odd power series a ,  a5 ,   a7 .....

 

To terminate the series at the nth power one must set 2n +(1-2E) = 0 , producing the quantization condition,

 

                                           En = n +1/2.

Example : Let n=0 , a0=1   ,  E =1/2   ; then   a2 = 0 . The function f(x) = a0 and the un normalized Ψ is

 

                                              Ψ 0 = exp( -x2/2)                  .        (54)

 

Let   n=1  ,  a1 =1 ,  E=3/2 ;then    a3 = 0  , f(x) = a1x = x  . The un normalized Ψ is

 

                                             Ψ 1= x exp( -x2/2)       .                (55)

 

 

Let   n=2  ,  a0 =1 ,  E=5/2 ;then    a2 = -2  , f(x) =  a0 + a2 x2 = 1-2x2 . The un normalized Ψ is

 

                                             Ψ 2=  ( 1-2x2 ) exp( -x2/2)       .                (56)

The first four Hermite polynomials (see e.g Ref 2) , in one of two possible conventions, are  written as

 

                                   H0(x) =1   ;   H1(x) = 2x     ; H2(x) = 4 x2 -2    ,  H3 (x)= 8 x3-12 x  (57)

which can be obtained from the formula

 

                            Hn(x) = (-1)n exp( x2 ) { dn exp(-x2) /dxn }             .  (58)

Our function f(x) of these examples is proportional to a Hermite polynomial. The wavefunctions of the harmonic oscillator are of the form                         

                                          Ψn (x)  = Nn Hn(x) exp(-x2/2)             ,    (59)

 where one can find by induction that the normalization constant is  Nn = (1/π1/4) {1/(2n n! )}1/2  .

 

Exercise : Justify the expression for Nn by normalizing  successively   Ψ 012 , Ψ 3  etc.

First calculate the integrals ∫ (Ψ n )2 dx to find out what pattern is established.

 

The normalization constant of  Ψ0 ,  is obtained from,

∫ (Ψ 0 )2 dx = ∫ (H0 exp(-x^2/2) )2 dx = π1/2  = 1/ (N0 )2  , thus   N0 = 1/π1/4

SAGE code

h0=1; psi0=h0*exp(-x^2/2);
s0=integral(psi0^2,x,-oo,oo);
n0=1/s0^(1/2);
s0 , n0

 

 

The normalization constant of  Ψ1 ,  is obtained from ,

 

  ∫ (Ψ 1 )2 dx = ∫ (H1 exp(-x^2/2) )2 dx π1/2 = S1 = 2 π1/2  =  1/(N1)2  ; thus N1 =(1/ π1/4 ) 1/(2)1/2 .

 

  The normalization constant of  Ψ2 ,  is obtained from ,                     

  ∫ (Ψ 2 )2 dx = ∫ (H2 exp(-x^2/2) )2 dx π1/2 = S2 = 8 π1/2   =  1/(N2)2 ; thus N2=(1/ π1/4 ) 1/( 22 *2)1/2 .

 

One last normalization constant .The normalization constant of  Ψ3 ,  is obtained from

  ∫ (Ψ 3 )2 dx = ∫ (H3 exp(-x^2/2) )2 dx π1/2 = S3 = 48 π1/2   =  1/(N3)2 ; thus N3 =(1/ π1/4 ) 1/( 23 *3!)1/2 .

Only the last two constants finally give away the formula

 

                                    Nn  = (1/ π1/4 ) 1/( 2n * n!)1/2                                  .     (58)

The normalized eigenfunctions in adimensional units are

 

                                  Ψ 2 (x)  =  {(1/ π1/4 ) 1/( 2n * n!)1/2  } Hn (x) exp(-x2/2)    . (59)  

 

 

      The following SAGE code can be used to plot some solutions (see http://www.sagenb.org/home/pub/590/  ). 

       x = var('x')
def u(n):
fac = ((2*factorial(n))^(-1/2))*(2^(1/4))
return fac * hermite(n, sqrt(2*pi)*x) * exp(-pi*x^2)

u_0(x) = u(0)
u_2(x) = u(2)
u_4(x) = u(4)

graph1 = plot(u_0(x), -3, 3, rgbcolor=(1, 0, 0))
graph2 = plot(u_2(x), -3, 3, rgbcolor=(0, 1, 0))
graph3 = plot(u_4(x), -3, 3, rgbcolor=(0, 0, 1))

show(graph1 + graph2 + graph3)

                           

 

 

 

 

Part c: Numerical solution ( last but not least)

 

Equation (23) becomes a finite difference equation when the second derivative is approximated by

 

                               d2 Ψ /dx2 ≈   (Ψn   -2Ψ n-1  + Ψ n-2  )/ (∆x)2                          .  (60)

The unit of length is Lscale  = ( h'2 /(mk) )1/4 ≡ 1 . In the integration one must choose ∆x <<  Lscale =1  .

 

The solution to

  - (1/2) d2 Ψ/dx2  + (1/2)  x2  Ψ = E Ψ   is,

 

    Ψn =  2Ψ n-1  - Ψ n-2  - 2 (∆x)2 {  E- (1/2) (xn-1)2  }Ψ n-1                            ,            (61)

 

where n is not an eigenvalue number but a supscript,   Ψn = Ψ (xn)   ,  Ψn-1 = Ψ (xn-1) = Ψ (xn - ∆x)  ,

Ψn-2 = Ψ (xn-2) = Ψ (xn - 2∆x) . 

The initial values are arbitrary ,

Ψ0 = Ψ(0)  = Ψ1 =Ψ( ∆x) =1   , ie.  (dΨ/dx)x=0 =0    for even functions   like the ground state  , the second excited state ,the fourth excited state, etc.

 

For odd functions Ψ0 = Ψ(0)  = 0 , (dΨ/dx)x=0 = 1 , Ψ1 = Ψ0 + ∆x .

 

A FORTRAN code is provided so that Eq (61) is integrated for a variable E until xasymtotic  >>  xclassical turning point . The classical turning point is defined from   E = V(x) = (1/2) x2  ;  xclassical turning point = (2 E) 1/2 . The "asymptotic "  Ψ (xasymtotic ) flips from positive to negative or vice versa when it passes through an energy eigenvalue.

 

 

 

Fig 8 . Even states  energy eigenvalues are found at E= 0.5 ,  2.5,  4.5,  6.5 .

 

 

 

Fig 8 . Odd states energy eigenvalues are found at E= 1.5 ,  3.5,  5.5  .

 

 

FORTRAN code

c Quantum harmonic oscillator
dimension psinf(100),energy(100)
data niter,nstep/80,2000/
v(x)=.5*x**2
e=.05
efinal=7.
de=(efinal-e)/float(niter)
do 110 iter=1,niter
xi=0.
xlim=2.5*sqrt(2.*e)
dx=xlim/float(nstep)
kp=int(float(nstep)/60.)
kount=kp
c IC for even functions
psi0=1.
psi1=psi0
c IC for odd functions
c psi0=0.
c psi1=psi0+dx
c print 200 ,0.,psi0,psi0/gauss(0.)
do 100 i=2,nstep
x=xi+dx*float(i)
psi2=2.*psi1-psi0 -dx**2*2.*(E-v(x-dx))*psi1
if(i.eq.kount)then
c print 200 ,x, psi2,psi2/gauss(x)
c print 200 ,-x,-psi2,-psi2/gauss(x)
kount=kount+kp
endif
psi0=psi1
psi1=psi2
100 continue
energy(iter)=e
psinf(iter)=psi2
20 continue
e=e + de
110 continue
do 30 i=1,niter
print 120 ,energy(i),psinf(i)/abs(psinf(i))
30 continue
120 format(3x,'energy,psinf/abs(psinf)=',2(3x,e10.3))
200 format(3x,'x,psi,h=',3(3x,e10.3))
stop
end

 

Finite Square Well Potential

A  finite square well (next figure) can be  defined by 

 

  V(x) =0     ,    -L < x < L

   V(x) = V0   ,           x>L    , x < -L.

The solutions can be classified according to their parity.

Starting with the ground state a set of solutions are symmetric  Ψ(x)  = Ψ(-x). Another set starting with the first excited state are anti symmetric  ,  Ψ(x)  = - Ψ(-x) .

The symmetric states ( even functions) will "resemble"  cos (k x) where k =(2E)1/2 and  k L ≈ (2n + 1)π/2 . It would be an equality if we had an infinite well. The anti symmetric states ( odd functions) will resemble sin(kx)  and kL ≈ nπ .In the regions /x/ >L the functions are decaying exponentials. Depending on the height V0  only a finite number of energy eigenvalues is permitted.

 

 

 

 

 

Fig 9 .Sample of class notes about the finite square potential.

 

For even states

 

                  Ψ 1 (x) = A exp (k1x)   ,    x< -L                   ,  k1 = ( 2 (V0 -E) )1/2   ,     (62)

                  Ψ 2 (x) = B cos (k2x)    ,       -L < x < L       ,   k2 = ( 2 E )1/2      ,          (63)

                  Ψ 2 (x) = A exp (-k1x)   ,    x > L    .

 

  From the continuity conditions on Ψ  and dΨ /dx   , at x= -L, we get , 

                                          A exp(-k1L) = B cos (k2L)                                           (64)

                                               k1 A exp(-k1L) = k2 B sin (k2L)                                       (65)

Dividing eq (65) by eq (64) gives the transcendental equation ,

 

                                  k1  =  k2 tan ( k2L)             .                                                (66)

The solutions can be obtained by graphical or numerical means .

 

 

For odd states the continuity equations ,(at x = -L) , lead to

 

                                                      A exp(-k1L) = - B sin (k2L) 

                                                   k1 A exp(-k1L) = k2 B cos (k2L) 

and the transcendental equation is now

 

                                                   k1  =  - k2  cot (k2L)                             .                     (67)

                                                          

Figures (10) and (11) show the allowed energies (three) for a finite square well with V0 = 8.

 

Fig 10 . Even state energies are found near E = 0.80  , 6.4  . 

 

A refinement of the energies is given by Newton's  formula for the roots. Define

         

                                            f (E) =     k1(E)  -  k2 (E)  tan ( k2(E)*L)  .

 

The roots are found from  the iteration formula

 

                                         En =  En-1  - f (En) / ( df/dE)n   .                         (68)

 

Ten iterations for the ground state E =0.784194

  i, e= 1 0.5
i, e= 2 0.939939141
i, e= 3 0.838132858
i, e= 4 0.790814519
i, e= 5 0.784336388
i, e= 6 0.784195423
i, e= 7 0.784194291
i, e= 8 0.78419435
i, e= 9 0.784194231
i, e= 10 0.784194291

 

Fortran code for Newton's roots formula

c roots of k1 = k2 tan ( k2*L)
real L , lscale ,k1,k2
data L ,V0 /1.,8. /
k1(e) = (2.* (V0 -E) )**(1./2.)
k2(e)=(2.*e)**(1./2.)
c f(e) for even states
f(e)=k1(e)-k2(e)*tan(k2(e)*L)
c f(e) for odd states
c f(e)= k1(e)+ k2(e)/tan(k2(e)*L)
df(e)= (f(e+epsi)-f(e))/epsi
Lscale=1/v0**.5
epsi=lscale/100.
c initial E must be selected near an eigenvalue
e=0.5
do 10 i=1,10
e1=e - f(e)/df(e)
print*,'i, e= ', i, e
e=e1
10 continue
stop
end


 

 

 

 

 

 

 

Fig 11. Odd states , only one energy is allowed ; E ≈ 3.0   .

 

The first three energy eigenvalues for the  infinite square well of size 2L  (letting L=1) are ,

 

En =  { π2 /(8 L2 ) } n2  = 1.23 n2  = 1.23   ,  4.92  , 11.07  . Since V0 is finite , the allowed penetration of  Ψ into the classically forbidden region represents a larger lambda and consequently a lower energy.

 

c energies for the finite square well
real k1,k2, L
data L,ne /1.,80/
k1(e) = (2.*(V0 -E))**(1./2.)
k2(e)=(2.*E)**(1./2.)
v0=8./L**2
print*,'Eisberg,.098*v0,.8*v0=',.098*v0,.8*v0
print*,' '
ei=0.01
ef=v0
de=(ef-ei)/float(ne)
e=ei
do 10 i=0,ne
even=k2(e)*tan(k2(e)*L)
odd= k2(e)/tan(k2(e)*L)
if(abs(even).lt.8.) print 100,e, k1(e) ,even
e=e+de
10 continue
100 format('e,k1,even=',3(3x,e10.3))
stop
end

A refinement of the energies

 

 

 

 

 

 

Finite square well  by numerical integration

 

Fig  12. Approximate energies of the first two and only even states  are  E ≈ 0.8   ,  6.76  . Compare with the estimate in figure 10.

 

 

   Fig 13 . Energy of the odd state ≈ 3.075  . The four digits are read from the data results.

 

Applying Newton's method  ,eq. (68) we get  E= 3.061.

 i, e= 1 2.5
i, e= 2 3.18810344
i, e= 3 3.06919932
i, e= 4 3.06180215
i, e= 5 3.06176519
i, e= 6 3.06176543
i, e= 7 3.06176496
i, e= 8 3.06176519
i, e= 9 3.06176543
i, e= 10 3.06176496
 

 

c numerical integration of finite square well
dimension psinf(100),energy(100)
real L
data L,v0 , niter,nstep/1.,8.,80,2000/
e=.05
efinal=v0
de=(efinal-e)/float(niter)
do 110 iter=1,niter
xi=0.
xlim=1.6*L
dx=xlim/float(nstep)
kp=int(float(nstep)/60.)
kount=kp
c IC for even functions
psi0=1.
psi1=psi0
c IC for odd functions
c psi0=0.
c psi1=psi0+dx
c print 200 ,0.,psi0,psi0/gauss(0.)
do 100 i=2,nstep
x=xi+dx*float(i)
psi2=2.*psi1-psi0 -dx**2*2.*(E-v(x-dx,L,v0))*psi1
c if(i.eq.kount)then
c kount=kount+kp
c endif
psi0=psi1
psi1=psi2
100 continue
energy(iter)=e
psinf(iter)=psi2
20 continue
e=e + de
110 continue
do 30 i=1,niter
print 120 ,energy(i),psinf(i)/abs(psinf(i))
30 continue
120 format(3x,'energy,psinf/abs(psinf)=',2(3x,e10.3))
200 format(3x,'x,psi,h=',3(3x,e10.3))
stop
end

function v(x,L,v0)
real L
if(x.lt.L)v=0.
if(x.ge.L)v=v0
return
end

 

 

END OF CHAPTER