Fig 1. A helium like atom with nuclear charge + Ze is an atom that has been stripped of all but two electrons.
The hamiltonian is ,
H = -(1/2m) ∆ 1 -(1/2m) ∆ 2 - Z k e2 /r1 - Z k e2 /r1 + k e2 /r12 + spin-orbit terms , ( 1 -a )
where ∆ 1 is the laplacian operating on the coordinates r1 ,∆ 2 operates on the coordinates of the second electron r2 .
Neglecting the spin-orbit terms (they are zero in the ground state since L=0 )
H = -(1/2)(p1)2 -(1/2)(p2)2 - Z /r1 - Z /r2 + 1 /r12 . (1-b)
As discussed in the introduction, the units of length and energy are (see the Introduction and chapter 8 of these Lectures),
Ls = 5.29E-011 meters , Es = 4.36E-018 joules = 27.2 eV .
A most simple approximation to the ground state writes the wave function as the product of two hydrogenic ground state wave functions ,
Ψ (r1 , r2) = φ1s (r1) φ1s ( r2)
= (Z3 /π) 1/2 exp(-Zr1) (Z3 /π) 1/2 exp(-Zr2) . (2)
Ψ (r1 , r2) is normalized since ∫ φ2 1s (r1) φ2 1s ( r2) dτ1dτ2 = 1 where dτ = 4π r2 dr .
Taking the average involves the integration
Eground state = ∫ Ψ H Ψ dτ1dτ2 . (3)
This separates in
E = ∫ φ1s (r1){-(1/2)(p1)2 - Z /r1 }φ 1s (r1) dτ1 + ∫ φ1s (r2){-(1/2)(p2)2 - Z /r2 }φ 1s (r2) dτ2
+ ∫ (1/r12) φ2 1s (r1) φ2 1s ( r2) dτ1dτ2 . (4)
The first two terms are the ground state energy of a hydrogenic atom with atomic number Z , i.e.
∫ φ1s (r1){-(1/2)(p1)2 - Z /r1 }φ 1s (r1) dτ1 = -Z2 /2 and also ∫ φ1s (r2){-(1/2)(p2)2 - Z /r2 }φ 1s (r2) dτ2 =-Z2 /2 , then
E = -Z2 + J (1s,1s) , (5)
where the coulombic repulsion between the two 1s electrons is,
J(1s,1s) = ∫ (1/r12) φ2 1s (r1) φ2 1s ( r2) dτ1dτ2 . (6)
Before calculating J we make a guess of its magnitude. Let the absolute values of the electron average positions be / r1 / =/ r2 / ~1/Z . Imagine that the elctrons are at opposite sides so that / r1 - r 2/ ~2/Z .
So the total energy is of the the order of E ~ -Z2 + 2/Z . If Z=2 , E ~ -3 au = -81.6 eV . We are not too far away ,since helium experimental ground state energy
is -2.90 au.
One can proceed to calculate (6) by dividing it in two parts. Define the potential V a(r1) and Vb (r1) by the integrals
V a(r1) = (1/r1) ∫0 r1 4π (r2 )2 φ2 1s ( r2) dr2 r1 ≥ r2 , (7)
Vb (r1) = ∫∞ r1 4π r2 φ2 1s ( r2) dr2 r1 ≤ r2 . (8)
The total potential is V(r1) = V a(r1) + Vb (r1) and the repulsion energy J simplifies to
J(1s,1s) = ∫ φ2 1s (r1) V(r1) dτ1 . (9)
Eq.(7) reflects the fact that outside a shell of radius r2 ,the potential at a distance r1 from the center is equal to qshell / r1 .
Eq.(8) reflects the fact that inside a shell of radius r2 , the potential at a distance r1 from the center is equal to qshell / r2 .
We have ,
Va (r1) = 1/r1 -( 2r1 Z2 + 2Z +1/r1) exp(-2r1 Z) , (10)
Vb(r1) = Z (2 r1 Z + 1) exp( - 2 r1 Z ) , (11)
and
V(r1) = 1/r1 - ( Z +1/r1 ) exp(-2Zr1) . (12)
Substitution of (10) -(11) in (9) gives ,
J(1s,1s) = ∫ φ2 1s (r1) { Va (r1) + Vb (r1 ) } dτ1 = 5Z/8 . (13)
E(Z) = - Z2 + (5/8)Z ~ au . (14-a)
For helium E(Z=2) = -2.75 au , (14-b)
(1au =27.2 eV) . The experimental value is -2.90 au while an important method,called Hartree -Fock gives -2.8617 au ( see ref. 8).
MAXIMA code
(%i1)
assume(r1>0);assume(r2>0);assume(Z>0);
(%o1) [r1 > 0]
(%o2) [r2 > 0]
(%o3) [Z > 0]
(%i4)
phi1(r1):=(Z^3/%pi)^(1/2)*exp(-Z*r1);
3
Z 1/2
(%o4) phi1(r1) := (---)
exp((- Z) r1)
%pi
(%i5)
phi2(r2):=(Z^3/%pi)^(1/2)*exp(-Z*r2);
3
Z 1/2
(%o5) phi2(r2) := (---)
exp((- Z) r2)
%pi
(%i6)
va(r1):=(1/r1)*integrate(4*%pi*r2^2*phi2(r2)^2,r2,0,r1);
1 2 2
(%o6) va(r1) := --
integrate(4 %pi r2 phi2
(r2), r2, 0, r1)
r1
(%i7) va(r1);
2 2 - 2 r1 Z
3 1 (2 r1 Z + 2 r1 Z + 1) %e
4 Z (---- -
----------------------------------)
3 3
4 Z 4 Z
(%o7)
------------------------------------------------
r1
(%i8)
(%i8)
vb(r2):=integrate(4*%pi*r2*phi2(r2)^2,r2,r1,inf);
2
(%o8) vb(r1) := integrate(4
%pi r2 phi2 (r2), r2, r1,
inf)
(%i9) vb(r1);
- 2 r1 Z
(%o9) Z (2 r1 Z + 1) %e
(%i10)
With the final integration
(SAGE code)
J=integrate(4*%pi*r1^2*(va(r1)+vb(r1))*phi1(r1)^2,r1,0,inf);
J = 5Z/8
Numerical integration of J(1s,1s) .
Suppose we want to find the coulombic repulsion of two electrons occupying orbitals p and q . An approximation to
J(p,q) can be found by taking first the spherical averages of φ2 p ( r1) and φ2 q ( r2) which is simply
φ2 p (ave) ( r) = (1/4π) ∫ φ2 p ( r) dΩ = (1/4π) ∫ ( R p ( r1) )2 Y* l,m Yl,m dΩ ( 15 ),
where the differential of the solid angle is dΩ = sin( θ ) dθ dφ' , ( φ' is the angle of r with the Z axis). The angular part integrates to 1.
The spherical average is the radial part of the wave function divided by the solid angle 4π ,
φ2 p (ave) ( r) = (1/4π) { R p ( r1) }2 . ( 16 )
The s (i.e. l =0) states are spherically symmetric. In this notes we take the spherical simplification for states with l ≠ 0.
A detailed treatment of the angular part is given in Ref. 2 , chapter IX .
A numerical integration of (9) is performed by the FORTRAN code given here.
Fortran code for J(1s,1s)
c spherical average for
J(p,q)
c let r1~ x and r2 ~ y
implicit real*8(a-h,o-z)
data z/2.d0/
data xi,yi,nstep/0.d0,0.d0
,4000/
equivalence (dx,dy),(xf,yf)
phi1s(r)=(z**3/pi)**.5d0*dexp(-z*r)
pi=2.d0*dasin(1.d0)
xf=12.d0/(2.d0*z)
dx=(xf-xi)/dfloat(nstep)
sum=0.d0
do 10 ix=1,nstep
x=xi+dx*dfloat(ix)
do 20 jy=1,nstep
y=yi+dy*dfloat(jy)
va=(1.d0/x)*4.d0*pi*y**2*phi1s(y)**2
vb=4.d0*pi*y*phi1s(y)**2
if(x.ge.y)sum=sum+dx*dy*(va*4.d0*pi*x**2*phi1s(x)**2)
if(x.lt.y)sum=sum+dx*dy*(vb*4.d0*pi*x**2*phi1s(x)**2)
20 continue
10 continue
print*,'sum=',sum
stop
end
J= sum= 1.24968122
Duplicating the parameter Z , to Z=4 in the numerical calculation gives J= 2.50 corroborating the analytical result that,
J(1s,2s) = (5/8) Z .
Perturbation theory
The approximation to the ground state of helium that was given above constitutes only the first term in a series of corrections.
Let the hamiltonian of helium be seprated in H = H0 + H1 , where
H0 = -(1/2)(p1)2 -(1/2)(p2)2 - Z /r1 - Z /r2 , ( 17 )
and H1 = 1/r12 . ( 18 )
We give a summary of time independent perturbation theory. Wave functions will be labeled with generic subscript n representing the associated quantum numbers.
Suppose the eigenfunctions of H0 are known. This means threre is a set of orthonormal wavefunction { Ψ n } , such that
< Ψ n / Ψ k >= δn,k , and
H0 Ψ n = E0 n Ψ n . ( 19 )
When the perturbation H1 is included the wave function and the energy have to be corrected .One writes for the n-th state
( H0 + H1 ) (Ψ n + ∑k≠n ck Ψ k ) = ( E0 n + E1 n + E2 n +...)( Ψ n + ∑k≠n ck Ψ k ) ( 20 )
Notice that the summation over k does not include the n wave function. The superscripts on the energy indicate , E1n first order correction , E2n is the second order correction and so on.
Equation ( 20 ) is separated by the different orders of correction . The wave function Ψ k are quantities to zero order, the
coefficients ck are assumed to be quantities of first order.
zero order H0 Ψ n = E0 n Ψ n (21 )
first order H0 ∑k≠n ck Ψ k + H1 Ψ n = E0 n ∑k≠n ck Ψ k + E1 n Ψ n ( 22 )
second order H1 ∑k≠n ck Ψ k = E1 n ∑k≠n ck Ψ k + E2n Ψ n (23 )
Using bra-ket notation , from (21) E0 n = < Ψ n / H0 / Ψ n > . Multiplying eq (22) by < Ψ n / gives the first order correction
E1 n = < Ψ n / H1 / Ψ n > . (24)
In the helium ground state example we have already worked , E1 n = J(1s,1s) , see equation (6) .
The coefficients ck are found using multiplying (22) by <Ψ p / , giving
< Ψ p /H0 / ∑k≠n ck Ψ k > + < Ψ p / H1 / Ψ n > = E0 n < Ψ p / ∑k≠n ck Ψ k > . Carrying out the multiplications and integrations
yields
E0p cp + < Ψ p / H1 / Ψ n > = E0 n cp or to first order (relabel p→ k ) ,
ck = < Ψ k / H1 / Ψ n > /( E0 n - E0p ) . (25)
Let us find the second order correction to the energy ,E2n , using (23) and the expression for ck . Multiplying (23) by <Ψ n / gives
E2n = ∑k≠n ck < Ψ n /H1 / Ψ k >
= ∑k≠n ( < Ψ n /H1 / Ψ k > )2 /(E0 n - E0k ) . (26)
If n represents the ground state, we conclude that E2n < 0 since the denominator in (26 ) , ( E0 n - E0k ) < 0.
Example: Calculate E 2 n for helium using (26) with only one term of (26). Estimate Eground state to second order.
We already know that E1 n = (5/8) Z .
We side step the important issue of the exclusion principle or the fact that electrons are indistinguishable . Let n=1 be a reference number representing the ground state, here it is not a quantum number . Let's label the ground state once more by
Ψ 1 = φ1s (r1) φ1s ( r2) = (Z3 /π) 1/2 exp(-Zr1) (Z3 /π) 1/2 exp(-Zr2) , (27)
as in eq .(1) .
For the the state Ψ 2 we postulate that the first electron is in state 1s and the second in state 2s and write
Ψ 2 = φ1s (r1) φ2s ( r2) = (Z3 /π) 1/2 exp(-Zr1) * {Z3/(32π ) }1/2 (2 - Zr2 )exp(-Z r2/2) .(28)
The terms in the denominator of (26) are E0 n = E0 1 = energy of 2 electrons in 1s hydrogenic state (without repulsion)
= 2( -Z2 /2) = -Z2 (29)
and E0k = E02 = energy of a hydrogenic atom with 1s electron plus energy of hydrogenic atom in 2s state
= -Z2 /2 - Z2 /8 . (30)
The integral in the numerator is
< Ψ n /H1 / Ψ k > = ∫ φ1s (r1) φ1s ( r2) (1/r12 ) φ1s (r1) φ2s ( r2) dτ1dτ2 . (31)
We adopt a new label for the integral (31) , wishing to be more explicit ,
J(1s,1s// 1s,2s) = ∫ φ1s (r1) φ1s ( r2) (1/r12 ) φ1s (r1) φ2s ( r2) dτ1dτ2 . (32)
This indicates that to the first electron are ascribed orbitals 1s ,1s and the second electron is represented by orbitals 1s,2s.
Adapting the results of (7)-(9) , one way the get the coulombic repulsion is given by
J(1s,1s// 1s,2s) = ∫ φ 1s (r1) φ 2s (r1) V(r1) dτ1 , (33)
with V(r1) = Va (r1) + Vb (r1) = { 1/r1 -( 2r1 Z2 + 2Z +1/r1) exp(-2r1 Z) }
+ { Z (2 r1 Z + 1) exp( - 2 r1 Z ) } ,
= 1/r1 - ( Z +1/r1 ) exp(-2Zr1) . (34)
as in (10) and (11).
Fortran code #2 is used to calculate (33) . The result is j(1s,1s//1s,2s) = 0.1840.
j1s2s,e1,e2,en1=
0.184006431 -4. -2.5 1.25
a,et,en2= 0.2000E+01
-0.2773E+01 -0.2257E-01
E0n = -4.000 , En1 = 1.25 , E 0 2 = -2.5 (one electron in 1s state and the second in 2s ) ,
En2 = J(1s,1s//1s,2s)2 /( -4 - (-2.5) ).
Etotal = -4.000 +
1.250 +(0.1840)2
/ [ -4.- (-2.5)]
(35)
= -2.77257061 , being a
small improvement from -2.75
, the first order
perturbation calculation.
Fortran code #2 for J(1s,1s//1s,2s)
c spherical average for
J(p,q)
c let r1~ x and r2 ~ y
implicit real*8(a-h,o-z)
real *8 J1s2s
data z/2.d0/
data xi,yi,nstep/0.d0,0.d0
,2000/
equivalence (dx,dy),(xf,yf)
phi1s(a,r)=(a**3/pi)**.5d0*dexp(-a*r)
phi2s(a,r)=(a**3/(32.d0*pi))**.5d0*(2.d0-a*r)*exp(-a*r/2.d0)
phi1s2s(a,r)=phi1s(a,r)*phi2s(a,r)
pi=2.d0*dasin(1.d0)
a=z
do 30 ia=1,1
xf=12.d0/(2.d0*a)
dx=(xf-xi)/dfloat(nstep)
sum=0.d0
do 10 ix=1,nstep
x=xi+dx*dfloat(ix)
do 20 jy=1,nstep
y=yi+dy*dfloat(jy)
va=(1.d0/x)*4.d0*pi*y**2*phi1s(a,y-dy/2.d0)**2
vb=4.d0*pi*y*phi1s(a,y-dy/2.d0)**2
if(x.ge.y)sum=sum+dx*dy*(va*4.d0*pi*x**2*phi1s2s(a,x-dx/2.d0))
if(x.lt.y)sum=sum+dx*dy*(vb*4.d0*pi*x**2*phi1s2s(a,x-dx/2.d0))
20 continue
10 continue
j1s2s=sum
e1=2.d0*(a**2/2.d0-z*a)
en1=5.d0*a/8.d0
e2=.25d0*(a**2/2.d0-z*a)
+(a**2/2.d0-z*a)
En2=sum**2/(e1-e2)
et=e1+en1+en2
print*,'j1s2s,e1,e2,en1=',j1s2s,e1,e2,en1
print*,' '
print 100,a,et,en2
a=a-.1d0
30 continue
100
format('a,et,en2=',3(3x,e11.4))
stop
end
Variational Calculation
Results of the perturbation theory (14) and (35) are markedly improved when a variational parameter α is introduced in the wave function , substituting for example the parameter Z such that Ψ (Z,r) → Ψ ( α , r) .
Suppose one introduces a normalized trial wave function Ψtrial for a given hamiltonian and attempts to estimate the ground state energy E0 from the intrgral ,
Etrial = ∫ ( Ψtrial )* H Ψtrial dτ . ( 36 )
We will show that Etrial ≥ E 0 . The proof goes as follows. Corresponding to the hamiltonian , there is a complete set of orthonormalized wave functions { Ψn } with energies E0 < E1 < E2 ....< Ei .
Ψtrial has the expansion , Ψtrial = ∑n=0 an Ψn and from the normalization condition , ∑n=0 an* an = 1.
From (36) we get , Etrial = ∑n=0 an* an En ≥ ∑n=0 an* an E0 = ( ∑n=0 an* an ) E0 = (1) E0 ;
thus Etrial ≥ E0 . (37)
The use of Ψ ( α , r) in (36) gives a functional Etrial ( α ). One finds the minima by taking the derivative d Etrial /dα and equating to zero yields α 0 . Then the minimum energy estimate is Etrial (α 0 ).
Combining the perturbation and variational methods give much improved values of the energy, as we will show.
Helium ground state energy by variation and perturbation method
Let
Ψ (α, r1 , r2) = φ1s (α , r1) φ1s ( α, r2) = (α 3 /π) 1/2 exp(- α r1) (α 3 /π) 1/2 exp(-α r2) . (38)
E0 n = < Ψ (α, r1 , r2) /H0 / Ψ (α, r1 , r2) >
=< Ψ (α, r1 , r2) / -(1/2)(p1)2 -(1/2)(p2)2 - Z /r1 - Z /r2 / Ψ (α, r1 , r2) > .(39)
The kinetic energy terms are
< φ1s (α , r1) / -(1/2)(p1)2 / φ1s (α , r1) > = < φ1s (α , r2) / -(1/2)(p2)2 / φ1s (α , r2) >
= a2 /2 . (40)
Sage code
var('a'); assume (a>0);
phi1s(r)=(a^3/pi)^(1/2)*exp(-a*r)
t1=
integral(phi1s(r)*(-(1/2)*diff(phi1s(r),r,2)
)*4*pi*r^2,r,0,oo);
t2=integral(phi1s(r)*(-diff(phi1s(r),r,1)
)*4*pi*r,r,0,oo);t1+t2
The potential energy terms are
< φ1s (α , r1) / -Z/r1 / φ1s (α , r1) > = < φ1s (α , r2) / -Z/r2 / φ1s (α , r2) >
= -Z a . (41)
Sage code
var('Z,a,r'); assume
(a>0);assume (Z>0);
phi1s(r)=(a^3/pi)^(1/2)*exp(-a*r);
V=integral(phi1s(r)^2*(-Z/r)*4*pi*r^2,r,0,oo);V
The repulsion term in the perturbation-variation calculations is
J(1s,1s//1s,1s) = < φ2 1s (α , r1) / 1/r12 / φ1s 2(α , r2) >
=(5/8) α , (42)
recall we had before (5/8) Z .
So for the ground state to first order (a = α )
En (a ) = 2 ( a2 /2 -Z a) + (5/8)a ,
= a2 + ( 5/8 -2Z ) a . (43)
d E /da= 2a + ( 5/8 -2Z ) = 0 , yields a0 = Z- 5/16 , and
Susbstitution in (43) gives Emin = E (a0) = (Z-(1/2)(5/8))2 - 2 (Z-(1/2)(5/8) )2
= - (Z - 5/16)2 . (43)
In the case of helium the estimate
E min = - 2.848 au , (44)
is a considerable improvement over the simple perturbation result given in (14-b).
Now we will consider a single term contributing to E2 n . We can use expression (33) and modify the FORTRAN code #2 by the
introduction of α instead of Z to calculate J(1s,1s// 1s,2s).
We find J with Z=1 to be 9.200E-2 and J when Z=2 to be 1.840E-1 ,thus we conclude that
J(1s,1s// 1s,2s) = ( 9.200E-2 ) α . hence J 2 = 8.464E-3 α2 . In our one term approximation ,
E2 n = 8.464E-3 a2 /{ (a2 /2-za) - (1/4)(a2 /2 -Za) } . (45)
Then
En =2 ( a2 /2 -Z a) + (5/8)a + 8.464E-3 a2 /{ (a2 /2-za) - (1/4)(a2 /2 -Za) } .(46)
The plot of E vs a shows a minima of - 2.86 au when a =1.69 .
Fig 2. Results of a variation perturbation calculation up to a second order term .
Energy value of 1s2s , 1S for the helium atom
We are still not mentioning the electron spin role in the energy levels. This will come shortly. But now let's examine an excited state. The wave function is to a first approximation
Ψ (r1 , r2) = φ1s (r1) φ2s ( r2) . (47)
Assume that Ψ (r1 , r2) contains the variation parameter a .
Then to first order
En (a) = (a2 /2 -Z a) +(1/4)(a2 /2 -Z a) + ∫ {φ2s ( r2)}2 (1/r12){φ1s (r1)}2 dτ1 dτ2 ,
=(5/4)(a2 /2 -Z a) + ∫ {φ2s ( r)}2 V(r) dτ , (48)
where V(r) = 1/r - ( a +1/r ) exp(-2ar) is given by eq. (12) .
Evaluating the integral in (48) gives (17/81) a . Hence
E n(a) = (5/4)(a2 /2 -Z a) + (17/81) a . (49)
a0 = {Z - (4/5)*(17/81) } For Z=2 , a0 = 1.83209872 and E(a0)= -2.0978663 ,the exact value is -2.15 au and the Hartree Fock value is -2.1698 , ( see ref . 8 , page 133)
Sage code to evaluate ∫ {φ2s ( r)}2 V(r) dτ = J(1s,2s) = (17/81) a
var('r,a'); assume (a>0);
V(r) = 1/r - ( a +1/r
)*exp(-2*a*r);
phi2s(r)=(a^3/(32*pi))^(1/2)*(2
- a*r)*exp(-a*r/2)
#phi2s(r)=exp(-a*r);
J= integral(phi2s(r)^2*V(r)*4*pi*r^2,r,0,oo);
The next correction to (49) will use two terms from the second order infinite series En2 .
We are reminded that the states are labeled as shown next.
state- n 1s2s Ψn = φ 1s (r1) φ 2s (r2) , En = (5/4)(a2 /2 -Za) (50-a)
state -1 1s1s Ψ1 = φ 1s (r1) φ 1s (r2) , E1= 2 (a2 /2 -Za) (50-b)
state 2 2s2s Ψ2 = φ 2s (r1) φ 2s (r2) , E2= 2(1/4) (a2 /2 -Za) (50-c)
If a=Z=2 ,
the energies would be
E1, En, E2= -4. -2.5 -1. in
au .
Our present En2 calculation will include one state below E1 and one state above E2 ,
En2 = < Ψn /1/r12 / Ψ1 > 2 /( E n - E1 ) + < Ψn /1/r12 / Ψ2 > 2 /( E n - E2 ) . (51)
The first term is positive since E n - E1 > 0 while the second is negative , since E n - E2 < 0.
We already have calculated the first term numerically that time. See eq (33) and the follow up.
Now we do it analytically ,
< Ψn /1/r12 / Ψ1 > = <φ 1s (r1) φ 2s (r2) /(1/r12 ) /φ 1s (r1) φ 1s (r2) >
= ∫ φ 1s (r) φ 2s (r) V(r) dτ
= {(4096 * 21/2 )/64827} a
= 0.0893550292 * a . (52)
To calculate < Ψn /1/r12 / Ψ2 > we define in analogy with (7)-(8)
V 2sa(r1) = (1/r1) ∫0 r1 4π (r2 )2 φ2 2s ( r2) dr2 r1 ≥ r2 , (53)
V2sb (r1) = ∫∞ r1 4π r2 φ2 2s ( r2) dr2 r1 ≤ r2 . (54)
Their values are
V2sa (r) = -(1/8)*((Z^4*r^4 + 4*Z^2*r^2 + 8*Z*r + 8)*e^(-Z*r)/Z^3 - 8/Z^3)*Z^3/r (55)
V2sb(r) = (1/8)*(Z^3*r^3 - Z^2*r^2 + 2*Z*r + 2)*Z*e^(-Z*r) . (56)
We have
<Ψn
/1/r12 /
Ψ2 > = ∫
φ 1s
( r) φ
2s
(r) V(r) 4π r2
dr
= {512*sqrt(2)/84375}
a = 0.00858165740960029*a
.
(59)
Sage code for V2sa(r) and
V2sb(r)
var('Z,r1,r2'); assume
(Z>0);assume(r1>0);
Sage code for
<Ψn
/1/r12 /
Ψ2 >
var('Z,r1,r2'); assume
(Z>0);assume(r1>0);
Two terms are added to (49)
E n(a)
= (5/4)(a2
/2 -Z a) +
(17/81) a + (8.936E-2
a)2 /(En
-E1)
+(8.582e-3*a)2
/(En -E2
) . (60)
sage code for the plot of of
En to third ordes
as given by (60)
var('a');Z=2;En=(5/4)*(a^2/2
-Z*a); E1= 2*(a^2/2
-Z*a);E2= 2*(1/4)*(a^2/2
-Z*a);
Fig 3 .
Energy , given by (80) vs the variation
parameter a .
At the value a=1.80 ,the
term (8.936E-2* a)^2
/(En -E1) = 1.74E-2
while the next term
(8.582e-3*a)^2/(En -E2)
equals -1.61E-4.
Energy value of 1s2p ,
1P for the
helium atom ( Eexact
= -2.1238 ref . 8 )
Using first order
variation-perturbation
theory and considering the
wave function to be a simple
product ,Ψ
=
φ1s (r1)
φ2p ( r2),
we have that
En0 =
< Ψ / H0 /Ψ
> = (5/4)(a2 /2
-Z a)
( 61 )
and the first order term is
En = < Ψ /
(1/r12 ) / Ψ
> = (59/243) a
( 62 )
We show below the code
employed for result
(62) .
En =
(5/4)(a2 /2 -Z a)
+ (59/243) a
The value that minimizes E
is a0 =
1.80576134 and
E( a0 ) =
-2.03798366 au representing
a 4 % error with respect to
the exact value.
The exclusion principle in
two electron atoms
An approximation for the total wave function for
an atom of N electrons and N
quantum state labels, is
given by the Slater's determinant
Fig 4. Slater'
determinant.
The most important
consequence of the
anti-symmetric wave function
is that the coulomb
repulsion , and the total
energy, are reduced by
negative energy terms
called exchange energy.
These are possible when two
electron in the atom have
parallel spins. This will be
illustrated below by
comparing helium excited
states energies for 1s2s
1S (electrons
have anti parallel spins)
and 1s2s 3S
(electrons have parallel
spins).
The label 1,2,3 N is a
shorthand for the state
label and all independent
variables including the spin
degree of freedom.
For example the first row
refers to electron number 1
. In φ1 (1)
occupies the quantum state
1, in φ1 (2 )
occupies the quantum state 2
and so on , until φ1
(N) represents the electron
in quantum state N.
Example write the wave
function for the ground
state of helium., i.e.
configuration 1s2
, 1S . (See fig
5)
One must now include spin in
an explicit way. The
notation φ1s
(1) α1 means that
electron number 1 occupies
the 1s state and its spin
function is α1 =
( 1,0)tranposed .
It is said that the spinor α1
, with a column of elements
1 and 0, represents the spin
"up: state.
The next element in the
first row, φ1s
(1) β1 ,
represents electron 1 in
state 1s with spin "down" .
The function β1 =
( 0,1)tranposed
is a column with elements 0
and 1.
We have also that φ1s
(2) α2 represents
electron number 2 in the 1s
state with spin up and φ1s
(2) β2 is
electron number 2 in the 1s
state with spin down.
Fig 5. Wave function of
helium ground state.
The complete wave
function is normalized
∫ {Ψ (1,2)}2
dτ1
dτ 2 =1
for ,
∫ {Ψ (1,2)}2
dτ1
dτ 2 = (1/2)∫ {
(φ1s
(1) α1 )2
(φ1s
(2) β2 )2
- 2 (φ1s
(1) α1 φ1s
(1) β1) ( φ1s
(2) β2 φ1s
(2) α2 )
+ ( φ1s
(1) β1 )2
(φ1s
(2) α2 )2
} dτ1 dτ
2
= (1/2) { 1
-2(0) +1 } = 1
.
(63)
For the ground state energy
, i.e. the configuration 1s2
, 1S , we have
the same results as before
E = < Ψ (1,2) / h01
+ h02 + 1/r12
/ Ψ (1,2) > =
< φ1s
(1)/ h01 /φ1s
(1) > + <φ1s
(2) / h02 / φ1s
(2)>
+ <φ1s
(1) 2/ 1/r12
/ φ1s
(2) 2 >
, (64)
where h01
= -(1/2) ∆1
-Z/r1
and
h02 = -(1/2) ∆2
-Z/r2
.The fact that <α1
/β1 > =
<α2 /β2
> = 0 is crucial in
the above results.
Fig 6. Wave function of state
1s2s 1S
(anti parallel spin
electrons)
For the state 1s2s 1S
the total wave function is
Ψ (1,2) = (1/21/2)
{ φ1s
(1) α1 φ2s
(2) β2 -
φ2s (1) β1
φ1s (2) α2
}
.
(65)
The energy is
E = <(1/21/2) {φ1s
(1) α1 φ2s
(2) β2 -
φ2s (1) β1
φ1s (2) α2
}/ h01 + h02
+ 1/r12 /
(1/21/2) { φ1s
(1) α1 φ2s
(2) β2
- φ2s
(1) β1 φ1s
(2) α2 } >
. (66)
Just as in the derivation of
(64) the orthogonality of
the spinors <α1
/β1 > = <α2
/β2 >
reduces (66) to
E = < φ1s
/ h0
/ φ1s
> +
< φ2s
/ h0
/ φ2s >
+ < φ
(1)1s 2
/ 1/r12
/φ (2) 2s 2
> .
(66)
A variation result for (66)
was given in eq. (49) .
Energy value of 1s2s ,
3S for the
helium atom
The Slater determinant is
obtained as in Fig . 7 .
Fig 7. Wave
function for state 1s2s
3S (
parallel spin electrons )
When comparing this state
with
1s2s , 1S
we see that there is an
extra negative term -K(1s,2s) contributing to
the total energy.
Ψ (1,2) = (1/21/2)
{ φ1s
(1) α1 φ2s
(2) α 2 -
φ2s (1) α1
φ1s (2) α2
}
.
(67)
Compare equation (67) with
(65) . In (67) all spinors
are of the α type , i.e. all
spins are up. This has an
important consequence for
the average < Ψ (1,2)
/ 1/r12 / Ψ (1,2)
> giving an extra negative
term that reduces the
coulombic repulsion.
The first term is the
coulombic repulsion
J(1s,2s) = ∫ φ
21s
(1) (1/r12) φ2 2s
(2) dτ1 dτ2
and the second is the
exchange term
K(1s,2s) = - ∫ φ 1s
(1) φ 2s
(1) (1/r12) φ 1s
(2) φ 2s
(2) dτ1 dτ2
.
(68)
we alreday have that
J(1s,2s)=(17/81) a
see the details after eq
(49). To find K(1s,2s)
we again use a MAXIMA code
integration.
The answer is K(1s,2s) = -
(16/729) z or
-(16/729) a
where a is the variation
parameter. Comparing with
J(1s,2s), we see that
K (1s,2s) ≈ - (1/9)
J(1s,2s) .
The total energy as a
function of one variation
parameter a is
E(a) = ( a2 /2 -Z
a) + (1/4) ( a2
/2 -Za) + (17/81) a
-(16/729) a
. (69)
Fig 8. Variational energy
for 1s2s ,
3S showing a
minima near E = -2.135 au
.The exact energy is -2.175
see Ref. 8 , page 173.
Maxima code for K(1s,2s)
%i15) Maxima 5.19.2 http://maxima.sourceforge.net
A short Table of coulombic
J(n,p) and exchange terms
K(n,p) is given when
hydrogenic wave functions
are used.
Table 1. Some terms for J
and K.
Coulombic terms
using hydrogenic wave
functions Exchange terms using
hydrogenic wave functions
Maxima code
(%i1) assume(z>0);
(%o1)
[z > 0] (%i2)
assume(r1>0);
(%o2)
[r1 > 0] (%i3)
assume(r2 >0);
(%o3)
[r2 > 0] (%i4)
phi1s:(z^3/%pi)^(1/2)*exp(-z*r1);
3/2
- r1 z
z %e
(%o4)
-------------
sqrt(%pi) (%i5)
phi2p:(z^3/(96*%pi))^(1/2)*z*r2*exp(-z*r2/2);
r2 z
- ----
5/2 2
r2 z %e
(%o5)
-------------------
4 sqrt(6) sqrt(%pi) (%i6)
va:(1/r1)*integrate(4*%pi*r2^2*phi2p^2,r2,0,r1);
Improper value assignment: 'va -- an
error. To debug this try
debugmode(true); (%i7)
va:(1/r1)*integrate(4*%pi*r2^2*phi2p^2,r2,0,r1);
4
4 3 3 2
2 - r1 z
5 24 (r1 z + 4 r1 z +
12 r1 z + 24 r1 z + 24) %e
z (-- -
-------------------------------------------------------)
5
5
z
z (%o7)
-----------------------------------------------------------------
24 r1 (%i8)
vb:integrate(4*%pi*r2*phi2p^2,r2,r1,inf);
3
3 2 2
- r1 z
z (r1 z
+ 3 r1 z + 6 r1 z + 6) %e
(%o8)
-------------------------------------------
24 (%i9)
j1s2p:integrate(phi1s^2*(va+vb)*4*%pi*r1^2,r1,0,inf);
59 z
(%o9)
----
243 (%i10)
Integration of
J(1s,1s) by Monte Carlo
method
answ, int(f(u))
= 0.625 ,
0.620409131
c Montecarlo integration of
J(1s,1s) -a < x < a
phi1s(r)=(Z^3/pi)^(1/2)*exp(-Z*r);
phi2s(r)=(Z^3/(32*pi))^(1/2)*(2
- Z*r )*exp(-Z*r/2);
#va=(1/r1)*integral(phi2s(r)^2*4*pi*r^2,r,0,r1)
#va
vb=integral(phi2s(r)^2*4*pi*r,r,r1,oo)
vb
phi1s(r)=(Z^3/pi)^(1/2)*exp(-Z*r);
phi2s(r)=(Z^3/(32*pi))^(1/2)*(2
- Z*r )*exp(-Z*r/2);
va=-(1/8)*((Z^4*r^4 +
4*Z^2*r^2 + 8*Z*r +
8)*e^(-Z*r)/Z^3 -
8/Z^3)*Z^3/r ;
vb= (1/8)*(Z^3*r^3 - Z^2*r^2
+ 2*Z*r + 2)*Z*e^(-Z*r);v=va+vb;
J=integral(phi1s(r)*phi2s(r)*v*4*pi*r^2,r,0,oo)
J
e=(5/4)*(a^2 /2 -Z*a) +
(17/81)*a + (8.936E-2* a)^2
/(En -E1)+
(8.582e-3*a)^2/(En -E2);
plot(e,a,1.4,2.2)
Using Lisp GNU Common Lisp (GCL)
GCL 2.6.8 (aka GCL)
Distributed under the GNU
Public License. See the file
COPYING.
Dedicated to the memory of
William Schelter.
The function bug_report()
provides bug reporting
information.
(%i1)assume(z>0);assume(r>0);assume(r1>0);
(%o1) [z > 0]
(%o2) [r > 0]
(%o3) [r1 > 0]
(%i4)
phi1s(u):=(z^3/%pi)^(1/2)*exp(-z*u);
3
z 1/2
(%o4) phi1s(u) := (---)
exp((- z) u)
%pi
(%i5)
phi2s(u):=(z^3/(32*%pi))^(1/2)*(2
- z*u )*exp(-z*u/2);
3
z 1/2 (- z) u
(%o5) phi2s(u) := (------)
(2 - z u) exp(-------)
32 %pi 2
(%i6)
va(r1):=(1/r1)*integrate(4*%pi*r^2*phi1s(r)*phi2s(r),r,0,r1);
1 2
(%o6) va(r1) := --
integrate(4 %pi r phi1s(r)
phi2s(r), r, 0, r1)
r1
(%i7)
vb(r1):=integrate(4*%pi*r*phi1s(r)*phi2s(r),r,r1,inf);
(%o7) vb(r1) := integrate(4
%pi r phi1s(r) phi2s(r), r,
r1, inf)
(%i8)
integrate(phi1s(r1)*phi2s(r1)*(va(r1)+vb(r1))*4*%pi*r1^2,r1,0,inf);
16 z
(%o8) ----
729
J(n,p) (n=row) ,(p=column)
1s
2s
2p
1s
(5/8)*Z
(17/81)*Z
(59/243)*Z
2s
(17/81)*Z
(77/512)Z
(83/512) Z
2p
(59/243)*Z
(83/512) Z
(93/512)Z
K(n,p) (n=row) ,(p=column)
1s
2s
2p
1s
-(16/729)Z
-(112/2187) Z
2s
-(16/729)Z
-(45/512) Z
c fix x1,x2, vol=scale**(dimensions)
, fix step 50 for number of
random
c arrays
dimension
xi(6),x1a(10000),y1a(10000),z1a(10000)
data z,nstep,nstepr/1.,5000,8000/
data seed,ncoord/2.5,3/
psi1s(r)=sqrt(z**3/pi)*exp(-z*r)
f(r1,r2)=psi1s(r1)**2*psi1s(r2)**2*(1./r12)
c f(u)=exp(-r**2/2.)
c f(u)=u
c f(u)=sin(u)**2
pi=2.*asin(1.)
x1=-3.2
x2= -x1
scale=(x2-x1)
vol=scale**6
slope=x2-x1
b=x1
kount=0
c creates one,two or three
arrays of rnd num x1a
,y1a,z1a
do 50 i=1,ncoord
x=sqrt(seed)*pi**(float(i+1))
do 55 j=1,nstepr
kount=kount+1
if(i.eq.1)then
pwr=3./5.
call
aleat(pwr,x,kount,slope,b,pi,zr)
x1a(j)=zr
x=abs(zr)*(float(j)+2.)+2.*pi
c if(x.eq.0.)x=float(kount)**.25
endif
if(i.eq.2)then
pwr=1./2.
call
aleat(pwr,x,kount,slope,b,pi,zr)
y1a(j)=zr
x=abs(zr)*(float(j)+2.5)+3.*pi
c if(x.eq.0.)x=float(kount)**.25
endif
if(i.eq.3)then
pwr=1./3.
call
aleat(pwr,x,kount,slope,b,pi,zr)
z1a(j)=zr
x=abs(zr)*(float(j)+2.5)+3.*pi
c if(x.eq.0.)x=float(kount)**.25
endif
55 continue
50 continue
sum=0.
do 10 i=1,nstep
xi(1)=x1a(i)
xi(4)=x1a(i+1)
xi(2)=y1a(i)
xi(5)=y1a(i+2)
xi(3)=z1a(i)
xi(6)=z1a(i+3)
c
print*,'x1,x2,x3=',xi(1),xi(2),xi(3)
r1=sqrt(xi(1)**2 + xi(2)**2
+ xi(3)**2)
r2=sqrt(xi(4)**2 + xi(5)**2
+ xi(6)**2)
r12=sqrt((xi(4)-xi(1))**2
+(xi(5)-xi(2))**2
+(xi(6)-xi(3))**2 )
c r=xi(1)
sum=sum + f(r1,r2)
10 continue
aint=vol*sum/float(nstep)
print*,'answ,int(f(u))=',
5.*Z/8. , aint
stop
end
subroutine
aleat(pwr,x,kount,slope,b,pi,zr)
y=x**pwr
rnd=y-int(y)
zr=slope*rnd + b
return
end
END OF CHAPTER