Reduction of Matrix

Leave a comment

In previous post, the number of combination between the core \Psi_{J_A} and the single particle state phi_{nlj} becomes massive as the number of m-states are huge. To simplify the coupling, assuming the interaction between the state \Psi_{J_AM_A}\phi_{nljm} and \Psi_{J_A'M_A'}\phi_{n'l'j'm'} are the same for all M_A, m, M'_A, m'. The Hamiltonian matrix can be reduced.

To demonstrate the idea, suppose the core only has a ground state. And the single-particle state has degenerated state with different m-state, say, j = 0,1

Lets the coupled states be

\displaystyle \begin{pmatrix}\Psi_{00}\phi_{00} \\ \Psi_{00}\phi_{1,-1} \\  \Psi_{00}\phi_{1,0} \\  \Psi_{00}\phi_{1,1} \end{pmatrix}

The Hamiltonian matrix is

\displaystyle H = \begin{pmatrix} E_1 & V & V & V \\ V & E_2 & 0 & 0 \\ V & 0 & E_2 & 0 \\ V & 0 & 0 & E2 \end{pmatrix} \Rightarrow \begin{pmatrix} E_1 & \sqrt{3}V \\ \sqrt{3}V & E_2 \end{pmatrix}

To give a concrete example, say E_1 = 0, E_2 = 1, V= 1

\displaystyle H = \begin{pmatrix} 0 & 1 & 1 & 1 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \end{pmatrix}

The eigenstates and eigenvalues are

\displaystyle e_1 = 2.3028, v_1 = ( 0.6011, 0.4614, 0.4614, 0.4614)

\displaystyle e_2 = -1.3028, v_2 = ( -0.7992, 0.3470, 0.3470, 0.3470)

\displaystyle e_3 = 1.000, v_3 = ( 0.0000, 0.4082, 0.4082, 0.8165)

\displaystyle e_4 = 1.000, v_4 = ( 0.0000, 0.7071, 0.7071, 0.0000)

We can see, there are 2 states the energies do not change and the coefficient just re-configured.

Let’s look at the eigen system of this matrix:

\displaystyle H = \begin{pmatrix} 0 & 1 \\ 1 & 1  \end{pmatrix}

\displaystyle e_1 = 2.3028, v_1 = ( 0.6011, 0.7992)

\displaystyle e_2 = -1.3028, v_2 = ( -0.7992, 0.6011)

I wish to prove it more general as a theorem, but I can just worked on few cases.


Spectroscopic factor (revisit)

Leave a comment

For transfer reaction, the wave function of nucleus B = A + 1 can be written as

\displaystyle \Psi_{J_B m_B}(A+1) = \sum_{A'nlj} \beta_{nlj}(B,A) A[\Psi_{J_A'}(A) \phi_{nlj}(r)]_{J_B m_B} \\ = \sum_{A' nlj} \beta_{nlj}(B,A') \frac{1}{N_A} \sum_{m_A' m}C_{J_A M_A j m}^{J_B M_B} \Psi_{J_A' m_A'}(A) \phi_{nljm}(r)

where the N_A = _{A+1}C_1 due to anti-symetrization between the single nucleon and the core nucleus A, the Clebsh-Gordon coefficient is come from the angular coupling, and \beta_{nlj}(B,A) is the square root of the spectroscopic factor for orbital nlj between nuclei B and A.

We can see that, to find out the spectroscopic factor, we can integrate out the core nucleus A

\displaystyle \beta_{nlj}(B,A) = \langle \Psi_{J_B mB} | A [\Psi_{J_A} \phi_{nlj}]_{J_B m_B} \rangle \\  = \sum_{m_A m} C_{J_A m_A j m}^{J_B m_B} \langle J_B m_B |a_{nljm}^\dagger |J_A m_A \rangle \\ = \sum_{m_A m} (-1)^{J_A - J + m_B} \sqrt{2J_B+1} \begin{pmatrix} J_A & j & J_B \\ m_A & m & -m_B \end{pmatrix} \langle J_B m_B |a_{nljm}^\dagger |J_A m_A \rangle

In the last step, we write the Clebsh-Gordon coefficient in Winger 3j-symbol. Using the Wigner-Eckart theorem

\displaystyle \langle J_B m_B |a_{nljm}^\dagger |J_A m_A \rangle \\ = \frac{(-1)^{J_B-m_B}}{\sqrt{2J_B+1}} \begin{pmatrix} J_B & j & J_A \\ -m_B & m &  m_A \end{pmatrix} \langle J_B || a_{nljm}^\dagger ||J_A \rangle \\ = \frac{(-1)^{-m_B - J_A - j }}{\sqrt{2J_B+1}} \begin{pmatrix} J_A & j & J_B \\ m_A & m &  -m_B \end{pmatrix} \langle J_B || a_{nljm}^\dagger ||J_A \rangle


\displaystyle \beta_{nlj}(B,A) = \sum_{m_A m} (-1)^{-2J} \sqrt{2J_B+1} \begin{pmatrix} J_A & j & J_B \\ m_A & m & -m_B \end{pmatrix}^2 \langle J_B|| |a_{nljm}^\dagger ||J_A\rangle

using the identity

\displaystyle  \sum_{m_A m} \begin{pmatrix} J_A & j & J_B \\ m_A & m & -m_B \end{pmatrix}^2 = \frac{1}{2J_B+1}

we have

\displaystyle \beta_{nlj}(B,A) = \frac{(-1)^{2J}}{\sqrt{2J_B+1}} \langle J_B|| |a_{nljm}^\dagger ||J_A\rangle

replacing the reduced matrix element in the Wigner-Eckart theorem with previous result,

\displaystyle \langle J_B m_B |a_{nljm}^\dagger |J_A m_A \rangle = \frac{(-1)^{J_B-m_B}}{\sqrt{2J_B+1}} \begin{pmatrix} J_B & j & J_A \\ -m_B & m &  m_A \end{pmatrix} \frac{\sqrt{2J_B+1}}{(-1)^{2j}} \beta_{nlj}(B,A)

multiple the adjoint

\displaystyle (2J_B+1)\begin{pmatrix} J_B & j & J_A \\ -m_B & m &  m_A \end{pmatrix}^2 \beta_{nlj}^2(B,A) \\=\langle J_A m_A |a_{nljm} |J_B m_B \rangle \langle J_B m_B |a_{nljm}^\dagger |J_A m_A \rangle

For A(d,p) B reaction, the nucleus A is in the ground state, we can sum all the m_B and m states, then sum all the J_B , notice that

\displaystyle \sum_{J_B m_B} |J_B m_B \rangle \langle J_B m_B | = 1

we have

\displaystyle \sum_{J_B} \frac{2J_B+1}{2J_A+1} \beta_{nlj}^2(B,A) =\langle J_A m_A |\sum_{m} a_{nljm}a_{nljm}^\dagger |J_A m_A \rangle


a_{nljm}a_{nljm}^\dagger = 1 - a_{nljm}^\dagger a_{nljm} \\ \langle J_A m_A | \sum_{m} a_{nljm}^\dagger a_{nljm}|J_A m_A \rangle = n_{nlj}(A)

The last equality mean the number of nucleon in orbital nlj in nucleus A. Thus, we obtain

\displaystyle \sum_{J_B} \frac{2J_B+1}{2J_A+1} \beta_{nlj}^2(B,A) = 2j+1  - n_{nlj}(A)

For B(p,d)A reaction, the nucleus B is in the ground state and we can sum all J_A, m_A, m,

\displaystyle \sum_{J_A} \beta_{nlj}^2(B,A) = n_{nlj}(B)

Thus, the adding A(d,p)B and removing A(p,d)C reaction from nucleus A, the sum of spectroscopic factors

\displaystyle  \sum_{J_B} \frac{2J_B+1}{2J_A+1} \beta_{nlj}^2(B,A) + \sum_{J_C} \beta_{nlj}^2(A,C) = 2j+1



Expressions of Total Nuclear Hamiltonian

Leave a comment

In the past, I always write the total Hamiltonian of N+1 system as

\displaystyle H = H_N + H_1 + V_{N1}

The 3 terms on the right are the Hamiltonian of the N-nucleon core, the Hamiltonian of the single-nucleon, and the interaction between the core and the single particle. And the wave function can be written as

\displaystyle \chi_k = \sum \beta_{ij} \Psi_i \phi_j


\displaystyle H_N \Psi_i = \epsilon^{N}_i \Psi_i,  H_1 \phi_j = \epsilon^{1}_j \phi_j

so far so good.

Breaking down the above Hamiltonian,

\displaystyle H = \sum_{i}^{N+1} K_i + \sum_{i\neq j}^{N+1} V_{ij}

\displaystyle H_N = \sum_{i}^{N} K_i + \sum_{i\neq j}^{N} V_{ij}


\displaystyle H = H_N + K_1 + \sum_{i}^{N} V_{i1}


\displaystyle H_1  = K_1 + \sum_{i}^{N} V_{i1}  - V_{N1}

Thus, we can see, the term V_{N1} is NOT really the sum of all interaction between the single-nucleon and the nucleons in the core. What is it? or What is H_1?

Ideally, in the single particle picture (or independent particle model) , we imagine the single-nucleon is bounded by an mean field created by the core, i.e.

\displaystyle H_{N1} = H_N + K + \bar{V} = H_N + H_1

\displaystyle H_{N1} \Psi_i \phi_j = \epsilon^{N1}_{ij} \Psi_i \phi_j

Adding a 2-body residual interaction, we restore the total Hamiltonian,

H = H_N + H_1 + R,   R = V_{1N} .

Thus, the V_{1N} is the residual interaction between the single-nucleon and the core. This residual interaction, is created after the mean field!

\displaystyle H = H_N + K_1 + \bar{V} +  \sum_{i}^{N} V_{i1} - \bar{V} = H_N + H_1 + V_{1N}

\displaystyle  H_1 = K_1 + \bar{V}

\displaystyle  V_{1N} = \sum_{i}^{N} V_{i1} - \bar{V}





Physics behind Woods-Saxon energy levels

Leave a comment

In the last post, I hope I explained how to find the Woods-Saxon energy levels for a given parameters. I just searched the best fit Woods-Saxon parameters to best fit the neutron single-particle energy for  209Pb. This is a double magic nucleus, and the there is no large fragmentation for the neutron excitation energy, thus, the outermost neutron in this nucleus can be considered as a good single-particle.

The best fit parameters are

Screen Shot 2019-03-21 at 01.47.03.png

The rms difference is just 78 keV!! For each level, the difference is not more than 50 keV!

So, how can we interpret this fitting result? It means, the energy levels can be well explained by WS mean field.

Recall that, the mean field actually included a lot things, it is the effective single particle potential that a nucleon is feeling. When we pushing the rms value to be minium. We are actually finding the best mean field. And the difference is due to the residual interaction.

But what contained in the mean field and residual interaction? The total Hamiltonian is

\displaystyle H = H_{N} + H_{NN} + H_{NNN} + ....

the mean field approach is to add an artificial potential V

\displaystyle H = (H_{N} + V) + (H_{NN} + H_{NNN} +... - V)  = H_{m} + H_{R}

such that $latex  H_{R} $ is minimum.

And by fitting energy levels using Woods-Saxon energy levels, we are basically doing the same thing! And the mean field is explicitly containing 2-body force, 3N-force and so on. So, can we say “because the Woods-Saxon can explain the energy level very well, the tensor force or other force is insignificant” ? the answer is no. Because the Woods-Saxon potential explicitly contains tensor force and other force.

Calculation of Woods-Saxon energy levels

Leave a comment

Hello everyone, I have been very busy lately, working on many things at the same time. Writing experiment proposal, testing a new detector array, flighting with the noise. Coding a program to read digitizer data for real time. Update the analysis code for coming experiments. Update/making some calculation codes, one of them is the calculation of Woods-Saxon energy levels.

Using RK4 method, a 2ndary 1-D differential potential can be solved easily (with some effort).

The Woods-Saxon potential with spin-orbital coupling is:

\displaystyle H = \frac{P^2}{2m} + V(r) + \frac{1}{r} \frac{d}{dr}V_{so}(r)\langle L\cdot S\rangle

\displaystyle V(r) = \frac{V_0 }{1+ \exp{\frac{r-R_0}{a_0}}}

\displaystyle V_{so}(r) = \frac{V_{so}}{1+\exp{\frac{r-R_{so}}{a_{so}}}}

There are 6 parameters, 3 for the central term, 3 for the spin-orbital term.

And the equation is

\displaystyle H\Psi(r, \theta, \phi) = E \Psi(r, \theta, \phi)

Since the radial wavefunction has to be zero when it is far away, we can search for energy, such that the wavefunction is close to zero at a certain range. ( this is the only method I know how to practically implemented. I can imagine that, there may be some other method, transform the whole equation with a fixed r and directly solve the energy. Something like least-square method for linear system can be transformed into Matrix formalism and be directly solved. But I don’t know how)

The code is available at github (currently only for neutron).

I also have an excel to solve a single energy. Feel free to leave in comment and I will email it to you.

3j, 6j, 9j-symbol for C++

Leave a comment

Since I cannot find a odd c++ code ( < c++11 ) for the calculation of those j-symbol, I made one for myself!

The code just using the formula which you can find in wiki or here. So, the code is not optimized, but for simple calculation, I guess it does not matter. Also, I checked, in my pc, calculate 1000 9-j symbols only take 1 or 2 sec.

The program first build

  1. factorial
  2. Clebsch–Gordan coefficient
  3. 3-j symbol
  4. 6-j symbol
  5. 9-j symbol

This is not a good place (or I don’t know how to) for code display…. any way.

#include <stdlib.h>
#include <cmath>

double factorial(double n){
if( n < 0 ) return -100.;

return (n == 1. || n == 0.) ? 1. : factorial(n-1) * n ;

double CGcoeff(double J, double m, double J1, double m1, double J2, double m2){
// (J1,m1) + (J2, m2) = (J, m)

if( m != m1 + m2 ) return 0;

double Jmin = abs(J1 - J2);
double Jmax = J1+J2;

if( J < Jmin || Jmax < J ) return 0;

double s0 = (2*J+1.) * factorial(J+J1-J2) * factorial(J-J1+J2) * factorial(J1+J2-J) / factorial(J+J1+J2 + 1.);
s0 = sqrt(s0);

double s = factorial(J +m ) * factorial(J -m);
double s1 = factorial(J1+m1) * factorial(J1-m1);
double s2 = factorial(J2+m2) * factorial(J2-m2);
s = sqrt(s * s1 * s2);

//printf(" s0, s = %f , %f \n", s0, s);

int kMax = min( min( J1+J2-J, J1 - m1), J2 + m2);

double CG = 0.;
for( int k = 0; k <= kMax; k++){
double k1 = factorial(J1+J2-J-k);
double k2 = factorial(J1-m1-k);
double k3 = factorial(J2+m2-k);
double k4 = factorial(J - J2 + m1 +k);
double k5 = factorial(J - J1 - m2 +k);
double temp = pow(-1, k) / (factorial(k) * k1 * k2 * k3 * k4 * k5);
if( k1 == -100. || k2 == -100. || k3 == -100. || k4 == -100. || k5 == -100. ) temp = 0.;

//printf(" %d | %f \n", k, temp);
CG += temp;

return s0*s*CG;


double ThreeJSymbol(double J1, double m1, double J2, double m2, double J3, double m3){

// ( J1 J2 J3 ) = (-1)^(J1-J2 - m3)/ sqrt(2*J3+1) * CGcoeff(J3, -m3, J1, m1, J2, m2) 
// ( m1 m2 m3 )

return pow(-1, J1 - J2 - m3)/sqrt(2*J3+1) * CGcoeff(J3, -m3, J1, m1, J2, m2);


double SixJSymbol(double J1, double J2, double J3, double J4, double J5, double J6){

// coupling of j1, j2, j3 to J-J1
// J1 = j1
// J2 = j2
// J3 = j12 = j1 + j2
// J4 = j3
// J5 = J = j1 + j2 + j3
// J6 = j23 = j2 + j3

//check triangle condition
if( J3 < abs(J1 - J2 ) || J1 + J2 < J3 ) return 0; 
if( J6 < abs(J2 - J4 ) || J2 + J4 < J6 ) return 0;
if( J5 < abs(J1 - J6 ) || J1 + J6 < J5 ) return 0;
if( J5 < abs(J3 - J4 ) || J3 + J4 < J5 ) return 0;

double sixJ = 0;

for( float m1 = -J1; m1 <= J1 ; m1 = m1 + 1){
for( float m2 = -J2; m2 <= J2 ; m2 = m2 + 1){
for( float m3 = -J3; m3 <= J3 ; m3 = m3 + 1){
for( float m4 = -J4; m4 <= J4 ; m4 = m4 + 1){
for( float m5 = -J5; m5 <= J5 ; m5 = m5 + 1){
for( float m6 = -J6; m6 <= J6 ; m6 = m6 + 1){

double f = (J1 - m1) + (J2 - m2) + (J3 - m3) + (J4 - m4) + (J5 - m5) + (J6 - m6);

double a1 = ThreeJSymbol( J1, -m1, J2, -m2, J3, -m3); // J3 = j12 
double a2 = ThreeJSymbol( J1, m1, J5, -m5, J6, m6); // J5 = j1 + (J6 = j23)
double a3 = ThreeJSymbol( J4, m4, J2, m2, J6, -m6); // J6 = j23
double a4 = ThreeJSymbol( J4, -m4, J5, m5, J3, m3); // J5 = j3 + j12

double a = a1 * a2 * a3 * a4;
//if( a != 0 ) printf("%4.1f %4.1f %4.1f %4.1f %4.1f %4.1f | %f \n", m1, m2, m3, m4, m5, m6, a);

sixJ += pow(-1, f) * a1 * a2 * a3 * a4;


return sixJ;

double NineJSymbol( double J1, double J2, double J3, double J4, double J5, double J6, double J7, double J8, double J9){

double gMin = min( min (min( abs(J1 - J2 ), abs(J4 - J5)) , abs( J4 - J6 )) , abs(J7 - J8));
double gMax = max( max (max( J1+J2, J4+J5), J3+J6), J7+J8);

//printf(" gMin, gMax = %f %f \n", gMin, gMax);

double nineJ = 0;
for( float g = gMin; g <= gMax ; g = g + 1){
double f = pow(-1, 2*g) * (2*g+1);
double s1 = SixJSymbol(J1, J4, J7, J8, J9, g);
if( s1 == 0 ) continue;
double s2 = SixJSymbol(J2, J5, J8, J4, g, J6);
if( s2 == 0 ) continue;
double s3 = SixJSymbol(J3, J6, J9, g, J1, J2);
if( s3 == 0 ) continue;
nineJ += f* s1*s2*s3;

return nineJ;

Single-particle structure of 19F

Leave a comment

About one year ago, I studied the nuclear structure of 19F. At that time, a lot of thing don’t understand and confused.  But now, I have a better understanding.

19F is always fascinating because the so complicated energy levels for such a relatively simple system with only 1 proton, 2 neutrons on top of a 16O core, which usually treated as double magic rigid core.

  • It ground state is 1/2+, which is unusual.
  • It is also well-deformed with \beta_2 \sim 0.4 , deduced from rotational band.
  • It also has very low lying negative parity state of 0.110 MeV. This state is also the band head of K=1/2- rotational band.
  • The rotational band of K=1/2± does not following J(J+1) curve (see the last picture from this post), this indicate the rigid rotor assumption is not so good.

3 of the rotational bands of 19F, from M. Oyamada et al., PRC11 (1975) 1578

When looking the region around 19F. From 16O a double magic core, to 20Ne very deformed nucleus (\beta_2 \sim = 0.7 ), and then go to 32Mg, the center of island of inversion. 16O, 18O, 18F, 20F are normal nuclei with normal shell order. 19F is a doorway nuclei for the complicated nuclear structure. That make the understanding the nuclear structure of 19F important.

  1. it helps to understand how deformation happen in small system,
  2. and how deformation can be described in particle-configuration.

The single-particle structure of the ground state of 19F has been studied using (p,2p) reaction. And found that the ground state to ground state transition only contain strength from 1s1/2 orbital. Thus, the wave function of 19F must bein the form

\displaystyle |^{19}F\rangle = \alpha |\pi 1s_{1/2} \rangle |^{18}O\rangle + \beta |\pi 0d_{5/2}\rangle |^{18}O^* \rangle + ...

The single-particle structures of the ground state and excited states are best studied using single-nucleon transfer. There are at least 4 directions, neutron-removal from 20F, neutron-adding from 18O, proton-removal from 20Ne, and proton-adding from 18O.  Also, a proton/neutron-removal from 19F itself to study the ground state properties.


In above figure, we also plotted the spectroscopic factors from 2 reactions, which is taken from G. Th. Kaschl et al., NPA155(1970)417, and M. Yasue et al., PRC 46 (1992) 1242. These are the most significant studies. It was suprised that the negative parity state of 0.110 meV can be populated from the 18O(3He,d) reaction. These suggest the 18O is not a good core that it has about 10% strength of two-nucleon hole in its ground state. From the 20Ne(d, 3He) reaction, the 0.110 MeV state is highly populated. Given that the ground state of 20Ne is mainly proton 1s1/2 strength due to deformation. These result indicated that the negative parity state of 19F is due to a proton-hole. Thus we have a picture for the 0.110 MeV state of 19F:


The 18O(3He,d) reaction put a 0p1/2 proton in 18O 2-proton hole state to form the negative parity state. And 20Ne(d, 3He) reaction remove a 0p1/2 proton from 20Ne.

But this picture has a problem that, how come it is so easy to remove the 0p1/2 proton from 20Ne? the p-sd shell gap is known to be around 6 MeV! or, why the proton hole in 19F has such a small energy?

In current understanding, the 2-neutron are coupled to J=0 pair, and no contribution to low-lying states. But is it true?

I suspect, the underneath reason for the proton 1s1/2 orbital is lower is because the proton 0d5/2orbital was repealed by the 2-neutron in 0d5/2 orbital due to the tensor force. And somehow, the tensor force becomes smaller in 20F when the 0d5/2 orbital is half filled.

And because the proton is in 1s1/2 and the 0d5/2 is just 0.2 MeV away, a huge configuration mixing occurred. and then, a Nilsson orbit is formed with beta = 0.4. This is an example of NN-interaction driving deformation.

Older Entries