Lindblad Master Equation Solvers

A Lindblad solver evolves a prescribed Markovian model while controlling integration and Hilbert-space truncation errors. Trace preservation and positivity are necessary checks, but neither proves that the chosen master equation describes the experiment. The Lindblad–GKSL equation explains the physical assumptions behind this model.

Helpful background. Vectorization constructs the generator; density-matrix diagnostics supplies state checks.

The Liouvillian initial-value problem

For collapse operators $L_\mu$ containing square-root rates,

$$\dot\rho=\mathcal L(\rho) =-\frac{i}{\hbar}[H,\rho] +\sum_\mu\left(L_\mu\rho L_\mu^\dagger -\frac12\{L_\mu^\dagger L_\mu,\rho\}\right).$$

A rate $\gamma$ belongs in $L=\sqrt\gamma\,A$, not $\gamma A$. In finite dimension with a time-independent generator, $\rho(t)=e^{t\mathcal L}\rho(0)$. An exponential action can avoid forming the full propagator. For time-dependent coefficients an adaptive ODE method evaluates the current generator at each stage.

Use sparse action or direct matrix products when an explicit $d^2\times d^2$ Liouvillian is too large. Multiple time scales create stiffness: fast decays can constrain an explicit integrator even when only slow populations are desired. Compare an implicit method or exponential method before tightening explicit tolerances indefinitely.

A solvable relaxation benchmark

Let $H=\hbar\omega|1\rangle\langle1|$ and $L=\sqrt\gamma|0\rangle\langle1|$. Then

$$\begin{aligned} \rho_{11}(t)&=\rho_{11}(0)e^{-\gamma t},\\ \rho_{01}(t)&=\rho_{01}(0)e^{(i\omega-\gamma/2)t},\\ \rho_{00}(t)&=1-\rho_{11}(t). \end{aligned}$$

Start from $(|0\rangle+i|1\rangle)/\sqrt2$ to test both coherent rotation and damping. Compare populations, real and imaginary coherences, trace, and the smallest eigenvalue over the entire time interval. A terminal population alone cannot detect a wrong Hamiltonian sign.

Forward Euler illustrates why a trace-preserving derivative is insufficient. Set $H=0$ and $x=\gamma\Delta t>0$. Euler damping multiplies the excited population by $1-x$ and the coherences by $1-x/2$. In the input-first Choi matrix, the corner block is

$$\begin{pmatrix}1&1-x/2\\1-x/2&1-x\end{pmatrix}, \qquad \det=-\frac{x^2}{4}.$$

Thus the step is not completely positive for any positive $x$, even when selected input states remain positive. A particular test state need not reveal the defect at every step size: the equal superposition has positive output for sufficiently small damping-only steps. Test the finite-step map or use convergence and state diagnostics together.

For this qubit, exact damping with $p=1-e^{-\gamma\Delta t}$ is a Kraus channel. With a time-independent noncommuting drive, composing a unitary half-step, an exact damping step, and a unitary half-step preserves complete positivity and trace at every step. Its global accuracy is still second order for smooth bounded generators; compare with a refined reference. Generic Runge–Kutta integration does not have this exact channel property.

Steady states as constrained null vectors

Solve

$$\mathsf Lr=0,\qquad \operatorname{vec}(I)^\dagger r=1.$$

One practical formulation replaces one redundant equation by the trace constraint; another augments the linear system. After reshaping, check the residual $\|\mathcal L(\rho_{\rm ss})\|$, Hermiticity, trace, and positivity. A normalized vector from a near-zero eigenvalue is not automatically a physical steady state.

The replaced equation must really be redundant, and the constrained system must have full rank. Inspect its conditioning when the relaxation gap is small. The steady-state residual alone does not bound the state error uniformly: on the trace-zero subspace, the inverse of the restricted generator determines how residuals are amplified. For a nonnormal generator, an eigenvalue gap alone is not a sufficient conditioning estimate.

If the kernel has dimension greater than one, the steady state is not unique. For pure dephasing, every diagonal density matrix is stationary. Initial conditions and conserved quantities determine the long-time state, and a single trace constraint cannot select it. Near degeneracy can also make steady-state observables sensitive to solver tolerances.

A small residual with a large state error

A finite decay chain makes the conditioning issue explicit. Take $N+1$ orthonormal levels $|0\rangle,\ldots,|N\rangle$, $H=0$, and the separate collapse operators

$$L_j=\sqrt\gamma\,|j-1\rangle\langle j|, \qquad j=1,\ldots,N,\qquad \gamma>0.$$

Every excited level feeds the one below it at the same rate. The unique stationary state is $\rho_{\rm ss}=|0\rangle\langle0|$. Consider the positive, normalized trial state

$$\tilde\rho=\frac1N\sum_{j=1}^N|j\rangle\langle j|, \qquad \mathcal L(\tilde\rho)=\frac\gamma N \left(|0\rangle\langle0|-|N\rangle\langle N|\right).$$

All interior population flows cancel in its residual. Nevertheless,

$$\frac12\|\tilde\rho-\rho_{\rm ss}\|_1=1, \qquad \|\mathcal L(\tilde\rho)\|_1=\frac{2\gamma}{N}.$$

The residual becomes small relative to the rate $\gamma$ as the dimension grows, while a ground-state measurement distinguishes the trial state perfectly from the answer. All nonzero population eigenvalues are $-\gamma$; ground–excited coherences decay at $\gamma/2$ and other coherences at $\gamma$. Thus the smallest nonzero decay rate of the full generator stays $\gamma/2$. For $N>1$, its population block contains a nontrivial Jordan chain, so that fixed eigenvalue gap does not control the inverse norm. On trace-zero operators the inverse amplifies this particular residual by $N/\gamma$ in trace norm.

For one fixed finite model the inverse bound remains finite. The example rules out an error certificate based only on a small residual and an eigenvalue gap across different models or increasing cutoffs. Supplement residuals with conditioning estimates and independently converged target observables.

A convergence protocol

First hold the Hilbert cutoff fixed and tighten time tolerances or refine the step until target observables stabilize. Then increase the cutoff independently, monitoring boundary occupation. Finally compare one small instance with exact diagonalization of the Liouvillian or a channel exponential. Keep absolute and relative tolerances meaningful for small coherences as well as order-one populations.

For an oscillator with $L=\sqrt\kappa\,a$, the vacuum must be stationary and $\langle n(t)\rangle=\langle n(0)\rangle e^{-\kappa t}$ when $H=\hbar\omega n$. Failure near a high-occupation cutoff can be a representation error rather than an ODE failure. See bosonic bases.

Exercise

Add excitation $L_+=\sqrt{\gamma_+}|1\rangle\langle0|$ to relaxation $L_-=\sqrt{\gamma_-}|0\rangle\langle1|$. Find the stationary excited population and population decay rate.

Solution

$\dot p_1=\gamma_+(1-p_1)-\gamma_-p_1$, so $p_1^{\rm ss}=\gamma_+/(\gamma_++\gamma_-)$ and the approach rate is $\gamma_++\gamma_-$. If both rates vanish there is no unique stationary population. The example tests upward/downward rate placement and the singular zero-rate limit.

For the decay chain above, find the mean time to reach $|0\rangle$ when the initial state is $\tilde\rho$. Explain why it can grow even though the decay eigenvalue gap stays fixed.

Solution

Starting at level $j$ requires $j$ successive jumps, each with mean waiting time $1/\gamma$, so the mean absorption time is $j/\gamma$. Averaging the uniform initial populations gives $(N+1)/(2\gamma)$. The cascade must transport population through many levels; its evolution contains polynomial factors multiplying decaying exponentials. A fixed eigenvalue gap alone does not bound this transient duration uniformly in $N$.

References

• Lindblad, G., “On the generators of quantum dynamical semigroups,” Communications in Mathematical Physics 48, 119–130 (1976), doi:10.1007/BF01608499.
• QuTiP developers, “Lindblad Master Equation Solver” and “Steady State Solvers”, QuTiP documentation — numerical implementations and solver choices.