Baretti LECTURES ON QUANTUM MECHANICS

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

Chapter 8 - The hydrogen atom

                        ∂ /∂x₂    =           ( m₂ /M) ∂/∂x_cm   +     ∂/∂r                                            .  (5)

Using (3) and (4) one gets

(1/m₁ ) ∂²  /∂x₁² = { (m₁/M²) (∂²  /∂x_cm² ) - (2/M) ∂² /∂x _cm ∂r +(1/m₁)∂² /∂r²  }          ,    (6)

(1/m₂ ) ∂²  /∂x₂² = { (m₂/M²) (∂²  /∂x_cm² ) + (2/M) ∂² /∂x _cm ∂r +(1/m₂)∂² /∂r²  }  .        (7)

Substitution of (6)-(7) in (3) gives two terms

 K_operator =    - (h'² /2 ) {(1/M ) ∂²  /∂x_cm²   +    (1/ μ )∂²  /∂r ² }          ,                     (8)

where μ = m₁ m₂ /(m₁ + m₂) is the reduced mass .  If m₁ >> m₂  ,  μ ≈ m₂ , the electron mass.

Let the total wave function be Ψ_total  = F(x_cm) Ψ(r) . The resultant wave equation is

(1/M ) Ψ(r)  ∂² F(x_cm) /∂x_cm²   +    F(x_cm)(1/ μ )∂² Ψ(r) /∂r ²

      + 2( E_total - U(r) ) F(x_cm) Ψ(r) = 0    .                                                            (9)

Let E_total = E_tranalation + E      and by separation of variables ,we have

                      (1/M ) ∂² F(x_cm) /∂x_cm²  = -2 E_translation F(x_cm )                        (10)

and the equation of interest is

                  (1/ μ )∂² Ψ(r) /∂r ²  = -2 ( E - U(r) ) Ψ(r)             .                       (11)

If we go to three dimensions , in spherical coordinates ,( r,θ, φ) eq. (11) goes over into

        + {(1/r² ) {∂² /∂ θ² + cotθ ∂/∂ θ +  (1/ sin²θ )  ∂² /∂φ² } Ψ ( r,θ, φ ) =   -2 μ (E - U(r) ) Ψ ( r,θ, φ )     .   (12)

It was shown in chapter 7 that the angular part has solutions

          -{∂² /∂ θ² + cotθ ∂/∂ θ +  (1/ sin²θ )  ∂² /∂φ² } Y_l,m (θ ,φ )

                  =    L² Y_l,m (θ ,φ )   = l(l+1) Y_l,m (θ ,φ )                   . (13)

Substitution of           Ψ ( r,θ, φ ) = R(r) Y_l,m (θ ,φ )     in (12) gives, 

             ( ∂² /∂r² +(2/r) ∂ /∂r ) R( r) - l(l+1) R(r) /r²    =     -2 μ (E - U(r) ) R (r)              .  (14)

We will solve (14) numerically.

The units and scales of magnitudes have to be established. Take the reduced mass μ =1 and as always h' (h bar)=1.

The coulombic potential energy is

                                                                 U(r)= -kZe²/r                                              .(15)

Set the factors       k Z e² =1 .

scale factors   the scale of length   L_s  is obtained using  (h'²/ μ ) (1/ L_s² ) ≡ kZe²/ L_s (see the Introduction chapter)

                                        L_s = {h'²/ (μ kZe²) }                                                    (16)

The scale of energy is obtained from    h'² / (μ L_s² ) = E_s . Therefore

                                                                        ^E_s = ke² /L_s                   .             (17)

 In the introduction we calculated to three digits these scales. but there we used m=9.11E-31 kg , whereas now we must employ

the reduced mass,

                                               μ = (m_e m_proton ) /( m_{e +} m_proton ) . Let m_proton = 1.67E-27 kg , then

                                                μ = 9.105E-31 kg   . 

So to three digits we may still employ the former calculations,

                           L_s   = 5.29E-011 meters ,     E_s = 4.36E-018  joules = 27.2 eV    .

The asymptotic form of  R(r)   is obtained from (14)  making   r→ ∞ .

         d ²R /dr² = -2 E R  = 2/E/ R ≡ k² R   . The physical solution must have -kr = -(2/E/)^1/2 r , as the argument of the exponetial.

Thus      

                                            lim R _{r→ ∞}  = exp {- (2/E/) ^1/2 r }                           . (17)

  It was shown in chapter 2, that the absolute value  of the n-th quantized state had a dependency

                      E_n ~  n^2α/(2+α) on the quantum number n, if the potential energy had the form U ~ r^α .    

For the coulombic potential α = -1  and   E_n ~ - 1/n² . The energy must be negative to have a bound state. Moreover from chapter 2 eq (9) ,

       E_n = -(1/2) (1/n²) . So

           lim R _{r→ ∞}  = exp {-  r /n }       where n is a quantum number equal to 1,2,3....

Our numerical results will justify such assertion.

Another limit to be examined is when r→0 .   

 As r→ 0    let    R~ r^s .The right hand side of (14) is of the order of  U r^s ~ r^s-1. Meanwhile the left side of the equation has three terms

proportional to r^s-2 .

They are  { s(s-1) +2 s - l (l+1) } r^s-2 .

The coefficient is equated to zero 

               s(s-1) +2 s - l (l+1) = 0 , and the solutions are s= l  , s = -(l +1) . Only the positive power , s= l ,can have physical meaning.

 One concludes that near the origin the wave function has the form  , R ~  r ^l  .

This leads to the following  IC at the origin ,

       R (0) =1     , d R(0)/dr = 0   , when l =0   ,             (18)

       R(0)=0       ,d R(0)/dr = 1   , if  l = 1

       R(0)=0       ,d R(0)/dr = 0   , if  l > 1 .

The difference equation solution to (14) is (setting μ = 1) ,

 R_n = 2 R_n-1 - R_n-2 + ( ∆r )² {  -(2/r_n-1) (R_n-1 -R_n-2 )/∆r + l(l+1) R_n-1 /( r_n-1)²

                                                                -2  (E - U(r_n-1 ) ) R_n-1   }     ,           (19)

where r_n = n ∆r .

    The following run has the energy running from -.6 au to -.0407 au.  Equation (19) is integrated ,for each energy value,well into the classically forbidden region  where E < U(r). The asymptotic sign given by Ψ(r_final) / abs( Ψ(r_final)    changes value when the equation goes through an energy eigenvalue.  From the run we see that the first three eigenvalues  with quantum numbers (l=0, m=0)  are,

                                                    -.505    <  E ₁ < -.496                              . (20)

                                                    -.126    <  E ₂ < -.117   

                                                     -.597    <  E ₃ < -.0502                                                                     

   This of course matches perfectly well with our expectation that E_n ~ -1/n² . We see that E_n = -1/(2n²) is the correct formula , for

E₁ = -.5    ; E₂ = -.125    , E₃ = -.0556 .

RUN

 enter atomic number
1.
atomic number, ang momentum 1. 0.

energy, psif/abs(psif)= -0.600E+00 0.100E+01
energy, psif/abs(psif)= -0.591E+00 0.100E+01
energy, psif/abs(psif)= -0.581E+00 0.100E+01
energy, psif/abs(psif)= -0.572E+00 0.100E+01
energy, psif/abs(psif)= -0.562E+00 0.100E+01
energy, psif/abs(psif)= -0.553E+00 0.100E+01
energy, psif/abs(psif)= -0.543E+00 0.100E+01
energy, psif/abs(psif)= -0.534E+00 0.100E+01
energy, psif/abs(psif)= -0.524E+00 0.100E+01
energy, psif/abs(psif)= -0.515E+00 0.100E+01

energy, psif/abs(psif)= -0.505E+00 0.100E+01
energy, psif/abs(psif)= -0.496E+00 -0.100E+01
 

energy, psif/abs(psif)= -0.486E+00 -0.100E+01
energy, psif/abs(psif)= -0.477E+00 -0.100E+01
energy, psif/abs(psif)= -0.467E+00 -0.100E+01
energy, psif/abs(psif)= -0.458E+00 -0.100E+01
energy, psif/abs(psif)= -0.448E+00 -0.100E+01
energy, psif/abs(psif)= -0.439E+00 -0.100E+01
energy, psif/abs(psif)= -0.429E+00 -0.100E+01
energy, psif/abs(psif)= -0.420E+00 -0.100E+01
energy, psif/abs(psif)= -0.410E+00 -0.100E+01
energy, psif/abs(psif)= -0.401E+00 -0.100E+01
energy, psif/abs(psif)= -0.391E+00 -0.100E+01
energy, psif/abs(psif)= -0.382E+00 -0.100E+01
energy, psif/abs(psif)= -0.373E+00 -0.100E+01
energy, psif/abs(psif)= -0.363E+00 -0.100E+01
energy, psif/abs(psif)= -0.354E+00 -0.100E+01
energy, psif/abs(psif)= -0.344E+00 -0.100E+01
energy, psif/abs(psif)= -0.335E+00 -0.100E+01
energy, psif/abs(psif)= -0.325E+00 -0.100E+01
energy, psif/abs(psif)= -0.316E+00 -0.100E+01
energy, psif/abs(psif)= -0.306E+00 -0.100E+01
energy, psif/abs(psif)= -0.297E+00 -0.100E+01
energy, psif/abs(psif)= -0.287E+00 -0.100E+01
energy, psif/abs(psif)= -0.278E+00 -0.100E+01
energy, psif/abs(psif)= -0.268E+00 -0.100E+01
energy, psif/abs(psif)= -0.259E+00 -0.100E+01
energy, psif/abs(psif)= -0.249E+00 -0.100E+01
energy, psif/abs(psif)= -0.240E+00 -0.100E+01
energy, psif/abs(psif)= -0.230E+00 -0.100E+01
energy, psif/abs(psif)= -0.221E+00 -0.100E+01
energy, psif/abs(psif)= -0.211E+00 -0.100E+01
energy, psif/abs(psif)= -0.202E+00 -0.100E+01
energy, psif/abs(psif)= -0.192E+00 -0.100E+01
energy, psif/abs(psif)= -0.183E+00 -0.100E+01
energy, psif/abs(psif)= -0.173E+00 -0.100E+01
energy, psif/abs(psif)= -0.164E+00 -0.100E+01
energy, psif/abs(psif)= -0.154E+00 -0.100E+01
energy, psif/abs(psif)= -0.145E+00 -0.100E+01
energy, psif/abs(psif)= -0.136E+00 -0.100E+01

energy, psif/abs(psif)= -0.126E+00 -0.100E+01
energy, psif/abs(psif)= -0.117E+00 0.100E+01
 

energy, psif/abs(psif)= -0.107E+00 0.100E+01
energy, psif/abs(psif)= -0.976E-01 0.100E+01
energy, psif/abs(psif)= -0.881E-01 0.100E+01
energy, psif/abs(psif)= -0.786E-01 0.100E+01
energy, psif/abs(psif)= -0.692E-01 0.100E+01

energy, psif/abs(psif)= -0.597E-01 0.100E+01
energy, psif/abs(psif)= -0.502E-01 -0.100E+01

energy, psif/abs(psif)= -0.407E-01 -0.100E+01

more ? =1 another run , =0 no more

********

If we run the program again with a very small spread in the energy say , from -.49 to -.51 , the exact energy can be pinpointed.

Here the number of steps - nstep- was increased to 5000 and  xlim=3.5*Z/abs(E).

 energy, psif/abs(psif)= -0.498571E+00 -0.100000E+01
energy, psif/abs(psif)= -0.498857E+00 -0.100000E+01
energy, psif/abs(psif)= -0.499143E+00 0.100000E+01

The estimate of E is now  -.49895 au .

Fig 2.  Hydrogenic ground state wave function. It is a decaying exponential.

We write for the ground state    Ψ_{n=1, l = 0 ,m=0}  = N exp(-r)   . It has to be normalized integrating its square over all space

∫  Ψ ² r² dr dΩ =1     , where the solid angle differential is  dΩ = sin(θ) dθ dφ   .  This gives     

                            N²  4 π ∫  exp(-2r) r² dr = N² π =1 ,

and therefore      N = 1/π^1/2 . The normalized ground state hydrogenic wave function is

                                               Ψ_{1, 0 ,0}  = (1/ π^1/2 ) exp(-r)                      . (21)

It may well happen that the choice of unit length may not include the factor Z ,like the length scale in (16).

Then the distance radial distance say u =h'²/ (μ ke²) }  is to our r  as

                                                u / r = Z                      

Then  exp(-u) = exp(-Zr)    ,  dr = du/Z     , r ² = (u²/Z²) .

The normalization becomes   N² ∫ exp(-2Zr) 4π(r²) dr    =N² ∫ exp(-2u) 4π(u²/Z²) du/Z = 1;

                                    or    N² π/Z³ = 1   , thus   N= (Z³ / π)^1/2 .

Eq (21) may be written as   Ψ_{1, 0 ,0}  = (Z³/ π )^1/2 exp(-Zr)  with the understanding that the definition of the scale of length does not include the atomic number Z.

FORTRAN # 8.1 - code for hydrogenic wave functions

c Hydrogenic atom 8 de octubre 2009
equivalence (dx,h)
dimension psi(0:10000) ,psinf(300),energy(300)
v(x)= -Z/x
g(x,i)=-2.*(e- v(x))*psi(i)+
$ al*(al+1.)*psi(i)/x**2
135 print*,' enter atomic number '
read*,Z
e=-.6*Z**2
c input angular momentum quantum number = al
al=0.
efinal= -z**2/(2.*4.**2)
nstep=1000
niter=60
de=(efinal-e)/float(niter)
c escale=27.2109
escale=1.
c print*,'escale,e=',escale,e
print*,'atomic number, ang momentum',z,al
print*,' '
do 110 iter=1,niter
xlim=3.0*Z/abs(e)
dx=xlim/float(nstep)
c xlim should go beyond the classical turning point where e=V(x)
c initial condtions psi(x)=x**al , psi'(x)= al*x**(al-1.)
c for al=0. psi(0)=1. psi(1)=psi(0)-dx
if(al.eq.0.)then
psi(0)=1.
psi(1)=psi(0)-dx
endif
c for al=1. psi(0)=0. psi'(0)=1. , psi(1) = psi(0)+dx
if(al.eq.1.)then
psi(0)=0.
psi(1)=psi(0)+dx
endif
c for al=2. psi(0)=0. psi(1)= (1./2.)*dx**2*al
if(al.eq.2.)then
psi(0)=0.
psi(1)= psi(0)+ (1./2.)*dx**2*al
endif
do 100 i=2,nstep
x=dx*float(i)
psi(i)=dx**2*g(x-dx,i-1)+2.*psi(i-1)-psi(i-2)-(2.*dx/(x-dx))*
$(psi(i-1)-psi(i-2))
c print*,'e,g,psi(i)=',e,g(x-dx,i-1),psi(i)
100 continue
energy(iter)=e
psinf(iter)=psi(nstep)
e=e + de
110 continue
do 20 i=1,niter
print 120,energy(i)*escale,psinf(i)/abs(psinf(i))
20 continue
120 format(3x,'energy, psif/abs(psif)=',2(4x,e10.3))
print*,' '
c call plotwave(psi,nstep,50,dx)
print*,'more ? =1 another run , =0 no more'
more=0
c read*,more
if(more-0)130,130,135
130 stop
end


subroutine plotwave(psi,nstep,ns,dx)
dimension psi(0:5000)
do 10 i=0,nstep,ns
x=dx*float(i)
print 100 , x, psi(i)
10 continue
100 format(2x,'x,psi=', 2(4x,e11.4))
return
end


 

Second eigenfunction

Set E= -(1/2) (1/2² ) = -1/8 in the above code with psi(0)=1. , psi(0+∆x) = 1.

Fig    3.   Ψ _{2, 0,0} =   Psi2s(r) = N ( 2-r)exp(-r/2).

To find  the normalization constant integrate use again the definition  ∫ { Ψ _{2, 0,0} }² 4π r² dr  =1 .

A useful definite integral ,for this calculations , is  ∫₀^∞  rⁿ exp(-ar) dr = Γ(n+1)/ aⁿ⁺¹ , or use a symbolic language e. g.

SAGE code

var('r');
f=( 2-r)*exp(-r/2);
integral(4*pi*r^2*f^2,r,0,oo)

32*pi .

 This gives    N² 32π = 1  , hence    N = (1/(32π ))^1/2  and

                   Ψ _{2, 0,0} =  (1/(32π ))^1/2 ( 2-r)exp(-r/2)                        .          (22)

Or one can select another choice of scale then as explained above , r → Z r , N → Z^3/2 N  , making

                                       Ψ _{2, 0,0} =  (Z³/(32π ))^1/2 ( 2-Zr)exp(-Zr/2)  .     (23)

To construct Ψ_3,0,0  we may set E = -(1/2) (1/3² ) in our code and plot R(r) finding the roots, just as we did with Ψ _{2, 0,0} .

However one may also obtain Ψ_3,0,0  using the orthogonality conditions.

Let    the un normalized  Ψ_3,0,0  be the polynomial of order (3-1) ,(with two constants a, b to be determined) ,

                                               Ψ_3,0,0  = (a + br + r² ) exp(-r/3)     .             (24)

Then we get two equations for a and b ,

       ∫ Ψ_3,0,0 Ψ _{1, 0,0} r² dr  = 0       ,   (25)

var('r');f1=exp(-r);
f2=( 2-r)*exp(-r/2);f3=(a+b*r+r^2)*exp(-r/3);
integral(r^2*f1*f3,r)

The indefinite integral is ,

I _1,3 =

-3/32*(8*r^2 + 12*r + 9)*a*e^(-4/3*r) - 3/128*(32*r^3 + 72*r^2 + 108*r +
81)*b*e^(-4/3*r) - 3/128*(32*r^4 + 96*r^3 + 216*r^2 + 324*r +
243)*e^(-4/3*r)

When evaluated at infinity  it is zero and at the lower limit r=0 , the only contributing terms are those not proportional to powers of r, with the signs changed ,

           I _1,3 =      (27/32)a + (243/128)b + 729/128 = 0

I _2,3    =    ∫ Ψ_3,0,0 Ψ _{2, 0,0} r² dr = 0                                                     .  (26)

c coeff hydrogenic radial wave functions
real lambda
dimension a(0:10)
a(0)=1.
5 print*,'anprinc , al? '
read*,anprinc , al
print *,'radial polynomial is of degree anprinc-1=',anprinc-1.
print*,'leading term is a0*r**(l)'
print*,' '
ifinal=int(anprinc)-1
lambda=1./anprinc
do 10 i=1,ifinal
p=float(i-1)
a(i)= a(i-1)*(2.*lambda*(p+al+1.) -2.)/( (p+al+1.)*(p+al+2.)
$-al*(al+1.) )
10 continue
do 20 i=0,ifinal
print 100,i ,a(i)
20 continue
100 format('i,a(i)=',i3,3x,e12.5)
stop
end
 

END OF CHAPTER