## Spectroscopic factor (revisit)

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$

Thus,

$\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)$

$\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$

recall

$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$

## 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

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

## 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, $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:

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

## Evaluation of mutual interaction

We are going to evaluation the integral

$\displaystyle G_{ij}^{hk} = \langle b_i(x) b_h(y) | \frac{1}{r_{xy}} | b_j(x) b_k(y) \rangle$

Recalling the multi-pole expansion,

$\displaystyle \frac{1}{r_{12}} = \sum_{l=0}^{\infty} \frac{4 \pi}{2l+1} \frac{r_{<}^l}{r_{>}^{l+1}} \sum_{m=-l}^{l} Y_{lm}^{*}(\Omega_1) Y_{lm}(\Omega_2)$

and the basis

$b_{nlm}(\vec{r}) = R_{nl}(r) Y_{lm}(\Omega)$

Set an averaged basis

$\displaystyle b_{nl}(\vec{r}) = R_{nl}(r) \frac{1}{\sqrt{2l+1}} \sum_{m=-1}^{l}Y_{lm}(\Omega)$

$\displaystyle \Gamma_{ij}^{hk}(l) = \frac{4 \pi}{2l+1} \sum_{m=-l}^{l} \frac{1}{\sqrt{(2l_i+1)(2l_j+1)(2l_h+1)(2l_k+1)}} \\ \sum_{m_i,m_j}\int Y_{l_i m_i}^{*}(\Omega_1) Y_{l m}^{*}(\Omega_1) Y_{l_j m_j}(\Omega_1) d \Omega_1 \\ \sum_{m_h,m_k} \int Y_{l_h m_k}^{*}(\Omega_2) Y_{l m}(\Omega_2) Y_{l_k m_k}(\Omega_2) d \Omega_2$

In the angular integrals, using Wigner 3-j symbol and the integral

$\displaystyle\int Y_{l_1 m_1}(\Omega) Y_{l m}(\Omega) Y_{l_2 m_2}(\Omega) d \Omega \\= \frac{\sqrt{(2l_1+1)(2l+1)(2l_2+1)}}{\sqrt{4\pi}} \begin{pmatrix} l_1 & l & l_2 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l & l_2 \\ m_1 & m & m_2 \end{pmatrix}$

$Y_{lm}^{*}(\Omega) = (-1)^m Y_{l(-m)}(\Omega)$

$\displaystyle \begin{pmatrix} l_1 & l & l_2 \\ m_1 & m & m_2 \end{pmatrix} = (-1)^{l_1+l+l_2} \begin{pmatrix} l_1 & l & l_2 \\ -m_1 & -m & -m_2 \end{pmatrix} \\ =(-1)^{l_1+l+l_2} \begin{pmatrix} l_1 & l_2 & l \\ m_1 & m_2 & m \end{pmatrix}$

we have

$\displaystyle\int Y_{l_i m_i}^{*}(\Omega) Y_{l m}^{*}(\Omega) Y_{l_j m_j}(\Omega) d \Omega \\= (-1)^{m_j} \frac{\sqrt{(2l_i+1)(2l+1)(2l_j+1)}}{\sqrt{4\pi}} \begin{pmatrix} l_i & l & l_j \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_i & l & l_j \\ m_i & m & -m_j \end{pmatrix}$

$\displaystyle\int Y_{l_h m_h}^{*}(\Omega) Y_{l m}(\Omega) Y_{l_k m_k}(\Omega) d \Omega \\= (-1)^{m_h} \frac{\sqrt{(2l_h+1)(2l+1)(2l_k+1)}}{\sqrt{4\pi}} \begin{pmatrix} l_h & l & l_k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_h & l & l_k \\ -m_h & m & m_k \end{pmatrix}$

Thus,

$\displaystyle \Gamma_{ij}^{hk}(l) = \sum_{m=-l}^{l} \sum_{m_i,m_j } (-1)^{m_j} \begin{pmatrix} l_i & l_j & l \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_i & l_j & l \\ m_i & -m_j & m \end{pmatrix} \\ \sum_{m_h,m_k} (-1)^{m_h} \begin{pmatrix} l_h & l_k & l \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_h & l_k & l \\ -m_h & m_k & m \end{pmatrix}$

The 3-j symbol restricted that

$l_i+l_j+l = even$

$l_h+l_k+l = even$

$m_i-m_j + m = 0, |m_i| \leq l_i , |m_j| \leq l_j$

$-m_h + m_k +m = 0, |m_h| \leq l_h , |m_k| \leq l_k$

I guess it is the most simplified formula for the angular part.

The total integral is

$\displaystyle G_{ij}^{hk} = \sum_{l=0}^\infty \Gamma_{ij}^{hk}(l) \langle R_i(x) R_h(y) | \frac{r_<^l}{r_>^{l+1}} |R_j(x) R_k(y) \rangle$

The angular integral imposes condition for $l$.

I am not sure this is a correct way to treat the problem.

First, the averaged basis is still an energy eigen state. It is not the eigen state for the angular part. So, this averaging could introduces an error and we should reminded that this is an approximation. But in the perturbation view point, this averaged basis is still valid.

Second thing is, the sum

$\displaystyle P_{(l_i m_i) (l_j m_j)}^{(l_h m_h)(l_k m_k)}(l) = \sum_{m=-l}^{l} \langle Y_{l_i m_i}|Y_{lm}^*|Y_{l_j m_j}\rangle \langle Y_{l_h m_h}|Y_{lm}|Y_{l_k m_k}\rangle$

is not symmetry for exchange of $i, j, h,k$ in general. For example,

$\displaystyle \frac{-1}{3}=P_{(1-1) (00)}^{(11)(00)}(1) {\neq} P_{(00)(1-1)}^{(11)(00)}(1) = 0$

This is a very uprising result that the mutual interaction dependent on the magnetic quantum number. Thus, in detail, we should use $n l m m_s$ as a basis.

Third, the sum $P_{ij}^{hk}(l)$ is depend on $l$. The mutual interaction require us to sum all possible $l$.

Fourth, the coupling between 1s2p triplet state, the total spin is $S = 1$, total L is $L = 1$, and the total angular momentum can be  $J = 0, 1, 2$. In our treatment, we did not coupled the angular momentum in the calculation explicitly. In fact, in the integral of the spherical harmonic, the coordinate are integrated separately, and the coupling seem to be calculated implicitly. I am not sure how to couple two spherical harmonics with two coordinates.

## Integration formulas of spherical harmonic

There are several important and useful formulas for the integration of spherical harmonic. We simplify the notation,

$\displaystyle \int_0^{\pi} \sin(\theta)d\theta\int_0^{2\pi}d\phi = \int d\Omega$

The first one is the average of spherical harmonic.

$\displaystyle \int Y_{lm} d\Omega = \sqrt{4\pi} \delta_{l0}\delta_{m0}$

The 2nd one is the orthonormal  condition.

$\displaystyle \int Y^{*}_{l'm'}Y_{lm} d\Omega = \delta_{l'l}\delta_{m'm}$

The 3rd one is triplet integral, we use the product of spherical harmonic,

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y^*_{l_3m_3} d\Omega \\ = \int \sum_{lm} \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l+1)}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} Y_{lm} Y^*_{l_3m_3} d\Omega \\= \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l_3+1)}} C_{l_10l_20}^{l_30} C_{l_1m_1l_2m_2}^{l_3m_3}$

The 4th one is another triple integral,

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y_{l_3m_3} d\Omega \\ = \int \sum_{lm} \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l+1)}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} Y_{lm} Y_{l_3m_3} d\Omega \\ = \int \sum_{lmLM} \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l+1)}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} \\ \sqrt{\frac{(2l+1)(2l_3+1)}{4\pi(2L+1)}} C_{l0l_30}^{L0} C_{lml_3m_3}^{LM}Y_{LM}d\Omega$

$\displaystyle = \sum_{lm} \sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} C_{l0l_30}^{00} C_{lml_3m_3}^{00}$

Notice that

$C_{lmLM}^{00} = (-1)^{L+M} \sqrt{\frac{1}{2L+1}} \delta_{Ll}\delta_{-m,M}$

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y_{l_3m_3} d\Omega = \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi (2l_3+1)}} C_{l_10l_20}^{l_30} C_{l_1m_1l_2m_2}^{l_3,-m_3} (-1)^{m_3}$

using Wigner 3-j symbol,

$C_{l_1m_1l_2m_2}^{l_3m_3} = (-1)^{l_1-l_2+m_3} \sqrt{2l_3+1} \begin{pmatrix} l_1 & l_2 & l_3 \\ m_1 & m_2 & -m_3 \end{pmatrix}$

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y_{l_3m_3} d\Omega \\= \sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi}} \begin{pmatrix} l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{pmatrix}$

For other integral, we can use

$Y^*_{lm}(\theta, \phi) = (-1)^{m}Y_{l(-m)}(\theta,\phi) = Y_{lm}(\theta, -\phi)$