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.

2-state system (an invert problem)

Leave a comment

Given 2 states with energy E_1 < E_2 , and interaction energy between the two states V , the Hamiltonian of the system is

\displaystyle H = \begin{pmatrix} E_1 && V \\ V && E_2 \end{pmatrix}

This system is already discussed in this post. The solution state in here again as:

The eigen-energies

\epsilon_{\pm} = \bar{E} \pm \sqrt{dE^2 + V^2 }

the eigen-states are

\Phi_{+} = \alpha \phi_1 + \beta \phi_2

\Phi_{-} = -\beta \phi_1 + \alpha \phi_2


\displaystyle \alpha = \frac{dE + \sqrt{dE^2 + V^2 }}{\sqrt{(dE+\sqrt{dE^2+V^2})^2 + V^2}}, \beta = \frac{V}{\sqrt{(dE+\sqrt{dE^2+V^2})^2 + V^2}}

The \phi_{1,2} are the wavefunction of pure state of E_{1,2} .

This can be easily extended to 3-state system or n-state system. The number of coupling constants is n(n-1)/2.

Experimentally, we observed the mixed states and the spectroscopic factors originated from one of the pure state. For example, in a neutron transfer reaction, the transferred neutron may coupled with the 2+ state of the core. The neutron is sitting in a pure orbital (\phi_1), say 1d5/2 orbital, it couple with the 2+ state (\phi_2) with interaction energy V . The final states will have energy \epsilon_{\pm} and we extract the spectroscopic factors (\alpha^2, \beta^2) for the 1d5/2 orbital $.

To state more clear, in a transfer reaction, we may observed two excited states A_+, A_- with transfer of |0p_{1/2}\rangle neutron. This neutron assumed to couple with a core state |C\rangle and a core excited state |C^+\rangle . And we supposed that the 2 observed states are:

|A_+ \rangle = \alpha |0p_{1/2}\rangle |C\rangle + \beta |0p_{1/2} \rangle |C^+\rangle  = \alpha |\phi_1\rangle + \beta |\phi_2\rangle

|A_- \rangle = \beta  |0p_{1/2}\rangle |C\rangle - \alpha |0p_{1/2} \rangle |C_+\rangle = \beta |\phi_1\rangle - \alpha |\phi_2\rangle

We then want to find out the un-perturbed energy E_{1,2} and the interaction V . And this is actually an easy problem by the diagonalization process.

The Hamiltonian was diagonalized into

\displaystyle H = \begin{pmatrix} E_1 && V \\ V && E_2 \end{pmatrix} =  P \cdot D \cdot P^{-1}


\displaystyle P = \begin{pmatrix} \alpha && -\beta \\ \beta && \alpha \end{pmatrix} , det(P) = 1

\displaystyle D = \begin{pmatrix} \epsilon_+ && 0 \\ 0 && \epsilon_- \end{pmatrix}


E_1 = \alpha^2 \epsilon_+ \beta^2 \epsilon_-

E_2 = \beta^2 \epsilon_+ \alpha^2 \epsilon_-

\displaystyle V = \alpha \beta (\epsilon_+ - \epsilon_-)  = \frac{1}{2}\sqrt{d\epsilon^2 - dE^2}


A thought on magnetic field

Leave a comment

Imagine there is a wire with some current on it. How do you know how the magnetic field look like without any knowledge but symmetry?


First of all, the field must be cylindrical, as the system is cylindrical.

Suppose we know it is rotate around the wire, but which direction?

of course, it is anti-clockwise, according to experiment.

But wait, does it somehow break the mirror symmetry?

Also, the so-called anti-clockwise is, in fact, a matter of notation, we DEFINED the magnetic field line is from north to south. In fact, there is nothing moving from north to south!!!

In my point of view, field line is not the fundamental quantity but the vector potential. The vector potential is same as the current direction and has the rotational symmetry and also mirror symmetry. The magnetic field is a kind of concept that simplify the picture.

So, when we put a small magnet, instead of the magnet will aligned with the field, may be think of the magnet align with the vector potential? just like the electric charge move along the potential. The magnetic charge moves perpendicular the the vector potential.

I don’t have a clear answer for this problem. But let me quote,

We think of the electric field as a condition in space set up by the system of point charges. . . .  The electric field is more than a calculational device. This concept enables us to avoid the problem of action at a distance . . .  We thus think of the force exerted on a charge q0 at point P as being exerted by the field at point P rather than by the charges, which are some distance away. Of course the field at point P is produced by the other charges, but not instantaneously.

— Paul Tipler

What exactly is an electric field? I have deliberately begun with what you might call the “minimal” interpretation of E, as an intermediate step in the calculation of electric
forces. But I encourage you to think of the field as a “real” physical entity, filling the space around electric charges. Maxwell himself came to believe that electric and magnetic fields are stresses and strains in an invisible primordial jellylike “ether”. Special relativity has forced us to abandon the notion of ether, and with it Maxwell’s mechanical interpretation of electromagnetic fields. (It is even possible, though cumbersome, to formulate classical electrodynamics as an “action-at-a-distance” theory, and dispense with the field concept altogether.) I can’t tell you, then, what a field is—only how to calculate it and what it can do for you once you’ve got it.

— David Griffiths

Rotational Band

Leave a comment

For deformed nuclei, it can be rotated in various angular momentum in Laboratory frame. Assume rigid body rotation, the energy is

\displaystyle E_J = \frac{1}{2}I\omega^2 = \frac{1}{2I}J^2

In QM, that becomes

\displaystyle H = \sum_{i=1}^{3} \frac{\hbar^2}{2I_i} J_i^2

For axial symmetry, I_1 = I_2 = I

\displaystyle H = \frac{\hbar^2}{2I} (J^2 - J_3^2) + \frac{\hbar^2}{2I_3}J_3^2

Remember, in deformed nuclei, the projection of J along the symmetry axis in the body-frame is K . The expected value of the Hamiltonian with state |Nn_z m_l K \rangle in the body-frame is proportional to J(J+1) for J^2 and K for J_3. i.e.

\displaystyle E_J = \frac{\hbar^2}{2I} J(J+1) + E_K

From body-frame to Lab-frame, we should apply the Wigner D-Matrix to the intrinsic wave function. ( I am not sure the following equation is correct, but the idea is rotating the body-frame wavefunction with Wigner D-Matrix to get the Lab-frame wave function. In Lab frame the total angular momentum must be a good Quantum number as rotational symmetry restored, so as J_z = M. The problem of the following equation is that the J is not a good Q-number in Nilsson wavefunction )

\displaystyle |JMK\rangle = \sum_{M} D_{MK}^{J} |Nn_zm_lK\rangle

However, the Wigner D-Matrix does not conserve parity transform:

\displaystyle D_{MK}^J \rightarrow (-1)^{J+K} D_{M-K}^{J}

In order to restored the parity, we need to include \pm K in the Lab-frame wave function.

\displaystyle |JMK\rangle = \sum_{M} \left( D_{MK}^J \pm (-1)^{J+K} D_{M-K}^J \right) |Nn_zm_lK\rangle

where + for positive parity, – for negative parity.

From the above equation, for K^\pi = 0^+ (0^-), J must be even (odd). For K > 0 , J = K, K+1, K+2, ... .

rotaional Band of 205Fm.png

rotational band of 253No.png

rotational band of 19F.png

We can see for K = 1/2 , the J = 5/2, 9/2, 11/2 are lower to the main sequence. This was explained by adding an extra term in the rotation Hamiltonian that connect \Delta K = 1 .

\displaystyle \langle JMK | H'(\Delta K = 1) |JMK \rangle

\displaystyle \rightarrow \langle D_{MK}^J | H' | D_{MK}^J \rangle+ \langle D_{M-K}^J |H'| D_{MK}^J \rangle + \langle D_{MK}^J | H' | D_{M-K}^J \rangle+ \langle D_{M-K}^J |H'| D_{M-K}^J \rangle

The term with \Delta K = 0 vanished. And since $\latex K = 1$, the only non-zero case is K = 1/2 .

A possible form of the H' (\Delta K = 1) = \frac{1}{2} \omega (J_+ + J_-) . These are the ladder operator to rise or lower the m-component by 1. In 19F case, we can think it is a single proton on top of 18O core.  A rotation core affect the proton with an additional force, similar to Coriolis force on earth.


some thoughts on mass

Leave a comment

Lets forget about the Higg mechanism.

Lets start from deuteron, it contains a proton and neutron, which bounded by strong nuclear force. The kinematics energy between the proton and neutron is about 35 MeV. Since the binding energy of deuteron is 2.2 MeV, thus the strong nuclear force is about -37 MeV. According to Einstein E = mc^2 . This 2.2 MeV is converted to 2.2 MeV/c2 of mass.

In this example, when a system contains internal motion, internal interaction, the “mass” is not just the mass of individual components, but also included the interaction energy.

Thus, what is mass? Only elementary particles had “pure” mass. The mass of all other particles must contains the energy of internal interaction and motion.

So, for the Earth, we know the mass from the orbital motion. But this mass, when the Earth exploded into tiny many pieces, some of the mass was gone because the internal interaction and motion no longer there.

What is the mass of the solar system? We can add all the masses of the sun, the planets, comets, ices, debris, etc, but we also need to add the internal motion and interaction. Although the gravitation force is weak, but still, we need to add this.

What is the mass of the Universe? Again, we can added all visible mass, but we also need to add the dark energy, the dark matter ( I suspect the dark matter is kind of illusion that we don’t fully understood the gravity, can the dark matter is due to gravitation energy, which act as mass?).

Argonne HELIOS spectrometer & its scientific success

Leave a comment

I have a chance to give a small seminar in Hong Kong University. Here I attached the ppt as pdf file.

argonne helios spectrometer &amp; its scientific success


The talk is an overview of the HELIOS spectrometer in Argonne National Laboratory ( which I am mainly working on). This talk also help me what is done and what is not explored. Hope you will find it interesting.

Older Entries