Solving Lindblad Equations

This notebook guide specifies a reproducible finite-dimensional calculation for solving Lindblad–GKSL master equations. The goal is not only to integrate an ordinary differential equation, but to verify the generator, compare solver representations, compute stationary states, inspect relaxation modes, and connect continuous-time dynamics with finite-time quantum channels.

As of this review, no executable notebook under notebooks/density-open-systems/lindblad-solvers/ is promoted as a reproduced artifact. This page is the admission contract for that notebook family: it states what the notebook should compute, what conventions it must declare, and which checks are required before any numerical output is cited.

Purpose

The notebook should demonstrate how to:

• build a finite-dimensional Lindblad generator from $H$, $L_k$, and rates;
• convert the generator into a matrix-vectorized Liouvillian;
• compare direct density-matrix right-hand sides with the vectorized Liouvillian;
• compare ODE integration with a matrix exponential for time-independent generators;
• compute steady states from the null space of the Liouvillian;
• inspect Liouvillian eigenvalues as decay rates and oscillation frequencies;
• validate trace preservation, Hermiticity, positivity, and complete positivity of finite-time maps;
• compare numerical results with analytic qubit models.

The first version should stay deliberately small: qubits and low-dimensional truncated oscillators are enough to expose the numerical issues.

Directory Plan

Use a dedicated directory:

notebooks/density-open-systems/lindblad-solvers/
  solving-lindblad-equations.ipynb
  README.md

The opening notebook cell or README.md should state:

• Python and package versions;
• basis ordering;
• density-matrix vectorization convention;
• Hamiltonian convention and units;
• whether rates are included in $L_k$ or kept as separate $\gamma_k$;
• ODE solver method and tolerances;
• norm used for residuals;
• random seed, if random states or random Hamiltonians are used;
• date and commit identifier when the notebook is promoted.

Mathematical Convention

Use the Lindblad–GKSL equation in the same convention as Lindblad–GKSL Equation:

$$\frac{d\rho}{dt} = - \frac{i}{\hbar}[H,\rho] + \sum_k\gamma_k \mathcal D[L_k]\rho, \qquad \gamma_k\ge0,$$

with

$$\mathcal D[L]\rho = L\rho L^\dagger - \frac12 \{L^\dagger L,\rho\}.$$

The notebook should decide in one place whether to represent each dissipator as $\gamma_k\mathcal D[L_k]$ or as $\mathcal D[\sqrt{\gamma_k}L_k]$. Both are standard, but mixing them is a common source of factor-of-two and factor-of-rate errors.

Column-Stacking Vectorization

For a $d\times d$ matrix $X$, define $\operatorname{vec}(X)$ by stacking columns:

$$\operatorname{vec}(X)_{a+bd} = X_{ab}, \qquad a,b=0,\ldots,d-1.$$

With this convention,

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

Therefore the time-independent master equation can be written as

$$\frac{d}{dt}\operatorname{vec}(\rho) = \mathbb L\operatorname{vec}(\rho),$$

where

$$\begin{aligned} \mathbb L =& - \frac{i}{\hbar} \left( I\otimes H - H^{\mathsf T}\otimes I \right) \\ &+ \sum_k\gamma_k \left[ L_k^*\otimes L_k - \frac12 I\otimes L_k^\dagger L_k - \frac12 (L_k^\dagger L_k)^{\mathsf T}\otimes I \right]. \end{aligned}$$

This formula should be tested against a direct function that computes $\dot\rho$ from matrix multiplications. The notebook should not rely only on a vectorization formula that has not been validated.

Baseline Models

Implement at least three analytically solvable generators. A later extension should add the driven two-level benchmark from Optical Bloch Equations to test detuning signs, saturation, and complex Liouvillian eigenvalues.

Pure dephasing qubit

Use

$$H = \frac{\hbar\omega_0}{2}\sigma_z,$$

and

$$\frac{d\rho}{dt} = - \frac{i}{\hbar}[H,\rho] + \frac{\Gamma_\phi}{2} (\sigma_z\rho\sigma_z-\rho).$$

The analytic solution is

$$\begin{aligned} \rho_{00}(t)&=\rho_{00}(0),\\ \rho_{11}(t)&=\rho_{11}(0),\\ \rho_{01}(t)&= e^{-i\omega_0t} e^{-\Gamma_\phi t} \rho_{01}(0). \end{aligned}$$

This model tests phase conventions, Hamiltonian signs, and coherence decay. Use Pure Dephasing Master Equation as the analytic reference.

Zero-temperature amplitude damping

Use

$$\sigma_-=\lvert g\rangle\langle e\rvert, \qquad \sigma_+=\lvert e\rangle\langle g\rvert,$$

and

$$\frac{d\rho}{dt} = - \frac{i}{\hbar}[H,\rho] + \Gamma_1\mathcal D[\sigma_-]\rho.$$

In the interaction picture,

$$\rho_{ee}(t) = e^{-\Gamma_1t}\rho_{ee}(0), \qquad \rho_{eg}(t) = e^{-\Gamma_1t/2}\rho_{eg}(0).$$

The finite-time channel has probability

$$p(t)=1-e^{-\Gamma_1t}.$$

This model tests nonunital relaxation and the relation between a master equation and a finite-time Amplitude-Damping Channel.

Finite-temperature two-level system

Use

$$\frac{d\rho}{dt} = - \frac{i}{\hbar}[H,\rho] + \Gamma_\downarrow\mathcal D[\sigma_-]\rho + \Gamma_\uparrow\mathcal D[\sigma_+]\rho.$$

The excited-state population satisfies

$$\dot p_e = - \Gamma_\downarrow p_e + \Gamma_\uparrow(1-p_e),$$

so the steady-state value is

$$p_e^{\mathrm{ss}} = \frac{\Gamma_\uparrow} {\Gamma_\uparrow+\Gamma_\downarrow}.$$

For a thermal bath, detailed balance additionally requires

$$\frac{\Gamma_\uparrow}{\Gamma_\downarrow} = e^{-\beta\hbar\omega_0}.$$

This model tests steady states and rate-equation reductions. Use Pauli Rate Equations, Thermal Master Equations, and Detailed Balance for the physical assumptions.

Optional truncated oscillator

A harmonic oscillator with loss is often written

$$\dot\rho = - \frac{i}{\hbar}[H,\rho] + \kappa(\bar n+1)\mathcal D[a]\rho + \kappa\bar n\,\mathcal D[a^\dagger]\rho.$$

If included, the notebook must state the Fock-space cutoff and test cutoff convergence. A truncated oscillator is finite-dimensional numerically but only approximates the infinite-dimensional model.

Numerical Workflow

The notebook should follow a validation-first workflow:

1. Define basis ordering and matrix constructors.
2. Define a direct function for $\mathcal L(\rho)$.
3. Define column-stacking vectorization and its inverse.
4. Build the matrix Liouvillian $\mathbb L$.
5. Compare $\operatorname{vec}(\mathcal L(\rho))$ with $\mathbb L\operatorname{vec}(\rho)$ on fixed test states.
6. Integrate the direct ODE in density-matrix form.
7. Integrate the vectorized ODE.
8. Compare both ODE results with $e^{t\mathbb L}\operatorname{vec}(\rho_0)$ for time-independent generators.
9. Compute the steady state from the null space of $\mathbb L$.
10. Inspect Liouvillian eigenvalues and identify relaxation modes.
11. Build finite-time maps at selected times and test them as quantum channels.

The first version should avoid random Hamiltonians as the main evidence. Fixed small systems make sign, transpose, and ordering errors easier to find.

Test States

Use the same small density-matrix set as Simulating Quantum Channels:

$$\rho_0=\lvert0\rangle\langle0\rvert, \qquad \rho_1=\lvert1\rangle\langle1\rvert,$$ $$\rho_+ = \lvert+\rangle\langle+\rvert, \qquad \lvert+\rangle = \frac{\lvert0\rangle+\lvert1\rangle}{\sqrt2},$$

and

$$\rho_y = \lvert +y\rangle\langle +y\rvert, \qquad \lvert +y\rangle = \frac{\lvert0\rangle+i\lvert1\rangle}{\sqrt2}.$$

Also test the maximally mixed state $I/2$. For amplitude damping it should move; for pure dephasing it should not.

Solver Comparisons

For time-independent generators, compare three equivalent computations:

$$\rho_{\mathrm{direct}}(t) \quad \rho_{\mathrm{vec}}(t) \quad \rho_{\mathrm{expm}}(t).$$

Here:

• $\rho_{\mathrm{direct}}(t)$ comes from an ODE solver acting on real and imaginary parts of the density matrix or on a complex state vector if the solver supports it;
• $\rho_{\mathrm{vec}}(t)$ comes from the vectorized ODE $\dot y=\mathbb L y$;
• $\rho_{\mathrm{expm}}(t)$ comes from $y(t)=e^{t\mathbb L}y(0)$.

The matrix exponential is not always the most efficient method for large systems, but it is an excellent validation tool for small time-independent examples.

For time-dependent Hamiltonians or rates, the notebook should not use $e^{t\mathbb L}$ unless $\mathbb L$ is constant. A time-ordered or stepwise method is needed when the generator changes during the interval.

Steady States

A steady state satisfies

$$\mathcal L(\rho_{\mathrm{ss}})=0, \qquad \operatorname{Tr}\rho_{\mathrm{ss}}=1, \qquad \rho_{\mathrm{ss}}\ge0.$$

For the conceptual distinctions among fixed points, attractors, conserved sectors, gaps, and metastable modes, see Steady States and Relaxation.

In vectorized form,

$$\mathbb L\,\operatorname{vec}(\rho_{\mathrm{ss}})=0.$$

The notebook should compute the null space and then impose trace normalization. For a unique steady state, the nullity should be one within numerical tolerance. If the nullity is larger, the model has multiple stationary states or conserved sectors, and the notebook should not silently pick an arbitrary vector. This null-space check is also a basic validation step for Reservoir Engineering, where uniqueness and the Liouvillian gap decide whether the intended target is actually stabilized.

For zero-temperature amplitude damping, the steady state is

$$\rho_{\mathrm{ss}} = \lvert g\rangle\langle g\rvert.$$

For the finite-temperature two-level system,

$$\rho_{\mathrm{ss}} = \begin{pmatrix} 1-p_e^{\mathrm{ss}}&0\\ 0&p_e^{\mathrm{ss}} \end{pmatrix}$$

in the $(\lvert g\rangle,\lvert e\rangle)$ basis.

Liouvillian Spectrum

The Liouvillian spectrum gives the relaxation-mode structure. If

$$\mathbb L v_\alpha=\lambda_\alpha v_\alpha,$$

then the corresponding mode evolves as

$$e^{\lambda_\alpha t}.$$

For a stable finite-dimensional Markovian semigroup, nonstationary modes should have

$$\operatorname{Re}\lambda_\alpha\le0.$$

The zero eigenvalues correspond to stationary operators. Negative real parts give decay rates. Nonzero imaginary parts give coherent oscillations or damped oscillatory modes.

The notebook should report:

• the eigenvalue closest to zero;
• the number of eigenvalues with magnitude below tolerance;
• the largest real part;
• the spectral gap when a unique steady state exists;
• whether eigenvectors reconstruct Hermitian modes or only complex mode pairs.

Eigenvectors of a nonnormal Liouvillian can be ill conditioned. The spectrum is useful, but residuals and physical checks remain necessary.

Finite-Time Channel Checks

For a time-independent generator, the finite-time map is

$$\Phi_t=e^{t\mathcal L}.$$

The notebook should build the superoperator matrix $e^{t\mathbb L}$ and convert it into the same Choi convention used by Choi Matrix. Then it should test:

$$\operatorname{Tr}\Phi_t(\rho)=\operatorname{Tr}\rho, \qquad J_{\Phi_t}\ge0.$$

For a correct Lindblad generator, $J_{\Phi_t}$ should be positive semidefinite for $t\ge0$, up to roundoff. This is a stronger check than testing positivity on a few system states.

The notebook should also compare special cases with known channels:

$$\lambda_{\mathrm{deph}}(t)=e^{-\Gamma_\phi t}, \qquad p_{\mathrm{amp}}(t)=1-e^{-\Gamma_1t}.$$

These comparisons link the master-equation notebook to the finite-channel notebook.

Validation Checks

The notebook should report pass/fail checks with numerical tolerances.

Check	Requirement
Direct versus vectorized generator	$\lVert\operatorname{vec}(\mathcal L(\rho))-\mathbb L\operatorname{vec}(\rho)\rVert\lt\epsilon$
Trace preservation of generator	$\lvert\operatorname{Tr}\mathcal L(\rho)\rvert\lt\epsilon$
Hermiticity preservation	$\lVert\mathcal L(\rho)-\mathcal L(\rho)^\dagger\rVert\lt\epsilon$ for Hermitian $\rho$
ODE versus exponential	$\lVert\rho_{\mathrm{ode}}(t)-\rho_{\mathrm{expm}}(t)\rVert\lt\epsilon$
Analytic solution	matrix elements match the linked model page
Steady-state residual	$\lVert\mathcal L(\rho_{\mathrm{ss}})\rVert\lt\epsilon$
Steady-state trace	$\lvert\operatorname{Tr}\rho_{\mathrm{ss}}-1\rvert\lt\epsilon$
Steady-state positivity	smallest eigenvalue $\ge-\epsilon$
Finite-time trace preservation	Choi partial-trace residual below tolerance
Finite-time complete positivity	Choi eigenvalues nonnegative within tolerance

Use a tolerance appropriate to double precision, such as $\epsilon=10^{-10}$ for ODE comparisons and $\epsilon=10^{-12}$ for exact algebraic identities, while keeping all tolerances declared in one place.

Expected Output

The notebook should produce tables rather than only plots:

• model name, parameters, and units;
• direct-versus-vectorized residuals;
• ODE-versus-exponential residuals at selected times;
• analytic-versus-numerical residuals;
• trace and minimum-eigenvalue checks for evolved states;
• steady-state residuals and null-space dimensions;
• Liouvillian eigenvalues sorted by real part;
• finite-time Choi eigenvalue minima.

Plots are helpful, especially for population relaxation and Bloch-vector trajectories, but the validation table is the evidence.

Suggested Parameter Sets

Use small parameter grids that include limiting cases:

$$\omega_0 \in \{0,1\}, \qquad \Gamma_\phi \in \{0,0.2,1\},$$ $$\Gamma_1 \in \{0,0.5,2\}, \qquad t \in \{0,0.1,1,5\}.$$

For finite-temperature damping, use one equilibrium example such as

$$\Gamma_\downarrow=1, \qquad \Gamma_\uparrow=0.25,$$

and verify that the steady excited-state population is $p_e^{\mathrm{ss}}=0.2$.

Common Numerical Failures

• Using row-major vectorization while applying column-stacking formulas.
• Forgetting the transpose in $\operatorname{vec}(AXB)=(B^{\mathsf T}\otimes A)\operatorname{vec}(X)$.
• Double-counting rates by using both $\gamma_k$ and $\sqrt{\gamma_k}L_k$.
• Treating an ODE solver’s small trace drift as physical.
• Clipping negative eigenvalues without reporting the original residual.
• Calling a nonunique null space a unique steady state.
• Using $e^{t\mathbb L}$ for a time-dependent generator.
• Testing positivity on a few states but never testing complete positivity of the finite-time map.
• Forgetting that a truncated oscillator requires cutoff checks.
• Reporting only visual trajectories, with no residual table.

Exercises

Vectorization identity

Using column stacking, show that $\operatorname{vec}(AXB)=(B^{\mathsf T}\otimes A)\operatorname{vec}(X)$.

Solution

Write

$$(AXB)_{ij} = \sum_{m,n} A_{im}X_{mn}B_{nj}.$$

In column stacking, the pair $(m,n)$ is mapped to the single index $m+nd$. The matrix $B^{\mathsf T}\otimes A$ has entries

$$(B^{\mathsf T}\otimes A)_{i+jd,\,m+nd} = B_{nj}A_{im}.$$

Multiplying by $\operatorname{vec}(X)_{m+nd}=X_{mn}$ gives the same sum:

$$\sum_{m,n} B_{nj}A_{im}X_{mn} = (AXB)_{ij}.$$

Pure dephasing modes

For the pure dephasing qubit with Hamiltonian $H=(\hbar\omega_0/2)\sigma_z$, what are the eigenvalues associated with the operators $\lvert0\rangle\langle1\rvert$ and $\lvert1\rangle\langle0\rvert$?

Solution

The coherences obey

$$\dot\rho_{01} = -(i\omega_0+\Gamma_\phi)\rho_{01},$$

and

$$\dot\rho_{10} = (i\omega_0-\Gamma_\phi)\rho_{10}.$$

Thus $\lvert0\rangle\langle1\rvert$ has Liouvillian eigenvalue

$$\lambda_{01}=-\Gamma_\phi-i\omega_0,$$

and $\lvert1\rangle\langle0\rvert$ has eigenvalue

$$\lambda_{10}=-\Gamma_\phi+i\omega_0.$$

They form a complex-conjugate pair, as expected for Hermiticity-preserving dynamics.

Thermal steady state

For a two-level system with $\Gamma_\downarrow=1$ and $\Gamma_\uparrow=0.25$, compute the steady excited-state probability.

Solution

The rate equation is

$$\dot p_e = - \Gamma_\downarrow p_e + \Gamma_\uparrow(1-p_e).$$

At stationarity,

$$0 = - (\Gamma_\downarrow+\Gamma_\uparrow)p_e^{\mathrm{ss}} + \Gamma_\uparrow.$$

Therefore

$$p_e^{\mathrm{ss}} = \frac{\Gamma_\uparrow} {\Gamma_\downarrow+\Gamma_\uparrow} = \frac{0.25}{1.25} = 0.2.$$

Trace-preserving left mode

Explain why a trace-preserving Liouvillian has a left zero mode corresponding to the identity operator.

Solution

Trace preservation means

$$\operatorname{Tr}\mathcal L(\rho)=0$$

for every $\rho$. With the Hilbert–Schmidt pairing, this is

$$\operatorname{Tr}\!\left[I\,\mathcal L(\rho)\right]=0.$$

Equivalently,

$$\mathcal L^\dagger(I)=0.$$

Thus the identity is a zero mode of the adjoint generator. In matrix-vectorized form, the vector representing the trace functional is a left eigenvector of $\mathbb L$ with eigenvalue zero.

Why check the Choi matrix?

If a numerical finite-time map sends the notebook’s five test states to positive density matrices, why is that not enough to prove complete positivity?

Solution

Positivity on a finite list of system states only checks those particular inputs. Complete positivity requires positivity after the map is extended by an identity operation on an arbitrary reference system. A map can look harmless on a small set of system states and still fail when applied to part of an entangled state.

For finite-dimensional channels, Choi positivity is the practical test:

$$J_{\Phi_t}\ge0.$$

That is why the notebook should construct the finite-time Choi matrix instead of relying only on evolved sample states.

Cross-Links

• Noise Simulation
• Lindblad–GKSL Equation
• Lindblad Operators
• Quantum Dynamical Semigroups
• Pure Dephasing Master Equation
• Amplitude Damping Master Equation
• Optical Bloch Equations
• Pauli Rate Equations
• Steady States and Relaxation
• Thermal Master Equations
• Simulating Quantum Channels
• Formula Sheet
• Approximation Checklist
• Validation Tests

References

• V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, “Completely positive dynamical semigroups of N-level systems,” Journal of Mathematical Physics 17, 821–825 (1976).
• G. Lindblad, “On the generators of quantum dynamical semigroups,” Communications in Mathematical Physics 48, 119–130 (1976).
• H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press (2002).
• C. W. Gardiner and P. Zoller, Quantum Noise, 3rd ed., Springer (2004).
• H. J. Carmichael, Statistical Methods in Quantum Optics 1, Springer (1999).
• E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I, 2nd ed., Springer (1993).
• P. Virtanen et al., “SciPy 1.0: fundamental algorithms for scientific computing in Python,” Nature Methods 17, 261–272 (2020).
• 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).