Topology by dissipation

Topological states of fermionic matter can be induced by means of a suitably engineered dissipative dynamics. Dissipation then does not occur as a perturbation, but rather as the main resource for many-body dynamics, providing a targeted cooling into a topological phase starting from an arbitrary initial state. We explore the concept of topological order in this setting, developing and applying a general theoretical framework based on the system density matrix which replaces the wave function appropriate for the discussion of Hamiltonian ground-state physics. We identify key analogies and differences to the more conventional Hamiltonian scenario. Differences mainly arise from the fact that the properties of the spectrum and of the state of the system are not as tightly related as in a Hamiltonian context. We provide a symmetry-based topological classification of bulk steady states and identify the classes that are achievable by means of quasi-local dissipative processes driving into superfluid paired states. We also explore the fate of the bulk-edge correspondence in the dissipative setting, and demonstrate the emergence of Majorana edge modes. We illustrate our findings in one- and two-dimensional models that are experimentally realistic in the context of cold atoms.


Introduction
Symmetries, their spontaneous breaking, and related order parameters were considered for a long time as the paradigm for understanding ordered states of matter. A paradigm shift was initiated in the late 1980s when the Landau-Ginzburg broken-symmetry theory of ordered phases-widely thought to be exhaustive-proved unable to characterize a new kind of phases with no local order parameter: topological phases, or phases with topological order [1,2]. Instead of being distinguished by symmetries, topological phases are characterized by distinct values of a non-local, topological order parameter, and phases transitions occur whenever the topology changes, signaled by discontinuities in this topological invariant. The existence of topological order may be conditioned on the existence of symmetries. However, as long as topological order is present, the underlying system generally exhibits topological features, i.e., features that are robust against arbitrary (symmetry-preserving) quasi-local perturbations. Spectral gap and ground-state degeneracy are typical topological properties which have been theoretically shown to be robust for wide classes of Hamiltonians [3][4][5][6]. Whereas the spectral gap is a property of the bulk-as topological order itself-the ground-state degeneracy generally depends on the boundary conditions imposed at the edges of the system and on the existence of topological defects in the bulk (e.g., vortices in a superfluid). Most importantly, the degeneracy can be traced to the existence of zero-energy modes localized at the edges or bound to topological defects, which are robust topological features as well. These objects can exhibit exotic behavior under spatial exchange (or "braiding") such as non-Abelian statistics [7][8][9], which opens up exciting possibilities for practical applications such as topological quantum memories and topological quantum computation [9][10][11].
The search for topological phases exhibiting quasiparticles with non-Abelian statistics has brought p-wave paired superfluids and superconductors to the forefront of theoretical and experimental condensed-matter research [8,9,[11][12][13]. Such systems have first been studied in two dimensions (2D), where they have been predicted to support topological phases with gapless chiral edges modes and quasiparticles known as Majorana zero modes, giving rise to Ising-type non-Abelian exchange statistics [8,12,14]. Following a seminal paper by Kitaev [15], the focus has moved more recently to networks of one-dimensional (1D) systems, which were shown to allow for similar topological features (non-Abelian statistics, in particular) as genuine 2D systems [16,17]. Recent proposals for solid-state [18,19] and cold atom [20] systems have made it possible for Majorana zero modes to enter the experimental stage, with promising first results in solid-state devices [21][22][23][24] and the perspective of increased future experimental efforts.
In recent years, the quest for topological states was extended to non-equilibrium systems, going beyond the Hamiltonian ground-state scenario. A first step in this direction was taken with periodically driven Hamiltonian systems [25][26][27], in which the time coordinate plays the role of an extra dimension, allowing for the realization of topological invariants with no equilibrium ground-state counterpart. In this work, we focus on a different paradigm in which Hamiltonian unitary dynamics is replaced by specifically designed dissipative dynamics described by a quantum master equation. Such a scenario was originally proposed as a means for quantum state preparation and quantum computation [28,29] and relies on the proper engineering of a coupling of the system to a suitable reservoir. In the context of cold atoms, such reservoir engineering may be seen as a natural extension of the more conventional Hamiltonian engineering, with similar advantages as compared to solid-state systems such as precise microscopic control and tunability. In previous works, we have shown how this concept can be utilized to "cool" or drive ensembles of atomic fermions into topologically ordered states in one [30] and two [31] dimensions in a targeted way, starting from an arbitrary initial state described a density matrix. The analysis of the many-body properties of the phases and phase transitions arising in these examples has revealed similarities but also differences between the physics of topological ground states of Hamiltonians and topological steady states resulting from a purely dissipative evolution.
In this work we put the results obtained in our two previous case studies into a broader theoretical perspective. We provide a framework for investigating nonequilibrium topological states that can be reached by means of engineered dissipation, developing a formalism and physical understanding that can also be used in situations where dissipation occurs as a perturbation. The natural object to study is the density matrix of the system, which does not necessarily correspond to a pure state described by a wave function alone. In the present article we focus on quadratic master equations with the aim of classifying topological states described by density matrices in analogy to the Hamiltonian ground-state scenario. All information contained in the density matrix is then equivalently encoded in the covariance matrix gathering all static singleparticle correlations. By identifying and exploiting the analogy between this object and a quadratic Hamiltonian in a "first-quantized" representation, we demonstrate how to classify topological phases in a non-equilibrium context where mixed states are allowed. Our analysis focuses on both bulk and edge properties.
As compared to a Hamiltonian ground-state scenario, key differences arise from the fact that the dynamics-or the spectral properties of the system-and the properties of its "ground" (steady) state-or the static correlation properties-are not as tightly related as in the Hamiltonian context. As far as the bulk is concerned, this crucial difference manifests itself in the fact that two independent spectral properties must be present to guarantee that the system is in a stable topological state: The first quantity that we identify is the dissipative gap, which corresponds to the slowest damping rate associated with modes belonging to the bulk of the system and is a direct counterpart of the excitation gap of a Hamiltonian spectrum. The second is the purity gap, which describes the purity of the mode belonging to the bulk which is most strongly mixed. Clearly, a purity gap is always present in the Hamiltonian context, since Hamiltonian ground states are by definition pure states. In our context, however, we argue that the system can undergo a topological phase transition if either (or both) of these two different gaps vanishes in a particular parameter regime.
The purity of the state plays a key role not only in the bulk, but also for the edge physics. In the Hamiltonian context, bulk-edge correspondence theorems describe a tight relation between the number of edge zero modes (i.e., modes that are decoupled from the Hamiltonian dynamics and thus do not evolve) found at the interface between two topologically distinct phases and the value of the topological invariant associated with each of the phases [13,[32][33][34][35]. We formulate a dissipative variant of such bulk-edge correspondence: Topological order ensures the existence, at the interface, of a fermionic subspace that is isolated from the bulk (with a dimension determined by the value of the topological invariant on both sides of the interface). However, in contrast to the Hamiltonian case, topology does not guarantee the decoupling of this subspace from the dynamics. As a result, the modes corresponding to this subspace can be either be zero-damping modes-i.e., modes that are decoupled from the dynamics similarly as in the Hamiltonian setting-or emerge as zero-purity modes-i.e., modes that are in a completely mixed state; in which no information can be stored. In the context of engineered dissipation, the simultaneous appearance of both zero-damping and zeropurity modes may give rise to intriguing physical effects, as we discuss in this work.
The fact that actual physical implementations of model Hamiltonians need often be properly modelled as open systems due to particle losses or dephasing, e.g., has been recognized in a number of theoretical works focusing on the stability of the edge modes [36][37][38][39] or on the very definition of topological order in such circumstances [40]. We emphasize that our approach is fundamentally different here, since dissipation does not occur as a perturbation but is rather harnessed as the main resource to generate the dynamics.
Our paper is organized as follows. In section 2, we discuss the dissipative framework of interest: We introduce the concept of "dark states" in a many-body context, and explain the main ideas behind the physical implementation of a dissipative counterpart of Kitaev's quantum wire, thereby illustrating how to engineer more general dissipative evolutions giving rise to superfluid paired states. We also provide both a second-and a first-quantized formulation of quadratic dissipative dynamics, and discuss the key properties that are necessary to understand the bulk and edge physics of Gaussian states in terms of either the corresponding density matrix or the associated covariance matrix. In section 3, we construct a symmetry-based topological classification of driven-dissipative systems using the covariance matrix, and identify relevant topological invariants in one and two dimensions. In section 4, we then apply this framework to identify the classes D and BDI of Altland and Zirnbauer as the symmetry classes that can be dissipatively targeted under physical constraints related to "typical" implementation schemes. As is well known, the edge modes of systems belonging to these two classes are Majorana fermions, which explains the potential of dissipatively induced superfluids to exhibit such modes. We also show that, in two dimensions, the quasi-locality of the dissipative operations acting on the system density matrix alone implies a vanishing Chern number. In section 5, we discuss the fate of the bulk-edge correspondence in the dissipative setting. We also show how to construct dissipative Majorana modes explicitly in a translation-invariant setting, and examine the robustness of such modes in the presence of typical imperfections. Section 6 is devoted to the discussion of the physical role of the dissipative gap for adiabatic manipulations-in particular, for the braiding of dissipative Majorana modes-showing that dissipative Majorana modes exhibit non-Abelian exchange statistics just as their Hamiltonian counterparts. The remainder of the paper deals with illustrative examples of our general framework and results: In section 7, we analyze a "zigzag" dissipative quantum wire exhibiting topological phase transitions of the three possible types allowed by the closure of the dissipative and/or purity gaps. In section 8, we illustrate in a two-dimensional model a dissipative mechanism that makes it possible to obtain unpaired Majorana modes in a topological phase characterized by an even-valued integer topological invariant.

Quantum master equations for many-body systems
The quantum master equation describing the time evolution of the reduced system density matrix ρ is given bẏ The commutator term familiar from the Heisenberg equation describes the coherent dynamics generated by a system Hamiltonian H. The second part, often referred to as Liouville operator or Liouvillian, describes the dissipative dynamics resulting from the interaction of the system with an environment, or "bath". In particular, the set of Lindblad operators (or "quantum jump" operators) i describe the coupling to that bath. The Liouville operator L has a characteristic Lindblad form: The anticommutator term describes dissipation and must be accompanied by fluctuations in order to conserve the "norm" Tr(ρ) of the system density matrix. The corresponding term, where the Lindblad operators act from both sides onto the density matrix, is called "recycling" or "quantum jump" term. Note that we have absorbed the rates κ i associated with each dissipative process into the definition of the Lindblad operators, making them carry dimension of square root of energy. The rates are non-negative, so that the density matrix evolution is completely positive, i.e., the eigenvalues of ρ remain positive under the combined dynamics generated by H and L [41]. The quantum master equation provides an accurate description of a system-bath setting with a strong separation of scales. More precisely, there must be a fast energy scale in the bath (as compared to the system-bath coupling) that justifies to integrate out the bath in second-order time-dependent perturbation theory. If, in addition, the bath has a broad bandwidth, the combined Born-Markov and rotating-wave approximations are appropriate, resulting in equation (1). Such a situation is generically realized in quantum optical few-level systems, e.g., for a laser-driven atom undergoing spontaneous emission. On the other hand, typical condensed matter systems do not display the necessary scale separations to justify a microscopic description in terms of a master equation. In systems with engineered dissipation [42] as investigated in this paper, we are interested in scenarios that share features from both quantum optical systems-in that they are coupled to Markovian quantum baths-and condensed matter systems-in that they dispose of a continuum of spatial degrees of freedom on a lattice. Using the manipulation tools of quantum optics, the validity of a Markovian master equation can be ensured, giving rise to a well-defined microscopic dissipative many-body dynamics. A similar level of microscopic control as obtained in Hamiltonian engineering in a cold atom context can be expected for this "Liouvillian engineering", which therefore is a natural extension of the former to a more general non-equilibrium situation. In this context, dissipation does not occur as a perturbation, but rather as the dominant resource of the many-body dynamics. In particular, here we will consider the case where the Hamiltonian is absent H = 0. Such a scenario can be useful from a practical point of view-for cooling systems into desired states-but also gives rise to interesting new many-body physics.

Dark states
In the long-time limit, a quantum system governed by equation (1) approaches a stationary or steady state ρ f = ρ(t → ∞) which generically corresponds to a mixed state. An interesting situation appears if, instead, the many-body density matrix is driven towards a pure stationary state, ρ f = |ψ D ψ| D [29,43]. In quantum optics, such pure states |ψ D that are obtained as a result of a dissipative evolution are called dark states. Mathematically, such dark states are zero modes of the Liouville operator. More precisely, a dark state is a zero mode shared by all Lindblad operators, The dark-state solution is the unique solution of the Liouville dynamics if (i) there exists exactly one normalized dark state |ψ D , and (ii) there are no stationary solutions other than this dark state [29,44]. In the specific case of interacting Liouville operators (higher than quadratic in the creation and annihilation operators) [29,44] or noninteracting Liouville operators (quadratic in the creation and annihilation operators), the fact that these conditions are satisfied can be proved rigorously. If present, the dynamics described by equation (1) for H = 0 corresponds to a directed motion-in the Hilbert space of the system-into the "sink" provided by the dark state, which is reached independently of the initial density matrix. In recent years, a number of theoretical [28,45] and experimental [46,47] studies have focused on how to construct Liouville operators such that, in the long-time limit, a quantum system reaches a welldefined, pure many-body steady state or exhibits novel phase transitions resulting from the competition between coherent and dissipative dynamics [48][49][50][51]. In particular, in the context of atomic fermions, it has been shown how to engineer number-conserving dissipative dynamics that drives the system into a pure BCS-type paired state in the absence of conservative forces [52,53]. The dissipative pairing mechanism forms a basis a) b)

Fermionic lattice
Auxiliary lattice Reservoir ⌦ +⌦ } Figure 1. (a) Unit cell for dissipation engineering. The lower potential wells correspond to the physical sites, whereas the upper site in between is an auxiliary site. Atoms in the lower sites are coherently coupled to the auxiliary site with opposite Rabi frequency ±Ω. Decay back to the lower sites occurs via spontaneous emission, where energy is released into a surrounding reservoir (see text). If the coupling to the upper level is sufficiently far detuned (∆ Ω), the latter can be integrated out, so that an effective dynamics for the lower sites is obtained. (b) In an optical lattice setup, this unit cell is repeated in a translation-invariant fashion multiple times in a natural way. The quantum wire corresponds to the lower sites-as anticipated above-which are populated by spinless (or spin-polarized) fermions in the cases discussed in this work. Rabi frequencies with alternating signs are realized by a commensurability condition on lattice and drive lasers (see text). Dissipation results from spontaneous phonon emission into a BEC reservoir (light grey).
for the targeted cooling into states with non-trivial topological properties far from thermodynamic equilibrium.

Physical realization
Here we briefly sketch the implementation idea common to the dissipative models studied in this paper. The basic setting consists of an atomic ensemble confined in an optical lattice (the system), which is driven coherently and immersed into a larger reservoir consisting of a different atomic species and playing the role of the dissipative bath. In the case of interest here, the constituents of the system are fermions. In cold atomic gases, the associated spin is realized in terms of hyperfine states, and thus both the cases of spinless and spinful fermions can be meaningfully considered. The bath constituents are chosen as bosonic atoms, so that the conservation of fermion parity in the system is guaranteed.
The working of the driven-dissipative mechanism is best illustrated by the unit cell Λ-configuration displayed in figure 1. The complete driven-dissipative process consists of two steps: The first step is a coherent excitation from the system (lower sites) to the auxiliary site in between. In the example of figure 1, we quasi-locally excite fermions on the system sites i and i + 1 (with annihilation and creation operators a i and a † i corresponding to each site) into an antisymmetric superposition ∼ a i − a i+1 , which can be controlled by a commensurability condition of the driving laser to the standing wave laser generating the optical superlattice (see reference [42] for details): If the driving laser has twice the wavelength of the lattice laser, there is phase shift of π in the effective Rabi frequency from one site to the next, leading to a relative minus sign. The auxiliary level is unstable if coupled to the reservoir. In this case, spontaneous phonon emission into the surrounding bath can occur, thereby giving rise to the second, dissipative step. The atoms are "brought back" to the lower sites in a quasi-local way ∼ a † i + a † i+1 ; since this process is isotropic, there is now a relative plus sign. If this driven-dissipative process is generated using a drive laser that is far detuned (Ω/∆ 1) from the auxiliary site resonance frequency, the auxiliary site can be integrated out and we obtain a Lindblad operator of the form A few remarks are in order: (i) For a driving laser superimposed over the extent of the whole optical superlattice, we obtain translation-invariant Lindblad operators as depicted in figure 1(b), up to system boundaries which are not shown. (ii) The quasilocality of the operators is controlled by the Wannier function overlap between the onsite wave functions involved in the combined excitation and de-excitation processes.
(iii) The Lindblad operators that can be realized in this setting have a generic form is a linear translation-invariant superposition of creation (annihilation) operators with generic properties: The excitation part (A i ) can be controlled to high accuracy-involving, in particular, the control of the relative phases in the superposition-since it is generated by a coherent laser beam. In two dimensions, for example, this allows to imprint angular momentum by choosing relative laser phases in different primitive directions of the lattice. On the other hand, the de-excitation or decay part (C † i ) results from spontaneous emission and is therefore unavoidably isotropic (or s-wave symmetric). (iv) The system particle number is conserved in our dissipative framework. This is reflected in the property [ i ,N ] = 0 for all Lindblad operators, whereN = in i the total particle number operator,n i = a † i a i . This exact microscopic property of the system, which implies an exact conservation of parity, is of importance to the discussion of the possible imperfections that may occur in the dissipative setting after performing approximations. Physically, this property originates from the fact that typical interactions in cold atomic systems are local density-density interactions. In particular, the system and bath constituents will interact via such coupling. On an even more elementary level, the fact that the bath is bosonic provides a further reason for fermion parity conservation. This aspect is crucially different from solidstate implementations which are not perfectly closed systems: There the environment is typically fermionic, which facilitates system-bath exchange processes affecting the parity of the system. (v) While no particle number exchange is possible between the system and the reservoir, energy can be exchanged. This enables the outflow of entropy from the system into the (infinite) reservoir, and the targeted cooling into pure many-body states. A crucial prerequisite for the entropy removal from the system is the coherent driving of the system. (vi) The fast energy scale ensuring the validity of the Born-Markov and rotating-wave approximations underlying our construction is provided by the band separation between the system and the auxiliary sites, which is the largest energy scale in the problem. (vii) The creation and annihilation part of the Lindblad operators is respectively symmetric and antisymmetric under the exchange i ↔ i + 1. Such property is important for reaching pure stationary states, as will become clear below.

From interacting Liouvillians to quadratic master equations
The physical implementation discussed above realizes a number-conserving microscopic dynamics, with the key advantage of conserving fermion parity as a consequence. The dynamics generated by the corresponding bilinear Lindblad operators i is described by an interacting (quartic) Liouvillian. The Lindblad operators are constructed in such a way that the stationary state is a unique dark state for a given even particle number 2N , characterized by a BCS pair wave function with fixed particle number |BCS, N that satisfies i |BCS, N = 0 for all i. The Liouville operator ensuring this property thus represents a parent Liouville operator for a given fixed number BCS pair wave function (see Appendix A for more details). Starting from the exact knowledge of the fixed-number dark-state wave function, we can switch in the thermodynamic limit from a fixed-number to a fixed-phase ensemble. In particular, the long-time evolution of the interacting master equation can be linearized based on this knowledge. The calculation presented in Appendix A can be summarized as That is, the product of creation and annihilation parts in the quadratic Lindblad operators transforms into a sum, giving rise to linear Lindblad operators. This relies on the property that C † i (A i ) is symmetric (antisymmetric). It provides a dynamical connection between the fixed-number and fixed-phase settings at the level of the equation of motion. The long-time dynamics is universal, in the sense that it is independent of the initial state.
We note that α = re iθ is a complex number in the above equation. Its phase is not fixed by the dynamics, but rather reflects spontaneous symmetry breaking in the dissipative setting of interest. The modulus r, on the other hand, is determined by the average particle number in the system (see Appendix A). In particular, for half filling and in the example of equation (3), we find r = 1; such that for θ = 0, without loss of generality, We recognize the quasi-local Bogoliubov quasiparticle operators of Kitaev's Hamiltonian quantum wire [15] (at half filling and with equal pairing and hopping amplitudes, up to normalization), emerging here naturally in the long-time limit of a dissipative dynamics. The ground-state condition of the Hamiltonian quantum wire, L i |BCS, θ = 0 for all i, is now interpreted as the dark-state condition of the linearized dissipative evolution. The corresponding wave function now has a fixed phase θ instead of a fixed number. Since the Lindblad operators form a complete Dirac algebra, i , L † j } = 0 for an infinite system with no boundaries, the uniqueness of the dark-state solution is obvious.
The quadratic dynamics obtained in the long-time limit makes the systems considered in this work amenable to a treatment analogous to the discussion of quadratic Hamiltonians when examining their topological properties. This dynamics was obtained by giving up exact particle number conservation, which is justified in the thermodynamic limit. The absence of exact particle number conservation thus emerges similarly as in the Hamiltonian scenario. There is, on a formal level, however, a crucial difference between dissipative and Hamiltonian settings. While a quadratic number non-conserving BCS Hamiltonian still conserves parity, formally such a property is not present in a quadratic Liouville evolution (for example, single fermions can be ejected into the environment, giving rise to quasiparticle poisoning [38,39]). However, remembering that the microscopic dissipative processes do conserve particle number and thus parity exactly, we can rule out parity non-conserving processes as potential imperfections arising in our scenario. The number non-conserving nature of the system is introduced "within the system only", and there is no exchange of particles with the reservoir. Physically, the number non-conserving processes describe pairwise creation and annihilation out of or into the mean field provided by the system itself.

Gaussian master equations
Having discussed how quadratic fermionic master equations naturally emerge in the longtime limit of interacting Liouville operators, we now summarize some general properties of such master equations. We do this in both a second-and a first-quantized formulation, working with operators or matrices, respectively, as familiar from the Hamiltonian setting [54]. For this discussion, it is useful to work in the real (or Majorana) basis of fermionic quadrature component operators. For a system with N sites, 2N real fermionic modes are introduced according to c 2n−1 = i a n − a † n , c 2n = a n + a † n , {c n , c m } = 2δ n,m , c † n = c n .
The fact that the master equation is quadratic in the fermion operators implies solutions in terms of Gaussian density operators. In the second-quantized formulation, this can be written as where c T = (c 1 , ..., c 2N ) is a column vector defined from the 2N Majorana operators and G is a real antisymmetric matrix (so that iG is Hermitian). Formally, ρ thus has the form of a canonical density matrix ρ c ∼ e −βH for a quadratic Hamiltonian. Instead of working in second quantization, we can move to a first-quantized formulation. As in the Hamiltonian scenario, the latter allows us to discuss symmetry and topological classifications in a more direct way. The key object here is the covariance or (equal-time) correlation matrix collecting the second moments, which is the only information contained in a Gaussian density operator. It is defined as We now look for an equation of motion for this object [49,55]. A straightforward way to derive such an equation is via the adjoint equation to a given master equation for the density operator [56], describing the evolution of an operator in the Heisenberg picture. For the correlation operatorΓ nm = i 2 2N j,k=1 c j G nm jk c k with real antisymmetric G nm jk = (δ j,n δ m,k − δ j,m δ n,k ), this reads Here we have written the linear quantum jump operators as where X = X T and Y = −Y T are real symmetric and antisymmetric matrices, respectively. Furthermore, X by construction is positive semidefinite. Taking the expectation value of equation (9), we readily find the fluctuation-dissipation equation describing the evolution of the real antisymmetric correlation matrix Γ, where we have suppressed the matrix indices. Denoting the steady-state correlation matrix asΓ, which satisfies the equation {X,Γ} = Y , we can give a clear physical meaning to the matrices X and Y . Writing Γ =Γ + δΓ, the approach to the steady state is governed by ∂ t δΓ = −{X, δΓ}; i.e., X alone governs the damping dynamics towards that steady state. The matrix Y describes fluctuations, which come along with dissipation in a probability preserving (∂ t Tr(ρ(t)) = 0) quantum mechanical evolution. The steady stateΓ depends on both X and Y . Finally, we remark that the correlation matrix is related to the density operator equation (7) by We may compare this to a Gaussian Hamiltonian equilibrium setting: Introducing the first-quantized, real and antisymmetric Hamiltonian matrix h in the Majorana basis via H = i 4 i,j h ij c i c j , we have at an arbitrary temperature T = 1/β which reduces to Γ (eq) = i sgn(ih) at T = 0.

Purity and purity gap
The purity of a Gaussian state defined by a particular correlation matrix Γ can be revealed by examining the spectrum of the Hermitian positive semidefinite matrix (iΓ) 2 , which we refer to as the purity spectrum. Pure Gaussian states are characterized by a "flat" purity spectrum with eigenvalues all equal to 1, whereas mixed Gaussian states exhibit eigenvalues smaller than 1, each zero eigenvalue indicating the existence of a completely mixed subspace. Intuition regarding the purity spectrum can be gained by constructing a fictitious quadratic Hamiltonian H Γ from the Hermitian matrix iΓ, namely, where the c i are the 2N Majorana basis operators introduced in the previous section. Since Γ is a real antisymmetric matrix, the spectrum of this Hamiltonian consists of real eigenvalues ± n (n = 1, 2, . . . , N ). Most importantly, the positive part of the spectrum of H Γ is the purity spectrum of the Gaussian state represented by the correlation matrix Γ (up to a square root). Exploiting this analogy further, we will identify the eigenvectors of (iΓ) 2 as "eigenmodes" or "modes" (of the fictitious Hamiltonian H Γ ). In particular, we will refer to modes of H Γ associated with zero eigenvalues as zero-purity modes and to the spectral gap of the latter as the purity gap. Such modes are defined in the mode space M consisting of operators of the form γ = v T c with v ∈ R 2N . Modes corresponding to a unit vector v will be referred to as "Majorana" modes since they satisfy γ † = γ and γ 2 = 1. Moreover, we will distinguish two types of zero-purity modes: intrinsic ones, which are determined by the dissipative dynamics, and extrinsic ones, which result from mixed initial conditions (and thus disappear when starting from pure initial states). From a topological perspective, Majorana zero-purity modes that have a topological origin will be referred to as genuine, as opposed to spurious ones. The purity of the steady state is determined by the dissipative dynamics and, if the steady state is not unique, by the purity of the initial state (i.e., by the initial conditions). In the case of interest in this work where the dissipative dynamics is quadratic, one can show (we refer to our previous work [31] for an explicit proof) that there exist initial conditions leading to a pure steady state whenever the corresponding Lindblad operators L i form a set of anticommuting operators, i.e., whenever {L i , L j } = 0 for all i, j * . In the matrix representation defined in the previous section, one can then establish a one-to-one correspondence between the matrices X and Y encoding the dynamics. Intuitively, this can be understood by examining the steady-state equation {X, Γ} = Y (where Γ now denotes the steady-state correlation matrix): if the steady state is pure, the spectrum of the associated correlation matrix Γ (i.e., the purity spectrum) essentially contains no information, since all of its eigenvalues are equal to 1. The information contained in X must therefore be exactly the same as that contained in Y , otherwise the steady-state equation would not be satisfied. In other words, X and Y both contain full information about the dissipative dynamics when the steady state is pure. In that case, one can construct yet another useful object encoding all information about the dynamics: the so-called parent Hamiltonian H parent naturally associated with the Hermitian matrix iY , defined as Clearly, the spectrum of H parent is directly related to that of Y and therefore to that of X as well for dissipative dynamics whose steady state is pure. Remembering the definition of the matrix Y in terms of the Lindblad operators, one can rewrite the parent Hamiltonian in the equivalent form This shows that pure steady states |ψ , which are "dark states" satisfying L i |ψ = 0 for all i (see equation (2)), can equivalently be seen as ground states of H parent . As opposed to the purely fictitious Hamiltonian H Γ that we have constructed above to quantitatively assess the purity of an arbitrary Gaussian state, independently of any dynamics, the parent Hamiltonian therefore describes features associated with the actual (dissipative) dynamics of the system-as expected from its definition from the matrix Y . In fact, we will argue in the next section that the spectrum of H parent encodes all information regarding pure steady states. We emphasize, however, that the parent Hamiltonian does not play such a prominent role in the more general case where the steady state of the dissipative dynamics is mixed (even when starting from pure initial states). As demonstrated in our previous work [31], the above discussion can be formalized and summarized as the following equivalent statements: (i) The steady state is pure (at least for pure initial states); (ii) {L i , L j } = 0 for all i, j (i.e., the Lindblad operators form a set of anticommuting operators); (iii) [X, Y ] = 0 and X 2 = − 1 4 Y 2 (in particular, the spectra of X and Y are directly related); (iv) The dissipative dynamics can be fully described using the parent Hamiltonian H parent = i L † i L i . This last point will be clarified in the next section.

Dissipative gap and Majorana zero-damping modes
In the case where the dissipative dynamics is quadratic, the dynamical approach to the steady state is governed by the associated matrix X, as mentioned in the previous section. This matrix, which is by construction real, symmetric and positive semidefinite, with eigenvalues κ j ≥ 0 and associated eigenvectors v j ∈ R 2N . The eigenvectors of X define "modes" in the mode space defined in the previous section. Assuming that they are normalized to unity, each eigenvector v j can be identified with a corresponding Majorana mode γ j = v T j c. Physically, the eigenvalues of X then correspond to particular damping rates associated with particular Majorana modes. Accordingly, we will refer to the spectrum of X as the damping spectrum and to Majorana modes γ j associated with a vanishing damping rate as Majorana zero-damping modes. One can show (for an explicit proof, see our previous work [31]) that a Majorana mode γ = v T c is a zero-damping mode whenever the following equivalent conditions are satisfied: . We remark, however, that Majorana zero-damping modes satisfying the above conditions need not be spatially localized. The topological nature of the system will play a crucial role in ensuring such localization, as we will demonstrate in section 5.2 below. In general, we will distinguish genuine Majorana zero-damping modes that have a topological origin from spurious ones which a mere artefacts of the dissipative dynamics. A simple example of spurious modes is provided by a lattice site on which no dissipative dynamics takes place. This gives rise to two Majorana zero-damping modes decoupled from the dynamics, which obviously do not have a topological origin.
The damping spectrum describes the dynamical separation (i.e., in time) between particular modes in a similar way as the spectrum of a Hamiltonian determines the energy separation between specific modes. Pushing the analogy further, one can see that the existence of a dissipative gap (or "damping gap") in the damping spectrum leads to the dynamical isolation of bulk and edge modes (through the quantum Zeno effect [57]), thereby providing a dissipative counterpart to the gap protection arising in the Hamiltonian context. Majorana zero-damping modes form a so-called decoherencefree subspace [58] unaffected by the dissipative dynamics and, in the presence of a finite dissipative gap, completely isolated from the rest of the system. The dissipative counterpart of a topological Hamiltonian ground-state degeneracy is then provided by the existence of a non-local decoherence-free subspace associated with zero-damping Majorana modes.
We finally clarify the role of the parent Hamiltonian defined in the previous section in light of the considerations above. When the steady state of the dissipative dynamics is pure, the spectrum of X (i.e., the damping spectrum) directly maps to the spectrum of Y , which in turn trivially maps to that of H parent . The parent Hamiltonian thus contains all information about the dissipative dynamics in that case. When the steady state is mixed (independently of the initial state), however, the tight relationship between X and Y (or H parent ) breaks down and the parent Hamiltonian becomes insufficient to describe the dynamics. In this more general case, one can show (we refer to our previous work [31]) that a zero mode of H parent (or, equivalently, of the matrix Y ) does not necessarily correspond to a zero mode of X (i.e., to a zero-damping mode), although the converse is always true. In other words, zero modes of the parent Hamiltonian need not be Majorana zero-damping modes of the dissipative dynamics. In fact, any zero mode of H parent which does not coincide with a zero mode of X gives rise, in steady state, to an intrinsic zero-purity mode. This crucial phenomenology will be key to understanding the non-equilibrium topological effects that will be exemplified below.
To conclude this section, we remark that Majorana zero-damping modes do not benefit from the protection mechanism featured by H parent due to the antisymmetry of the matrix Y . While the antisymmetry of Y forces H parent to have an even number of Majorana zero modes, such that spatially isolated modes cannot be affected by local perturbations, the fact that X is symmetric does not lead to such restriction. Although this is potentially harmful for Majorana zero-damping modes, we will demonstrate in the sections below that this can also lead to intriguing physics with no Hamiltonian ground state counterpart.

Topological properties of the bulk
In this section, we focus on the topological properties of the bulk of driven-dissipative fermionic systems with Gaussian steady states. In particular, we identify the correlation matrix, which fully describes such states, as a fictitious first-quantized Hamiltonian and use the latter to construct a topological classification in complete analogy to the conventional Hamiltonian scenario.
The topological classification of gapped states of non-interacting fermions can be achieved on the basis of symmetry properties of the corresponding Hamiltonian under time-reversal, charge-conjugation (or particle-hole) and chiral (or sublattice) transformations, as was proposed by Schnyder et al [59] following the classification of random matrices developed by Altland and Zirnbauer [54]. Ten symmetry classes were proved to be sufficient for an exhaustive classification of topological phases in any spatial dimension, and an alternative approach was later introduced by Kitaev in the powerful framework of topological K-theory [60][61][62][63] (see reference [64], e.g., for a self-contained review). We argue below that all classification schemes developed in the Hamiltonian setting can be automatically applied to classify the Gaussian steady states of a dissipative dynamics. We do not provide an exhaustive classification, however, since this can be done straightforwardly based on the references cited above. Instead, we construct an explicit topological classification for two symmetry classes of particular interest for this work, namely, for dissipative systems belonging to the symmetry classes D and BDI.

Steady-state symmetries
We first study how symmetries of the Lindblad operators translate into symmetries of the correlation matrix. To this end, we consider a Gaussian dissipative dynamics with unique steady state, i.e., the corresponding correlation matrix Γ is a unique solution of with matrices X and Y defined as in section 2.5 (see equation (10), in particular). We assume that the Lindblad operators L i are invariant, up to a phase factor, under some symmetry group G: where g ∈ G. In order to preserve the linearity of the Lindblad operators in the fermionic operators, G must act linearly on the 2N Majorana operators c j introduced in section 2.5 above. These operators form an orthonormal basis of the mode space . Clearly, any symmetry g ∈ G (unitary or antiunitary) must act linearly on M and transform the Majorana basis defined by the operators c j into another Majorana basis. In mode space, a symmetry g ∈ G must therefore be represented by a real orthogonal matrix S g ∈ O(2N ), Note that this formula allows to analyze the symmetry properties of the state also in the more general case where the dynamics is governed by both a Liouvillian and a Hamiltonian, see equation (24) below. The Lindblad operators are defined in the Nambu space N ∼ = C 2N of operators L = l T c with l ∈ C 2N , which can be viewed as a complexification of M. In Nambu space, the relevant representation of g ∈ G is given by S g if g is unitary and S g K (= KS g ) if g is antiunitary. Here K is the complex conjugation operator defined such that Ki = −iK. The Lindblad operators L i = l T i c are therefore invariant (up to a phase factor) under the symmetry g ∈ G if and only if where φ i ∈ [0, 2π) (note that the phase factors e iφ i do not affect the form of the Liouvillian). Using equation (10), we then find that the matrices X and Y encoding the dissipative dynamics have the properties where the positive and negative signs corresponds to the cases where g is unitary or antiunitary, respectively. The steady-state equation can then be written as and, since we have assumed that the steady state is unique, we obtain This shows that the symmetries of the Lindblad operators are naturally reflected in symmetries of the steady-state correlation matrix Γ * . In turn, this implies that the matrix Γ can be used to construct a symmetry-based topological classification of Gaussian steady states. Of particular importance for this purpose are the two discrete symmetries corresponding to particle-hole (PHS) and time-reversal (TRS) symmetry, respectively. The former corresponds to the "+" sign in equation (24) and is trivially satisfied as S g = 1 in our case, as we argue below, while the latter corresponds to the "-" sign and depends on the specific form of the Lindblad operators. We remark that chiral symmetry (defined as the combination of PHS and TRS [33]) is automatically satisfied whenever TRS is present, since the system always has PHS, by construction. In that case, there exists matrices S C and S T corresponding to PHS and TRS, respectively, such that equation (24) is satisfied (with a "+" sign for S C and a "-" sign for S T ), and one can easily verify that the combination of PHS and TRS (i.e., chiral symmetry) leads to

Topological classification and topological invariants
Let us consider an arbitrary Gaussian steady state ρ, fully characterized by its Using that the matrix Γ is real and antisymmetric, we construct a corresponding fictitious free-fermion Hamiltonian as in section 2.6 where the purity spectrum was defined (not to be confused with the parent Hamiltonian introduced in the same section), thereby establishing a one-to-one correspondence between the set of Gaussian steady states and the set of free-fermion Hamiltonians. It is clear that the symmetries of the dissipative system-embedded in Γ-are the same as that of the Hamiltonian system defined by H Γ , since Γ can be viewed as the "first-quantized" Hamiltonian corresponding to the "second-quantized" Hamiltonian H Γ . The problem of classifying steady states according to topological properties is therefore equivalent to that of classifying Hamiltonian systems of non-interacting fermions, which is the main message of this section. Consequently, both the symmetry-based classification of references [54,59] and the K-theory approach of references [60][61][62] can be directly applied in the dissipative framework. Note that the Hamiltonian H Γ always takes a Bogoliubov-de Gennes form when expressed in terms of the original fermionic operators a i , a † i , and is therefore automatically particle-hole symmetric.
The topological classification crucially relies on the existence of a bulk spectral gap and is essentially based on the mathematical concept of homotopy equivalence, or equivalence under continuous deformations * . More specifically, two gapped Hamiltonians H Γ and H Γ are considered as "topologically equivalent" if they can be continuously deformed into each other without closing the gap. The dissipative counterpart of this equivalence is provided by the mapping of equation (25): two Gaussian steady states corresponding to correlation matrices Γ and Γ will be considered as topologically equivalent and referred to as belonging to the same topological phase if and only if they can be continuously deformed into each other without closing the bulk purity gap †. The existence of a bulk purity gap is therefore the key ingredient required to define topological order in the dissipative setting.
The mapping defined by equation (25) and the results of references [59][60][61][62] provide us, in principle, with a general topological classification of all possible (purity) gapped Gaussian steady states according to symmetries and to the spatial dimension of the system. We now sketch this construction focusing on the symmetry classes that are most relevant for the dissipative systems considered in this work.
An arbitrary 2N × 2N steady-state correlation matrix Γ (N being the number of fermionic modes, or the number of sites for systems of spinless fermions defined on a lattice) can be brought to a block diagonal form where Q is an orthogonal matrix and ± n are the real eigenvalues forming the spectrum of the Hermitian matrix iΓ. The purity spectrum is defined by the spectrum of the real symmetric matrix (iΓ) 2 , which is doubly degenerate with positive eigenvalues 2 n . Assuming that it is gapped, such that | n | > 0 for all n, the matrix Γ can be continuously deformed into a topologically equivalent matrixΓ with a "flat" purity spectrum Since (iΓ) 2 = 1, this "spectrally flattened" correlation matrix defines a pure Gaussian state which is topologically equivalent to the not necessarily pure original steady state of interest. The matrixΓ allows us to construct an orthogonal spectral projection operator P (see, e.g., reference [13]) defined as which projects onto the N -dimensional subspace associated with eigenvectors of iΓ with negative eigenvalues ‡. This operator plays a crucial role in the topological * In the framework of K-theory, stable equivalence also plays a crucial role in the definition of "topological equivalence" (see reference [60]). †Note that the spectrum of H Γ is defined by that of the Hermitian matrix iΓ, and is therefore in one-to-one correspondence with the spectrum of (iΓ) 2 , i.e., with the purity spectrum.
‡From the point of view of the Hamiltonian H Γ (see equation (25)), P projects onto the subspace of eigenstates of H Γ with negative energy.
classification of Gaussian steady states as well as in the construction of associated topological invariants, as we will demonstrate below.
We first consider the case of spinless fermions on a d-dimensional lattice with periodic boundary conditions evolving under a translation-invariant dissipative dynamics. It will be convenient to label the local Majorana operators c i as c i,λ , where i refers to a particular lattice site at position r i and λ = 1, 2 distinguishes the two local Majorana operators associated with the corresponding fermionic creation and annihilation operators a i and a † i , i.e., In momentum space, the steady-state correlation matrix Γ then takes the form of a 2 × 2 complex antihermitian matrix Γ(k) with components (λ, µ = 1, 2) where N is the total number of lattice sites and c k,λ = 1 Since the matrix Γ(k) satisfies the condition Γ(k) * = Γ(−k), the spectrum of the Hermitian matrix iΓ(k) is given by eigenvalues ± (k) with (k) = (−k) > 0. This allows us to introduce the "spectrally flattened" momentum-space correlation matrixΓ(k) (the eigenvalues of iΓ(k) being ±1) and the associated spectral projection operator projecting onto the subspace C 2 − (k) of the 2-dimensional complex vector space C 2 spanned by the (complex) eigenvectors of iΓ(k) associated with the negative eigenvalue −1. The operators P (k) form a family of operators labeled by the wavevectors k belonging to the Brillouin zone, i.e., to the d-dimensional torus T d = [−π, π] d . They define a fiber bundle ⊕ k∈T d C 2 − (k) on the Brillouin zone manifold, with fibers C 2 − (k) assigned to each point k ∈ T d . The problem of classifying Gaussian steady states according to topology therefore reduces to that of classifying fiber bundles over a torus T d (the Brillouin zone). In the present case, the fibers are vector spaces and the fiber bundle is thus a vector bundle, for which a complete topological classification can be constructed using the approach of reference [59] or K-theory [60][61][62]. We present below the results pertaining to the two types of systems that are most relevant to this work (see section 4), namely: 2D dissipative systems without TRS (symmetry class D of Altland and Zirnbauer) and 1D dissipative systems with TRS (symmetry class BDI).
The topological classification generally reduces to the identification of all possible homotopy classes of continuous maps k → P (k) from the Brillouin zone manifold to the space of spectral projection operators. Since the 2×2 complex matrix iΓ(k) is Hermitian and unitary, with (iΓ(k)) 2 = 1, the spectral projection operator can be expressed as where n(k) is a unit vector on the 2-sphere (S 2 ) and σ is a vector of Pauli matrices. In the absence of additional symmetries, the spectral projection operator thus describes a mapping from the Brillouin zone torus to the sphere * .
In 2D dissipative systems where the steady state belongs to the symmetry class D, the only symmetry thatΓ(k) (or P (k)) has is the PHS automatically embedded in Gaussian steady states. The homotopy classes of maps k → n(k) then form a group which can be identified with the homotopy group π 2 (S 2 ) = Z † (which is non-trivial) and one can have non-trivial topological steady states and distinguish them via an integer topological invariant known as the Chern number, defined as where in the last line we have used the fact thatΓ(k) = i(n(k) · σ) (see equation (31)).
In 1D dissipative systems where PHS is the only symmetry, the Brillouin zone corresponds to a circle S 1 and the homotopy classes of maps k → n(k) form a homotopy group π 1 (S 2 ) which is trivial. This means that the vector n(k) cannot "twist" in a way that cannot be continuously untwisted as we move along the Brillouin zone circle, such that all steady states are necessarily topologically trivial. The situation can change in the presence of additional symmetries, however, since additional constraints on P (k) (or, equivalently, onΓ(k) or n(k)) can potentially restrict the set of allowed continuous deformations. If the 1D dissipative system has TRS, then it also has chiral symmetry and, hence, belongs to the symmetry class BDI. In that case, as argued below equation (24), one can find a unitary matrix Σ such that Σ 2 = 1 and for all k. (In our case, TRS takes the form The matrix Σ has eigenvalues ±1 and can be expressed as Σ = a·σ, where a is a real unit vector that does not depend on k. The above condition then takes the form a · n(k) = 0 for all k. In other words, chiral symmetry restricts the vector n(k) from the sphere to a circle in the plane perpendicular to a. As a result, the relevant homotopy group becomes non-trivial, given by π 1 (S 1 ) = Z. The corresponding topological invariant can be constructed by choosing a basis where with θ(k) ∈ [0, 2π) and e iθ(k) = n y (k) + in x (k). The U (1) phase θ(k) therefore fully characterizesΓ(k) and contains all required information to construct a topological * P (k) can alternatively be seen as an element of U (2)/U (1) × U (1) G(1, 2), G(m, n) being the set of all n-dimensional subspaces in C m (also known as a complex Grassmann manifold) [59]. †π 2 (S 2 ) formally classifies maps from S 2 to S 2 , but since π 1 (S 2 ) is trivial, it equivalently classifies maps from T 2 to S 2 .
invariant distinguishing between different topological classes of steady states. The relevant topological invariant, the winding number, takes the form where we have introduced projection operators P ± = (1 ± Σ)/2 defined so that We remark that the above construction of topological classes and topological invariants crucially relies on translational symmetry. However, one expects bulk topological properties to be robust against weak dissipative perturbations that preserve the purity gap as well as the symmetries of the system. This can be rigorously demonstrated in the K-theoretic framework of reference [60] and in the more general framework of noncommutative geometry developed in references [13,66], allowing to extend the above construction to disordered (and possibly finite) systems. Although the details are beyond the scope of this work, we emphasize that the existence of quantized (integer) topological invariants in such systems crucially relies on the quasi-local nature of the spectral projection operator. Here the quasi-local nature of the corresponding steadystate correlation matrix is ensured in the presence of a dissipative gap, which is a spectral property of the model in contrast to the purity gap, which relates to mode occupations. More specifically, the existence of a finite dissipative gap ∆ implies an exponentially fast relaxation to the steady state (on a characteristic time scale τ ∼ 1/∆) and leads to the exponential decay of all spatial correlations, such that the spectral projection operator (see equation (28)) satisfies for some constants c, α > 0 * . Note that the converse is also true [50]. The operator P ij defined by (P ij ) λµ = P iλ,jµ (see reference [13]) is thus quasi-diagonal, i.e., there exists some constants c , c > 0 and α > d (d being the spatial dimension of the system) such that where ||.|| HS is the Hilbert-Schmidt norm and ||.|| the usual operator norm. This "localization" condition satisfied by the spectral projection operator in the presence of a finite dissipative gap crucially allows for the definition of the Chern number and winding number topological invariants in disordered finite systems. Such a construction can be found in the appendix C of reference [13], for example. * Note that the purity gap of Γ ensures that the exponential decay of all spatial correlations remains under the continuous "spectral flattening" transformation Γ →Γ.

Phenomenology of non-equilibrium topological phase transitions
The above discussion shows that the existence of two gaps is necessary to identify the topological properties of the bulk (Gaussian) steady state: (i) the purity gap, which allows to identify a particular topological class (e.g., Z) to which the steady state belongs, and (ii) the dissipative gap, which ensures that the steady-state correlations are quasi-local (as defined above), so that in particular the Chern and winding number topological invariants defined above (see equations (32) and (35)) remain quantized in the presence of disorder. Since the dissipative and purity spectra are in general independent quantities-which can be seen from the fact that two independent matrices X and Y are required to describe the dissipative dynamics-these two gaps can close independently of each other. However, the closure of the dissipative gap alone gives rise to a critical behavior in steady state. As a consequence, non-equilibrium topological phase transitions can occur in three distinct ways: (i) Via the closure of the bulk dissipative gap only (criticality); (ii) Via the closure of the bulk purity gap only (no criticality); (iii) Via the closure of the both the bulk dissipative and purity gaps (criticality).
These three possibilities will be exemplified in explicit models in section 7 below. Clearly, the same phenomenology must appear at interfaces between distinct non-equilibrium topological phases. We will investigate such physics in detail in section 5 below.

Physical constraints
Our discussion so far was very general, with the only assumption that the steady state of the dissipative dynamics is a Gaussian fermionic state. We now discuss restrictions that arise in "typical" experimental realizations, as sketched in section 2.3 above. More specifically, we consider quasi-local dissipative processes in one-and two-dimensional lattice systems; assuming that the corresponding Lindblad operators act on a quasilocal subset of sites and have a spatially isotropic (or s-wave symmetric) creation part resulting from spontaneous emission processes. We additionally assume that the system is translation-invariant and that the steady state of the dissipative dynamics is unique and pure. Under this purity assumption, the steady state of the system can be viewed as the ground state of the parent Hamiltonian associated with the dissipative dynamics (see equation (15) and the discussion of section 2.6), and the existence of topological order can be assessed in the exact same way as in the Hamiltonian setting. In particular, as for the topological classification of section 3.2 above, one can rely on the general Hamiltonian classification of gapped topological phases of non-interacting fermions (according to symmetry and spatial dimension) developed in references [33,54,59,67]. Although we assume a pure steady state, we remark that the conclusions drawn below will also hold in the more general case where the steady state is not necessarily pure but has a gapped purity spectrum. We have argued in section 3.2 that an arbitrary state exhibiting a finite purity gap is topologically equivalent to a pure state exhibiting a completely "flat" purity spectrum. To identify the possible classes of topological steady states that can be reached when taking account the "typical" physical constraints described above, one can therefore focus exclusively on pure states, without loss of generality. We emphasize that in that case the (automatically "flat") correlation matrix iΓ describing the pure steady state precisely corresponds to the parent Hamiltonian in its "first-quantized" form, namely, H parent = i i,jΓ ij c i c j = HΓ, where HΓ is the fictitious Hamiltonian that can always be constructed from the (Gaussian) steady state (see equations (25) and (28) and discussion thereof).
Translational symmetry makes it most convenient to express the Lindblad operators in the momentum-space form L k = u k a k + v k a † −k . In order for our assumption of a pure and unique steady state to be satisfied, the Lindblad operators must form a complete set of anticommuting operators (i.e., {L k , L k } = 0 for all k, k ), which implies that Translational symmetry additionally restricts the possible lattices for the system to Bravais lattices. In that case, the set of Lindblad operators is complete whenever there are as many Lindblad operators as lattice sites. As argued above, the pure steady state can equivalently be described as the ground state of the parent Hamiltonian H parent = k L † k L k (see section 2.6), which takes here the explicit form Dropping the constant term k |v k | 2 , it becomes where we have defined Since the Lindblad operators are generally defined up to a phase, we can assume v k to be real, without loss of generality. Using equation (39), we thus find Our assumption that the creation part of the Lindblad operators is isotropic implies that v k = v −k * . Since equation (39) must be satisfied, we conclude that the condition u k = −u −k must be imposed. In practice, this can easily be achieved by tuning the phases of the driving lasers, as argued in our previous work [31], and will therefore be assumed in the following. We then have

Symmetry classes
Equations (41), (43), and (44) show that the parent Hamiltonian takes the generic form of a Bogoliubov-de Gennes (BdG) Hamiltonian and therefore exhibits, by construction, particle-hole (or charge-conjugation) symmetry (PHS). In the scheme of reference [33], further classification can be achieved by considering time-reversal symmetry (TRS). In the Nambu representation, the parent Hamiltonian can be written as where H k is the "first-quantized" parent Hamiltonian and Ψ † k = (a † k , a −k ) the Nambu field operator. Particle-hole and time-reversal symmetries are then present whenever there exists 2 × 2 unitary matrices U C and U T , respectively, such that the following conditions on the Hamiltonian are satisfied: The symmetry class of the system is determined by whether or not such matrices can be identified and, if so, on the sign of U * C U C and U * T U T (which can only take values ±1, as argued in reference [33]). In the present case, the parent Hamiltonian describes a spin-polarized superfluid and exhibits-again, by construction-PHS with U C = σ x , such that U * C U C = +1. It is clear from the form of H parent (see equations (41), (43), and (44)) that the system additionally exhibits TRS if and only if ∆ * −k = ∆ k which, remembering equation (42), is satisfied provided that ∆ k is real up to a global phase * . In that case, equation (47) holds for U T = 1 (such that U * T U T = +1) and the system exhibits chiral symmetry.
In summary, two symmetry classes of steady states can be realized through "typical" quantum-reservoir engineering (see section 2.3) in driven-dissipative systems of spinpolarized fermions: if the Lindblad operators break TRS, the steady state belongs to the symmetry class D of Altland and Zirnbauer; otherwise it belongs to the symmetry class BDI. According to the topological classification of reference [33], one must therefore consider either (i) 2D systems in the symmetry class D (with broken TRS) or (ii) 1D systems in the symmetry class BDI (with TRS) in order to have the (a priori) possibility of reaching topologically non-trivial steady states. In these systems, distinct topological states are characterized by distinct values of an integer-valued (Z) topological invariant: * One must keep in mind that the gap function ∆ k is defined up to a global phase, since H parent is invariant under the global U (1) gauge transformation a † k → e iθ a † k .
the Chern and winding number, respectively. Most importantly, the zero-energy edge modes appearing in such systems in the Hamiltonian setting are Majorana fermions. This opens up the possibility to generate spatially isolated Majorana zero-damping modes in the dissipative setting of interest in this work.

Chern number
Let us now examine more closely the case of a 2D translation-invariant dissipative system with broken TRS, with a pure and unique steady state that accordingly belongs to the symmetry class D of Altland and Zirnbauer and is characterized by the Chern number topological invariant of equation (32). We show below that quasi-local Lindblad operators generating such dissipative dynamics do not allow to obtain phases with a nonzero Chern number, in contrast to the Hamiltonian setting. Although this result seems to be a no-go statement for Majorana zero modes in such systems, we will demonstrate in section 8 that unpaired Majorana zero-damping modes can exist in dissipative systems with vanishing Chern number provided that the steady state is not pure but mixed.
As we have discussed above, a pure steady state can equivalently be viewed as the ground state of the parent Hamiltonian H parent associated with the dissipative dynamics, given by equation (45) (see also equations (41), (43) and (44)), and its bulk topological properties can be characterized by means of the spectral projector P k = 1 2 (1 + n k · σ), where n k · σ = H k /N k (N k being a normalization factor defined such that ||n k || = 1). The corresponding Chern number is determined by the vector n k (see equation (32)), which takes the explicit form where N k = ξ 2 k + |∆ k | 2 = |u k | 2 + |v k | 2 * . Note that N k corresponds to the spectrum of H parent and thus defines the dissipative gap (see section 2.7). In order for the vector n k to be well defined, the functions u k and v k must not vanish simultaneously anywhere in the Brillouin zone. Introducing ϕ k = |ϕ k |e iθ k = v k /u k and a "Fermi surface" F = {k : |ϕ k | > 1/ } (defined for some arbitrary > 0) which generally † takes the form of a finite union of piecewise smooth, simple closed curves F ,λ (i.e., F = λ F ,λ ), the Chern number can be expressed (see reference [31]) as a sum ) * We assume that the system is infinite such that the components of n k are continuous functions of k in the Brillouin zone, in which case the Chern number given by equation (32) is a well-defined quantity.
†More pathological Fermi surfaces can in principle be encountered, but the Chern number is a mere definition in that case. of winding numbers defined by Choosing 1, the above expressions tell us that the value of the Chern number can be inferred from the behavior of ϕ k = v k /u k around the zeroes of u k -which do not coincide with zeroes of v k since N k must be non-zero. Since the norm of ϕ k diverges upon approaching such points, it is the phase winding of ϕ k around these points that encodes all information about the Chern number. More specifically, each of the zeroes of u k contributes to the Chern number by an integer value corresponding to the winding number of the phase θ k around the latter (in a conventional, e.g., counter-clockwise direction). Since the function v k can always be chosen as real (see discussion below equation (43)), θ k corresponds to the phase of u k (up to a sign that can be absorbed into the definition of the Chern number by a change of convention ν 2D → −ν 2D ) and the Chern number simply reduces to the sum of the phase windings of the function u k around its zeroes.
In the Hamiltonian setting, the phase of u k is defined by the phase of the gap function ∆ k , whereas the norm of u k is defined by the norm of ∆ k and the value of the independent quantity ξ k (the single-particle dispersion) [12], such that the phase and the modulus of u k are essentially independent. In our dissipative setting, however, u k is a single complex function defined on the Brillouin zone and, therefore, the phase and the norm of u k are closely related. This crucially restricts the possible values of the Chern number: the function u k defines a smooth vector field (Re(u k ), Im(u k )) on a torus and, according to the Poincaré-Hopf theorem, the sum of its phase windings around its zeroes must vanish. We thus conclude that the only value of the Chern number that can be achieved in the above dissipative setting is zero.
Note that the vanishing of the Chern number is a direct consequence of the smoothness of the vector field corresponding to u k , which in turn is due to the quasi-local nature of the Lindblad operators. The function u k defines a smooth (C ∞ ) vector field if and only if the norm of the coefficients u j−i = u(r j − R i ) encoding the annihilation part of the Lindblad operators (see equation (52)) decays faster than any power of 1/|r j −R i |. A power-law decay of the Lindblad operators would therefore be necessary to be able to engineer phases characterized by a non-zero Chern number in our dissipative setting.

Dissipative edge physics
Thus far we have concentrated on bulk properties of purely dissipative systems described by a quadratic Lindblad master equation. It is clear, however, that the existence of nontrivial topological bulk properties must somehow be reflected in the physics occurring at physical edges or in topological defects-which we both refer to as "edges" in the following. In the Hamiltonian setting, this intimate connection has been rigorously demonstrated and formalized-at least, for non-interacting systems-as a bulk-edge correspondence (or bulk-boundary correspondence) relating the topological nature of the bulk to the number of gapless modes occurring at an edge [13,[32][33][34][35]. Of particular interest here is the bulk-edge correspondence pertaining to (i) 1D systems characterized by a winding number topological invariant ν 1D (see equation (35)) and (ii) 2D systems characterized by a Chern number topological invariant ν 2D (see equation (32)), which tells us that a number m = |ν (1) − ν (2) | of gapless edge modes must be present at the interface between two topological phases characterized by integer topological invariants (winding or Chern numbers, as appropriate) ν (1) and ν (2) , respectively. In general, such modes are robust (e.g., against disorder) as long as the Hamiltonian system remains in the same topological class. Since the purity gap can close in the dissipative setting of interest in this work, it is highly unclear whether similar universal signatures of bulk topological properties can be observed. We thus proceed, in the next sections, to investigate the edge physics that may arise in systems described by a quadratic dissipative dynamics.

Dissipative bulk-edge correspondence
Let us consider a generic dissipative dynamics characterized by a finite dissipative gap and described by Lindblad operators L i or, equivalently, by matrices X and Y as defined in section 2.5. It is clear that the bulk-edge correspondence pertaining to Hamiltonian systems must also be satisfied when the steady state of the dissipative dynamics is pure, since such a state can equivalently be regarded as a ground state of the parent In that case, the spectrum of H parent defines the damping spectrum, and at least m = |ν (1) − ν (2) | Majorana zero-damping modes must be present at the interface between two nonequilibrium topological phases characterized by Chern or winding number topological invariants ν (1) and ν (2) . In the more general situation where the steady state exhibits a finite purity gap but is not pure, the Hamiltonian bulk-edge correspondence obviously still applies to the parent Hamiltonian, such that H parent still exhibits a minimum of m = |ν (1) − ν (2) | Majorana zero modes at the interface, but none of these modes is guaranteed to be Majorana zero-damping modes of the dissipative dynamics. As argued in section 2.7, Majorana zero mode of H parent must either correspond to Majorana zerodamping modes or give rise to intrinsic Majorana zero-purity modes in steady state. The dissipative counterpart of the Hamiltonian bulk-edge correspondence can therefore be stated as follows: As long as the topological nature of the steady state remains the same in the bulk-which is true if the purity gap remains finite and if the symmetries of the Lindblad operators are preserved-a minimum total number of m = |ν (1) − ν (2) | ≡ ∆ν Majorana zero-damping and intrinsic zero-purity modes must be present at the interface between two non-equilibrium topological phases characterized by topological invariants ν (1) and ν (2) , respectively. With obvious notations, this result can be summarized in the form We remark that such a dissipative bulk-edge correspondence is only guaranteed to hold in the presence of a finite dissipative gap, since the existence of such a gap is necessary to ensure the exact quantization of the topological invariants, as discussed in section 3.2.
Most importantly, the inequality appearing in equation (51) originates from the fact that "spurious" Majorana zero modes may arise "accidentally", as opposed to Majorana zero modes which we refer to as "genuine" that have topological origins. A strict equality is therefore only obtained when restricting to genuine Majorana zero-damping and zeropurity modes. Despite its important predictive power, the above result does not provide specific information regarding the exact number of genuine Majorana zero-damping modes that appears at an edge-as we will argue below, such knowledge must come from additional considerations unrelated to the topological nature of the system. One may therefore wonder whether this number can be a robust quantity and, if so, under what conditions. This will be the focus of the next sections: we will first investigate systems whose edge dissipative dynamics emerges in a very natural way by extending the bulk dissipative dynamics as close as possible to a physical edge. Since the topological nature of generic Gaussian states is only well-defined in the presence of a finite purity gap and since all Gaussian states featuring such a gap can be continuously deformed to pure states (see section 3.2), we will focus on dissipative dynamics driving the system into a pure steady state. In that case, the purity of the system obviously forces the number (m purity ) edge of Majorana zero-purity modes to vanish in equation (51), and one can try to construct Majorana zero-damping modes explicitly. Having examined this simple situation where the number (m damping ) edge of genuine Majorana zero-damping modes is solely determined by the topological nature of the bulk, we will then extend our discussion to the more general case where the dissipative dynamics can affect the purity of the system while maintaining a finite purity gap in the bulk, and conclude by clarifying the role of topology in the dissipative setting central to this work.

General form of Majorana zero-damping modes for translation-invariant dissipative processes
We have introduced above a dissipative bulk-edge correspondence which crucially relies on the existence of two gaps: the dissipative and purity gaps. Both these gaps must be finite in the bulk in order for equation (51) to hold and-as argued in section 3.2-for the bulk to have a well-defined topological nature with associated topological invariants that are quantized. Any of these two gaps can in principle close at an edge, which gives rise to the variety of possibilities allowed by equation (51). In the following section, however, we focus on the case where the purity gap remains open at the edge, thus forbidding the appearance of intrinsic Majorana zero-purity modes. In that case, the topological properties of the steady state are equivalent to that of a pure state-as argued in section 3.2-and the bulk-edge correspondence describes the behavior of a single gap at the edge (namely, the dissipative gap) exactly as in the Hamiltonian setting. From a e ? S i ⌘ i L i Figure 2. Left: Translation-invariant dissipative dynamics terminated at the edge in the case of "cross-shaped" Lindblad operators acting on a subset of 5 sites (as in the explicit example of section 8). Lindblad operators acting within the system (in blue) are taken into account into the dissipative dynamics, whereas Lindblad operators requiring truncation at the edge (in red) are not considered. Right: Lindblad operator L i associated with a lattice site i that coincides with its center of symmetry S i , and example of an edge whose orientation is according to a vector e ⊥ (see main text).
topological point of view, all steady-state properties are the same as if the steady state were completely pure. We therefore restrict ourselves, without loss of generality, to pure (Gaussian) steady states. In addition we assume translation invariance in the following in order to make the explicit construction of Majorana zero-damping modes analytically tractable.
We consider a d-dimensional lattice system with a (d−1)-dimensional physical edge and assume that the system evolves under a translation-invariant dissipative dynamics consisting of a periodic repetition, on the lattice, of a single quasi-local dissipative process (or Lindblad operator) everywhere in the bulk and as close as possible to the edge (see figure 2), so that every lattice site i of the bulk becomes associated with a Lindblad operator L i whose form is independent of i. In order to ensure that the steady state is pure we further assume that the Lindblad operators form a set of anticommuting operators. For simplicity (and in order to incorporate the typical physical constraints discussed in section 4), we moreover restrict ourselves to Lindblad operators L i that possess a center of symmetry, which we denote as S i (the index i used to distinguish Lindblad operators then corresponds, by convention, to the index of the lattice site located closest to S i * ). We define a position vector R i corresponding to S i and similarly define the position of the lattice sites j as r j . In the translation-invariant setting defined above, the Lindblad operators then take the generic form where I(i) denotes the subset of sites j onto which L i acts in a non-trivial way (see figure 2), and u j−i , v j−i , and a j are shorthand notations for u(r j − R i ), v(r j − R i ), and a(r j ), respectively * . As argued in section 2.7 above, a necessary and sufficient condition for a Majorana mode γ to correspond to a Majorana zero-damping mode of the dissipative dynamics is {L i , γ} = 0 (for all i). Using this condition as well as the translation-invariant form of the Lindblad operators, one can show (see Appendix A.3) that any Majorana zero-damping mode γ must take the generic form m n e n ) + h.c. , (53) where N > 0 is a normalization factor, φ i ∈ [0, 2π) a phase, {e n } n=1,2,...,d a set of primitive vectors associated with the d-dimensional Bravais lattice on which the system is defined, and {m 1 , m 2 , . . . , m d } a set of integers defined such that the vectors r i + d n=1 m n e n span the positions of all sites in the system. Most importantly, {β en } n=1,2,...,d is a set of real factors determining the increase of the "spatial wave function" corresponding to γ in each of the directions defined by the primitive vectors e n . Note that β j−i = (β i−j ) −1 is satisfied owing to the translation invariance discussed above (β j−i being a shorthand notation for β(r j − r i )) †.
Equation (53) states that γ-or, more precisely, its spatial "wave function"-is completely delocalized when |β en | = 1 for all n. In the less pathological case where |β en | = 1 for some n, γ grows exponentially in the direction σ n e n with σ n = ±1 defined such that |β σnen | > 1 (or, equivalently, such that |β −σnen | = |β σnen | −1 < 1). In that case, γ can be normalized if and only if the system is finite in the direction of "steepest" exponential increase, i.e., as long as the edge of the system is perpendicular to the direction defined by the unit vector e ⊥ ∝ d n=1 log (|β σnen |)σ n e n , as shown in figure 2. The corresponding Majorana zero-damping mode γ is then exponentially localized along that edge, decaying in every direction −σ n e n away from the edge (into the bulk) on a characteristic length scale ξ n = 1/ log (|β σnen |).
In general, the existence of Majorana zero-damping modes γ of the form (53) can be assessed from the knowledge of the Lindblad operators by looking for solutions of the equation which can be expressed solely in terms of {β en } n=1,2,...,d using the relations (i) β k−i = β k−j β j−i (for any triple of indices (i, j, k)); where u j−i and v j−i are fixed coefficients defining the Lindblad operators (see equation (52)) and (j, j) denotes pairs of sites located symmetrically around the center of symmetry S i of the Lindblad operators. If there exists a set {β en } n=1,2,...,d of real factors satisfying equation (54) for some value of the phase φ i , the system supports at least a single Majorana zero-damping modes of the form (53). We remark that the phase φ i has a direct physical meaning in 2D systems where time-reversal symmetry is broken (symmetry class D), where it is simply defined by the orientation of the edge. We refer to Appendix A.3 for further details as well as for an explicit derivation of the above results.

Robustness of Majorana zero-damping modes and role of topology
We now investigate the robustness of Majorana zero-damping modes against weak quasi-local dissipative perturbations in the general case where the steady state is not necessarily pure but is characterized by a finite bulk purity gap, so that its topological nature is well-defined. We will assume that the dissipative dynamics giving rise to such modes consists of n independent dissipative processes (described by Lindblad operators L i (i = 1, 2, . . . , n)) and will distinguish two types of perturbations: (i) weak quasilocal perturbations affecting each of these n dissipative processes and (ii) perturbations that introduce additional weak quasi-local dissipative processes (or Lindblad operators). Distinguishing spurious from genuine Majorana zero-damping modes as in section 2.7, we will demonstrate that topology does not guarantee the robustness of genuine Majorana zero-damping modes against perturbations, although it is crucial to be able to isolate such modes spatially. Robustness may additionally be ensured by geometrical properties of the system imposed through quantum reservoir engineering. Majorana modes are mathematically defined in the mode space M ∼ = R 2N consisting of real linear combinations of local Majorana basis operators c j (j = 1, 2, . . . , 2N ), where N is the number of lattice sites in the system. In order for a particular Majorana mode γ = γ † to be an exact zero-damping mode of the dissipative dynamics, one must have {L i , γ} = {L † i , γ} = 0 for all i (see section 2.7), which translates, in mode space, as two independent conditions {L i + L † i , γ} = 0 and {L i − L † i , γ} = 0 for each i. Every Lindblad operator therefore imposes two independent constraints for Majorana zerodamping modes in R 2N as long as it is neither Hermitian nor antihermitian, which can safely be assumed, e.g., in the presence of disorder. Consequently, the dynamics generated by n < N dissipative processes necessarily gives rise to 2(N − n) exact Majorana zero-damping modes, independently of the form-quasi-local or not, with or without specific symmetries-of the corresponding Lindblad operators, as illustrated in figure 2. We emphasize, however, that such modes which trivially do not take part in the dynamics need not be spatially localized and isolated from each other (as is required for them to exhibit non-Abelian exchange statistics; see section 6). As embedded in the dissipative bulk-edge correspondence, it is the topological nature of the bulk that plays a crucial role in isolating Majorana zero-damping modes away from each other, allowing to "fractionalize" the fermionic degrees of freedom of the system. Robust isolated Majorana zero-damping modes may therefore arise from the interplay between topology and the "incomplete" nature of the dissipative dynamics (i.e., n < N ), as we will demonstrate through explicit examples in section 8 below. Robustness in that case stems from the fact that the number n of dissipative processes (or Lindblad operators) taking part in the dynamics is a well-controlled quantity in typical physical implementations. Intuitively, this can be understood from the fact that driven-dissipative processes considered here require a finite (typically large) amount of energy to occur and therefore do not arise unless they are deliberately introduced via quantum reservoir engineering (see section 2.3). In light of this physical constraint, dissipative perturbations of the type (i) discussed above are the most relevant ones to consider in assessing the robustness of Majorana zero-damping modes. We will focus on the latter.
We remark that there is a priori no guarantee to find Majorana zero-damping modes in the case where the system is driven by n ≥ N Lindblad operatorsindependently of the topological nature of the bulk-since, as discussed above, every generic dissipative process introduced in the dynamics leads to two additional constraints in mode space and further reduces the subspace available for Majorana zero-damping modes. Whereas spurious Majorana zero-damping modes can generally be removed by introducing additional dissipative processes-leaving no trace in the purity or damping spectra-genuine Majorana zero-damping modes, in contrast, must either survive or give rise to intrinsic Majorana zero-purity modes as dictated by the dissipative bulk-edge correspondence. Which of these two outcomes actually occurs cannot be determined from topological properties but rather depends on the dissipative boundary conditions, i.e., on the specific form of the edge dissipative dynamics resulting for the combined properties of bulk engineered dissipative dynamics and the physical edge delimiting it spatially. As we will demonstrate using explicit examples in sections 7,8 below, the dissipative dynamics can be constrained by e.g. geometric properties of the edge in such a way that the purity gap closes in a robust, controlled manner at that edge, thereby giving rise to an odd number of Majorana zero-damping modes in a phase whose bulk topology would naively imply the existence of an even number of such modes. Interesting phenomena with no Hamiltonian counterparts may therefore arise from Liouvillian dynamics as far as the edge physics is concerned, with intriguing effects such as the transformation of Abelian vortices into non-Abelian ones (see section 8).

Non-Abelian statistics of dissipative Majorana modes
A key property of Majorana modes in the Hamiltonian context is their non-Abelian statistics under spatial exchange. An important question is therefore whether or not this property is also present in the case of dissipative Majorana modes (more precisely, of Majorana zero-damping modes). Here we provide a simple affirmative argument following reference [30]: Physically, this may be understood from the fact that Majorana modes in the context of topological insulators or superfluids are not dynamical, but are purely static degrees of freedom. The actual bulk dynamics is irrelevant; it can be either unitary or dissipative.
In order to verify non-Abelian exchange statistics, we first explain in more detail how, in a dissipative setting, the edge subspace of the density matrix-containing the Majorana modes-is isolated from the bulk subspace of the latter. We then consider the effect of adiabatic parameter changes of the Liouville operator or, more generally, of the operator generating the system dynamics.

Dissipative isolation of the edge mode subspace
As argued in sections 2.6, 2.7 and 5, an important prerequisite for the presence of stable dissipative Majorana modes is the existence of a zero-mode subspace as well as its isolation from the bulk. Here we formulate the conditions for such a situation in a general framework that is valid beyond the quadratic setting introduced in section 2.4. To this end, we introduce projectors p and q = 1 − p on the edge and bulk subspaces, respectively. A decoupled edge subspace then appears if the Lindblad operators i are block diagonal in this projection, i.e., i,pq = i,qp = 0, with an edge block identical to zero, i.e., i,pp = 0. The dissipative evolution then reads Here, the bulk dissipative evolution L qq [ρ qq ] has a Lindblad form and ρ pp is a constant of motion, by construction. Crucially, the density matrix elements ρ qp and ρ pq coupling these sectors damp out according to ρ qp = e − i † i,qq i,qq t ρ qp (t = 0). In the presence of a dissipative gap, this process is exponentially fast. In steady state, the density matrix then takes a factorized form ρ = ρ edge ⊗ ρ bulk , with ρ edge = ρ pp and ρ bulk = ρ qq . Clearly, the absence of dynamics in the decoherence-free subspace does not allow for a controlled initialization of its occupation. This property is in complete analogy with a Hamiltonian ground-state scenario, where the zero-energy edge subspace is decoupled from Hamiltonian dynamics. The preparation and detection ideas developed in [68] for the ground states of fermionic atoms can be directly applied to the dissipative setting.
Moving from a second-to a first-quantized representation (for a quadratic setting and using a Majorana basis as in section 2.5), we obtain the evolution showing explicitly that the fadeout of bulk-edge coherences is indeed controlled by the dissipative gap, i.e., by the slowest rate in X qq .

Adiabatic parameter changes and braiding
We consider a steady state with Majorana zero-damping modes whose corresponding eigenvectors |α span a non-local decoherence-free subspace described by the density matrix elements (ρ edge ) αβ = α|ρ|β = ρ αβ . The decoherence-free subspace has the propertyρ αβ = 0. We now study the time evolution of the density matrix in a comoving basis |a(t) = U (t)|a(0) that follows the decoherence-free subspace of the edge modes, i.e., that preserves the propertyρ αβ = 0. Without specifying the actual dynamics generating the physical evolution, the time evolution in that basis is given by We define the vector potential A ab = a(0)|U †U |b(0) = i 2 ( a|ḃ − ȧ|b ), which is antisymmetric and Hermitian, by construction. Taking into account the normalization ∂ t b(t)|a(t) = 0, we then obtain whereρ ab ≡ a(t)|∂ t ρ|b(t) is the time evolution in the instantaneous basis. The Heisenberg commutator clearly reflects the appearance of a gauge structure [69][70][71][72][73] in the density matrix formalism. Crucially, this structure emerges independently of whether the physical dynamics encoded inρ ab -inserted in the last inequality-is unitary or dissipative. The transformation applied in the zero-mode subspace of either a Hamiltonian or a Liouvillian starting from the initial condition ρ αβ (0) is given by ρ αβ (t) = (V (t)ρ(0)V (t) † ) αβ , with time-ordered V (t) = T exp (−i t 0 dτ A(τ )) where A(t) αβ is the projection of the vector potential onto the decoherence-free subspace. The adiabatic change of the parameters is the key feature for such a state transformations to be realized while preserving the zero-mode subspace. Here one crucially requires the ratio between the rate of parameter changes encoded in A and the dissipative gap to be very small. Due to this separation of time scales-enabled by the non-evolving subspacethe decoherence-free subspace is never left. This phenomenon is sometimes referred to as the Quantum Zeno effect [57].
In the first-quantized formulation, the eigenvalues and eigenvectors of X now depend on time. The transformation into the instantaneous basis |a(t) = U (t)|a 0 , where |a 0 is the initial reference basis, is now orthogonal. The vector potential is given by the Hermitian antisymmetric matrix A ab = i 2 ( a(t)|ḃ(t) − ȧ(t)|b(t) ), and the correlation matrix evolves in this basis according to (h is defined above equation (13)) Starting form this basic understanding, one can construct adiabatic local parameter changes of the Lindblad operators to perform elementary dissipative Majorana moves [30]. Such procedure can then be applied sequentially in order to exchange two particular modes, keeping them sufficiently far apart from each other during the process (in networks of 1D wires, such exchange can be made using a T-junction geometry [16]). The unitary braiding matrix describing the complete process is B ij = exp π 4 γ i γ j for two Majorana modes i, j. Since [B ij , B jk ] = 0 for i = j, non-Abelian statistics is obtained.

Illustrative examples in 1D
In this section, we apply our results to several dissipative models in one dimension. In particular, we investigate the interplay between criticality and purity at the onset of nonequilibrium topological phase transitions, exemplifying the results of section 3.3. We focus on 1D lattice systems of spinless fermions with translational symmetry evolving under a dissipative dynamics described by Lindblad operators L n = m v n−m a † m + u n−m a m , where a † m and a m are creation and annihilation operators associated with a lattice site m and v n−m and u n−m are translation-invariant functions. The system is then most conveniently described in momentum space with Fourier-transformed Lindblad operators L k = m e ikm L m of the form Here we assume that as will be satisfied in the examples examined below. The dissipative dynamics obviously occurs in decoupled momentum sectors corresponding to the fermionic modes a † k and a −k . Defining a corresponding basis of Majorana operators (c 2k−1 , c 2k , c 2(−k)−1 , c 2(−k) ), the matrices X k and Y k describing the dynamics in each of these sectors take the following form (see section 2.5): where σ µ=x,y,z denotes the usual Pauli matrices. The system exhibits critical behavior if at least one of the eigenvalues λ (1,2) k = ±|u * k ± v k | 2 of X k vanishes for some k, i.e., if the dissipative gap closes. Note that in that case the matrix Y k identically vanishes, since its eigenvalues are given by ±(|u * k + v k |)(|u * k − v k |). All matrices X k of the form (62) can be diagonalized through the same unitary transformation. This allows us to obtain the steady-state correlation matrix in a generic form (solving the steady-state equation {X k , Γ k } = Y k ; see section 2.5); namely, The symmetry class to which the steady state belongs can be readily determined from this expression. We recall that one can associate a fictitious free-fermion Hamiltonian H Γ = i i,j Γ ij c i c j with the steady-state correlation matrix Γ (see section 3.2). Here we write H Γ = 2 k Ψ † k H k Ψ k in a momentum-space Nambu representation with Ψ † k = (a † k , a −k ) and obtain where N 2 k = |v k | 2 + |u k | 2 . Equation (61) implies that ξ −k = ξ k and ∆ * −k = ∆ k , and one can easily verify that time-reversal and particle-hole symmetries are satisfied as H * −k = +H k and σ x H * −k σ x = −H k , respectively (see equations (46) and (47)). Consequently, all (Gaussian) steady states resulting from a translation-invariant dissipative dynamics described by matrices X k and Y k of the form (62) and (63) have chiral symmetry and accordingly belong to the symmetry class BDI of Altland and Zirnbauer. In what follows, we present three examples of dissipative systems belonging to the symmetry class BDI exhibiting a topological phase transition due to the change of some external parameter κ. We first consider an example where the dissipative dynamics gives rise to a pure steady state in the whole parameter range of κ and the system exhibits critical behavior at the point where the topological phase transition occurs. We then discuss an example in which the topological phase transition is also driven by criticality but accompanied by the closure of the purity gap. Finally, we provide a third example in which the topological phase transition does not lead to any critical behavior. The three examples presented below therefore illustrate the three possibilities for non-equilibrium topological phase transitions identified in section 3.3.
We illustrate our results by plotting, as a function of κ, the dissipative gap ∆ d (given by the minimum eigenvalue of X) and the purity gap which we define here as ∆ p = min k {Tr(Γ 2 k )/4}. We also consider finite systems with open boundary conditions but translation-invariant Lindblad operators and investigate the behavior of the edge zero-damping modes in the vicinity of the topological phase transition.

Example 1: topological phase transition of a pure state with criticality
We consider translation-invariant Lindblad operators L n that act on three consecutive sites: where κ is a real parameter controlling the coupling to the central site n * . In momentum space, these Lindblad operators take the form L k = v k a † k − u k a −k with v k = κ + 2 cos k and u k = −2i sin k (note that equation (61) is satisfied). The corresponding matrix X k is diagonal, namely, X k = (8 + 2κ 2 + 8κ cos k)/(4 + κ 2 )I, and one can readily see that * Note that the overall multiplicative factor is introduced for convenience and does not affect the physical properties of the system.  (66) exhibits critical behavior at κ = ±2, as indicated by the closure of the dissipative gap ∆ d . At these points, the system undergoes a transition from a topologically trivial (ν 1D = 0) to a topologically ordered (ν 1D = 2) phase. Since the purity gap ∆ p = 1, the steady state of the dissipative dynamics is always pure. (b) Divergence of the localization length associated with one of the two Majorana zero-damping modes found at the edges upon approaching the topological phase transition at κ = 2.
Expressing the momentum-space correlation matrix in the form Γ k = i(n k · σ) with n k ∈ R 3 similarly as in section 3, we obtain n k = g k (0, 4κ(sin k + sin 2k), κ 2 + 4κ cos k + 4 cos 2k). As expected, this vector n k is non-zero for all k in the whole parameter range of κ where the dissipative gap is finite, i.e., for |κ| = 2. In that case one can "spectrally flatten" the steady-state correlation matrix Γ k by normalizing the vector n k to unity in order to calculate the corresponding winding number topological invariant ν 1D (see equation (35) and discussion above). Here we find ν 1D = 2 for |κ| < 2 and ν 1D = 0 for |κ| > 2 (for an explicit calculation, see our previous work [31]). We therefore identify a topological phase transition at |κ| = 2 where the dissipative gap closes and the system becomes critical. These results are illustrated in figure 3(a).
We now examine the case of a finite system (or "chain") of N sites with open boundary conditions. We first notice that only N − 2 three-site Lindblad operators of the form (66) can be "placed" on such an open chain. Hence, according to our discussion of section 5.3, four exact Majorana zero-damping modes (two at each edge, by symmetry) must be present independently of the topological properties of the system (in particular, for any value of the parameter κ). Obviously, this number which remains constant across the topological phase transition cannot be used to identify the latter. However, the topological phase transition can be revealed in two other ways: (i) One can remove spurious Majorana zero-damping modes with no topological origin by introducing additional arbitrary Lindblad operators acting at the edges; assuming that these extra Lindblad operators preserve chiral symmetry * , the system must then exhibit a total of ν 1D = 2 Majorana zero-damping and/or zero-purity modes in the parameter range |κ| < 2, as dictated by the dissipative bulk-edge correspondence (see section 5.1). (ii) One can examine the behavior of the localization length associated with the Majorana zero-damping modes as the parameter κ is tuned across the topological phase transition; owing to criticality, the characteristic length scale associated with at least one of the two modes found at an edge must diverge at the phase transition point |κ| = 2.
Here we follow strategy (ii): Since the Lindblad operators L n are translationinvariant and have a symmetric form around their central site n †, we use the results of section 5.2 and assess the existence of Majorana zero-damping modes that decay into the bulk by solving equation (54). Choosing φ i = 0, the latter equation takes the explicit form 0 = κ + 2β ex (e x being a unit vector joining neighboring sites of the chain) and we readily obtain β ex = −κ/2, showing that a least one of the two Majorana zero-damping modes at the edge must decay (exponentially) into the bulk when |κ| < 2 (as required by |β ex | < 1) on a characteristic length scale ξ = −1/ log (|β ex |) = −1/ log (|κ|/2). As expected, ξ diverges at |κ| = 2 and thus reveals the topological phase transition. Figure 3(b) illustrates numerical results corroborating this behavior.

Examples 2 and 3: topological phase transition of a mixed state with and without criticality
We now consider a dissipative dynamics that involves two types of translation-invariant Lindblad operators: acting on neighboring and next-to-nearest neighboring sites, respectively. The Lindblad operators L n generate a dissipative dynamics that drives the system into a pure steady state corresponding to the ground state of a topologically non-trivial Kitaev chain [15], as demonstrated in our previous work [30]; two Majorana zero-damping modes are found in that case (one at each edge). The Lindblad operators L (2) n describe two decoupled Kitaev chains, with four Majorana zero-damping modes (two at each edge). Below we discuss two scenarios: (i) We first assume that both dissipative processes L (1) n and L (2) n occur * Note that time-reversal symmetry only must be ensured since particle-hole symmetry is automatically satisfied in our dissipative framework.
†Namely, the creation and annihilations parts of L n (see equation (66)) have s-wave and p-wave symmetry, respectively.   coherently, so that the relevant Lindblad operators are L n = (L (1) n +κL (2) n )/(2(κ 2 +κ+1)) with κ ∈ R. (ii) We then consider the case where both dissipative processes "compete" against each other, so that the relevant Liouvillian to describe the dissipative dynamics is L = (L 1 + κL 2 )/(1 + κ), where L α denotes the Liouvillian corresponding to a single Lindblad operator L n as well as their combinations (i) and (ii) can be realized in a "zigzag" geometry, e.g., as depicted in Fig. 4.
To conclude our investigation of case (i), we consider a finite chain of N sites with open boundary conditions. Since only N − 2 Lindblad operators of the form L n ∝ L (1) n + κL (2) n can be applied on such a chain, we again obtain four exact Majorana zero-damping modes that are either genuine or spurious. A numerical analysis below the topological phase transition point κ = 1 reveals that two of these modes (one at each edge) decay exponentially into the bulk on a characteristic length scale ξ = −1/ log(κ) which diverges at κ = 1, as expected (see figure 5(b)).
We now turn to the investigation of case (ii) in which the two dissipative processes L (1) n and L (2) n compete against each other. In that scenario, one can readily verify that X k = 2I, such that the damping spectrum is gapped (and "flat") for all k. The steady- n , respectively, exhibits a topological phase transition at κ = 1 from a state with winding number ν 1D = 1 to a state with winding number ν 1D = 2 that results from the closure of the purity gap ∆ p only. The system never becomes critical, since the dissipative gap ∆ d remains constant in the whole parameter range of κ. state correlation matrix is then given by Γ k = g (cos k + κ cos 2k)iσ y −(sin k + κ sin 2k)σ z (sin k + κ sin 2k)σ z (cos k + κ cos 2k)iσ y , with g = 1/(1 + κ) and eigenvalues ± √ 1 + κ 2 + 2κ cos k/(1 + κ), implying that the steady state is pure if and only if κ = 0 or κ → ∞. The corresponding vector n k is non-zero (for all k) except at κ = 1 in which case Γ k=±π = 0. For κ = 1, one finds n k = g(0, −(sin k + κ sin 2k), −(cos k + κ cos 2k)) and the corresponding winding number topological invariant is ν 1D = 1 for κ < 1 and ν 1D = 2 for κ > 1. These results, which are presented in figure 6, demonstrate that case (ii) provides an example of a topological phase transition that is not driven by criticality but rather by the closure of the purity gap only.
Examining case (ii) on a finite chain with open boundary conditions, we find four Majorana zero-damping modes (two at each edge) in the whole parameter range of κ. In accordance with the dissipative bulk-edge correspondence (see section 5.1), all of these modes must be genuine for κ > 1 where ν 1D = 2, while two of them (one at each edge) must be spurious for κ < 1 where ν 1D = 1. Crucially, the spatial wave function of these four modes does not exhibit any critical behavior * upon approaching the topological phase transition point κ = 1.

Illustrative example in 2D
In the previous section, we have exemplified the phenomenology of Majorana zerodamping modes at physical edges in the context of 1D systems. Here we focus instead on the physics that arises when topological defects are introduced in the bulk of the system, away from physical edges. We consider topological defects in the form of dissipative vortices in infinite 2D systems in order to illustrate such physics in the simplest possible scenario.

Dissipative vortices
Dissipative vortices are most conveniently defined starting from an infinite 2D lattice system evolving under a dissipative dynamics generated by translation-invariant Lindblad operators L i = j u j−i a j +v j−i a † j (using notations as in section 5.2 above; see equation (52), in particular). We assume that the corresponding Liouvillian exhibits a finite dissipative gap and that its steady state is characterized by a finite purity gap, so that the topological nature of the system is well-defined (i.e., the steady state belongs to a specific topological class and is characterized by quantized topological invariants); in particular, the dissipative bulk-edge correspondence introduced in section 5 is satisfied * . We then introduce a dissipative vortex at a position defined by the vector R 0 † by modifying the translation-invariant coefficients u j−i ≡ u ij defining the annihilation part of the Lindblad operators in the following way: where (r j , ϕ j ) are polar coordinates defining the position of each site j with respect to R 0 , f (r) is a real and positive function that vanishes as r → 0 and reaches a constant value as r δ (δ being a characteristic length scale associated with the dissipative vortex, defining the vortex core), and e −i ϕ j describes the vortex phase winding times around the origin R 0 ( defining the so-called vorticity). While vortices exhibit similar properties in Hamiltonian systems, the specific form of a dissipative vortex defined above is motivated by its natural realization in typical implementation schemes with cold atoms based on optical vortex imprinting (we refer to our previous work [31] for further details).
We remark that the qualitative features of a dissipative vortex do not depend on the explicit form of the function f (r) describing its core. Of crucial importance is the fact that f (r) vanishes as r → 0, so that the vortex core can be identified with a region of space in which the system is topologically equivalent to the vacuum, with all sites essentially occupied ‡. If the topological nature of the bulk is non-trivial, a small circular * One can assume, without loss of generality, that the Lindblad operators form a complete set of anticommuting operators, so that the system is driven into a pure steady state independently of the initial conditions (see discussion of section 3.1).
†Note that this position need not coincide with a lattice site. ‡This can be understood from the fact that the Lindblad operators acting in the vortex core have an annihilation part that essentially vanishes, namely, edge (or domain wall) of radius ∼ δ must then appear around the vortex core, separating this vacuum-like region from the surrounding bulk which is potentially topologically nontrivial. Remembering the dissipative bulk-edge correspondence discussed in section 5, one then naively expects Majorana zero-damping modes and/or intrinsic Majorana zeropurity modes to appear at this domain wall. We argue below that the possibility of having such modes, however, additionally and crucially depends on the phase winding of the dissipative vortex, as expected by analogy to the Hamiltonian setting. The crucial role of the vortex phase winding in prohibiting or allowing the existence of Majorana zero modes in the vortex core can be understood in analogy to the Hamiltonian setting. In the Hamiltonian context, circular domain walls typically trap low-energy modes with quantized angular momentum m and a corresponding energy E ∼ m [12]. When a flux of the U (1) gauge field associated with the fermions corresponding to n quanta of angular momentum is inserted in the region enclosed by the domain wall, the angular momentum of these low-energy modes shifts to values m+n and the energy of the latter shifts accordingly, thereby prohibiting or allowing for zeroenergy modes. We argue that the same argument applies here in the case of a dissipative vortex, namely: the phase winding defines a flux π threading the vortex core and shifts the eigenvalues corresponding to low-damping or low-purity modes bound to the vortex core (in the damping or purity spectrum, respectively). This can be understood from the fact that the Hamiltonian phenomenology discussed above automatically applies to the parent Hamiltonian associated with the dissipative dynamics, whose zero modes either correspond to Majorana zero-damping modes or to intrinsic Majorana zero-purity modes (see discussion of section 2.7). In analogy with the Hamiltonian setting (see, e.g., [12]), we then expect Majorana zero-damping modes and/or intrinsic Majorana zero-purity modes to appear in dissipative vortices with odd vorticity only, as we corroborate below in an explicit model.
We remark that the physics associated with a dissipative vortex remains qualitatively the same upon insertion or removal of 2π fluxes in the vortex core since such modifications can be seen as gauge transformations [12]. The physical properties of a dissipative vortex are therefore determined by the vorticity modulo 2, and one can restrict oneself to = 0 or = 1, without loss of generality. Since the case of a dissipative vortex with = 0 is physically equivalent to that of a small circular physical edge, we will focus exclusively on dissipative vortices with = 1. Remarkably, the physics arising in a single dissipative vortex with a π flux ( = 1) in an infinite system on the plane can be shown to be equivalent that occurring at the edge of a semi-infinite system with a cylinder geometry and no flux (see figure 7 and our previous work [31] for further details). This equivalence-which will be useful in the sections below-can intuitively be understood from the fact that the π flux of the U (1) gauge field introduced by a vortex with = 1 naturally arises in cylinder geometry due to the (extrinsic) curvature of the latter. . Left: Semi-infinite plane "wrapped up" into a semi-infinite cylinder in the e y direction (shown in red). Right: Topological equivalence between the semi-infinite cylinder geometry and an infinite planar geometry with a single dissipative vortex "puncturing" the plane.

Explicit 2D model and bulk properties
We now consider, in the framework defined above, an explicit model defined on a square lattice with primitive vectors e x and e y , with translation-invariant Lindblad operators L i = j u j−i a j + v j−i a † j whose only non-zero coefficients are v 0 = β; v ±ex = 1; v ±ey = 1; where β is a real parameter which can be used to tune the system across phase transitions. Every site i therefore has an associated Lindblad operator L i of the form where the creation part C † i is isotropic or s-wave symmetric with respect to the central site i (as naturally arises in typical implementation schemes; see section 2.3) and the annihilation part A i has a p-wave symmetry, i.e., A i = (∇ x + i∇ y )a i with ∇ λ a i ≡ a i+e λ − a i−e λ . In momentum space, these translation-invariant Lindblad operators take the form L k = u k a k + v k a † −k , where u k = 2i(sin (k x ) + i sin (k y )); v k = β + 2(cos (k x ) + cos (k y )). (76) Since the above functions satisfy u k v k = −u −k v −k , the Lindblad operators form a complete set of anticommuting operators and the system is driven into a pure Gaussian steady state corresponding to the ground state of the parent Hamiltonian H parent = k L † k L k (see equation (39) and discussion thereof). The damping spectrum, which coincides with the spectrum of H parent in that case, is then given by the norm Remembering the explicit form of the functions u k and v k defined in equation (76) above, one can easily verify that the corresponding dissipative gap closes at the parameter values β = 0, ±4. We thus expect the existence of four distinct-possibly topologicalphases depending on the value of β. We remark that the Lindblad operators are not invariant under time-reversal symmetry due to the p-wave symmetry embedded in their annihilation part. The system thus belongs to the symmetry class D of Altland and Zirnbauer (see section 3.2) and is characterized by the Chern number topological invariant defined in equation (32). One can easily verify that, in accordance with the discussion of section 4.2, the Chern number vanishes in the whole β parameter range where the system exhibits a finite dissipative gap (we refer to our previous work [31] for an explicit proof). As dictated by the dissipative bulk-edge correspondence introduced in section 5, we therefore naively do not expect to find Majorana zero-damping modes in the system if edges are introduced. Despite this conclusion, let us now try to construct such modes explicitly in the translation-invariant setting of section 5.2 anyway.

Construction of Majorana zero-damping modes
The possibility of having Majorana zero-damping modes at a physical edge when the dissipative dynamics is generated by translation-invariant Lindblad operators as in section 5.2 above can be assessed from the explicit form of the Lindblad operators (see equation (74)) by looking for a solution of equation (54). Choosing φ i = 0 (which corresponds to choosing a particular orientation of the edge; see section 5.2), the real and imaginary parts of equation (54) here take the explicit form 0 = β + 2β ex + β ey + (β ey ) −1 , 0 = β ey − (β ey ) −1 .
One thus finds β ey = ±1 and β ex = −(β/2 ± 1), which implies the existence of a solution with |β ex | < 1 if and only if 0 < |β| < 4. Consequently, there exists at least one Majorana zero-damping mode γ of the form defined by equation (53) when 0 < |β| < 4 if the system possesses an edge in the −e x direction * , in which case γ is uniformly spread along the edge (since |β ey | = 1) and exponentially localized in the +e x direction away from the edge on a characteristic length scale ξ = −1/ log (|β ex |) = −1/ log (|β/2 ± 1|). The semi-infinite planar geometry of the system in that case is depicted in figure 7(a). Remarkably, the special points β = 0, ±4 where the dissipative gap closes coincide with those at which |β ex | = 1, i.e., at which the localization length associated with γ diverges. This suggests the possible existence of topological phase transitions occurring at the gap-closing points |β| = 0, ±4 which are not identified by the (vanishing) Chern number topological invariant.

Topological origin
We have demonstrated above that the system can support Majorana zero-damping modes at an edge (in the parameter range 0 < |β| < 4) despite its vanishing Chern number. We argue below that this apparent contradiction stems from the fact that we have assumed translational symmetry, and establish the topological origin of Majorana zero-damping modes in the above translation-invariant setting. It is clear that translational symmetry is not robust against disorder. However, assuming such a symmetry will allow us to demonstrate, in the simple model above, a key mechanism with no Hamiltonian counterpart through which spatially isolated Majorana zerodamping modes can be obtained in a topological phase characterized by an even integer topological invariant, in contradiction with the Hamiltonian bulk-edge correspondence.
We first remark that the construction of section 8.3 above remains unchanged if we assume that the system is finite in the e y direction, "wrapped up" into a cylinder (as depicted in figure 7(b)). Since such a semi-infinite cylinder geometry is topologically equivalent to that of an infinite plane "punctured" by a single dissipative vortex (figure 7(c)), Majorana zero-damping modes of the form found in section 7 must therefore similarly appear in the core of a dissipative vortex on the plane * . In these two equivalent geometries, the topological origin of Majorana zero-damping modes can be identified by moving to momentum space in the e y direction associated with translational (or rotational, in the vortex case) symmetry (see figures 7(b) and (c)). The system then reduces to a "stack" of semi-infinite 1D wires (in the e x direction) associated with different momenta k y . In particular, the Fourier-transformed counterpart of the Lindblad operators defined in equation (75) take the form in the momentum sectors corresponding to k y = 0 or π, where i indexes the lattice sites in the remaining spatial direction e x and κ = β + 2 for k y = 0 and β − 2 for k y = π, respectively. In these two sectors, the system thus reduces to the 1D wire with chiral symmetry investigated in section 7 above (up to an overall dissipation rate √ 4 + κ 2 which does not affect any of the topological properties of the system). Owing to chiral symmetry, such a system belongs to the symmetry class BDI of Altland and Zirnbauer and is characterized by the winding number topological invariant ν 1D defined by equation (35). As discussed in section 7 and depicted in figure 3 thereof, one finds ν 1D = 2 for |κ| < 2 and ν 1D = 0 for |κ| > 2. Discontinuities in the winding number at κ = ±2 (or, equivalently, at β = 0, ±4) clearly establish the occurrence of nonequilibrium topological phase transitions at the dissipation gap-closing points identified in section 8.3.
In summary, we have identified a non-trivial topological invariant ν 1D = 2 in the parameter range 0 < |β| < 4 † by taking advantage of translation invariance to reduce * Note that the translation symmetry along the circumference of the cylinder translates as a rotational symmetry around the vortex core.
the original 2D problem in semi-infinite cylinder geometry or, equivalently, in planar geometry with a single = 1 dissipative vortex, to a 1D wire defined along the axis of the cylinder or in the radial direction away from the vortex core (e x direction shown in figures 7(b) and (c)). This provides a topological explanation as to why we were able to demonstrate the existence of at least one Majorana zero-damping mode for 0 < |β| < 4 in section 8.3. In fact, invoking the dissipative bulk-edge correspondence of section 5.1, we can now argue that a total of ν 1D = 2 Majorana zero-damping modes and/or Majorana zero-purity modes must be present at the edge of the semi-infinite system in cylinder geometry or, equivalently, in the core of the = 1 dissipative vortex on the (infinite) plane, as long as the symmetries of the system are preserved * .
The precise number (m damping ) edge ≤ 2 (see equation (51) and discussion thereof) of genuine Majorana zero-damping modes that is found at the edge of the cylinder or in the vortex core depends on the specific dissipative boundary conditions that appear there. Crucially, physically distinct boundary conditions naturally emerge in these two systems, although their geometries are topologically equivalent. In the case of the dissipative vortex, in particular, one can demonstrate explicitly that a single Majorana zero-damping mode is found in the vortex core (alongside with a single Majorana zero-purity mode such that the dissipative bulk-edge correspondence is satisfied; see equation (51)). The key mechanism through which the dissipative vortex isolates exactly one Majorana zero-damping mode despite the fact that the bulk is characterized by an even (integer) topological invariant (ν 1D ) originates from the geometry of the vortex core alone and is therefore generic. We refer to our previous work [31] for an explicit proof of this statement.

Physical properties of dissipative vortices carrying single Majorana zero-damping modes
In the Hamiltonian context, π flux vortices carrying single Majorana zero-energy modes have been demonstrated to exhibit non-Abelian exchange statistics [8]. The underlying proof exclusively relies on the algebraic properties of Majorana modes and, crucially, does not depend on the dynamics giving rise to them (see section 6.2). As a result, the same conclusions apply in the dissipative setting of interest in this work. In particular, π flux (i.e., = 1) dissipative vortices carrying single Majorana zero-damping modes must exhibit non-Abelian exchange statistics (when moved around through adiabatic parameter changes as discussed in section 6.2). In light of the mechanism illustrated in the above explicit model, we thus conclude that vortices can have, in our dissipative setting, non-Abelian exchange statistics in a phase characterized by an even integer topological invariant, in stark contrast to what can occur in the Hamiltonian context. . Note that, due to the finite size of the system, vortices "interact" more strongly (via the edge) as they come closer to the edge, giving rise to the symmetric behavior that can be seen around d ≈ 16. (c) and (d): Low-lying part of the damping (respectively purity) spectrum for β = 2 featuring a gap with two quasi-zero (∼ 10 −15 ) eigenvalues (indicated by a red arrow) corresponding to Majorana zero-damping (zero-purity) modes trapped in each of the vortex cores. All results were obtained for vanishingly small vortex cores, i.e., with f (r) = 1 everywhere except at r = 0 where it vanishes (see equation (73)).
In addition to their non-Abelian statistics, vortices carrying single Majorana zero modes behave, in our dissipative setting, in a fully analogous way as in the Hamiltonian context. In particular, two such vortices "interact" as the distance d between the latter becomes small, acquiring an exponentially small residual damping rate ∼ e −d/ξ , where ξ is the localization length associated with the Majorana zero-damping mode that they carry. We have verified such phenomenology through numerical simulations. Figure 8 summarizes the results obtained in the case of a finite system on the plane with two "interacting" = 1 dissipative vortices.

Conclusions
In this paper we have provided a discussion of non-equilibrium superfluid topological states generated by means of engineered dissipation, highlighting analogies and differences to the more conventional Hamiltonian setting. Similarities mainly stem from the fact that topological properties are properties of the state of the systemencoded in the static correlations-with some care to be taken regarding the purity of the state. In particular, identifying the correlation matrix as a tool for the symmetry-based topological classification of arbitrary Gaussian states analogous to a first-quantized quadratic Hamiltonian, it is clear that the classification of topological states according to the symmetry classes of Altland and Zirnbauer extends to a general non-equilibrium steady-state situation, highlighting the universality of this classification. Differences, on the other hand, mostly stem from the fact that the spectral or dynamical properties of the system and the properties of the steady state are not as tightly related as in the Hamiltonian ground-state scenario. The interplay between the well-known characteristics of topological states and the new elements brought in by the dissipative nature of the dynamics can give rise to effects with no Hamiltonian counterpart. Examples are topological phase transitions which do not require the dissipative spectral gap to close, or the trapping of unpaired Majorana modes in a bulk with vanishing Chern number.
We have presented in this work a detailed discussion of dissipatively induced topological superfluid states in 1D and 2D systems which is a first step towards a complete picture of non-equilibrium topological order. In this respect, various directions remain to be explored in future work: From the viewpoint of dissipation engineering, we may ask ourselves whether other symmetry classes than the ones identified in this work can be achieved in realistic setups. So far we have concentrated on topological superfluids of spinless fermions. Adding spin degrees of freedom opens up the perspective of reaching topological insulating phases with potentially more robust edge modes. It will also be interesting to investigate whether going beyond quasi-local Lindblad operators is possible in realistic cold atom experiments, allowing to reach phases characterized by a non-vanishing Chern number.
From a many-body perspective, an interesting direction will be to explore the nature of topological phases and phase transitions in more detail. This concerns, in particular, the interplay of dissipative and Hamiltonian dynamics, which we did not examine in this work. While this intentional omission is not a severe limitation in the context of cold atomic gases where the Hamiltonian dynamics can be suppressed, including the effects of such dynamics would give a more complete picture of the robustness of topological phases under general combined unitary and dissipative dynamics. From the viewpoint of critical phenomena, topological phase transitions accompanied by a closure of the dissipative gap offer intriguing examples of criticality in fermionic systems. Going beyond a mean-field approach by including the effects of Hamiltonian and dissipative interactions will be necessary to reach a detailed understanding of this issue.
From a quantum information point of view, it will be important to investigate in more detail the analogies and differences between dissipatively induced non-local decoherence-free subspaces and more conventional Hamiltonian edge mode subspaces. To answer these questions, a more systematic classification of harmful dissipative perturbations as well as a more detailed understanding of the dynamics of excitations on top of the dissipatively stabilized state will be necessary. The propagation of excitations or defects on top of a state which is entirely induced by dissipation can be expected to be vastly different from a Hamiltonian ground state. Also, manipulation and readout tools tailored to the dissipative setting will have to be developed. and dark state there is a fixed number setting with * The reverse statement holds as well. Here, as in the main text the creation (annihilation) parts C † i (A i ) are linear in the spinless fermion operators and obey the following requirements: The Fourier transforms for the creation and annihilation part are local in momentum space, The Fourier transformed fixed phase and fixed number Lindblad operators then read, respectively (ii) Antisymmerty -The functions have the property u −k = ±u k , v −k = ∓v k , such that the related "wave function" ϕ is antisymmetric, This implies that the set L k form a full Dirac algebra up to a normalization, for all k, k . We may introduce normalized operators viā . The generator of the state is bilinear and reads in terms of the momentum space wave function The fixed number operators satisfy [ i ,N ] = 0 and thus conserve total particle numberN = i a † i a i . The corresponding dark states |BCS, N contain 2N fermions. Instead, in the fixed phase setting, [L i ,N ] = 0. α = re iθ is an arbitrary complex number. Within the mean field described below, the modulus r can be related to the average particle number. θ has an interpretation in terms of a fixed, superfluid phase. This is made explicit from the expansion of the coherent state wave function where N is a global normalization. k is restricted to k such that k λ > 0 in all primitive lattice directions λ.
With these preparations, we can prove the above statement. We work in momentum space and start with the number conserving setting. The Lindblad operators k = q C † q−k A q are normal ordered, such that k |vac = 0. The dark state property i |BCS, N = k |BCS, N = 0 for all i or k is therefore equivalent to following relation for all k, This is true if and only if vqu q−k uqv q−k = ϕq ϕ q−k holds for all k, i.e., for the wave function ϕ q = v q /u q up to a momentum independent complex number. The latter can be absorbed into the constant α.
It is easy to show the uniqueness of the dark state in the fixed phase ensemble. The coherent state wave function (A.10) can be written as |BCS, θ ∝ k L k L −k |vac , so that it is indeed the unique wave function annihilated by the full set of Dirac operators L k up to normalization. We are not aware of an analogous argument for the fixed number case (for each given particle number 2N ), since the bilinears k do not obey a simple algebra. Numerical simulations in the context of spinful Lindblad operators [52,53] however suggest a unique steady dark state.
The two wave functions |BCS, N , |BCS, θ are equivalent in the thermodynamic limit N → ∞: We consider the relative number fluctuations of the fixed phase state, whereN = kn k ,n k = c † k c k is the total particle density, the average is taken in |BCS, θ , and n k = n k = |v k | 2 . Here we have usedn 2 k =n k , implying n 2 k = n k , as well as n qnk = n q n k for k = q due to the product nature of the BCS state. Since 0 ≤ n k ≤ 1 and n k = 0 only for a non-extensive set of modes, the scaling of the numerator is ∼ N , and that of the denominator ∼ N 2 in the thermodynamic limit. It is the mean-field (product) nature of the state which is responsible for this scaling.
The fixed phase wave function is an exact solution of the number conserving equation (if the particle number is not fixed to a specific N ), since |BCS, θ ∝ exp(αG † )|vac = N (αG † ) N /N !|vac . The above discussion shows that for fixed N → ∞, the fixed phase wave function approximates the fixed number one in a controlled way. We will rely on this fact, which is solely related to large N , in the mean-field theory to be discussed next. many reorderings are necessary: where in the last equality we have chosen the normalization Tr k ρ k = 1 and the cyclic invariance of the trace. In the thermodynamic limit, q =p (.) = q (.). Close to the steady state we can evaluate the correlation functions using the stationary form equation (A.10) for the state, which take a factorized form, i.e., a † q a q = |v q | 2 , a q a † q = |ū q | 2 , a −q a q = −ū * qv q , a † q a † −q = −v * qū q . We can thus factorize a real positive number κ 0 from the last equation, such that mean-field Liouvillian operator for the mode pair {k, −k} takes the form In addition to the effective dissipative rate, we can also calculate the amplitude of the relative coefficient α in dependence of the particle filling n = q a † q a q from meanfield theory. While for the number conserving equation the particle number operator is a constant of motion due to [ i ,N ] = 0 (i.e., the expectation values of any power of the total particle number are conserved, ∂ t N m = 0), no such property is found for the mean-field dynamics. However, for any initial average particle filling n = N /L d there is an α = re iθ such that the first moment, i.e., the average particle density, is conserved, and ∂ t N = 0 holds. In other words, the value of α fixes the average particle number, which is intuitive since α measures the relative strength of creating fermions vs. annihilating fermions. The number equation of state for which n(t) = const. is given by (A. 16) While, therefore, the modulus |α| = r is fixed by an additional physical condition from the initial state, no such condition exists which fixes the phase θ. This is physically meaningful, since the original number conserving operators are invariant under global U (1) phase rotations. The interpretation of α = re iθ therefore is the following: (i) the macroscopic phase is fixed spontaneously by the dynamics, i.e., it indicates spontaneous symmetry breaking of U (1) phase rotations in the dissipative setting. (ii) The amplitude relates to the average particle number in the sample, which can be chosen such as to match the average particle number of the number conserving setting n = N/L d in the initial state. In this section, we provide an explicit derivation of equation (54) in the framework introduced in section 5.2 (using the same notations).
An arbitrary Majorana operator γ = γ † can be expressed, in terms of the fermionic creation and annihilation operators a † i and a i corresponding to distinct lattice sites i of the system, in the form where the sum runs over all sites j belonging to the system and α j ∈ C with j |α j | 2 = 1 (so that {γ, γ † } = 2 as required for a Majorana operator). The operator γ then corresponds to a Majorana zero-damping mode of the dissipative dynamics described by the translation-invariant Lindblad operators L i = j∈I(i) u j−i a j +v j−i a † j (equation (52) where γ i denotes the restriction of the Majorana operator γ to the sites j ∈ I(i) where L i acts non-trivially. These anticommutation conditions obviously constrain the possible form of γ i in a way that solely depends on the form of L i . Since the form Lindblad operator L i does not depend on i due to translation invariance, γ i must be of the form where β j−i are complex factors (β j−i being a shorthand notation for β(r j −r i ), as defined in section 5.2) which, for consistency with equation (A.17), must satisfy the following set of equations: leading to a single "consistency relation" which is valid for any triple of indices (i, j, k) associated with lattice sites belonging to the system, with β i−i ≡ β 0 = 1. Since β j−i = (β i−j ) −1 and β * j−i = (β * i−j ) −1 , the factors β j−i must be real.
Equation (A.22) implies that any Majorana zero-damping mode γ can be fully constructed from the only knowledge of d real factors β en , where e n (n = 1, 2, . . . , d) are primitive vectors associated with the d-dimensional Bravais lattice on which the system is defined. Such a construction can be made starting from any lattice site i, choosing a phase φ i ∈ [0, 2π) for the coefficient α i of γ (see equation (A.17)) corresponding to that particular site. This allows us to express γ in the generic form presented in the main text (see section 5. where in the last step we have symmetrized the expression using the fact that I(i) contains pairs of sites (j, j) located symmetrically around the center of symmetry S i of L i (see section 5.2). The fact that the Lindblad operators are symmetric around S i and that the dissipative dynamics leads to a pure steady state as assumed in the main text implies that u j−i = −u j−i and v j−i = v j−i (see also equation (39) and discussion thereof). Equation (A.24) therefore reduces to which is the result presented in the main text (see equation (54)). We note that this equation can be further simplified when the center of symmetry S i of L i coincides with a lattice site. In that case, β j−i = β i−j = (β j−i ) −1 and we find real factors β en satisfying equation (A.25) for some particular phase φ i ∈ [0, 2π), the system can support at least one Majorana zero-damping mode of the form (53). Otherwise there can be no such mode. In order to assess the possibility of having Majorana zero-damping modes, one therefore has to solve a polynomial equation in d real variables ({β en } n=1,2,...,d ) with coefficients that may be real or complex depending on the coefficients e −iφ i u j−i , since v j−i can always be chosen as real (see section 4). The existence of solutions thus depends on the phase φ i that is chosen as an "initial condition" for the construction of γ starting from site i according to equation (A.23). We finally discuss the physical meaning of the phase φ i in the two types of dissipative systems that are most relevant to this work (see section 4): (i) 1D dissipative systems with TRS (symmetry class BDI) and (ii) 2D dissipative systems without TRS (symmetry class D) * . In the 1D case, TRS constrains the coefficients u j−i to be real up to a global phase and the phase φ i crucially allows to compensate for such a phase so that all coefficients e −iφ i u j−i are real. Equation (A.25) then reduces to a univariate polynomial equation with real coefficients which, as opposed to a univariate polynomial equation with complex coefficients (as would generically be obtained in the absence of TRS), can have robust solutions (i.e., solutions that survive if the coefficients u j−i and v j−i are slightly modified). In the 2D case, v j−i typically exhibits an isotropic (or s-wave symmetric) form due to physical constraints (see section 4) and TRS can only be broken if u j−i ≡ u(r j − r i ) transforms under rotations as an eigenfunction of angular momentum with odd integer eigenvalue = 0, namely, u(r) = u(r, ϕ) ∼ f (r)e −i ϕ with polar coordinates (r, ϕ), f (r) being an arbitrary positive function of r. We then find e −iφ i u j−i = e −iφ i u(r, ϕ) = u(r, ϕ + φ i / ), showing that modifying the phase φ i simply corresponds to rotating the reference frame or, equivalently, to modifying the orientation of the edge of the system. If there exists a solution of equation (54) for some φ i , there must exist a solution for any φ i ∈ [0, 2π) or, equivalently, for any orientation of the edge. TRS breaking thus crucially allows for the existence of Majorana zero-damping modes at edges that are curved arbitrarily, as expected physically.