setup the system

This loads the package.

<<AtomicDensityMatrix`

We define an atomic system consisting of two even-parity lower states and two odd-parity upper states. We apply a light field with components at frequencies (Subscript[ω, 1]) (near resonant with the |1⟩->|3⟩ transition), (Subscript[ω, 2]) (near resonant with the |2⟩->|3⟩ transition), and  (Subscript[ω, c]) (near resonant with the |2⟩->|4⟩ transition).

Define the atomic system.

system=Sublevels[{
AtomicState[1,J->0,L->0,S->0,NaturalWidth->0,Energy->0,Polarizability[0]->0],
AtomicState[2,J->1,L->1,S->0,NaturalWidth->Γ,Energy->ω0,Polarizability[0]->-2/3,Polarizability[2]->2/3,BranchingRatio[1]->1]
}];

Define the optical field

field=((2(Sqrt[6]))/(ReducedME[1,{Dipole,1},2])){(E^(I (ϕx-t ω))) ΩRx ,(E^(I (ϕy-t ω))) ΩRy ,0}

The Hamiltonian for the system subject to the optical field.

MatrixForm[H=Expand@Hamiltonian[system,ElectricField->field,MagneticField->{Ωrf Sin[ωrf t],0,ΩL}/BohrMagneton]+Hamiltonian[system,ElectricField->{0,0,(Sqrt[Δ0])},Interaction->{Polarizability}]//FullSimplify]

The level diagram for the system.

LevelDiagram[system,H/.{ω0->1,Ωrf->.5,ΩL->.3,Δ0->.1}]

Apply the rotating-wave approximation to the Hamiltonian.

(Hrwa=RotatingWaveApproximation[system,H,ω]/.ω->ω0+Δ)//MatrixForm

IntrinsicRelaxation and TransitRelaxation supply the relaxation matrices.

MatrixForm[relax=IntrinsicRelaxation[system]+TransitRelaxation[system,γt]]

Remove explict time dependence from the density matrix.

SetOptions[DensityMatrix,TimeDependence->False,ComplexExpandVariables->False]

OpticalRepopulation and TransitRepopulation supply the repopulation matrices.

MatrixForm[repop=OpticalRepopulation[system]+TransitRepopulation[system,γt]/.BranchingRatio[a_,b_]->(Subscript[R, a,b])]

Here are the evolution equations.

TableForm[eqs=LiouvilleEquation[system,Hrwa,relax,repop]//Expand,TableHeadings->{DMVariables[system],None}]

Convert to c form.

Collect[eqs[[All,2]]/.{(E^(I ϕx)) ->Ex/ΩRx,(E^(-I ϕx)) ->Exc/ΩRx,(E^(I ϕy)) ->Ey/ΩRy,(E^(-I ϕy)) ->Eyc/ΩRy, γt->gt,Γ->g0,ΩL->deltaL, Δ0->delta0z[z],Δ->delta},DMElementPattern[],FullSimplify]/. Sin[t ωrf]->(Sqrt[2])rf/ Ωrf/.ρ->r;
MapThread[Equal,{DMVariables[system]/.ρ->dr,%}]/.Complex[0,a_]->a i/.Table[{Label[system[[i]]],M[system[[i]]]}->i,{i,Length[system]}];
DeleteCases[%,(Subscript[dr, a_,b_])==_/;b<a];
StringJoin[StringReplace[ToString@CForm[#]<>";\n",{"=="->"=","Subscript(dr,"~~a_~~","~~b_~~")":>"dr"<>a<>b<>"_dt","Subscript("~~r_~~","~~a_~~","~~b_~~")":>r<>a<>b}]&/@%]
Export[ToFileName[NotebookDirectory[],"code.txt"],%]

Find medium polarization

TensorForm[dcart=Most@ToCartesian@WignerEckart[system,{Dipole,1}]]

4π I ω0 n ((ReducedME[1,{Dipole,1},2])/(2(Sqrt[6]))) Sum[DensityMatrix[system][[i,j]]#[[j,i]],{i,Length[system]},{j,i}]&/@dcart//Simplify
ExpandDipoleRME[system,% ]/.ω0->2π/λ//Simplify
% eta/(((3 n Γ (λ^2) )/(16 π)))/.Table[{Label[system[[i]]],M[system[[i]]]}->i,{i,Length[system]}]/.ρ->r