Operators

Sx ** K[l_, m_] ^= (Sp ** K[l, m] + Sm ** K[l, m])/
          2 ;
Sy ** 
          K[l_, m_] ^= (Sp ** K[l, m] - 
              Sm ** K[l, m])/(2 I);
Sp ** K[l_, m_] ^= 
        Sqrt[l (l + 1) - m (m + 1)] K[l, 
            m + 1];
Sm ** K[l_, m_] ^= 
        Sqrt[l (l + 1) - m (m - 1)] K[l, 
            m - 1];
Sz ** K[l_, m_] ^= m K[l, m];

List of states

Klist = {K[1, 1], K[1, 0], K[1, -1]}

Tomat[e_] := 
      Table[Coefficient[e, Klist[[j]]], {j, 1, 3}]

Matrix  < m | H | n>

Hmat = mB Tomat[Thread[ Sx ** Klist]] // MatrixForm

List of density matrix elements for solution of differential equation

rhomat = Table[rho[m, n][t], {m, 1, 3}, {n, 1, 3}]

Initial conditions

rhoinit = 
      Table[If[m == 1 && n == 1, 1, 0], {m, 1, 3}, {n, 1, 
          3}]

Formal defintion of time derivative

rhodot = Table[rho[m, n]'[t], {m, 1, 3}, {n, 1, 3}]

Liouville equaiton

LiouvEq = 
      Thread[Flatten[
            rhodot] == -I/h Flatten[
              Hmat . rhomat - rhomat . Hmat]]

Initial conditions in the form needed for solution of differential equation

InitCond = 
      Thread[Flatten[rhomat /. t -> 0] == 
          Flatten[rhoinit]]

Solution of differential equation

sol = DSolve[Union[LiouvEq, InitCond], Flatten[rhomat], t]

rhomat /. sol[[1]] // MatrixForm

Expressions for expectation values of spin components in terms of the density matrix

SzExp = Tr[Tomat[Thread[ Sz ** Klist]] . rhomat]

SxExp = Tr[Tomat[Thread[ Sx ** Klist]] . rhomat]

SyExp = Tr[Tomat[Thread[ Sy ** Klist]] . rhomat]

Spin oscillates as one would expect

Plot[SzExp /. sol[[1]] /. (mB/h) -> 1, {t, 0, 
        10}]

Plot[SyExp /. sol[[1]] /. (mB/h) -> 1, {t, 0, 
        10}]

Plot[SxExp /. sol[[1]] /. (mB/h) -> 1, {t, 0, 
        10}]

Density matrix after 1/4 of rotation

(rhomat /. sol[[1]] /. t -> Pi h/(2 mB)) // 
      MatrixForm

Spherical tensor components of density matrix element  < J, m | ρ | Jp,  
      mp >  
(see for example Appendix G  .4 of Budker et al book)

SpherCom[J_, m_, Jp_, mp_] := 
      Sum[T[L, m - mp, J, Jp] (-1)^(mp - Jp) 
          ClebschGordan[{J, m}, {Jp, -mp}, {L, m - mp}], {L, 0, 
          2 J}]

Representation of initial density matrix in terms of spherical tensor components

SpherInit = 
      Sum[If[m == 1 && n == 1, 1, 0] SpherCom[1, m, 1, 
            n], {m, -1, 1}, {n, -1, 1}]

Rotation matrix U describing coupling between m and n components of spherical tensor of rank J under rotation by an anlge al  around a rotation axis specified by angles th and ph in polar coordinates 
(these are given for example by Sec .  4.5  .2,  
        Eq .  (8) and Sec .  
            4.15  .2  Eq .  ( 
              2) in Varshilovich book on angular momentum,  
        although with slightly different normalizations)

RotU[J_, m_, n_, th_, ph_, al_] := 
        4 Pi/Sqrt[2 (2*J + 1)] 
          Sum[(-1)^(m - n) ksi[L, J, al] 
              ClebschGordan[{J,  n},  {L,  m - n},  {J, 
                  m}] SphericalHarmonicY[L, n - m, th, ph], {L, 0, 
              2  J}];
ksi[L_, J_, m_] := 
      1/Sqrt[2 Pi] 
        Sum[Exp[-I m n] 
            ClebschGordan[{J,  n},  {J,  -n},  {L, 
                0}] (-1)^(n - J), {n, -J, J}]

Transformation of the initial density matrix under rotation by Pi/2 around x axis

RotSpher = 
      Simplify[SpherInit /. (T[L_, m_, J_, Jp_] :> 
              Sum[T[L, k, J, Jp] RotU[L, m, k, Pi/2, 0, Pi/2], {k, -L, 
                  L}])]

Inverse transformation of the spherical tensors into a density matrix

ZeemCom[L_, M_, J_, Jp_] := 
      Sum[Outer[Times, Tomat[K[J, m]], 
            Tomat[K[Jp, m - M]]] (-1)^(m - M - Jp) 
          ClebschGordan[{J, m}, {Jp, M - m}, {L, M}], {m, -J, 
          J}]

Resulting rotated density matrix

Simplify[RotSpher /. T[L_, M_, J_, Jp_] :> ZeemCom[L, M, J, Jp]] // 
      MatrixForm