General setup

Mathematica 8 defines WignerD, so we need to unprotect it.

If[$VersionNumber≥8.,Unprotect[WignerD],Remove[WignerD]]

Load the package.

<<AtomicDensityMatrix`

Set options.

SetOptions[DensityMatrix,ComplexExpandVariables->Subscript,TimeDependence->True];

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

SetOptions[NDSolve,PrecisionGoal->1,AccuracyGoal->2,MaxSteps->(10^6)];

SetOptions[LevelDiagram,JLabel->True,MLabel->True,Arrowheads->{-.07,.07},MLabelPosition->{Top}];

Atomic system setup

Define the atomic system.

system=Sublevels[{
AtomicState[1,J->0,L->0,S->0,NaturalWidth->0,Energy->0,Parity->Odd],
AtomicState[2,J->1,L->1,S->0,NaturalWidth->Γ2,Energy->ω2,BranchingRatio[1]->1,BranchingRatio[3]->0,Parity->Even],
AtomicState[3,J->0,L->0,S->0,NaturalWidth->Γ3,Energy->ω3,BranchingRatio[1]->0,BranchingRatio[2]->0,Parity->Odd]
}];
TableForm[system=DeleteStates[system,Label->2,M->0]]

Define the optical field.

field=ToCartesian[Transpose@{{(((E^(I(-t ω))) ((Subscript[Ω, Re,1])+I (Subscript[Ω, Im,1])))/(ReducedME[1,{Dipole,1},2])),0,(((E^(I (-t ω))) ((Subscript[Ω, Re,-1])+I (Subscript[Ω, Im,-1])))/(ReducedME[1,{Dipole,1},2]))}}]//FullSimplify

The Hamiltonian for the system subject to the optical field.

MatrixForm[Ho=Expand@Hamiltonian[system,ElectricField->field,MagneticField->{0,0,ΩL/BohrMagneton}]/. ReducedME[2,{Dipole,1},3]->0]

The Hamiltonian for the system subject to the coupling field

MatrixForm[HE=Expand@Hamiltonian[system,ElectricField->{Ec(Sqrt[6])/ReducedME[2,{Dipole,1},3],0,0},Interaction->{E1}]/. ReducedME[1,{Dipole,1},2]->0]

The total Hamiltonian

MatrixForm[H=Ho+HE//TrigToExp]

The level diagram showing field detunings. Dashed arrows indicate dc fields. The level diagram is set to not evaluate, because it needs an updated version of the package that I haven’t released yet.

LevelDiagram[system/.(Parity->_)->Sequence[],H/.{ω2->3.,ω3->3.5,ΩL->1.,ω->3.}]

Apply the rotating-wave approximation to the Hamiltonian.

(Hrwa=
RotatingWaveApproximation[system,H,{ω},
TransformMatrix->MatrixExp[-I t DiagonalMatrix[{0,ω,ω,ω}]]]/.ω->ω2+Δ/.ω3->ω2+Δ23)//MatrixForm //Simplify

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

Density-matrix and field variables at one point i

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

Density matrix evolution equations at one point i

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]/.(Subscript[Ω, a_,b_])->(Subscript[Ω, a,b,i])[t]/.Δ23->(Subscript[Δ23, i])]

Initial conditions for density matrix:

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

Here we find the field evolution equations with finite difference approximation (second-order upwind scheme) for co-propagating beams. h is the the grid spacing in the spatial dimension.

polarizationcomponents=ExprToReIm[ExpandDipoleRME[system,4π I ω2 n0 ReducedME[1,{Dipole,1},2]Components@ ECT@Sum[DensityMatrix[system][[j,i]]WignerEckart[system[[i]],{Dipole,1},system[[j]]],{i,Length[system]},{j,i}]/.ReducedME[2,{Dipole,1},3]->0/.(Subscript[ρ, r_,a_,b_])[t]-> (Subscript[ρ, r,a,b,i])[t]]/.ω2->2π/λ/.n0 (λ^2)Γ2->4 π η/3]

fieldcomponents=ExprToReIm[Expand[ReducedME[1,{Dipole,1},2]Components@ToContravariant[field]]/.t->0]/.(Subscript[Ω, a_,b_])->(Subscript[Ω, a,b,i])[t]

UFDWeights[m_, n_, s_] := CoefficientList[Normal[Series[(x^s)((Log[x])^m),{x,1,n}]/(h^m)],x]

fieldeqs=Thread[D[fieldcomponents,t]==c polarizationcomponents-c UFDWeights[1,2,2].Table[(fieldcomponents/.i->i+j),{j,-2,0}]]

Initial conditions for fields (assume uniform in space).

initfields={(Subscript[Ω, Re,1,i])[t0]==(Subscript[Ω0, 1]),(Subscript[Ω, Re,-1,i])[t0]==(Subscript[Ω0, -1]),(Subscript[Ω, Im,1,i])[t0]==0,(Subscript[Ω, Im,-1,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,-1,0])[t]==(Subscript[Ω0, -1]),(Subscript[Ω, Im,1,0])[t]==0,(Subscript[Ω, Im,-1,0])[t]==0}

Results

setup

Choose number of spatial points.

n=40;

All system variables for all spatial points.

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

Equations for all points for the co-propagating case.

alleqs=(#[[1]]-#[[2]]==0&)/@(Flatten@Table[Join[eqs,fieldeqs],{i,0,n}]/.{(Subscript[Ω, r_,s_,i_])[t_]/;i<0->(Subscript[Ω, r,s,0])[t],(Subscript[Ω, r_,s_,i_])[t_]/;i>n->(Subscript[Ω, r,s,n])[t]}/.(((D[#[[1]],t]==_)->#&)/@(boundaryconds/.n0->n)));

allinits=Flatten@Table[Join[inits,initfields],{i,0,n}];

Relabel variables to enforce efficient ordering in NDSolve.

allvars1=Table[(Subscript[y, i])[t],{i,Length[allvars]}];

varreps=MapThread[Rule,{allvars,allvars1}];

alleqs1=alleqs/.Dispatch[varreps]/.Dispatch[D[varreps,t]];

allinits1=allinits/.(varreps/.t->t0);

plots

Here we choose parameters and integrate the equations. Blue is (SuperMinus[σ]) field intensity, magenta is (SuperPlus[σ])field intensity.

Parameters:
(Subscript[Ω0, -1]): input Rabi frequency for (SuperMinus[σ]) (rad/s) 
(Subscript[Ω0, 1]):  input Rabi frequency for (SuperPlus[σ])(rad/s) 
Γ2: natural width of J=1 level (rad/s) 
Γ3: natural width of upper J=0 level (rad/s)
γt: transit relaxation rate (rad/s)
Ec: static electric field coupling strength  (rad/s)
Δ: Light detuning from J=1 level (rad/s)
ΩL: Zeeman shift (rad/s)
(Subscript[Δ23, i_]): frequency splitting between upper J=1 and J=0 level as a function of the position i, as i varies from zero to n (rad/s)
h: distance between grid points (m)
c: light speed (m/s)
η: light coupling constant (rad/s/m)

params={(Subscript[Ω0, -1])->0.,(Subscript[Ω0, 1])->1. 2π (10^6),Γ2->1. 2π (10^6),Γ3->0,γt->.001 (10^6),Ec->6. 2π (10^6),Δ->0. 2π (10^6),ΩL->50. 2π (10^6),c->3. (10^8),η->3. (((794.7 (10^(-9))))^2)(.5 (10^15))(2 π (10^6))/(4.π),h->5./n,(Subscript[Δ23, i_])->(-2.+4i/n) 2π  (10^6),t0->0};
AbsoluteTiming[sol=NDSolve[Evaluate[N[{alleqs1,allinits1}/.params]],allvars1,{t,0,8. (10^(-6))}];]
ListLinePlot[Transpose@Table[Evaluate[{{i h,((Abs[((Subscript[Ω, Re,-1,i])[t]+I (Subscript[Ω, Im,-1,i])[t])])^2)},{i h,((Abs[((Subscript[Ω, Re,1,i])[t]+I (Subscript[Ω, Im,1,i])[t])])^2)}}/.params],{i,0,n}]/.varreps/.sol[[1]]/.t->8. (10^(-6)),FrameLabel->{"Position (m)","Intensity (Subscript[Ω, R]^2) (rad/s)^2"}]

Same plot with flipped magnetic field

params={(Subscript[Ω0, -1])->0.,(Subscript[Ω0, 1])->1. 2π (10^6),Γ2->1. 2π (10^6),Γ3->0,γt->.001 (10^6),Ec->6. 2π (10^6),Δ->0. 2π (10^6),ΩL->-50. 2π (10^6),c->3. (10^8),η->3. (((794.7 (10^(-9))))^2)(.5 (10^15))(2 π (10^6))/(4.π),h->5./n,(Subscript[Δ23, i_])->(-2.+4i/n) 2π  (10^6),t0->0};
AbsoluteTiming[sol=NDSolve[Evaluate[N[{alleqs1,allinits1}/.params]],allvars1,{t,0,8. (10^(-6))}];]
ListLinePlot[Transpose@Table[Evaluate[{{i h,((Abs[((Subscript[Ω, Re,-1,i])[t]+I (Subscript[Ω, Im,-1,i])[t])])^2)},{i h,((Abs[((Subscript[Ω, Re,1,i])[t]+I (Subscript[Ω, Im,1,i])[t])])^2)}}/.params],{i,0,n}]/.varreps/.sol[[1]]/.t->8. (10^(-6)),FrameLabel->{"Position (m)","Intensity (Subscript[Ω, R]^2) (rad/s)^2"}]