1 Introduction↩︎

No system we make measurements on is ever exactly isolated. To the best of our knowledge, the quantum density matrix \(\varrho_{\mathrm{world}}\) of everything evolves according to a von Neumann equation of the form \[\label{world} \frac{\partial \varrho_{\mathrm{world}}}{\partial t} = - i [ H_{\mathrm{world}}, \varrho_{\mathrm{world}} ]\tag{1}\] for some Hermitian Hamiltonian \(H_{\mathrm{world}}\). In any given measurement, one however only has access to a subset of all available degrees of freedom and so in practice one tracks the quantum state of an open system defined by \[\rho := \underset{\mathrm{env}}{\mathrm{Tr}}\left[ \varrho_{\mathrm{world}} \right]\] where we have taken a partial trace over the unobserved degrees of freedom living in the environment, giving rise to a reduced density matrix \(\rho\) of size \(D \times D\), and obeying \[\label{DMcond} \rho^{\dagger}= \rho \;, \qquad \mathrm{Tr}[\rho]=1 \;, \qquad \rho \geq 0 \;.\tag{2}\] Unlike Eq. (1 ), the time evolution of \(\rho\) can be non-unitary. In practice, when the map taking \(\rho(0) \to \rho(t)\) is a quantum channel¹, one finds the generic form [3] \[\begin{align} \label{nonMark} \frac{\partial \rho}{\partial t} = - i [ H(t) , \rho(t) ] + \sum_{\alpha\beta} c_{\alpha\beta}(t) \Big[ A_{\alpha}(t) \rho(t) A_{\beta}^{\dagger}(t) - \frac{1}{2} \big\{ A_{\beta}^{\dagger}(t) A_{\alpha}(t), \rho(t) \big\} \Big] \end{align}\tag{3}\] for a set of time-dependent jump operators \(A_{\alpha}(t)\) which are not necessarily Hermitian, a Hermitian system Hamiltonian \(H(t)\) and a set of time-dependent real coefficients \(\mathfrak{c}(t) = [c_{ij}(t)]\). When the time evolution is a quantum channel, this equation is guaranteed to preserve the conditions (2 ) defining a density matrix, although how these conditions manifest themselves in terms of the structure of \(A_{\alpha}(t)\), \(H(t)\) and \(\mathfrak{c}(t)\) can be highly non-trivial in the general (non-Markovian) case [4], [5].

Things simplify significantly in the Markovian limit, where the system evolves slowly compared to the timescales of the environment, allowing the dynamics to be treated as memoryless². In this regime, the general evolution given in Eq. (3 ) simplifies: the jump operators \(A_{\alpha}\) become time-independent constants, and the coefficient matrix \(\mathfrak{c} = [c_{\alpha\beta}]\) becomes constant and positive semi-definite. Diagonalizing the coefficient matrix such that \(\mathfrak{c} = U^{\dagger} \gamma U\), one may rotate the jump operators to \(B_\alpha := \sqrt{\gamma_\alpha} \sum_{\beta} U_{\beta \alpha} A_{\beta}\) (with \(\gamma_\alpha \geq 0\) the positive semi-definite eigenvalues of \(\mathfrak{c}\)), resulting in the standard form of the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) [6], [7] equation, in which the infinitesimal time evolution of \(\rho\) is generated by the Lindbladian \(\mathcal{L}\): \[\begin{align} \label{Lindblad95sch} \frac{\partial \rho}{\partial t} = \mathcal{L}(\rho)\equiv - i [ H , \rho(t) ] + \sum_{\alpha} \Big[ B_{\alpha} \rho B_{\alpha}^{\dagger} - \frac{1}{2} \big\{ B_{\alpha}^{\dagger} B_{\alpha}, \rho(t) \big\} \Big], \end{align}\tag{4}\] where there are at most \(D^2 - 1\) independent jump operators \(B_{\alpha}\) in the second dissipative term of the Lindblad equation responsible for information loss into the environment. The Lindblad equation turns out to be invariant under unitary reshuffling of the jump operators \(B_{\alpha} \to \sum_{\beta} U_{\alpha\beta} B_{\beta}\) as well as the inhomogeneous shift \(B_{\alpha} \to B_{\alpha} + z_\alpha \mathbb{I}\) with \(H \to H - \frac{i}{2} \sum_{\beta} (z^{\ast}_\beta B_{\beta} - z_\beta B_{\beta}^{\dagger} ) + x \mathbb{I}\) for complex \(z_\alpha\) and real \(x\). As done in the seminal work [7], one may use this ambiguity in the choice of \(H\) and \(B_{\alpha}\) to pick the jump operators to be orthonormal with \(\mathrm{Tr}[ B^{\dagger}_{\alpha} B_{\beta} ] = 2 \delta_{\alpha\beta}\) as well as traceless with \(\mathrm{Tr}[B_\alpha] =0\). It turns out that choosing traceless jump operators minimizes the size of the dissipative term under several natural super-operator norms, which provides a physical motivation for this choice and selects a canonical form of the Lindblad equation [8]. This minimality property also holds more generally in certain non-Markovian extensions.

Lindblad equations are highly valuable for modelling a wide range of dissipative processes, particularly in fields such as quantum optics, where they describe photon loss, spontaneous emission, and other cavity QED effects. Beyond quantum optics, they have been applied to condensed matter physics, quantum information theory, quantum error correction, holography and cosmology [9].

One of the features of Lindblad equations is that they generically take pure states to mixed states (inducing decoherence), via the non-unitary dissipative term in Eq. (4 ). To emphasize this point, one may pass to the Dirac (or interaction) picture, so that the unitary terms involving the Hamiltonian are absent. One defines \[\begin{align} \overline{\rho}(t) := \mathcal{U}^{\dagger}(t) \rho(t) \mathcal{U}(t) \qquad \qquad \mathrm{and} \qquad \qquad \overline{B}_{\alpha} = \mathcal{U}^{\dagger}(t) B_\alpha \mathcal{U}(t) \end{align}\] as well as the unitary time evolution operator \(\mathcal{U}(t) = \mathcal{T}\exp\left[ - i \int_0^{t} {\rm d}t'\; H(t') \right]\) using the system Hamiltonian \(H\) (assuming the evolution starts at time \(t=0\)), which converts the Lindblad equation to \[\begin{align} \label{Lindblad} \frac{\partial \overline{\rho}}{\partial t} = \mathcal{L}(\overline{\rho})\equiv \sum_{\alpha} \Big[ \overline{B}_{\alpha}(t) \overline{\rho}(t) \overline{B}_{\alpha}^{\dagger}(t) - \frac{1}{2} \big\{ \overline{B}_{\alpha}^{\dagger}(t) \overline{B}_{\alpha}(t), \overline{\rho}(t) \big\} \Big] \;. \end{align}\tag{5}\] The goal of this work is to develop a simple framework for studying the flow of the above Lindblad equation, centered around the concept of gradient flow³.

1.1 Relation to prior literature↩︎

In the classical setting, Markov processes are well-known (especially in the context of optimal transport) to admit a gradient flow formulation with respect to classical entropy functions. This perspective has been extended to the quantum setting in recent years, where Lindblad equations — describing quantum Markovian dynamics — have likewise been shown to arise as gradient flows of suitable quantum generalizations of entropy [12]–[15]. Note however that these works rely on a gradient fixed by a non-trivial distance metric, and a scalar entropy function often depending on a known steady-state density matrix \(\rho_{\mathrm{eq}}\) and/or dynamics assuming detailed balance. In particular, when the Lindbladian \(\mathcal{L}\) (where \(\mathcal{L}(\rho)\) is the RHS of Eq. (4 )) admits a stationary state \(\rho_{\mathrm{eq}}\) (i.e.\(\mathcal{L}(\rho_{\mathrm{eq}}) = 0\)) and satisfies detailed balance with respect to it, one can define a gradient flow structure for the evolution. The detailed balance condition requires that \(\mathcal{L}\) is self-adjoint with respect to the KMS inner product \(\langle A,B \rangle_{\mathrm{eq}} := \mathrm{Tr}[A^\dagger B \rho_{\mathrm{eq}}]\) so that, for all operators \(A\) and \(B\), \[\langle \mathcal{L}[A], B \rangle_{\mathrm{eq}} = \langle A, \mathcal{L}[B] \rangle_{\mathrm{eq}},\] which ensures the dynamics respects a form of reversibility. Gibbs states, with \(\rho_{\mathrm{eq}} \propto e^{-\beta H}\), are a standard example of such equilibrium states. In the simplest cases, detailed balance implies that rates between energy eigenstates obey \(\Gamma_{n \to m} / \Gamma_{m \to n} = e^{\beta (E_n - E_m)}\) so there is a notion of reversibility.

The works [12], [13], [15] formalize Lindblad evolution with detailed balance as a gradient flow of the von Neumann (in the case of self-adjoint \(\mathcal{L}\)) and quantum relative entropy (relative to the steady state \(\rho_{\mathrm{eq}}\)). In addition, [14] generalizes these constructions to sandwiched Rényi divergences assuming GNS-detailed balance.

In contrast to this body of work, our approach exhibits Lindblad dynamics as a gradient flow with respect to a flat (Euclidean) distance metric, without relying on detailed balance or the previous knowledge of a equilibrium state \(\rho_{\mathrm{eq}}\). This simpler structure allows us to frame the dynamics as ordinary gradient descent in operator space, uncoupled from thermodynamic notions of reversibility.

1.2 Main results↩︎

The main result of the paper is to show that the Lindblad equation (5 ) can usefully be rewritten as \[\label{Lindblad95PhiR95intro} \frac{\partial \rho}{\partial t} = \mathcal{L}(\rho)=- \frac{\partial}{\partial \rho^{T}} \Phi(\rho) + R(\rho)\tag{6}\] with the definitions \[\begin{align} \Phi(\rho) & := & \frac{1}{4} \sum_{\alpha} \bigg( \big| [B_\alpha, \rho ] \big|^2_{\mathrm{HS}} - \mathrm{Tr}\big( \rho^2 [B_{\alpha},B_{\alpha}^{\dagger}]\big) \bigg) \label{scalar95intro} \\ R(\rho) & := & \frac{1}{2} \sum_{\alpha} \Big( B_{\alpha} \rho B_{\alpha}^{\dagger} - B_{\alpha}^{\dagger} \rho B_{\alpha} \Big), \end{align}\tag{7}\] where \(|O|_{\mathrm{HS}} = \sqrt{ \mathrm{Tr|}[ O^{\dagger} O ]}\) is the Hilbert-Schmidt norm, and where \(\frac{\partial}{\partial \rho^{T}} \Phi\) is the matrix gradient of the scalar potential \(\Phi\).

In §2 we consider the case of Hermitian jump operators with \(B_{\alpha}^{\dagger} = B_{\alpha}\), where the Lindblad equation simplifies to pure gradient flow, as the operator \(R(\rho)\) vanishes in this special (degenerate) case, and uniquely determines \(\Phi\) (up to constant shifts). This makes the structure of the evolution particularly transparent: when there are enough orthogonal jump operators, the system tends to flow towards the maximally mixed state, which acts as a late-time attractor. This gradient structure also clarifies why the purity of the state always decreases over time in a particularly simple way.

In the more general case where the jump operators are non-Hermitian, \(R(\rho)\) does not vanish in general (though it does for qubits), resulting in a more intricate form of evolution. We examine this in §3, where we show that the right-hand side of Eq. (6 ) naturally admits a Helmholtz decomposition familiar from vector calculus. In this setting, \(R(\rho)\) is a solenoidal matrix field (divergence-free in \(\rho\)) and crucially, the decomposition into \(R\) and \(\Phi\) is not unique. Nevertheless, there are two ways to interpret the evolution as a kind of gradient flow. The first, discussed in §3.2, involves a complexification of the density matrix. The second, presented in §3.3, exploits the ambiguity in \(R\) and \(\Phi\) by choosing a so-called orthogonal Helmholtz-Hodge decomposition such that \(\frac{\partial}{\partial \rho^{T}} \Phi\) is orthogonal to \(R\). This condition guarantees that \(\Phi\) serves as a Lyapunov function, with the late-time steady states of \(\mathcal{L}\) determined solely by the potential \(\Phi\).

To make the dynamics more tractable, we also additionally express the Lindblad equation in terms of generalized Bloch vectors for both Hermitian and non-Hermitian jump operators, rewriting the \(D \times D\) density matrix \(\rho(t)\) in terms of a real \(D^2 - 1\) dimensional Bloch vector \(\mathbf{a}(t)\). This reformulation is particularly helpful for geometric interpretations and for finding solutions (especially for qudits with \(D \geq 3\)). To illustrate the physical content of Eq. (6 ), we explore explicit examples in the Hermitian case in §2 and in the non-Hermitian case in §3.

In what follows we work exclusively in the Dirac (interaction) picture, and so drop the overlines on the operators when the context is clear.

2 Hermitian jump operators↩︎

To begin, we assume that the jump operators appearing in Eq. (5 ) are Hermitian⁴, so that \[\begin{align} B_\alpha^{\dagger} = B_\alpha \;\forall \alpha. \end{align}\] The main reason we start out with this case is for simplicity. Nevertheless there are many physical situations in which Hermitian jump operators appear in Lindblad equations. For example, in pure dephasing channels the jump operator is Hermitian (a particle number operator), leading to the decay of coherences in the density matrix. Also in the Caldeira-Leggett model of quantum Brownian motion [21], the jump operator takes the schematic form \(B = \sqrt{T} Q + \frac{i}{\sqrt{T}} P\), where \(Q\) and \(P\) are the system position and momentum operators, with \(T\) the temperature. While this operator is generally non-Hermitian, in the high-temperature \(T\to \infty\) limit one finds \(B \simeq \sqrt{T} Q\) (which is Hermitian) showing how Hermitian jump can emerge approximately in certain regimes. In this limit, the dynamics is purely dissipative, consistent with the physical picture that the system primarily loses coherence and energy into a hot bath. More generally, whether Hermitian jump operators appear depends on the specific system-environment model and the approximations made in deriving the Lindblad equation.

2.1 Gradient flow in matrix space↩︎

In this case, the Lindblad equation (5 ) in the Dirac picture simplifies to \[\begin{align} \label{Lindblad95H} \frac{\partial \rho}{\partial t} = - \frac{1}{2} \sum_{i} \big[ B_\alpha, [ B_\alpha, \rho ] \big] \;. \end{align}\tag{8}\] Defining the scalar function as in Eq. (7 ) with \(B_{\alpha}^{\dagger} = B_{\alpha}\), \[\begin{align} \label{Vdef95comm} \Phi(\rho) = \frac{1}{4} \sum_{\alpha} \big| [B_\alpha, \rho ] \big|^2_{\mathrm{HS}}, \end{align}\tag{9}\] one finds generically that \(\Phi \geq 0\). The key observation⁵ is that the gradient of this scalar function with respect to the reduced density matrix is given by \[\frac{\partial}{\partial \rho^{T} } \Phi(\rho) = \frac{1}{2} \sum_{\alpha} \big[ B_{\alpha}(t), [ B_{\alpha}(t), \rho(t) ] \big]\] which means that it simply expresses the Lindblad equation Eq. (8 ) as a gradient flow equation \[\label{Herm95gradientflow} \frac{\partial \rho}{\partial t} = \mathcal{L}(\rho) = - \frac{\partial}{\partial \rho^{T}} \Phi(\rho) \;.\tag{10}\] This formulation is remarkable as it casts the non-unitary evolution of the reduced density matrix in a highly specialized form. One can interpret \(\Phi\) as a potential which the time evolution minimizes — indeed the flow of the state is always in the direction of steepest descent of the potential \(\Phi\). This is clearest when studying the purity of the reduced density matrix, defined as \[\label{purity} \gamma := \mathrm{Tr}[\rho^2] \;.\tag{11}\] which is a simple basis-independent scalar which takes the value \(1/D \leq \gamma \leq 1\) depending on how mixed \(\rho\) is. In the case that \(\rho\) is pure (in which case \(\rho^2 = \rho\) is a projection) it takes the value 1, and in the maximally mixed case (where \(\rho = \mathbb{I}/D\)) takes the value \(1/D\).

By differentiating the purity, one finds that \[\label{diff95purity} \frac{\partial \gamma}{\partial t} = 2 \mathrm{Tr}\Big[ \rho \frac{\partial \rho}{\partial t} \Big] = - 4 \Phi(\rho) \;.\tag{12}\] which shows in a particularly simple way that purity is always non-increasing over the entire evolution (since \(\Phi \geq 0\)). The potential determines the instantaneous rate of change of the purity, or equivalently the rate at which the Rényi-2 entropy \(S_{2} =- \log(\gamma)\) increases monotonically (similar statements apply to the von Neumann entropy, although the relation is more involved). In the Hermitian case, the potential is therefore the simple quantity that drives entropy production.

2.1.1 Qubit flows↩︎

As an illustration of the gradient flow picture, let us consider a simplified scenario in which \(D=2\) and the jump operators are constant in time. We study the reduced density of a qubit with evolving under the influence of 3 constant Hermitian jump operators \(\{ B_{\alpha} \}_{\alpha=1}^3\) (assumed to be distinct).

The gradient flow in Eq. (10 ) suggests that the evolution tends toward minima of the potential. Minimizing the potential means to find the zeros of the right hand side of Eq. (10 ), and so one seeks the steady states \(\rho_{\mathrm{SS}}\) where \[\label{steady95offlow} \frac{\partial}{\partial \rho^{T} } \Phi\; \bigg|_{\rho = \rho_{\mathrm{SS}}} = \frac{1}{2} \sum_{\alpha} \big[ B_\alpha, [ B_{\alpha}, \rho_{\mathrm{SS}} ] \big] = 0.\tag{13}\] In the generic case of three linearly independent jump operators, the steady state must therefore be proportional to the identity matrix, and is thus given by the maximally mixed state \(\rho_{\mathrm{SS}} = \frac{1}{2}\mathbb{I}\).

To better understand the dynamics of the solution, we rewrite⁶ \[- \frac{1}{2} \sum_{\alpha} \big[ B_\alpha, [ B_\alpha, B_{\beta} ] \big] = \;- \sum_{\alpha} \bigg( \mathrm{Tr}[B_{\alpha}^2] B_{\beta} - \mathrm{Tr} [B_{\alpha}B_{\beta}] B_{\alpha} \bigg) \;.\] If one furthermore picks the jump operators (see the text below Eq. (4 )) so they are orthogonal under the Hilbert-Schmidt inner product with \(\mathrm{Tr} [B_{\alpha}B_{\beta}]= 2 \delta_{\alpha\beta}\) then one gets \[\label{eigenmatrices} - \frac{1}{2} \sum_{\alpha} \big[ B_\alpha, [ B_\alpha, B_{\beta} ] \big] = \;- \Gamma_{\beta} B_{\beta} \qquad \mathrm{with} \qquad \Gamma_{\beta} = \sum_{\alpha \neq \beta} \mathrm{Tr}[B_{\alpha}^2]\tag{14}\] which shows that \(B_{\beta}\) are eigenmatrices of the operator \(- \frac{\partial}{\partial \rho^{T} } \Phi(\rho) = - \frac{1}{2} \sum_{\alpha} \big[ B_\alpha, [ B_{\alpha}, \rho ] \big]\) appearing on the RHS of the Lindblad equation with eigenvalues \(-\Gamma_{\beta}\). Consequently, the manifestly positive constants \(\Gamma_{\beta}\) are the rates with which the solution relaxes to the maximally mixed state.

2.2 Bloch Vectors↩︎

Quantum mechanics is inherently invariant under changes of bases, and the Lindblad equation of course respects this principle. The Lindblad equation remains unchanged under any unitary transformation \(\rho \to U^{\dagger} \rho U\), where \(U\) is a fixed element of the unitary group \(U(D)\), provided that the operators \(B_i\) and \(H\) are also rotated accordingly. This suggests that \(\rho\) naturally resides within the adjoint representation of \(U(D)\).

Given this, it is natural to express \(\rho\) in terms of the basis elements of the corresponding Lie algebra. Since \(U(D) = U(1) \times SU(D)\) with \(SU(D)\) the special unitary group, the Lie algebra \(\mathfrak{u}(D)\) is insensitive to the \(U(1)\) phase. This means that \(\rho\) can be expressed in terms of the Lie algebra \(\mathfrak{su}(D)\), which consists of traceless, Hermitian \(D \times D\) matrices.

The vector space \(\mathfrak{su}(D)\) is \(D^{2} - 1\) dimensional, which is consistent with the (real) degrees of freedom contained within a density matrix once the conditions of Hermicity and unit trace are accounted for. Taking \(\require{upgreek} \{ \uplambda_{i} \}_{i=1}^{D^2 - 1}\) to be an orthonormal basis for \(\mathfrak{su}(D)\) such that \[\require{upgreek} \label{lambdaconditions} \uplambda_j^{\dagger} = \uplambda_j \;, \qquad \mathrm{Tr}\big[ \uplambda_j \big] = 0\;, \qquad \mathrm{Tr}\big[ \uplambda_j \uplambda_k \big] = 2 \delta_{jk} \;,\tag{15}\] one can then express the density matrix as \[\require{upgreek} \label{blochN} \rho = \frac{1}{D} \; \mathbb{I} + \sqrt{ \text{\scalebox{.85}{\frac{D-1}{2D}}} } \; a_j \uplambda_{j},\tag{16}\] where we use the Einstein summation convention for repeated indices, and the vector \(\mathbf{a} \in \mathbb{R}^{D^2 -1}\) is the generalized Bloch vector [22]–[28] corresponding to the density matrix \(\rho\). Given the conditions (15 ), the inversion of Eq. (16 ) is given by \[\require{upgreek} a_j = \frac{1}{\sqrt{2 - 2/D}} \; \mathrm{Tr}[ \rho \uplambda_{j} ] \;.\] The Bloch vector is a convenient object to study, as it completely eliminates the redundant degrees of freedom present in the \(D \times D\) complex-valued density matrix. Moreover, it enables a direct analysis of the density matrix within the more intuitive Euclidean space, providing a clearer geometric interpretation of the behaviour of the density matrix.

Important to this geometric interpretation is the expression of the purity (11 ) as \[\begin{align} \label{purityN} \gamma := \mathrm{Tr}[ \rho^2 ] = \frac{1 + (D-1) |\mathbf{a}|^2}{D} = \begin{cases} \; \frac{1}{D} \qquad & \mathrm{if} \; \mathbf{a} = \mathbf{0} \\ \; \; 1 \qquad & \mathrm{if} \; |\mathbf{a}| = 1 \end{cases} \end{align}\tag{17}\] This means that \(\mathrm{Tr}[\rho^2] \leq 1\) implies that Bloch vectors lie within a hyperball of unit radius \(|\mathbf{a}| \leq 1\), with the boundary of the ball corresponding to pure states, and \(\mathbf{a} = \mathbf{0}\) corresponding to the maximally mixed state.

The conditions of Hermiticity and unit trace are built into the Bloch vector definition in Eq. (16 ). The tradeoff is that positivity (\(\rho \geq 0\)) is generally hard to characterize in terms of the Bloch vector except for the special case of qubits. When \(D = 2\), there is a one-to-one correspondence between physical density matrices and the unit ball in \(\mathbb{R}^3\), since in that case positivity is equivalent to \(\gamma \leq 1\).

For \(D \geq 3\), this no longer holds: many vectors \(\mathbf{a}\) defined by Eq. (16 ) correspond to non-positive matrices and therefore the set of physical Bloch vectors is a proper subset of the unit ball in \(\mathbb{R}^{D^2 - 1}\), and it is non-trivial to describe the structure fo this set [22]–[27].

Crucial to the study of generalized Bloch vectors is the identity [29] \[\require{upgreek} \label{mult} \uplambda_j \uplambda_k = \frac{2}{D} \delta_{jk} \mathbb{I}_D + ( i f_{jk\ell} + d_{jk\ell} ) \uplambda_\ell,\tag{18}\] where \(f_{jk\ell}\) are the (completely antisymmetric) structure constants and \(d_{jk\ell}\) are the (completely symmetric) “\(d\)-coefficients” (which turn out to vanish identically for the case \(D=2\) and otherwise are non-zero). This means that \[\require{upgreek} \label{fcomm95def} \left[ \uplambda_j , \uplambda_k \right] = 2 i f_{jk\ell} \uplambda_\ell \qquad \mathrm{and} \qquad \left\{ \uplambda_j , \uplambda_k \right\} = \frac{4}{D} \delta_{jk} \mathbb{I}_D + 2 d_{jk\ell} \uplambda_\ell \;.\tag{19}\] Given a basis \(\require{upgreek} \{ \uplambda_j \}_{j=1}^{N^2 - 1}\) satisfying (15 ) one can determine \(f_{jk\ell}\) and \(d_{jk\ell}\) via the relationships \[\require{upgreek} \begin{align} f_{jk\ell} = - \tfrac{i}{4} \mathrm{Tr}\Big( \left[ \uplambda_j , \uplambda_k \right] \lambda_{\ell} \Big) \qquad \mathrm{and} \qquad d_{jk\ell} = \tfrac{1}{4} \mathrm{Tr}\Big( \left\{ \uplambda_j , \uplambda_k \right\} \uplambda_\ell \Big) \;. \end{align}\] Given these definitions, given two vectors \(\mathbf{a},\mathbf{b} \in \mathbb{R}^{D^2 - 1}\) one may define two vector operations [22], [30] \[\begin{align} (\mathbf{a} \wedge \mathbf{b})_{j} = f_{jk \ell} a_{j} b_{k} \qquad \mathrm{and} \qquad (\mathbf{a} \star \mathbf{b})_{j} = d_{jk \ell} a_{j} b_{k} \;. \end{align}\] Since these operations depend on a choice of \(\require{upgreek} \{\uplambda_{j}\}_{j=1}^{D^2 - 1}\) they are basis-dependent operations. There is a sense in which the antisymmetric \(\wedge\)-product is a generalization of the cross product in 3-dimensions, although there isn’t a clear analog for the symmetric \(\star\)-product. Neither of these operations are associative i.e.\(( \mathbf{a} \star \mathbf{b} ) \star \mathbf{c} \neq \mathbf{a} \star ( \mathbf{b} \star \mathbf{c} )\), although they are distributive such that \(\mathbf{a} \star ( \mathbf{b} + \mathbf{c} ) = \mathbf{a} \star \mathbf{b} + \mathbf{a} \star \mathbf{c}\) and \((\mathbf{a} \star \mathbf{b} ) \cdot \mathbf{c} = ( \mathbf{b} \star \mathbf{c} ) \cdot \mathbf{a}\) and satisfy for example [30] \[\begin{align} (\mathbf{a} \wedge \mathbf{b}) \cdot (\mathbf{c} \wedge \mathbf{d}) =\big( (\mathbf{a} \wedge \mathbf{b} ) \wedge \mathbf{c} \big) \cdot \mathbf{d} \end{align}\] Since these operations are only sparsely discussed in the literature, we summarize several of their properties in Appendix 5.

2.3 Gradient flow in the Bloch sphere↩︎

Our goal now is to express the earlier Lindblad equation on the Bloch sphere to facilitate an easier geometric interpretation [31], [32], and to easier find solutions. Assuming the Hermitian jump operators \(B^{\alpha}\) are traceless, one can parametrize them in terms of “jump vectors” \(\mathbf{b}^{(\alpha)}\) such that \[\require{upgreek} \begin{align} B_{\alpha} = \sqrt{ \text{\scalebox{.85}{\frac{D-1}{2D}}} } \; b^{(\alpha)}_j \uplambda_{j} \;. \end{align}\] Using Eq. (19 ) the Lindblad equation (8 ) is trivially found to be⁷ \[\begin{align} \frac{\partial \mathbf{a}}{\partial t} = \frac{D-1}{D} \; \mathbf{b} \wedge \big( \mathbf{b} \wedge \mathbf{a} \big). \end{align}\] It may be equivalently expressed as a gradient flow equation \[\begin{align} \frac{\partial \mathbf{a}}{\partial t} = - \frac{D}{D-1} \; \nabla_{\mathbf{a}} \Phi(\mathbf{a}) \qquad \mathrm{with} \qquad \Phi(\mathbf{a}) = \frac{(D-1)^2}{2D^2} \sum_{\alpha} | \mathbf{b}^{(\alpha)} \wedge \mathbf{a} |^2 \;. \end{align}\] This complements the earlier gradient flow equation, by mapping it to real Euclidean space. This linear equation can also be written as the matrix equation \[\begin{align} \label{Mmatrix95def} \frac{\partial \mathbf{a}}{\partial t} = \mathbb{M}\mathbf{a}(t) \qquad \mathrm{with} \qquad \mathbb{M}_{ij} = \frac{D-1}{D} \sum_{\alpha} f_{i k \ell } f_{\ell m j} b^{(\alpha)}_{k} b^{(\alpha)}_{m} \;. \end{align}\tag{20}\] Note furthermore that the matrix \(\mathbb{M}\) is closely related the scalar potential \(\Phi\) in the sense that it is directly proportional to its Hessian such that \(\frac{\partial^2 \Phi}{\partial a_i \partial a_j} = - \frac{D-1}{D} \mathbb{M}_{ij}\).

2.3.1 Qubit flow in the Bloch sphere↩︎

Let us return to the \(D=2\) qubit problem from §2.1.1 now phrased in terms of Bloch vectors. Here the basis elements of \(\mathfrak{su}(2)\) are of course the Pauli matrices, where the structure coefficients are the Levi-Cevita symbol \(f_{ijk} = \epsilon_{ijk}\), and the \(d\)-ceofficients vanish identically \(d_{ijk}=0\) (in a sense much of the simplicity that arises for qubits is due to this fact).

One finds now in this language that the matrix \(\mathbb{M}\) from earlier is \[\label{M95qubit1} \mathbb{M} = \frac{1}{2} \sum_{\alpha} \big( \mathbf{b}^{(\alpha)}\mathbf{b}^{(\alpha)T} - (\mathbf{b}^{(\alpha)} \cdot \mathbf{b}^{(\alpha)} ) \mathbb{I} \big) \;.\tag{21}\] Since \(\mathrm{Tr}[B_{\alpha} B_{\beta}] = \frac{1}{2} \mathbf{b}^{(\alpha)} \cdot \mathbf{b}^{(\beta)}\) here (since the \(d\)-coefficients vanish for \(\mathfrak{su}(2)\)), the orthogonality condition from §2.1.1 implies that the Bloch vectors are orthogonal with respect to the standard dot product with \(\mathbf{b}^{(\alpha)} \cdot \mathbf{b}^{(\beta)} \propto \delta_{\alpha \beta}\). In this case, the matrix is easily diagonalized, since one now finds \[\mathbb{M} \mathbf{b}^{(\beta)} = \;- \Gamma_{\beta} \mathbf{b}^{(\beta)} \qquad \mathrm{with \;} \Gamma_{\beta} = \tfrac{1}{2} \sum_{\alpha \neq \beta} \mathbf{b}^{(\alpha)} \cdot \mathbf{b}^{(\alpha)}\] with the same rates \(\Gamma_{\beta}\) as given in Eq. (14 ), now written in terms of Bloch vectors. The solution (for constant jump operators) is then given by \[\mathbf{a}(t) = e^{\mathbb{M} t} \mathbf{a}(0) = \sum_{\beta} e^{- \Gamma_{\beta} t} \frac{\mathbf{b}^{(\beta)} \cdot \mathbf{a}(0) }{ \mathbf{b}^{(\beta)} \cdot \mathbf{b}^{(\beta)}} \; \mathbf{b}^{(\beta)}\] This shows how the solution sinks towards the steady solution at \(\mathbf{0}\) at late times \(t \to \infty\), corresponding to the maximally mixed state (consistent with the text surrounding Eq. (13 )).

In the case that we have a single jump operator the picture slightly changes: consider \(B^{(1)} = B\) and \(B^{(2)} = B^{(3)} = 0\) so that the earlier matrix (21 ) becomes \[\begin{align} \mathbb{M} = \frac{1}{2} \big( \mathbf{b} \mathbf{b}^{T} - (\mathbf{b} \cdot \mathbf{b} ) \mathbb{I} \big) \end{align}\] with Bloch vector \(\mathbf{b}\) corresponding to the traceless \(B\). The spectrum of this matrix is different since \(\mathbb{M} \mathbf{b} = \mathbf{0}\) and \(\mathbb{M} \mathbf{b}_{\perp} = - \frac{1}{2} ( \mathbf{b} \cdot \mathbf{b} ) \mathbf{b}_{\perp}\) for any vector \(\mathbf{b}_{\perp}\) perpendicular to \(\mathbf{b}\), where of course the eigenspace corresponding to eigenvalue \(- \frac{1}{2} ( \mathbf{b} \cdot \mathbf{b} )\) is of dimension 2. The solution in this case is \[\begin{align} \label{singlejump95qubit} \mathbf{a}(t) & = & \frac{ \mathbf{a}(0) \cdot \mathbf{b} }{ |\mathbf{b}|^2 } \mathbf{b} + \frac{e^{- \frac{1}{2} | \mathbf{b} |^2 t} }{|\mathbf{b}|^2} \; \big( \mathbf{b} \times \mathbf{a}(0) \big) \times \mathbf{b} \end{align}\tag{22}\] where \(\wedge = \times\) is the usual cross product in \(\mathbb{R}^3\). We see that the state sinks towards the projection of the initial state vector \(\mathbf{a}(0)\) onto \(\mathbf{b}\).

2.3.2 Qutrit flow on the Bloch sphere↩︎

All of the earlier qubit calculations can be carried out just as easily using either Bloch vectors or matrices. Here, we demonstrate how Bloch vectors can significantly simplify the solution of Lindblad equations in cases where solving them directly in matrix form would be considerably more challenging. Specifically, we consider the case of \(D=3\) qutrits with a single jump operator, so that the Lindblad equation takes the form (c.f. Eq. (20 )). \[\begin{align} \frac{\partial \mathbf{a}}{\partial t} = \mathbb{M}\mathbf{a}(t) \qquad \mathrm{with} \qquad \mathbb{M}_{ij} = - \frac{2}{3} \sum_{\alpha} f_{\ell i k } f_{mjk} b^{(\alpha)}_{m} b^{(\alpha)}_{\ell} \;. \end{align}\] where \(\mathbf{a},\mathbf{b} \in \mathbb{R}^{8}\) now. To actually perform calculations in this case it is useful to go an specific basis of \(\mathfrak{su}(3)\), where the obvious choice are the well-known Gell-mann matrices⁸. In this basis, one has the Cartan-Killing metric tensor \[\begin{align} f_{ \ell j k} f_{ m j k } = 3 \delta_{\ell m}\; , \end{align}\] which can be used to show for example that \(\mathrm{Tr}[\mathbb{M}] = - 2 |\mathbf{b}|^2\). Repeated application of this and other identities found in Appendix 5 shows that the characteristic equation for the eigenvalues \(-\Gamma\) of the matrix \(\mathbb{M}\) is given by \[\begin{align} (-\Gamma)^2 \Big( \; (-\Gamma)^3 + |\mathbf{b}|^2 (-\Gamma)^2 + \tfrac{1}{4} |\mathbf{b}|^4 (-\Gamma) + \tfrac{1}{54} |\mathbf{b}|^6 - \tfrac{1}{18} \big[ (\mathbf{b} \star \mathbf{b}) \cdot \mathbf{b} \big]^2 \; \Big)^2= 0 \;. \end{align}\] The cubic equation in the above always has real roots (which can be seen by studying the cubic discriminant), and one finds that there are 8 eigenvalues \[\begin{align} (-\Gamma) \in \{ \; 0, \; 0 , \; r^{(0)}, \; r^{(0)}, \; r^{(1)}, \; r^{(1)}, \; r^{(2)}, \; r^{(2)} \; \} \end{align}\] with the cubic roots \[\begin{align} r^{(n)} = - \frac{ 2 }{ 3 } |\mathbf{b}|^2 \sin^{2}\bigg( \; \frac{1}{6} \; \mathrm{arccos}\left( \frac{ 6 \; [ \mathbf{b} \star \mathbf{b} \cdot \mathbf{b} ]^2 }{ | \mathbf{b} |^6 } - 1 \right) - \frac{\pi n}{3} \; \bigg) \qquad \mathrm{for} \; n = 0 , 1, 2 \;. \end{align}\] One can see trivially that the eigenvalues \(-\Gamma\) are negative semi-definite, which means that the rates \(\Gamma\) are positive and so the solution sinks towards the steady state. Furthermore, the eigenspace corresponding to the eigenvalue \(0\) is spanned by the vectors \(\mathbf{b}\) and \(\mathbf{b} \star \mathbf{b}\), which means that the steady state is a projection onto this subspace analagous to the qubit case in Eq. (22 ). The degeneracy seen in the above eigenvalues is in fact intimately related to the star product: defining a matrix \(\mathbb{S}_{\mathbf{b}}\) taking \(\mathbb{S}_{\mathbf{b}} \mathbf{a} := \mathbf{b} \star \mathbf{a}\), it turns out that \([\mathbb{M}, \mathbb{S}_{\mathbf{b}}] = 0\). This shows that the single-jump Hermitian operator Lindblad equation is invariant under star products.

3 Non-Hermitian jump operators↩︎

In this section we study the full Lindblad equation with non-Hermitian jump operators \[\begin{align} \label{nonHerm95Lind2} \frac{\partial \rho}{\partial t} = \mathcal{L}(\rho)=\sum_{\alpha} \Big[ B_{\alpha} \rho B_{\alpha}^{\dagger} - \frac{1}{2} \big\{ B_{\alpha}^{\dagger} B_{\alpha}, \rho \big\} \Big] \;. \end{align}\tag{23}\] We see that it does not quite undergo gradient flow in the sense described in §2. Nevertheless we show that it exhibits several features which prove useful and similar to gradient flow: in particular, the steady states are still those states which minimize a scalar potential.

3.1 Lack of gradient flow↩︎

Let us begin by analogy with the earlier Eq. (12 ) which expresses the evolution of the purity in terms of a scalar function \(\Phi\). We find in the non-Hermitian case that \[\frac{\partial \gamma}{\partial t} = 2 \mathrm{Tr}\Big[ \rho \frac{\partial \rho}{\partial t} \Big] = - 4 \Phi(\rho),\] with the definition \[\Phi(\rho) := \frac{1}{2} \sum_{\alpha} \bigg( \mathrm{Tr}\big( \rho^2 B_\alpha^{\dagger} B_{\alpha}\big) - \mathrm{Tr}\big( \rho B_{\alpha}^{\dagger} \rho B_{\alpha} \big) \bigg) = \frac{1}{4} \sum_{\alpha} \bigg( \big| [B_\alpha, \rho ] \big|^2_{\mathrm{HS}} - \mathrm{Tr}\big( \rho^2 [B_{\alpha},B_{\alpha}^{\dagger}]\big) \bigg).\] Notice how \(\Phi(\rho)\) reduces to Eq. (9 ) in the Hermitian case where \(B_{\alpha}^{\dagger} = B_{\alpha}\). The potential \(\Phi(\rho)\) now can have any sign (contrary to the Hermitian case); however, this point will not matter too much, as the potential is in fact no longer unique in this case (for reasons that will become clear below).

Importantly, notice that the gradient of this scalar potential does not yield the RHS of the Lindblad equation since \[- \frac{\partial}{\partial \rho^{T}} \Phi(\rho) = \frac{1}{2} \sum_{\alpha} \bigg( B_{\alpha} \rho B^{\dagger}_{\alpha} + B^{\dagger}_{\alpha} \rho B_{\alpha} - \left\{ B_{\alpha}^{\dagger} B_{\alpha}, \rho \right\} \bigg),\] c.f. Eq. (10 ). Indeed, the gradient of \(\Phi\) generates a string of operators which is necessarily symmetric under the interchange of \(B_{\alpha}\) and \(B_{\alpha}^{\dagger}\). The anti-symmetric nature of the Lindblad equation under this interchange is why it cannot be expressed as gradient flow. We instead rewrite Eq. (23 ) as \[\label{Lindblad95PhiR} \frac{\partial \rho}{\partial t} = \mathcal{L}(\rho)= - \frac{\partial}{\partial \rho^{T}} \Phi(\rho) + R(\rho),\tag{24}\] with the definition \[R(\rho) := \frac{1}{2} \sum_{\alpha} \Big( B_{\alpha} \rho B_{\alpha}^{\dagger} - B_{\alpha}^{\dagger} \rho B_{\alpha} \Big).\] The second term is anti-symmetric under \(B_{\alpha} \to B_{\alpha}^{\dagger}\), which is why it cannot be generated from gradient flow. We can interpret the Lindblad evolution as the combination of a driving force \(R(\rho)\) and a conservative force generated by the potential \(\Phi\) .

3.2 Gradient flow via complexification↩︎

We can generate the full Lindbladian as a gradient by complexifying \(\rho\), i.e.by considering a flow in the larger space of matrices with \(\rho\) not necessarily Hermitian. Consider the scalar potential \[\mathcal{F}(\rho,\rho^{\dagger}) = - \frac{1}{2} \sum_{\alpha} \mathrm{Tr}\big[ B_{\alpha} \rho^\dagger B_{\alpha}^\dagger \rho - B_{\alpha}^\dagger \rho^\dagger B_{\alpha} \rho \big],\] which treats \(\rho\) and \(\rho^\dagger\) as independent variables.

Differentiating \(\mathcal{F}\) with respect to \(\rho\) while holding \(\rho^\dagger\) fixed yields \[\frac{\partial \mathcal{F}}{\partial \rho^T} = - \frac{1}{2} \sum_{\alpha} \bigg( B_{\alpha} \rho^\dagger B_{\alpha}^\dagger - B_{\alpha}^\dagger \rho^\dagger B_{\alpha} \bigg). \;\] If we then proceed to set \(\rho^{\dagger} = \rho\), we arrive at \[- \frac{\partial \mathcal{F}}{\partial \rho^T}\Big|_{\rho=\rho^\dagger} = \frac{1}{2} \sum_{\alpha} \bigg( B_{\alpha} \rho B_{\alpha}^\dagger - B_{\alpha}^\dagger \rho B_{\alpha} \bigg) = R(\rho) \;.\] There is thus a sense in which the RHS of the Lindblad equation is the gradient of \(\Phi + \mathcal{F}\), but on the space of complexified density matrices. Therefore \[\frac{\partial \rho}{\partial t}= - \frac{\partial }{\partial \rho^T}\Big( \Phi(\rho)+\mathcal{F}(\rho,\rho^{\dagger}) \Big)\; \Big|_{\rho=\rho^\dagger} \;,\] We note that this equation must be interpreted with some care, because the condition \(\rho=\rho^\dagger\) is imposed after taking the derivative \(\frac{\partial}{\partial \rho^T}\). Thus this is not a gradient flow in the traditional sense, rather a flow with respect to a gradient in a larger space (the space of complexified density matrices) which is then projected onto the desired subspace (the subspace of positive-semidefinite operators with unit trace) at each step along the flow.

3.3 Orthogonal Helmholtz-Hodge decomposition↩︎

We first observe that the matrix divergence of the operator \(R(\rho)\) in 24 vanishes: \[\mathrm{div}_{\mathbf{\rho}}R(\rho) := \sum_{ij} \frac{\partial R_{ij}}{\partial \rho_{ij}} = 0 \;.\] This means that \(R(\rho)\) is a solenoidal or incompressible matrix field, in the standard sense of multivariable calculus. Hence \(R(\rho)\) has no sources or sinks as a function of \(\rho\). Thus it is natural to expect that the steady states of the Lindblad equation are determined by the scalar potential alone, with \(R\) instead acting analogously to a magnetic field in electromagnetism: it induces rotations but does not by itself create or remove steady states (specifically in the case when \(R\) and \(\nabla \Phi\) are orthogonal as we soon discuss). While \(R(\rho)\) does not contribute simply to purity loss, it may play a role in redistribution processes that could be connected to notions of energy or heat currents [17] in specific system-bath models. A precise identification would require embedding the Lindblad equation into a microscopic setting, which lies beyond our present scope.

The decomposition of the Lindblad equation into a gradient term and a solenoidal term is not unique: consider any scalar function \(\varphi(\rho)\) which satisfies \[\label{LaplacianVarphi} \Delta_{\rho} \varphi(\rho) := \mathrm{div}_{\mathbf{\rho}}\frac{\delta}{\delta \rho^{T}} \varphi(\rho) = \sum_{ij} \frac{\partial^2 \varphi}{\partial^2 \rho_{ij}} = 0\tag{25}\] i.e.where the operation \(\Delta_{\rho}\) is the matrix Laplacian, so that the scalar \(\varphi\) is a harmonic function of \(\rho\). One can then shift \(\Phi\) and the solenoidal field \(R\) such that \[\label{gauge1} \Phi \to \widetilde{\Phi} = \Phi + \varphi \qquad\qquad \mathrm{and} \qquad\qquad R \to \widetilde{R} = R + \text{\scalebox{.85}{\frac{\delta \varphi}{\delta \rho^{T}}}}\tag{26}\] where \(\widetilde{R}\) is also a solenoidal field satisfying \(\mathrm{div}_{\mathbf{\rho}} \widetilde{R}(\rho)=0\) because of Eq. (25 ), and the Lindblad equation remains invariant under this transformation, since \[\label{Lindblad95PhiRTILDE} \frac{\partial \rho}{\partial t} = \mathcal{L}(\rho)= - \frac{\partial}{\partial \rho^{T}} \widetilde{\Phi}(\rho) + \widetilde{R}(\rho) \;.\tag{27}\] cf. Eq. (24 ). Of course, in taking the transformation (26 ) all we’ve done is added and subtracted \(\frac{\delta \varphi}{\delta \rho^{T}}\), so there should be no surprise that the Lindblad equation remains the same. As we describe next, we can exploit this redundancy to obtain a type of generalized gradient flow.

Generalizing gradient flow↩︎

We saw above that the Lindblad equation for non-Hermitian jump operators is, strictly speaking, not gradient flow. When interpreted correctly as an orthogonal Helmholtz-Hodge decomposition [36]–[38], there is a sense in which the Lindblad can nonetheless be interpreted as a generalization of gradient flow.

To begin, we compute the orbital derivative of \(\Phi\) along solutions \(\rho(t)\) of the Lindblad equation as \[\frac{\partial }{\partial t} \Phi\big(\rho(t) \big) = \sum_{jk} \frac{\partial \Phi}{\partial \rho_{jk}} \frac{\partial \rho_{jk}}{\partial t} = \mathrm{Tr}\bigg( \frac{\partial \Phi}{\partial \rho^{T}} \frac{{\rm d}\rho}{{\rm d}t} \bigg) \;.\] Using the Lindblad equation Eq. (24 ) one finds that this is equivalently \[\label{orbPhi} \frac{\partial }{\partial t} \Phi\big(\rho(t) \big) =- \mathrm{Tr}\bigg( \frac{\partial \Phi}{\partial \rho^{T}} \frac{\partial \Phi}{\partial \rho^{T}} \bigg) + \mathrm{Tr}\bigg( \frac{\partial \Phi}{\partial \rho^{T}} R \bigg) \;.\tag{28}\] The strategy now is to eliminate the second term in the above equation. As emphasized earlier, the Helmholz-Hodge decomposition is not unique, and one may always shift \(\Phi\) and \(R\) by a harmonic function \(\varphi\) satisfying \(\Delta_{\rho} \varphi =0\). This means one can select \(\varphi\) so that \(\widetilde{\Phi} = \Phi + \varphi\) and \(\widetilde{R} = R + \frac{\partial \varphi}{\partial \rho^T}\) give rise to \[\label{orthHH} \mathrm{Tr}\Big( \frac{\partial \widetilde{\Phi}}{\partial \rho^{T}} \widetilde{R} \Big) =0,\tag{29}\] i.e. so that the (Hermitian) matrices \(\frac{\partial \widetilde{\Phi}}{\partial \rho^{T}}\) and \(\widetilde{R}\) are orthogonal under the Hilbert-Schmidt inner product. Note that since the Lindblad equation is linear in \(\rho\) it is always true that one can find such a decomposition [38]. Picking such \(\widetilde{\Phi}\) and \(\widetilde{R}\) constitutes an orthogonal Helmholtz-Hodge decomposition, where one now finds that Eq. (28 ) gets translated to \[\frac{\partial }{\partial t} \widetilde{\Phi}\big(\rho(t) \big) =- \bigg| \frac{\partial \widetilde{\Phi}}{\partial \rho^{T}} \bigg|^2_{\mathrm{HS}} \;.\] This is useful because it shows that solutions \(\rho(t)\) ensure that the orbital derivative of \(\Phi'\) is non-increasing with \(\frac{{\rm d}}{{\rm d}t } \widetilde{\Phi} \leq 0\). This means that the behavior of solutions is governed by the level sets of \(\widetilde{\Phi}(\rho(t))\). In particular, we find that the potential is minimized at the equilibrium point, and \(\widetilde{\Phi}\) can be interpreted as a Lyapunov function [39].

Perhaps the most useful property of orthogonal Helmholtz decompositions is that at a fixed point \(\rho_{\mathrm{SS}}\) (defined by \(\mathcal{L}(\rho_{\mathrm{SS}}) = 0\)) one finds that \(0 =- \frac{\partial \widetilde{\Phi}}{\partial \rho^{T}} + \widetilde{R} \; |_{\rho_{\mathrm{SS}}}\). It then trivially follows from the orthogonality property (29 ) that \[\mathcal{L}(\rho_{\mathrm{SS}}) = 0 \qquad \iff \qquad \frac{\partial \widetilde{\Phi}}{\partial \rho^{T}} \Big|_{\rho_{\mathrm{SS}}} = 0\] since the orthogonality property (29 ) means that \(| \frac{\partial \widetilde{\Phi}}{\partial \rho^{T}} |^2_{\mathrm{HS}} \big|_{\rho_{\mathrm{SS}}} = 0\) at the steady state. Hence we have established that the fixed points of the evolution are determined purely by the potential \(\widetilde{\Phi}\).

In practice, however, it appears that it is non-trivial to find potentials satisfying Eq. (29 ), although there are simple ways of proceeding. For example, a simple harmonic function \(\varphi\) that is quadratic in \(\rho\) and can be used to shift \((\Phi,R) \to (\widetilde{\Phi},\widetilde{R})\) via Eq. (26 ) is \(\varphi = \mathrm{Tr}[ \rho^2 X ]\) with the constraint \(\mathrm{Tr}[X] =0\). This is because \(\Delta_{\rho} \mathrm{Tr}[ \rho^2 X ] = 2 D \mathrm{Tr}[X]\) which vanishes when \(X\) is traceless and so makes \(\varphi\) harmonic in the sense of Eq. (25 ). In principle, one can pick \(X\) to so that \((\widetilde{\Phi},\widetilde{R})\) are orthogonal.

It turns out however that this approach fails to capture the generic case, so we instead reformulate the discussion in terms of Bloch vectors, where the orthogonality conditions are most transparent.

3.4 Orthogonality in terms of Bloch vectors↩︎

We now pass to the Bloch vector picture, where the somewhat unwieldly matricial calculus gets recast as familiar vector calculus. This is particularly helpful for constructing the conditions required for the potentials \(\widetilde{\Phi}\) and \(\widetilde{R}\) to be orthogonal in Eq. (29 ). In the case of non-Hermitian jump operators, the only difference from before is that we must parametrize the “jump vectors” \(\mathbf{b}^{(\alpha)}\) as complex vectors in \(\mathbb{C}^{D^2 - 1}\), such that \[\require{upgreek} B_{\alpha}= \text{\scalebox{.85}{\frac{\mathrm{Tr}[B_{\alpha}]}{D}}} \; \mathbb{I}_{D} + \sqrt{ \text{\scalebox{.85}{\frac{D-1}{2D}}} } \; b^{(\alpha)}_j \uplambda_{j} \quad \qquad \mathrm{or} \qquad \quad B_{\alpha}^{\dagger} = \text{\scalebox{.85}{\frac{\mathrm{Tr}[B_{\alpha}]^{\ast}}{D}}} \; \mathbb{I}_{D} + \sqrt{ \text{\scalebox{.85}{\frac{D-1}{2D}}} } \; b^{(\alpha)\ast}_j \uplambda_{j} \;.\] By rewriting the Lindblad equation Eq. (23 ) as a string of commutators and anti-commutators, \[\frac{\partial \rho}{\partial t} = \tfrac{1}{4} \big[ [B_{\alpha}^{\dagger},\rho], B_{\alpha} \big] + \tfrac{1}{4} \big[ [B_{\alpha},\rho], B_{\alpha}^{\dagger} \big] - \tfrac{1}{4} \big[ \{ B_{\alpha}^{\dagger},\rho\}, B_{\alpha} \big] + \tfrac{1}{4} \big[ \{ B_{\alpha},\rho\}, B_{\alpha}^{\dagger} \big] \label{better} \;,\tag{30}\] one finds that

upon assuming that the jump operators are traceless⁹. This reduces to the earlier Eq. (20 ) when the jump operator is Hermitian (i.e. when \(\mathbf{b}^{(\alpha)}\) is real). This can be compactly written as [40] \[\label{pqohtgnj} \frac{\partial \mathbf{a}}{\partial t} = \mathbb{M}\mathbf{a}(t) + \mathbf{v} \; ,\tag{32}\] where the \((D^2 - 1) \times (D^2 - 1)\) matrix \(\mathbb{M}\) and vector \(\mathbf{v} \in \mathbb{R}^{D^2-1}\) have components

Decomposing Eq. (20 ) into a gradient and solenoidal term analogous to Eq. (27 ) one writes \[\frac{\partial \mathbf{a}}{\partial t} = - \text{\scalebox{.85}{\frac{D}{D-1}}} \nabla_{\mathbf{a}} \widetilde{\Phi}(\mathbf{a}) + \widetilde{\boldsymbol{\mathcal{R}}}(\mathbf{a}),\] where \(\widetilde{\Phi}\) is precisely the potential (just expressed in terms of the Bloch vector), and has the general form \[\begin{align} \label{mathcalPR} \widetilde{\Phi}(\mathbf{a}) := \text{\scalebox{.85}{\frac{D}{D-1}}} \left[ \; \text{\scalebox{.85}{\frac{1}{2}}} \mathbf{a}^{T} \mathbb{P} \mathbf{a} + \mathbf{p}^{T} \mathbf{a} \right] \end{align}\tag{33}\] where \(\mathbb{P}\) is an unknown \((D^2-1) \times (D^2 - 1)\) symmetric matrix and \(\mathbf{p}\) is an unknown \((D^2 -1)\)-dimensional vector. The solenoidal vector function \(\widetilde{\boldsymbol{\mathcal{R}}}\) is analogous to \(\widetilde{R}\) in Eq. (27 ). With this definition one then finds that \[\begin{align} \label{mathcalPR2} - \text{\scalebox{.85}{\frac{D}{D-1}}} \nabla_{\mathbf{a}} \widetilde{\Phi}(\mathbf{a}) = - \mathbb{P} \mathbf{a} - \mathbf{p} \qquad \mathrm{and} \qquad \widetilde{\boldsymbol{\mathcal{R}}}(\mathbf{a}) = ( \mathbb{M} + \mathbb{P}) \mathbf{a} + \mathbf{v} + \mathbf{p} \end{align}\tag{34}\] in order to ensure that Eq. (20 ) continues to hold. Enforcing that \(\widetilde{\boldsymbol{\mathcal{R}}}\) is solenoidal with \(\nabla_{\mathbf{a}} \cdot \widetilde{\boldsymbol{\mathcal{R}}} = 0\) additionally implies that \[\label{Tr95cond} \mathrm{Tr}[\mathbb{P}] = - \mathrm{Tr}[\mathbb{M}] = - (D-1) \sum_{\alpha} \; \mathbf{b}^{(\alpha)\ast} \cdot \mathbf{b}^{(\alpha)} \; ,\tag{35}\] where we used the Cartan-Killing form \(f_{njk}f_{mjk} = D \delta_{nm}\) of \(\mathfrak{su}(D)\) and the relation \(f_{njk}d_{mjk}=0\). The orthogonality condition¹⁰ \(\nabla_{\mathbf{a}} \widetilde{\Phi} \cdot \widetilde{\boldsymbol{\mathcal{R}}} =0\) then implies \[\begin{align} \big[ \mathbb{P} \mathbf{a} + \mathbf{p} \big] \cdot \big[ ( \mathbb{M} + \mathbb{P}) \mathbf{a} + \mathbf{v} + \mathbf{p} \big] = 0 \qquad \forall \mathbf{a} \;. \end{align}\] This translates to the set of equations: \[\begin{align} \mathbb{M}^{T} \mathbb{P} + \mathbb{P}\mathbb{M} + 2 \mathbb{P}^2 & =& 0 \tag{36} \\ \mathbb{M}^{T} \mathbf{p} + \mathbb{P} ( \mathbf{v} + 2\mathbf{p} ) & = & \mathbf{0} \tag{37} \\ \mathbf{p} \cdot (\mathbf{v} + \mathbf{p}) & = & 0 \tag{38} \end{align}\] As pointed out in [38], the first equation (for the unknown matrix \(\mathbb{P}\)) is a special case of the algebraic Riccati equation, which often arises in the context of optimal control theory. In Appendix 6, we outline a general algorithm for developing solutions to Eqs. (36 )-(38 ) with the constraint \(\mathrm{Tr}[\mathbb{P}] = - \mathrm{Tr}[\mathbb{M}]\) from Eq. (35 ), and provide an explicit example in §3.5.2.

3.5 Examples with non-Hermitian jump operators↩︎

We here provide two explicit examples. First we solve for qubit gradient flow in the case of a non-Hermitian jump operator, and then we consider the orthogonal Helmholtz decomposition for a particular qutrit problem.

3.5.1 Qubits with non-Hermitian jump operators↩︎

Let us consider the qubit \(D=2\) case again, in the case of a single jump operator which we take to be traceless. Note that qubits are special because the symmetric \(d\)-coefficients (from Eq. (18 )) vanish in this special case. The terms in Eq. (31 ) involving \(\star\) therefore vanish and \(\wedge = \times\) in this case, giving \[\begin{align} \mathbf{a}'(t) & = & \tfrac{1}{2} \mathrm{Re} \left[ \mathbf{b} \times (\mathbf{b}^{\ast} \times \mathbf{a} ) \right] + \tfrac{1}{4} \mathbf{b}^{\ast} \times (\mathbf{b} \times \mathbf{a} ) + \tfrac{1}{2} \mathrm{Im}[ \mathbf{b}^{\ast} \times \mathbf{b} ] \;. \end{align}\] In this simple case, one finds that this is exactly gradient flow, where \[\begin{align} \mathbf{a}'(t) & = & - \nabla_{\mathbf{a}} \Phi \qquad \mathrm{with}\;\Phi = \Big( \tfrac{1}{4} \big[ \mathbf{b} \times \mathbf{a}(t) \big] \cdot \big[ \mathbf{b}^{\ast} \times \mathbf{a}(t) \big] + \tfrac{1}{2} \mathrm{Im}[ \mathbf{b}^{\ast} \times \mathbf{b} ] \cdot \mathbf{a} \Big) \;. \end{align}\] By inspection, one clearly sees that the potential is minimized at a point proportional to \(\mathrm{Im}[\mathbf{b}^{\ast} \times \mathbf{b}]\) as we verify more explicitly below.

The above can be written as a matrix equation \[\begin{align} \mathbf{a}'(t) = \mathbb{M} \mathbf{a} + \tfrac{1}{2} \mathrm{Im}[\mathbf{b}^{\ast} \times \mathbf{b} ] \qquad \qquad \mathrm{with} \;\mathbb{M} = \tfrac{1}{4} \mathbf{b}^{\ast} \mathbf{b}^{T} + \tfrac{1}{4} \mathbf{b} \mathbf{b}^{\dagger} - \tfrac{1}{2} ( \mathbf{b}^{T} \mathbf{b}^{\ast} ) \mathbb{I} \end{align}\] Diagonalizing the matrix, one finds eigenvalues \(-\Gamma_j\) and eigenvectors \(\mathbf{r}_j\) such that \[\begin{align} -\Gamma_{1} & = - \tfrac{1}{2} \mathbf{b} \cdot \mathbf{b}^{\ast} \\ -\Gamma_{2} & = - \tfrac{1}{4} \mathbf{b} \cdot \mathbf{b}^{\ast} - \tfrac{1}{4} \sqrt{ ( \mathbf{b} \cdot \mathbf{b} ) ( \mathbf{b}^{\ast} \cdot \mathbf{b}^{\ast} ) } \\ -\Gamma_{3} & = - \tfrac{1}{4} \mathbf{b} \cdot \mathbf{b}^{\ast} + \tfrac{1}{4} \sqrt{ ( \mathbf{b} \cdot \mathbf{b} ) ( \mathbf{b}^{\ast} \cdot \mathbf{b}^{\ast} ) } \end{align}\begin{align} \mathbf{r}_{1} & = \mathbf{b} \times \mathbf{b}^{\ast} \\ \mathbf{r}_{2} & = \mathbf{b} - \tfrac{\mathbf{b} \cdot \mathbf{b} }{ \sqrt{ ( \mathbf{b} \cdot \mathbf{b} ) ( \mathbf{b}^{\ast} \cdot \mathbf{b}^{\ast} ) } } \mathbf{b}^{\ast} \\ \mathbf{r}_{3} & = \mathbf{b} + \tfrac{\mathbf{b} \cdot \mathbf{b} }{ \sqrt{ ( \mathbf{b} \cdot \mathbf{b} ) ( \mathbf{b}^{\ast} \cdot \mathbf{b}^{\ast} ) } } \mathbf{b}^{\ast} \\ \end{align}\] Notice that in the case of a Hermitian operator, the first two eigenvalues both become \(- |\mathbf{b}|^2/2\) and the last eigenvalue becomes 0. Meanwhile the first two eigenvectors formulae stop working as the above formulae predict the eigenvectors are the zero vector, while the last eigenvector is just \(\mathbf{b}\).

The determinant of \(\mathbb{M}\) is non-zero and the matrix is invertible, this means that the steady state Bloch vector \(\mathbf{a}_{\mathrm{SS}}\) is given by \[\begin{align} \mathbf{0} = \mathbb{M} \mathbf{a}_{\mathrm{SS}} + \tfrac{1}{2} \mathrm{Im}[ \mathbf{b}^{\ast} \times \mathbf{b}] \qquad \implies \qquad \mathbf{a}_{\mathrm{SS}} = - \mathbb{M}^{-1} \; \tfrac{1}{2} \mathrm{Im}[ \mathbf{b}^{\ast} \times \mathbf{b}] = \frac{ \mathrm{Im}[ \mathbf{b}^{\ast} \times \mathbf{b}] }{\mathbf{b} \cdot \mathbf{b}^{\ast} } \end{align}\] The length of the steady state (related to the purity via Eq. (17 )) is \[\begin{align} |\mathbf{a}_{\mathrm{SS}} |^2 = 1 - \frac{ ( \mathbf{b} \cdot \mathbf{b}) ( \mathbf{b}^{\ast} \cdot \mathbf{b}^{\ast} ) }{ ( \mathbf{b} \cdot \mathbf{b}^{\ast})^2 } \;. \end{align}\] The general solution to \(\mathbf{a}'(t) = \mathbb{M} \mathbf{a} + \tfrac{1}{2} \mathrm{Im}[\mathbf{b}^{\ast} \times \mathbf{b} ]\) is given by \[\begin{align} \mathbf{a}'(t) & = & \mathbf{a}_{\mathrm{SS}} + e^{ \mathbb{M} t } \big( \mathbf{a}(0) - \mathbf{a}_{\mathrm{SS}} \big) \;= \; \frac{ \mathrm{Im}[\mathbf{b}^{\ast} \times \mathbf{b} ] }{\mathbf{b} \cdot \mathbf{b}^{\ast} } + \sum_{j=1}^3 e^{ - \Gamma_j t } \frac{ \mathbf{r}_j^{\ast} \cdot \mathbf{a}(0) }{ \mathbf{r}^{\ast}_j \cdot \mathbf{r}_j } \mathbf{r}_j \end{align}\] this shows that the solution relaxes towards the steady state, with timescales \(1/\Gamma_j\). It is interesting that one can arbitrarily tune the jump operator \(B\) to tune the final steady state — this has potential application for thermal state preparation.

3.5.2 Qudits with non-Hermitian jump operators↩︎

In this subsection, we demonstrate the use of the framework in the calculation of a qudit undergoing Lindblad evolution with the single jump operator \(B\) with components \(B_{jk} = \gamma \sqrt{j} \delta_{j,j+1}\), where \(\gamma > 0\) is a constants with units of energy. This corresponds to amplitude damping [41], where the jump operator is the annihilation operator for a \(D\)-level system. Using the standard Gell-mann matrices one finds for example with \(D=3\) that \[B = \gamma \left[\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & \sqrt{2} \\ 0 & 0 & 0 \end{matrix} \right] \qquad \implies \qquad \mathbf{b} = \frac{\sqrt{3}}{2} \gamma \; \left( 1, i, 0, 0, 0, \sqrt{2}, \sqrt{2} i, 0 \right) \;.\] The Lindblad equation for this choice of jump operator can be written in terms of Bloch vectors as the equation \[\frac{{\rm d}\mathbf{a}}{{\rm d}t} = \mathbb{M} \mathbf{a}(t) + \mathbf{v} \quad \mathrm{with} \;\mathbb{M} = \left[ \begin{smallmatrix} -\gamma ^2/2 & 0 & 0 & 0 & 0 & \sqrt{2} \gamma ^2 & 0 & 0 \\ 0 & -\gamma ^2/2 & 0 & 0 & 0 & 0 & \sqrt{2} \gamma ^2 & 0 \\ 0 & 0 & -\gamma ^2 & 0 & 0 & 0 & 0 & \sqrt{3} \gamma ^2 \\ 0 & 0 & 0 & -\gamma ^2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -\gamma ^2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - 3 \gamma ^2/2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 3 \gamma ^2/2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 \gamma ^2 \\ \end{smallmatrix} \right] \;\mathrm{and} \; \mathbf{v} = \left[\begin{smallmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \gamma ^2 \end{smallmatrix} \right] \;.\] Note that \(\det \mathbb{M} = \frac{9}{8} \gamma^{16} \neq 0\) which means the matrix \(\mathbb{M}\) is invertible: the steady-state is clearly given by \[\label{dampSS} \mathbf{a}_{\mathrm{SS}} = - \mathbb{M}^{-1} \mathbf{v} = \left[\begin{smallmatrix} 0 \\ 0 \\ {\sqrt{3}}/{2} \\ 0 \\ 0 \\ 0 \\ 0\\ 1/2 \end{smallmatrix} \right]\tag{39}\] which using Eq. (16 ) means that the steady-state is given by the density matrix \(\require{upgreek} \rho_{\mathrm{SS}} = \frac{1}{3} \mathbb{I} + \tfrac{1}{2} \uplambda_{3} + \tfrac{1}{2\sqrt{3}} \uplambda_{8}=\mathrm{diag}(1,0,0)\) corresponding to a pure state, in which the system ends up in the ground state of the computational basis.

Note however, that one may use instead the orthogonal Helmholtz-Hodge decomposition for studying the behavior of the solution. Indeed one finds that \[\mathbb{P} = \left[ \begin{smallmatrix} {\gamma ^2}/{3} & 0 & 0 & 0 & 0 & -{\gamma ^2}/{3 \sqrt{2}} & 0 & 0 \\ 0 & {\gamma ^2}/{3} & 0 & 0 & 0 & 0 & -{\gamma ^2}/{3 \sqrt{2}} & 0 \\ 0 & 0 & {3 \gamma ^2}/{4} & 0 & 0 & 0 & 0 & -{ \sqrt{3} \gamma ^2 }/{4} \\ 0 & 0 & 0 & \gamma ^2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \gamma ^2 & 0 & 0 & 0 \\ - {\gamma ^2}/{3 \sqrt{2}} & 0 & 0 & 0 & 0 & {5 \gamma ^2}/{3} & 0 & 0 \\ 0 & -{\gamma ^2}/{3 \sqrt{2}} & 0 & 0 & 0 & 0 & {5 \gamma ^2}/{3} & 0 \\ 0 & 0 & - {\sqrt{3} \gamma ^2}/{4} & 0 & 0 & 0 & 0 & {9 \gamma ^2}/{4} \end{smallmatrix} \right] \qquad \mathrm{and} \qquad \mathbf{p} = \left[ \begin{smallmatrix} 0 \\ 0\\ -{\sqrt{3} \gamma ^2}/{4} \\ 0 \\ 0 \\ 0 \\ 0 \\ - {3 \gamma ^2}/{4} \end{smallmatrix} \right]\] solves Eqs. (36 )-(38 ) for this particular choice of jump operator. Using Eq. (33 ) this means that the potential is \[\begin{align} \widetilde{\Phi}(\mathbf{a}) & := & \text{\scalebox{.85}{\frac{3}{2}}} \left[ \; \text{\scalebox{.85}{\frac{1}{2}}} \mathbf{a}^{T} \mathbb{P} \mathbf{a} + \mathbf{p}^{T} \mathbf{a} \; \right] \\ &=& \frac{\gamma ^2}{16} \bigg(4 a_1^2-4 \sqrt{2} a_6 a_1+4 a_2^2+9 a_3^2+12 a_4^2+12 a_5^2+20 a_6^2\\ && +20 a_7^2+27 a_8^2-6 \sqrt{3} a_3-4 \sqrt{2} a_2 a_7-6 \sqrt{3} a_3 a_8-18 a_8\bigg) \notag \end{align}\] It is easy to check that this potential is minimized by exactly the steady state \(\mathbf{a}_{\mathrm{SS}}\) given in Eq. (39 ), as required.

4 Conclusions↩︎

In this work, we have reformulated the finite-dimensional Lindblad equation as gradient flow. For Hermitian jump operators, the evolution admits a simple gradient flow structure with respect to a scalar potential \(\Phi\), offering clear insight into the dynamics and steady-state behaviour. In the non-Hermitian case, the steady-state structure is still fixed by a potential \(\Phi\), but the presence of a solenoidal term \(R(\rho)\) complicates the dynamics. We demonstrated that one can still cast the equation in a gradient-like form via (i) a complexication of the potential or (ii) an orthogonal Helmholtz-Hodge decomposition. There is still work to be done in better understanding the more complicated dynamics of the non-Hermitian jump operator case.

The algorithm described in Appendix 6 for constructing an orthogonal Helmholtz-Hodge decomposition is generic, making it well suited for numerical implementation, but must still be applied on a case-by-case basis. In practice, finding explicit expressions for the orthogonal choices of \(\widetilde{\Phi}\) and \(\widetilde{R}\) is difficult for generic jump operators. We leave the development of a general analytic method for this task to future work.

Throughout this work, we have phrased everything in the interaction picture, absorbing the Hamiltonian into the definition of the jump operators. In realistic settings, interactions with an environment generically shift the free Hamiltonian, and incorporating these shifts into the interaction-picture jump operators is not always straightforward. As a result, the gradient flow structure may be modified by an additional von Neumann term, which can always be absorbed into the \(R(\rho)\) contribution.

There are several other potentially interesting directions to explore:

• Thermal State Preparation: It would be valuable to understand whether thermal steady states can always be engineered by appropriately choosing the jump operators, and how special this feature is to qubits.

• Numerical Applications: Since gradient-flow equations underlie many optimization methods (like gradient descent widely used in machine learning), they may provide useful intuition for numerical strategies. One could speculate that this structure might help in approximating or truncating the infinite hierarchy of equations that are often encountered in studies of adjoint master equations and mean-field approximations [42]. We leave this as an interesting direction for future work.

• Alternative Metrics: We assumed a flat distance metric throughout in this work. Extending our framework to incorporate more physically natural metrics — such as the Bures or Fubini-Study metrics — could more clearly connect our approach to the entropic gradient flows studied in [12]–[15]. The connection to the exact renormalization group (ERG) flows studied in [43] could also be interesting. There, the ERG was shown to be a Lindbladian quantum channel, and monotonicity of distinguishability followed from the data processing inequality. It would be interesting to relate this to our scalar potential \(\Phi\) and its interpretation as a Lyapunov function.

• Eigenstate Thermalization Hypothesis (ETH): In [44], Lindblad dynamics were derived from coupling to a finite-sized environment obeying ETH (see also [45]). It would be interesting to explore whether our gradient flow formalism applies to such settings or to broader non-equilibrium problems related to ETH [46].

• Non-Markovian extensions: We suspect that our approach can be generalized to non-Markovian systems. Non-Markovian master equations are notoriously difficult to solve, and in particular we suspect that the study of late-time states can be simplified by use of the formalism presented here. This approach could complement recent results such as those in [40].

• Infinite-dimensional systems: Extending our formulation to infinite-dimensional systems involving quantum field theories would be very interesting as understanding steady states in these problems tends to be very challenging (for example in inflationary cosmology [47]–[49]). It may be tractable to export the tools presented here, as well as those involving Bloch vectors to this case using known results about \(\mathfrak{su}(\infty)\) [50].

We hope that our perspective on Lindblad evolution provides a useful framework for further analytical and numerical investigations across quantum information and open systems.

Acknowledgements↩︎

We thank Patrick Hayden, Alejandro Kunold, and Arnab Pradhan for useful discussions. The authors are grateful to William Kretschmer and Shozab Qasim for comments on a preliminary version of the manuscript.

References↩︎