In this computational essay, we explore a quantum particle whose position is monitored continuously. This is what is usually called continuous weak measurement. It starts from the stochastic Schrödinger equation (the SDE that drives a state while it is being measured), turns the measurement step into ordinary Bayesian updating, and runs the whole thing on a position grid with a Trotter-Suzuki splitting scheme. We focus on two examples in particular. A free particle, where the measurement narrows a spread-out state to a packet of finite steady width and the noisy measurement record alone tells you where it settled. And a harmonic trap, where continuous measurement squeezes the position variance below the ground state's, and the quantum Kalman filter that levitated-nanoparticle and membrane experiments actually run rebuilds the conditional mean position from the noisy signal.
Part I: The Equation of Motion, and Measurement as Bayes' Rule
Part I: The Equation of Motion, and Measurement as Bayes' Rule
Start with the equation of motion. A particle with Hamiltonian whose position is measured continuously and weakly obeys a stochastic Schrödinger equation (SSE), the diffusive quantum-state-diffusion form,
H
In words: the first term is ordinary unitary evolution under ; the last is the stochastic measurement backaction, proportional to the Wiener increment (an independent Gaussian of mean zero and variance in each step), which nudges the state toward a sharper position; and the middle term is the nonlinear drift whose only job is to keep normalized. The strength sets how fast position information arrives, and its dimensions (inverse length squared, per time) say what kind of rate it is: what accrues is not position but statistical distinguishability: the log-likelihood ratio between two candidate positions apart grows at the rate (equivalently their Kullback-Leibler rate), a fact Part III will watch in action. (Factor conventions for and for the record differ across the literature by twos; translate parameters before comparing formulas across sources.) The detector's own output over the same increment is the measurement record
H
dW
dt
|ψ〉
λ
Δx
2λ
2
(Δx)
λ
which is the conditional mean position buried in white noise. Equivalently, the record is signal (the mean) plus shot noise (the same ). Because that one increment drives both the state and the record, an observer who knows the initial state can feed the record into the same dynamics and track the conditional state as the experiment runs, never seeing the wavefunction directly. That is quantum filtering, and it is what a real experiment computes, in software, from the detector's streaming output. Part IV does it concretely, and there the record fixes the conditional mean while the width follows a rule that never looks at the record.
x
dW
Keep one piece of bookkeeping straight from the start. In an experiment only the record is data; the conditional state, its mean, and its width are inferences from that data through an assumed model. In this essay even the record is synthetic, generated by the same equations that propagate the state, so whenever two of our computed objects agree below, the agreement is an internal consistency check of the model, never evidence for it. What evidence for the model itself looks like in a laboratory, whiteness tests on the filter's surprises and independent out-of-loop calibration, appears at the close of Part IV.
The above SSE is the limit of a concrete two-step update, and it is the update, not the differential form, that the simulation runs. Over one step the detector returns a single weak-measurement outcome , a reading of the position blurred by Gaussian noise, and the state is revised by Bayes' rule. The reading has likelihood
dt0
x
a Gaussian in centered on the candidate position . (The dropped constant is ; it cancels in the normalized update below, and needs restoring only if is read as a normalized density rather than a reweighting factor.) Because the reading is noisy, the outcome is not drawn from
directly but from the prior predictive , the density
convolved with the detector noise. And the conditional (posterior) density is the prior reweighted by the likelihood,
x
x
2λdt/π
P()
x
2
ψ(x) |
P()=∫P(x)
dx
x
x
2
ψ(x) |
2
ψ(x) |
which is Bayes' rule verbatim: posterior likelihood prior. At the level of amplitudes this reweighting is a Kraus operator ()∝, the square root of the likelihood, applied and renormalized; multiplying each amplitude by carries the phase along. The particle then evolves unitarily for the step, so one full step is Bayes' rule followed by Hamiltonian evolution:
∝
×
M
x
-λdt
2
(-)
x
x
e
P(x)
x
In short: measure by Bayes' rule, then evolve. (The ordering is a discretization, not physics: over a real interval the two proceed together, and splitting them costs an error first order in ; Part III puts a number on it.) Notice the one fact we will lean on for the rest of the notebook: the measurement half of this update knows nothing about . The outcome distribution, the Bayesian reweighting, and the renormalization are identical for a free particle or a trap; only the unitary carries the potential. So generalizing the simulation to any will mean generalizing that one factor and leaving the measurement untouched, which is exactly what the next Part does.
dt
H
-idtℏ
H
e
V(x)
Glossary
Glossary
The symbols of this Part sort into what a detector records and what an observer infers from it. What a detector actually returns is finite, one number per time step , so here we write those recorded increments with ; the continuous -forms of the equations above (, ) are their idealizations, living only inside the model. The last column gives each symbol's standing: recorded, calibrated from the record, or inferred through the model.
Δt
Δ
d
d
Y
t
d
W
t
Δt0
Symbol | Meaning | Standing, and how it is checked | |
Δ Y k | the measurement record over one step, Δ Y k λ x W k | Recorded: the detector's output each step, the only raw data here. | |
x k | one weak-measurement outcome per step, x k Y k λ Δt) | Calibrated from the record, not a separate reading: the same Δ Y k | |
x | the conditional mean position, ∫x 2
| Inferred: the filter of Part IV reconstructs it from the record, and feedback acts on it. | |
Δ W k | the shot-noise increment (variance Δt Δ Y k λ x | Inferred: the record minus the estimated mean. If this leftover is not pure noise, the model is wrong. |
Only the record is recorded; the outcome is that same increment calibrated into position units, and the conditional mean is inferred from the whole record by filtering, which is exactly what Part IV shows a real experiment doing in real time.
Δ
Y
k
x
k
x
Part II: Simulating Any Potential on a Grid
Part II: Simulating Any Potential on a Grid
Part I gave the update as an equation; to run it we have to make one decision that everything else follows from. A wavefunction is defined at every point of the infinite line, which is more than a computer can hold, so we keep it only on a finite grid: equally spaced points across a window of length , at spacing . The state becomes a list of complex numbers =ψ(), and every operation below is arithmetic on that list. The care is all in choosing the window and the spacing, and in carrying, next to the positions, the momenta we will need. Let me explain each piece before writing the grid.
ψ(x)
n
L
dx=L/n
n
ψ
j
x
j
Why a periodic box. Cutting space to a window of length leaves two loose ends, and we glue them: the point is identified with , so the segment becomes a ring of circumference . We choose this periodic boundary for one reason. Momentum states are the plane waves , and on a ring such a wave is allowed only if it closes up smoothly after going once around, which forces the momentum onto the discrete ladder
L
x=-L/2
x=+L/2
L
ipx/ℏ
e
In other words, periodicity is exactly what turns momentum from a continuum into an evenly spaced list we can store beside the positions, and it makes those plane waves an exact basis. It costs nothing as long as
stays negligible at the edges, so the packet never feels the seam where the ends are glued.
2
ψ |
How to choose the length and the spacing. The length must be large enough to hold the whole wavefunction, tails included, for the entire run: the packet spreads, drifts, and (under measurement) heats, and must contain all of it, or a tail that reaches one edge reappears at the other. The spacing must be small enough to resolve the fastest wiggle of , which is set by its largest momentum: a component of momentum oscillates with de Broglie wavelength (we keep the symbol for the measurement strength), and capturing a wave needs at least two samples across it, so small enough that the biggest momentum you expect stays below . The two ends are mirror images: the box length fixes the momentum spacing (fine momentum needs a big box), the sampling fixes the momentum reach (high momentum needs fine sampling). So choosing the position grid already fixes the momentum grid.
L
L
dx
ψ
p
2πℏ/p
λ
dx
πℏ/dx
dp=2πℏ/L
πℏ/dx
What the Fourier transform is, and why we want it. We will keep switching between two descriptions of the same state: how much amplitude sits at each position, the list , and how much sits at each momentum, the coefficients when we write as a sum of the plane waves above, . The recipe that turns the first list into the second is the Fourier transform: a reversible change of basis, the fixed sum = that correlates the state against each plane wave, with an inverse that rebuilds . (The built-in refers phases to the first grid point rather than to , which multiplies by a per-mode sign; it cancels in every round trip below and in
.)
ψ
j
c
m
ψ
ψ(x)=
∑
m
c
m
ix/ℏ
p
m
e
c
m
1
n
∑
j
ψ
j
-iℏ
p
m
x
j
e
ψ
Fourier
x=0
c
m
2
c m |
How a wavefunction sits on the grid. After all that, representing is the easy part: evaluate it at the grid points and normalize, with
=1. The normalization makes each entry carry a factor : the list stores ), so
is probability per grid point rather than a density, and every expectation below is a plain dot product for exactly that reason. That list is the state in the position description; its Fourier transform is the same state in the momentum description. First fix units so that :
ψ(x)
ψ(ψ(),…,ψ())
x
1
x
n
∑
j
2
ψ j |
dx
dx
ψ(x
j
2
ψ j |
ℏ=m=1
In[]:=
hbar=1;mass=1;
Define the grid given the space and the size of the box:
n
l
In[]:=
grid[n_,l_]:=With[{dx=l/n},<|"n"->n,"dx"->dx,"x"->dxRange[-n/2,n/2-1],"p"->(2Pihbar/l)RotateLeft[Range[-n/2,n/2-1],n/2]|>];
The grid itself returns an association with , the spacing , the positions =dx{-n/2,…,n/2-1} (centered on the origin, so the box is ), and the quantized momenta =2πℏm/L. The momenta are stored rotated, nonnegative first and the negative half moved to the end, because that is the order returns them in, so that later one multiplication lands each amplitude on its own momentum. Fix the transform's sign convention once, as a pair used everywhere below:
n
dx
x
j
[-L/2,L/2)
p
m
m
Fourier
In[]:=
ft[ψ_]:=Fourier[ψ,FourierParameters->{0,-1}];ift[c_]:=InverseFourier[c,FourierParameters->{0,-1}];
A convenient initial state is a Gaussian wave packet centered at with mean momentum and position width . Following the recipe just above, we write its formula on the grid points and normalize:
x
0
p
0
s
g["x"]
In[]:=
gaussian[g_,x0_,p0_,s_]:=Normalize[Exp[-(g["x"]-x0)^2/(4s^2)+Ip0g["x"]/hbar]];
The two numbers we will watch most are the conditional mean position and the conditional variance =-. On the grid each is a dot product against the probability density
:
=∫x
dx
x
2
ψ(x) |
Σ
x
2
x
2
x
2
ψ |
In[]:=
xMean[g_,ψ_]:=g["x"].Abs[ψ]^2;xVar[g_,ψ_]:=g["x"]^2.Abs[ψ]^2-xMean[g,ψ]^2;
For the measurement step, we implement it in the two pieces: first the outcome , then the Bayesian update of the state.
x
The outcome is a sample from the prior predictive . Recall that is
convolved with Gaussian noise of variance . So for the sampling of , we can treat it as =X+ξ, with distributed as
and distributed as a Gaussian with zero mean and standard deviation :
x
P()
x
P()
x
2
ψ |
1/(4λdt)
x
x
X
2
ψ(x) |
ξ
1/
4λdt
In[]:=
drawSyntheticOutcome[g_,λ_,dt_][ψ_]:=RandomChoice[Abs[ψ]^2->g["x"]]+RandomVariate[NormalDistribution[0,1/Sqrt[4λdt]]];
The intermediate draw is a sampling device for that convolution, nothing more: it is discarded on the spot, and no position is thereby assigned to the particle; only leaves the detector. ( is exact in that sense: it samples the finite-step outcome law without approximation.)
X
x
drawSyntheticOutcome
Second the Bayesian update itself: reweight the amplitudes by the square root of the likelihood, , and renormalize. That is the whole measurement step:
-λdt
2
(x-)
x
e
In[]:=
measUpdate[g_,λ_,dt_,xb_][ψ_]:=Normalize[Exp[-λdt(g["x"]-xb)^2]ψ];
(One numerical caveat for reuse: a wildly off-center , or a very strong measurement, can drive every weight below machine range before , which would then return a zero vector; at the strengths here that never bites, but a hardened version subtracts the peak exponent before exponentiating.)
x
Normalize
Is this really Bayes' rule? It should be: the updated density
ought to equal the prior
times the likelihood , normalized. Confirm the two constructions agree, for a test state and outcome:
2
ψ post |
2
ψ(x) |
-2λdt
2
(x-)
x
e
In[]:=
With[{gDemo=grid[512,40.],λ=.75,dt=.01,xb=-.5},With[{ψt=gaussian[gDemo,1.,0.5,2.]},With[{post=Abs[measUpdate[gDemo,λ,dt,xb][ψt]]^2,bayes=Normalize[Exp[-2λdt(gDemo["x"]-xb)^2]Abs[ψt]^2,Total]},Max@Abs[post-bayes]]]]
Out[]=
1.73472×
-17
10
Before the runs themselves, two one-line helpers they will need.
Two more checks catch the two ways a finite grid can quietly lie. The box is a ring, so probability that reaches the wall wraps around to the far side; and the grid holds only a limited span of momenta, so probability that reaches the fastest ones can no longer be told from slow motion. Both stay harmless only while the packet keeps clear of those edges, so we measure how much probability has reached each: the outer tenth of the box, and the fastest tenth of the momenta. Both are the same sum over the outer tenth of an axis, so define that once as a helper and apply it to the position and to the momentum:
Part III: The Free Particle: No Potential, Localizing and Heating
Part III: The Free Particle: No Potential, Localizing and Heating
Run one trajectory from a wide packet of width 2, with the noise seeded so the run repeats exactly:
Take the conditional width at each step, the square root of the position variance:
Then the momentum leakage: the largest share that ever reaches the fastest tenth of the momenta, which stays negligible only while the packet never reaches the top of the momentum ladder.
The offset is first order in the step, so halving the step should roughly halve it: the offset at the full step divided by the offset at the half step should come out near two.
The surviving offset halves with the step: pure splitting bias, first order exactly as Part II promised, and none of it physics.
At unit detection efficiency the conditional state stays pure, which is what lets us carry it as a single state vector; the state only becomes mixed when we average over the noise, which is the master equation, coming next.
Now the second known fact, the measurement backaction. Work on a small 32-point grid over a box of width 16:
Use a fine step and a short total time:
Start from a single narrow packet at the origin:
Build the cat as the normalized sum of a left packet and a right packet:
Look at its two lobes:
Then raise the density over position and time:
The two ridges compete and one grows while the other fades, with no potential in sight.
Each outcome contributes a log-odds increment, left lobe against right; define it:
Accumulate those increments and turn the running total into the posterior probability for the left lobe:
First look at the raw stream the detector logs, the signal buried in noise:
So the outcomes say left lobe and the grid agrees. The road there is the instructive part: the belief builds only over many steps, and can even lean the wrong way first. For the comparison, first read the exact left-lobe mass the grid state carries at each step, which the observer never sees:
The runs split between the two lobes, each driven to a definite verdict by its own outcomes. This is the outcome stream doing the observer's inference, the same accumulation a real experiment performs in software beside the detector. Now we turn on a potential and ask what survives.
Part IV: The Harmonic Trap: Position Sharper Than the Ground State
Part IV: The Harmonic Trap: Position Sharper Than the Ground State
Feed the trap Hamiltonian into the stochastic Schrödinger equation of Part I and apply it to each moment with Itô calculus. The equations of motion split cleanly in two. The two means follow the classical Hamilton drift plus a measurement gain acting on the record's noise:
Solve for the fixed point where all three rates vanish:
Exactly one root passes it:
Select that physical branch:
Both check out: the steady state is a pure Gaussian that reduces to the free particle when the trap is switched off.
Define a function, and the strength beside it:
And confirm the two forms agree:
Second, that its slope is negative everywhere, so watching harder always squeezes the position further:
Its weak-measurement limit:
And its strong-measurement limit:
Confirm it:
Evolve one trajectory from the ground state and take its conditional width over the run:
Plot it in units of the vacuum, against the closed-form Riccati value (red):
The grid variance relaxes to the Riccati steady value and sits well below the zero-point; the small remaining offset is the same first-order splitting bias Part III measured, independent of the grid spacing, and it shrinks with the step. The relaxation is the same in every run, because the conditional variance obeys a deterministic (noiseless) equation while only the mean is driven by the record.
Reading the Position Out of the Noise
Reading the Position Out of the Noise
The point of watching the particle is to know where it is. But the detector hands back something discouraging: a stream of position readings, each one almost pure noise. Can you recover the particle's mean position from that? You can, and the surprise is that you don't need the full wavefunction to do it, only five numbers, updated as each reading arrives. Those five numbers are themselves the conditional quantum state, its two means and three covariances; what you skip is the 256-amplitude grid, not the state. This five-number recipe is the quantum Kalman filter, the quantum form of the Kalman filter that turns a noisy signal into a running state estimate throughout engineering. It is what a trapped-nanoparticle or vibrating-membrane experiment does, inside a richer noise model: it never holds the full wavefunction, it watches the measured signal and runs this small calculation to know the mean position moment by moment.
Also keep the conditional mean the full grid computes, the reference the filter is checked against later:
Look at the readings on their own. Each scatters far wider than the conditional mean it carries, so the readings are mostly noise and the signal cannot be picked out by eye:
The gain is nonzero only on the two means, so the readings move them and never the spreads. The last three rows of the drift are the covariance block, exactly the symbolic Riccati of this Part with the numbers substituted, so one definition serves both the algebra and the filter:
The surprise itself, the innovation, is the newest reading minus what the estimate predicted, and it reduces to the outcome and the step alone:
One step is that update:
Run the filter along the whole record from the ground state, and keep its estimate of the conditional mean:
Compare that estimate with the grid's conditional mean, taking the largest gap between them as a fraction of the full range that mean covers over the run:
The gap is a tiny fraction of that range: the small price of finite time steps, not a real disagreement between the five numbers and the grid. From the noisy signal alone, five numbers recover the conditional mean that the full 256-point grid computed. Overlay them:
The plot overlays two curves for the average position over time, and as with the lobe plot only one of them is something an experiment could actually draw. The red dashed curve is the estimate the five-number filter reports from nothing but the noisy stream of readings, so it is observable, the very quantity a real detector tracks. The gray curve is the average position carried by the full 256-amplitude wavefunction, and since no experiment holds that wavefunction it is not observable: it is drawn only as the reference to check the filter against, and the red lies right on top of it. The filter reaches it without ever building the wavefunction: it starts from the known initial state and, at each reading, nudges its running estimate by the surprise, how far the reading fell from what the estimate predicted, adjusting only its five numbers. So five numbers, driven by the raw signal alone, reproduce the average position that 256 complex amplitudes compute. That is what a real experiment does: to keep track of the particle's average position it never holds the wavefunction, it runs these five numbers in real time from the measured light, as a levitated nanoparticle or a vibrating membrane experiment does. Part III's outcomes answered a single yes-or-no question, which packet won; this one hands back the whole moving conditional mean.
The Price of an Imperfect Detector
The Price of an Imperfect Detector
The first cost is purity. A perfect measurement held the state exactly at the uncertainty limit; a lossy one holds it above that limit, by exactly the missing fraction:
and the curve between the two ends prices every detector:
Part V: Averaging Over Runs: A Trapped State That Never Settles
Part V: Averaging Over Runs: A Trapped State That Never Settles
Now build the run from the ground state and tabulate. Start the density matrix from the trap ground state:
Choose a handful of times to sample:
The tabulated moments sit fifty-fold or more below both bounds across this window; run long enough and any finite grid must eventually betray the ramp when the packet meets a cutoff. Watch that betrayal happen, on Part III's deliberately small free grid run far past its comfort. Run the free master equation for many thousands of steps:
By the end enough has piled into that momentum tail for the grid to break down, exactly where the two curves parted:
From there on the small grid is lying about the continuum, the honest content of every cutoff warning above. The trap's ramp, plotted before this caveat, sits inside its safe window, so it carries no such lie.
A genuine thermal steady state needs a mechanical damping channel to balance that input, which is exactly why the optomechanics experiments cool with feedback rather than watch the bare measurement; adding that channel to the master equation is the natural next extension.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Continuous quantum measurement, from Bayes' rule to the quantum Kalman filter
by Mohammad Bahrami
Wolfram Community, STAFF PICKS, August 5, 2026
https://community.wolfram.com/groups/-/m/t/3774755
by Mohammad Bahrami
Wolfram Community, STAFF PICKS, August 5, 2026
https://community.wolfram.com/groups/-/m/t/3774755