Quantum Jump Simulation

A photon-counting trajectory describes a state conditioned on detected emissions. For an undriven excited atom, the entire trajectory is determined by one random decay time. Sampling that time exactly gives a useful first simulation: there is no integration-step bias, so discrepancies from the master equation come from finite sampling.

Helpful background. Quantum-jump trajectories explain conditional evolution; amplitude damping supplies the ensemble answer. The code uses Python 3 and NumPy.

One emission from an excited atom

Work in the interaction picture with $H=0$ and $L=\sqrt\Gamma|g\rangle\langle e|$. Assume ideal photon counting, unit efficiency, and no other channels. The unnormalized no-count state is

$$|\widetilde\psi(t)\rangle=e^{-\Gamma t/2}|e\rangle.$$

Its squared norm gives the survival probability $S(t)=e^{-\Gamma t}$. Differentiation gives the first-emission density $w(t)=\Gamma e^{-\Gamma t}$. After emission the state is $|g\rangle$ and stays there, so there is at most one jump.

If $u$ is uniform on $(0,1)$, then $\tau=-\log u/\Gamma$ has precisely this waiting-time distribution. NumPy’s exponential sampler performs the equivalent sampling directly. A time grid is used only to display the trajectories, not to decide when an emission occurs.

Run the simulation and quantify sampling error

Set $\Gamma=1$ in the chosen time unit. For $N$ independent trajectories, the fraction surviving at time $t$ estimates $p_e(t)$. Its standard deviation as an estimator is

$$\operatorname{SE}[\widehat p_e(t)] =\sqrt{\frac{S(t)[1-S(t)]}{N}}.$$

This uncertainty is statistical; reducing a plot’s time spacing does not reduce it.

import numpy as np

rng = np.random.default_rng(20261013)
Gamma = 1.0
N = 100000
tau = rng.exponential(scale=1/Gamma, size=N)
for t in [0.25, 0.5, 1.0, 2.0]:
    estimate = np.mean(tau > t)
    exact = np.exp(-Gamma*t)
    se = np.sqrt(exact*(1-exact)/N)
    assert abs(estimate-exact) < 6*se
    print(t, 'survival:', estimate, 'exact:', exact, 'SE:', se)
assert abs(tau.mean()-1/Gamma) < 6/(Gamma*np.sqrt(N))
print('mean waiting time:', tau.mean())

# A finite observation window censors late emissions.
T = 1.0
observed = tau[tau <= T]
conditional_mean = 1/Gamma-T/np.expm1(Gamma*T)
assert abs(observed.mean()-conditional_mean) < 0.01
print('mean among clicks before T:', observed.mean())

# General initial state (|g> + i|e>)/sqrt(2).
# Half the trajectories never emit; finite times have density
# (Gamma/2)*exp(-Gamma*t).
u = rng.random(N)
first = np.full(N, np.inf)
finite = u < 0.5
first[finite] = -np.log1p(-2*u[finite])/Gamma
no_count = first > T
v = np.array([1, 1j*np.exp(-Gamma*T/2)]) / np.sqrt(2)
v /= np.linalg.norm(v)
rho_no = np.outer(v, v.conj())
rho_g = np.diag([1.0, 0.0])
mean_rho = no_count.mean()*rho_no+(1-no_count.mean())*rho_g
exact_rho = np.array([[1-0.5*np.exp(-Gamma*T),
                       -0.5j*np.exp(-Gamma*T/2)],
                      [0.5j*np.exp(-Gamma*T/2),
                       0.5*np.exp(-Gamma*T)]])
assert np.linalg.norm(mean_rho-exact_rho) < 0.01
assert np.linalg.eigvalsh(mean_rho).min() >= -1e-12
print('ensemble state:', mean_rho, sep='\n')

At $t=1$, $S=e^{-1}\approx0.367879$ and the standard error for $N=100000$ is about $0.001525$. A different seed changes the estimate within statistical fluctuations. Error need not decrease monotonically for every new seed or sample size; the $N^{-1/2}$ law describes its scale across repeated samples.

Two traps in interpreting the record

A finite window changes the observed mean. If you discard trials with no click before $T$, the remaining times have conditional density $w(t)/(1-e^{-\Gamma T})$ on $[0,T]$. Their mean is

$$\mathbb E[\tau\mid\tau\le T] =\frac1\Gamma-\frac{T}{e^{\Gamma T}-1}.$$

At $\Gamma T=1$ this is about $0.418023/\Gamma$, not $1/\Gamma$. Keep the no-click trials when estimating a survival curve; do not normalize a truncated histogram and compare it with an untruncated exponential density.

A no-count state can change continuously. For $|\psi_0\rangle=a|g\rangle+b|e\rangle$, the unnormalized no-count state is $a|g\rangle+b e^{-\Gamma t/2}|e\rangle$. Its survival probability is $|a|^2+|b|^2e^{-\Gamma t}$, leaving probability $|a|^2$ of never emitting. Conditioning on no count gradually favors the ground component. The last part of the code includes this possibility and averages density matrices, not state vectors.

With inefficient detection, “no observed click” can also include an undetected emission. Its conditioned state is generally mixed. Deleting a random subset of the displayed clicks while retaining ideal pure conditional states would give the wrong state estimate.

Extend to a driven atom

For $H=\hbar\Omega X/2$, the atom can emit repeatedly. Evolve the unnormalized state with $H_{\mathrm{eff}}=H-i\hbar L^\dagger L/2$ until its squared norm reaches a sampled uniform threshold, apply the jump, then restart. A numerical implementation must resolve threshold-crossing times and compare its ensemble with a Lindblad solver.

A fixed-step method instead uses $p_{\rm jump}\simeq\Delta t\langle L^\dagger L\rangle$. It introduces time-step bias; sufficiently many trajectories cannot remove that bias. In the simple excited-state example its survival probability after $n$ steps is $(1-\Gamma\Delta t)^n$, rather than $e^{-\Gamma n\Delta t}$.

Exercises

Waiting-time mean. Integrate $tw(t)$ and find the variance.

Solution

The first two moments are $1/\Gamma$ and $2/\Gamma^2$, so the variance is $1/\Gamma^2$. The sample mean over $N$ independent complete waiting times has standard error $1/(\Gamma\sqrt N)$.

A no-count update. Start in $(|g\rangle+i|e\rangle)/\sqrt2$. At $t=(\log2)/\Gamma$, find the no-count probability and the conditional excited population.

Solution

The no-count probability is $(1+1/2)/2=3/4$. Its unnormalized excited population is $1/4$, so its normalized excited population is $1/3$. The unconditional excited population remains $1/4$ after averaging the no-count and jump branches.

A sampling budget. How many excited-state trajectories give standard error at most $0.001$ at $\Gamma t=1$?

Solution

Require $N\ge e^{-1}(1-e^{-1})/10^{-6}\approx232544.16$, hence at least $232545$. This is a standard-error target, not a guarantee that every observed error is smaller than $0.001$.

Step bias. Expand the fixed-step survival probability at fixed $t=n\Delta t$.

Solution

$\log S_{\Delta t}=n\log(1-\Gamma\Delta t)=-\Gamma t-\Gamma^2t\Delta t/2+O(\Delta t^2)$. Thus this method underestimates survival, with leading relative error $-\Gamma^2t\Delta t/2$.

References

• H. J. Carmichael, An Open Systems Approach to Quantum Optics, Springer (1993).
• J. Dalibard, Y. Castin, and K. Mølmer, “Wave-function approach to dissipative processes in quantum optics,” Physical Review Letters 68, 580–583 (1992). doi:10.1103/PhysRevLett.68.580.
• K. Mølmer, Y. Castin, and J. Dalibard, “Monte Carlo wave-function method in quantum optics,” Journal of the Optical Society of America B 10, 524–538 (1993). doi:10.1364/JOSAB.10.000524.
• M. B. Plenio and P. L. Knight, “The quantum-jump approach to dissipative dynamics in quantum optics,” Reviews of Modern Physics 70, 101–144 (1998). doi:10.1103/RevModPhys.70.101.