Diffusive Trajectory Simulation

A continuous measurement produces a noisy record and a state conditioned on that record. This example measures a qubit’s $Z$ component without a Hamiltonian. Because the measured observable is conserved, the conditional state can be calculated exactly from the accumulated record. This gives a positive-state benchmark for more general stochastic solvers.

Helpful background. Diffusive trajectories explain the state and record equations; measurement backaction distinguishes conditional information from ensemble disturbance. The code uses Python 3 and NumPy.

Continuous Z measurement and its record

Take $Z=\operatorname{diag}(1,-1)$, $\kappa>0$, and detection efficiency $0\le\eta\le1$. The Itō equations are

$$\begin{aligned} d\rho_c&=\kappa(Z\rho_cZ-\rho_c)dt +\sqrt{\eta\kappa}(Z\rho_c+\rho_cZ-2z_c\rho_c)dW,\\ dY&=2\sqrt{\eta\kappa}\,z_cdt+dW, \qquad z_c=\operatorname{Tr}(Z\rho_c). \end{aligned}$$

Here $dW$ is a Wiener increment with mean zero and variance $dt$. Ignoring the record gives $\rho_{01}(t)=e^{-2\kappa t}\rho_{01}(0)$ for every efficiency. Lower efficiency loses information; it does not reduce the total disturbance specified by $\kappa$.

For either initial state $I/2$ or $|+\rangle\langle+|$, the initial $Z$ populations are equal. The distribution of the record can be sampled as follows: choose $s=+1$ or $-1$ with equal probabilities and generate

$$Y(t)=2\sqrt{\eta\kappa}\,s t+W(t).$$

This is a convenient sampling identity for this commuting measurement’s record statistics. It does not assert that an initially coherent qubit secretly had a definite $Z$ value. The state conditioned on the record must still be computed with the quantum measurement update.

An exact conditional update

The likelihood ratio for $s=+1$ versus $s=-1$ is $e^{4\sqrt{\eta\kappa}Y(t)}$. Therefore

$$p_+(t)=\frac{1+\tanh[2\sqrt{\eta\kappa}Y(t)]}{2}, \qquad z_c=2p_+-1.$$

For an initial $|+\rangle$, the off-diagonal entry is

$$\rho_{01,c}(t)=\sqrt{p_+(t)[1-p_+(t)]}\,e^{-2(1-\eta)\kappa t}.$$

The square-root factor describes conditioning on the observed measurement. The exponential describes the unobserved dephasing. For initial $I/2$, the off-diagonal entry instead remains zero; the same record distribution and diagonal update apply.

The determinant is $p_+(1-p_+)[1-e^{-4(1-\eta)\kappa t}]\ge0$, so the coherent initial state’s update preserves positivity at every sampled time. At $\eta=1$ it remains pure. This exact formula depends on $H=0$ and a single conserved measured observable; it is not a general-purpose stochastic integrator.

import numpy as np

def simulate(eta, N=20000, steps=200, T=1.0, kappa=1.0):
    rng = np.random.default_rng(20261013)
    s = rng.choice([-1.0, 1.0], size=N)
    dt = T/steps
    Y = np.zeros(N)
    sample_record = []
    for _ in range(steps):
        dY = 2*np.sqrt(eta*kappa)*s*dt + np.sqrt(dt)*rng.normal(size=N)
        Y += dY
        sample_record.append(dY[0])
    z = np.tanh(2*np.sqrt(eta*kappa)*Y)
    pp = (1+z)/2
    coherence = np.sqrt(pp*(1-pp))*np.exp(-2*(1-eta)*kappa*T)
    purity_plus = pp**2+(1-pp)**2+2*coherence**2
    purity_mixed = pp**2+(1-pp)**2
    determinants = pp*(1-pp)-coherence**2
    assert determinants.min() >= -1e-12
    assert np.max(purity_plus) <= 1+1e-12
    assert abs(z.mean()) < 6*z.std(ddof=1)/np.sqrt(N)+1e-12
    target = 0.5*np.exp(-2*kappa*T)
    se = coherence.std(ddof=1)/np.sqrt(N)
    assert abs(coherence.mean()-target) < 6*se+1e-12
    print('eta:', eta, 'mean coherence:', coherence.mean(),
          'target:', target, 'SE:', se,
          'mean purity from I/2:', purity_mixed.mean())
    return z, coherence, purity_plus, purity_mixed, np.array(sample_record)

for eta in [0.0, 0.3, 1.0]:
    z, coherence, purity, mixed_purity, record = simulate(eta)
    if eta == 0:
        assert np.allclose(z, 0)
        assert np.allclose(mixed_purity, 0.5)
    if eta == 1:
        assert np.allclose(purity, 1)

# Why naive Euler-Maruyama can fail even at a small step:
# start from |+>, so x=1, z=0, and choose dW=2*sqrt(dt).
dt, kappa, eta = 0.01, 1.0, 1.0
dW = 2*np.sqrt(dt)
x = 1-2*kappa*dt
z = 2*np.sqrt(eta*kappa)*dW
minimum_eigenvalue = (1-np.hypot(x, z))/2
assert minimum_eigenvalue < 0
print('Euler step minimum eigenvalue:', minimum_eigenvalue)

The mean coherence should agree with $e^{-2}/2\approx0.067668$ within its sampling uncertainty at all three efficiencies. The ensemble of initially mixed states stays $I/2$, although the conditional states become more definite when $\eta>0$. In this ideal quantum nondemolition model, any fixed positive efficiency eventually distinguishes $Z=\pm1$; inefficiency slows purification rather than imposing a permanently mixed limit. Additional dynamics can change that conclusion.

Comparing exact records with a stochastic integrator

For Euler–Maruyama, use $\Delta W=\sqrt{\Delta t}\,\xi$ with independent unit normal $\xi$, and evaluate the coefficients at the beginning of each step. However, a trace-preserving Euler update can still leave the Bloch ball. The final code lines deliberately demonstrate this failure for a two-standard-deviation noise increment. Symmetrization or normalization does not repair a negative eigenvalue.

A general solver should be compared with this exact commuting example using both smaller steps and larger ensembles. They address different errors. For the exact update above, changing the sampling grid alone introduces no state-propagation error; it merely changes which points of the continuous record are displayed.

The innovation is $dW=dY-2\sqrt{\eta\kappa}z_cdt$. In the continuous-time limit it is Wiener noise relative to the observed past. If you reconstruct it from finite record bins using only the state at the bin’s beginning, finite-bin corrections remain. Do not demand an exact Gaussian identity at arbitrarily large bin size. Also, the raw current $\Delta Y/\Delta t$ has noise variance $1/\Delta t$; compare currents with the same averaging bandwidth.

Exercises

Likelihood ratio. Derive the exponent in the update for $p_+$.

Solution

For $a=2\sqrt{\eta\kappa}$, the two record densities are normal with means $\pm at$ and variance $t$. Their logarithmic ratio is $[-(Y-at)^2+(Y+at)^2]/(2t)=2aY=4\sqrt{\eta\kappa}Y$. Equal prior probabilities give the stated logistic, or hyperbolic-tangent, update.

Purity and missing information. Show that a pure initial state remains pure for $\eta=1$, and explain the case $\eta=0$.

Solution

At unit efficiency, $|\rho_{01,c}|^2=p_+(1-p_+)$, so the determinant is zero and the trace-one state is pure. At zero efficiency the record is independent noise, $p_+=1/2$, and the coherence is $e^{-2\kappa t}/2$. The state follows unconditional dephasing.

Positivity of an Euler step. Find the smallest eigenvalue in the deliberate failure above.

Solution

The updated Bloch vector is $(0.98,0,0.4)$, so its length is $\sqrt{1.1204}>1$. The smaller eigenvalue is $(1-\sqrt{1.1204})/2\approx-0.02924$ despite exact trace preservation.

Does every trajectory purify monotonically? Interpret a short interval in which a mixed conditional state becomes less pure.

Solution

A noisy record can temporarily move the posterior probabilities back toward equality. Individual purity need not increase monotonically; information gain is a statistical statement. In this conserved-$Z$ model, the long-time record nevertheless selects an eigenstate almost surely when $\eta>0$.

References

• H. J. Carmichael, An Open Systems Approach to Quantum Optics, Springer (1993).
• K. Jacobs, Quantum Measurement Theory and its Applications, Cambridge University Press (2014).
• K. Jacobs and D. A. Steck, “A straightforward introduction to continuous quantum measurement,” Contemporary Physics 47, 279–303 (2006). doi:10.1080/00107510601101934; author manuscript.
• H. M. Wiseman and G. J. Milburn, Quantum Measurement and Control, Cambridge University Press (2010).