Phase shift of elastics scattering

Leave a comment

I found that the derivation of most “google result” is not clear enough. So here is my derivation. Before process, people may need to review the pervious post.

Most people start from,

u_l = A( \hat{h}_l^- - S_l \hat{h}_l^+ )

that annoying me because it is somehow not “natural”. Why there is a “minus” sign? Why the \hat{h}_l^- is the first term? For my self, a more natural way is,

u_l = a \hat{h}_l^+ + b \hat{h}_l^-

where a, b are complex numbers, but that is still not so natural, because in numerical calculation, for simplicity, there is no complex number, we only have,

u_l = \alpha \hat{j}_l + \beta \hat{n}_l

The first term is alway there as it is the free solution and bounded at r = 0. the second term is caused by the potential.

The goal is to find a solution take the form

\displaystyle \psi = A \left( e^{i \vec{k} \cdot \vec{r}} + f(\theta) \frac{e^{ikr}}{r} \right)

where the first term is free wave and the second term is scattered wave. The solution for elastics scattering is

\displaystyle \psi = \sum C_l P_l (\cos\theta) \frac{u_l}{kr} \rightarrow \sum C_l P_l(\cos\theta) (\alpha \hat{j}_l + \beta \hat{n}_l)

we used the substitution,

\displaystyle R_l(r) = \frac{u_l(\rho)}{\rho},  \rho = kr.

The radial function can be solved using Rungu-Kutta method on the equation,

\displaystyle \frac{d^2}{d\rho^2} u_l = \frac{2 m_\mu}{\hbar^2} (V-E) u_l + \frac{l(l+1)}{\rho^2}

and the solution of u_l at far away is,

u_l \rightarrow  \alpha \hat{j}_l + \beta \hat{n}_l.

the arrow means r \rightarrow \infty. So, the problem is how to rewrite the solution. In the way, we will see how the phase shift or the S-matrix was found.

The free solution is the spherical wave,

\displaystyle e^{i \vec{k} \cdot \vec{r}} = \sum_l (2l+1) i^l P_l(\cos\theta) j_l(kr)

The spherical Bessel function j_l(kr) cna be express as Heankel function

h_l^{\pm} = n_l \pm i j_l \rightarrow e^{\pm i (kr - l \frac{\pi}{2})}

The + sign is outgoing wave.

\displaystyle u_l \rightarrow (\alpha \hat{j}_l + \beta \hat{n}_l)

\displaystyle = \frac{\alpha}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \frac{\beta}{2}(\hat{h}_l^{+} + \hat{h}_l^{-})

\displaystyle = \frac{\alpha + i \beta}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \beta \hat{h}_l^{+}

\displaystyle = (\alpha - i \beta ) \left( \frac{\hat{h}_l^+ - \hat{h}_l^-}{2i} + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

\displaystyle = (\alpha - i \beta ) \left( \hat{j}_l + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

Since the u_l should be normalized, we can se \alpha = \cos \delta and \beta = \sin\delta.

\displaystyle \frac{\beta}{\alpha - i \beta } = \sin(\delta) e^{i\delta}

We put u_l back

\displaystyle \psi \rightarrow \sum_l C_l P_l (cos\theta)(\alpha - i \beta ) \left( \hat{j}_l + \sin(\delta) e^{i\delta} \hat{h}_l^+\right)

By setting

\displaystyle C_l = A i^l \frac{2l+1}{\alpha - i \beta} ,

we have the first term is the free wave function. In the second term, \hat{h}_l^+ \rightarrow e^{i(kr - l \frac{\pi}{2}}) / kr . Notice that

e^{i l \frac{\pi}{2}} = i^{-l}

That cancel the i^l term in C_l . And we have

 \displaystyle f(\theta) = \sum (2l+1) P_l (\cos\theta) \frac{\sin(\delta) e^{i\delta}}{k}

some people will write the u_l as \hat{h}_l^{\pm} and the S-matrix,

\displaystyle u_l = \frac{\alpha + i \beta} {2i} \hat{h}_l^+ - \frac{\alpha - i \beta}{2i} \hat{h}_l^-

\displaystyle = -\frac{\alpha - i \beta}{2i} \left( \hat{h}_l^- - \frac{\alpha + i \beta}{\alpha - i \beta} \hat{h}_l^+ \right)

\displaystyle = A' (\hat{h}_l^- - S_l \hat{h}_l^+)


\displaystyle S_l =\frac{\alpha + i \beta}{\alpha - i \beta} = e^{2i\delta} .

Remember that this is the S-matrix for elastics scattering.


Solving radial SE numerically


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.


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.


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.


Cross seciton of Coulomb Scattering

Leave a comment

The coulomb scattering looks very easy, the formula of the differential cross section in CM frame is,

\frac{d\sigma}{d\Omega} = (\frac{Z1 Z2 e}{4 E_{cm}})^2 \frac{1}{sin^4(\theta_{cm}/2)},

where e = 1.44 MeV fm and 1 fm^2 = 10 mb. The tricky point is, in most experiment, we are working in Laboratory frame that require frame transformation.

The relationship of the energy in CM frame and the energy in the Lab frame can be found by Lorentz transform, and use the total kinematic energy (both particle 1 and particle 2). In the CM frame, we can image we have a fixed virtual target on the center of mass, and there is only 1 object moving at energy of the total kinematic energy.

For example, we have a target of mass m_2, a projectile with mass m_1 and energy T_1, a classical energy in CM frame is

E_{cm} = \frac{m_2}{m_1+m_2} T_1

In fact the E_{cm} has only 5% difference between relativistic and non-relativistic even up to 500 MeV

When calculating the integrated cross section, we can do it in the CM frame, but it is more intuitive to do it in the Lab frame. In this case, we need to transform the differential cross section from the CM frame to the Lab frame, mathematically, the transformation is done by a factor called Jacobian.

We can compare the result using the kinematic calculator in LISE++.

Screenshot from 2016-04-22 01:08:55.png

In the above plot, the blue line is the d.s.c. in CM frame, and the red line is d.s.c. of the 9Be in Lab frame. Jacobian was added, therefore, the zero degree d.s.c. of 9Be is larger than the 180 degree d.s.c. in the CM frame.

Screenshot from 2016-04-22 01:25:11.png

The grazing angle of the scattering, can be calculated by the shorted distance between the target and the projectile. In the Lab frame, the target is not fixed, so it is not easy to know the shortest distance. But in the CM frame, the virtual target is fixed, and we can calculate the distance using the E_{cm} and \theta_{cm}.

Very short introduction to Partial-wave expansion of scattering wave function

Leave a comment

In a scattering problem, the main objective is solving the Schrödinger equation


where H is the total Hamiltonian of the scattering system in the center of momentum, K is the kinetic energy and V is the potential energy. We seek for a solution \psi,

\displaystyle \psi_{k}^{+}(r)=e^{i\vec{k}\cdot \vec{r}}+f(\theta)\frac{e^{ikr}}{kr}

The solution can be decomposed

\displaystyle \psi_{k}^{+}(r)=R_{l}(k,r)Y_{lm}(\theta,\phi)=\frac{u_{l}(k,r)}{kr}Y_{lm}(\theta,\phi)

The solution of u_{l}(k,r) can be solve by Runge-Kutta method on the pdf

\displaystyle \left(\frac{d^2}{d\rho^2} + 1 - \frac{l(l+1)}{\rho^2} \right)u_{l}(k,\rho)=U(\rho)u_{l}(k,\rho)

where \rho=kr, k=\sqrt{2\mu E}/\hbar, \mu=(m_1+m_2)/(m_1 m_2) and U=V/E.

For U = 0, the solution of u_l is

\displaystyle u_{l}(k,r)=\hat{j}_l(\rho) \xrightarrow{r\rightarrow \infty} \sin(r') = \frac{e^{ir'}-e^{-ir'}}{2i}

where r' = kr-l\pi/2 and \hat{j}_l is the Riccati-Bessel function. The free wave function is

\displaystyle \phi_k(r)=e^{i\vec{k}\cdot\vec{r}}=\sum\limits_{l=0} P_l(\cos(\theta)) \frac{2l+1}{2ikr}i^l (e^{ir'}-e^{-ir'})

where P_l(x) is the Legendre polynomial.

Note that, if we have Coulomb potential, we need to use the Coulomb wave instead of free wave, because the range of coulomb force is infinity.

For U\neq 0, the solution of u_l(r<R) can be found by Runge-Kutta method, where R is a sufficiency large that the potential V is effectively equal to 0.  The solution of u_l(r>R) is shifted

\displaystyle u_{l}(k,r>R)=\hat{j}_l(\rho)+\beta_l \hat{n}_l(\rho) \xrightarrow{r\rightarrow \infty} \frac{1}{2i}(S_l e^{ir'}-e^{-ir'})

where S_l is the scattering matrix element, it is obtained by solving the boundary condition at r = R. The scattered wave function is

\displaystyle \psi_k(r)=\sum\limits_{l=0} P_l(\cos(\theta)) (2l+1) i^l \frac{u_l(r)}{kr}

put the scattered wave function and the free wave function back to the seeking solution, we have the f(\theta)

 \displaystyle f(\theta) = \sum\limits_{l=0} P_l(\cos(\theta)) \frac{2l+1}{2ik} (S_l - 1)

and the differential cross section

\displaystyle \frac{d\sigma}{d\Omega}=|f(\theta)|^2.

In this very brief introduction, we can see

  • How the scattering matrix S_l is obtained
  • How the scattering amplitude f(\theta) relates to the scattering matrix

But what is scattering matrix? Although the page did not explained very well, especially how to use it.

GEANT4 installation – Concept

Leave a comment

for the moment, i think GEANT4 is a set of libraries (of physical process, meterial) and a tool (Monte Carlo method).

As you can see, Visualization is not included in GEANT4. I mean, GEANT4 can output some files for other Visualization programs, but those programs have to be install separately.

A complete environment = GEANT4 + Visualization 

//——————– GEANT4 installation

To install GEANT4, you can go to to download to source code.

The installation method can be found in

basically, unzip the source code.

tar -zxvf geant4-source
mkdir geant4-build
cd geant4-build
cmake -DGEANT4_INSTALL_DATA=ON -DGEANT4_USE_NETWORKVRML=ON -DCMAKE_INSTALL_PREFIX=/path/to/geant4-install  geant4-source
make install


There are two options, 1 is the GEANT4_INSTALL_DATA, this will download and install the physical process for radiation-matter interaction. 2 is the GEANT4_USE_NETWORKVRML, this will enable network VRML visuallization. This is a local install, all the installed files are located at /path/to/geant4-install. After the installation, better to write this in .bashrc

cd ~/geant4.10.02-install/bin
cd ~/


The only difficulty is the CMAKE 3.0 or above. But all errors i met can be easily found in google.

//————— Visualization program.

I skipped the openGL, as it is well known for unfriendly for NVidia display card.

I tried DAWN and freeWRL. DAWN is a static visualization. It read *.prim files and output an *.eps file for fixed angle. freeWRL is a VRML display using java-script, an alternative is openVRML. But I cannot install openVRML.

DAWN can be installed from, I have zero problem.

FreeWRL is a bit tricky. I can only install the version, not earlier, not later. I don’t know why. It require javascript engine, apt-get install libmozjs185-dev. when it said any thing missing, just look for libXXXX.

//——————————– Example.

I took the example form the geant4-source/example/basic/B1, copy it to ~/geant4Works/B1.

mkdir B1-build
cd B1-build
cmake -DGeant4_DIR=/path/to/geant4-install/lib/Geant4-10.2.0  B1

to run, you simple ./exampleB1

the macro files are *.mac. vis.mac is visualization macro, and will be called by init_vis.mac. You should go to vis.mac and modify it.

To run,

Idle> /control/execute run1.mac


I modified the last line be /run/beamOn 100, the result is

Screenshot from 2016-01-27 13:46:48.png

you can see there are 100 protons passing through. I don’t understand the example B1, I’m just showing you what is a proper result.


Scattering Matrix

Leave a comment

at the point of scattering ( t = 0 ), the wave function and the incoming and out-going wave function can be related as:

\left| \psi_{in} \right> \rightarrow \left| \psi \right> \leftarrow \left| \psi_{out}\right>

where the incoming and out-going wavefunction is very far away from the field of the scatter, and thus, they are free and we say they are asymptotic.

Let a time propagator with a full Hamiltonian be U(t), and a time propagator with free Hamiltonian be U^0(t) . thus the time behavior of the scattering wave functions can be related as:

U(t) \left| \psi \right> \rightarrow U^0(t) \left| \psi_{in} \right> for t \rightarrow - \infty

U(t) \left| \psi \right> \rightarrow U^0(t) \left| \psi_{out} \right> for t \rightarrow + \infty

in equation:

\left| \psi \right> = lim_{t \rightarrow - \infty} U^\dagger (t) U^0(t) \left| \psi_{in} \right> = \Omega_+ \left| \psi_{in} \right>

\left| \psi \right> = lim_{t \rightarrow + \infty} U^\dagger (t) U^0(t) \left| \psi_{out} \right> = \Omega_- \left| \psi_{out} \right>

Thus, we have the scattering matrix or the S-matrix.

\left| \psi_{out} \right> = \Omega_-^\dagger \Omega_+ \left|\psi_{in} \right> = S\left| \psi_{in} \right>

Now assume for a particular state generated by an accelerator is Φ, and a particular out-going asymptotic state is χ. we have:

\left| \phi \right> = \Omega_+ \left| \phi_+ \right>

\left| \chi \right> = \Omega_+ \left| \chi_- \right>

thus, the probability amplitude for the scattering between these 2 states is:

\omega( \chi \leftarrow \phi) = | \left<\chi_- | \phi_+ \right> |^2 =| \left<\chi|\Omega_-^\dagger \Omega_+| \phi \right> |^2 = |\left< \chi | S | \phi \right>|^2

if we expand a wave function  in momentum basis:

\left| \psi \right> = \int d^3p \left| p \right> \left< p | \psi \right>

\psi_{out}(p) = \int d^3p' \left<p | S| p' \right> \psi_{in}(p')

now, we are going to show the energy conservation of the scattering operator or matrix, by showing that the scattering operator S commutes with the Hamiltonian. from

Exp(\frac{i}{\hbar} H \tau) \Omega_\pm =\Omega_\pm Exp( \frac{i}{\hbar} H^0 \tau )

differential it then we have :

H \Omega_\pm = \Omega_\pm H^0


H^0 = \Omega_\pm^\dagger H \Omega_\pm

S H^0 = H^0 S

together with the wavefunction:

0 = \left<p'| [H^0,S] | p\right> = (E_{p'} -E_p ) \left< p'|S|p\right>

thus implies,

\left<p'|S|p\right> = \delta(E_{p'}-E_p) g( p' \leftarrow p)

since, at the forward direction,  the change of momentum is zero, we can write S = 1 + R, then,

\left<p'|S|p\right> =\delta(p-p') - 2 \pi i \delta(E_{p'}-E_p) t( p' \leftarrow p)

the t(p' \leftarrow p) is called on-shell T-matrix. since the energy must be equal, required by the delta function, thus, the momentum magnitude must be equal, therefore, the 2 momentums s on a shell. The T-matrix also related to the scattering amplitude by:

f(p' \leftarrow p) = - (2 \pi)^2 m t(p' \leftarrow p)

then the S-matrix becomes,

\left<p'| S|p \right> = \delta(p - p') + \frac{i}{2\pi m} \delta(E_{p'} - E ) f(p' \leftarrow p)

Schrödinger Equation and Scatter Equation

Leave a comment

The Spatial part of the Time-Independent Schrödinger Equation (TISE) is

\left(- \frac {\hbar^2}{2m} \nabla^2 +V(r) \right) \psi(r) = E \psi(r)

by setting

k^2 = 2m E / \hbar^2   and    U(r) = 2 m V(r) /\hbar^2

the equation becomes a wave equation with a source, or scattering equation.

(\nabla^2 +k^2) \psi(r) = U(r) \psi(r)

for solving it, we have to find the Green function such that

(\nabla^2 + k^2 ) G(r,r') = \delta(x-x')

The solution is easy ( i will post it later, you can think the Green function is the inverse of the Operator)

G(r,r') = - \frac{1}{4\pi} \frac {Exp( \pm i k \cdot (r-r') )} { |r - r'| }

the particular solution is

\psi_p(r) = \int {G(r,r') U(r') \psi(r') dr'}

plus the homogeneous solution

\psi(r) = \phi(r) +\int {G(r,r') U(r') \psi(r') dr'}

but it is odd that the solution contain itself!

Older Entries