The Quantum Stochastic Master Equation (SME) describes the evolution of an open quantum system interacting with its environment in a probabilistic manner. It extends the Lindblad master equation by incorporating stochastic noise, often modeled using quantum Wiener or Poisson processes. QSMEs are crucial for understanding decoherence, quantum measurement dynamics, and feedback control in quantum information processing. Implementation is used from: Pierre Rouchon, Jason F. Ralph (2015), Efficient Quantum Filtering for Quantum Feedback Control. Phys. Rev. A91, 012118. https://doi.org/10.1103/PhysRevA.91.012118
Dynamical equation of the quantum system:
ρ=t-i[H,ρ]+
∑
i
Γ
i
V
i
ρ
†
V
i
-
1
2

†
V
i
V
i
,ρ+
∑
i
γ
i
L
i
ρ
†
L
i
-
1
2

†
L
i
L
i
,ρ+
∑
i
η
i

w
i
γ
i

L
i
ρ+ρ
†
L
i
-Tr
L
i
ρ+ρ
†
L
i
ρ
There are jump/Lindblad operators due to the environment, denoted by
V
i
, that contribute only as decoherence term in the dynamical equation (deterministic) with the decoherence rate
Γ
i
.
There are jump/Lindblad operators due to the monitoring of the system (with a corresponding readout), denoted by
L
i
, that contribute as decoherence terms and also stochastic terms in the dynamical equation, with the decoherence rate of
γ
i
and the measurement channel efficiency of
η
i
. Their corresponding output signal reads as

R
j
=
η
i
tTr
L
i
ρ+
†
L
i
+
w
i
We have implemented the efficient numerical scheme proposed in this paper: https://arxiv.org/pdf/1410.5345​
The main function for simulation is the following one:
[
ρ
0
,H,L,η,V,δt,
t
f
]
The function for manual evolution, with ρ the initial state,
H
the Hamiltonian,
L
the list of jump operators being monitored (and the evolution is conditioned upon their corresponding readout currents), η the readout efficiency,
V
the list of jump operator not being monitors (e.g. environmental noise; by default it is None), δt the time increment and
t
f
the final time.
​
We call it manual evolution because we do not use numerical features of Mathematica such as NDSolve or ItoProcess. The main reason for that is we were not able to using those functionalities while preserving important features of a Bloch sphere dynamics.
The readout currents are sowed in .
Below examples are taken from Prof. Gabriel T. Landi’s melt notebook.

ℒ[ρ]=-i[
Ωσ
x
,ρ]+γD[
σ
z
]
and γ >> Ω

Hamiltonian, jump operators and damping rates:
ℒ[ρ]=-i[
Ωσ
x
,ρ]+γD[
σ
z
]
In[]:=
SeedRandom[1];​​Ω=1.0;H=ΩQuantumOperator["X"];
Jump operators and damping rates (note they are List):
In[]:=
γs={2.0};​​Ls={QuantumOperator["Z"]};
Time increment and final time:
In[]:=
δt=0.005;​​tf=20.;
Initial state:
In[]:=
ρ0=QuantumState["0"];
Solve the Lindblad master equation for the density matrix:
∂
t
ρ=-i[H,ρ]+
∑
i
γ
i
L
i
ρ
†
L
i
-
1
2

†
L
i
L
i
,ρ
In[]:=
ρt=QuantumEvolve[H,Ls->γs,ρ0,{t,0,tf}]
Out[]=
In[]:=
Plot[Evaluate@Re@ρt[t]["BlochVector"],{t,0,tf},PlotRange->All]
Out[]=
In[]:=
trajectory=ρ0,H,
γs
Ls,δt,tf;//AbsoluteTiming
Out[]=
{0.115558,Null}
Check if any un-physical state:
In[]:=
PositiveSemidefiniteMatrixQ/@trajectory//Tally
Out[]=
{{True,4001}}
In[]:=
bloch=Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&/@trajectory;
In[]:=
ListLinePlot[Transpose[bloch]]
Out[]=
1000
2000
3000
4000
-1.0
-0.5
0.5
1.0
For each trajectory, find Bloch vector:
In[]:=
Show[QuantumState["UniformMixture"]["BlochPlot"],ListLinePlot3D[bloch]]
Out[]=
Compare stochastic averaged trajectory with Lindbladian for 200 realizations (they must match):
In[]:=
Module{trjs,bloch,average,tf},​​tf=5;​​trjs=Tableρ0,H,
γs
Ls,δt,tf,200;​​bloch=Map[Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&]/@trjs;​​average=Table[Re[Evaluate@ρt[t]["BlochVector"]],{t,0,tf,δt}]//Transpose;​​Table[ListLinePlot[{average[[i]],Mean[(Transpose/@bloch)[[All,i,All]]]},ImageSize->200],{i,3}]
Out[]=

,
,


ℒ[ρ]=-i[
Ωσ
x
,ρ]+γD[
σ
z
]
and γ << Ω

Hamiltonian, jump operators and damping rates:
ℒ[ρ]=-i[
Ωσ
x
,ρ]+γD[
σ
z
]
In[]:=
SeedRandom[1];​​Ω=1.0;H=ΩQuantumOperator["X"];
Jump operators and damping rates (note they are List):
In[]:=
γs={.1};​​Ls={QuantumOperator["Z"]};
Time increment and final time:
In[]:=
δt=0.005;​​tf=20.;
Initial state:
In[]:=
ρ0=QuantumState["0"];
Solve the Lindblad master equation for the density matrix:
∂
t
ρ=-i[H,ρ]+
∑
i
γ
i
L
i
ρ
†
L
i
-
1
2

†
L
i
L
i
,ρ
In[]:=
ρt=QuantumEvolve[H,Ls->γs,ρ0,{t,0,tf}]
Out[]=
Plot the evolution:
In[]:=
Plot[Evaluate@Re@ρt[t]["BlochVector"],{t,0,tf},PlotRange->All]
Out[]=
In[]:=
trajectory=ρ0,H,
γs
Ls,δt,tf;//AbsoluteTiming
Out[]=
{0.116894,Null}
Check if any un-physical state:
In[]:=
PositiveSemidefiniteMatrixQ/@trajectory//Tally
Out[]=
{{True,4001}}
Bloch vector evolution
In[]:=
bloch=Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&/@trajectory;
Visualize them:
In[]:=
ListLinePlot[Transpose[bloch]]
Out[]=
1000
2000
3000
4000
-1.0
-0.5
0.5
1.0
Compare stochastic averaged trajectory with Lindbladian for 200 realizations (they must match) :
In[]:=
Module{trjs,bloch,average,tf},​​tf=20;​​trjs=Tableρ0,H,
γs
Ls,δt,tf,200;​​bloch=Map[Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&]/@trjs;​​average=Table[Re[Evaluate@ρt[t]["BlochVector"]],{t,0,tf,δt}]//Transpose;​​Table[ListLinePlot[{average[[i]],Mean[(Transpose/@bloch)[[All,i,All]]]},ImageSize->200],{i,3}]
Out[]=

,
,


Diffusion for spontaneous emission

Hamiltonian, jump operators and damping rates:
In[]:=
SeedRandom[1];​​Δ=1.;Ω=2.;​​H=1/2(ΩQuantumOperator["X"]+ΔQuantumOperator["Z"]);
Jump operators and damping rates (note they are List):
In[]:=
γs={.2};​​Ls={QuantumOperator["J-"]};
Time increment and final time:
In[]:=
δt=0.01;​​tf=10.;
Initial state:
In[]:=
ρ0=QuantumState["+"];
Solve the Lindblad master equation for the density matrix:
∂
t
ρ=-i[H,ρ]+
∑
i
γ
i
L
i
ρ
†
L
i
-
1
2

†
L
i
L
i
,ρ
In[]:=
ρt=QuantumEvolve[H,Ls->γs,ρ0,{t,0,tf}]
Out[]=
Plot the evolution:
In[]:=
Plot[Evaluate@Re@ρt[t]["BlochVector"],{t,0,tf},PlotRange->All]
Out[]=
In[]:=
trajectory=Tableρ0,H,
γs
Ls,δt,tf,200;//AbsoluteTiming
Out[]=
{6.77819,Null}
Check if any un-physical state:
In[]:=
Map[PositiveSemidefiniteMatrixQ]/@trajectory//Flatten//Tally
Out[]=
{{True,200200}}
In[]:=
bloch=Map[Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&]/@trajectory;
In[]:=
ListLinePlot[Transpose[bloch[[1]]]]
Out[]=
Compare stochastic averaged trajectory with Lindbladian for 200 realizations (they must match) :
In[]:=
Module[{average},​​average=Table[Re[Evaluate@ρt[t]["BlochVector"]],{t,0,tf,δt}]//Transpose;​​Table[ListLinePlot[{average[[i]],Mean[(Transpose/@bloch)[[All,i,All]]]},ImageSize->200],{i,3}]]
Out[]=

,
,


Reproducing Fig. 4.6 of Wiseman and Milburn

Exp #1

Hamiltonian, jump operators and damping rates:
In[]:=
SeedRandom[1];​​Δ=0;Ω=3.;​​H=1/2(ΩQuantumOperator["Y"]+ΔQuantumOperator["Z"]);
Jump operators and damping rates (note they are List):
In[]:=
γs={1};​​Ls={QuantumOperator["J-"]};
Time increment and final time:
In[]:=
δt=0.01;​​tf=10.;
Initial state:
In[]:=
ρ0=QuantumState["+"];
Solve the Lindblad master equation for the density matrix:
∂
t
ρ=-i[H,ρ]+
∑
i
γ
i
L
i
ρ
†
L
i
-
1
2

†
L
i
L
i
,ρ
In[]:=
ρt=QuantumEvolve[H,Ls->γs,ρ0,{t,0,tf}]
Out[]=
Plot the evolution:
In[]:=
Plot[Evaluate@Re@ρt[t]["BlochVector"],{t,0,tf},PlotRange->All,PlotLegends->{"x","y","z"}]
Out[]=
x
y
z
Manual simulation:
In[]:=
trajectory=ρ0,H,
γs
Ls,δt,tf;
Check if any un-physical state:
In[]:=
PositiveSemidefiniteMatrixQ/@trajectory//Tally
Out[]=
{{True,1001}}
In[]:=
bloch=Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&/@trajectory;
In[]:=
ListLinePlot[Transpose[bloch]]
Out[]=
In[]:=
Module{trjs,bloch,average,tf},​​tf=10;​​trjs=Tableρ0,H,
γs
Ls,δt,tf,500;​​bloch=Map[Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&]/@trjs;​​average=Table[Re[Evaluate@ρt[t]["BlochVector"]],{t,0,tf,δt}]//Transpose;​​Table[ListLinePlot[{average[[i]],Mean[(Transpose/@bloch)[[All,i,All]]]},ImageSize->200,PlotRange->All],{i,3}]
Out[]=

,
,

In[]:=
ListLinePlot[Transpose[{ArcTan[#1,#2],#3}&@@@bloch],PlotLegends->{"ϕ","cos[θ]"}]
Out[]=
ϕ
cos[θ]
In[]:=
Show[QuantumState["UniformMixture"]["BlochPlot"],ListLinePlot3D[bloch]]
Out[]=

Exp #2

Hamiltonian, jump operators and damping rates:
In[]:=
SeedRandom[1];​​Δ=0;Ω=3.;​​H=1/2(ΩQuantumOperator["Y"]+ΔQuantumOperator["Z"]);
Jump operators and damping rates (note they are List):
In[]:=
γs={1};​​Ls={Exp[π/2]QuantumOperator["J-"]};
Time increment and final time:
In[]:=
δt=0.01;​​tf=10.;
Initial state:
In[]:=
ρ0=QuantumState["+"];
Solve the Lindblad master equation for the density matrix:
∂
t
ρ=-i[H,ρ]+
∑
i
γ
i
L
i
ρ
†
L
i
-
1
2

†
L
i
L
i
,ρ
In[]:=
ρt=QuantumEvolve[H,Ls->γs,ρ0,{t,0,tf}]
Out[]=
Plot the evolution:
In[]:=
Plot[Evaluate@Re@ρt[t]["BlochVector"],{t,0,tf},PlotRange->All]
Out[]=
Manual
In[]:=
trajectory=ρ0,H,
γs
Ls,δt,tf;
Check if any un-physical state:
In[]:=
PositiveSemidefiniteMatrixQ/@trajectory//Tally
Out[]=
{{True,1001}}
In[]:=
bloch=Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&/@trajectory;
In[]:=
ListLinePlot[Transpose[bloch]]
Out[]=
In[]:=
Module{trjs,bloch,average,tf},​​tf=10;​​trjs=Tableρ0,H,
γs
Ls,δt,tf,500;​​bloch=Map[Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&]/@trjs;​​average=Table[Re[Evaluate@ρt[t]["BlochVector"]],{t,0,tf,δt}]//Transpose;​​Table[ListLinePlot[{average[[i]],Mean[(Transpose/@bloch)[[All,i,All]]]},ImageSize->300,PlotRange->All],{i,3}]
Out[]=

,
,

In[]:=
ListLinePlot[Transpose[{ArcTan[#1,#2],#3}&@@@bloch],PlotLegends->{"ϕ","cos[θ]"}]
Out[]=
ϕ
cos[θ]
In[]:=
Show[QuantumState["UniformMixture"]["BlochPlot"],ListLinePlot3D[bloch]]
Out[]=

Two-point function for the stochastic homodyne current

Hamiltonian, jump operators and damping rates:
In[]:=
SeedRandom[1];​​Δ=0;Ω=3.5;​​H=1/2(ΩQuantumOperator["Y"]+ΔQuantumOperator["Z"]);
Jump operators and damping rates (note they are List):
In[]:=
γs={.7};​​Ls={QuantumOperator["J-"]};
Time increment and final time:
In[]:=
δt=0.01;​​tf=4000.;
Initial state:
In[]:=
ρ0=QuantumState["0"];
Generate trajectory and output currents:
In[]:=
{trajectory,{Is}}=Reap@ρ0,H,
γs
Ls,δt,tf;
Check if any un-physical state:
In[]:=
PositiveSemidefiniteMatrixQ/@trajectory//Tally
Out[]=
{{True,400001}}
In[]:=
ListLinePlot[Transpose[bloch]]
Out[]=
Since only one jump operator being monitored, only one current:
In[]:=
I1=Re@Transpose[Is][[1]];
Two point correlations:
In[]:=
ℱ1=TwoPointCorrelation[Chop@Flatten@Is,Floor[8/δt],δt,4];​​ℱ2=TwoPointCorrelation[butterworthFilter[Chop@Flatten@Is,50,δt],Round[8/δt],δt,4];
In[]:=
ListLinePlot[{ℱ1,ℱ2}]
Out[]=

Initialization

In[]:=

Install quantum paclet

In[]:=
PacletInstall["https://wolfr.am/DevWQCF",ForceVersionInstall -> True]​​<<Wolfram`QuantumFramework`
Out[]=
PacletObject
Name: Wolfram/QuantumFramework
Version: 1.4.9

In[]:=

Manual evolution

In[]:=
ClearAll[,ℛ]
Implementation of https://arxiv.org/pdf/1410.5345 to ensure positivity of ρ
In[]:=
ℛ[ρ_,L_,δt_,δw_,η_]:=MapThread[
#3
Tr[(#1+ConjugateTranspose[#1]).ρ]δt+#2&,{L,δw,η}]
In[]:=
ClearAll[]​​[ρ_,H_,L_,η_:None,V_:None,δt_,tf_]:=​​Module{δw,Lmat,Vmat,Hmat,ρ0,Heff,ηm},​​(*η are efficiencies; if not given they are all 1*)​​ηm=If[MatchQ[η,None],ConstantArray[1,Length[L]],η];​​(*turning initial operators or states into matrices*)​​ρ0=ρ["Computational"]["DensityMatrix"];​​Lmat=#["Computational"]["Matrix"]&/@L;​​Vmat=#["Computational"]["Matrix"]&/@V;​​Hmat=H["Computational"]["Matrix"];​​(*given Eq.7, above, this is the term containing everything that is not random*)​​Heff=IdentityMatrix[2]-δt(I Hmat+.5Sum[ConjugateTranspose[l].l,{l,Lmat}]+If[MatchQ[V,None],0,.5Sum[ConjugateTranspose[l].l,{l,Vmat}]]+.5Total[ηm Table[l.l,{l,Lmat}]]);​​(*for each L element, generate a random process*)​​δw=RandomVariateNormalDistribution0,
δt
,{Length[L],Floor[tf/δt]};​​FoldList​​Module{M,st,Υ},​​(*for each step, create the stochastic/random terms to be added to M in Eq.7*)​​Υ=Total[
ηm
Lmat Sow@ℛ[#1,Lmat,δt,#2,ηm]];​​M=Heff+.5Υ.Υ+Υ;​​(*evolve*)​​st=M.#1.ConjugateTranspose[M]+If[MatchQ[V,None],0,δt Sum[l.#1.ConjugateTranspose[l],{l,Vmat}]]+δt Sum[(1-ηm[[l]])Lmat[[l]].#1.ConjugateTranspose[Lmat[[l]]],{l,Range[Length@Lmat]}];​​stTr[st]&,ρ0,Transpose[δw]

Two point correlation and Butterworth filter

Gabriel’ s function :
In[]:=
Clear[TwoPointCorrelation];​​TwoPointCorrelation[data_,hmax_,dt_:None,steps_:1]:=Module{=Length@data,ave = Mean[data],vec1,M,corr},​​vec1 = data1;;-hmax;​​M = Table[data〚1+i;;-hmax+i〛,{i,1,hmax,steps}];​​corr=
1

(M.vec1-
2
ave
);​​If[dt===None,corr,{dt Range[1, hmax,steps],corr}]​​
https://mathematica.stackexchange.com/questions/46500/data-filtering-using-butterworth-filter
In[]:=
Clear[butterworthFilter];​​butterworthFilter[data_,ωc_,dt_:1,order_:2]:=Module[{f,b},​​f=RecurrenceFilter[ToDiscreteTimeModel[ButterworthFilterModel[{"Lowpass",order,ωc }],dt],data];​​b=RecurrenceFilter[ToDiscreteTimeModel[ButterworthFilterModel[{"Lowpass",order,ωc }],dt],Reverse[data]];​​(f+Reverse[b])/2​​]

CITE THIS NOTEBOOK

Quantum stochastic master equation (SME)​
by Mohammad Bahrami​
Wolfram Community, STAFF PICKS, January 31, 2025
​https://community.wolfram.com/groups/-/m/t/3368265