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

During the pulse of Bx field we can ignore other terms since they are small and pulse duration is small relative to T.

Hx = Bx Tomat[Thread[ Sx ** Klist]]

List of density matrix elements for solution of differential equation

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

Initial conditions :   all atoms are in m = 
        0 state

rhoinit = 
      Table[If[m == 2 && n == 2, 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[
              Hx . rhomat - rhomat . Hx]]

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

Now we want all atoms in m = ±1 state, this happens when   2  Bx t/h=π

rhoinit1 = rhomat /. sol[[1]] /. 
        t -> (Pi h/(2  Bx))

New initial conditions

InitCond1 = 
      Thread[Flatten[rhomat /. t -> 0] == 
          Flatten[rhoinit1]]

The Hamilitonian for Bz,  dE and (AE^2) interaction

H1 = {{Bz + dE + AE, 0, 0}, {0, 0, 0}, {0, 
          0, -Bz - dE + AE}}

LiouvEq1 = 
      Thread[Flatten[
            rhodot] == -I/h Flatten[
              H1 . rhomat - rhomat . H1]]

sol1 = 
      DSolve[Union[LiouvEq1, InitCond1], Flatten[rhomat], t]

Note that (AE^2) term disappears

rhoinit2 = rhomat /. sol1[[1]] /. t -> T

Use this density matrix for the initial conditions before the second pulse

InitCond2 = 
      Thread[Flatten[rhomat /. t -> 0] == 
          Flatten[rhoinit2]]

LiouvEq2 = 
      Thread[Flatten[
            rhodot] == -I/h Flatten[
              Hx . rhomat - rhomat . Hx]]

sol2 = 
      DSolve[Union[LiouvEq2, InitCond2], Flatten[rhomat], t]

Take the length of the second pulse to be the same as the first

rhomat /. sol2[[1]] /. 
      t -> (Pi h/(2  Bx))

Pop = 
      FullSimplify[rho[2, 2][t] /. sol2[[1]] /. 
          t -> (Pi h/(2  Bx))]

Good, diagonal elements have to be real

Plot[
      Pop /. {dE -> 0, h -> 1, T -> 1}, {Bz, -5, 
        5}]

For maximum sensitivity to dE term we need  Bz T/h = 
      Pi/4

In a real experiment one has to consider  what happens if we did not ignore (AE^2) term during Bx pulse,  
      if there is some By field,  
      and if the Bx pulse is not exactly the right length to find out to what extend the sensitivity to (AE^2) can be suppressed.