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

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

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

using namespace std;

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

Short-range interaction on same nucleons pair

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

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 .
3. Most parallel states are either

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 then or via versa. If is formed by pure S = 1 state, thus, the most affected state is as 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, ,

The J = 3 state is unaffected by δ-interaction

Sum rule of 9-j symbol

The 9-j symbol is the coupling coefficient when combining 2 nucleons with state and  to from the state and

or in simpler form

or reverse.

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

we can also see that

---

For example, when and coupled to   . We have

Since for each , it is degenerated. Thus, the sum rule is

The sum is equal the a fraction, because the left-hand side is degenerated for the state .

Winger 6-j and 9-j symbol

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, .

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

the other way is

The 6-j symbol is

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

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

In Mathematica, there is a build in function

---

The 9-j symbol is the coupling between 4 angular momenta, .

The 9-j symbol can be used in coupling 2 nucleons, .

The 9-j symbol is

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

Where 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 and to from a state is