ClearAll["Global`*"]
(*Computing EIT susceptibility for Zeeman degenerate atoms*)
(* author: B. Gouraud*)
(*Equations from my thesis, section 3.1.3 *)
(*Update and important bug solved 2015-08-12, updating sign conventions on 2015-10-02, 2015-11-26: a factor 1/2 was missing in the B vector*)
(*Usage: Define the atom and polarizations in the PARAMETERS part, execute the rest of the code.*)

(*//////////////////////////////////*)
(*////////PARAMETERS////////////////*)
(*//////////////////////////////////*)

Fg=0;  (*Angular momentum of the g,s and e states*)
Fs=1;
Fe=1;

ϵpp={0,0,1}; (*Polarization of the probe. Will be normalized below. Format: {σ+,π,σ-}*)
ϵcc={0,1,0}; (*Polarization of the control*)

(*CHOSE HERE which Zeeman sublevels are populated*)
(*ρ0 defines which Zeeman sublevels are initially populated*)
ρ0[μ_,ν_]:=Piecewise[{{1/(2*Fg+1),μ==ν}},0]; (*all g states equally populated: unpolarized atom*)

(*ρ0[μ_,ν_]:=Piecewise[{{1/(2*Fg+1),μ==ν&&μ==2}},0];*)
(* |g,+2> state only is populated: atom pumped if mf=+2 *)

(*The following should not necessarily be defined, if they are notn Mathermatica will give explicite expressions.*)
(*Γ=1;*) (*Γ=1, all frequencies expressed in units of Γ*)
(*Δ=0;*)(*Control beam detuning*)
Ωp=-3; (*Rabi frequencies of the probe and control*)(*Since Ω=-d.E, Ω is chosen equal to -1 to obtain the susceptibility. Times a factor 3=2L+1 where L is the photonic momentum.*)
(*Ωc=0;*)

(*//////////////////////////////////*)
(*//////END PARAMETERS//////////////*)
(*//////////////////////////////////*)

J=-(δ-I*Γ/2);
c=((Ωc^2)/(4(δ-Δ)));

NormS[ϵ_]:=Sqrt[Abs[ϵ[[1]]]^2+Abs[ϵ[[2]]]^2+Abs[ϵ[[3]]]^2];
ϵp=ϵpp/NormS[ϵpp]; (*Normalized polarizations*)
ϵc=ϵcc/NormS[ϵcc];
s[ϵ_,q_]:=ϵ[[q+2]];(* component q (q in {-1,0,1}) of a vector in expressed in the spherical basis *)

Clebsch=Piecewise[{{ClebschGordan[{FE,i2},{1,i1-i2},{FG,i1}],-FG≤i1≤FG&&-FE≤i2≤FE&&-1≤i2-i1≤1 }},0];
cC[μ_,ν_]:=Clebsch/.{FG->Fs,FE->Fe,i1->μ,i2->ν}(*Clebsch-Gordan coefficient for the control field. Noted Cμ,ν in the thesis.*)
cP[μ_,ν_]:=Clebsch/.{FG->Fg,FE->Fe,i1->μ,i2->ν}(*Clebsch-Gordan coefficient for the probe field. Noted Pμ,ν in the thesis.*)

(*Defining the M coefficients (intermediate steps: f and Mf).*)
f[ν_,q_]:=Piecewise[{
{Sum[cC[ν-i,ν]^2*Abs[s[ϵc,-i]]^2,{i,-1,1}],q==0},
{cC[ν+1,ν]*cC[ν+1,ν+1]*Conjugate[s[ϵc,1]]*s[ϵc,0]+cC[ν,ν]*cC[ν,ν+1]*Conjugate[s[ϵc,0]]*s[ϵc,-1],q==1},
{cC[ν,ν]*cC[ν,ν-1]*Conjugate[s[ϵc,0]]*s[ϵc,1]+cC[ν-1,ν]*cC[ν-1,ν-1]*Conjugate[s[ϵc,-1]]*s[ϵc,0],q==-1},
{cC[ν+1,ν]*cC[ν+1,ν+2]*s[ϵc,-1]*Conjugate[s[ϵc,1]],q==2},
{cC[ν-1,ν]*cC[ν-1,ν-2]*s[ϵc,1]*Conjugate[s[ϵc,-1]],q==-2}}
,0
]

Mf:=SparseArray[{{k_,l_}/;Abs[k-l]≤2->f[k-Fe-1,l-k]},{2*Fe+1,2*Fe+1}];
M:=c*Mf+J*IdentityMatrix[2*Fe+1]; (*Matrix to be inversed*)

(*Defining the B vector.*)
B[μ_,ν_]:=-Ωp/2*Sum[cP[ν-q,ν]*Conjugate[s[ϵp,-q]]*ρ0[μ,ν-q],{q,-1,1}];
Bv:=Table[B[μ,k-Fe-1],{k,2*Fe+1}];

(*Solve the linear system*)
inv:=LinearSolve[M] (*Inverse the M matrix*)
sols=Table[inv[Bv],{μ,-Fg,Fg}];

Print["EIT susceptibility"]
su=Table[Sum[Piecewise[{{cP[m,m-q]*Conjugate[(sols[[m+1+Fg]])[[m-q+Fe+1]]],1≤m-q+1+Fe≤2*Fe+1}},0],{m,-Fg,Fg}],{q,-1,1}]
(*Piecewise necessary?*)(*The EIT susceptibility*)
(*The susceptibility projected on the probe polarization.*)
suProjected=Simplify[Sum[su[[i]]*Conjugate[ϵp[[i]]],{i,1,3}],{Γ,Ωc,δ,Δ}∈Reals](*A dot Conjugate[B]*)

constants={Ωc->1,Δ->0.2,Γ->1}
(*Plot all imaginary and real parts of the susceptibility.*)
GraphicsRow[Plot[{Im[su[[#]]]/.constants,Re[su][[#]]/.constants},{δ,-2,2},PlotRange->All]&/@{1,2,3},ImageSize->Large]
Plot[{Im[suProjected]/.constants,Re[suProjected/.constants]},{δ,-2,2}]