Super-operator structures and no-go theorems for dissipative quantum phase transitions

In the thermodynamic limit, the steady states of open quantum many-body systems can undergo nonequilibrium phase transitions due to a competition between coherent and driven-dissipative dynamics. Here, we consider Markovian systems and elucidate structures of the Liouville super-operator that generates the time evolution. In many cases of interest, an operator-basis transformation can bring the Liouvillian into a block-triangular form, making it possible to assess its spectrum. The spectral gap sets the asymptotic decay rate. The super-operator structure can be used to bound gaps from below, showing that, in a large class of systems, dissipative phase transitions are actually impossible and that the convergence to steady states follows an exponential temporal decay. Furthermore, when the blocks on the diagonal are Hermitian, the Liouvillian spectra obey Weyl ordering relations. The results apply, for example, to Davies generators and quadratic systems, and are also demonstrated for various spin models.

Due to practical constraints and our aim of manipulating these systems efficiently, they are inevitably open in the sense that they are coupled to the environment. For the trapped atom and ion systems, the employed electronic state space is not entirely isolated. One has, for example, heating of motional degrees of freedom, photoionization, and spontaneous emission, as well as fluctuations in laser beam intensities and pulse timings. Such environment couplings naturally lead to dissipation and decoherence.
As an open system is generally entangled with its environment, it has to be described by a density operator ̺(t). In Markovian systems, the dynamics is then governed by a Lindblad master equation ∂ t ̺ = L(̺) with [15][16][17][18][19] In addition to the unitary Hamiltonian part, the Liouville super-operator L contains dissipative terms, where Lindblad operators L k describe the environment couplings with strengths γ k ≥ 0. At long times, systems with a finitedimensional Hilbert space converge to steady states ̺ ss obeying L(̺ ss ) = 0 or evolve unitarily in an asymptotic subspace [20]. All eigenvalues of L have non-positive real parts [18], and the largest nonzero eigenvalue real part sets the asymptotic decay rate -the so-called Liouville gap ∆. The exceptional control in the aforementioned experimental systems can be used to design and tune the environment couplings. In this way, one could stabilize useful states [21][22][23] such as cluster or more general graph states [24,25] for measurement-based quantum computation [26,27], or suitable initial states for quantum phase estimation [28,29] and quantum simulation. Moreover, in a competition of coherent and driven-dissipative dynamics, one could drive the system into (novel) phases of matter that may not be accessible in equilibrium.
Here, we investigate conditions under which phase transition in the steady state can or cannot occur on a fundamental level. In the literature, the theoretical analysis is often based on mean-field approximations: For many-body lattice systems, interaction terms are decoupled in the mean field theory to arrive at nonlinear equations of motion for local observables. The mean-field solutions may suggest rich phase diagrams and can include multi-stable and oscillatory phases [22,[30][31][32][33][34][35][36][37][38][39][40][41][42][43]. The validity of such results is often unclear as mean field theory is notoriously unreliable for dynamical problems, even for closed systems. In several cases, qualitative discrepancies have been found in comparison with other approaches [36,[44][45][46][47][48][49][50][51][52][53].
Like quantum phase transitions in closed systems [54], a dissipative phase transition is characterized by a non-analytic dependence of steady-state expectation values on system parameters. This, in turn, requires a non-analytic change in the steady-state density matrix and, hence, a level crossing [55]. So, the spectral gap ∆ of the Liouvillian needs to close at the transition point [56,57].
In this paper, we describe classes of Liouvillians that can be brought into a block-triangular form by operator-basis transformations. These transformations are in general nonorthogonal and can, for example, be found using symmetries and other dynamical constraints. The Liouvillian spectrum can then be assessed through the spectra of the blocks on the diagonal. If the blocks are all Hermitian, the spectrum is completely real. The eigenvalues of sums of such Liouvillians obey Weyl ordering relations, which can result in bounds on the gap and, hence, preclude steady-state phase transitions (Sec. III). This class of open systems encompasses, for example, the Davies generators [58][59][60][61] that arise in the context of thermalization for systems that are weakly coupled to thermal baths. We then generalize further to Liouvillians that can be transformed to block-triangular form but have non-Hermitian blocks on the diagonal, and we discuss ways to bound their spectral gaps (Secs. IV, V, and VI). All results are illustrated with exemplary many-body lattice models. The block-triangular structures can make (seemingly) perturbative treatments exact in the sense that effective Liouvillians co-incide with specific diagonal blocks in the block-triangular representation and, hence, yield exact Liouvillian eigenvalues. Exploiting the block structure also makes numerical approaches more efficient and can lead to integrable subproblems.

II. A DISSIPATIVE MODEL WITH Z2 SYMMETRY BREAKING
Let us start the discussion with a concrete spin-1/2 system on a cubic lattice of N sites in d dimensions. It serves as a motivation and example for the general discussion in Sec. III. We can try to mimic the physics of the non-dissipative transverse Ising model at zero temperature. The paramagnetic state |⇒ := |→→ . . . → with σ x i |⇒ = |⇒ can be stabilized by Lindblad operators L x i := σ x+ i . The ferromagnetic states |⇑ := |↑↑ . . . ↑ and |⇓ are stabilized by controlled spin-flip operators acting on bonds (i, j). Here, σ α i are the Pauli matrices on site i, σ x+ = |→ ←|, and P s i projects onto σ z or σ x eigenstates |s on site i with s ∈ {↑, ↓, →, ←}. For simplicity, we include no Hamiltonian terms such that and The sum in Eq. (2c) runs over all lattice bonds, and The term D z has been included to make |⇑ ⇑| and |⇓ ⇓| the only steady states for γ x → 0. The system has a Z 2 symmetry generated by U x = i σ x i .

A. The disordered phase
For γ x ≫ γ f , γ z , the system is in the disordered (paramagnetic) phase. At γ f = γ z = 0, the non-degenerate steady state is |⇒ ⇒| = i P → i . The first excitations can be generated by flipping one spin in the bra or ket and have eigenvalue λ = −γ x /2. Noting that σ y and σ z are linear combinations of |← →| and |→ ←|, we choose the operators as a basis for the subspace of first excitations such that The corresponding left eigenvectors of L are B z i = σ z i /2 and B y i = σ y i /2. They obey the biorthonormality conditions Taking D ± and D z into account perturbatively, the extensive degeneracy of the first excitations is lifted. For a ddimensional cubic lattice, the only nonzero matrix elements of the effective Liouvillian So, with site positions r j ∈ Z d and α ∈ {z, y}, the excitations become A α k = 1 N j e ik·rj A α j with crystal momentum k and dispersions up to higher-order corrections. Thus, to leading order, the Liouvillian gap stays fixed with −λ z k=0 = γ x /2.
Numerically, one also finds an isolated eigenvalue λ b inbetween the first excitation and the single-magnon continuum. For small γ x , it obeys −γ x /2 > λ b > −2γ z − γ x /2. The nature of this excitation will be clarified in Sec. III D.

C. Further observations and interpretation
Beyond perturbation theory, it turns out that all eigenvalues of the Liouvillian (2) are real and that, surprisingly, the gap is fixed at γ x /2 in the entire phase diagram. At this point, it is evident that the analogy to the Hamiltonian transverse Ising model does not carry very far here: We have an isolated transition point at γ x = 0 without a continuum of gapless excitations. There is Z 2 symmetry breaking at this point, but we do not actually have an extended ferromagnetic phase. Also, the excitations above the ferromagnetic steady states are magnons instead of domain walls. While this can be rectified in a modified model, the gap being fixed at γ x /2 in this system is remarkable and lead to an interesting general finding: The unapparent cause for the observed spectral properties is that one can apply a non-orthogonal basis transformation such that D x , D ± , and D z assume a block-triangular form with Hermitian blocks on the diagonal. This has far-reaching consequences, as discussed in the following.

III. LIOUVILLIANS WITH REAL SPECTRUM
The reason for the gap in the model (2) being finite for all nonzero γ x and further properties of its Liouvillian spectrum are closely related to the fact that the eigenvalues of the competing terms γ x D x and γ f D f + γ z D z are all real. Let us discuss, more generally, conditions under which Liouvillians have an entirely real spectrum.

A. Hermiticity
The simplest scenario are Liouvillians that are self-adjoint L = L † and, hence, have a real spectrum. Here the adjoint of a super-operator M is defined through the Hilbert-Schmidt inner product in Eq. (4) such that B|M(A) = M † (B)|A . In particular, the adjoint of the Liouvillian (1) generates the time evolution in the Heisenberg picture, ∂ t B = L † (B). To construct a self-adjoint Liouvillian, one can for example choose H = 0 and dissipators such that for every Lindblad operator L k with rate γ k , there exists a second term with Recalling that the spectrum of a matrix is invariant under similarity transformations, a more general class of models is captured by the following.

Lemma 1. L has a real spectrum if we can find an (in general, non-orthogonal) basis in which it is Hermitian, i.e., if there exists an invertible super-operator G such that
For the orthonormal local operator basis (|→ →|, |← ←|, |← →|, |→ ←|), it is non-Hermitian; hence, it is non-Hermitian in all orthonormal bases. In particular, it has the off-diagonal element In contrast, choosing the non-orthogonal right basis it attains the diagonal form and is hence Hermitian. Here, [D x i ] Bx denotes the matrix with elements B m |D x i |A n for A n ∈ B x and corresponding left basis elements B m ∈ (½, −|← ←|, σ z /2, σ y /2) that obey biorthonormality B m |A n = δ m,n .
A well-known class of Liouvillians, covered by Lemma 1, are Davies generators [58][59][60][61][62][63][64][65][66][67]. They arise when a system is weakly coupled to the environment through Hermitian system operators S k and bath operators, and when the bath is in thermal equilibrium. The Lindblad operators L k (ω) are then given by the Fourier coefficients of S k , i.e., e iHt S k e −iHt =: with the temperature β −1 being set by the environment state. As L k (ω) connects H eigenstates with energy difference ω, it obeys e −βH L k (ω) = e βω L k (ω)e −βH . It follows that Davies generators can describe thermalization in the sense that the canonical ensemble ̺ β = e −βH /Z is a steady state and that the so-called detailed-balance condition is obeyed. This shows that Lemma 1 applies with G = G 1/2 β . Note that our previous example, D x , is not a Davies generator; for example, its steady state does not have full rank.

B. Block-triangular with Hermiticity
In order to understand the observations made for the spin model in Sec. II, we need to generalize further. Proposition 1. If we can find an (in general, non-orthogonal) operator basis in which L is block-triangular with Hermitian blocks M i on the diagonal (BTH), i.e., if there exists an invertible super-operator G such that with M i = M † i , then the spectrum of L is real, consisting of the eigenvalues of the matrices M i .
on a single bond with respect to the operator basis B ′ x on both sites. The basis elements, given in Eq. (14), are here denoted by ∅, x, z, and y. This is due to the fact that the spectrum of a triangular matrix is given by its diagonal elements and, in Eq. (13), GLG −1 is transformed to lower or upper triangular form when diagonalizing the blocks M i . For the typical many-body systems that we consider, it is of course practically impossible to transform the Liouvillian to Jordan normal form (or to diagonalize all blocks M i ). As we will see in examples, the idea behind Proposition 1 is to use symmetries or other properties of the system to attain a BTH form (13).
As a simple example we can again invoke the dissipator D . For a single site, it is upper triangular in the basis (|→ →|, |← ←|, |← →|, |→ ←|) with the diagonal elements (0, −1, −1/2, −1/2) giving the real spectrum. Hence, the full dissipator D x = i D x i is also upper triangular when using the corresponding tensor product basis.
A more intricate example is the dissipator D f = D + + D − of the spin system in Eq. (2). D f can be brought into a lower block-triangular form when using B x [Eq. (10)] as the singlesite operator basis. And if we modify the basis slightly to all blocks on the diagonal become Hermitian, such that showing that its spectrum is real. The nonzero matrix elements are given in Table I. The global BTH structure of the many-body system will be discussed in more detail in Sec. III D.

C. Weyl ordering and a no-go theorem for phase transitions
Let the Liouvillians L = L (1) + L (2) consist of two terms that are BTH with the same block structure and have descendingly ordered eigenvalues 0 = λ Then, the eigenvalues 0 = λ 1 ≥ λ 2 ≥ . . . of L are also real and obey the ordering relations This follows by application of Weyl's theorem [68] to the blocks on the diagonal. In particular, if we have two Hermi-tian n × n matrices M and N , and denote their descendingly ordered eigenvalues by λ k (M) and λ k (N ), then one can use the variational characterization of eigenvalues according to the minmax principle [68] to show that λ k (M) + λ n (N ) ≤ λ k (M + N ) ≤ λ k (M) + λ 1 (N ) for all k ∈ [1, n]. For Liouvillians with a real spectrum, the largest eigenvalue is necessarily λ 1 = 0, corresponding to a steady state, and we thus obtain Eq. (15). With this, we arrive at the following no-go theorem for dissipative phase transitions. Proposition 2. Consider L(γ) = γ 0 L (0) + ν γ ν L (ν) with a basis in which the Liouvillians L (k) are all BTH with the same block structure. The spectra of L(γ) and its components then obey the Weyl ordering relations (15). Specifically, if L (0) has a unique steady state and finite gap ∆, then the gap of L(γ) is at least γ 0 ∆, and the system can only be critical for γ 0 = 0. As discussed in the following section, this result explains, for example, why the gap in the spin-1/2 model (2) cannot be smaller than γ x /2 in the entire phase diagram.
(2)], which features Z 2 symmetry breaking at the ferromagnetic point γ x = 0, and let us explain its properties in the light of Proposition 2. In the non-orthogonal single-site basis (14), the dissipators D x = i D x i and D z = i D z i become diagonal and lower triangular, respectively, where Also, D f assumes lower BTH form with the nonzero elements given in Table I. With the dissipator D x in mind, let us refer to the elements (|→ →|, σ x , 2 1/4 σ z , σ y ) of the operator basis (14) as the vacuum ∅, and x, z, and y particles, respectively. We can then characterize the actions of the three dissipators in that basis as follows. D x is diagonal. D z is diagonal except for terms that generate a z particle from the vacuum (∅ → z). The non-diagonal elements of the reflection-symmetric dissipator D f can create one x particle or a particle pair of type x, y, or z from the vacuum (∅∅ → ∅x, xx, yy, zz). For a bond with an x particle and one vacuum site, D f can let the particle hop or create a particle pair (∅x → x∅, xx, yy, zz). When we have a bond with a y or z particle and one vacuum site, D f can let the particle hop or create an x particle (∅y → y∅, yx; ∅z → z∅, zx). When two particles of equal flavor meet on a bond, D f can transform them into a pair of the other two flavors (xx ↔ yy ↔ zz ↔ xx). When two particles of different flavor meet on a bond, D f can swap their positions (αβ ↔ βα).
Thus, the total number of particles never decreases under the action of L, and the particle number parities P α = N α mod 2 are conserved. So, we can structure the operator space into blocks, ordered by increasing particle number N x + N y + N z and can further decompose into sectors with different parities (P x , P y , P z ). All blocks on the diagonal are Hermitian, and L is lower block-triangular.
The dissipator D x has the spectral gap 1/2 and, according to Proposition 2, the gap of L can hence never fall below γ x /2. This excludes steady-state phase transitions at nonzero γ x . The perturbative treatment in Sec. II A suggested that the gap remains fixed at γ x /2 for all γ f and γ z ; it was associated with the k = 0 momentum state of a single z particle. According to the block-triangular structure described above, the perturbatively determined eigenvalues (6) of L are actually exact: The effective Liouvillian B α i |L|A β j in Eq. (5) simply consists of the single yand z-particle blocks on the diagonal in our chosen basis B ′ x . As pointed out in Proposition 1, their eigenvalues are exact L eigenvalues. (This is, of course, not true for the corresponding eigenvectors.) The largest block eigenvalue −γ x /2 is an exact eigenvalue of L for all γ f and γ z . Finally, we can convince ourselves that it is always the overall largest nonzero L eigenvalue: The single x-particle block on the diagonal also features a cosine dispersion with the largest eigenvalue −γ x − 2γ z corresponding to a state with momentum k a = π for a = 1, . . . , d. For the blocks with total particle number N > 1, an upper spectral bound of −N γ x /2 follows from the diagonal elements due to D x . Thus, the spectral gap of L is indeed exactly γ x /2 for all γ x , γ f , γ z .
Finally, let us comment on the isolated second excitation with eigenvalue λ b above the single-magnon continuum found in the perturbative analysis for small γ x in Sec. II B. It turns out to be the largest eigenvalue of the two-particle block on the diagonal with even parities P x = P y = P z = 0 and corresponds to a two-particle bound state. . For the study of spectral gaps, it would of course be sufficient for the block that contains the first excitation to be Hermitian. If this form is difficult or impossible to attain, other options are to study the singular values of M i or to bound the gap. The latter can, e.g., be done using the eigenvalues of the Hermitian components or Gershgorin's circle theorem, as discussed in the following. Proposition 3. If we can find a (basis) transformation G such that L assumes the block-triangular form

              
, and if we can find negative upper bounds on the largest eigenvalue real parts of the block matrices M i ( = M † i ), then no phase transitions can occur.
Again, the idea is to exploit symmetries or other properties of L to attain such a block-triangular form. The cost for transforming into Jordan normal form generally grows exponentially in the system size. To achieve that the zero eigenvalue is located in the first 1 × 1 block of GLG −1 , the first right basis element should be chosen as the actual steady state or, if it is not known, at least the first left basis element should be the identity. In these cases, the spectral gap ∆ of L is determined by the largest real part of the eigenvalues of the blocks M i>0 , i.e., If the asymptotic subspace (Re λ k = 0) cannot be separated in the block structure from further excitations, one needs to additionally assess the largest nonzero eigenvalue real part of any block that contains an element of the asymptotic subspace. While this generalization is straightforward, for conciseness, we assume the former scenario in the following.
It can also be obtained variationally as Let λ = −∆ + iϕ be the M i eigenvalue with largest real part. If the imaginary part ϕ is nonzero, the corresponding singular value ν = ∆ 2 + ϕ 2 would only provide an upper bound on ∆. But if it can be established that ϕ = 0, then ν = ∆.
In general, it is unfortunately difficult to assess whether λ is real. Even if the smallest singular value is unique, we only know that it corresponds to a real eigenvalue, and, in general, it only gives an upper bound on ∆.
Using the Hermitian component. -A more straightforward approach is to study the spectrum of the Hermitian component (17). The largest and smallest eigenvalues µ max and µ min of the Hermitian component M H i of a Liouvillian block M i bound the real parts of its eigenvalues such that Clearly, if M i (A) = λA with normalized A, then λ = A|M i |A and Re λ = A|M H i |A ∈ [µ min , µ max ]. The quality of the bound depends on the chosen basis, i.e., the transformation G in Proposition 3, and one can try to apply non-orthogonal transformations in the blocks to decrease µ max . In the best case, the Hermitian and anti-Hermitian components of M i commute such that the M H i eigenvalues coincide with the real parts of the M i eigenvalues.
Gershgorin circle theorem. -A well-known theorem for bounding the eigenvalues of matrices is Gershgorin's circle theorem [69][70][71]. It states that each eigenvalue of a square n × n matrix M is contained in one of n disks in the complex plane, centered at the diagonal elements M i,i with radii Alternatively, one can choose the radii as the column sums instead of the row sums such that R i = j =i |M j,i |. Again, the quality of the resulting bound on the largest eigenvalue real part depends on the chosen transformation G. In our case, the circle theorem is only useful for diagonally dominant matrices.
All three approaches allow us to assess the gap of Liouvillians that can be transformed to a block-triangular form (13) through properties of the blocks M i on the diagonal instead of the entire Liouvillian. The issue may then reduce to a simple few-particle problem, allow for the application of techniques for integrable systems, or substantially reduce the cost for numerical methods. One can apply exact diagonalization for individual blocks or employ variational techniques like tensor network states [72][73][74][75][76][77][78][79]. Deriving in such a way lower bounds on the spectral gap, one can exclude dissipative phase transitions.

V. SPECTRAL BOUNDS FOR SPIN-1/2 SYSTEMS WITH SPONTANEOUS EMISSION
Let us illustrate how the spectral gap can be bounded by bringing the Liouvillian into a block-triangular form as in Proposition 3 and by then assessing the largest eigenvalues of the Hermitian components of the blocks M i on the diagonal; cf. Eq. (21). In particular, we consider some simple Markovian lattice models L = H + γD e of two-level systems (spins-1/2, qubits) with various Hamiltonians and spontaneous emission D e = i D e i corresponding to Lindblad operators L i = σ − i , i.e., Such models are, for example, naturally realized in Rydberg atom experiments [8][9][10], where σ − i corresponds to a spontaneous deexcitation of atom i that is accompanied by photon emission. In the following, we set γ = 1.

A. Magnetic fields
If the Hamiltonian only contains a magnetic field term H α := − i h α i σ α i , the issue reduces to a single-site problem. In the Pauli operator basis the dissipator takes the block-triagonal form The single-site Liouvillians H z i + D e i and H y i + D e i for fields in the z and y directions are in basis B given by respectively. They are lower block-triangular, but the occurring 2 × 2 blocks are neither Hermitian nor anti-Hermitian.
In the case of the z field, the eigenvalues of the bottom right 2 × 2 block are −1/2 ± 2ih z i . Their real part agrees with the eigenvalues of the Hermitian component (17), which is simply diag(−1/2, −1/2).
In the case of the y field, the eigenvalues of the central 2 × 2 block are For |h y i | ≥ 1/8, they are complex with real part −3/4. For |h y i | < 1/8, the eigenvalues are real, going to λ ± = −1/2, −1 in the limit h y i → 0. The latter two values agree with the spectral bounds set by the Hermitian component diag(−1, −1/2). In fact, for any field in the xy plane, the eigenvalues of the single-site Liouvillian are given by the expressions in Eq. (26) in addition to λ = 0, −1. This is due to the z rotation symmetry of D e .
For a magnetic field in a generic direction with Thus, in all cases, the gap is bounded by ∆ ≥ 1/2.

B. XX and YY interactions with field
Next, let us discuss systems with XX interactions and fields in x direction on more or less arbitrary lattices, i.e., systems with the dissipator D e as in Eq. (23) and Hamiltonians of the form with scalar coupling coefficients J i,j . We again employ the Pauli single-site basis (24), labeling the elements by (e, z, x, y). The field term −i[σ x i , · ] only generates transitions y ↔ z on site i. The XX interaction term −i[σ x i σ x j , · ] only generates transitions ez ↔ xy and ey ↔ xz (and reflections) on bond (i, j). Thus, the number N y + N z of y and z sites is conserved under the Hamiltonian. This implies that the Liouvillian is lower block-triangular if we choose the identity as the first basis element and order all further basis elements by increasing N y + N z .
The first block M 0 on the diagonal is 1 × 1, giving the zero eigenvalue of the steady state. The basis for the second block M 1 on the diagonal consists of all product operators with zero y and z sites (N y + N z = 0) and N x ≥ 1 x sites. As M 1 does not contain any Hamiltonian matrix elements, it is diagonal, giving eigenvalues −N x /2. Here, each x site contributes with −1/2 according to Eq. (25).
The basis for the third block M 2 consists of all product operators that are dynamically connected to the operator with N y + N z = 1 and N x = 0. The Hamiltonian cannot move the y/z excitation, but sites that are connected to the y/z site through the XX interaction can oscillate between e and x as discussed above. The Hermitian component M H 2 only contains D e terms and is hence diagonal. The largest M H 2 eigenvalue is −1/2, corresponding to a single y site. With the lattice coordination number ζ, the smallest M H 2 eigenvalue is −1 − ζ/2, corresponding to one z site and ζ connected x sites. For small ζ, the exact spectrum of the blocks with N y + N z = 1 can be determined analytically.
The blocks with N y + N z > 1 contain multiple y/z sites. These are still immobile but can now interact via their clouds of surrounding x sites if the "distance" is less than twice the interaction range. The Hermitian components (17) of the block matrices are again diagonal, with the largest eigenvalue being −(N y + N z )/2.
In conclusion, the largest eigenvalue real part for every Liouvillian block M n>0 is bounded from above by −1/2. Thus, irrespective of lattice geometry and (in)homogeneity of couplings J i,j and fields h i , the gap is always ∆ = 1/2, corresponding to a single x site. Due to the symmetry of D e , this result readily generalizes to interactions and fields that are rotated by any angle around the z axis, such as YY interactions and y fields.
In this basis, the dissipator is upper block-triangular with The Hamiltonian term H = −i[H, · ] conserves the numbers of up sins N k and N b in the ket and bra spaces, where the elements of the site basis (28) have quantum numbers (N k , N b ) = (0, 0), (1, 1), (1, 0), (0, 1). As the off-diagonal term in [D e i ] Bz simultaneously decreases both quantum numbers, the full Liouvillian (23), assumes upper block-triangular form if we order the operator basis by increasing N k + N b . For fixed N k + N b , we can for example order by increasing N k . The block matrices on the diagonal are where the second term comprises the Hamiltonian matrix elements with N k and N b up spins. The excitations (block-matrix eigenvectors) are operators |E k E b | with H eigenstates |E k and |E b from the N k and N b up-spin subspaces. Their eigenvalues are Thus, the spectral gap is again ∆ = 1/2, irrespective of the detailed Hamiltonian. This simple scenario, of a Hamiltonian with a conserved charge and charge-decreasing Lindblad operators can also be used to construct exactly solvable open quantum systems [80,81].

VI. SPECTRAL BOUNDS FOR QUADRATIC FERMION AND BOSON SYSTEMS
References [82,83] discuss another class of open manybody systems for which L can be brought into a blocktriangular form. Consider a Markovian system of identical fermions or bosons with annihilation and creation operators a j and a † j that obey the canonical (anti-)commutation relations. If the Hamiltonian is bilinear in the ladder operators and Lindblad operators are either linear or bilinear and Hermitian in the ladder operators, we call such a system quadratic. In circuit QED systems [13,84,85], for example, linear Lindblad operators arise naturally from photon loss and pump process, while the coupling of cavities can lead to bilinear Lindblad operators [33,86].
As detailed in Ref. [82], one can always express the Liouville super-operator of a quadratic system in terms of a suitable set of ladder super-operators [82,[87][88][89][90] such that it assumes a block-triangular form in the corresponding biorthogonal operator Fock basis. One can efficiently determine or bound the spectra of the smallest blocks, e.g., to establish that a model is gapless. Examples for gapless fermionic and bosonic systems are given in Ref. [83]. Conversely, one can exploit the block-triangular structure to establish that the dissipative gap is at least κ if we add linear Lindblad operators L j = √ 2κ a j for bosons or for fermions. This implies, for example, that, without symmetry constraints beyond invariance under single-particle basis transformations and fermionic particle-hole symmetry, all gapped quadratic systems belong to the same phase [83].

VII. DISCUSSION
In Sec. II, we discussed a non-trivial spin-1/2 model L = γ x D x + γ f (D + + D − ) + γ z D z , where dissipator D x tries to polarize all spins in the +x direction and the D ± dissipators stabilize the all-up and all-down states. The competition between the different terms leads to Z 2 symmetry breaking at the ferromagnetic point γ x = 0. Perturbation theory suggested that the Liouville gap is ∆ = γ x /2, independent of the other system parameters. In Sec. III, we introduced the general class of Liouvillians that, after an appropriate transformation, become block-triangular with Hermitian blocks on the diagonal (BTH). This implies that all eigenvalues are real and that eigenvalues obey Weyl ordering relations. A special case are Davies generators, which arise in the context of thermalization; they have full-rank steady states and are actually entirely Hermitian after an appropriate basis transformation. The model of Sec. II is an example that does not fall in the Davies class and has a non-trivial block-triangular structure. A very simple non-Davies example is a spin-1/2 with Lindblad operator σ − and pure steady state |↓ ↓|. When considering the competition between several BTH Liouvillian terms, the Weyl ordering relations imply that the gap of each individual term is a lower bound for the gap of the full Liouvillian. For a large class of models, this excludes phase transitions. On the other hand, for local couplings, the nonzero gap enables an efficient preparation of the corresponding steady states [91] and makes them stable with respect to perturbations [92].
In Sec. IV, we considered Liouvillians that can be transformed into a block-triangular form but have non-Hermitian blocks on the diagonal. There are several paths to bound the Liouvillian gaps of such systems, e.g., through the block singular values, Gershgorin's circle theorem, or, maybe most effectively, through the Hermitian components of the block matrices. In Sec. V, we saw how the latter approach can be used to establish nonzero gaps for classes of lattice spin models with spontaneous emission, e.g., for Hamiltonians with mag-netic fields and XX or YY interactions with almost no restrictions on lattice geometry, interaction range, and spatial variations.
For now, our methods for finding operator transformations that block-triangularize Liouvillians were based to the use of symmetries and further dynamical constraints as demonstrated in the examples of Secs. II, III D, V, and VI. It would be great to develop more sophisticated approaches.
Of course, there exist models, where the discussed conditions that exclude dissipative phase transitions are not met. A classic example is the phase transition in the driven-dissipative Kerr oscillator [93][94][95][96]. Here, the block-triangular structure of quadratic models (Sec. VI) is spoiled by the quartic Hamiltonian term a † a † aa, which is due to photon-photon interaction. A similar interaction term spoils the block-triangular structure in the driven-dissipative Bose-Hubbard model, and a steadystate phase transition is expected to occur in two dimensions [50,51].
As seen in the examples, the block-triangular structures may imply that perturbative treatments actually yield exact Liouvillian eigenvalues. In general, reducing the spectral problem to the blocks on the diagonal substantially decreases the costs for numerical approaches like exact diagonalization and can make variational techniques like tensor network state methods [72][73][74][75][76][77][78][79] applicable. It is also conceivable that the block structure may, in some models, allow for the application of methods for integrable systems like the Bethe ansatz [97,98], Luttinger liquid theory [99,100], or conformal field theory [101,102].

ACKNOWLEDGMENTS
We gratefully acknowledge helpful discussions with Jianfeng Lu, Hendrik Weimer, and Xin Zhang, and support through U.S. Department of Energy grant DE-SC0019449.