Lindblad equation: quantum jumps

Graphics[{
Line[{{0,0},{1,0}}],
Line[{{0,2},{1,2}}],
Text["|1⟩ = |g⟩",{1.3,0}],
Text["|0⟩ = |e⟩",{1.3,2}],
Arrow[{{0.3,0.1},{0.3,1.9}}],
Arrow[{{0.83,1.9},{0.83,0.1}}],
Text["γN",{0.19,1.}],
Text["γ(N+1)",{1.03,1.}]

},BaseStyle->{FontFamily->"Times",20},
PlotRange->{{-0.1,1.6},{-0.3,2.3}}]

Export["/Users/gtlandi/Dropbox/Z QInfo notes/Figures/Open/LindbladQubit_Example1_drawing.pdf",%];

LoadPauliMatrices[]

Clear[ψ0,ρ,γ,f];
ψ0={Cos[θ],Exp[I (π/4)]Sin[θ]}/.θ->π/3;
γ=0.01;
𝒩=5.;
sol=First@NDSolve[
{D[ρ[t],t]==-(I/2)(σz.ρ[t]-ρ[t].σz)+ γ(𝒩+1)Lindbladian[σm][ρ[t]]+γ 𝒩 Lindbladian[σp][ρ[t]],ρ[0]==out[ψ0,ψ0]},ρ,{t,0,tf=300}];

s[ρ_]:=Chop@2{((ρ[[2,1]]+ρ[[1,2]])/2),((ρ[[2,1]]-ρ[[1,2]])/(2I)), ρ[[1,1]]-(1/2)}
Show[BlochSphere[],Graphics3D[{Black,Thickness[0.0066],Line[Chop@Table[s[(ρ/.sol)[t]],{t,0,tf,0.1}]]}]]

(*Export["/Users/gtlandi/Dropbox/Z QInfo notes/Figures/Open/LindbladQubit_Example1.png",%,ImageResolution->200];*)

Clear[ψ0,ρ,γ,f];
ψ0={Cos[(π/8)],Exp[I (π/4)]Sin[(π/8)]};
γ=0.3;
H = (1/2)(σz + γ σx);
𝒩=2.;
sol=First@NDSolve[
{D[ρ[t],t]==-(I/2)(H.ρ[t]-ρ[t].H),ρ[0]==out[ψ0,ψ0]},ρ,{t,0,tf=1000}];

s[ρ_]:=Chop@2{((ρ[[2,1]]+ρ[[1,2]])/2),((ρ[[2,1]]-ρ[[1,2]])/(2I)), ρ[[1,1]]-(1/2)}
Show[BlochSphere[],Graphics3D[{Black,Thickness[0.0066],Line[Chop@Table[s[(ρ/.sol)[t]],{t,0,100,0.1}]]}]]

LoadPauliMatrices[]

ℒ = UnitaryToVec[(Ω/2)σz + (λ/2)σx] + γ(𝒩+1)LindbladToVec[σm]+γ 𝒩 LindbladToVec[σp]//Simplify//mf;

ρ=AnalyticalSteadyState[ℒ]//mf;

sx=Tr[σx.ρ]//Simplify
sy=Tr[σy.ρ]//Simplify
sz=Tr[σz.ρ]//Simplify

Limit[{sx,sy,sz},λ->0]

Limit[{sx,sy,sz},λ->∞]

LoadPauliMatrices[]

ψ0=Cos[(π/8)]{1,0}+Exp[I (π/4)]Sin[(π/8)]{0,1};
δt=0.1;
𝒩 = 1.0;
γ = 0.03;
H = (1/2)σz;
Heff=H-I (γ/2)((𝒩+1)σp.σm+𝒩 σm.σp);
K=MatrixExp[-I Heff δt];
Sm = σm;
Sp = σp;

Pm[x_]:=Chop[γ (𝒩+1) δt Conjugate[x].Sm\[ConjugateTranspose].Sm.x ]
Pp[x_]:=Chop[γ 𝒩 δt Conjugate[x].Sp\[ConjugateTranspose].Sp.x ]
jumpm[x_]:=((Sm.x)/(Sqrt[Conjugate[x].Sm\[ConjugateTranspose].Sm.x]))
jumpp[x_]:=((Sp.x)/(Sqrt[Conjugate[x].Sp\[ConjugateTranspose].Sp.x]))
evol[x_]:=((K.x)/(Sqrt[Conjugate[x].K\[ConjugateTranspose].K.x]))

traj[x_,M_]:=Module[{tr=Table[x,M],r= RandomReal[1,M]},
Do[

tr[[j]]=tr[[j-1]];

Which[
r[[j]]>(Pm[tr[[j]]]+Pp[tr[[j]]]),tr[[j]]=evol[tr[[j]]],
r[[j]]≤Pp[tr[[j]]], tr[[j]]=jumpp[tr[[j]]],
(Pp[tr[[j]]]<r[[j]]≤(Pm[tr[[j]]]+Pp[tr[[j]]])),tr[[j]]=jumpm[tr[[j]]]
];


,{j,2,M}];

tr
]

trajExample=traj[ψ0,700];
ListLinePlot[((Abs[trajExample[[All,2]]])^2)]

Show[BlochSphere[],Graphics3D[{Black,Thickness[0.0045],Line[Chop@Table[s[out[ψ,ψ]],{ψ,trajExample}]]}],ImageSize->400]

ntraj = 100;
trajAve=(1/ntraj)Total@Table[out[#,#]&/@ traj[ψ0,700],{i,ntraj}];
s[ρ_]:=Chop@2{((ρ[[2,1]]+ρ[[1,2]])/2),((ρ[[2,1]]-ρ[[1,2]])/(2I)), ρ[[1,1]]-(1/2)}

Show[BlochSphere[],Graphics3D[{Black,Thickness[0.0045],Line[Chop@Table[s@rho,{rho,trajAve}]]}],ImageSize->700]

Tr[σp.σm.{{f,0},{0,1-f}}]