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),  (Subscript[ω, 3]) (near resonant with the |2⟩->|4⟩ transition), and (Subscript[ω, 4]) (near resonant with the |1⟩->|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]), (Subscript[ω, 3]), and (Subscript[ω, 4]).

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

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/.{ω3->ω2+ω43,ω4->ω1+ω41},{ω1,ω2,ω43},
TransformMatrix->MatrixExp[-I t DiagonalMatrix[{-ω1,-ω2,0,ω43}]]]/.
{ω1->-Energy[1]+δ1,ω2->-Energy[2]+δ2,ω43->Energy[4]+δ3-δ2,ω41->Energy[4]+δ4-δ1})//MatrixForm //Simplify

Assume that it is degenerate four-wave mixing  (Subscript[ω, 1])-(Subscript[ω, 2]) =  (Subscript[ω, 4])-(Subscript[ω, 3]).

Hrwa=Hrwa/. δ4->δ1-δ2+δ3

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

Atomic and thus fixed parameters

phparams={ϕ1->0,ϕ2->0,ϕ3->0,ϕ4->0};

decayparams={γt->2π .01,(Subscript[Γ, 3])->2 2π 2.7,(Subscript[Γ, 4])->2 2π 3,(Subscript[R, 4,1])->0.5,(Subscript[R, 4,2])->0.5,(Subscript[R, 3,_])->.5};

From here tunable parameters

detunparams={δ1->2 d,δ2->0,δ3->0};

omegaparams={Ω1->2π .005,Ω2->2π 15.,Ω3->2π 10.,Ω4->2π 0.001};

Make plots from paper.

params=Join[phparams, decayparams, detunparams, omegaparams];

Plot[{Re[(Subscript[ρ, 1,3])/.NSolve[eqs/.params,DMVariables[system]][[1]]],Im[(Subscript[ρ, 1,3])/.NSolve[eqs/.params,DMVariables[system]][[1]]]},{d,-60,60},Frame->True,PlotPoints-> 50,MaxRecursion->1,PlotRange->{{-60,60},Full}]

Plot[{Re[(Subscript[ρ, 1,4])/.NSolve[eqs/.params,DMVariables[system]][[1]]],Im[(Subscript[ρ, 1,4])/.NSolve[eqs/.params,DMVariables[system]][[1]]]},{d,-60,60},Frame->True,PlotPoints-> 50,MaxRecursion->1,PlotRange->{{-60,60},All}]