Create Presentation
Download Presentation

Download Presentation

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

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

655 Views

Download Presentation
Download Presentation
## Solving the radial Schr ö dinger equation of hydrogen atom for l = 0

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -

**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**Time dependent Schrodinger Equation**• With separation of variable • The time-dependent part is decoupled, resulting in time-independent Schrodinger equation**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**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]**Physical boundary conditions necessary for the s-state**hydrogen atom • R(r) 0 as r0 • R (r) 0 as r∞ Or equivalently • u(r) 0 as r0 • u (r) 0 as r∞ • From theory we also expect that as r∞, u(r) rexp (-r)**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**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**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**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**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**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**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?**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**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**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**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**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.**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”