General setup

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={
AtomicState[1,NaturalWidth->0,Parity->Even],
AtomicState[2,NaturalWidth->0,Parity->Even],
AtomicState[3,Energy->0,NaturalWidth->(Subscript[Γ, 3]),Parity->Odd],
AtomicState[4,NaturalWidth->(Subscript[Γ, 4]),Parity->Odd]
};

Define the optical field with three frequencies, (Subscript[ω, 1]), (Subscript[ω, 2]), and (Subscript[ω, c]).

field=OpticalField[ω1,{Ω1/ReducedME[1,{Dipole,1},3],ϕ1}]+OpticalField[ω2,{Ω2/ReducedME[2,{Dipole,1},3],ϕ2}]+OpticalField[ωc,{Ωc/ReducedME[2,{Dipole,1},4],ϕc}]

The Hamiltonian for the system subject to the optical field. Each field is assumed to interact with only one transition—the other terms are set to zero.

MatrixForm[H=Expand@Hamiltonian[system,ElectricField->field]/. Cos[_] ReducedME[_,{Dipole,1},_]->0]

The level diagram for the system.

LevelDiagram[system,H/.{Energy[1]->-1.5,Energy[2]->-1,Energy[4]->.5}]

Apply the rotating-wave approximation to the Hamiltonian.

(Hrwa=RotatingWaveApproximation[system,H/.ωc->ω2+ω4,{ω1,ω2,ω4},TransformMatrix->MatrixExp[-I t DiagonalMatrix[{-ω1,-ω2,0,ω4}]]]/.{ω1->-Energy[1]+δ1,ω2->-Energy[2]+δ2,ω4->Energy[4]+δc-δ2})//MatrixForm

IntrinsicRelaxation and TransitRelaxation supply the relaxation matrices.

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

Version using complex DM variables

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

Make plots from paper.

Plot[Re[(Subscript[ρ, 1,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 5.9,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 5.,Ω1->2π .005,Ω2->2π 5.,δ2->0,δ1->2 d},DMVariables[system]][[1]]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-5 (10^(-4)),5 (10^(-4))}}]

Plot[Re[(Subscript[ρ, 2,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 5.9,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 5.,Ω2->2π .005,Ω1->2π 5.,δ1->0,δ2->2 d},DMVariables[system]][[1]]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-3 (10^(-4)),3 (10^(-4))}}]

Plot[Im[(Subscript[ρ, 1,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 3.,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 3.,Ω1->2π .005,Ω2->2π 5.,δ2->0,δ1->2 d},DMVariables[system]][[1]]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-8 (10^(-4)),0 (10^(-4))}}]

Plot[Im[(Subscript[ρ, 1,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 3.,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 6.,Ω1->2π .005,Ω2->2π 5.,δ2->0,δ1->2 d},DMVariables[system]][[1]]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-6 (10^(-4)),0 (10^(-4))}}]

MapThread[Equal,{DMVariables[system]/.ρ->dr,Collect[eqs[[All,2]]/.Complex[0,a_]->a i/.{(E^(i ϕ1))Ω1->2 E1,(E^(-i ϕ1))Ω1->2 E1c,(E^(i ϕ2))Ω2->2 E2,(E^(-i ϕ2))Ω2->2 E2c,(E^(i ϕc))Ωc->2 Ec,(E^(-i ϕc))Ωc->2 Ecc,δ1->d1,δ2->d2,δc->d3,γt->gt,(Subscript[Γ, 3])->G3,(Subscript[Γ, 4])->G4}/.ρ->r,i,FullSimplify]}]
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"],%]

Version using real DM variables

Remove explict time dependence from the density matrix.

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

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

Make plots from paper.

Plot[(Subscript[ρ, Re,1,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 5.9,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 5.,Ω1->2π .005,Ω2->2π 5.,δ2->0,δ1->2 d},DMVariables[system]][[1]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-5 (10^(-4)),5 (10^(-4))}}]

Plot[Re[(Subscript[ρ, 2,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 5.9,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 5.,Ω2->2π .005,Ω1->2π 5.,δ1->0,δ2->2 d},DMVariables[system]][[1]]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-3 (10^(-4)),3 (10^(-4))}}]

Plot[Im[(Subscript[ρ, 1,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 3.,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 3.,Ω1->2π .005,Ω2->2π 5.,δ2->0,δ1->2 d},DMVariables[system]][[1]]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-8 (10^(-4)),0 (10^(-4))}}]

Plot[Im[(Subscript[ρ, 1,3])/.NSolve[eqs/.{ϕ1->0,ϕ2->0,ϕc->0,γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 3.,(Subscript[R, 4,1])->0,(Subscript[R, 4,2])->1,(Subscript[R, 3,_])->.5,δc->0,Ωc->2π 6.,Ω1->2π .005,Ω2->2π 5.,δ2->0,δ1->2 d},DMVariables[system]][[1]]],{d,-60,60},Frame->True,PlotRange->{{-60,60},{-6 (10^(-4)),0 (10^(-4))}}]

eqs[[All,2]]/.{Ω1->2 E1, Ω2->2E2,Ωc->2Ec,δ1->d1,δ2->d2,δc->d3,γt->gt,(Subscript[Γ, 3])->G3,(Subscript[Γ, 4])->G4,Complex[0,a_]->a i,ρ->r,Cos->cos,Sin->sin}

MapThread[Equal,{DMVariables[system]/.ρ->dr,Collect[eqs[[All,2]]/.{Ω1->2 E1, Ω2->2E2,Ωc->2Ec,δ1->d1,δ2->d2,δc->d3,γt->gt,(Subscript[Γ, 3])->G3,(Subscript[Γ, 4])->G4,Complex[0,a_]->a i,ρ->r,Cos->cos,Sin->sin},i,FullSimplify]}]
StringJoin[StringReplace[ToString@CForm[#]<>";\n",{"=="->"=","Subscript(dr,"~~f1_~~f2_~~","~~a_~~","~~b_~~")":>"d"<>f1<>f2<>"r"<>a<>b<>"_dt","Subscript(r,"~~f1_~~f2_~~","~~a_~~","~~b_~~")":>f1<>f2<>"r"<>a<>b}]&/@%]
Export[ToFileName[NotebookDirectory[],"code.txt"],%]