Examples of the use of 6-j symbols

Many of these sections rely on the previous sections, so you will need to step through, evaluating cells as you go for the later cells to work

Turn off the annoying and not helpful warnings

Off[General::spell1];
Off[General::spell];
Off[ClebschGordan::phy];
Off[ClebschGordan::tri];

Verify the definition of the 6-j symbols in terms of sums over Clebsch-Gordan coefficients

Mathematica doesn't deal well with algebraic forms of 6-j's and C-G coefficients, so take a couple of specific examples.

Example 1:

Take specific values of F2, F1, J2, J1, I1, and assume the operator is a vector (k = 1). (I called it I1, because otherwise mathematic thinks its i=√-1) (This example is appropriate, e.g., for the 5p3/2 -> 5s1/2 transition in 87Rb, which has I=3/2.)

F2 = 2;
F1 =  3;
J1 = 3/2;
J2  =  1/2;
I1  =  3/2;
k  =  1;

From class, we derived the form of  <F2 ||(T^k)|| F1>  in terms of  <J2 ||(T^k)|| J1>  in terms of a sum of a product of three C-G's.  The reduced matrix element <F2 ||(T^k)|| F1> should be independent of q, mF1, and mF2, so take q=mF1=mF2=0, when evaluating the sum over three C-G's from class:  <F2 ||(T^k)|| F1> /  <J2 ||(T^k)|| J1> =

Simplify[
(((Sqrt[2  F2 + 1])/(Sqrt[2  J2 + 1]))) 
        (1/ClebschGordan)[{F1, 0}, {k, 0}, {F2, 0}]
Sum[
ClebschGordan[{J2, mJ2}, {I1, mI2}, {F2, 
                0}]
ClebschGordan[{J1, mJ1}, {k, 0}, {J2, 
                mJ2}] 
ClebschGordan[{J1, mJ1}, {I1, 
                mI1}, {F1, 0}], 
{mJ1, -J1, 
            J1}, {mI1, -I1, I1}, {mJ2, -J2, J2}, {mI2, -I1, 
            I1}]]

This is also defined in terms of 6-j coefficients, and the two agree.

(-1)^(F1 + J2 + k + I1) Sqrt[2 F1 + 1] Sqrt[2 F2  + 1] 
      SixJSymbol[{J2, J1, k}, {F1, F2, I1}]

Note that in this particular definition (<F2 ||(T^k)|| F1> /  <J2 ||(T^k)|| J1>) = (<F1 ||(T^k)|| F2> /  <J1 ||(T^k)|| J2>)  (i.e. flipping F2 <->F1 and J2<->J1 gives the same answer)

F2 = 3;
F1 =  2;
J1 = 1/2;
J2  =  3/2;
(-1)^(F1 + J2 + k + I1) Sqrt[2 F1 + 1] Sqrt[2 F2  + 1] 
      SixJSymbol[{J2, J1, k}, {F1, F2, I1}]

Example 2:

Take specific values of F2, F1, J2, J1, I1, and assume the operator is a vector (k = 1).

F1 = 2;
F2  =  2;
J1 = 1/2;
J2  =  1/2;
I1  =  3/2;
k  =  1;

From class, we derived the form of  <F2 ||(T^k)|| F1>  in terms of  <J2 ||(T^k)|| J1>  in terms of a sum of a product of three C-G's.  The reduced matrix element <F2 ||(T^k)|| F1> should be independent of q, mF1, and mF2, but in this case, the selection rules don't allow mF, q = 0, so we have to take another value of F's. Take q=0, mF2 = 2, mF1 = 2.

Simplify[
(((Sqrt[2  F2 + 1])/(Sqrt[2  J2 + 1]))) 
        (1/ClebschGordan)[{F1, F1}, {k, 0}, {F2, F2}]
Sum[
ClebschGordan[{J2, mJ2}, {I1, mI2}, {F2, 
                F2}]
ClebschGordan[{J1, mJ1}, {k, 0}, {J2, 
                mJ2}] 
ClebschGordan[{J1, mJ1}, {I1, 
                mI1}, {F1, F1}], 
{mJ1, -J1, 
            J1}, {mI1, -I1, I1}, {mJ2, -J2, J2}, {mI2, -I1, 
            I1}]]

This is also defined in terms of 6-j coefficients, and the two agree.

(-1)^(F1 + J2 + k + I1) Sqrt[2 F1 + 1] Sqrt[2 F2  + 1] 
      SixJSymbol[{J2, J1, k}, {F1, F2, Ival}]

A Different definition, used in the rest of this Mathematica notebook

Note that a different definition of the reduced matrix element which doesn't have the √(2F2+1) and √(2J2 +1) would be as below:

F1 = 3;
F2  =  2;
J1 = 3/2;
J2  =  1/2;
I1  =  3/2;
k  =  1;

In this particularly definition for <F2 ||(T^k)|| F1>, we don't have the (Sqrt[2  F2 + 1]) and (Sqrt[2  J2  + 1]):

Simplify[
(1/ClebschGordan)[{F1, 0}, {k, 0}, {F2, 0}]
Sum[
ClebschGordan[{J2, mJ2}, {I1, 
                mI2}, {F2, 0}]
ClebschGordan[{J1, mJ1}, {k, 
                0}, {J2, mJ2}] 
ClebschGordan[{J1, 
                mJ1}, {I1, mI1}, {F1, 0}], 
{mJ1, -J1, 
            J1}, {mI1, -I1, I1}, {mJ2, -J2, J2}, {mI2, -I1, 
            I1}]]

This <F2 ||(T^k)|| F1> is also defined in terms of 6-j coefficients, and the two agree, as long as (Sqrt[2  F2 + 1]) is replaced with (Sqrt[2  J2  + 1]) )

(-1)^(F1 + J2 + k + I1) Sqrt[2 F1 + 1] Sqrt[2 J2  + 1] 
      SixJSymbol[{J2, J1, k}, {F1, F2, I1}]

Note that in this particular definition (<F2 ||(T^k)|| F1> /  <J2 ||(T^k)|| J1>) =  \(TraditionalForm\`\((\@\(\(2\ F2\)\(+\)\(1\)\(\ \)\)/\@\(\(2\ F1\)\(+\)\(1\)\(\ \)\))\)\ \((\ \@\(\(2\ J1\)\(+\)\(1\)\(\ \)\)/\@\(\(2\ J2\)\(+\)\(1\)\(\ \ \)\))\)\) (<F1 ||(T^k)|| F2> /  <J1 ||(T^k)|| J2>)  (i.e. flipping F2 <->F1 and J2<->J1 gives a different answer by the amount  \(TraditionalForm\`\((\@\(\(2\ F2\)\(+\)\(1\)\(\ \)\)/\@\(\(2\ F1\)\(+\)\(1\)\(\ \)\))\)\ \((\ \@\(\(2\ J1\)\(+\)\(1\)\(\ \)\)/\@\(\(2\ J2\)\(+\)\(1\)\(\ \ \)\))\)\) )

F1 = 2;
F2  =  
        3;
J1 = 
        1/2;
J2  =  
        3/2;
(-1)^(F1 + J2 + k + I1) 
      Sqrt[2 F1 + 1] Sqrt[2 J2  + 1] 
      SixJSymbol[{J2, J1, k}, {F1, F2, I1}]
((Sqrt[2 F2+1 ]) (Sqrt[2 J1+1 ])/((Sqrt[2 J2+1 ]) (Sqrt[2 F1+1 ])))(-1)^(F1 + J2 + k + I1)Sqrt[
        2 F1 + 1]Sqrt[
        2 J2  + 1]SixJSymbol[{J2, J1, k}, {F1, F2, I1}]

Application to photon scattering

Consider an atom with a single valence electron with a ground state in an s-orbital, e.g. the configuration is ns, where n is the principle quantum number. (This could be any of the alkalais Li, Na, K, Rb etc., or any ion of the alkaline earths, Be+ etc.)  Assume that the atom is well described by LS coupling. Consider second order scattering excitations through the  np configuration back to ns.  Assume the nucleus has angular momentum Ival=1/2. The three J-states to consider are (D1 and D2 are common labels for the two excied transitions)
| ground > = | ns 2S1/2>  
| excitedD1 > = | np 2P1/2>
| excited D2> = | np 2P3/2>
Including the hyperfine with Ival =1/2, there will be six states, with Fground = 0,1   FexcitedA = 0,1 FexcitedB = 2,1, for a total of six F states (each of which has 2F+1 mF states.) The resonant frequency of the p states is ~ω2p, but the two p states are split by spin orbit coupling which (fine structure) shifts the levels according to 
(ΔFS/2)(J(J+1)-S(S+1)-L(L+1)), 
Explicitly, this is  P3/2 is shifted up by ΔFS/2 and P1/2 down by -ΔFS.  (Total splitting is 3/2 ΔFS)
The hypefine interaction further splits the levels accordint to
(ΔHFS/2)(F(F+1)-J(J+1)-I(I+1)),

Set up and define the six states of the problem:

S  =  1/2;

    Inuc  =  1/2;

    Lgrnd  =  0;

    Jgrnd  =  1/2;

    Lp  =  1;

    JD1  =  1/2;

    JD2  =  3/2;

    SO[j1_, j2_, 
          j3_] :=  (1/2) (j3 (j3 + 1) - j2 (j2 + 1) - 
              j1 (j1 + 1)) ;

Label the states by {S,L,J,I,F, ω}, where ω is the resonant frequency, and group them according excited or ground:

GndStateF0  =  {S, Lgrnd, Jgrnd, Inuc, 
        0, ΔHFSgrnd SO[Jgrnd, Inuc, 
            0]}  
GndStateF1  =  {S, Lgrnd, Jgrnd, Inuc, 
        1, ΔHFSgrnd SO[Jgrnd, Inuc, 
            1]}
D1StateF0  = {S, Lp, JD1, Inuc, 
        0, ω2p + ΔFS SO[Lp, S, 
              JD1] + ΔHFSD1 SO[JD1, Inuc, 
              0]  }
D1StateF1  = {S, Lp, JD1, Inuc, 
        1, ω2p + ΔFS SO[Lp, S, 
              JD1] + ΔHFSD1 SO[JD1, Inuc, 
              1]  } 
D2StateF1  = {S, Lp, JD2, Inuc, 
        1, ω2p + ΔFS SO[Lp, S, 
              JD2] + ΔHFSD2 SO[JD2, Inuc, 
              1]  } 
D2StateF2  = {S, Lp, JD2, Inuc, 
        2, ω2p + ΔFS SO[Lp, S, 
              JD2] + ΔHFSD2 SO[JD2, Inuc, 
              2]  } 
GndStates  =  {GndStateF0, 
          GndStateF1};
D1States  =  {D1StateF0, D1StateF1};
D2States  =  {D2StateF1, D2StateF2};

Define some functions for extracting S, L, J, F, ω:

Sval[state_ ] := state[[1]];
Lval[state_ ] := state[[2]];
Jval[state_ ] := state[[3]];
Ival[state_ ] := state[[4]];
Fval[state_ ] := state[[5]];
ωval[state_]  :=  
        state[[6]];
HighIndex[
          state_]  :=  (Fval[state])^2 - (Ival[state] + 
                Jval[state])^2 + 1;
LowIndex[
          state_]  :=  (Fval[state])^2 - (Ival[state] - 
                Jval[state])^2 + 1;
MidIndex[
          state_]  :=  (Fval[state])^2 - (Ival[state] - 
                Jval[state])^2 + 1 + 
          Fval[state];
NumStates[States_] :=  
        HighIndex[
            States[[(Dimensions[States])[[1]]  ]]  ] -  
          LowIndex[States[[1]]  ];

First: Consider first order transition elements between ground and excited states

<gnd FmF | dq | F' mF' >  Note that these elements are block off-diagonal, since the dipole doesn't connect ground to ground or excited to excited. These elements can be written as a reduced matrix element  < J || d || J' > times a coefficient involving 6-j symbols Note: I'm using the √(2J'+2) √(2F+1) normalizaton:

Mcoef[A_, mA_, dλ_, B_, 
          mB_] :=  
(-1)^(Fval[A] + 1 + 
                Jval[B] + Ival[A])*
Sqrt[
            2 Jval[B] + 1]*
Sqrt[
            2  Fval[A] + 1]*
SixJSymbol[{Jval[B],  
              Jval[A],  1}, {Fval[A],  Fval[B],  
              Ival[A]}]*
ClebschGordan[{Fval[A], mA}, {1, 
              dλ}, {Fval[B], mB}];

Here, dλ is the q value for the ploarization of the light that is causing the transition, -1, 0, +1 for σ-, π, σ+ light., which is represented as a vector {1,0,0}, {0,1,0} and {0,0,1} for σ-, π, σ+ respectively.

Make a table to look at the relative strengths of  the first order transitions elements for all polarizations

Ignore the complicated syntax, this just makes a table of all the possible F,mF  to F'mF' via dλ combinations.

DMatrix[GndHFStates_, ExHFStates_, d_] := 
Flatten[
          Table[
Table[
Flatten[
                Table[
Table[
                    Mcoef[GndHFStates[[j]], mi, d, ExHFStates[[n]], 
                      mj], {mi, -Fval[GndHFStates[[j]]], 
                      Fval[GndHFStates[[j]]] } 
], {j, 1, Dimensions[
                        GndHFStates][[1]]}
]], {mj, -Fval[ExHFStates[[n]] ], 
                Fval[ExHFStates[[n]] ]} 
],   {n, 
              1, Dimensions[ExHFStates][[1]]}
], 
          1];
D1Matrix  = 
        Transpose[ 
          Table[DMatrix[GndStates, D1States, d], {d, -1, 1}], {3, 1, 
            2}];
D2Matrix  =  
        Transpose[
          Table[DMatrix[GndStates, D2States, d], {d, -1, 1}], {3, 1, 
            2}];

The following matrices are the block off-diagonal part of the full coupling Hamiltonian, and should be interpreted as follows: The D1-> ground is 4x4:  4 excited states and 4 ground states. F', mF' = |0,0>  |1,-1>   |1,0>  |1,1>, similarly for the ground state. The D2 -> ground is 4x8, with 8 excited states (2F+1 for F=1 and 2F+1 for F=2, running from |1,-1>  up to |2,2>). (see notes). For each state coupling, the vector inside tells which polarization couples.  For example, in the D1 transition, the excited |1,-1>  state is coupled to the ground |0,0> state by 1/3.  The fully stretched transition Fgnd=1 mFgnd=1 to Fexc=2 mFexc=2 has coupling strength 1. (which is why i used this normalization convention)

MatrixForm[D1Matrix^2]
MatrixForm[D2Matrix^2]

Now consider second order transitions from ground to ground states (but only considering a single set of levels J excited (D1 or D2 only)  for now).

This is a sum over   <gnd FmF | d | i >< i | d |gnd F' mF' >/(ω-ωi), and represents the scattering rate for incoming polarization pol =   For the two level problem, we were able to average over scattered polarizations to get |d.Esc|^2 = (1/3) |d|^2 .  Here we can do the same thing, but  the vector matrix element (called MD1 here) has to be summed over all intermediate states first. (There was only one intermediate state in the 2 level problem.) For now I have left out the reduced matrix elements, which provide an overal constant.

MD1[ω_, pol_, Astate_, mFA_, Bstate_, 
          mFB_]  :=  (∑_(n = 1)^(Dimensions[D1States][[1]]))(∑_(j = -Fval[D1States[[n]]])^(Fval[D1States[[n]]]))((1/(ωval[D1States[[n]]] - ωval[
                        Astate] - ω)) (D1Matrix[[MidIndex[
                            D1States[[n]]] + j]][[MidIndex[
                          Bstate] + 
                        mFB]]) pol . (D1Matrix[[MidIndex[
                              D1States[[n]]] + j]][[MidIndex[
                            Astate] + 
                          mFA]]) );
ScattD1[ω_, pol_, Astate_, mFA_, Bstate_, mFB_]  :=  
        MD1[ω, pol, Astate, mFA, Bstate, mFB] . 
            MD1[ω, pol, Astate, mFA, Bstate, mFB]/3;

MD2[ω_, pol_, Astate_, mFA_, Bstate_, 
          mFB_]  := (∑_(n = 1)^(Dimensions[D2States][[1]]))(∑_(j = -Fval[D2States[[n]]])^(Fval[D2States[[n]]]))((1/(ωval[D2States[[n]]] - ωval[
                        Astate] - ω)) (D2Matrix[[MidIndex[
                            D2States[[n]]] + j]][[MidIndex[
                          Bstate] + 
                        mFB]]) pol . (D2Matrix[[MidIndex[
                              D2States[[n]]] + j]][[MidIndex[
                            Astate] + 
                          mFA]])  );
ScattD2[ω_, pol_, Astate_, mFA_, Bstate_, mFB_]  :=  
        MD2[ω, pol, Astate, mFA, Bstate, mFB] . 
            MD2[ω, pol, Astate, mFA, Bstate, mFB]/3;

The Hyperfine splittings are typically very small compared to the transitions frequencyies. (Since we are considering only individual fine structure levels, the splitting is not too important, it just shifts the two levels.) Plot this for different input states, with different polarizations, using example values of 
ω2p = 10^15, 
ΔFS = 10^13
ΔHFS = 10^10 for grnd, 10^9 for excited.
In particular, compare the Δm = 2 scattering rates to the Δm=1, or Δm=0 rates: Notice in the zoomed in version, one can easily get reasonable Δm =2, but that far from resonance, the Δm=2 rate can be orders of magnitude smaller than the Δm=1,0 rates, which stay on the same order of magnitude throughout. (Also note the interesting interference effects giving zero Δm=1 rates at some point, but the total rate is never zero). (Also note that for frquencies ω-> 0, we are probably not justified in making the rotating wave approximation....but the general idea is right.)

vals  =  {ω2p ->  10^15, ΔFS  ->  
            10^13,  ΔHFSgrnd  ->  
            10^10, ΔHFSD2 ->  
            10^9};
d = {1, 0, 0};      

Astate  =  GndStates[[2]];      
mFA  =  
      1;
Bstate  =  GndStates[[2]];     
mFB  =  -1;
scattΔM2  =  
        ScattD2[ω, d, Astate, mFA, Bstate, mFB] /.  
          vals;

mFB  =  0;    
scattΔM1  =  
        ScattD2[ω, d, Astate, mFA, Bstate, mFB] /.  
          vals;

mFB  =  1;   
scattΔM0  =  
        ScattD2[ω, d, Astate, mFA, Bstate, mFB] /. 
          vals;

scatttotal  =  
        scattΔM2 + scattΔM1 + 
          scattΔM0;

Plot[{Log[10, scattΔM2], Log[10, scattΔM1], 
        Log[10, scattΔM0], Log[10, scatttotal]}, {ω,  
        0 10^15, 1.6 10^15}, 
      PlotStyle ->  {RGBColor[1, 0, 0], RGBColor[0, 0, 1], 
          RGBColor[0, 1, 0], RGBColor[0, 1, 1]}, 
      AxesLabel ->  "Wide frequency range, red = \
ΔM=2"]
Plot[{Log[10, scattΔM2], Log[10, scattΔM1], 
        Log[10, scattΔM0], Log[10, scatttotal]}, {ω,  
        1.00499 10^15, 1.00501 10^15}, 
      PlotStyle ->  {RGBColor[1, 0, 0], RGBColor[0, 0, 1], 
          RGBColor[0, 1, 0], RGBColor[0, 1, 1]}, 
      AxesLabel ->  "Zoomed in, the two peaks are the F=2, F=1 \
excited states"]

Where does the suppression come from? Analytically, the rates look like: (the fine structure offset is the same, and the ground state hyperfine is the same, so, let δ =  ω - (ω 2p+ΔFS/2-ΔHFSgrnd/4),  )  Notice that the Δm=2 has a minus sign between the two terms, which for causes the terms to cancel if ΔHFSD2=0. This is the result of a sum rule on the 6-j coefficients.

d = {1, 0, 0};       
Astate  =  
      GndStates[[2]];      
mFA  =  
      1;
Bstate  =  GndStates[[2]];    
mFB  =  -1;
scattΔM2  = 
      ScattD2[δ  + ω2p + ΔFS/
            2 - ΔHFSgrnd/4 , d, Astate, mFA, Bstate, 
        mFB]
mFB  =  0;    
scattΔM1  = 
      ScattD2[δ  + ω2p + ΔFS/
            2  - ΔHFSgrnd/4, d, Astate, mFA, Bstate, 
        mFB]
mFB  =  1;    
scattΔM0  = 
      ScattD2[δ  + ω2p + ΔFS/
            2  - ΔHFSgrnd/4, d, Astate, mFA, Bstate, 
        mFB]

Doing a Taylor series expansion reveals that the Δm=2 has a different power law, namely (1/δ^2)(ΔHFSD2/δ)^2 ∝ 1/δ^4:

Series[
      scattΔM2, {ΔHFSD2, 0, 
        2}]
Series[
      scattΔM1, {ΔHFSD2, 0, 
        2}]
Series[
      scattΔM0, {ΔHFSD2, 0, 2}]

Combining the two excited J levels,  (D1 and D2) transitions:

Unfortunately I didn't have time to do this, but the general idea is:  express < J || d || J' > in terms of < L || d || L'> in exactly the same way that < F || d || F' > was expressed in terms of < J || d || J' >,   then when everything is added up, detuning far from the fine structure will make the Δm =1 transitions get really small.