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

Example: L=1,S=1/2, I = 0 Zeeman structure. (e.g. He+, 2p states)

L = 1;
S = 1/2;

Spin-Orbit In the J=3/2 (4 states) and J=1/2 (2 states) mJ basis, the spin orbit is diagonal, and has the form

a = (β/2) (J (J + 1) - L (L + 1) - 
              S (S + 1))  /.  
        J -> 1/2
b = (β/2) (J (J + 1) - L (L + 1) - 
              S (S + 1))  /.  
        J -> 3/2
Hso[β_]  = 
        DiagonalMatrix[{b, b, b, b, a, a}];
MatrixForm[Hso[β]]

B-field:  The magnetic field can couple states with the same mJ, but cannot couple states that differ in mJ.  That leads to diagonal terms, and the only off diagonal terms are for the   3/2,1/2 <--> 1/2,1/2 and  3/2,-1/2 <--> 1/2,-1/2.  The coupling elements are a product of Clebsch-Gordans with (mL+2mS).  Another way to view the diagonal terms is that they should just be gJ  mJ.

Couple[J1_, J2_, mJ_]  :=  
        Sum[ClebschGordan[{L, mL}, {S, mS}, {J1, mJ}] 
            ClebschGordan[{L, mL}, {S, mS}, {J2, mJ}] (mL + 
                2  mS), {mL, -L, L}, {mS, -S, 
            S}];
g32  =  3/
          2 + (S (S + 1) - 
              L (L + 1))/(2 3/2 (3/2 + 1))    
g12  =  
      3/2  +  (S (S + 1) - 
              L (L + 1))/(2 1/2 (1/2 + 1))

Couple[3/2, 3/2, 
      3/2]     
g32  3/2    
Couple[3/2, 3/2, 
      1/2]      
g32   1/2    
Couple[1/2, 1/2, 1/2]     
g12   1/2    
Couple[1/2, 3/2, 
      1/2]

HZeeman[B_] =( ({{2, 0, 0, 0, 0, 0}, {0, 2/3, 0, 0, -(Sqrt[2])/3, 0}, {0, 0, -2/3, 0, 0, -(Sqrt[2])/3}, {0, 0, 0, -2, 0, 0}, {0, -(Sqrt[2])/3, 0, 0, 1/3, 0}, {0, 0, -(Sqrt[2])/3, 0, 0, -1/3}}))B;

Diagonalize, and take a look:

A = Eigensystem[Hso[β] + HZeeman[B]];

Plot[Evaluate[A[[1]] /. β ->  1], {B, 0, 8}]   
Plot[
      Evaluate[D[A[[1]] /. β ->  1, B]], {B, 0, 3}]

MatrixForm[
      Normal[Series[
          A[[1]] /. β ->  1, {B, 0, 2}]]]

Example: L=0,S=1/2, J=1/2, I = 3/2 Hyperfine Zeeman structure. (e.g. 87Rb, ground state)
The Breit-Rabi formula

J = 1/2;
Ival = 3/2;

Hyperfine interaction: In the F=2 (5 states) and F=1 (3 states) mF basis, the hyperfine interactions (Ahf I.J) is diagonal, and has the form (for J=1/2, there is no quadrupole field gradient, and therefore no quadrupole term).

a = (Ahf/2) (F (F + 1) - J (J + 1) - 
              Ival (Ival + 1))  /.  
        F -> 1
b = (Ahf/2) (F (F + 1) - J (J + 1) - 
              Ival (Ival + 1))  /.  
        F -> 2
Hhf[Ahf_]  = DiagonalMatrix[{b, b, b, b, b, a, a, a}];

B-field:  The magnetic field can couple states with the same mF, but cannot couple states that differ in mF.  That leads to diagonal terms, and the  off-diagonal terms are for the   2,1 <--> 1,1 and  2,0 <--> 1,0 and 2,-1 <--> 1,-1.  The coupling elements are a product of Clebsch-Gordans with mF (since μN ~0).  Another way to view the diagonal terms is that they should just be gF  mF.  Note that F=2 and F=1 have opposite g-factors: the same atomic angular momentum gives rise to opposite Zeeman shifts.

Couple[F1_, F2_, mF_]  :=  
        Sum[ClebschGordan[{J, mJ}, {Ival, mI}, {F1, mF}] 
            ClebschGordan[{J, mJ}, {Ival, mI}, {F2, mF}] mJ, {mJ, -J, 
            J}, {mI, -Ival, Ival}];

    g2  =  1/
          2 + (J (J + 1) - 
              Ival (Ival + 1))/(2   2 (2 + 
                  1))    

    g1  =  1/
          2 + (J (J + 1) - 
              Ival (Ival + 1))/(2  1 (1 + 1))    

    gJ  =  2;

Couple[2, 2, 2] gJ     
g2 gJ  2   
Couple[2, 2, 1]  gJ    
g2 gJ   1    
Couple[2, 2, 0]  gJ    
g2 gJ   0    
Couple[1, 1, 1]  gJ   
g1 gJ  1    
Couple[1, 2, 
        1]  gJ    
Couple[1, 2, 
        0]  gJ    
Couple[1, 
        2, -1]  gJ

HZeemanhf[B_] =( ({{1, 0, 0, 0, 0, 0, 0, 0}, {0, 1/2, 0, 0, 0, (Sqrt[3])/2, 0, 0}, {\[Placeholder], 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, -1/2, 0, 0, 0, (Sqrt[3])/2}, {0, 0, 0, 0, -1, 0, 0, 0}, {0, (Sqrt[3])/2, 0, 0, 0, -1/2, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, (Sqrt[3])/2, 0, 0, 0, 1/2}}))B;

Note that each state couples with at most one other state.  These can be diagonalized easily,  individually as 2x2 matrices. Let Mathematica do it: Two states (stretched states) are uncoupled to anything.

A = 
        Eigensystem[Hhf[Ahf] + HZeemanhf[B]];
MatrixForm[A[[1]]]

Plot[Evaluate[A[[1]] /. Ahf ->  1], {B, 0, 3}, 
      PlotStyle ->  {Thickness[0.004]}]   
Plot[
      Evaluate[D[A[[1]] /. Ahf ->  1, B], {B, 0, 3}, 
        PlotStyle ->  {Thickness[0.004]}]]

Example: L=1,S=1/2, J=3/2, I = 3/2 Hyperfine Zeeman structure. (e.g. 87Rb, excited 2p3/2 state)

J = 3/2;
Ival = 3/2;
S  =  1/2;
L = 1;

Hyperfine interaction: In the F=3 (7 states) F=2 (5 states) F=1 (3 states) and F=0 (1 state) mF basis, the hyperfine interactions (Ahf I.J) is diagonal, and has the form (for J=3/2, there is a quadrupole field gradient, but for simplicity here I am ignoring it).

a = (Ahf/2) (F (F + 1) - J (J + 1) - 
              Ival (Ival + 1))  /.  
        F -> 0
b = (Ahf/2) (F (F + 1) - J (J + 1) - 
              Ival (Ival + 1))  /.  
        F -> 1
c = (Ahf/2) (F (F + 1) - J (J + 1) - 
              Ival (Ival + 1))  /.  
        F -> 2
d = (Ahf/2) (F (F + 1) - J (J + 1) - 
              Ival (Ival + 1))  /.  
        F -> 3
Hhf[Ahf_]  = 
        DiagonalMatrix[
Join[
Table[
              d, {j, 1, 2 3 + 1}], 
Table[
              c, {j, 1, 2 2 + 1}], 
Table[
              b, {j, 1, 2 1 + 1}], 
{a}]];

B-field:  The magnetic field can couple states with the same mF, but cannot couple states that differ in mF.  That leads to diagonal terms, and multiple  off-diagonal terms are coupled.  The coupling elements are a product of Clebsch-Gordans with mJ (since μN ~0).  Another way to view the diagonal terms is that they should just be gF  mF. Since J = I, the gF factor is particularly simple: 1/2.

Couple[F1_, F2_, mF_]  :=  
        Sum[ClebschGordan[{J, mJ}, {Ival, mI}, {F1, mF}] 
            ClebschGordan[{J, mJ}, {Ival, mI}, {F2, mF}] mJ, {mJ, -J, 
            J}, {mI, -Ival, Ival}];
gF  =  1/2     
gJ  =  3/
          2 + (S (S + 1) - 
              L (L + 1))/(2 3/2 (3/2 + 1))

Note that more things couple than before: e.g. the F=2, mF = 0 states couple to F=1 and F=3.  Why doesn't F=2 couple to F=0?  (the SparseArray command is merely putting all the correct off-diagonal elements in place)

Nsize  =  
        2 3 + 1  + 2 2 + 1  + 2 1  + 1  + 
          1;
HZeemanhf[
          B_]  = (DiagonalMatrix[
Join[
Table[
                    gF gJ j, {j,  
                      3, - 3, -1}], 
Table[
                    gF gJ j, {j,  
                      2, - 2, -1}], 
Table[
                    gF gJ j, {j, 
                      1, -1, -1}], 
{0}]
] + 
SparseArray[
{{2, 
                      2 3 + 2} ->  
                    Couple[3, 2, 2]  gJ , {2 3 + 2, 2} -> 
                    Couple[3, 2, 2]  gJ  , 
{3, 
                      2 3 + 3} ->  
                    Couple[3, 2, 1]  gJ , {2 3 + 3, 3} -> 
                    Couple[3, 2, 1]  gJ, {4, 2 3 + 4} ->  
                    Couple[3, 2, 0]  gJ , {2 3 + 4, 4} -> 
                    Couple[3, 2, 0]  gJ, 
{5, 
                      2 3 + 5} ->  
                    Couple[3, 2, -1]  gJ , {2 3 + 5, 5} -> 
                    Couple[3, 2, -1]  gJ , 
{6, 
                      2 3 + 6} ->  
                    Couple[3, 2, -2]  gJ , {2 3 + 6, 6} -> 
                    Couple[3, 2, -2]  gJ , 
{2 3 + 
                        3, 2 3 + 1 + 2 2 + 2} ->  
                    Couple[2, 1, 1]  gJ , {2 3 + 1 + 2 2 + 2, 
                      2 3 + 3} -> 
                    Couple[2, 1, 1]  gJ, 
{2 3 + 4, 
                      2 3 + 1 + 2 2 + 3} ->  
                    Couple[2, 1, 0]  gJ , {2 3 + 1 + 2 2 + 3, 
                      2 3 + 4} -> 
                    Couple[2, 1, 0]  gJ, 
{2 3 + 5, 
                      2 3 + 1 + 2 2 + 4} ->  
                    Couple[2, 1, -1]  gJ , {2 3 + 1 + 2 2 + 4, 
                      2 3 + 5} -> 
                    Couple[2, 1, -1]  gJ, 
{Nsize - 
                        2, Nsize} ->  
                    Couple[1, 0, 0]  gJ , {Nsize, Nsize - 2} ->  
                    Couple[1, 0, 0]  gJ }
, {Nsize, 
                  Nsize}]) B;

MatrixForm[HZeemanhf[B]]

Diagonalize, and take a look:

A = Eigensystem[Hhf[1] + HZeemanhf[B]];

Plot[Evaluate[A[[1]]], {B, 0, 2}, 
      PlotStyle ->  {Thickness[0.004]}]