This document shows how to simulate a continuously monitored qubit with
NDSolve
. We have implemented this approach for a model that follows four quantities: three numbers that describe the qubit state (i.e., the components of the Bloch vector) and the detector's accumulated readout. The dynamical equations for these four quantities are usually written as a nonlinear stochastic differential equation in Itô form. The obstacle is that
NDSolve
cannot use noise as an ordinary function of time, or, more generally, any variable that must be sampled from a distribution. Following the essence of the Wong–Zakai theorem, we therefore convert the original Itô equation to its equivalent Stratonovich form, which is suited to smooth approximations of the noise, and replace the rough Wiener path with a continuous piecewise-cubic path.
NDSolve
can then solve the resulting random ordinary differential equation. As the noise grid is refined, this regularized equation approaches the original stochastic model under the conditions of the Wong–Zakai theorem. Symbolic checks recover the standard master-equation evolution exactly, while finite-sample simulations compare the mean and spread with independent references. These numerical comparisons support the method but do not by themselves prove convergence.

The trajectory equations: the state and the readout in one SDE

A continuously monitored quantum state can be understood schematically as
d
ρ
c
=ℒ(
ρ
c
)dt+ℳ(
ρ
c
)dW
: the drift
ℒ(
ρ
c
)dt
describes deterministic Hamiltonian evolution, relaxation, dephasing, and the unconditional disturbance caused by measurement, while
ℳ(
ρ
c
)dW
is the stochastic, record-dependent measurement update. The detector record obeys
dQ=s(
ρ
c
)dt+dW
, where
s(
ρ
c
)
is the expected signal, so
dW=dQ-s(
ρ
c
)dt
is the innovation (the unpredictable part of the observed result) and the same Wiener increment updates both the record and the state. Because the stochastic average
𝔼[dW]=0
and
d
2
W
=dt
,
dW
is of order
dt
; thus
d
ρ
c
is not an ordinary time derivative but a deterministic prediction plus a random correction conditioned on what the detector reports. Averaging over all possible records removes the stochastic term and recovers the master equation.
Continuous homodyne monitoring of a driven qubit produces two coupled equations driven by the same noise: one for the observer's running estimate of the Bloch vector
{x,y,z=

σ
x
,

σ
y
,

σ
z

, in the convention

σ
z
=|g〉〈g|-|e〉〈e|
that puts the ground state at
z=+1
(some literature uses the opposite sign), and one for the readout
Q
the detector accumulates while measuring

σ
z
. Stacking the readout as a fourth component, the pair is a single scalar-noise Itô stochastic differential equation, a drift

a
times
dt
plus a diffusion

b
times one Wiener increment
dW
:
d
x
y
z
Q
=
-
eff
Γ
2
x
-
eff
Γ
2
y-
Ω
x
z
γ
1
(1-z)+
Ω
x
y
Γ
CI
z

a
dt+
-
Γ
CI
xz-
Γ
BA
y
-
Γ
CI
yz+
Γ
BA
x
Γ
CI
(1-
2
z
)
1

b
dW.
This is the target: the process whose mean we will check against the Lindblad master equation and whose spread we will check against a native integrator. The top three rows, the Bloch state, close on their own, nothing in them depends on
Q
; the fourth row is the record, an integral whose drift
Γ
CI
z
carries the

σ
z
signal and whose noise is the same
dW
that kicks the state.
The six parameters have units of inverse time:
Γ
CI
is the informative (coherent-information) rate, which sharpens the conditional state toward a

σ
z
eigenstate;
Γ
BA
is the backaction rate of the orthogonal homodyne quadrature, which rotates the transverse phase;
Γ
d
is the measurement-induced ensemble dephasing;
γ
ϕ
is intrinsic pure dephasing;
γ
1
=1/
T
1
is relaxation; and
Ω
x
is the Rabi frequency. One combination recurs, the effective transverse decay
eff
Γ
2
=
γ
1
2+
γ
ϕ
+
Γ
d
. When
Γ
d
>0
, the detector efficiency is
η=(
Γ
CI
+
Γ
BA
)/(2
Γ
d
)
, with a quantum-limited detector satisfying
0≤η≤1
. The benchmark ensembles below use
η=1
; one deliberately unphysical example violates the admissibility condition. Numerically, time and all rates are nondimensionalized by a chosen reference rate (the representative example takes that reference to be
γ
1
). Throughout, the code takes a rate vector
r={
Γ
CI
,
Γ
BA
,
Γ
d
,
γ
ϕ
,
γ
1
,
Ω
x
}
.
For a single scalar noise, the Itô equation above is identical to a Stratonovich equation with the same diffusion and a drift lowered by one exact term (the Itô-Stratonovich conversion).
d
x
y
z
Q
=
Γ
BA
+
Γ
CI
2
-
eff
Γ
2
-
Γ
CI
2
z
x-
Γ
BA
Γ
CI
yz
Γ
BA
+
Γ
CI
2
-
eff
Γ
2
-
Γ
CI
2
z
y+
Γ
BA
Γ
CI
xz-
Ω
x
z
γ
1
(1-z)+
Γ
CI
z(1-
2
z
)+
Ω
x
y
Γ
CI
z
dt+
-
Γ
CI
xz-
Γ
BA
y
-
Γ
CI
yz+
Γ
BA
x
Γ
CI
(1-
2
z
)
1
∘dW.
We feed the Stratonovich form to the random ODE because continuous, piecewise-smooth approximations of a one-dimensional Wiener path converge under the usual Wong-Zakai hypotheses to the Stratonovich solution. The fact that the driving noise is scalar matters: with several noncommuting noise channels, the limiting equation can depend on how the missing iterated integrals are approximated. At any nonzero grid spacing the ODE is only a regularized model; the equality with the target SDE is a limit statement, not an identity at finite resolution. The readout row is the same in both calculi because its diffusion coefficient is the constant one, so its conversion correction vanishes. After nondimensionalizing time,
Q
is dimensionless; before nondimensionalization it carries the same units as
W
, namely
time
.

Encoding the Stratonovich vector field

Encode the Stratonovich drift, all four rows and name the entries of the rate vector
r={
Γ
CI
,
Γ
BA
,
Γ
d
,
γ
ϕ
,
γ
1
,
Ω
x
}
:
In[]:=
ClearAll[driftStrat];​​driftStrat[{x_,y_,z_,q_},r_]:=​​With[{ΓCI=r[[1]],ΓBA=r[[2]],Γd=r[[3]],​​γϕ=r[[4]],γ1=r[[5]],Ωx=r[[6]]},​​With[{Γ2=γ1/2+γϕ+Γd},​​{((ΓBA+ΓCI)/2-Γ2-ΓCIz^2)x-Sqrt[ΓBAΓCI]yz,​​((ΓBA+ΓCI)/2-Γ2-ΓCIz^2)y+Sqrt[ΓBAΓCI]xz-Ωxz,​​γ1(1-z)+ΓCIz(1-z^2)+Ωxy,​​Sqrt[ΓCI]z}]];
Encode the Stratonovich diffusion, the same

b
as the Itô form:
In[]:=
ClearAll[diffStrat];​​diffStrat[{x_,y_,z_,q_},r_]:=​​With[{ΓCI=r[[1]],ΓBA=r[[2]]},​​{-Sqrt[ΓCI]xz-Sqrt[ΓBA]y,​​-Sqrt[ΓCI]yz+Sqrt[ΓBA]x,​​Sqrt[ΓCI](1-z^2),​​1}];
To read these back as formulas, write the six rates as symbols, using string subscripts so the labels stay
Γ
CI
,
Γ
BA
,…
and never collide with the Bloch
x
or fuse into a
C·I
product. The readbacks and the symbolic conversion below use
x,y,z,q,Q,t
, a formal noise
w
, the starting values, and the rate symbols themselves as unassigned variables, so clear them first:
In[]:=
ClearAll[x,y,z,q,Q,t,w,x0,y0,z0,q0,Γ,γ,Ω];​​ratesSym={Subscript[Γ,"CI"],Subscript[Γ,"BA"],​​Subscript[Γ,"d"],Subscript[γ,"ϕ"],​​Subscript[γ,"1"],Subscript[Ω,"x"]};
Verify that the encoded drift is the Stratonovich drift written above, by reading it back on the symbolic rates:
In[]:=
driftStrat[{x,y,z,Q},ratesSym]//MatrixForm
Out[]//MatrixForm=
-yz
Γ
BA
Γ
CI
+x-
γ
1
2
-
γ
ϕ
-
2
z
Γ
CI
+
1
2
(
Γ
BA
+
Γ
CI
)-
Γ
d
xz
Γ
BA
Γ
CI
+y-
γ
1
2
-
γ
ϕ
-
2
z
Γ
CI
+
1
2
(
Γ
BA
+
Γ
CI
)-
Γ
d
-z
Ω
x
(1-z)
γ
1
+z(1-
2
z
)
Γ
CI
+y
Ω
x
z
Γ
CI
As one can see, the four rows are term for term the drift in the display equation, with Wolfram ordering each sum its own way: the state rows carry the
-
Γ
CI
2
z
sharpening and the
Γ
BA
Γ
CI
cross term, and the readout row is the bare
Γ
CI
z
signal. Confirm the diffusion the same way:
In[]:=
diffStrat[{x,y,z,Q},ratesSym]//MatrixForm
Out[]//MatrixForm=
-y
Γ
BA
-xz
Γ
CI
x
Γ
BA
-yz
Γ
CI
(1-
2
z
)
Γ
CI
1
As expected, the three state rows are

b
and the readout row is the constant one, the additive noise that makes
Q
the same in both stochastic calculi. The code and the equations are one object, so from here on we compute with the code.
Fix a representative normalized benchmark, rescaled by
γ
1
; its readout is quantum-limited (
η=1
),
Γ
CI
+
Γ
BA
=2
Γ
d
exactly:
In[]:=
ratesB={0.252,0.1,0.176,1.,1.,3.};
Start from the ground state with the readout zeroed, a four-vector: the three Bloch numbers followed by the record
Q
:
In[]:=
initB={0.,0.,1.,0.};
Fix a short time window and the trajectory count:
In[]:=
tfB=3.;nB=600;
For the trajectory to stay a density matrix, the measurement must fit inside the total dephasing,
Γ
CI
+
Γ
BA
≤2
Γ
d
+
γ
ϕ

, an admissibility condition derived below. For
Γ
d
>0
it reads
η≤1+
γ
ϕ

Γ
d
; a physical detector with
η≤1
and nonnegative intrinsic dephasing satisfies it. Confirm this set passes:
In[]:=
ratesB[[1]]+ratesB[[2]]<=2(ratesB[[3]]+ratesB[[4]])
Out[]=
True
The quantum-limit claim is one division away; define the efficiency and read this set's:
In[]:=
ClearAll[efficiency];​​efficiency[r_]:=(r[[1]]+r[[2]])/(2r[[3]]);​​efficiency[ratesB]
Out[]=
1.
The efficiency is one at the input precision, by construction.
The reference for the mean is the Lindblad master equation. Because the Itô drift

a
is affine in
{x,y,z}
, the ensemble mean obeys that drift as a closed linear ODE. Encode

a
's three state rows once; the symbolic conversion check below must land on this same object, so one encoding serves both:
In[]:=
ClearAll[driftIto];​​driftIto[{x_,y_,z_},r_]:=​​With[{Γd=r[[3]],γϕ=r[[4]],γ1=r[[5]],Ωx=r[[6]]},​​With[{Γ2=γ1/2+γϕ+Γd},​​{-Γ2x,-Γ2y-Ωxz,γ1(1-z)+Ωxy}]];
Calling that encoding "the Lindblad master equation" requires a derivation. Build the generator itself, the drive
H=
Ω
x
2
σ
x
with relaxation
γ
1
𝒟[
σ
-
]
toward the
z=+1
ground state and total dephasing
γ
ϕ
+
Γ
d
2
𝒟[
σ
z
]
, where
𝒟[L]ρ=Lρ
†
L
-
1
2
(
†
L
Lρ+ρ
†
L
L)
; project it onto the Pauli components and subtract
driftIto
:
In[]:=
With[{sx=PauliMatrix[1],sy=PauliMatrix[2],sz=PauliMatrix[3],​​sm={{0,1},{0,0}},id=IdentityMatrix[2]},​​Module[{ρ,diss,gen},​​ρ=(id+xsx+ysy+zsz)/2;​​diss[L_]:=L.ρ.ConjugateTranspose[L]-​​(ConjugateTranspose[L].L.ρ+ρ.ConjugateTranspose[L].L)/2;​​gen=-I(ratesSym[[6]]/2)(sx.ρ-ρ.sx)+​​ratesSym[[5]]diss[sm]+((ratesSym[[4]]+ratesSym[[3]])/2)diss[sz];​​Simplify[(Tr[gen.#]&/@{sx,sy,sz})-driftIto[{x,y,z},ratesSym]]]]
Out[]=
{0,0,0}
Use forty target intervals across the shortest estimated timescale as the working mesh. This is a reproducible default, not a certified tolerance:

Regularizing the noise: a Wiener path NDSolve can integrate

On that grid, draw independent Gaussian increments with variance equal to the realized spacing, accumulate them, and interpolate the sampled path:
Visualize one realization of the path and its piecewise-defined derivative, the classical forcing that stands in for white noise:

The ensemble: matching the mean and the spread

For the harder regimes and the step-size study below we will want only the final values, so reduce each generator to them. The native final Bloch vectors:
The plot compares both sampled ensemble means with the exact Lindblad curve over the whole stored grid. In the white-noise theory the ensemble average obeys the master equation exactly; finite ensembles fluctuate around it.
This is the largest pointwise standardized gap on the stored grid. Because it is a maximum over many correlated times, ordinary one-, two-, or three-standard-error thresholds do not calibrate it as a simultaneous test. Use it as a scale diagnostic; a formal whole-curve test would require the sampling distribution of the maximum, for example from an independent bootstrap.
The mean sees only the drift, so it is a weak test. The measurement backaction lives in the standard deviation of the final Bloch vector, its spread across the ensemble, compared across the two routes:
Compare the Stratonovich and regularized spreads by dividing their absolute difference by their combined standard error,
Read each component's gap in those units:
The returned numbers express each spread difference in approximate combined-standard-error units. They are useful diagnostics, but they are not a distribution-free hypothesis test and do not justify the phrase “statistically indistinguishable” without a calibrated sampling analysis.
Verify that closed form, substituting a point on the sphere and subtracting the claim:
The radius runs past one. At the chosen initial boundary point the radial diffusion is zero while the radial generator points outward, so the model immediately violates the invariance condition. The equation still integrates; what fails is its interpretation as a density-matrix trajectory.
The radius sticks to the sphere at solver error, the sharpest test of the diffusion's scale in the document: a noise term wrong by any factor would walk the path off the sphere at once. On the boundary, positivity is not an inequality respected but an equality enforced.

The readout: the average signal and the shared noise

Overlay the ensemble-averaged record on that prediction:
As expected, the mean record follows the integrated mean signal. Each trajectory satisfies
The result measures the largest pointwise discrepancy in standard-error units. Because the maximum is taken over correlated times, pointwise Gaussian thresholds do not give its false-alarm probability. It is a diagnostic of scale, not a calibrated simultaneous test and not a claim about unsampled times.
The final-record spread supplies the corresponding diffusion-sensitive comparison for the fourth row. Here the two finite samples differ on the scale estimated above.
Rather than trust the line width, put a number on it, the largest gap between the stripped record and the driving path across the window; define the residual once, its solver goals adjustable, and read it at the default tolerance:
The residual decreases when the requested goals are tightened, as expected for integration error. The exact regularized ODE satisfies the identity; the reported nonzero residual is numerical.
Tighten the goals and the violation must follow the solver, not the noise grid:

The two time steps: which one sets the answer

Check the native integrator at three step sizes. The runs use the same integer seed, but changing the step changes how random numbers are consumed, so they should be treated as separate Monte Carlo samples rather than pathwise-coupled trajectories:
Read the two refined rows as distances from the working-step row, in units of the estimates' sampling errors in quadrature:
The refined estimates show no resolved trend relative to their Monte Carlo uncertainty. That supports using the working step as a practical reference at this sample size, but it does not prove convergence.
The table shows how the estimated mean and spread vary across the sampled noise meshes. Put the mean column in pointwise Monte Carlo standard-error units relative to the exact Lindblad value:

A strong readout: the same route where the drift is large

Read its efficiency:
Run the native reference for the strong case:
Run the regularized ensemble for the strong case:
Compare the mean of the final Bloch vector across the three routes:
As expected, the spreads agree closely where the Stratonovich correction is largest. Put the sampling scale on all three components again:
Repeat the finite mesh-sensitivity check in the strong-readout case, where the conversion correction is larger:
Now read each row as a distance from the native strong-case spread, in units of the two estimates' sampling errors in quadrature:
These three estimates show no resolved trend beyond their Monte Carlo scale. That supports the working mesh for the displayed observable and sample size; it does not certify the forty-interval heuristic for other rates, initial states, observables, or accuracy targets.

What is now true

The exact algebra is stronger than the numerical evidence: the Itô conversion, Lindblad generator, readout identity, and Bloch-ball flux all reduce symbolically to the claimed formulas. The simulations are finite-resolution checks. Their means and spreads agree with the stated references at the displayed scales, but regularization bias, stochastic-integrator bias, ODE error, and Monte Carlo uncertainty remain distinct error sources.

Acknowledgment

Thanks to Jose M. Martin-Garcia (Wolfram R&D) for suggesting this approach as a way to simulate random processes.

CITE THIS NOTEBOOK

Wong–Zakai regularization for continuous quantum measurement with NDSolve​
by Mohammad Bahrami​
Wolfram Community, STAFF PICKS, August 15, 2026
https://community.wolfram.com/groups/-/m/t/3781131