Transcript Project Presentation

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
mn
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
mn
, 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
l0
all ml
(1)  (1)  0
BOUNDARY CONDITION
l0
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 )
l0
l0
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 , u0
(ii) r   , u  re  r
(i)
(i)
r 0
,R0
(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 )
l0
l0
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  cos1 x)  cos
Z (r , x)  r  cos2  x
  [0,  ]
X (r, x,  )  r  sin(cos1 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