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)

Work with real and imaginary parts of the density matrix variables.

SetOptions[DensityMatrix,ComplexExpandVariables->Subscript];

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]).

SetOptions[OpticalField,CartesianCoordinates->{x,y,z}]

field=OpticalField[{(Subscript[ω, 1]),(Subscript[k, 1])},{(Subscript[Ω, 1])/ReducedME[1,{Dipole,1},3],(Subscript[ϕ, 1])}]+OpticalField[{(Subscript[ω, 2]),(Subscript[k, 2])},{(Subscript[Ω, 2])/ReducedME[2,{Dipole,1},3],(Subscript[ϕ, 2])}]+OpticalField[{(Subscript[ω, 3]),(Subscript[k, 3])},{(Subscript[Ω, 3])/ReducedME[2,{Dipole,1},4],(Subscript[ϕ, 3])}]+
OpticalField[{(Subscript[ω, 4]),(Subscript[k, 4])},{(Subscript[Ω, 4])/ReducedME[1,{Dipole,1},4],(Subscript[ϕ, 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/.{(Subscript[ω, 3])->(Subscript[ω, 2])+ω43,(Subscript[ω, 4])->(Subscript[ω, 1])+ω41}/.{(Subscript[k, 3])->(Subscript[k, 2])+k43,(Subscript[k, 4])->(Subscript[k, 1])+k41},{(Subscript[ω, 1]),(Subscript[ω, 2]),ω43},
TransformMatrix->MatrixExp[-I t DiagonalMatrix[{-(Subscript[ω, 1]),-(Subscript[ω, 2]),0,ω43}]+I z DiagonalMatrix[{-(Subscript[k, 1]),-(Subscript[k, 2]),0,k43}]]]/.
{(Subscript[ω, 1])->-Energy[1]+δ1,(Subscript[ω, 2])->-Energy[2]+δ2,ω43->Energy[4]+δ3-δ2,ω41->Energy[4]+δ4-δ1}/.{k43 ->(Subscript[k, 3])-(Subscript[k, 2]),k41 ->(Subscript[k, 4])-(Subscript[k, 1])})//MatrixForm //Simplify

Set (Subscript[ω, 1])-(Subscript[ω, 2]) =  (Subscript[ω, 4])-(Subscript[ω, 3]), and phase-matching: -(Subscript[k, 1])+(Subscript[k, 2])-(Subscript[k, 3])+(Subscript[k, 4])=0

(Hrwa=Hrwa/. δ4->δ1-δ2+δ3/. (Subscript[k, 4])->(Subscript[k, 1])-(Subscript[k, 2])+(Subscript[k, 3])/.(E^(I (Subscript[ϕ, j_]))) (Subscript[Ω, j_])->(Subscript[Ω, Re,j,i])[t]+I (Subscript[Ω, Im,j,i])[t]/.(E^(-I (Subscript[ϕ, j_]))) (Subscript[Ω, j_])->(Subscript[Ω, Re,j,i])[t]-I (Subscript[Ω, Im,j,i])[t])//MatrixForm

IntrinsicRelaxation and TransitRelaxation supply the relaxation matrices.

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

OpticalRepopulation and TransitRepopulation supply the repopulation matrices.

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

Density-matrix and field variables for one point

vars=Join[DMVariables[system]/.(Subscript[ρ, r_,a_,b_])[t]-> (Subscript[ρ, r,a,b,i])[t],Table[{(Subscript[Ω, Re,j,i])[t],(Subscript[Ω, Im,j,i])[t]},{j,4}]]//Flatten

Density matrix evolution equations:

TableForm[eqs=Expand@LiouvilleEquation[system,Hrwa,relax,repop]/.{(((Subscript[ρ, r_,a_,b_])))'[t]->(((Subscript[ρ, r,a,b,i])))'[t]}/.(Subscript[ρ, r_,a_,b_])[t]-> (Subscript[ρ, r,a,b,i])[t]]

Initial conditions for density matrix:

inits=InitialConditions[system,TransitRepopulation[system,1],t0]/. (Subscript[ρ, r_,a_,b_])[t_]-> (Subscript[ρ, r,a,b,i])[t]

Field evolution equations with finite difference approximation (first-order upwind scheme) for co-propagating beams. h is the the grid spacing in the spatial dimension.

fieldeqs={
(((Subscript[Ω, Re,1,i])))'[t]==η c ((Subscript[ρ, Im,1,3,i])[t])-c((Subscript[Ω, Re,1,i])[t]-(Subscript[Ω, Re,1,i-1])[t])/h,(((Subscript[Ω, Im,1,i])))'[t]==η c  (Subscript[ρ, Re,1,3,i])[t]-c((Subscript[Ω, Im,1,i])[t]-(Subscript[Ω, Im,1,i-1])[t])/h,
(((Subscript[Ω, Re,2,i])))'[t]==η c ((Subscript[ρ, Im,2,3,i])[t])-c((Subscript[Ω, Re,2,i])[t]-(Subscript[Ω, Re,2,i-1])[t])/h,(((Subscript[Ω, Im,2,i])))'[t]==η c  (Subscript[ρ, Re,2,3,i])[t]-c((Subscript[Ω, Im,2,i])[t]-(Subscript[Ω, Im,2,i-1])[t])/h,
(((Subscript[Ω, Re,3,i])))'[t]==η c ((Subscript[ρ, Im,2,4,i])[t])-c((Subscript[Ω, Re,3,i])[t]-(Subscript[Ω, Re,3,i-1])[t])/h,(((Subscript[Ω, Im,3,i])))'[t]==η c  (Subscript[ρ, Re,2,4,i])[t]-c((Subscript[Ω, Im,3,i])[t]-(Subscript[Ω, Im,3,i-1])[t])/h,
(((Subscript[Ω, Re,4,i])))'[t]==η c ((Subscript[ρ, Im,1,4,i])[t])-c((Subscript[Ω, Re,4,i])[t]-(Subscript[Ω, Re,4,i-1])[t])/h,(((Subscript[Ω, Im,4,i])))'[t]==η c  (Subscript[ρ, Re,1,4,i])[t]-c((Subscript[Ω, Im,4,i])[t]-(Subscript[Ω, Im,4,i-1])[t])/h}

Field evolution equations with finite difference approximation (first-order upwind scheme) for fields 1 and 2 forward propagating, fields 3 and 4 backward propagating.

counterfieldeqs={
(((Subscript[Ω, Re,1,i])))'[t]==η c ((Subscript[ρ, Im,1,3,i])[t])-c((Subscript[Ω, Re,1,i])[t]-(Subscript[Ω, Re,1,i-1])[t])/h,(((Subscript[Ω, Im,1,i])))'[t]==η c  (Subscript[ρ, Re,1,3,i])[t]-c((Subscript[Ω, Im,1,i])[t]-(Subscript[Ω, Im,1,i-1])[t])/h,
(((Subscript[Ω, Re,2,i])))'[t]==η c ((Subscript[ρ, Im,2,3,i])[t])-c((Subscript[Ω, Re,2,i])[t]-(Subscript[Ω, Re,2,i-1])[t])/h,(((Subscript[Ω, Im,2,i])))'[t]==η c  (Subscript[ρ, Re,2,3,i])[t]-c((Subscript[Ω, Im,2,i])[t]-(Subscript[Ω, Im,2,i-1])[t])/h,
(((Subscript[Ω, Re,3,i])))'[t]==η c ((Subscript[ρ, Im,2,4,i])[t])+c((Subscript[Ω, Re,3,i+1])[t]-(Subscript[Ω, Re,3,i])[t])/h,(((Subscript[Ω, Im,3,i])))'[t]==η c  (Subscript[ρ, Re,2,4,i])[t]+c((Subscript[Ω, Im,3,i+1])[t]-(Subscript[Ω, Im,3,i])[t])/h,
(((Subscript[Ω, Re,4,i])))'[t]==η c ((Subscript[ρ, Im,1,4,i])[t])+c((Subscript[Ω, Re,4,i+1])[t]-(Subscript[Ω, Re,4,i])[t])/h,(((Subscript[Ω, Im,4,i])))'[t]==η c  (Subscript[ρ, Re,1,4,i])[t]+c((Subscript[Ω, Im,4,i+1])[t]-(Subscript[Ω, Im,4,i])[t])/h}

Initial conditions for fields (assume uniform in space).

initfields={(Subscript[Ω, Re,1,i])[t0]==(Subscript[Ω0, 1]),(Subscript[Ω, Re,2,i])[t0]==(Subscript[Ω0, 2]) Exp[-(((t0/tp))^2)],(Subscript[Ω, Re,3,i])[t0]==(Subscript[Ω0, 3]),(Subscript[Ω, Re,4,i])[t0]==(Subscript[Ω0, 4]) Exp[-(((t0/tp))^2)],(Subscript[Ω, Im,1,i])[t0]==0,(Subscript[Ω, Im,2,i])[t0]==0,(Subscript[Ω, Im,3,i])[t0]==0,(Subscript[Ω, Im,4,i])[t0]==0}

Boundary conditions for co-propagating fields. Fields 1 and 3 are constant, fields 2 and 4 are pulsed. The first point is 0, and the last point is n0.

boundaryconds={(Subscript[Ω, Re,1,0])[t]==(Subscript[Ω0, 1]),(Subscript[Ω, Re,2,0])[t]==(Subscript[Ω0, 2]) Exp[-((((t)/tp))^2)],(Subscript[Ω, Re,3,0])[t]==(Subscript[Ω0, 3]),(Subscript[Ω, Re,4,0])[t]==(Subscript[Ω0, 4])Exp[-((((t)/tp))^2)],(Subscript[Ω, Im,1,0])[t]==0,(Subscript[Ω, Im,2,0])[t]==0,(Subscript[Ω, Im,3,0])[t]==0,(Subscript[Ω, Im,4,0])[t]==0}

Boundary conditions for counter-propagating fields. Fields 1 and 3 are constant, fields 2 and 4 are pulsed.

counterboundaryconds={(Subscript[Ω, Re,1,0])[t]==(Subscript[Ω0, 1]),(Subscript[Ω, Re,2,0])[t]==(Subscript[Ω0, 2]) Exp[-((((t)/tp))^2)],(Subscript[Ω, Re,3,n0])[t]==(Subscript[Ω0, 3]),(Subscript[Ω, Re,4,n0])[t]==(Subscript[Ω0, 4])Exp[-((((t)/tp))^2)],(Subscript[Ω, Im,1,0])[t]==0,(Subscript[Ω, Im,2,0])[t]==0,(Subscript[Ω, Im,3,n0])[t]==0,(Subscript[Ω, Im,4,n0])[t]==0}

Results

More setup

Choose number of spatial points.

n=100;

All system variables for all spatial points.

allvars=Flatten@Table[vars,{i,0,n}];

Equations for all points for the co-propagating case.

TableForm[alleqs=Join[DeleteCases[Flatten@Table[Join[eqs,inits,fieldeqs,initfields],{i,0,n}],Alternatives@@(D[#[[1]],t]==_&/@(boundaryconds/.n0->n))],(boundaryconds/.n0->n)]];

Equations for all points for the counter-propagating case.

TableForm[allcountereqs=Join[DeleteCases[Flatten@Table[Join[eqs,inits,counterfieldeqs,initfields],{i,0,n}],Alternatives@@(D[#[[1]],t]==_&/@(counterboundaryconds/.n0->n))],(counterboundaryconds/.n0->n)]];

Choose tolerance for NDSolve.

SetOptions[NDSolve,PrecisionGoal->4,AccuracyGoal->4];

SetOptions[Plot,PlotRange->All,Frame->True];
SetOptions[ListDensityPlot,PlotRange->All];

Here we choose parameters and integrate the equations for the co- and counter-propagating cases.

case I (f=0)

params={γt->0 2π 0.01 (10^6),(Subscript[Γ, 3])->2π 3.0 (10^6),(Subscript[Γ, 4])->2π 3.0 (10^6),(Subscript[R, 3,1])->.5,(Subscript[R, 4,1]) ->.5,(Subscript[R, 3,2])->.5,(Subscript[R, 4,2]) ->.5,δ1->0,δ2->0,δ3->0,c->3. (10^8),η->2 3. (((794.7 (10^(-9))))^2)(10^15) 6 2 π (10^6)/(8.π),t0->-15. (10^(-6)),tp->4. (10^(-6)),(Subscript[Ω0, 1])->3. 2π  (10^6),(Subscript[Ω0, 2])->1. 2π (10^3),(Subscript[Ω0, 3])->6. 2π (10^6),(Subscript[Ω0, 4])->0 1. 2π  (10^3),h->1.5 (10^(-2))/n};

alleqs1=Expand[alleqs/.params];
sol=NDSolve[alleqs1,allvars,{t,-15. (10^(-6)),15. (10^(-6))}];

TableForm[allcountereqs1=Expand[allcountereqs/.params]];
countersol=NDSolve[allcountereqs1,allvars,{t,-15. (10^(-6)),15. (10^(-6))}];

Time dependence before and after the cell.  Black is (Subscript[Ω, 2]) pulse before the medium, blue is  (Subscript[Ω, 2]) pulse after the medium, red is (Subscript[Ω, 4]) pulse after the medium. Solid lines are co-propagating, dashed lines are fields 1 and 2 forward propagating, fields 3 and 4 backward propagating.

Plot[Evaluate[Join[{((((Subscript[Ω, Re,2,0])[t]+I (Subscript[Ω, Im,2,0])[t]))^2),((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,n])[t]+I (Subscript[Ω, Im,4,n])[t]))^2)}/.sol[[1]],{((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,0])[t]+I (Subscript[Ω, Im,4,0])[t]))^2)}/.countersol[[1]]]],{t,-5. (10^(-6)),5. (10^(-6))},PlotStyle->{Black,Blue,Red,Directive[Blue,Dashed],Directive[Red,Dashed]},FrameLabel->{"Time (μs)","Intensity ((Subscript[Ω, R]^2))"}]

Intensity of fields 2 and 4 for co-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.sol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.sol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]

Intensity of fields 2 and 4 for counter-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.countersol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.countersol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]

case II (f=0.5)

params={γt->0 2π 0.01 (10^6),(Subscript[Γ, 3])->2π 3.0 (10^6),(Subscript[Γ, 4])->2π 3.0 (10^6),(Subscript[R, 3,1])->.5,(Subscript[R, 4,1]) ->.5,(Subscript[R, 3,2])->.5,(Subscript[R, 4,2]) ->.5,δ1->0,δ2->0,δ3->0,c->3. (10^8),η->2 3. (((794.7 (10^(-9))))^2)(10^15) 6 2 π (10^6)/(8.π),t0->-15. (10^(-6)),tp->4. (10^(-6)),(Subscript[Ω0, 1])->3. 2π  (10^6),(Subscript[Ω0, 2])->1. 2π (10^3),(Subscript[Ω0, 3])->6. 2π (10^6),(Subscript[Ω0, 4])->.5 1. 2π  (10^3),h->1.5 (10^(-2))/n};

alleqs1=Expand[alleqs/.params];
sol=NDSolve[alleqs1,allvars,{t,-15. (10^(-6)),15. (10^(-6))}];

TableForm[allcountereqs1=Expand[allcountereqs/.params]];
countersol=NDSolve[allcountereqs1,allvars,{t,-15. (10^(-6)),15. (10^(-6))}];

Time dependence before and after the cell.  Black is (Subscript[Ω, 2]) pulse before the medium, blue is  (Subscript[Ω, 2]) pulse after the medium, red is (Subscript[Ω, 4]) pulse after the medium. Solid lines are co-propagating, dashed lines are fields 1 and 2 forward propagating, fields 3 and 4 backward propagating.

Plot[Evaluate[Join[{((((Subscript[Ω, Re,2,0])[t]+I (Subscript[Ω, Im,2,0])[t]))^2),((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,n])[t]+I (Subscript[Ω, Im,4,n])[t]))^2)}/.sol[[1]],{((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,0])[t]+I (Subscript[Ω, Im,4,0])[t]))^2)}/.countersol[[1]]]],{t,-5. (10^(-6)),5. (10^(-6))},PlotStyle->{Black,Blue,Red,Directive[Blue,Dashed],Directive[Red,Dashed]},FrameLabel->{"Time (μs)","Intensity ((Subscript[Ω, R]^2))"}]

Intensity of fields 2 and 4 for co-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.sol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.sol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]

Intensity of fields 2 and 4 for counter-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.countersol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.countersol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]

case III (f=1)

params={γt->0 2π 0.01 (10^6),(Subscript[Γ, 3])->2π 3.0 (10^6),(Subscript[Γ, 4])->2π 3.0 (10^6),(Subscript[R, 3,1])->.5,(Subscript[R, 4,1]) ->.5,(Subscript[R, 3,2])->.5,(Subscript[R, 4,2]) ->.5,δ1->0,δ2->0,δ3->0,c->3. (10^8),η->2 3. (((794.7 (10^(-9))))^2)(10^15) 6 2 π (10^6)/(8.π),t0->-15. (10^(-6)),tp->4. (10^(-6)),(Subscript[Ω0, 1])->3. 2π  (10^6),(Subscript[Ω0, 2])->1. 2π (10^3),(Subscript[Ω0, 3])->6. 2π (10^6),(Subscript[Ω0, 4])-> 1. 2π  (10^3),h->1.5 (10^(-2))/n};

alleqs1=Expand[alleqs/.params];
sol=NDSolve[alleqs1,allvars,{t,-15. (10^(-6)),15. (10^(-6))}];

TableForm[allcountereqs1=Expand[allcountereqs/.params]];
countersol=NDSolve[allcountereqs1,allvars,{t,-15. (10^(-6)),15. (10^(-6))}];

Time dependence before and after the cell.  Black is (Subscript[Ω, 2]) pulse before the medium, blue is  (Subscript[Ω, 2]) pulse after the medium, red is (Subscript[Ω, 4]) pulse after the medium. Solid lines are co-propagating, dashed lines are fields 1 and 2 forward propagating, fields 3 and 4 backward propagating.

Plot[Evaluate[Join[{((((Subscript[Ω, Re,2,0])[t]+I (Subscript[Ω, Im,2,0])[t]))^2),((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,n])[t]+I (Subscript[Ω, Im,4,n])[t]))^2)}/.sol[[1]],{((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,0])[t]+I (Subscript[Ω, Im,4,0])[t]))^2)}/.countersol[[1]]]],{t,-5. (10^(-6)),5. (10^(-6))},PlotStyle->{Black,Blue,Red,Directive[Blue,Dashed],Directive[Red,Dashed]},FrameLabel->{"Time (μs)","Intensity ((Subscript[Ω, R]^2))"}]

Intensity of fields 2 and 4 for co-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.sol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.sol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]

Intensity of fields 2 and 4 for counter-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.countersol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.countersol[[1]]],{t,-5. (10^(-6)),5. (10^(-6)),.2 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]

Slow light

params={γt->2π 0.01 (10^6),(Subscript[Γ, 3])->2π 3.0 (10^6),(Subscript[Γ, 4])->2π 3.0 (10^6),(Subscript[R, 3,1])->.5,(Subscript[R, 4,1]) ->.5,(Subscript[R, 3,2])->.5,(Subscript[R, 4,2]) ->.5,δ1->0,δ2->0,δ3->0,c->3. (10^8),η->2 3. (((794.7 (10^(-9))))^2)(10^15) 6 2 π (10^6)/(8.π),t0->-1.4 (10^(-6)),tp->0.1 (10^(-6)),(Subscript[Ω0, 1])->2 1.9 (10^7),(Subscript[Ω0, 2])->2 3.1 (10^5),(Subscript[Ω0, 3])->1.9 (10^7),(Subscript[Ω0, 4])->0,h->1.5 (10^(-2))/n};

alleqs1=Expand[alleqs/.params];
sol=NDSolve[alleqs1,allvars,{t,-1.4 (10^(-6)),1.4 (10^(-6))}];

TableForm[allcountereqs1=Expand[allcountereqs/.params]];
countersol=NDSolve[allcountereqs1,allvars,{t,-1.4 (10^(-6)),1.4 (10^(-6))}];

Time dependence before and after the cell.  Black is (Subscript[Ω, 2]) pulse before the medium, blue is  (Subscript[Ω, 2]) pulse after the medium, red is (Subscript[Ω, 4]) pulse after the medium. Solid lines are co-propagating, dashed lines are fields 1 and 2 forward propagating, fields 3 and 4 backward propagating.

Plot[Evaluate[Join[{((((Subscript[Ω, Re,2,0])[t]+I (Subscript[Ω, Im,2,0])[t]))^2),((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,n])[t]+I (Subscript[Ω, Im,4,n])[t]))^2)}/.sol[[1]],{((((Subscript[Ω, Re,2,n])[t]+I (Subscript[Ω, Im,2,n])[t]))^2),((((Subscript[Ω, Re,4,0])[t]+I (Subscript[Ω, Im,4,0])[t]))^2)}/.countersol[[1]]]],{t,-.5 (10^(-6)),.5 (10^(-6))},PlotStyle->{Black,Blue,Red,Directive[Blue,Dashed],Directive[Red,Dashed]},FrameLabel->{"Time (μs)","Intensity ((Subscript[Ω, R]^2))"}]

Intensity of fields 2 and 4 for co-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.sol[[1]]],{t,-.5 (10^(-6)),.5 (10^(-6)),.02 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.sol[[1]]],{t,-.5 (10^(-6)),.5 (10^(-6)),.02 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]

Intensity of fields 2 and 4 for counter-propagating case.

GraphicsRow[{ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,2,i])[t]+I (Subscript[Ω, Im,2,i])[t]))^2)}/.params/.countersol[[1]]],{t,-.5 (10^(-6)),.5 (10^(-6)),.02 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 2])((|)^2))"],ListDensityPlot[Chop@Flatten[Table[Table[Evaluate[{t (10^6),i h (10^2),((((Subscript[Ω, Re,4,i])[t]+I (Subscript[Ω, Im,4,i])[t]))^2)}/.params/.countersol[[1]]],{t,-.5 (10^(-6)),.5 (10^(-6)),.02 (10^(-6))}],{i,0,100,5}],1],FrameLabel->{"t (μs)","z (cm)"},PlotLabel->"Intensity (|(Subscript[Ω, 4])((|)^2))"]},ImageSize->8 72]