Solving the radial Schrцdinger equation of hydrogen atom for l = 0

1
Solving the radial Schrödinger equation of
hydrogen atom for l 0

• Yoon Tiem Leong
• Presentation at the
• Weekly coffee meeting,
• The Theory Group
• School of Physics, USM
• 21 Nov 2008

2
Time dependent Schrodinger Equation

• With separation of variable

• The time-dependent part is decoupled, resulting
  in time-independent Schrodinger equation

3
Time independent Schrodinger Equation in
spherical coordinate

• The spatial function eigen function and the
  energy eigen value E are determined by solving an
  eigen value problem based on the time-independent
  Schrodinger equation

4
Solving TISE with separation of variables
5
Radial part simplified
6
Atomic units
7
Radial part simplified
8
The second order differential equation to solve
numerically

• In Mathematica, the differential equation is
  expressed as
• (-r/2)ur ur-ure r ur

9
Physical boundary conditions necessary for the
s-state hydrogen atom

• R(r) ?0 as r?0
• R (r) ? 0 as r?8
• Or equivalently
• u(r) ?0 as r?0
• u (r) ? 0 as r?8
• From theory we also expect that as
• r?8, u(r) ? rexp (-r)

10
Numerical boundary conditions

• To solve a second order differential equation
  numerically, two boundary conditions are
  necessary
• Since from theory we expect as r?8, u(r) ? r exp
  r
• Hence, in the numerical program, we set
  u(rmax)rmaxExprmax
• For the second B.C
• u(rmax-epsilon)(rmax-epsilon)Exprmax-epsilon
• rmax has to be chosen (with some trial and error)
  such that it simulates a cut-off such that u
  effectively drops to zero when r gt rmax

11
Eigen value of E

• The radial TISE for hydrogen is actually an eigen
  value problem, with discrete eigen value E that
  has to be solved for
• The numerical behavior of solution to u(r) is
  dependent on the value of E
• If E takes on the true value, u(r) will behave
  properly, i.e. u(r 0)0
• Else the boundary condition u(r0) 0 will not
  be satisfied
• These behavior will show up in the Mathematica
  solution

12
NDSolve command line

• soln_NDSolve
• -1u'r-r/2u''r-urr enur,
• urmaxrmaxExp-rmax,
• urmax-epsilon(rmax-epsilon)Exp-(rmax-epsilon
  ),
• u,
• r,epsilon/5,rmax
• Note that the range cannot be starting from r 0
  but only at r epsilon, or else the numerical
  code will not work due to computational artifact
  effect

13
Input values

• eneguess nepsilon
• Input values
• nlow-410500
• nup-390000
• eguess-0.1
• Epsilon N10(-6)
• tolN10(-5) ? determine the accuracy of the
  calculation

14
NDSolve with a trial E en

• The radial equation
• (-r/2)ur ur-ure r ur
• is solved with boundary conditions
• u(rmax)rmaxExprmax
• u(rmax-epsilon)(rmax-epsilon)Exprmax-epsilon
• with a trial value of eneguess nepsilon
• The energy E en, is parametrised in terms of
  n
• The result is a list of numerical values of ur
  as a function of r from repsilon to rrmax

15
The zero of u0

• Once the numerical solution of ur is obtained
  we can the check the value of u0 correspond to
  the value of en eguess nepsilon
• u0 is energy-dependent (controlled by n)
• We then plot a graph of u0 versus n to locate
  the interval of n within which the zero of u0
  occur
• To investigate the zero of u0 we have to tell
  the program the range of n, nlow, nup.
• Have to choose nlow, nup wisely

16
Root finding

• It should occur around n400,000, corresponds to
  eneguessepsilonn
• What is the exact value n with u0 zero?

17
Bisection method to find root

• Set two end values, n1, n2, such that
  un1,r0un2,r0 lt0, so that the root lies
  between n1, n2
• n1 -410500 usolzeron1 1.823106
• n2 -390000 usolzeron2 -1.30187106
• Then define nave(1/2)(n1n2), and evaluate
  unave,r0 to determine whether the sign of
  unave,r0un1,r0 or unave,r0un2,r0 is
  negative

18
Bisection method to find root (cont)
nave
n2
n1

• If unave,r0un1,r0gt0, the root must lies
  between (nave,n2), then set n1? nave
• If unave,r0un2,r0gt0, the root must lies
  between (n1,nave), then set n2? nave
• n1,n2 will be updated in every step
• After n1 or n2 has been updated, then update
  nave?(1/2)(n1 n2)
• The interval n1,n2 becomes narrower and
  narrower
• Stop until the criteria of either
• en1-en2lttol or unave,0 lt tol is met

19
Result

• Drop out from the loop once the tol criteria is
  met
• The most updated nave is the value of n of which
  u0 is zero
• In the example, nave -400000
• enave -0.5, c.f. theoretical expectation,
• E -0.5

20
Profile of u(r) vs r
21
Whats next

• To generalise to non-zero l for hydrogen atom
• To generalise the program to treat Helium atom as
  perturbation on the hydrogen atom, by including
  additional effects (apart from the Coulombic
  potential from nucleus) coming from corrections
  due to coulombic interaction between the
  electrons (Hatree pontential), exchange and
  correlation effect
• In particular, the Hartree interaction for He
  atom, due to the electrostatic potential
  generated by the charge distribution of one of
  the two electrons on the other one, can be
  calculated as

22
Hatree-Fock scheme

• The inclusion of Hatree interaction and
  exchange-correlation effect in He calculation has
  to be implemented in a self-consistent manner.
• The full program is called Hatree-Fock
  calculation, requiring extensive programming
  scheme that iterates the eigen energies of the
  multi-electron atom until the eigen-energies
  become convergent.

23
Conclusion

• The program developed in this talk calculates the
  eigen energy and the radial wave function of
  s-state hydrogen atom, and is readily expanded to
  treat cases of higher l
• In addition, the preliminary program presented
  here is the starting point to enter the full
  Hatree-Fock calculation for atoms with higher
  number of electrons

The pdf version of the Mathematica code of this
presentation can be found at the file
schrodinger4.pdf