Solving radial SE numerically

Leave a comment

The time-independent Schrödinger equation is

(-\frac{\hbar^2}{2m}\nabla^2 + V ) \Psi = E \Psi

Using the Laplacian in spherical coordinate. and Set \Psi = R Y

\nabla^2 R Y - \frac{2m}{\hbar^2}(V-E) R Y = 0

\nabla^2 = \frac{1}{r^2}\frac{d}{dr}(r^2 \frac{d}{dr}) - \frac{1}{r^2} L^2

The angular part,

L^2 Y = l(l+1) Y

The radial part,

\frac{d}{dr}(r^2\frac{dR}{dr}) - l(l+1)R - \frac{2mr^2}{\hbar^2}(V-E) R = 0

To simplify the first term,

R = \frac{u}{r}

\frac{d}{dr}(r^2 \frac{dR}{dr})= r \frac{d^2u}{dr^2}

A more easy form of the radial function is,

\frac{d^2u}{dr^2} + \frac{l(l+1)}{r^2} u - \frac{2m}{\hbar^2} (V-E) u = 0

The effective potential U

U = V + \frac{\hbar^2}{m} \frac{l(l+1)}{r^2}

\frac{d^2u}{dr^2} + \frac{2m}{\hbar^2} (E - U) u = 0

We can use Rungu-Kutta method to numerically solve the equation.

RK4.PNG

The initial condition of u has to be 0. (home work)

I used excel to calculate a scattered state of L = 0 of energy 30 MeV. The potential is a Wood-Saxon of depth 50 MeV, radius 3.5 fm, diffusiveness 0.8 fm.

e2.PNG

Another example if bound state of L = 0. I have to search for the energy, so that the wavefunction is flat at large distance. The outermost eigen energy is -7.27 MeV. From the radial function, we know it is a 2s orbit.

ex1.PNG

Electromagnetic multi-pole moment

5 Comments

Electromagnetic multipole comes from the charge and current distribution of the nucleons.

Magnetic multipole in nucleus has 2 origins, one is the spin of the nucleons, another is the relative orbital motion of the nucleons.  the magnetic charge or monopoles either not exist or very small. the next one is the magnetic dipole, which cause by the current loop of protons.

Electric multipole is solely by the proton charge.

From electromagnetism, we knew that the multipole has  different radial properties, from the potential of the fields:

\Psi(r) = \frac{1}{4\pi\epsilon_0} \int\frac{\rho(r')}{|r-r'|}d^3r'

A(r) = \frac{\mu_0}{4\pi}\int\frac{J(r')}{|r-r'|}d^3r'

and expand them into spherical harmonic by using:

\frac{1}{|r-r'|} = 4\pi\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l} \frac{1}{2l+1}\frac{r_{<}^l}{r_>^{l+1}} Y_{lm}^*(\theta',\phi')Y_{lm}(\theta,\phi)

we have

\Psi(r) = \frac{1}{\epsilon_0} \sum_{l,m}\frac{1}{2l+1}\int Y_{lm}^*(\theta',\phi') r'^l\rho(r')d^3r' \frac{Y_{lm}(\theta,\phi)}{r^{l+1}}

A(r)=\mu_0 \sum_{l,m}\frac{1}{2l+1}\int Y_{lm}^*(\theta',\phi') r'^l J(r') d^3r' \frac{Y_{lm}(\theta,\phi)}{r^{l+1}}

we can see the integral give us the required multipole moment. the magnetic and electric are just different by the charge density and the current density. we summarize in this way :

q_{lm} = \int Y^*_{lm}(\theta',\phi') r'^l \O(r') d^3 r'

where O can be either charge or current density. The l determine the order of multipole. and the potential will be simplified :

M(r)=\sum_{l,m}\frac{1}{2l+1} q_{lm} \frac{Y_{lm}(\theta,\phi)}{r^{l+1}}

were M can be either electric or magnetic potential, and i dropped the constant. since the field is given by 1st derivative, thus we have:

  1. monopole has 1/r^2 dependence
  2. dipole has 1/r^3
  3. quadrapole has 1/r^4
  4. and so on

The above radial dependences are same for electric or magnetic. for easy name of the multipole, we use L-pole, which L can be 0 for monopole, 1 for dipole, 2 for quadrapole, etc.. and we use E0 for electric monopole, M0 for magnetic monopole.

Since the nucleus must preserver parity, and the parity for electric and magnetic moment are diffident.the different come from the charge density and current density has different parity. The parity for charge density is even, but for the current density is odd. and 1/r^2 has even parity, 1/r^3 has odd parity. therefore

  • electric L-pole — (-1)^{L}
  • magnetic L-pole — (-1)^{L+1}

for easy compare:

  • E0, E2, E4… and M1,M3, M5 … are even
  • E1,E3,E5…. and M0, M2, M4…. are odd

The expectation value for L-pole, we have to calculate :

\int \psi^* Q_{lm} \psi dx

where Q_{lm} is multipole operator ( which is NOT q_{lm}), and its parity is follow the same rule. the parity of the wave function will be canceled out due to the square of itself. thus, only even parity are non-Zero. those are:

  • E0, E2, E4…
  • M1,M3, M5 …

that make sense, think about a proton orbits in a circular loop, which is the case for E1, in time-average, the dipole momentum should be zero.

Deuteron

Leave a comment

The deuteron is the nucleus that contains 1 proton and 1 neutron. The spin and isospin of proton and neutron are the same, both are equal to half.  It is the only stable state for 2 nucleons. Deuteron provides an unique aspect to study the inter nuclear force. The strong force are believed to be charge independent. Thus, the strong force can be more easily to study on deuteron due to the absent of other force or eliminate from the Coulomb force, which is understood very much.

The mass of deuteron is 1876.1244MeV. The binding energy is then 2.2245MeV. It was determined by the slow neutron capture of a proton. The emitted gamma ray is approximately equal to the binding energy and the deuteron mass was calculated.

Deuteron has no excited state. It is because any excitation will easily to make the system break apart.

The parity is positive from experiment. If we separate the deuteron wavefunction into 3 parts. The proton wavefunction, neutron wavefunction and the orbital wavefunction. Under the only force, the strong force in this system, proton and neutron are the same nucleon with different state. Thus, the parity are the same for proton and neutron. So, the product of these 2 wavefunction always has positive parity. The total parity then is solely given by the angular orbital.

Any orbital wave function can be represented by the spherical harmonic, Y(l,m) .

The parity transform is changing it to

Y(l,m) \rightarrow (-1)^l Y(l,m)

So, the experimental face of positive parity fixed the angular momentum must be even.

Ok, we just predicted the possible angular momentum from parity.

The experimental fact on spin is 1. Since J = L + S, and the value of J can take every integer from |L-S| to L + S. and L must be even. The spin of proton and neutron is 1/2. Thus the possible S is 0 or 1 ( we are using L-S coupling scheme ). J = 1 = L + S , that tell us S must be odd to give out 1 for an even L. Thus S=1. So, the only possible L is 0 and 1. Thus, the possible state of deuteron is (L,S) = (0,1) or (2,1). Therefore, a deuteron is a mixed state, if without any further argument.

Now, 2 out of 3 parts of the wave function symmetry were determined by symmetry argument. The isospin can now be fixed by the 2 fermions state must be antisymmetry. The spatial state symmetry is even by L = 0 or 2. And for the state (L , S) = ( 0, 1 ), the spin state is symmetric. Thus, the isospin must be antisymmetric. Since the algebra for isospin and spin are the same. We use T = 0 for the isospin. Thus a complete wavefunction is ( L , S , T ) = ( 0 , 1, 0 ). For the other possible state (L , S) = ( 2 , 1 ) , we can use same argument for isospin state. And for the degenerated state with Ms = +1, 0, -1. By the symmetry of the raising and lowering ladder operator, they all preserved the symmetry. Thus, the Ms = 0 state can only be the + state.

So, we now have 2 possible states of deuteron. If the hamiltonian is commute with L^2 and  S^2, both L and L is a good quantum number and those states are eigen state. And the deuteron ground state must be one of them.

Laplacian in spherical coordinate

Leave a comment

the Momentum operator in spherical coordinate

\nabla^2 = \frac {1}{r^2}\frac {\partial } { \partial r} \left ( r^2 \frac {\partial} {\partial r} \right ) - \frac {1}{r^2} L^2

where L is the Reduced angular momentum operator. the minus sign is very important for giving a correct sign. the original angular momentum operator J is related by:

J=\hbar^2 L

by compare the Laplacian in spherical coordinate, the L is

L^2 = - \frac {1}{sin(\theta)} \frac {\partial}{\partial \theta} \left( sin(\theta) \frac {\partial}{\partial \theta} \right ) - \frac {1}{sin(\theta)} \frac{\partial^2} {\partial \phi ^2}

But this complicated form is rather useless, expect you are mathematic madman.

we can start from classical mechanic

\vec{L} = \vec {r} \times \vec{p}

L_x = y \frac {\partial} {\partial z} - z \frac {\partial}{\partial y }

L_y = z \frac {\partial} {\partial x} - x \frac {\partial}{\partial z }

L_z = x \frac {\partial} {\partial y} - y \frac {\partial}{\partial x }

with the change of coordinate

\begin {pmatrix} x \\ y \\ z \end{pmatrix} = \begin {pmatrix} r sin(\theta) cos(\phi) \\ r sin(\theta) sin(\phi) \\ r cos(\theta) \end{pmatrix}

and the Jacobian Matrix M_J , which is used for related the derivatives.

since

\frac {\partial}{\partial x} = \frac {\partial r}{\partial x} \frac {\partial} {\partial r} +\frac {\partial \theta}{\partial x} \frac {\partial} {\partial \theta}+\frac {\partial \phi}{\partial x} \frac {\partial} {\partial \phi}

\frac {\partial}{\partial y} = \frac {\partial r}{\partial y} \frac {\partial} {\partial r} +\frac {\partial \theta}{\partial y} \frac {\partial} {\partial \theta}+\frac {\partial \phi}{\partial y} \frac {\partial} {\partial \phi}

\frac {\partial}{\partial z} = \frac {\partial r}{\partial z} \frac {\partial} {\partial r} +\frac {\partial \theta}{\partial z} \frac {\partial} {\partial \theta}+\frac {\partial \phi}{\partial z} \frac {\partial} {\partial \phi}

which can be simplify

\nabla_{(x,y,z)} = M_J^T \nabla_{(r, \theta, \phi )}

M_J = \frac {\partial ( r, \theta, \phi) }{\partial (x,y,z)}

M_J^{\mu\nu} = \frac {\partial \mu}{\partial \nu}

then, we have

L_x = i sin(\phi) \frac {\partial }{\partial \theta} +i cot(\theta) cos(\phi) \frac { \partial }{\partial \phi}

L_y =-i cos(\phi) \frac {\partial }{\partial \theta} + i cot(\theta) sin(\phi) \frac { \partial }{\partial \phi}

L_z = - i \frac {\partial }{\partial \phi}

However, even we have the functional form, it is still not good.  we need the ladder operator

L_+ = L_x + i L_y = Exp(i \phi) \left( \frac {\partial }{\partial \theta} + i cot(\theta) \frac { \partial }{\partial \phi} \right)

L_- = L_x - i L_y = Exp(-i \phi) \left( \frac {\partial }{\partial \theta} - i cot(\theta) \frac { \partial }{\partial \phi} \right)

notice that

L_+^\dagger = L_-

so, just replacing i \rightarrow -i .

when we looking for the Maximum state of the spherical Harmonic Y_{max}(\theta, \phi)

L_+ Y_{max}(\theta,\phi) = 0 *)

use the separable variable assumption.

Y_{max}(\theta, \phi) = \Theta \Phi

L_+ \Theta \Phi = 0 = - Exp(i \phi) \left( \frac {d\Theta}{d \theta} \Phi + i cot(\theta) \frac { d\Phi}{d\phi} \right) \Theta

\frac {tan(\theta)}{\Theta} \frac { d \Theta} {d \Theta } = - \frac {i}{\Phi} \frac {d \Phi} {d \phi} = m

the solution is

Y_{max}(\theta,\phi) = sin^m(\theta) Exp(i m \phi )

L^2 Y_{max}(\theta, \phi) = m(m+1) Y_{max}(\theta,\phi)

an application on Hydrogen wave function is here.

Hydrogen Atom (Bohr Model)

Leave a comment

OK, here is a little off track. But that is what i were learning and learned. like to share in here. and understand the concept of hydrogen is very helpful to understand the nuclear, because many ideas in nuclear physics are borrow from it, like “shell”.

The interesting thing is about the energy level of Hydrogen atom. the most simple atomic system. it only contains a proton at the center, um.. almost center, and an electron moving around. well, this is the “picture”. the fact is, there is no “trajectory” or locus for the electron, so technically, it is hard to say it is moving!

why i suddenly do that is because, many text books said it is easy to calculate the energy level and spectrum for it. Moreover, many famous physicists said it is easy. like Feynman, Dirac, Landau, Pauli, etc… OK, lets check how easy it is.

anyway, we follow the usual say in every text book. we put the Coulomb potential in the Schrödinger equation, change the coordinate to spherical. that is better and easy for calculation because the coulomb potential is spherical symmetric. by that mean, the momentum operator (any one don’t know what is OPERATOR, the simplest explanation is : it is a function of function.) automatically separated into 2 parts : radial and angular part. The angular part can be so simple that it is the Spherical harmonic.

Thus the solution of the “wave function” of the electron, which is also the probability distribution of  the electron location, contains 2 parts as well. the radial part is not so trivial, but the angular part is so easy. and it is just Y(l,m) .

if we denote the angular momentum as L, and the z component of it is Lz, thus we have,

L^2 Y(l,m) = l(l+1) \hbar^2 Y(l,m)

L_z Y(l,m) = m \hbar Y(l,m)

as every quadratic operator, there are “ladder” operator for “up” and “down”.

L_\pm Y(l,m) =\hbar \sqrt{l(l+1) - m(m\pm 1)} Y(l,m \pm 1)

which means, the UP operator is increase the z-component by 1, the constant there does not brother us.

it is truly easy to find out the exact form of the Y(l,m) by using the ladder operator. as we know, The z component of the a VECTOR must have some maximum. so, there exist an Y(l,m) such that

L_+ Y(l,m) =0

since there is no more higher z-component.

by solve this equation, we can find out the exact form of Y(l,m) and sub this in to L2, we can knowMax(m) = l . and apply the DOWN operator, we can fins out all Y(l,m) , and the normalization constant is easy to find by the normalization condition in spherical coordinate, the normalization factor is sin(\theta) , instead of 1 in rectangular coordinate.

\int_0^\pi \int_0^{2 \pi} Y^*(l',m') Y(l,m) sin(\theta) d\theta d \psi = \delta_{l' l} \delta_{m' m}

more on here