Baretti LECTURES ON QUANTUM MECHANICS

by Reinaldo Baretti Machín   (reibaretti2004@yahoo.com)

  Free counter and web stats

Chapter 7 - Angular Momentum

Lectures on Quantum Mechanics -Introduction  

Chapter 1 - Some experimental facts

 Chapter 2 - The Old Quantum theory

Chapter 3 - The Postulates of Quantum Mechanics  

Chapter 4- Potential Barriers

Chapter 5 - Potential Wells

Chapter 6 -  Periodic Potentials

Chapter 7 - Angular momentum

Chapter 8 - The Hydrogen Atom

FREE DOWNLOAD FORTRAN-FORCE 2.0.8

 

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. Applied mathematics for engineers and physicists  , by Louis Albert Pipes

5.  Theoretical physics by A. S Kompaneets 

6. Classical Mechanics by Herbert Goldstein

 

Angular momentum in classical mechanics

It would be instructive to review the classical problem of a body moving in a central filed. Consider the gravitational problem of a particle of mass m  in a plane X-Y  (see Fig  ) moving under the influence of  a central potential  V(r)= -GM/r .  G is the universal gravitational constant ( G =6.67E-11Nm2/kg2 )  and M is the mass of the attractor e.g. M = 5.98E24 kg .

The equations of motion in classical mechanics require that  two initial conditions be given  , the position and the velocity. another option is the specification of the total energy ,E , and the angular momentum L.

The Lagrangian function , labeled here Lg , is the difference between the kinetic and potential energies (see ref. 5 -6)

                                 Lg = T - V  = (1/2) m ( vr2 + vt2 )  + GMm/r                         

                                                  = (1/2)m( r'2 + r2 θ'2) + GMm/r                    (1) 

where r'=dr/dt (the derivative of the magnitude of r)   ,   θ' =d θ/dt . 

The total energy is another constant of motion ,

                                       E = T + V =  (1/2) m ( vr2 + vt2 )  - GMm/r        .  (2)       

 

 

Fig 1. Motion in a plane under  gravitational attraction.

The equation of motion for each independent coordinate qi is

     d (  ∂Lg/∂q'i ) /dt  =  ∂Lg/∂q i                                                                                   (3)

where   q'i  = dqi/dt .

Let qi = θ . Notice that ∂Lg/∂θ = 0  . Therefore  d (  ∂Lg/∂ θ' ) /dt =0     and 

          ∂Lg/∂ θ' = constant in time with dimensions of Lg - energy -time~ angular momentum.

          ∂Lg/∂ θ' = mr2 θ' ≡ L (angular momentum ) = constant of the motion.                (4)

The vector L is also defined by the cross product , L = r x (m v) = r x p.

Substituting θ' = L /( mr2 ) in (2) gives the total energy in terms of  r and r' ,

  E = (m/2) r'2 + {L2/(2mr2 )  - GMm/r }                                 .                               (5)

E and L ( the two constants of the motion) can be set to arbitrary values independently. This is somewhat similar to setting r0 and v0 as initial arbitrary conditions.

The term in parenthesis is called the effective potential

                                       Vef (r) = L2/(2mr2 )  - GMm/r                                        (6)

and  L2/(2mr2 ) is the centrifugal potential . An example of Vef(r) is shown in Fig 2 in arbitrary units where L=1.The particle m can not fall to the centre due the repulsive wall  L2/(2mr2 ) .                                      

var('r'); L=1 ; f= (L^2/(2*r^2) -1/r) ;
y=plot(f,r,0.4,4);
show(y)

 

Fig 2. L=1 , with E ≈ -.3 the motion traces an ellipse. The two turning points are r1 ≈ 0.5 , r2 ≈ 3 .

 

The equation of motion for r is

               m d2r/dt2 = ∂Lg/ ∂r =  m r θ'2 - GMm/r2                          (7)

When θ' = L /( mr2 ) is substituted in (7) we get                    

 

                                               r'' =    L2/(m2 r3 ) - GMm/r2     .       (8)

One can solve it by finite differences

                                rn = 2rn-1 -rn-2 + (∆t)2 { L2/(m2 rn-1 3 ) - GM/rn-1 2   }    . (9)

 

 

 


             Fig 3.Initial conditions are r0 and  v0 =(v0 cos(α ) , v0 sin( α )  ) .  The total energy is Etotal = (1/2) mv02 -GMm/r0 and the angular momentum   L = m r0 v0 sin(α ). 

 

 

Fig 4. A circular path is obtained when  v0 = (G*M/r0)1/2 , with v0 perpendicular to r0  (α = π/2),  L = r0 m v0 ,

Another example is given with E(rel) =E r0/(G*mearth*m) = -.5 , but α = 3π/4 .

Fig 5. The starting point is x0=rearth, y0 =0  and  v0  makes an angle 3π/4 with the +X axis. The plot shows the center of attraction at the origin  (x=0,y=0). It is one focal point of the ellipse.

 

 

 

 

FORTRAN code

 

c motion under gravitational potential /
c Solution using d^2r/dt^2= L**2/(m**2*r**3) -G*M/r**2
implicit real*8(a-h,o-z)
real*8 Mearth,m,L
data rearth,Mearth,m/6.36d6,5.98d24,1.d0/
data G/6.67d-11/
v(r)=-G*mearth*m/r
pi=2.d0*dasin(1.d0)
c alfa inital velocity angle with respect to +X axis
alfa =pi/2.d0
r0=rearth
escale=G*mearth*m/r0
v0=dsqrt(G*mearth/r0)
L=r0*m*v0*sin(alfa)
tscale=r0/v0
dt=tscale/500.d0
tf=6.5d0*tscale
nstep=int(tf/dt)
kp=int(dfloat(nstep)/80.d0)
kount=kp
E=.5d0*m*v0**2 -G*mearth*m/r0
print*,'E/escale=',E/escale
theta0=0.d0
x0=r0*dcos(theta0)
y0=r0*dsin(theta0)
x1=x0+dt*v0*dcos(alfa)
y1=y0+dt*v0*dsin(alfa)
r1=dsqrt( x1**2 +y1**2)
theta1=theta0+dt*L/(m*r0**2)
print 100 ,0.d0 ,theta0,x0,y0
do 10 i=2,nstep
t=dt*dfloat(i)
r2=2.d0*r1-r0+dt**2*( L**2/(m**2*r1**3)-G*mearth/r1**2 )
theta2=theta1+dt*L/(m*r1**2)
if(i.eq.kount)then
x=r1*dcos(theta2)
y=r1*dsin(theta2)
print 100 ,t,theta2,x,y
c print*,' '
kount=kount+kp
endif
r0=r1
r1=r2
theta1=theta2
10 continue
100 format('t,theta,x,y=',4(3x,d10.3))
stop
end

 

  

 

The Laplacian in orthogonal curvilinear coordinates

From the energy equation (5) ,  E = (m/2) r'2 + {L2/(2mr2 )  - GMm/r }   it is not hard to divinate that the hamiltonian in a problem with spherical coordinates would be

                        H  =  pr2 /(2m)  L2/(2mr2 ) + V(r)                                         (10)

where      pr2   +  L2/r2 are differential operators . These two should equal

                                               -h'2 2 = pr2   +  L2/r2                               .  (11)

Not surprisingly the classical angular momentum has become a differential operator.

The laplacian operating on  Ψ , in spherical coordinates , is given by

2Ψ  = (1/r2 ) ∂ (r2 ∂Ψ /∂r )/∂r + [1/(r2 sin(θ) )] ∂ (sin θ ∂ Ψ/∂θ)/∂θ + [1/(r2 sin2 (θ))] ∂2 Ψ /∂φ2    . (12)

 

We will derive this expression for the laplacian operator by two different methods.

Recall first that the  laplacian operator in cartesian coordinates is obtained ,easy enough, from the "dot" of the divergence operating on a gradient.

The gradient of a scalar V(x,y,z) is a vector  (see Ref. 4) ,

           grad  V(x,y,z) = i ∂ V/∂x  + j ∂V /∂ y   +   k  ∂V /∂z                                  . (1)

Finding the divergence , which is the scalar product of del "dot" del ,  we have,

div (grad V) = del   ( i ∂ V/∂x  + j ∂V /∂ y   +   k  ∂V /∂z )

                              =  ( i ∂ V/∂x  + j ∂V /∂ y   +   k  ∂V /∂z )  (i ∂ V/∂x  + j ∂V /∂ y   +   k  ∂V /∂z ).

But the unit vectors   i ,j k are independent of x,y,z  and they are orthonormal ,  i ∙ i =1    ,  i ∙j = 0 ,i ∙ k = 0 , etc.

So the  the Laplacian operator in cartesian coordinates is

                   ∆ V = ∂2V/∂x2  +  2 V /∂ y2   +    2 V /∂z2                           .    (2)     

 

  In spherical coordinates the unit vectors are labeled r u , θ u ,   φu  .They are shown in Fig 6. Their components along X , Y , Z and the dependency on the angles is also given.

 

Fig 6. Orientation of spherical unit vectors  runit , θ unit  ,φ unit  and their expressions in terms of the cartesian unit vectors  i ,j k.A few  derivatives of the unit vectors with respect  to θ and φ    are necessary

runit /∂ θ =  i cosθ cosφ+ j cosθ sin φ   - k sin θ = θunit   ,        (3)  

∂  runit /∂φ = - i sinθ sin φ+ j sin θ cos φ = sin(θ) φ unit      (4)

θ unit /∂ θ = -i sinθ cosφ- j sinθ sin φ - k cosθ = - runit     ,           (5)  

θ unit  /∂ φ = (-i sin φ   + j cosφ  ) cosθ  = cos (θ)  φ unit   ,  (6)

φ unit /∂ θ =  0           ,                                                 (7)

 ∂φ unit /∂ φ =  -cosφ -j sinφ         .                               (8)                          

 

To get the laplacian in spherical coordinates , start with the gradient

                    ∆  = runit  ∂/∂r +     unit 1/r) ∂/∂ θ  +   φ unit /(r sinθ ) ∂ /∂φ          .  (9)

The dot product

∆∙ ∆ = ( runit  ∂/∂r +   unit 1/r) ∂/∂ θ  +  φ unit /(r sinθ ) ∂ /∂φ) ∙ (runit  ∂/∂r +   unit 1/r) ∂/∂ θ  +  φ unit /(r sinθ ) ∂ /∂φ ) (10)

has in principle 9 terms., but most of them are zero.

Now bear in mind that the unit vectors are orthonromal and the relations (3) - (8) ; then

runit  ∂/∂r { runit  ∂/∂r +   unit 1/r) ∂/∂ θ  +  φ unit /(r sinθ ) ∂ /∂φ}  =  ∂2 /∂r2  ,       (10)

unit 1/r) ∂/∂ θ { runit  ∂/∂r  +   unit 1/r) ∂/∂ θ  +  φ unit /(r sinθ ) ∂ /∂φ}

= (1/r) ∂ /∂r  + ((1/r2 )2 /∂ θ2    ,                                                                       (11)

( φ unit /(r sinθ ) ∂ /∂φ)∙{ runit  ∂/∂r  +   unit 1/r) ∂/∂ θ  +  φ unit /(r sinθ ) ∂ /∂φ} =

=  (1/r) ∂ /∂r  + (cotθ /r2 )∂/∂ θ + (1/(r2 sin2 θ) ) ∂2 /∂φ2                         .            (12)

Adding (10),(11) and (12) gives the laplacian in spherical coordinates.

2 = ∂2 /∂r2 +(2/r) ∂ /∂r + (1/r2 ) {∂2 /∂ θ2 + cotθ ∂/∂ θ(1/ sin2θ )  ∂2 /∂φ2 } .    (13)

Another method employs Gauss divergence theorem in conjunction with orthogonal curvilinear coordinates (Ref. 4)

The distance squared  between two points in orthogonal curvilinear coordinates is given by

 

                        ds2 =  (h1 du1)2  + ( h2 du2 )2  + (  h3 du3)2                     .   (13)

For example in spherical coordinates ,  du1=dr    ,  du2 = dθ  , du3= dφ  and

                      ds2 =   dr2  + r2 dθ2  +  (r sinθ )2 dφ2                            .     (14)

An elementary volume element is

                                         ∆τ =  h1 h2 h3 du1du2 du3                           .        (15)

 

Comparing with(13) ity follows that    the factors  h1 =1   , h2 = r  , h3 =  r sinθ

Let the vector A  be equal to the gradient of the scalar function  V(u1 , u2 ,u 3)  and write

       A     e1 (1/h1) ∂V/∂u1     +   e2 (1/h2) ∂V/∂u2  +  e3 (1/h3) ∂V/∂u3       

             =         e1 A1       +  e2 A2  +        e3 A3                                               .  (16)

where e1  ,  e2  ,  e3  are units vectors along the coordinates  u1 ,  u2  and u3 .

Of course in cartesian coordinates,  e1= i  ,  e2= j  ,  e3 = k , u1 = x,  u2 = y  ,  u3 = z .

Consider the definition of  Gauss divergence theorem (see ref. 4)

                            div A =   lim ∆τ → 0  (  ∫ A ∙ d S ) /∆τ                      .    (17)

The flux along the u1 direction is made of two contributions ,

 ∆Φ1 = -A1 h2 h3 du2 du3  + { A1 h2 h3 du2 du3  + [∂ (A1 h2 h3) /∂u1] du1 du2 du3 }

         = [∂ (A1 h2 h3) /∂u1] du1 du2 du3                                                          . (18)

 

 Dividing by ∆τ we get  

         ∆Φ1 /∆τ =  ( 1/(h1 h2 h3 ) ) {∂ (A1 h2 h3) /∂u1}                   .       (19)

 

Introduce now  A1 in (19) ,

∆Φ1 /∆τ  =  ( 1/(h1 h2 h3 ) ) {∂ [ (h2 h3 / h1)  ∂V/∂u1 ] /∂u1}          .  (20)

The other terms are obtained from (20)  by permuting the indices .Therefore the laplacian operator applied to V is

 

2 V =(1/(h1 h2 h3 ) ){ ∂ [ (h2 h3 / h1)  ∂V/∂u1 ] /∂u1  + {∂ [ (h1 h3 / h2)  ∂V/∂u2 ] /∂u2}

 

                                                             + {∂ [ (h1 h2 / h3)  ∂V/∂u3 ] /∂u3}}          .     (21)

Example :

a ) spherical coordinates . spherical coordinates Fig 7 Spherical coordinates

 

We have  h1 =1   , h2 = r  , h3 =  r sinθ .

The first term in (21) is  (1/(r2 sin(θ)) ∂ [ (r2 sin(θ) )  ∂V/∂r ] /∂r , which gives

2 V/∂r2  +(2/r)∂V/∂r  . 

The second term is (1/(r2 sin(θ)) ∂ [ ( sin(θ) )  ∂V/∂θ ] /∂θ , which equals

(1/r2) ∂2 V/∂θ2  + (cot(θ)/r2 ) ∂V/∂θ   .

The last term is (1/(r2 sin(θ)) ∂ [ (1/ sin(θ) )  ∂V/∂φ ] /∂φ  , or

                          {1/(r sin(θ))2 }∂2 V/∂φ2       .

 

b) cylindrical coordinates :  

Cylindrical Coordinates Fig 8. Cylindrical coordinates.

Beware of the coordinates labels. In this figure  r2 = x2 + y2 and θ corresponds to φ in the spherical

coordinates drawing.

ds2 =   dr2  + r22  +  dz2  with   h1 =1  ,  h2 =r ,   h3 =1 . From  eq (21)

2 V =(1/(h1 h2 h3 ) ){ ∂ [ (h2 h3 / h1)  ∂V/∂u1 ] /∂u1  + {∂ [ (h1 h3 / h2)  ∂V/∂u2 ] /∂u2}   + {∂ [ (h1 h2 / h3)  ∂V/∂u3 ] /∂u3}}

2 V =(1/r) { ∂ [ r  ∂V/∂r ] /∂r  + ∂ [ (1/r) ∂V/∂θ] /∂θ + ∂ [ r  ∂V/∂ z] /∂z  }

       = ∂2 V/∂r2 + (1/r)∂V/∂r + (1/r2) ∂ 2V/∂θ2 + ∂ 2V/∂ z2 .

 

Returning to Physics, the kinetic energy operator in spherical coordinates has the form

(-h'2/(2m)) *( ∆2 ) = 1/(2m)) *(-h'2  ∂2 /∂r2 +(2/r) ∂ /∂r ) 

                              +((-h'2/(2m))*{(1/r2 ) {∂2 /∂ θ2 + cotθ ∂/∂ θ(1/ sin2θ )  ∂2 /∂φ2 } . (22)

The  angular momentum operator is defined by its square

 

        L2 = -h'2 {∂2 /∂ θ2 + cotθ ∂/∂ θ(1/ sin2θ )  ∂2 /∂φ2 }                                       .  (23)

To find the characteristic  functions of L2  let  ,

Y(θ ,φ ) satisfy 

{∂2 /∂ θ2 + cotθ ∂/∂ θ(1/ sin2θ )  ∂2 /∂φ2 }Y(θ ,φ ) = -C Y(θ ,φ )                   ,  (24)

where -C is an eigenvalue. Assuming separation of variables  Y(θ ,φ ) = Θ(θ )  Φ(φ) ,one has from (24),

 

 (1/Θ) sin(θ ) d {sin(θ) dΘ/dθ }/dθ + C (sinθ )2 + (1/Φ) d2 Φ/dφ2 = 0          .    (25)

The first step in finding the solution is to set   (1/Φ) d2 Φ/dφ2  = -m2 , m=0, +- 1 ,+-2, .. i.e.  Φ( φ) =  A exp(-i mφ) . The separation constant m must

be an integer since   Φ( φ + 2π ) = Φ( φ ) ,thus  m2π has to be an integral number of revolutions.

We can identify the z component of the angular momentum Lz with the operators

 

 

                                                                           Lz = -ih' ∂ /∂φ                   ,   (25)

such that                                             Lz2 Y ( θ,φ )  = m2h'2 Y.                        (26)

  

Eq (25) is now

(1/Θ) sin(θ ) d {sin(θ) dΘ/dθ }/dθ + C (sinθ )2  - m2 = 0  .                                (27)

A change of variable is introduced  with  x = cos(θ)  and label Θ = y(x) . We have dx/dθ = - sin(θ) = -(1-x2 )1/2 .

So   dΘ/dθ =dy/dθ = (dy/dx)(dx/dθ) = -(1-x2 )1/2 (dy/dx)                                      ,      (28)

and

d2 Θ/dθ2  =d2 y/dθ2 = (1-x2) (d2 y/dx2 ) - x (dy/dx)                                    .  (29)

 

Eq(26) can be cast in the form (see ref. 8),

(1- x2 ) y'' -2x y ' + { C - m2/(1-x2) } y  =  0 , where y'=dy/dx , y''= d2 y/dx2   .    (30)

The range for x is  -1 ≤  x ≤ +1   since x= cos(θ)  and     0 ≤  θ ≤ π  . 

Set m=0 initially in (29) . Letting  y = xs in (29)  near the origin gives s(s-1) xs-2 = 0 , implying either  (s=0 even solutions) , (s = 1 ,odd solutions) .

Though we seek a polynomial solution for y(x), we start with a series solution ,   y = n=0 an xn . The derivatives are y' = ∑n=0 n an xn-1   , y '' =  n=0 n(n-1) an xn-2 . Substituting in (29) with m=0 gives for the x n power  gives

                                               [ (n+2)(n+1)an+2 +{ C - n(n+1) } an =  0  ,

producing the rercusion relation 

                                               an+2 =  - { (C- n(n+1) )/( (n+2)(n+1) ) } an              .   (31) 

A polynomial of order  n= l  is generated by equating C = l ( l+1) , then

 

 an+2 =  - { (l(l+1)- n(n+1) )/( (n+2)(n+1) ) } an                                               ,       (32)

 

and   an+2 = 0 .

The solutions of  

(1- x2 ) y'' -2x y ' + { l(l+1) } y  =  0         are the Legendre polynomials labeled Pl (x) , with an l subscript ,the denoting the order of the polynomial. They are normalized to have the value  unity at x= +1.

Example:  Find the first three polynomials.

Let l=0  and a0 =1 .From (31)  a2 = -0 a0 = 0 and  P0(x) = 1.

Let l=1  and a1=1 . From (31)  a3 = -0 a1 =0  and  P1(x) =a1x = x .

Let l=2 and a0 ≠ 0 . From (31)  a2 =  -(6/2) a0 = - 3a0  . The un normalized polynomial is  P2(x) = (-3x2 +1)a0 .

Imposing  P2(x=1) =1 , gives a0 = -1/2  and the normalized polynomial is

P2(x) = (1/2) ( 3x2 - 1)  .

Example : Integrate  (1- x2 ) y'' -2x y ' + { l(l+1) } y  =  0  , numerically for l=5.

Let the inital conditions be y(x=1) =1. From the DE itself  , ư'(1)= (l(l+1))y(1)/(2 xi) = (1/2)l(l+1).

Writing the DE as a difference equation we get

                  yn = 2 yn-1 - yn-2 + [(∆x)2 /(1-x2 n-1 ) ]{  2(xn-1) ( yn-1 - yn-2 )/∆x - l (l+1) yn-1 }       .   (33)

 

 

Fig 9.  P5(x)  obtained numerically employing the Fortran code.

FORTRAN code

c Legendre polynomials
real L
data L /5.0/
data p0,xi,xf,nstep/1.0,1.0,-1.0,6000/
f(x,p0,p1)=(1.0/(1.0-x**2))*( 2.0*x*(p1-p0)/dx -L*(L+1.0)*p1)
dydx= (L*(L+1.0))*p0/(2.0*xi)
xi=1.0
xf=-1.0
nstep=1000
kp=int(float(nstep)/100.0)
kount=kp
dx=(xf-xi)/float(nstep)
p1=p0+dx*dydx
print*,' L=',L
print*,' '
print 100, xi, p0
do 10 i=2,nstep
x=xi+dx*float(i)
p2=2.*p1-p0+dx**2*f(x-dx,p0,p1)
if(i.eq.kount)then
print 100, x, p2
kount=kount+kp
endif
p0=p1
p1=p2
10 continue
print 100, x, p2
100 format(1x,'x,P(x)=',2(4x,e10.3))
stop
end


We direct now our attention toward the solutions of 

                             (1- x2 ) y'' -2x y ' + { l(l+1) - m2/(1-x2) } y  =  0                                       (33)

which contains two eigenvalues  l and m . They are called associated Legendre functions and labeled  Plm (x) , with the m superscript.

 

From (23) and (24) after identifying C and m we have that the eigenfunctions of the angular momentum satisfy

 -h'2 {∂2 /∂ θ2 + cotθ ∂/∂ θ(1/ sin2θ )  ∂2 /∂φ2 } Yl m (θ,φ) =  L2 Yl m (θ,φ) = h'2 l(l+1) Yl m (θ,φ)    .   (34)

In units of h'    the average 

              <  L2    > = <  Lx2 >   + < Ly2   >  + < Lz2   >  = <  Lx2 >   + < Ly2   >  + m2   ,

or                    l(l+1) - m2 ≥ 0 implying       abs(m) ≤ l  .

For a given l   ,   m is allowed 2l+1 values,  from +l , l-1, l-2,...0 ,-1, ...-l   .

 

The analytic solution is given by (see ref. 8) ,

 P_\ell^{m}(x) = \frac{(-1)^m}{2^\ell \ell!} (1-x^2)^{m/2}\  \frac{d^{\ell+m}}{dx^{\ell+m}}(x^2-1)^\ell.      for m ≥ 0   .                                     (35)

For negative m 

 


P^{-m}_\ell(x) = (-1)^m \frac{(\ell-m)!}{(\ell+m)!} P^{m}_\ell(x).
                                                                                        . (36)

Employing Dirac notation, the angular momentum L2 and Lz  operate on / l , m> ,

 

 L^2 | l, m \rang = {\hbar}^2 l(l+1) | l, m \rang
 L_z | l, m \rang = \hbar m | l, m \rang

where

 \lang \theta , \phi | l, m \rang = Y_{l,m}(\theta,\phi)   .

 

 

 

Reverting to our numerical approach if (33)is multiplied by (1-x2 )  as x→ 1 this leaves  -m2 y(1) =0 ,  implying y(1)=0.

Using the FORTRAN code  (Associated Legendre functions) we solve (33) numerically to get the un normalized associated Legendre functions with L=4 ,

m=1,2,3,4.

 

Fig 10 . Plot P41 (x)   obtained numerically.     IC  y(0)=0 , y'(0) =1 ,   L=4  , m=1.

 

 

 

Fig 11. P42 (x)   obtained numerically.     IC  y(0)=1 , y'(0) =0 ,   L=4  , m=2.

 

Fig 12.  Plot P43 (x)   obtained numerically.     IC  y(0)=0 , y'(0) =1 ,   L=4  , m=3.

 

 

 

Fig 13.  Plot P44 (x)   obtained numerically.     IC  y(0)=1 , y'(0) =0 ,   L=4  , m=4.

From the DE ( eq(33)) we can not discriminate numerically  between positive and negative m.




FORTRAN code (Associated Legendre functions)

c associated legendre functions (explore finding the value of C )
c DE (1-x**2)*y''-2.*x* y' +(L(L+1.) - m**2/(1.-x**2) ) y=0.
c parity L-m ~even , L-m ~ odd
implicit real*8(a-h,o-z)
real*8 m ,L
data nx ,xi,xf/2000,0.d0,1.d0/
data L, m/4.d0,4.d0/
data ip/1/
f(x,y0,y1)= (1.d0/(1.d0-x**2))*( 2.d0*x*(y1-y0)/dx -( L*(L+1.d0)
$ -m**2/(1.d0-x**2))*y1)
dx=(xf-xi)/dfloat(nx)
kp=int(dfloat(nx)/60.d0)
kount=kp
c inital conditions y0=1., dydx=0. or y0=0,dydx=1. parity ~ (L-m)
y0=1.d0
dydx=0.d0
y1=y0+dx*dydx
print*,'L,m ,y0,dydx=',L , m , y0, dydx
print*,' '
if(ip.eq.1)print 110 ,xi,y0
do 10 i=2,nx
x=xi+dx*dfloat(i)
y2=2.d0*y1-y0+dx**2*f(x-dx,y0,y1)
if(i.eq.kount)then
if(ip.eq.1)print 110 ,x,y2
kount=kount+kp
endif
y0=y1
y1=y2
if(i.eq.nx)print 110 ,x,y2
10 continue
110 format('x,y=',2(4x,d10.3))
stop
end
 

The associated Legendre function -squared -is used to represent in rather qualitative fashion the spatial configurations of the complete

wavefunction.  Let the atomic wavefunction be of the form Ψ (r, θ ,φ ) ~ R(r) Plm (θ) exp(i mφ) . Then

Ψ Ψ * ~ R(r)2 (Plm (θ) )2 .

In fig 14 we reproduce the pz orbital. We take z coordinate = (P10 (θ) )2 and the X-Y plane radius to be

(P10 (θ) )2 tan ( θ ). 

 

Fig 14. A plot of l=1 ,m=0 or  pz orbital.  Notice the different scales employed in the Z axis and the X-Y plane.

 

 

Fig 15. Plot of l=1 ,m=1 or -1 . The x-y  plane value ≡ {P11 (θ) }2 , zcoordinate = {P11 (θ) }2 /tan (θ ).

 

A more artistic rendition of  orbitals is available on many sites for example this one

taken from  http://boomeria.org/chemlectures/ch9orbitals1.jpg . The symbols s, p, d, f correspond to l=0,1,2,3.

The drawings are a function of    ( Plm (θ) )2 .

 

FORTRAN code for the orbitals

c orbitals (P l,m (theta))^2 phi^2 = Plm^2
data thetai, ntheta/0.,40/
p10(u)=u
P11(u)=-sqrt(1.-u**2)
pi=2.*asin(1.)
thetaf=pi
dtheta=(thetaf-thetai)/float(ntheta)
do 10 i=0,ntheta
theta=dtheta*float(i)
u=cos(theta)
if(theta.ne.pi/2.)x=p10(u)**2*tan(theta)
print 100 ,x ,p10(u)**2
print 100 ,-x ,p10(u)**2
print 100 ,-x ,-p10(u)**2
print 100 ,x ,-p10(u)**2
10 continue
100 format('x,P10*2=',2(4x,e10.3))
stop
end

 

Table of associated Legendre functions

P_{0}^{0}(x)=1
P_{1}^{-1}(x)=-\begin{matrix}\frac{1}{2}\end{matrix}P_{1}^{1}(x)
P_{1}^{0}(x)=x
P_{1}^{1}(x)=-(1-x^2)^{1/2}
P_{2}^{-2}(x)=\begin{matrix}\frac{1}{24}\end{matrix}P_{2}^{2}(x)
P_{2}^{-1}(x)=-\begin{matrix}\frac{1}{6}\end{matrix}P_{2}^{1}(x)
P_{2}^{0}(x)=\begin{matrix}\frac{1}{2}\end{matrix}(3x^{2}-1)
P_{2}^{1}(x)=-3x(1-x^2)^{1/2}
P_{2}^{2}(x)=3(1-x^2)
P_{3}^{-3}(x)=-\begin{matrix}\frac{1}{720}\end{matrix}P_{3}^{3}(x)
P_{3}^{-2}(x)=\begin{matrix}\frac{1}{120}\end{matrix}P_{3}^{2}(x)
P_{3}^{-1}(x)=-\begin{matrix}\frac{1}{12}\end{matrix}P_{3}^{1}(x)
P_{3}^{0}(x)=\begin{matrix}\frac{1}{2}\end{matrix}(5x^3-3x)
P_{3}^{1}(x)=-\begin{matrix}\frac{3}{2}\end{matrix}(5x^{2}-1)(1-x^2)^{1/2}
P_{3}^{2}(x)=15x(1-x^2)
P_{3}^{3}(x)=-15(1-x^2)^{3/2}
P_{4}^{-4}(x)=\begin{matrix}\frac{1}{40320}\end{matrix}P_{4}^{4}(x)
P_{4}^{-3}(x)=-\begin{matrix}\frac{1}{5040}\end{matrix}P_{4}^{3}(x)
P_{4}^{-2}(x)=\begin{matrix}\frac{1}{360}\end{matrix}P_{4}^{2}(x)
P_{4}^{-1}(x)=-\begin{matrix}\frac{1}{20}\end{matrix}P_{4}^{1}(x)
P_{4}^{0}(x)=\begin{matrix}\frac{1}{8}\end{matrix}(35x^{4}-30x^{2}+3)
P_{4}^{1}(x)=-\begin{matrix}\frac{5}{2}\end{matrix}(7x^3-3x)(1-x^2)^{1/2}
P_{4}^{2}(x)=\begin{matrix}\frac{15}{2}\end{matrix}(7x^2-1)(1-x^2)
P_{4}^{3}(x)= - 105x(1-x^2)^{3/2}
P_{4}^{4}(x)=105(1-x^2)^{2}