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. A 91, 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[]:=
MatchQ[{1},None]
Out[]=
False
trajectory=ρ0,H,
γs
Ls,δt,tf;//AbsoluteTiming
Out[]=
{10.7171,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[]=

ℒ[ρ]=-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[]:=
average=Table[Evaluate@Re[ρt[t]["BlochVector"]],{t,0,tf,δt}]//Transpose;
In[]:=
trajectory=ρ0,H,
γs
Ls,δt,tf;//AbsoluteTiming
Out[]=
{0.116607,Null}
Check if any un-physical state:
In[]:=
PositiveSemidefiniteMatrixQ/@trajectory//Tally
Out[]=
{{True,4001}}
In[]:=
bloch=Table[Re@Tr[#.PauliMatrix[i]],{i,3}]&/@trajectory;
Visualize them:
In[]:=
Table[ListLinePlot[Transpose[{Range[0,tf,δt],#}]&/@{average[[i]],Transpose@bloch[[All,i]]}],{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:
Initial state:
Plot the evolution:
Check if any un-physical state:

Reproducing Fig. 4.6 of Wiseman and Milburn

Exp #1

Hamiltonian, jump operators and damping rates:
Jump operators and damping rates (note they are List):
Time increment and final time:
Initial state:
Plot the evolution:
Manual simulation:
Check if any un-physical state:

Exp #2

Hamiltonian, jump operators and damping rates:
Jump operators and damping rates (note they are List):
Time increment and final time:
Initial state:
Plot the evolution:
Manual
Check if any un-physical state:

Two-point function for the stochastic homodyne current

Hamiltonian, jump operators and damping rates:
Jump operators and damping rates (note they are List):
Time increment and final time:
Initial state:
Generate trajectory and output currents:
Check if any un-physical state:
Since only one jump operator being monitored, only one current:
Two point correlations:

Initialization

In[]:=

Install quantum paclet

Manual evolution

Implementation of https://arxiv.org/pdf/1410.5345​
to ensure positivity of ρ

Two point correlation and Butterworth filter

Gabriel’ s function :
https://mathematica.stackexchange.com/questions/46500/data-filtering-using-butterworth-filter

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