Find evolution equations

Load the package.

<<AtomicDensityMatrix`

Use density matrix variables with explict time dependence.

SetOptions[DensityMatrix,TimeDependence->True]

Pull quantum numbers and other basic info about the states from a database.

s12data=AtomicData[{"Rb",87,{"Kr",{5,"s"}},{2,"S",(1/2)}},{Energy,J,L,S,NuclearSpin,NaturalWidth,Parity}]

p12data=Append[AtomicData[{"Rb",87,{"Kr",{5,"p"}},{2,"P",(1/2)}},{J,L,S,NuclearSpin,Parity}],BranchingRatio[0]->1]

p32data=Append[AtomicData[{"Rb",87,{"Kr",{5,"p"}},{2,"P",(3/2)}},{J,L,S,NuclearSpin,Parity}],BranchingRatio[0]->1]

Using the quantum numbers, create a list of all hyperfine and Zeeman sublevels of the J states, which are labeled 0, 1, and 2 for reference. This list will be passed to the functions from the ADM package that create the DM evolution equations.

system0=Sublevels[{AtomicState[0,s12data],AtomicState[1,p12data],AtomicState[2,p32data]}];

For simplicity, we delete some excited states from the system, leaving only (Subscript[P, 1/2]) F=1 and (Subscript[P, 3/2]) F=2.

system=DeleteStates[system0,(Label==1&&F≠1)||(Label==2&&F≠2)]

Define the optical field with four frequencies, (Subscript[ω, 1]), (Subscript[ω, 2]), (Subscript[ω, 3]), and (Subscript[ω, 4]). Here (Subscript[Ω, 1]) and (Subscript[Ω, 1]) are left-circularly polarized, and (Subscript[Ω, 3]) and (Subscript[Ω, 4]) are right-circularly polarized.

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

Parameters for the OpticalField function are {frequency, wavenumber}, {electric amplitude, phase}, {rotation angle (relative to PolarizationVector), ellipticity}. Here (Subscript[Ω, i]) are the Rabi frequencies defined in terms of the dipole reduced matrix elements.

field=OpticalField[{(Subscript[ω, 1]),(Subscript[k, 1])},{(Subscript[Ω, 1])/ReducedME[0,{Dipole,1},1],(Subscript[ϕ, 1])},{0,π/4}]+
OpticalField[{(Subscript[ω, 2]),(Subscript[k, 2])},{(Subscript[Ω, 2])/ReducedME[0,{Dipole,1},1],(Subscript[ϕ, 2])},{0,π/4}]+
OpticalField[{(Subscript[ω, 3]),(Subscript[k, 3])},{(Subscript[Ω, 3])/ReducedME[0,{Dipole,1},2],(Subscript[ϕ, 3])},{0,-π/4}]+
OpticalField[{(Subscript[ω, 4]),(Subscript[k, 4])},{(Subscript[Ω, 4])/ReducedME[0,{Dipole,1},2],(Subscript[ϕ, 4])},{0,-π/4}]

The Hamiltonian for the system subject to the optical field. Each field is assumed to interact with only one transition—the replacement rule (Cos[_]|Sin[_]) ReducedME[_,{Dipole,1},_]->0 causes other terms to be set to zero.

MatrixForm[H=Expand@Hamiltonian[system,ElectricField->field,MagneticField->{0,0,ΩL/BohrMagneton}]/. (Cos[_]|Sin[_]) ReducedME[_,{Dipole,1},_]->0]

The level diagram for the system, showing optical couplings. Note that both resonant and off-resonant (counter-rotating) couplings are shown, because we have not yet performed the rotating-wave approximation.

LevelDiagram[system,H/.{Energy[0]->0,Energy[1]->2,Energy[2]->4,HyperfineA[_]->.2,HyperfineB[_]->.1,ΩL->.1},ParityOffset->False]

Here we apply the rotating-wave approximation to the Hamiltonian. Here we construct a list of frequency shifts for the list of Zeeman sublevels that will have the effect of removing the optical frequencies from the Hamiltonian. I.e., we hold (Subscript[S, 1/2])F=2 fixed, shift (Subscript[P, 1/2]) down by (Subscript[ω, 1]), shift (Subscript[S, 1/2]) F=1 down by (Subscript[ω, 1])-(Subscript[ω, 2]), and shift (Subscript[P, 3/2]) down by (Subscript[ω, 1])-(Subscript[ω, 2])+(Subscript[ω, 3]).

{Label[#],F[#]}&/@system
shifts=%/.{{0,2}->0,{1,_}->(Subscript[ω, 1]),{0,1}->(Subscript[ω, 1])-(Subscript[ω, 2]),{2,_}->(Subscript[ω, 1])-(Subscript[ω, 2])+(Subscript[ω, 3])}

Make a unitary transformation matrix to implement the chosen frequency shifts.

MatrixForm[transmat=MatrixExp[-I t DiagonalMatrix[shifts]+I z DiagonalMatrix[shifts/.ω->k]]]

Write (Subscript[ω, 4]) and (Subscript[k, 4]) in terms of the non-degenerate detuning and phase-mismatch parameters δω and δk, which we set here to zero for simplicity. We then apply the transform matrix to the Hamiltonian and set off resonant terms oscillating at harmonics of (Subscript[ω, 1]), (Subscript[ω, 2]), and (Subscript[ω, 3]) to zero.

MatrixForm[H1=RotatingWaveApproximation[system,H/.{(Subscript[ω, 4])->(Subscript[ω, 1])-(Subscript[ω, 2])+(Subscript[ω, 3])-δω,(Subscript[k, 4])->(Subscript[k, 1])-(Subscript[k, 2])+(Subscript[k, 3])-δk}/.{δω->0,δk->0},{(Subscript[ω, 1]),(Subscript[ω, 2]),(Subscript[ω, 3])},TransformMatrix->transmat]]

There are remaining fast-oscillating terms in the Hamiltonian at the difference frequency between (Subscript[ω, 1]) and (Subscript[ω, 2]), due to the fact that each optical field can interact with both ground-state hyperfine levels. We set these terms to zero, which is equivalent to assuming that each field interacts with only the F=1 or the F=2 ground-state hyperfine sublevel.

MatrixForm[H2=H1/.{(E^(ϕ_-I t (Subscript[ω, 1])+I t (Subscript[ω, 2])))->0,(E^(ϕ_+I t (Subscript[ω, 1])-I t (Subscript[ω, 2])))->0}]

Find the frequencies at which each field is assumed to be resonant with its F->F’ transition.

ω1res=Hamiltonian[{SelectState[system, Label==1&&F==1]}][[1,1]]-Hamiltonian[{SelectState[system, Label==0&&F==2]}][[1,1]]

ω2res=Hamiltonian[{SelectState[system, Label==1&&F==1]}][[1,1]]-Hamiltonian[{SelectState[system, Label==0&&F==1]}][[1,1]]

ω3res=Hamiltonian[{SelectState[system, Label==2&&F==2]}][[1,1]]-Hamiltonian[{SelectState[system, Label==0&&F==1]}][[1,1]]

Rewrite the frequencies  (Subscript[ω, 1]), (Subscript[ω, 2]), and (Subscript[ω, 3]) in terms of a detuning  (Subscript[δ, 1]), (Subscript[δ, 2]), or (Subscript[δ, 3]) from resonance with the appropriate transition between ground-state and excited state hyperfine levels, and subtract a constant term off of the diagonal to simplify the appearance.

MatrixForm[Hrwa=(H2/.{(Subscript[ω, 1])->ω1res+(Subscript[δ, 1]),(Subscript[ω, 2])->ω2res+(Subscript[δ, 2]),(Subscript[ω, 3])->ω3res+(Subscript[δ, 3])})-((3 HyperfineA[0])/4)IdentityMatrix[Length[system]]]

The level diagram showing resonant (co-rotating) optical couplings.

diagram=LevelDiagram[system,Hrwa/.{(Subscript[δ, 1])->-ω1res,(Subscript[δ, 2])->-ω2res,(Subscript[δ, 3])->-ω3res}/.{Energy[0]->0,Energy[1]->2,Energy[2]->3,HyperfineA[_]->.2,HyperfineB[_]->.1,ΩL->.1},ParityOffset->False,Epilog->{Text[Style["( ^2)(Subscript[S, 1/2])",SingleLetterItalics->False],{-2.9,0}],Text[Style["( ^2)(Subscript[P, 1/2])",SingleLetterItalics->False],{-2.9,1.6}],Text[Style["( ^2)(Subscript[P, 3/2])",SingleLetterItalics->False],{-2.9,2.7}]},ImageSize->3.5 72,ImagePadding->{{35,0},{10,10}}]

IntrinsicRelaxation and TransitRelaxation supply the matrices describing relaxation due to spontaneous decay and atomic transit, respectively. (Subscript[γ, t]) is the transit rate.

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

OpticalRepopulation and TransitRepopulation supply the matrices describing repopulation of the ground state due to spontaneous decay and atomic transit.

MatrixForm[repop=OpticalRepopulation[system]+TransitRepopulation[system,γt]]

Here are the evolution equations.

TableForm[eqs0=LiouvilleEquation[system,Hrwa,relax,repop]]

Here we pull useful numerical atomic data for Rb out of the database. These numbers are in omega units, so that the unit “Hertz” actually corresponds to rad/s.

atomicdata=Join[AtomicData[{"Rb",87,{"Kr",{5,"s"}},{2,"S",(1/2)}},{HyperfineA},Label->0],AtomicData[{"Rb",87,{"Kr",{5,"p"}},{2,"P",(1/2)}},{HyperfineA,NaturalWidth},Label->1],AtomicData[{"Rb",87,{"Kr",{5,"p"}},{2,"P",(3/2)}},{HyperfineA,HyperfineB,NaturalWidth},Label->2],{Wavelength[1]->1/Wavenumber[{"Rb",87,{"Kr",{5,"p"}},{2,"P",(1/2)}}],Wavelength[2]->1/Wavenumber[{"Rb",87,{"Kr",{5,"p"}},{2,"P",(3/2)}}]}]

Here we find DM elements that are always identically zero, so we can remove them from the evolution equations. We put in sample values for all of the parameters, use the values for the atomic data from above, and set time derivatives to zero with (Subscript[ρ, s1_,s2_])'[t]->0 to find the steady state.

sol=NSolve[Evaluate[eqs0/.{(Subscript[ϕ, _])->1.,δk->0,δω->0,(Subscript[Ω, 1])->(10^11),(Subscript[Ω, 2])->2. (10^10),(Subscript[Ω, _])->(10^10),(Subscript[δ, _])->1. (10^7),γt->10.,ΩL->.1}/.atomicdata/.Hertz->1/.Mega->1/.(Subscript[ρ, s1_,s2_])'[t]->0],DMVariables[system]][[1]];

Find the position in the solution list of all of the DM elements that are zero.

delpos=Position[sol,_->0.+0. I]

Find the list of DM elements that correspond to these zero positions.

delvars=Extract[DMVariables[system],delpos];

Create a list of rules to set these DM elements to zero.

delreps=((#->0)&/@delvars);

Find the list of variables corresponding to nonzero DM elements.

vars=Delete[DMVariables[system],delpos]

Remove all of the zero elements from the evolution equations.

TableForm[eqs=Delete[eqs0/.delreps,delpos]]

Initial conditions for the time-dependent case.

inits=Delete[InitialConditions[system,TransitRepopulation[system,1],0],delpos]

Find equations for the steady state by setting time derivatives to zero. We also substitute in numerical values for natural widths in units of (10^6) rad/s.

steadyeqs=eqs/.(Subscript[ϕ, _])->0/.atomicdata/.Hertz->1/.Mega->1/.(Subscript[ρ, s1_,s2_])'[t]->0

Here we find expressions for the fractional absorption (change in electric-field amplitude) and phase shift of the probe light per unit path length and unit atomic density, i.e., {(1/(Subscript[ℰ, 0]))((Subscript[dℰ, 0])/((Subscript[n, 0]) dℓ)),((dϕ)/((Subscript[n, 0]) dℓ))}. These correspond to the imaginary and real parts of the index of refraction, respectively. Since the probe light does not interact with the D2 transition, we set DM elements involving the (Subscript[P, 3/2]) state to zero. We also set the vanishing DM elements found above to zero. We find the expressions for arbitrary ellipticity and then set ellipticity to π/4 after simplification in order to avoid an artificial divide by zero problem.

Observables[system,Energy[1],(Subscript[Ω, 2])/ReducedME[0,{Dipole,1},1],{0,ϵ}][[1]];
obs0=Simplify[%/.DMElementPattern[]/;Label2==2||F1==2->0/.delreps]/.ϵ->π/4

Calculate linear probe absorption with weak probe light on resonance and no pump to find reference level

params={(Subscript[Ω, 1])->0,(Subscript[Ω, 2])->.00001,(Subscript[Ω, 3])->0,(Subscript[Ω, 4])->0,(Subscript[δ, 2])->0,(Subscript[δ, 1])->0,(Subscript[δ, 3])->0,γt->.0001,ΩL->0};
linabs=-obs0[[1]]/.Chop@NSolve[Evaluate[steadyeqs/.params],vars][[1]]/.params

Divide fractional absorption and phase shift by the linear absorption per unit length to find observables per absorption length.

obs={-1,1}obs0/linabs

Plot steady-state results

SetOptions[ListLinePlot,Frame->True,Axes->False,PlotRange->All];

Here we set parameters and plot the fractional absorption and phase shift of the probe light per absorption length, corresponding to imaginary and real parts of the index of refraction, respectively, as a function of probe detuning.

params={(Subscript[Ω, 1])->100.,(Subscript[Ω, 2])->.01,(Subscript[Ω, 3])->10.,(Subscript[Ω, 4])->0.,(Subscript[δ, 2])->d,(Subscript[δ, 1])->0,(Subscript[δ, 3])->0,γt->.1,ΩL->0.};
steadyeqs1=steadyeqs/.params;
table=Table[{d,obs}/.Chop@NSolve[steadyeqs1,vars][[1]]/.params,{d,-200,200,2}];
ListLinePlot[Transpose[{table[[All,1]],table[[All,2,1]]}]]
ListLinePlot[Transpose[{table[[All,1]],table[[All,2,2]]}]]

Plot atomic populations for a given set of parameters.

params={(Subscript[Ω, 1])->10,(Subscript[Ω, 2])->.01,(Subscript[Ω, 3])->1.,(Subscript[Ω, 4])->0.,(Subscript[δ, 2])->0,(Subscript[δ, 1])->0,(Subscript[δ, 3])->0,γt->.1,ΩL->0.};
steadyeqs1=steadyeqs/.params;
sol=Chop@NSolve[steadyeqs1,vars][[1]]/.params;
LevelDiagram[system,Hrwa/.{(Subscript[δ, 1])->-ω1res,(Subscript[δ, 2])->-ω2res,(Subscript[δ, 3])->-ω3res}/.params/.{Energy[0]->0,Energy[1]->2,Energy[2]->3,HyperfineA[_]->.2,HyperfineB[_]->.1,ΩL->0},DensityMatrix[system]/.sol/.DMElementPattern[]->0,ParityOffset->False,Epilog->{Text[Style["( ^2)(Subscript[S, 1/2])",SingleLetterItalics->False],{-2.9,0}],Text[Style["( ^2)(Subscript[P, 1/2])",SingleLetterItalics->False],{-2.9,1.6}],Text[Style["( ^2)(Subscript[P, 3/2])",SingleLetterItalics->False],{-2.9,2.7}]},ImageSize->3.5 72,ImagePadding->{{35,0},{10,10}},PopulationStyle->PointSize[.018]]

Eliminate conjugate DM variables

delvars1=Intersection[vars,Cases[vars,DMElementPattern[]/;!OrderedQ[{State2,State1}]->(Subscript[ρ, State2,State1])[t]]]
delpos1=Flatten[Position[vars,#]&/@delvars1,1]
vars1=Delete[vars,delpos1]

delreps1=delvars1/.(Subscript[ρ, s1_,s2_])[t]->((Subscript[ρ, s1,s2])[t]->Conjugate[(Subscript[ρ, s2,s1])[t]])

eqs1=Delete[eqs,delpos1]/.delreps1

inits1=Delete[inits,delpos1]

Convert equations to XMDS form

Join[atomicdata,
Table[eta[s]==3((Wavelength[s])^2)NaturalWidth[s]/(4 π),{s,Union[Label/@ExcitedStates[system]]}]/.atomicdata,{rt6==N[(Sqrt[6])],rt3==N[(Sqrt[3])],rt2==N[(Sqrt[2])]}
]/.Mega->(10^6)/.Hertz->1/.Centimeter->.01/.Rule->Equal/.{NaturalWidth->g,HyperfineA->ha,HyperfineB->hb,Wavelength->lambda}
StringReplace["//---------------- Constants.txt starts ------------------
"<>StringJoin["const double "<>ToString@CForm[#]<>";\n"&/@%]<>"//---------------- Constants.txt ends   ------------------
",{"=="->"=",Shortest["dr("~~a__~~","~~b__~~")"]:>"dr"<>a<>b<>"_dt",Shortest["r("~~a__~~","~~b__~~")"]:>"r"<>a<>b,"ha("~~a_~~")":>"ha"<>a,"hb("~~a_~~")":>"hb"<>a,"lambda("~~a_~~")":>"lambda"<>a,"eta("~~a_~~")":>"eta"<>a,"g("~~a_~~")":>"g"<>a,"W("~~a_~~")":>"E"<>a,"Wc("~~a_~~")":>"Ec"<>a,"d("~~a_~~")":>"d"<>a}]

Export[ToFileName[NotebookDirectory[],"Constants.txt"],%]

Convert equations to C form

{Label[#],F[#],M[#]}&/@system
labelrep=Table[%[[i]]->ToString@PaddedForm[i,2,NumberPadding->{"0",""},NumberSigns->{"",""}],{i,Length@%}]

eqs1/.{(Subscript[δ, s_])->(Subscript[δ, s])+(Subscript[k, s])v,δω->δω+δk v}/.{NaturalWidth->g,HyperfineA->ha,HyperfineB->hb,(Subscript[ρ, s1_,s2_])'[t]:>dr[s1,s2],(Subscript[ρ, s1_,s2_])[t]:>r[s1,s2],γt->gt,ΩL->WL,(Subscript[δ, i_])->d[i],(Subscript[k, i_])->k[i],(E^(a_.+I (Subscript[ϕ, i_])))->(E^a)W[i]/(Subscript[Ω, i]),(E^(a_.-I (Subscript[ϕ, i_])))->(E^a)Wc[i]/(Subscript[Ω, i])}/.labelrep/.Power[E,b_]->exp[b]/.Conjugate->conj/.Complex[0,a_]->a i/.{1/(Sqrt[6])->1/rt6,1/(Sqrt[3])->1/rt3,1/(Sqrt[2])->1/rt2}
StringReplace[
"//---------------- RbEquations.txt starts ------------------
"<>StringJoin[ToString@CForm[#]<>";\n"&/@%]<>
"//---------------- RbEquations.txt ends  ------------------
",{"=="->"=",Shortest["dr(\""~~a__~~"\",\""~~b__~~"\")"]:>"dr"<>a<>b<>"_dt",Shortest["r(\""~~a__~~"\",\""~~b__~~"\")"]:>"r"<>a<>b,"ha("~~a_~~")":>"ha"<>a,"hb("~~a_~~")":>"hb"<>a,"g("~~a_~~")":>"g"<>a,"W("~~a_~~")":>"E"<>a<>"a","Wc("~~a_~~")":>"conj(E"<>a<>"a)","d("~~a_~~")":>"delta"<>a,"k("~~a_~~")":>"Kvec"<>a}]

Export[ToFileName[NotebookDirectory[],"RbEquations.txt"],%]

List of variables

StringReplace[
"<!-- ############### RbChosenRho.txt starts ################ -->
"<>StringJoin[ToString@CForm[#]<>" "&/@(vars1/.labelrep/.(Subscript[ρ, s1_,s2_])[t_]:>r[s1,s2])]<>
"
<!-- ############### RbChosenRho.txt ends  ################# -->
",{"=="->"=",Shortest["dr(\""~~a__~~"\",\""~~b__~~"\")"]:>"dr"<>a<>b<>"_dt",Shortest["r(\""~~a__~~"\",\""~~b__~~"\")"]:>"r"<>a<>b,"ha("~~a_~~")":>"ha"<>a,"hb("~~a_~~")":>"hb"<>a,"g("~~a_~~")":>"g"<>a,"W("~~a_~~")":>"E"<>a,"Wc("~~a_~~")":>"Ec"<>a,"d("~~a_~~")":>"d"<>a,"k("~~a_~~")":>"k"<>a}]

Export[ToFileName[NotebookDirectory[],"RbChosenRho.txt"],%]

Convert initial conditions to c form.

inits1/.{(Subscript[δ, s_])->(Subscript[δ, s])+(Subscript[k, s])v,δω->δω+δk v}/.{NaturalWidth->g,HyperfineA->ha,HyperfineB->hb,(Subscript[ρ, s1_,s2_])'[t_]:>dr[s1,s2],(Subscript[ρ, s1_,s2_])[t_]:>r[s1,s2],γt->gt,(Subscript[δ, i_])->d[i],(Subscript[k, i_])->k[i],(E^(a_.+I (Subscript[ϕ, i_])))->(E^a)W[i]/(Subscript[Ω, i]),(E^(a_.-I (Subscript[ϕ, i_])))->(E^a)Wc[i]/(Subscript[Ω, i])}/.labelrep/.Power[E,b_]->exp[b]/.Complex[0,a_]->a i/.{1/(Sqrt[6])->1/rt6,1/(Sqrt[3])->1/rt3,1/(Sqrt[2])->1/rt2}
StringReplace[
"//---------------- RbInits.txt starts ------------------
"<>StringJoin[ToString@CForm[#]<>";\n"&/@%]<>
"//---------------- RbInits.txt ends  ------------------
",{"=="->"=",Shortest["dr(\""~~a__~~"\",\""~~b__~~"\")"]:>"dr"<>a<>b<>"_dt",Shortest["r(\""~~a__~~"\",\""~~b__~~"\")"]:>"r"<>a<>b,"ha("~~a_~~")":>"ha"<>a,"hb("~~a_~~")":>"hb"<>a,"g("~~a_~~")":>"g"<>a,"W("~~a_~~")":>"E"<>a,"Wc("~~a_~~")":>"Ec"<>a,"d("~~a_~~")":>"d"<>a,"k("~~a_~~")":>"k"<>a}]

Export[ToFileName[NotebookDirectory[],"RbInits.txt"],%]

Find propagation equations

Here we find the propagation equations.

First, for each field 1--4, extract the list of sublevels that the field is assumed to interact with.

transitioncriteria={Label==0&&F==2||Label==1&&F==1,Label==0&&F==1||Label==1&&F==1,Label==0&&F==1||Label==2&&F==2,Label==0&&F==2||Label==2&&F==2}

subsys=Table[SelectStates[system,transitioncriteria[[k]]],{k,Length[transitioncriteria]}]

The reduced matrix elements corresponding to the transition for each field.

rmes={ReducedME[0,{Dipole,1},1], ReducedME[0,{Dipole,1},1], ReducedME[0,{Dipole,1},2], ReducedME[0,{Dipole,1},2]};

The first-order wave equation is given in terms of the rotating-frame density matrix by
	((∂(Subscript[E, 0]))/(∂z))+(1/c)((∂(Subscript[E, 0]))/(∂t))=((4π i ω n)/c)(∑_(m<n))(Subscript[ρ, n,m])(Subscript[d, m,n])

First we calculate the right-hand side. Here are the covariant spherical components of (∑_(m<n))(Subscript[ρ, n,m])(Subscript[d, m,n]) for each of the four transitions, in terms of the reduced dipole matrix element ||d||.

polarizationcomponents1=Table[Sum[Boole[BranchingRatio[m][n]=!=0&&BranchingRatio[n][m]===0]DMElement[SublevelLabel@n,SublevelLabel@m,t]WignerEckart[m,{Dipole,1},n],{n,s},{m,s}],{s,subsys}]/.delreps

Multiply through by 4π i ω n ||d||, and rewrite ||d((||)^2) in terms of the natural width Γ of the transition.

polarizationcomponents2=Table[Expand[ExpandDipoleRME[subsys[[k]],4π I Energy[Last[rmes[[k]]]] n0 rmes[[k]] polarizationcomponents1[[k]]/.delreps1]],{k,Length[subsys]}]

Make the substitution (Subscript[η, s])=((3(Subscript[λ, s]^2)(Subscript[Γ, s])n)/(4π)).

polarizationcomponents3=polarizationcomponents2/.Energy[s_]->2π/λ[s]/.NaturalWidth[s_]->4 π η[s]/(n0 ((λ[s])^2)3)//Simplify

Cartesian components of the amplitude of each of the 4 fields.

fieldcomponents1=Table[field/.{(Subscript[k, _])->0,t->0,(Subscript[ϕ, _])->0,(Subscript[Ω, s_])/;i≠s->0},{i,4}]

Convert to covariant spherical components and multiply by ||d||.

fieldcomponents2=ToCovariant/@fieldcomponents1 rmes/.(Subscript[Ω, s_])->(Subscript[Ω, s])[z]

Find propagation equations using the spherical components of the polarization and fields. Here we leave out the time-derivative term, which we will add in below.

propeqs=Solve[D[fieldcomponents2,z]==Simplify@polarizationcomponents3,Table[(Subscript[Ω, i])'[z],{i,4}]][[1]]/.Rule->Equal

Convert propagation equations to c form.

propeqs/.{η[s_]->Ndens eta[s],NaturalWidth->g,(Subscript[ρ, s1_,s2_])[t]:>r[s1,s2],(Subscript[Ω, i_])'[z]->dW[i]}/.labelrep/.Complex[0,a_]->a i/.Conjugate->conj/.a:Power[_?NumericQ,_]:>N[a]//Simplify
StringJoin[Table[ToString@CForm[%[[j]]]<>" - Lt[E"<>ToString@j<>"];\n",{j,4}]]
"//---------------- RbPropEquations.txt starts ------------------
"<>StringReplace[%,{"=="->"=","eta("~~a_~~")":>"eta"<>a,Shortest["dr(\""~~a__~~"\",\""~~b__~~"\")"]:>"dr"<>a<>b<>"_dt",Shortest["r(\""~~a__~~"\",\""~~b__~~"\")"]:>"r"<>a<>b,"ha("~~a_~~")":>"ha"<>a,"hb("~~a_~~")":>"hb"<>a,"g("~~a_~~")":>"g"<>a,"dW("~~a_~~")":>"dE"<>a<>"_dz","Wc("~~a_~~")":>"Ec"<>a,"d("~~a_~~")":>"d"<>a}]<>
"//---------------- RbPropEquations.txt ends  ------------------
"
Export[ToFileName[NotebookDirectory[],"RbPropEquations.txt"],%]