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)u’’[r] – u’[r]-u[r]=e r u[r]

9. Physical boundary conditions necessary for the s-state hydrogen atom • R(r) 0 as r0 • R (r)  0 as r∞ Or equivalently • u(r) 0 as r0 • u (r)  0 as r∞ • From theory we also expect that as r∞, 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∞, u(r)  r exp r • Hence, in the numerical program, we set u(rmax)=rmax*Exp[rmax] • For the second B.C: u(rmax-epsilon)=(rmax-epsilon)*Exp[rmax-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 > 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(r=0) = 0 will not be satisfied • These behavior will show up in the Mathematica solution

12. NDSolve command line • sol[n_]:=NDSolve[ • {-1u'[r]-r/2u''[r]-u[r]==r e[n]u[r], • u[rmax]==rmax*Exp[-rmax], • u[rmax-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 • e[n]=eguess + n*epsilon • Input values: • nlow=-410500; • nup=-390000; • eguess=-0.1 • Epsilon= N[10^(-6)] • tol=N[10^(-5)]  determine the accuracy of the calculation

14. NDSolve with a trial E = e[n] • The radial equation • (-r/2)u’’[r] – u’[r]-u[r]=e r u[r] • is solved with boundary conditions • u(rmax)=rmax*Exp[rmax] • u(rmax-epsilon)=(rmax-epsilon)*Exp[rmax-epsilon] • with a trial value of e[n]=eguess + n*epsilon • The energy E = e[n], is parametrised in terms of n • The result is a list of numerical values of u[r] as a function of r from r=epsilon to r=rmax

15. The zero of u[0] • Once the numerical solution of u[r] is obtained we can the check the value of u[0] correspond to the value of e[n] = eguess + n*epsilon • u[0] is energy-dependent (controlled by n) • We then plot a graph of u[0] versus n to locate the interval of n within which the zero of u[0] occur • To investigate the zero of u[0] 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 n~400,000, corresponds to e[n]=eguess+epsilon*n • What is the exact value n with u[0] zero?

17. Bisection method to find root • Set two end values, n1, n2, such that u[n1,r=0]*u[n2,r=0] <0, so that the root lies between n1, n2 • n1= -410500 usolzero[n1]= 1.823×10^6 • n2= -390000 usolzero[n2]= -1.30187×10^6 • Then define nave=(1/2)(n1+n2), and evaluate u[nave,r=0] to determine whether the sign of u[nave,r=0]*u[n1,r=0] or u[nave,r=0]*u[n2,r=0] is negative

18. Bisection method to find root (cont) • If u[nave,r=0]*u[n1,r=0]>0, the root must lies between (nave,n2), then set n1 nave • If u[nave,r=0]*u[n2,r=0]>0, 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 • e[n1]-e[n2]<tol or u[nave,0] < tol is met nave n2 n1

19. Result • Drop out from the loop once the tol criteria is met • The most updated nave is the value of n of which u[0] is zero • In the example, nave= -400000 • e[nave]= -0.5, c.f. theoretical expectation, E = -0.5

20. Profile of u(r) vs r

21. What’s 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”