Steady States and Relaxation

A steady state of a time-independent Markovian master equation is a density operator that the generator leaves unchanged. If

$$\dot\rho = \mathcal L(\rho),$$

then a steady state satisfies

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

Steady states are fixed points of the finite-time evolution:

$$\Phi_t(\rho_{\mathrm{ss}}) = e^{t\mathcal L}\rho_{\mathrm{ss}} = \rho_{\mathrm{ss}}.$$

Finding a stationary state is only the beginning: its uniqueness, attracting basin, and relaxation rates determine what an experiment approaches. A dark subspace can protect information, while a slowly decaying mode can make a state look stationary before it changes again.

Helpful background. Quantum dynamical semigroups relate the generator to finite-time evolution. Solving Lindblad equations supplies the numerical methods. The distinct closed-system thermalization problem concerns restricted observations of globally unitary dynamics.

Liouvillian Kernel

Assume a time-independent finite-dimensional GKSL generator throughout the spectral discussion. The linear map $\mathcal L$ is often called the Liouvillian. Its kernel is

$$\ker\mathcal L = \{X:\mathcal L(X)=0\}.$$

Physical steady states are the positive, trace-one elements of this kernel. The kernel may also contain non-Hermitian or nonpositive operators; those are useful for linear algebra, but they are not states by themselves.

In finite dimensions, a time-independent completely positive trace-preserving semigroup has at least one stationary state. A practical way to see this is to average a trajectory over a long time:

$$\bar\rho_T = \frac{1}{T} \int_0^T dt\,\Phi_t(\rho_0).$$

The state space is compact, and

$$\mathcal L(\bar\rho_T)=\frac{\Phi_T(\rho_0)-\rho_0}{T}\longrightarrow0.$$

Any accumulation point is therefore stationary. This argument proves existence even when an individual trajectory keeps oscillating.

Trace preservation is a statement about the adjoint generator:

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

Thus the identity is always a left zero mode. Right zero modes are stationary operators. Left zero modes beyond the identity are conserved observables or sector labels, and they often explain why the steady state is not unique.

Relaxation Modes

If $\mathcal L$ is diagonalizable, one can expand operators in right eigenoperators:

$$\mathcal L(R_a)=\lambda_a R_a.$$

Then

$$e^{t\mathcal L}R_a = e^{\lambda_a t}R_a.$$

For every finite-dimensional CPTP semigroup, nonstationary modes have

$$\operatorname{Re}\lambda_a\le0.$$

The real part gives decay or growth. Physical stability rules out positive real parts. The imaginary part gives oscillation frequency:

$$e^{\lambda_a t} = e^{\operatorname{Re}\lambda_a t} e^{i\,\operatorname{Im}\lambda_a t}.$$

Zero eigenvalues describe stationary operators. Nonzero purely imaginary eigenvalues describe persistent oscillations, not slow decay. Modes with small negative real parts relax slowly. Jordan blocks for decaying eigenvalues can give $t^m e^{\lambda t}$ factors; blocks on the imaginary axis must be trivial because the semigroup is bounded.

Many Liouvillians are nonnormal, so right eigenoperators need not be orthogonal and may be ill conditioned. Norm-dependent transient amplification and large spectral prefactors are possible, but CPTP evolution still contracts the trace distance between physical states, including distance to a stationary state. Compare spectral predictions with the observable and norm actually being used.

Liouvillian Gap

When the steady state is unique and all other modes decay, the Liouvillian gap is commonly defined as

$$\Delta_{\mathcal L} = \min_{\lambda_a\ne0} \left[ - \operatorname{Re}\lambda_a \right],$$

where the minimum is taken over nonstationary modes. The associated relaxation time scale is roughly

$$\tau_{\mathrm{rel}} \sim \frac{1}{\Delta_{\mathcal L}}.$$

This statement is only a first estimate. Prefactors, nonnormality, initial-state overlap with slow modes, and time-dependent observables can matter. Still, the gap is the main first diagnostic for how quickly a Markovian system forgets its initial condition.

Multiple zero modes must be counted separately: the displayed minimum excludes them and therefore need not become zero when stationary states are nonunique. Pure dephasing, for example, has eigenvalues $0,0,-\Gamma_\phi,-\Gamma_\phi$; its decaying modes have gap $\Gamma_\phi$, although its stationary manifold is a line. Nonzero purely imaginary modes do make that minimum zero. Always state whether “gap” means decay toward a stationary manifold or toward one unique state.

Unique Versus Multiple Steady States

A unique steady state means there is exactly one density operator $\rho_{\mathrm{ss}}$ satisfying $\mathcal L(\rho_{\mathrm{ss}})=0$. For a finite-dimensional continuous-time GKSL semigroup, uniqueness implies convergence from every initial density operator (Schirmer and Wang, 2010). It does not imply rapid convergence. This implication must not be transferred without qualification to arbitrary discrete channels or infinite-dimensional systems.

Multiple steady states can arise from:

• conserved quantities;
• disconnected classical sectors;
• decoherence-free subspaces;
• symmetry-protected dark subspaces;
• absorbing states in a rate equation;
• missing weak processes that would connect sectors in a more complete model.

For a Pauli Rate Equation, this is graph theory: disconnected communicating classes or absorbing states give multiple stationary distributions. For a quantum Liouvillian, the same idea appears as a nontrivial zero-mode structure of $\mathcal L$ and $\mathcal L^\dagger$.

In experiments and simulations, an apparent family of steady states can also be a metastable manifold: the system looks stationary for a long time, but an extremely slow mode eventually selects a unique state.

Thermal and Nonequilibrium Steady States

A thermal steady state has the Gibbs form

$$\rho_\beta = \frac{e^{-\beta H_S}} {\operatorname{Tr}(e^{-\beta H_S})}$$

for the standard weak-coupling secular generator with thermal rates. Its commuting Lamb shift does not automatically change the bare Gibbs weights. Strong-coupling equilibrium can instead involve a Hamiltonian of mean force; use a consistent derivation. Detailed balance explains the thermal rate condition.

Stationarity alone is weaker:

$$\mathcal L(\rho_\beta)=0$$

does not by itself prove detailed balance. A generator can have a stationary distribution while supporting irreversible cycles, currents, or coherences. Detailed balance rules out equilibrium currents by relating forward and reverse processes.

With multiple baths, drives, measurement feedback, biased reservoirs, or engineered dissipation, the steady state is generally a nonequilibrium steady state. It can be time independent while carrying heat, particle, photon, or information currents.

Dark States and Attractors

A pure state $\lvert\psi\rangle$ is a simple dark steady state if

$$L_k\lvert\psi\rangle=0 \quad \text{for all }k, \qquad H\lvert\psi\rangle=E\lvert\psi\rangle.$$

Then

$$\rho_\psi = \lvert\psi\rangle\langle\psi\rvert$$

satisfies $\mathcal L(\rho_\psi)=0$.

This condition is sufficient, not necessary. Lindblad operators are not unique, and pure stationary states can sometimes be written in forms where the jumps do not literally annihilate the vector. The invariant object is the full Liouvillian.

A dark state is also not necessarily an attractor. If several states are dark, the generator may protect a subspace rather than prepare one unique state. Reservoir Engineering is largely the art of making the desired state or subspace stationary while ensuring all unwanted components are drained toward it; Dissipative state preparation explains how to test attraction as well as stationarity.

Metastability

Metastability means the system first approaches a long-lived low-dimensional manifold and only later relaxes to the final steady state. Spectrally, it appears when a few eigenvalues lie very close to zero while the rest are separated by a larger decay scale.

Schematically, suppose

$$0 = \lambda_0, \qquad |\operatorname{Re}\lambda_1|, \ldots, |\operatorname{Re}\lambda_m| \ll |\operatorname{Re}\lambda_{m+1}|.$$

Then after the fast modes die out, the state is well approximated by a combination of the stationary state and the $m$ slow modes. On intermediate times it can look as if there are several steady states, even if the exact steady state is unique.

A concrete qubit example is

$$\mathcal L=\gamma\mathcal D[Z]+\epsilon\mathcal D[X],\qquad 0<\epsilon\ll\gamma.$$

The Bloch components obey $\dot r_x=-2\gamma r_x$, $\dot r_y=-2(\gamma+\epsilon)r_y$, and $\dot r_z=-2\epsilon r_z$. After time of order $1/(2\gamma)$, the state is nearly diagonal; it then drifts much more slowly toward the unique state $I/2$. At $\epsilon=0$ the diagonal states are exactly stationary. For $\gamma=1000\,\mathrm{s}^{-1}$ and $\epsilon=0.01\,\mathrm{s}^{-1}$, the two characteristic times are 0.5 ms and 50 s. This separates an approximate plateau from an exact family of fixed points.

For $\gamma=1000\,\mathrm{s}^{-1}$ and $\epsilon=0.01\,\mathrm{s}^{-1}$, coherence and population imbalance have decay times 0.5 ms and 50 s. The shaded interval illustrates a nearly diagonal intermediate state whose populations have barely changed. The horizontal axis is logarithmic; neither plateau nor the eventual $I/2$ state requires changing the generator.

Examples

Pure dephasing

For qubit pure dephasing with $\Gamma_\phi>0$ in the $\sigma_z$ basis,

$$\dot\rho_{00}=0, \qquad \dot\rho_{11}=0, \qquad \dot\rho_{01}=-\Gamma_\phi\rho_{01}$$

in the rotating frame. Every diagonal state is stationary. The steady-state manifold is not a single point:

$$\rho_{\mathrm{ss}} = \begin{pmatrix} p&0\\ 0&1-p \end{pmatrix}, \qquad 0\le p\le1.$$

Pure dephasing destroys coherence but does not choose a thermal population. See Pure Dephasing Master Equation.

Amplitude damping

For zero-temperature amplitude damping,

$$\dot\rho = \Gamma\mathcal D[\sigma_-]\rho,$$

with $\Gamma>0$, the unique steady state is

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

The population mode decays with rate $\Gamma$, while coherences decay with rate $\Gamma/2$ before adding Hamiltonian phase rotation. Thus different observables can relax on different time scales. See Amplitude Damping Master Equation.

Two-state rate equation

For upward and downward rates,

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

for a positive total rate the stationary population is

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

and the nonzero eigenvalue is

$$\lambda = - (\Gamma_\uparrow+\Gamma_\downarrow).$$

This is the population-sector relaxation mode of the corresponding two-level thermal master equation.

Driven two-level system

The Optical Bloch Equations have a unique steady state under ordinary radiative decay and continuous drive. The nonzero modes can be real or complex depending on drive strength and detuning. Complex modes describe damped Rabi oscillations, while their real parts determine decay envelopes.

Computing the stationary set

Vectorize in a declared convention and solve $\mathbb L\,\mathrm{vec}(\rho)=0$ together with Hermiticity, trace one, and positivity. Check the residual on the original matrix equation. If the kernel is multidimensional, a normalized arbitrary null vector need not even be positive; determine the physical convex set and use conserved observables to select the limit compatible with the initial state.

A small eigenvalue alone does not prove a physically observable plateau. Check mode amplitudes, separation from faster decay rates, and direct evolution. In an oscillator truncation, occupation accumulating at the highest retained level can manufacture a false steady state. Increase the cutoff and compare physical moments and boundary population before interpreting the result.

Exercises

1. Fixed point from the generator. Show that if $\mathcal L(\rho_{\mathrm{ss}})=0$, then $e^{t\mathcal L}\rho_{\mathrm{ss}}=\rho_{\mathrm{ss}}$.

Solution

Use the power series:

$$e^{t\mathcal L}\rho_{\mathrm{ss}} = \sum_{n=0}^{\infty} \frac{t^n}{n!} \mathcal L^n(\rho_{\mathrm{ss}}).$$

The $n=0$ term is $\rho_{\mathrm{ss}}$. Every term with $n\ge1$ contains at least one $\mathcal L(\rho_{\mathrm{ss}})$ and therefore vanishes. Thus

$$e^{t\mathcal L}\rho_{\mathrm{ss}} = \rho_{\mathrm{ss}}.$$

2. Pure-dephasing steady states. For a qubit with $\dot\rho_{01}=-\Gamma_\phi\rho_{01}$ and fixed populations, describe the steady-state set.

Solution

At stationarity, $\rho_{01}=0$ if $\Gamma_\phi>0$, while $\rho_{00}$ and $\rho_{11}$ are unchanged except for normalization. Therefore

$$\rho_{\mathrm{ss}} = \begin{pmatrix} p&0\\ 0&1-p \end{pmatrix}, \qquad 0\le p\le1.$$

There is a continuum of steady states, not a unique thermal state.

3. Two-state relaxation gap. The rate matrix for $p=(p_g,p_e)^{\mathsf T}$ is

$$K = \begin{pmatrix} -\Gamma_\uparrow & \Gamma_\downarrow\\ \Gamma_\uparrow & -\Gamma_\downarrow \end{pmatrix}.$$

Find its eigenvalues.

Solution

Because columns sum to zero, one eigenvalue is $0$. The trace is

$$\operatorname{Tr}K = - (\Gamma_\uparrow+\Gamma_\downarrow).$$

For a $2\times2$ matrix, the sum of eigenvalues is the trace, so the other eigenvalue is

$$\lambda_1 = - (\Gamma_\uparrow+\Gamma_\downarrow).$$

The relaxation gap is $\Gamma_\uparrow+\Gamma_\downarrow$ when their sum is positive, including one-way relaxation. If both vanish, every population vector is stationary.

4. Dark but not unique. Two qubits undergo collective decay with $J_-=\sigma_-^{(1)}+\sigma_-^{(2)}$. Both $\lvert gg\rangle$ and the singlet $\lvert\psi_-\rangle=(\lvert ge\rangle-\lvert eg\rangle)/\sqrt2$ are annihilated by $J_-$. Does collective decay alone uniquely prepare the singlet?

Solution

No. Both states are dark under the same jump operator:

$$J_-\lvert gg\rangle=0, \qquad J_-\lvert\psi_-\rangle=0.$$

Since the dark space is at least two-dimensional, collective decay alone does not select the singlet as a unique attractor. Additional drives, dissipators, measurement conditioning, or symmetry breaking are needed to make the target uniquely attracting.

5. Thermal versus stationary. A rate equation on three states has a stationary distribution but also a nonzero clockwise probability current around the cycle. Is that an equilibrium detailed-balance steady state?

Solution

No. Stationarity only requires the net inflow and outflow at each state to balance. Detailed balance requires each edge current to vanish pairwise:

$$W_{n\leftarrow m}p_m^{\mathrm{ss}} = W_{m\leftarrow n}p_n^{\mathrm{ss}}.$$

A circulating current violates this pairwise condition, so the state is a nonequilibrium steady state rather than an equilibrium detailed-balance state.

References

• V. V. Albert and L. Jiang, “Symmetries and conserved quantities in Lindblad master equations,” Physical Review A 89, 022118 (2014). DOI.
• R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications, revised and expanded edition, Springer (2007). Publisher.
• H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press (2002). Publisher.
• 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.
• G. Lindblad, “On the generators of quantum dynamical semigroups,” Communications in Mathematical Physics 48, 119–130 (1976). DOI.
• K. Macieszczak, M. Guță, I. Lesanovsky, and J. P. Garrahan, “Towards a theory of metastability in open quantum dynamics,” Physical Review Letters 116, 240404 (2016). DOI.
• S. G. Schirmer and X. Wang, “Stabilizing open quantum systems by Markovian reservoir engineering,” Physical Review A 81, 062306 (2010). Author manuscript.