Solving Lindblad Equations

A Lindblad equation is a linear differential equation for a density matrix, but a generic numerical step need not preserve positivity. This example solves a thermal qubit in two representations, compares both with an exact solution, and checks the stationary state. Halving the step size separates numerical error from physical decay.

Helpful background. The Lindblad equation defines the generator, and thermal master equations explain its rates. The complete example uses Python 3 and NumPy.

A thermal qubit with an exact answer

Use $(|0\rangle,|1\rangle)=(|g\rangle,|e\rangle)$, set $\hbar=1$, and take $H=\omega|1\rangle\langle1|$. Rates are included in the jump operators:

$$L_\downarrow=\sqrt{\gamma_\downarrow}|0\rangle\langle1|, \qquad L_\uparrow=\sqrt{\gamma_\uparrow}|1\rangle\langle0|.$$

Thus $\mathcal L(\rho)=-i[H,\rho]+\sum_k\mathcal D[L_k]\rho$, with no extra rate multiplying a dissipator. For $G=\gamma_\downarrow+\gamma_\uparrow>0$,

$$\begin{aligned} p_e(t)&=p_*+[p_e(0)-p_*]e^{-Gt}, & p_*&=\frac{\gamma_\uparrow}{G},\\ \rho_{01}(t)&=e^{(i\omega-G/2)t}\rho_{01}(0). \end{aligned}$$

The positive phase of $\rho_{01}$ follows from $E_0-E_1=-\omega$. With $\gamma_\downarrow=1$, $\gamma_\uparrow=0.25$, and $\omega=1$, the stationary state is $\operatorname{diag}(0.8,0.2)$. These rates describe a thermal bath at $\beta\hbar\omega=\log4$ when the weak-coupling thermal model applies.

Build the Liouvillian without changing conventions

Column stacking means A.reshape(-1, order='F'). It obeys $\operatorname{vec}(AXB)=(B^{\mathsf T}\otimes A)\operatorname{vec}(X)$, so

$$\mathbb L={}-i(I\otimes H-H^{\mathsf T}\otimes I) +\sum_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].$$

A transpose and an adjoint are different operations here. Test this formula on matrix units, including non-Hermitian ones, before trusting a trajectory. The following fourth-order Runge–Kutta method works on either matrices or vectors because NumPy applies its additions componentwise.

import numpy as np

I = np.eye(2, dtype=complex)
H = np.diag([0.0, 1.0])
down, up = 1.0, 0.25
sm = np.array([[0, 1], [0, 0]], dtype=complex)
Ls = [np.sqrt(down)*sm, np.sqrt(up)*sm.conj().T]

def rhs(A):
    out = -1j*(H@A-A@H)
    for L in Ls:
        R = L.conj().T@L
        out += L@A@L.conj().T-(R@A+A@R)/2
    return out

LL = -1j*(np.kron(I, H)-np.kron(H.T, I))
for L in Ls:
    R = L.conj().T@L
    LL += np.kron(L.conj(), L)-np.kron(I, R)/2-np.kron(R.T, I)/2

vec = lambda A: A.reshape(-1, order='F')
mat = lambda v: v.reshape(2, 2, order='F')
for i in range(2):
    for j in range(2):
        E = np.zeros((2, 2), dtype=complex)
        E[i, j] = 1
        assert np.allclose(vec(rhs(E)), LL@vec(E))
assert np.linalg.norm(vec(I).conj()@LL) < 1e-12

def rk4(f, y, T, steps):
    y = y.copy()
    h = T/steps
    for _ in range(steps):
        a = f(y)
        b = f(y+h*a/2)
        c = f(y+h*b/2)
        d = f(y+h*c)
        y += h*(a+2*b+2*c+d)/6
    return y

rho0 = np.array([[0.5, -0.5j], [0.5j, 0.5]])
T = 2.0
pe = 0.2+0.3*np.exp(-1.25*T)
coh = -0.5j*np.exp((1j-0.625)*T)
exact = np.array([[1-pe, coh], [coh.conjugate(), pe]])
errors = []
for n in [20, 40, 80]:
    direct = rk4(rhs, rho0, T, n)
    vectorized = mat(rk4(lambda v: LL@v, vec(rho0), T, n))
    assert np.linalg.norm(direct-vectorized) < 1e-12
    assert abs(np.trace(direct)-1) < 1e-12
    assert np.linalg.norm(direct-direct.conj().T) < 1e-12
    assert np.linalg.eigvalsh(direct).min() >= -1e-12
    errors.append(np.linalg.norm(direct-exact))
    print(n, 'Frobenius error:', errors[-1])
assert 14 < errors[0]/errors[1] < 18
assert 14 < errors[1]/errors[2] < 18

# Add trace = 1 to the homogeneous stationary equation.
A = np.vstack([LL, vec(I).conj()])
b = np.array([0, 0, 0, 0, 1], dtype=complex)
rho_ss = mat(np.linalg.lstsq(A, b, rcond=None)[0])
assert np.allclose(rho_ss, np.diag([0.8, 0.2]))
print('stationary state:', rho_ss, sep='\n')
print('Liouvillian eigenvalues:', np.linalg.eigvals(LL))

The eigenvalues should be $0,-1.25,-0.625+i,-0.625-i$. The error decreases by approximately $2^4=16$ when the step is halved, as expected for this smooth fourth-order calculation. Agreement between direct and vectorized integration checks the representation; agreement with the analytic answer checks the propagation. Neither comparison alone covers both errors.

What a solver result does and does not establish

Runge–Kutta is not a completely positive integrator for arbitrary step size. Trace preservation can hold while positivity fails. Even for a purely Hamiltonian qubit, forward Euler sends a unit Bloch vector to one longer than unity after a perpendicular infinitesimal rotation. Do not repair such errors by silently clipping eigenvalues: reduce the step or change the numerical method, then compare results.

For small constant generators, a reliable matrix-exponential routine provides another comparison with $e^{t\mathbb L}$. Avoid constructing that exponential through eigenvector inversion when the Liouvillian is defective or ill conditioned. For time-dependent generators, a single exponential of $t\mathbb L(t)$ is generally wrong; integrate the changing generator or use a controlled time-ordered approximation.

The least-squares stationary calculation above works because this model has a unique stationary state. Pure dephasing has a whole family of diagonal stationary states. A numerical null-space basis is then a basis of stationary operators, not automatically a list of density matrices; trace and positivity still matter. See steady states and relaxation.

To extend the example to an oscillator, use both $\kappa(\bar n+1)\mathcal D[a]$ and $\kappa\bar n\mathcal D[a^\dagger]$. Increase the Fock cutoff separately from decreasing the integration step. These test truncation error and time-discretization error, respectively.

Exercises

Column stacking. Derive the vectorization identity used above.

Solution

The matrix element at row $i+jd$, column $m+nd$ of $B^{\mathsf T}\otimes A$ is $B_{nj}A_{im}$. Multiplying by $X_{mn}$ and summing gives $(AXB)_{ij}$. The transpose on $B$ is required by the index order.

Which coherence rotates which way? Change the Hamiltonian to $H=\omega Z/2$, keeping the same basis labels. What replaces the coherence phase?

Solution

Now $E_0-E_1=\omega$, so $\rho_{01}(t)$ contains $e^{-i\omega t}$. With this Hamiltonian $|0\rangle$ has the higher energy for $\omega>0$; keeping the original downward jump would no longer describe decay toward its ground state.

A left zero mode. Why does the code test vec(I).conj() @ LL?

Solution

The row vector computes the trace. Its product with $\mathbb L$ vanishes exactly when $\operatorname{Tr}\mathcal L(A)=0$ for every matrix $A$, or equivalently $\mathcal L^\dagger(I)=0$.

Predict an error. If the error at $80$ steps is $\epsilon$, estimate it at $160$ steps and explain when the estimate fails.

Solution

In the fourth-order convergence regime it is about $\epsilon/16$. This estimate can fail before the step is small enough, near nonsmooth time dependence, or after roundoff dominates. The ratio must be observed, not assumed from the method’s name.

References

• H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press (2002).
• 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). doi:10.1063/1.522979.
• E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd edition, Springer (1993).
• G. Lindblad, “On the generators of quantum dynamical semigroups,” Communications in Mathematical Physics 48, 119–130 (1976). doi:10.1007/BF01608499.