Computational Notebooks

Numerical work becomes trustworthy when the physics, representation, algorithm, and validation evidence can all be inspected independently. A plot that looks plausible is not sufficient. A mature computational result should survive analytic limiting cases, physicality checks, resolution changes, alternative implementations, and a clean rerun from recorded inputs.

The pages in this chapter are notebook specifications. As of this review, they do not claim that corresponding executable artifacts have been reproduced or promoted. Each page defines an admission contract: the model, conventions, calculations, tests, metadata, and accepted outputs required before numerical results can be cited as reproduced artifacts.

The Notebook Contract

Every notebook in this chapter should contain nine layers.

1. Physics goal: State the question and the observable to be computed.
2. Mathematical model: Give the Hamiltonian, channel, generator, stochastic equation, or work protocol.
3. Units and parameters: Record dimensions, scales, basis order, frames, and numerical values.
4. Analytic limiting case: Identify at least one result known independently of the numerical method.
5. Numerical method: Name the integrator, matrix representation, stochastic convention, optimizer, and approximation.
6. Convergence checks: Vary time step, tolerance, truncation, trajectory count, or sampling resolution as appropriate.
7. Reproducible implementation: Pin dependencies, record random seeds, and separate source inputs from generated outputs.
8. Validation test: Compare against an exact solution, structural identity, second method, or independently computed benchmark.
9. Extensions: Mark optional investigations separately from the validated baseline.

The contract makes a notebook more than an illustrated derivation. It states what would falsify the calculation.

Reading Path

Read this page	Use it for
Simulating Quantum Channels	Constructing Kraus and Choi representations, applying channels, and testing complete positivity and trace preservation.
Bloch Vector Noise Models	Visualizing dephasing, depolarization, and amplitude damping as affine maps of a qubit Bloch vector.
Solving Lindblad Equations	Comparing ODE and matrix-exponential solvers, finding steady states, and inspecting Liouvillian spectra.
Quantum Jump Simulation	Generating jump records, no-jump evolution, waiting-time histograms, and ensemble averages.
Diffusive Trajectory Simulation	Integrating stochastic master equations, generating noisy records, and checking conditional purification.
Non-Markovian Toy Models	Comparing exact enlarged-system evolution with Markovian reductions and trace-distance revivals.
Decoherence Timescale Estimation	Computing coherence envelopes from spectra and filter functions and defining fitted timescales.
Optimal Control Toy Problems	Optimizing small pulse families under dephasing and testing fidelity and robustness.
Quantum Thermodynamics Toy Models	Building work distributions and testing Jarzynski and Crooks relations for a driven two-level system.

A productive first route is channels, Bloch vectors, and Lindblad equations. The two trajectory notebooks then add stochastic records. The remaining notebooks test memory, spectral estimation, control, and thermodynamic bookkeeping.

From Physics Statement to Numerical Artifact

A reliable workflow proceeds in a fixed order.

Specify the observable

Begin with the quantity that will answer the physics question. Examples include

$$\langle \sigma_z(t)\rangle, \qquad \operatorname{Tr}[\rho(t)^2], \qquad g^{(2)}(\tau),$$

or a work probability $P(W)$. This prevents the implementation from expanding into an unstructured simulation of every available variable.

Fix representation conventions

State basis ordering, tensor-factor ordering, vectorization, complex-conjugation conventions, and units. If column stacking is used, then

$$\operatorname{vec}(A\rho B) = (B^{\mathsf T}\otimes A) \operatorname{vec}(\rho).$$

Using row stacking changes the matrix representation of the superoperator. A solver may still run after a convention mismatch while producing a physically incorrect evolution.

Nondimensionalize when useful

Choose a reference frequency $\Omega_0$ or time $t_0$ and expose the dimensionless combinations that control the dynamics:

$$\tau=\Omega_0t, \qquad \tilde\gamma=\frac{\gamma}{\Omega_0}, \qquad \tilde\Delta=\frac{\Delta}{\Omega_0}.$$

Nondimensionalization improves conditioning and makes parameter regimes easier to compare. Physical units should still be recoverable from recorded metadata.

Derive a baseline before implementing

Every notebook needs a case whose answer is known. For zero-temperature qubit relaxation,

$$\rho_{ee}(t) = \rho_{ee}(0)e^{-\gamma t}.$$

For pure dephasing under the convention used by the notebook, the off-diagonal element must follow the corresponding analytic exponential. The convention must be checked rather than inferred from the name of a collapse operator.

Implement the smallest validated case

Start with the minimum Hilbert space, one parameter set, and one observable. Add sweeps, optimization, or additional modes only after that baseline passes its tests. Complexity should enter in auditable increments.

Structural Physicality Checks

Numerical state evolution should be tested against properties that do not depend on the plotted observable.

Density operators

For every reported state, monitor

$$\epsilon_{\mathrm{tr}} = \left| \operatorname{Tr}\rho-1 \right|,$$ $$\epsilon_{\mathrm{H}} = \left\| \rho-\rho^\dagger \right\|,$$

and the smallest eigenvalue $\lambda_{\min}(\rho)$. Small negative eigenvalues can arise from roundoff, but persistent or tolerance-dependent negativity can signal a bad integrator, an invalid generator, or an inconsistent approximation.

Trace, Hermiticity, and positivity tests serve different purposes. Renormalizing the trace does not repair non-Hermiticity or negative eigenvalues.

Quantum channels

For Kraus operators $\{K_\alpha\}$, a trace-preservation residual is

$$\epsilon_{\mathrm{TP}} = \left\| \sum_\alpha K_\alpha^\dagger K_\alpha -I \right\|.$$

For a Choi matrix $J(\Phi)$, check Hermiticity, positive semidefiniteness, and the appropriate partial-trace identity. A few positive test states cannot establish complete positivity.

Master equations

A finite-dimensional Lindblad solver should preserve trace and Hermiticity and should keep positive initial states positive to numerical tolerance. It should also reproduce known stationary states and decay rates. If the generator is represented as a matrix $\mathbb L$, verify

$$\frac{d}{dt} |\rho\rangle\!\rangle = \mathbb L |\rho\rangle\!\rangle$$

against direct evaluation of the operator equation on randomly selected test states.

Stochastic trajectories

Conditional states require normalization, record statistics, and ensemble checks. A trajectory algorithm can preserve each state’s norm while using the wrong jump probabilities or stochastic drift. Structural tests must therefore include both state and record.

Convergence Is Part of the Result

A numerical value without a resolution study is an estimate with an unknown numerical error.

Time-step and solver tolerance

Let $O_h(t)$ denote an observable computed at step size $h$. A simple refinement residual is

$$R_h = \max_t \left| O_h(t)-O_{h/2}(t) \right|.$$

Report $R_h$ over the interval actually used in the analysis. One matching endpoint does not establish convergence of the trajectory.

Adaptive solvers require absolute and relative tolerances, an error norm, and any maximum-step restriction. Tightening only one tolerance may not probe the dominant error.

Hilbert-space truncation

Bosonic modes require a cutoff $n_{\max}$. Convergence should be checked by increasing the cutoff and monitoring both observables and boundary population. A small occupation of the highest retained level is useful evidence but not sufficient when dynamics repeatedly reaches the boundary.

Quadrature and frequency grids

Noise-spectrum and filter-function calculations require stated integration intervals, infrared and ultraviolet cutoffs, grid spacing, and quadrature rules. Apparent convergence on a fixed window can hide cutoff dependence.

Optimization

An optimizer’s termination flag is not a physics validation. Repeat from multiple initial guesses, compare with a baseline control, enforce amplitude and bandwidth constraints, and test the final pulse on a finer propagation grid than the optimization grid.

Statistical Error in Trajectory Simulations

For $N$ independent trajectories with conditional states $\rho_c^{(n)}(t)$, the ensemble estimator is

$$\bar\rho_N(t) = \frac{1}{N} \sum_{n=1}^{N} \rho_c^{(n)}(t).$$

For an observable $A$, define samples

$$X_n(t) = \operatorname{Tr} \left[ A\rho_c^{(n)}(t) \right].$$

The estimated standard error of the sample mean is

$$\operatorname{SE} \left[ \bar X_N(t) \right] = \frac{s_X(t)}{\sqrt N},$$

where $s_X$ is the sample standard deviation. Quadrupling the number of independent trajectories reduces ordinary Monte Carlo error by only a factor of two.

For waiting-time distributions or rare jumps, mean curves can converge before tails do. Report binning rules, censoring, number of events, and confidence intervals. Reusing correlated random streams across nominally independent trajectories invalidates the usual standard-error estimate.

Records and Stochastic Conventions

Diffusive equations must state whether they use Itō or Stratonovich calculus. For an Itō Wiener increment,

$$\mathbb E[dW_t]=0, \qquad dW_t^2=dt.$$

A discrete implementation should verify that normalized increments have approximately zero mean, variance $dt$, and negligible unintended temporal correlations.

For a counting process,

$$dN_t\in\{0,1\}, \qquad dN_t^2=dN_t.$$

The time step must keep the probability of two or more unresolved jumps negligible when the algorithm assumes at most one. Event-driven methods remove that particular restriction but introduce their own root-finding and propagation tolerances.

The random seed belongs in run metadata, not as a substitute for uncertainty analysis. A trustworthy result should not depend qualitatively on one unusually favorable seed.

Cross-Method Validation

Whenever feasible, compare representations that fail differently.

Calculation	Primary method	Independent comparison
channel action	Kraus sum	superoperator or Choi reconstruction
Lindblad dynamics	ODE integration	matrix exponential for a small system
steady state	long-time propagation	null space of the Liouvillian
jump ensemble	Monte Carlo trajectories	direct master-equation solution
diffusive ensemble	stochastic integration	unconditional dephasing solution
non-Markovian toy model	enlarged unitary evolution	analytic one-excitation solution
coherence envelope	numerical spectral integral	white-noise or quasistatic limit
control pulse	optimization grid	finer independent propagation grid
work distribution	sampled histogram	exact transition-probability sum

Agreement between two methods is strongest when they do not share the same implementation path. Two wrappers around the same erroneous matrix construction are not independent validation.

Reproducibility Metadata

Each accepted run should record:

• a stable identifier for the notebook source;
• the execution environment and dependency versions;
• operating-system and hardware details when numerically relevant;
• all physical parameters with units;
• basis, tensor ordering, and vectorization convention;
• solver, tolerances, step restrictions, and truncations;
• random-number generator and seeds;
• optimizer configuration and initial guesses;
• input-data checksums where external data are used;
• run timestamp and elapsed time;
• generated data files separate from display figures;
• validation-test results and acceptance thresholds.

The Reproducibility Status page owns promotion labels across the site. A notebook specification should not silently upgrade itself from planned to reproduced because it ran once on one machine.

Data and Figure Integrity

Plots should be generated from saved numerical data rather than being the only output. Axes need units, parameter values, and normalization conventions. Logarithmic plots must expose zeros, clipped values, and sign handling.

For each figure, preserve enough information to answer:

1. Which run generated it?
2. Which data columns were plotted?
3. Which transformations or smoothing operations were applied?
4. Which validation tests passed?
5. Can the figure be regenerated without manual editing?

Interpolation may aid visualization but must not be confused with improved simulation resolution. Smoothing a noisy trajectory does not increase the detector bandwidth or trajectory count.

Acceptance Levels

Level	Evidence	Permitted claim
specification	model and required tests are documented	the calculation is planned and auditable
executable	code runs in a recorded environment	the implementation executes
validated	analytic, structural, and convergence tests pass	baseline outputs are numerically supported
reproduced	an independent clean run regenerates accepted outputs	the artifact is reproduced under stated conditions
benchmarked	independent methods or trusted data agree within tolerance	performance or accuracy comparison is supported

These levels describe evidence, not scientific importance. A small validated two-level calculation is more trustworthy than an elaborate untested simulation.

Diagnosing Failures

Symptom	Likely causes	First checks
trace drift	wrong generator, loose integration, vectorization mismatch	evaluate trace derivative and reduce step size
negative populations	invalid equation, unstable method, excessive step	inspect generator and convergence under refinement
wrong steady state	sign, rate, or basis error	solve the analytic rate balance and Liouvillian null space
trajectory average mismatch	incorrect jump probability or stochastic drift	test one-step expectation and record moments
cutoff dependence	inadequate Hilbert space	inspect boundary occupation and increase truncation
false revival	finite-size recurrence or sampling noise	vary environment size, cutoff, and ensemble count
optimizer reports implausible fidelity	objective or propagator mismatch	recompute with an independent fine-grid solver
fluctuation relation fails	work-sign, partition-function, or reverse-protocol error	test normalization and convention reversal

Failure is useful information. It should remain visible in validation logs rather than being removed from the final notebook.

Canonical Boundaries

This chapter owns computational specifications and validation workflows for the open-system examples listed here.

• Quantum Channels and Noise owns channel definitions and derivations.
• Markovian Master Equations owns GKSL dynamics and physical interpretation.
• Continuous Measurement and Quantum Trajectories owns stochastic measurement theory.
• Non-Markovian Dynamics owns memory diagnostics and their limitations.
• Quantum Control and Feedback owns control theory and objective design.
• Quantum Thermodynamics owns definitions of work, heat, entropy production, and fluctuation relations.
• Notebooks owns the site-wide notebook registry.
• Validation Tests owns site-wide validation policy.

Notebook pages should link to canonical theory rather than becoming alternate derivations of the same result.

Common Mistakes

• Treating a successful execution as validation.
• Omitting basis, tensor-order, vectorization, or unit conventions.
• Testing positivity only on a few hand-picked input states.
• Renormalizing an invalid state and hiding the original residual.
• Reporting solver tolerances without a refinement study.
• Choosing a Hilbert-space cutoff from convenience rather than convergence.
• Showing a trajectory ensemble without statistical uncertainty.
• Using one random seed as evidence of robustness.
• Fitting a decoherence time without reporting the fit window and model.
• Optimizing and validating a control pulse on the same coarse grid.
• Comparing methods that share the same faulty matrix construction.
• Saving figures without the numerical data and run metadata that generated them.
• Describing an unpromoted notebook specification as a reproduced artifact.

Exercises

Trace-preservation residual

A qubit amplitude-damping channel uses

$$\begin{aligned} K_0 &= \begin{pmatrix} 1&0\\ 0&\sqrt{1-p} \end{pmatrix}, \\[4pt] K_1 &= \begin{pmatrix} 0&\sqrt p\\ 0&0 \end{pmatrix}. \end{aligned}$$

Verify analytically that the trace-preservation residual vanishes for $0\le p\le1$.

Solution

Direct multiplication gives

$$K_0^\dagger K_0 = \begin{pmatrix} 1&0\\ 0&1-p \end{pmatrix},$$

and

$$K_1^\dagger K_1 = \begin{pmatrix} 0&0\\ 0&p \end{pmatrix}.$$

Their sum is $I$, so

$$\epsilon_{\mathrm{TP}} = \left\| K_0^\dagger K_0 + K_1^\dagger K_1 -I \right\| =0.$$

The condition $0\le p\le1$ also keeps the square roots real and gives the physical damping family.

Time-step order estimate

An observable at a fixed time is $O_h=0.7310$, $O_{h/2}=0.7270$, and $O_{h/4}=0.7260$. Assuming the leading error is proportional to $h^q$, estimate the observed order $q$.

Solution

The refinement differences are

$$\Delta_h = |O_h-O_{h/2}| =0.0040$$

and

$$\Delta_{h/2} = |O_{h/2}-O_{h/4}| =0.0010.$$

For leading error $Ch^q$,

$$\frac{\Delta_h}{\Delta_{h/2}} \approx 2^q.$$

The ratio is $4$, so $q\approx2$. This estimate is meaningful only if the three resolutions are already in the asymptotic convergence regime.

Trajectory count

A trajectory estimate has standard error $0.02$ using $N=2500$ independent trajectories. Approximately how many trajectories are required for standard error $0.005$ if ordinary Monte Carlo scaling applies?

Solution

Since standard error scales as $N^{-1/2}$, reducing it by a factor of four requires increasing $N$ by a factor of sixteen:

$$N' = 16N = 40000.$$

This estimate assumes independent samples and approximately unchanged variance. Rare-event observables may converge more slowly in practice.

Bosonic cutoff check

A damped oscillator calculation at cutoff $n_{\max}=20$ places probability $3\times10^{-3}$ in level $n=19$. At cutoff $n_{\max}=30$, the observable of interest changes by $4\%$. Is the smaller cutoff validated?

Solution

No. The non-negligible boundary population already warns that the state reaches the truncation edge, and the $4\%$ observable shift under refinement confirms material cutoff dependence. The cutoff should be increased until both the boundary diagnostic and the observable change satisfy a predeclared tolerance.

Independent validation

Why is comparing a matrix-exponential solver with an ODE solver useful only if their Liouvillian construction is also tested independently?

Solution

Both propagators may act perfectly on the same incorrectly constructed matrix. Their agreement would then validate numerical exponentiation and integration of that matrix, not its correspondence with the intended operator equation. Testing the Liouvillian action against direct operator evaluation separates representation errors from propagation errors.

Cross-Links

• Simulating Quantum Channels
• Bloch Vector Noise Models
• Solving Lindblad Equations
• Quantum Jump Simulation
• Diffusive Trajectory Simulation
• Non-Markovian Toy Models
• Decoherence Timescale Estimation
• Optimal Control Toy Problems
• Quantum Thermodynamics Toy Models
• Notebook Index
• Notebooks
• Validation Tests
• Reproducibility Status

References

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM (2002).
• E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd rev. ed., Springer (1993).
• E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd rev. ed., Springer (1996).
• H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press (2002).
• H. J. Carmichael, An Open Systems Approach to Quantum Optics, Springer (1993).
• J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP: An open-source Python framework for the dynamics of open quantum systems,” Computer Physics Communications 183, 1760–1772 (2012).
• J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP 2: A Python framework for the dynamics of open quantum systems,” Computer Physics Communications 184, 1234–1240 (2013).
• G. K. Sandve, A. Nekrutenko, J. Taylor, and E. Hovig, “Ten simple rules for reproducible computational research,” PLoS Computational Biology 9, e1003285 (2013).