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;
}
Advertisements

Short-range interaction on same nucleons pair

Leave a comment

From “Nuclear structure from a simple perspective” by Richard F. Casten, Chapter 4. The first half discusses the coupling between pp or nn T=1 isovector pair under δ-interaction, which is a good approximation to short-range interaction.

I found that the book makes it a bit complicated. The only 3 cases matter

Capture.PNG

The logic is follow. As the spatial part of the total wavefunction must be symmetric for the δ-interaction to be effective. The spin-isospin part must be antisymmetric as the total wavefunction must be antisymmetric due to identical fermions system. Also, the isospin part must be isovector or symmetric. Thus, the spin part must be antisymmetric or S = 0.

  1. In other word, the δ-interaction has no effect on the S =1 state.
  2. Also, the book said that l_1 + l_2 - J = even.
  3. Most parallel states are either J = J_{min} , J_{max}

Using these 3 reasons, and consider these 3 cases, which actually covered all possible combinations and all cases. Point 2 also means that only half of the coupled states are affected. Since the number of state are even, when J_{max} = even then J_{min} = odd or via versa. If J_{max} is formed by pure S = 1 state, thus, the most affected state is J_{min} as J_{min} is the other most parallel states, in which the overlap of the wavafunctions is maximum.

In the book, the example is d5/2 f7/2, J = 1, 2, 3, 4, 5, 6. This is belong to case 1, thus, the lowest state is J = 1. another example is d5/2 g7/2, which is case 2, the lowest state is J = 6.

For equivalent orbit, either case 1 or case 3. The lowest state is J = 0.


There are one interesting case. when coupling the s1/2. Since the s-orbit is isotropic, thus, no matter how the other orbit, the overlap is always maximum. Since pure S = 1 state is unaffected, thus, we can somehow, imagine the nucleon pair is forming S = 0 pair. In fact, for example, when s1/2 couples with d5/2, using 9-j symbol, L = l_1+l_2, S = s_1+s_2 ,

|J = 2 \rangle = \frac{1}{10} |L S = 20\rangle + \frac{1}{15\sqrt{2}} |21\rangle

|J = 3 \rangle = \frac{1}{6\sqrt{5}} |21\rangle

The J = 3 state is unaffected by δ-interaction

Sum rule of 9-j symbol

Leave a comment

The 9-j symbol is the coupling coefficient when combining 2 nucleons with state |l_1 s_1 j_1 \rangle and |l_2 s_2 j_2 \rangle to from the state |j_1 j_2 J M \rangle and |L S J M \rangle

|L S J M \rangle = \sum_{j_1 j_2 } |j_1  j_2 J M \rangle \begin{pmatrix} l_1 & s_1 & j_1 \\ l_2 & s_2 & j_2 \\ L & S & J \end{pmatrix}

or in simpler form

|L S J M \rangle = \sum_{j_1 j_2 } |j_1  j_2 J M \rangle \langle j_1 j_2 | L S \rangle

or reverse.

|j_1 j_2 J M \rangle = \sum_{L S } |L  S J M \rangle \langle L S | j_1 j_2 \rangle

or to say, this is the coefficient between LS coupling and jj coupling scheme.

we can also see that

\langle L S | j_1 j_2 \rangle = \begin{pmatrix} l_1 & l_2 & L \\ s_1 & s_2 & S \\ j_1 & j_2 & J \end{pmatrix}


For example, when p_{1/2} and p_{3/2} coupled to J = 1  . We have

\displaystyle|p_{1/2} p_{3/2} (J = 1) \rangle = \frac{1}{9} |01\rangle -  \frac{1}{6\sqrt{6}} |10\rangle + \frac{1}{12} |11\rangle + \frac{1}{36} |21\rangle

Since for each |LS\rangle, it is (2L+1) (2S+1) degenerated. Thus, the sum rule is

\displaystyle \sum_{LS} \left(\langle L S | j_1 j_2 \rangle\right)^2 (2L+1) (2S+1) = \frac{1}{(2j_1+1)(2j_2+1)}

The sum is equal the a fraction, because the left-hand side is (2j_1+1) (2j_2+1) degenerated for the state |p_{1/2} p_{3/2} (J = 1) \rangle .

Winger 6-j and 9-j symbol

Leave a comment

The meaning of 3-j symbol is same as Clebsch-Gordan coefficient. So, we skip in here.

I am not going to construct the 6-j symbol from 3-j symbol. In here, I just state the meaning and usage in Mathematica.


The 6-j symbol is the coupling between 3 angular momenta, j_1, j_2, j_3 .

There are 2 ways to couple these 3 angular momenta. First,

j_1 + j_2 + j_3 \rightarrow j_{12} + j_3 \rightarrow J

the other way is

j_ 1 + j_2 + j_3 \rightarrow j_1 + j_{23}  \rightarrow J

The 6-j symbol is

\begin{pmatrix} j_1 & j_2 & j_{12} \\ j_3 & J & j_{23} \end{pmatrix}

We can see that there are 4 vector-sum must satisfy.

\Delta(j_1, j_2, j_{12})

\Delta(j_2, j_3, j_{23})

\Delta(j_1, j_{23}, J)

\Delta(j_{12}, j_3, J)

If we draw a line to connect these 4 vector-sum, we have:

Screen Shot 2018-05-16 at 15.34.47.png

In Mathematica, there is a build in function

\textrm{SixJSymbol}[ \left\{j_1, j_2, j_{23} \right\}, \left\{j_3, J , j_{23} \right\}]


The 9-j symbol is the coupling between 4 angular momenta, j_1, j_2, j_3, j_4 .

The 9-j symbol can be used in coupling 2 nucleons, l_1, s_1, l_2, s_2 .

The 9-j symbol is

\begin{pmatrix} l_1 & s_1 & j_1 \\ l_2 & s_2 & j_2  \\  L & S & J \end{pmatrix}

We can see, each row and column must satisfy the vector-sum.

Unfortunately, there is no build in function in Mathematica. The formula for 9-j symbol is

\displaystyle \begin{pmatrix} l_1 & s_1 & j_1 \\ l_2 & s_2 & j_2  \\  L & S & J \end{pmatrix} \\ = \sum_{g} (-1)^{2g} (2g+1) \begin{pmatrix} l_1 & s_1 & j_1 \\ j_2 & J & g \end{pmatrix} \begin{pmatrix} j_2 & s_2 & j_2 \\ s_1 & g & S \end{pmatrix} \begin{pmatrix} L & S & J \\ g & l_1 & l_2 \end{pmatrix}

Where g sum all possible value, which can be calculate using the 6 couplings inside the 3 6-j symbols.To check your result, the coupling between d_{5/2} and f_{7/2} to from a L = 5, S = 0, J = 0 state is

\begin{pmatrix} 2 & 1/2 & 5/2 \\ 3 & 1/2 & 7/2  \\  5 & 0 & 5 \end{pmatrix} = \frac{1}{2\sqrt{770}}