Quantum Stochastic Master Equation

Way back when I was a postdoc, the Lindblad equation was where all the action was. People had creaky Fortran codes which they used to generate their results. I was amazed that they worked, being full of strange constants and fudge factors. Here is a Mathematica/Wolfram Language notebook–taken from Mohammed Bahrami at Wolfram Research–with a few minor modifications that solves the Quantum Stochastic Master Equation and plots the results in an easily understandable self-contained form.

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.161403,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[]=

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.137536,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[]=

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[]=

{7.43195,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

In[]:=

​​