Report

Time independent Schrödinger equation(TISE) in spherical form The first thing to do is to rewrite the Schrödinger equation in the spherical coordinate system. Solving TISE by separation of variable In our project we separate the code in to 4 parts 1st part finding the energy 2nd part solving the angular and radial function 3rd part get the probability density 4th part visualization of the probability density Part 1 Finding the energy Solve the Radial Part Energy spectrum is the important property in quantum mechanics, it can be obtain only by solving the radial part of Schrödinger equation. The solution for radial function will be obtain only if the corresponding energy is the Eigen energy for the radial function. Hence, in this part, we will keep guessing the energy until the solution obtain fit the boundary condition. 1st we have to substitute all constant in Radial function in to atomics unit to simplify our simulation. 2nd simplify our equation and solving it using the command DSolve or NDSolve . In my code, I choose the command DSolve for some reason To solve the equation using DSolve or NDSolve boundary condition is important: B.C. When r 0, R[r] 0; When r ∞, R[r] ∞; From theory, r ∞ , R[r] rexp (-r) ; In our code, we change the R[r] term to U[r], where U[r] = r*R[r] to avoid some singularity when the DSolve executed. 3rd we plug in a series of guess energy in radial function. Eguess[n]=Eguessinitial + n*epsilon. Epsilon is a very small interval which can be adjust to change the accuracy. Plot the graph U[0] versus n, it is equivalent to U[0] versus Eguess, when the graph intercept with the xaxis then corresponding energy is the Eigen energy we want. The intersection point can be found using bisection method. U[0] versus n, when l=0 U[0] versus n, when l=1 Part 2 Solving the angular and radial function 1 2 1 1 2 2m (r ) 2 (sin ) 2 2 ( E V ) 0 r 2 r r r sin r sin 2 2 For Hydrogen atom, the potential is given as: - V e2 4 0 r and therefore, the Schrödinger equation for Hydrogen atom is given by, 1 2 1 1 2 2m e2 (r ) 2 (sin ) 2 2 (E ) 0 r 2 r r r sin r sin 2 2 4 0 r The wavefunction (r, , ) can be dissolved into three components by assuming all of these three functions are independent of each another. By using separation of variables, (r , , ) R(r )( )( ) For convenience, we just write ( ) ( ) R R(r ) Finally, we will have for : (i) Equation for (ii) Equation for (ii) Equation for R d 2 2 m 0 : l d 2 0 ml l : 1 d (sin d ) [l (l 1) ml2 ] 0 sin d d sin 2 0 l (n 1) 2ml 1 d 2 dR e2 l (l 1) ( r ) [ ( E ) ]R 0 : 2 2 2 r dr dr 4 0 r r by universal scaling, let m e2 1 and 2l 1 4 0 0 l (n 1) Orthogonality Functions y1 ( x) y,2 ( x) ,……. defined on some interval a x b are called orthogonal on this interval with respect to the weight function r ( x) 0 if for all n and all m different from n . b a r ( x) ym ( x) yn ( x)dx 0 The norm ym ym m n of ym is defined by b a 2 r ( x) ym ( x)dx mn The functions y1 , y2 , …… are called orthonormal on a x b if they are orthogonal on this interval and all have norm 1. Why is the orthogonality property so important? It is to normalise the wavefunction. Why do we need to normalise the wavefunction? It is because we need the total area under the curve to be equal to 1. To do normalisation, for b a r mn , say r is independent of x (r ) ym ( x)dx 1 2 b a ym ( x)dx 1 2 ⇒ r 1 b a But, r A2 A r 1 b a 2 ym ( x)dx 2 ym ( x)dx Solving for ( ) :d 2 2 m l 0 d 2 (i) Periodic BCs : (0) (2 ) (ii) Periodic BCs : ' (0) ' (2 ) BOUNDARY CONDITION However, MATHEMATICA NDSolve function can only sustained a maximum number of 2 constraints for the 2nd order ODE. Indeed, it should be A , we set it to be 1 because NDSolve can take numerics only. Let (0) 1; (2 ) 1 Then, we have to normalise the function 2 0 2 d 1 NORMALISING CONDITION Solving for ( ) :First, we have to make use of substitution to change the function to (x) x cos ml d 2 d (1 x ) 2 2 x [l (l 1) ] 0 2 dx dx (1 x ) The boundary conditions are:2 2 l0 all ml (1) (1) 0 BOUNDARY CONDITION l0 ml is even ml is odd (1) (1) 1 BOUNDARY CONDITION W (1) W (1) 1 BOUNDARY CONDITION W ( x) x ( x) 2 2 ml 2 dW 2 2 d W (1 x ) 2 [ 2 l (l 1) ]W 0 dx x dx x (1 x 2 ) This is because theoretically, when l 0 , ml is even : x 1, (1) 1 ;x 1, (1) 1 l 0 , ml is odd : x 1, (1) 1 ; x 1, (1) 1 BUT x 1, x (1) 1 ; x 1, x (1) 1 Solving for R(r ) :There are 2 different cases for the radial function: R(r ) R(r ) l0 l0 r u (r ) r u ( r ) r R( r ) d 2u 2 l (l 1) 2 u u 2 En u 2 dr r r Because with this substitution, at r 0, d 2 R 2 dR 2 l (l 1) 2 R R 2 En R 2 dr r dr r r 0 R(0) 0 The boundary conditions are then set to be:- r 0 , u0 (ii) r , u re r (i) (i) r 0 ,R0 (ii) r , R re r For l 0 , when we transform the substitution function , we will have singularity at r 0 . u (r ) back into R(r ) To solve this, we have no other choices but to do some improper ‘shifting’, that is, When r 0 , (r ) R Since our plot range is u (r ) r r 0,3 , then we will set u (r ) r to be 0.1 R(r ) l0 l0 r With this modification, the singularity will be “hidden away”. r The same thing again, we have to normalise both the theta angular function and radial function as well. 0 0 2 dx 1 2 r R dr 1 2 NORMALISING CONDITION NORMALISING CONDITION Finally, we obtain the total wavefunction, (r, , ) to be the multiplication of all these 3 interpolating functions. Part 3 Obtaining the probability density Now, we have to plot the density plot of the wavefunction. It is given as (r , , ) 2 Now, the question arises. We are not able to visualise a function of R 3 . We have to choose 2 parameters to plot. Which 2 parameters are to be plotted? ( ) depends on ml ( ) depends on ml R(r ) depends on ml , l , En ,l In this situation, the energy level is dependent on different excited state, n , and for different excited energy state, different orbital, l will have different energy as well. Thus, we will plot the (r , , 0) 2 by varying only and r at a fixed 0 . First, we need to change the wavefunction from spherical coordinate back to the Cartesian coordinate, and plot all the points discretely. X (r, , ) r sin cos Y (r , , ) r sin sin Z (r , ) r cos But, we have substitution of x cos , and fixed r [0, rm ax] [0, ] [0,2 ] 0 , then sin 0 Thus, our plot is independent on Y (r , x, ) 0 However, the range of from 0 to is not complete. Therefore, we need to have another mirror image into the range of another [ ,2 ] . [2 , ] X (r, x, ) r sin(2 cos1 x) cos Z (r , x) r cos2 x [0, ] X (r, x, ) r sin(cos1 x) cos Z ( r , x) r x changing of plot interval, changing of plot interval, [0, ] [0, ] x [1,1] x [1,1] Part 4 visualising the probability density plot Refer to source code