Symmetry classes of open fermionic quantum matter

We present a full symmetry classification of fermion matter in and out of thermal equilibrium. Our approach starts from first principles, the ten different classes of linear and anti-linear state transformations in fermionic Fock spaces, and symmetries defined via invariance properties of the dynamical equation for the density matrix. The object of classification are then the generators of reversible dynamics, dissipation and fluctuations featuring in the generally irreversible and interacting dynamical equations. A sharp distinction between the symmetries of equilibrium and out of equilibrium dynamics, respectively, arises from the different role played by `time' in these two cases: In unitary quantum mechanics as well as in `micro-reversible' thermal equilibrium, anti-linear transformations combined with an inversion of time define time reversal symmetry. However, out of equilibrium an inversion of time becomes meaningless, while anti--linear transformations in Fock space remain physically significant, and hence must be considered in autonomy. The practical consequence of this dichotomy is a novel realization of antilinear symmetries (six out of the ten fundamental classes) in non-equilibrium quantum dynamics that is fundamentally different from the established rules of thermal equilibrium. At large times, the dynamical generators thus symmetry classified determine the steady state non-equilibrium distributions for arbitrary interacting systems. To illustrate this principle, we consider the fixation of a symmetry protected topological phase in a system of interacting lattice fermions. More generally, we consider the practically important class of mean field interacting systems, represented by Gaussian states. This class is naturally described in the language of non-Hermitian matrices, which allows us to compare to previous classification schemes in the literature.


I. INTRODUCTION
The distinction between different unitary and antiunitary symmetries [1][2][3][4] is a powerful organizing principle in the classification of quantum matter. It has been spectacularly successful in the description of gapped fermionic matter, where the identification of topologically twisted ground states on the background of ten fundamental symmetry classes [5] culminated in the periodic table of topological insulators and superconductors [6,7].
In this paper we ask how the concept of symmetry classifications can be generalized to fermion matter pushed out of thermal equilibrium by an external environment. This question is motivated in part by recent experimental progress in the physics of condensed matter, atomic condensates, and optics, which led to the realization of novel phases of quantum matter in engineered environments [8][9][10][11][12][13][14][15][16]. These developments call for the classification of symmetries and topologies of open quantum matter, in extension of existing frameworks for closed system quantum ground states.
Earlier work in this direction has put the emphasis on the most apparent consequence of environmental coupling, lossy dynamics and its description in terms of non-Hermitian matrix operators [17][18][19][20][21]. We here take a more general perspective and note that a comprehensive description of out of equilibrium symmetries must account for the interplay of dissipation and environmental fluctuations. Our approach to the full problem starts with the realization just how straightforward the description of symmetries in the complementary case of isolated sys-tems actually is. There, the full information is stored in the symmetries of a single Hermitian operator, the system Hamiltonian,Ĥ. The latter encodes the symmetries on the microscopic level via the definition ofĤ from Fock space operators, it describes the symmetries of state evolution via the evolution operatorÛ = exp(−iĤt) ( = 1), and those of long time stationary states through projectors onto the many body ground state ofĤ, or, slightly more generally, a thermal distributionρ ∼ exp(−Ĥ/T ).
Out of equilibrium, the situation becomes distinctly more complex. How does an exhaustive set of symmetries describing a system out of equilibrium look like? Which elements of the theory assume the role of the Hamiltonian in the description of these symmetries? And on the basis of what physical principles should they be described in mathematical terms? The formulation of concrete, and surprisingly simple answers to these questions is the mission of the research reported in this paper.
To get warmed up to the subject, consider a system whose dynamics is subject to damping, external driving, and an intrinsic Hamiltonian. In its evolution, these influences manifest themselves in the competition of dissipation, fluctuations, and unitary dynamics (cf. Fig. 1). We aim to understand in which ways these three are constrained by symmetries. On top of the concrete challenges formulated above, we are immediately facing a fundamental issue, the status of time reversal. Six out of ten of the fundamental symmetries of equilibrium quantum matter make reference to this symmetry. However, out of equilibrium, time reversal looses its physical meaning. Does this mean that the number of symmetries reduces  (47).) Out of equilibrium, the dynamics is specified in terms of three generators,Ĥ,D,P , representing unitary evolution, dissipation, and fluctuations, respectively. In either case, initial Fock space states evolve into stationary states whose symmetries are determined by the dynamical generators, but not the initial configurations. b) In equilibrium dynamics, fluctuation and dissipation are locked by the fluctuation-dissipation relation, Eq. (47), reducing the 'operator space' to two freedoms, H andD controlling the relaxation into an equilibrium configuration. Anti-unitary symmetries appear in combination with an inversion of time. The representation of anti-unitary symmetries, T, depends on whether the dynamics is unitarŷ D =P = 0 (T in combination with time reversal, t → −t), irreversible equilibriumD ∝P (T in combination with thermal time reversal, t → −t + iβ), or irreversible out of equilibrium (T in autonomy).
from ten to four? The answer is no. To understand what is happening, recall that time reversal in quantum mechanics is implemented by a combined operation inverting time, t → −t, and subjecting operators to an anti-unitary transformation. Out of equilibrium, the former looses its meaning, but the latter does not. We are thus led to investigate the status of anti-linear operations in autonomy.
This realization, surprisingly, appears to be novel, and it leads to unexpected structures. To mention one example, we call a quantum Hamiltonian 'chiral', stabilizing a chiral equilibrium phase, if it anti-commutes with a Pauli matrix σ 3 in some representation. However, it turns out that a Hamiltonian participating in the stabilization of a chiral out of equilibrium phase must commute with σ 3 . Such constraints have practical bearings, for example for the engineered preparation of symmetry enriched out of equilibrium phases as we will show.
Returning to the general topic, our discussion below thus focuses on the manifestations of symmetries in the dynamical evolution of quantum states. The approach starts from first principles with a representation of states via density operatorsρ({a i , a † i }) in Fock spaces spanned by creation and annihilation operators {a i , a † i }. Symmetry operations X are realized in ten families of unitary and anti-unitary transformations of these operators, nowadays mostly labeled by the Cartan symbols, A, AIII, AI,. . ., CI (see Appendix A for a review). These symmetries are well documented, meaning that the microscopically formulated starting point of the theory is under control. In particular, the ten-fold classification implies the existence of a unique symmetry label assigned to individual quantum states at each instance of time. For example, in the case of thermodynamic equilibrium, ρ eq ({a i , a † i }) = e −βĤ({ai,a † i }) , the state inherits the symmetries of the HamiltonianĤ. More generally, however, the classification statement remains formal before the state is described in the context of its dynamical evolution.
The dynamics of quantum states is described by evolution equations of the structure ∂ tρ (t) =Xρ(t). (1) This equation describes both Markovian state evolution, whereXρ(t) =X (t)ρ(t) is time local and non-Markovian cases whereXρ(t) = R + dsX (s)ρ(t − s) contains a retarded convolution over time. In either case, the symmetries of the linear dynamical generatorX and their manifestations in the solutionsρ are the central topic of this work. Our discussion will be general in that it treats interacting and non-interacting systems on the same footing. We first note that Fock space symmetry transformations X affect the operatorsρ({a i , a † i }) andX ({a i , a † i }) featuring in the equation asρ →ρ X andX →X X via their representation on the creation and annihilation operators. Importantly, the realization of a symmetry in the evolution equation may include an additional transformation of the dynamical parameter 'time', t, itself. A transformation is called a symmetry, if it leaves the evolution equation invariant, i.e. if the equation looks the same before and after the transformation. If this condition is met,ρ andρ X solve the same equation and are therefore identical; the symmetry of the dynamically evolving state is established, including in the limitρ(t → ∞), where it defines the stationary phases.
The moment, the dynamics includes elements of irreversibility induced by the coupling to an environment, the inversion of time t → −t in Eq. (1) is no longer physical -movies of irreversible processes don't make sense in reverse. The generic situation in this case is that of non-equilibrium dynamics, which will also be the focus of this paper. With inversion of time out of the picture, the invariance condition readsX X = +X . More specifically, the generators,X =X (Ĥ,D,P ) can generally be expressed in terms of three subordinate operators, describing the contributions of unitary evolution, dissipation and fluctuations to the dynamics. In this representation, the symmetry criterion splits into three, individually for these operators.
The symmetry classification of the generator of dynam-icsX defines the criteria required for the stabilization of a (stationary) state of definite symmetry. We already mentioned that the criteria obtained in this way differ from those of equilibrium systems. For example, the conditions for the Hamiltonian contribution to a dynamical evolution assume a form opposite to those in the equilibrium case. Given the scarcity of general principles characterizing out of equilibrium quantum distributions, the specification of symmetry criteria universally described in terms of the generatorsĤ,D,P is an important contribution of this work.
At this point, we have mentioned two settings, the limit of closed system unitary dynamics, and that of nonequilibrium irreversible dynamics. However, sandwiched between these two, we have a third major class, that of irreversible equilibrium dynamics, commonly associated to the physics of thermalization. Thermalization is irreversible, and as in the non-equilibrium case, a naive inversion of time in Eq. (1) is not physical. However, unlike in the non-equilibrium case, there still applies a principle of 'micro-reversibility'. In essence, it states that the rates of microscopic processes are determined by those of their time reversed inverse processes. Micro-reversibility implies a symmetry under shift inversion ∆t → −∆t + iβ, where β is inverse temperature, and ∆t the difference entering the correlation of observables at different times. In order for this condition to hold, fluctuations, P , and dissipation, D must be locked to each other via the fluctuation dissipation theorem Eq. (47) below. From the larger perspective of out of equilibrium dynamics, thermal equilibrium thus defines a 'fine tuned' case, much as unitary dynamics (D = P = 0) is an even stronger confined limit. Either limit comes with symmetry principles specific to its constraints, and different from the general case. The hierarchy of different settings is illustrated in Fig. 1 b).

A. Synopsis and summary of results
We now turn to a more concrete level and summarize the main findings of our work. Proceeding in a bottom up manner, we first discuss the realization of symmetries in fermionic Fock space, before turning to their representation in different descriptions of effective dynamics. We also compare our results to related work on symmetry classification of open quantum systems. Symmetries -Our starting point is the representation of symmetries as unitary or anti-unitary transformations in Fock space through their action on fermion operators {a i , a † i }. These operators represent the system after tracing over the environmental degrees of freedom, and no reference to a particular type of dynamics is made yet. Following the reasoning of Refs. [5,22], all we can say at this level is that modulo unitary equivalence ten different classes of transformations need to be distinguished.
More specifically, the basic operations from which these symmetries are built by composition are an anti-unitary transformation (T i T −1 = −i) acting as Ta i T −1 = u T,ij a j , where u T is a unitary matrix, and a unitary operation exchanging annihilators and creators Ca j C −1 = u C,ij a † j (see Eq. (4) for a more detailed representation). For bookkeeping purposes, we also define the combination S = T • C as an anti-linear (S i S −1 = − i ) operation exchanging annihilation and creation operators, Sa j S −1 = u S,ij a † j . Invariance of evolution equations -Consider Eq. (1), for the case whereρ =ρ({a i , a † i }) is the reduced density operator describing the state of an open quantum system, andX =X ({a i , a † i }) the dynamical generator governing its out of equilibrium evolution. The transformations X = C, T, S affecting the fermion operators as a i → Xa i X −1 ≡ (a i ) X define induced operationŝ X →X X andρ →ρ X . We understand an operation X as a physical symmetry if the transformed evolution equation ∂ tρX (t) =X XρX (t) remains invariant in the sense thatX In this case, the equation looks the same before and after the transformation, which impliesρ X (t) =ρ(t) for the solutions: the density operator does not change under the transformation and inherits the symmetry. Anti-unitary symmetry: equilibrium vs. non-equilibrium -Let us now discuss in more specific terms how the anti-unitary transformation T is represented in the evolution equation. Specifically, we will distinguish between the three different settings mentioned previously, unitary state evolution, equilibrium evolution, and general out-of equilibrium evolution.
In the first (textbook) case,X The generalization to the case of thermal equilibrium is complicated somewhat by the fact that non-Markovianity becomes essential for fermion systems. Referring for a detailed discussion to section V, the reason is that the frequency dependent Fermi-Dirac distribution coupling dissipation and fluctuation generators via Eq. (47) introduces retardation in the time evolution. Within this setting, the previously mentioned micro-reversibility principle manifests itself via the Kubo-Martin-Schwinger (KMS) relation [23,24] stating invariance of expectation values of any two-time correlators (with time difference ∆t) under an operation ∆t → −∆t − i/T , where T is temperature. This invariance motivates the definition of a generalized time reversal operation acting on functions of time as Note that like conventional time reversal E 2 β = 1 is an involutory operation. Combined with the Fock space symmetry T, it defines the extension of quantum mechanical time reversal to irreversible systems at thermal equilibrium. The ensuing thermal time reversal operation is compatible with the presence of a global time arrow in irreversible equilibrium dynamics. However, out of equilibrium dynamics excludes changes of the time variable. T now acts in autonomy as an anti-linear symmetry realized in Fock space. We are thus led to the conclusion that the watershed distinguishing between anti-linear symmetries with and without time inversion is the boundary between unitary or equilibrium dynamics on the one side, and non-equilibrium dynamics on the other.
What are the consequences of this finding? We first note that both, in and out of equilibrium the stationary states,ρ ≡ρ(t → ∞) satisfy identical symmetries: ρ =ρ X is uniquely fixed by the action of symmetries in Fock space (we will formulate the ensuing symmetry classes in more concreteness later in the text). This implies a high level of universality for a system's state, extending the ten-fold classification of equilibrium phases to the full realm of non-equilibrium stationary states. However, both, the dynamical generators, and the dynamical processes leading to stationarity satisfy opposite symmetry principles. For example, a Hamiltonian contribution,Ĥ, to a non-equilibrium generator must satisfŷ H T = −Ĥ to stabilize a T-symmetric stationary state, opposite to the equilibrium case. These differences are crucially important to the realization of stationary states by dissipative protocols, a point we will illustrate on the example of a topological phase in a chiral symmetry class. An equally important consequence is that in either case the language in which the respective symmetry criteria are articulated involves no more than the ten fundamental Fock space symmetries.
Below, we will discuss the above invariance principle for different realizations of dynamical evolution, including 'interacting'X s of quartic order, and non-Markovian generators required to include the case of thermal equilibrium. In either case, the dynamical generators contain three operators,Ĥ,D,P describing the effect of unitary state evolution, dissipation and quantum fluctuations, respectively. These operators are to the non-equilibrium system what the Hamiltonian is to an isolated quantum system. As we will discuss in detail, an immediate consequence of this statement is that individual (non-Hermitian) operators cannot define a phase, it takes the more structured information contained in the three op-eratorsĤ,D,P to do that.
We note that the representation of the anti-linear symmetry T as a pure Fock space symmetry (out of equilibrium), or in connection with time-reversal as T β (equilibrium) has consequences for the combined chiral transformation S = T • C as well. It does not affect, however, the charge conjugation transformation, which is unrelated to time altogether. Gaussian dynamics -While our approach works for general interacting systems, in the second part of the paper, we specialize to Gaussian evolutions, the irreversible generalization of free fermion systems. This setting allows us to present the general framework in more concrete terms, and it defines the non-equilibrium counterpart of the free fermion ground state classification [5]. The operatorX is now quadratic in a i , a † i , and we may turn to a first quantized representation in terms of matrices. The deterministic generator,K, is represented by a non-Hermitian matrix, K = H − iD, with Hermitian contribution H, and semi-positive Hermitian damping matrix, D, and the fluctuation generatorP by an anti-Hermitian matrix P , see Sec. III B for the precise relation between secondand first quantized representation. (Throughout, we use carets to label all operators in second quantized representation other than a i , a † i , while first quantized operators come without.) The symmetry criteria discussed previously now assume the form of matrix symmetries under transposition or complex conjugation, and give rise to a table with 40 entires, see table II: 10 symmetry conditions for Hamiltonians in equilibrium, and 3 × 10 for the Hamiltonian, dissipation, and fluctuation generators out of equilibrium, where the number 10 refers to the universal Cartan labels. Coming from a Hamiltonian perspective, the symmetry relations assume an unfamiliar form, which, again, has its origin in the absence of time inversion in the present setting.
For a given initial state, the generator of dynamics determines the system's state via a Gaussian density operator. Specifically, we represent the stationary long time limit (prior to normalization) as ρ ≡ exp(−Θ), where the Hermitian matrix Θ defines the effective Hamiltonian. In equilibrium, Θ = βH is determined by the Hamiltonian, implying that state and Hamiltonian transform identically under symmetry operations. Out of equilibrium, Θ = Θ(H, D, P ) is a function of the dynamical generator. However, the realization of symmetries on Θ is determined by the underlying Fock space operations and does not depend on whether the dynamical evolution was in or out of equilibrium (see the last column in table I.) For example, the state of a free fermion topological insulator with 'chiral symmetry' is described by a state anticommutative with an involutory matrix such as the Pauli matrix, σ 3 , i.e. σ 3 Θσ 3 = −Θ. In equilibrium Θ = βH and the Hamiltonian H obeys the same rule. However, if Θ defines an out of equilibrium distribution defined by a protocol with Hamiltonian participation, the latter must be commutative, σ 3 Hσ 3 = +H. The origin of this perhaps counterintuitive result can be traced back to Eq. (1): the presence or absence of an explicit time inversion in the realization of S = T • C in the equation accounts for the relative sign.
Irrespective of the different realizations of individual symmetries for the generator of Gaussian dynamics, a main conclusion for the Gaussian state classification is then that it leads to the definition of ten matrix symmetry classes, and that these classes are in one-to-one relation to the Cartan classes A, AI, AII, . . . , D defining the '10-fold way', cf. Appendix A for an overview. This correspondence becomes of key importance when we proceed to the classification of state topologies.
Topology -The objects entering the topological classification are macroscopic Slater determinants defined by the single particle eigenstates |γ of Gaussian density operators (or, equivalently, their effective Hamiltonians, Θ), where γ indicates the dependence on parameters such as a lattice momentum or spin-degrees of freedom. We consider partially filled systems and assume the existence of a subset of states occupied with high probability, p γ > 1 2 . These states assume a role analogous to that of quantum ground states, and deviations off p γ = 1 describe the residual effects of heating in a system with large excitation gap. Whether or not the parametric dependence of the states |γ admits the definition of a topological invariant depends on the symmetries ofρ, and on the dimensionality of the parameter space γ. A corollary of the Cartan classification of symmetries is that this information can be lifted from the periodic table of topological insulators and superconductors [6,7]. For example, a Θoperator chiral in the above sense, σ 3 Θσ 3 = −Θ, admits the definition of a state topology in odd but not in even parameter dimension, etc.
Notice that the effective Hamiltonian plays a double role in this context. It defines topological structures via the parametric dependence of its ground states, and controls the occupancy of these states via the probabilities p γ = f ( γ ), where Θ|γ = γ |γ and f ( ) = 1/(e + 1). The stable occupancy of the ground state of Θ requires the presence of a global purity gap ∆ p [25][26][27] defined in Eq. (22), avoiding totally mixed modes with occupation probability p γ = 1/2.
As a second condition, we require that the stationary stateρ = exp(−Θ) is attained with a finite rate, and that excitations out of it relax back in finite time. This permits slow variations of system parameters in time without leaving the instantaneous stationary state. The dynamical approach towards Θ is controlled by the deterministic generator, K = H − iD, and the minimal rate set by the spectral gap, ∆ s , which will be defined in Eq. (21). In passing, we note that the complex eigenvalue spectra of dissipation generators are interesting objects in their own right (see Refs. [20,28] for review). Their singularities in the complex plane, dubbed 'exceptional points', can be classified in terms of topological principles and leave signatures in specific dynamical response functions. However, the stationary density matrix implies a long time limit/integral over frequencies effectively averaging over these structures.
We will also address the formation of edge states at the boundaries between bulk topological phases. By definition, an edge is defined by the change in a topological invariant which in turn requires the closure of the purity gap [26] (unless the invariant is destroyed by violation of a symmetry condition). Ultimately, one would like to apply such edges as resources for the realization and manipulation of topologically protected edge states. In this regard, the closure of the purity gap appears to be bad news. By definition, a closing purity gap implies fully mixed configurations, e.g. Majorana edge qubits forced into an equal probability configuration of up and down states. However, as we will discuss, a way out of this dilemma is a simultaneous closure of purity and spectral gap at the edge. We will see how this happens, for example, in systems with Lindbladian state evolution where the topological state is the 'dark state' of the dynamics, and how this defines manipulable edge spaces.
Finally, we remark that topological classification generally requires more structural input than symmetry classification. For example, for a given fermion Hamiltonian -classifiable according to the tenfold symmetry scheme -the objects of topological classification can be zero temperature ground states [29] (as in the physics of topological insulators), unitary time evolution operators [30,31] (as relevant for the topology of Floquet systems), or group cohomology structures [32] (as relevant for the classification of interacting symmetry protected topological phases). The exploration of the full scope of topological classifications based on the non-equilibrium symmetry classification discussed here is a subject transcending the scope of this work.
Relation to other work in the field -We already mentioned that previous work on symmetry classifications emphasized the non-Hermiticity of non-interacting dissipatively damped dynamics in first quantized matrix representations (encoded inĤ andD, but discarding fluctu-ationsP ). Building on symmetry classifications of such non-Hermitian matrices [17], it has been reasoned that the Hermitian adjoint η : X → X † becomes a symmetry operation in addition to the standard unitary and anti-unitary symmetries of quantum mechanics [18][19][20]. The inclusion of this operation defined an extended classification, now containing 38 rather than the 10 classes governing the Hermitian generators of unitary time evolution [18,19]. However, causality arguments led Ref. [33] to the conclusion that only ten of these -defined as combinations of η, time reversal, T and charge conjugation, C -are physical in dissipative evolution.
Our analysis of symmetries is different from these works in that (i) it starts from a representation of symmetries in Fock space (as we will see, the passage from second quantized operators to matrices is not an innocent one), and (ii) system dynamics is fundamental to the symmetry classification. The latter requires consideration of all three generatorsĤ,D andP on equal footing, and on all levels of the description. For example, P contains the system's distribution function. If this information is not kept track of, the system may appear to be in different classes, depending on whether parts of it are integrated out or not. Erasure of this information may thus spoil the unambiguous assignment of systems to classes. (The sensitivity of time reversal to the realization of system partitions has also been noted in [34].) It is important to point out that the symmetry classes of stationary states depend on the symmetries of the above generators but not on those of initial states, unless multiple stationary states exist. For example, the initial state of a quantum system prepared with definite center of mass momentum breaks T. One may consider fine tuned protocols for which this breaking is preserved including under the evolution by T invariant generators. However, this situation is not generic. Rather, generic long time evolution will eventually eradicate the memory of the initial state and stabilize a T-symmetric phase. The disclaimer that we are focusing on universal symmetry classes is important inasmuch as non-equilibrium dynamics starting from specific initial states leaves ample room for dynamical or transient manifestations of symmetries and/or topology not captured by our analysis. As examples, we mention work on the classification of quantum quenches in closed systems governed by topologically non-trivial Hamiltonians [35][36][37], and on the structure of the complex spectra of dissipative generators, and topological signatures caused by the presence of singular 'exceptional points' in open system's state evolution [38][39][40][41][42][43][44][45][46][47][48][49][50], see [28] for review. Such structures show in observables probing the approach to stationarity, but not in the stationary phases themselves.
Plan of the paper -The focus in this paper is on general structures. However, for illustrative purposes we have included the discussion of one exemplary case study, namely the physics of an interacting model in class BDI (possessing a symmetry under both T and C), reducing to a variant of a Majorana chain in the Gaussian limit. This example is deliberately chosen to illustrate all concepts introduced below in the simplest possible scenario.
We will start in section II with a discussion of symmetries in Fock space. In section III we discuss the interplay of symmetries and topology in Markovian dynamics. This will introduce the material required to discuss the aforementioned case study in section IV. Section V goes beyond the Markovian limit and lifts the discussion of symmetries to the framework of the Keldysh path integral. This extension is required, e.g., to include the important limit of thermal equilibrium states. We conclude in section VI. Technical details are largely relegated to several Appendices. The inclusion of non-Markovian processes (required, e.g., to describe the relaxation towards a quantum thermal distribution) is beyond the scope of the Lindbladian description but can be addressed within the Keldysh path integral framework.

II. SYMMETRIES IN DRIVEN OPEN QUANTUM DYNAMICS
In this section, we start out from a precise definition of all symmetry operations X relevant to our discussion in fermionic Fock space. This part of the discussion makes no reference to the concrete realization ofX . In the remaining parts of the paper, we will then explore the manifestations of invariance for the different settings listed in table I.

A. Symmetry operations in Fock space
The symmetry operations relevant to our discussion are defined by U, T, C, S (for a review, see [29,51]), where U : The action of the symmetry operations on the creation operators is obtained by taking the adjoint of the above, e.g., Ta † i T −1 ≡ū Tij a † j , where here and throughout the overbar denotes complex conjugation.
Eq. (4) exhausts the list of unitary (U, C) and antiunitary (T, S) transformations of the operator algebra exchanging (C, S) or not (U, T) creators and annihilators. The transformations are defined to include an optional purely unitary operation, u X , where here and in the following X = T, C, S unless refined otherwise. Defined as it is, the list contains redundancy. For example, the combination of T and C is a transformation of type S = T • C, two different T, T combine to U = T • T , etc. However, due to the distinct physical meaning of the transformations, it pays to consider them separately.
Specifically, T is in the class of anti-unitary operations required to describe time reversal in unitarily evolving systems. Although the direct meaning of a 'time reversing' operation gets lost in irreversible dynamics as anticipated above, anti-unitary transformations continue to play an important role.
The 'charge conjugation' transformation, C, exchanges the role of particle and holes via a unitary operation in Fock space, C i C −1 = +i [52]. As mentioned above, the definition of S = T • C is technically redundant, but considered here for its role in the description of 'particlehole' exchange, and the ensuing 'chiral' symmetries.
The group, U, includes the familiar unitary symmetries such as number conservation, U = exp(iαN ), (2), (σ, σ =↑, ↓ are spin indices), etc. In the definition of symmetry classes, the complementary set of operations T, C, S must be considered in relation to the unitaries, U. More precisely, these operations define meaningful symmetry classes only if they commute with the unitary symmetries of a system, and in the consequence act within the irreducible Fock subspaces defined by them -sectors of conserved angular momentum, particle number, lattice symmetry, etc. For a more detailed discussion we refer to Refs. [22,29,51], and to Appendix B, where we discuss this point on a few illustrative examples.
Finally, we note that all operations introduced above are compatible with operator products in that X(ÔÔ )X −1 = (ÔÔ ) X =Ô XÔ X . This feature will become important when applying transformations to evolution equations such as First quantized representation -To efficiently describe the above transformations in the language of matrices, we define the Nambu operators A i ≡ (a i , a † i ) T . The operations as defined in Eq. (4) act on these opera- For later reference, we note the relations where here and in the following the Pauli matrices σ i act in Nambu space. A general bilinear free fermion operator O has the Nambu representation defined in terms of a first quantized matrix operator O = {O ij }. SubjectingÔ to the above transformations aŝ (where only (TO ij T −1 ) =Ō ij acts non-trivially on matrix elements), we obtain the first quantized version of the symmetry operations (for the relation in the third row, see also App. B, Eq. (B1)), We say that an operatorÔ has X-symmetry, ifÔ X =Ô, or O X = O in first quantization. The operations defined by Eq. (9) will be the basis for the definition of matrix symmetry classes. Symmetry tabulation -Referring to Refs. [5,22] for a more comprehensive discussion, we now briefly identify the symmetry classes defined by application of T, C, and S = C • T within the sectors of definite unitary symmetry: two-fold application of either X = C, T to an X-symmetric operator leads to The commutativity of V X with all operators O in the representation space [53] requires V X = ±1 N , and these two options define the symmetry classes X = ±1. The absence of X-symmetry is called X = 0. Similarly, the presence/absence of symmetry under S is defined as S = 0 or S = 1 [54]. Counting all options, we obtain the list of ten symmetry classes tabulated in Appendix A for the convenience of the reader.
Already at this stage it is evident that, irrespective of the specific type of the irreversible dynamics (Markovian/non-Markovian and equilibrium/nonequilibrium) we have a maximum number of ten symmetry classes. However, as we will demonstrate below, the conditions under which a given class emerges depend on whether equilibrium or non-equilibrium dynamics is considered.

III. SYMMETRIES IN MARKOVIAN DYNAMICS
In this section we discuss how the symmetry principles introduced above materialize in the case of Markovian dynamics, i.e. that of systems coupled to a bath with vanishingly short memory. We will start out from a representation of Markovian state evolution in the form of a general Lindblad equation, then specialize to Gaussian state evolution, and finally discuss topological structures.

A. Symmetries in Lindbladian dynamics
To study general Fock space transformations introduced above in the case of Markovian dynamics, we consider a quantum master equation of Lindblad form [55][56][57][58][59] Here,L α are the 'jump operators' coupling the system to a bath, and the sum extends over distinct types of coupling. Typical realizations include jump operatorŝ L α ∼ a i linear in fermion operators, or number conserving couplingsL α ∼ a † i a j . Number conserving systembath couplings imply quartic (or even higher order) Lindblad operators,L † αLα = O(a 4 ), and this is why we need to include nonlinearities even if the focus will ultimately be on Gaussian states.
We now apply the strategy Eq. (5) to the Lindblad equation, withρ →ρ X and (iĤ) Under these conditionsρ andρ X satisfy the same irreversible evolution equation and are thus identical,ρ =ρ X . (In cases where the stationary state depends on the initial state, symmetry of the latter becomes an additional condition [60].) At this level, the analysis includes interacting settings. For illustration, see the case study of section IV, where we consider an interacting one-dimensional chain with BDI-symmetry, (T, C, S) = (+1, +1, 1).

B. Symmetries of Gaussian states
While non-linear Lindbladian equations are as complex as their reversible (von Neumann) cousins, one often has situations where the Hamiltonian is quadratic, i.e.Ĥ = 1 2 A †T HA, and Hartree-Fock mean field approximations may be applied to non-linear operatorsL † αLα to define quadratic approximations (cf. Refs. [25,26,61]). Under this condition [62][63][64], we have a reduction [65] αL † with a semi-positive and Hermitian matrix M , and Hermitian H. These conditions reflect the complete positivity of the Lindblad generator [55,56], and the dynamical conservation of Hermiticity ofρ. The Lindblad equation then assumes the form (similarly to Eq. (7), . Due to fermion exchange symmetry, the kernel, O, defining any bilinear formÔ = 1 2 A †T OA satisfies O T = −σ x Oσ x (cf. Eq. (B1)). We here have assumed without loss of generality that O is traceless. This feature suggests a decomposition   from ψ|M |ψ ≥ 0 and ψσ x |M |σ x ψ ≥ 0 for any state |ψ ). With these definitions we have In addition, the conservation of Hermiticity of the density matrixρ by the Lindblad generator in each time step implies These three building blocks fully describe the dynamical generator, and they all have individual physical meaning: H describes the reversible contribution, D the dissipative damping generator, and P generates fluctuations.
Eq. (B1) applied to the bilinear forms Eq. (11) implies a reductionM = −iP in the anti-commutator term of the Lindblad equation. However, due to the presence of ρ the full content of M remains in the jump term: with the positive and real constant c Proceeding as before, we subject Eq. (16) to the operations in Eq. (4), and deduce the symmetry relations the generator compounds H, D, P must obey to obtain form invariance of the evolution equation, and symmetrŷ ρ =ρ X of the density operator. Referring for the details of the straightforward calculation to Appendix C, this leads to the transformation criteria summarized in table I.
Notice the similarities and differences to the symmetry classification in unitary dynamics. Similarly to that case, the ten different combinations are applied to operators which are individually Hermitian (H, D), or i-times Hermitian (P ). However, unlike with unitary state evolutionρ → e −iĤtρ e iĤt , where just one Hermitian operator determines the symmetry of the classification object, we here have a situation, where the three operators, H, D, P , are essential to the description ofρ, or better to say of the Hermitian effective Hamiltonian Θ governing the Gaussian statesρ = exp(− 1 2 A †T ΘA). In view of the fact that the focus in the literature (cf. Refs. [8-16, 33, 66]) often is on the dissipative damping operator K ≡ H − iD, it is worth pointing out that the symmetry of all three, H, D, P is essential to the invariance of the state,ρ.
Second, we note that due to the lack of the extraneous time inversion t → −t implied in the definition of unitary time reversal symmetry, the symmetries T and T • C = S are realized differently than in Hamiltonian (or equilibrium) dynamics. For example, T-symmetry of the density operator requires H = −U TH U † T , sign different from the standard time reversal condition of quantum mechanics. Similarly, the 'chiral symmetry', S, requires commutativity H = +U T S HŪ S , again with opposite sign. In Sec. V B, we will compare this symmetry requirement to that in the case of unitary evolution.
Finally, notice the absence of a Hermitian adjoint η : O → O † applied to non-Hermitian operators in the present framework. To see why, notice that our first quantized T acts as T : (9)). By contrast, Ref. [33] T . The two transformations differ by a relative application of η, T = Tη (up to unitaries). The T transformation was motivated by a criterion of causality -namely to avoid a sign change in the imaginary contribution to the damping generator K = H − iD -and by requiring a smooth connection of the irreversible dynamics to purely unitary evolution.
In our approach, the causality/sign issue does not arise since the transformation t → −t is avoided. This is seen by noting that K couples to the dynamics (symbolically) as exp(−iKt) = exp((−iH − D)t). The absence of the t → −t operation in our understanding of symmetries in irreversible non-equilibrium dynamics makes η never appear in any symmetry operation.
Rather than probing Θ directly, we here consider the stationary covariance matrix, Γ = lim t→∞Γ (t) [62][63][64], which carries the same information but is more directly accessible. The covariance matrix is a 2N ×2N Hermitian matrix defined asΓ where a is a composite index comprising the Hilbert space index, i, and the two-component Nambu index. Inspection of Γ in the eigenbasis shows that i.e. the covariance matrix and the effective Hamiltonian carry identical information, and in particular share the same ground state. Substitution of Eq. (17) into Eq. (10) shows that Γ is obtained as the stationary solution of the equation K ≡ H − iD. The long time limit of the solution is obtained as Eq. (20) reveals much of the physics of the covariance matrix, and the dynamical processes stabilizing it. Specifically, the equation shows that the stationary state is obtained by retarded/advanced propagation of the fluctuation matrix P , as described by the retarded and advanced Green functions (ω − K) −1 and (ω − K † ) −1 . The non-negative matrix D contained in K = H − iD defines the relaxation rate at which the stationary state is attained. We require a finite minimal rate, as defined by the spectral gap, This requirement permits slow (adiabatic) changes of system parameters in time such that the system stays at all times in the instantaneous steady state. Note that only the imaginary part of the eigenvalues, i.e. the lifetime of the slowest decay rate in the system, enters the damping gap; thus the Hamiltonian alone may be zero, or have zero modes without corrupting an open damping gap. Finally, it is straightforward to check either by inspection of the representation Eq. (20), or of the evolution equation Eq. (17), that the symmetries of K, H, P listed in table II extrapolate to symmetries of Γ and Θ as indicated in the last column.
On this basis, we now discuss the topological structures defined by Θ. This operator plays a role analogous to the gapped single particle Hamiltonians describing topological insulators or superconductors. The dissipative analog of the single particle gap in these systems is the purity gap [26,27], which can be defined via the eigenvalues of Γ = tanh(Θ/2) as The purity gap sets a lower bound for the negative part of Θ's spectrum relative to 0. In terms of Θ, a vanishing purity gap is similar to a 'metal', where a chemical potential intersecting a band invalidates the definition of a topological ground state. The occupation probability of an Θ-eigenstate ψ γ is given by the Fermi function p γ = f ( γ ) = 1/(1 + e γ ). This occupation balance demonstrates that negative eigenvalues of finite minimal modulus are required to stabilize a fully occupied ground state in the thermodynamic limit containing a macroscopically large number of states. In contrast, a state with vanishing purity gap features at least one mode with γ * = 0, representing a totally mixed fermion state with occupation probability p γ * = 1/2. For systems with finite spectral and purity gap, the ground state {ψ γ } of Θ defined by the Slater determinant of all negative eigenvalue Θ eigenstates becomes the subject of topological classification. Since the symmetries of Θ are realized identically to those describing Hermitian Hamiltonians, the classification of dissipative topological phases becomes equivalent to that described by the periodic table of topological insulators (cf. table III). The same goes for physical information obtained via topological principles from the periodic table. For example, a two-dimensional dissipative Chern insulator in class A ((T, C, S) = (0, 0, 0)) [61,[76][77][78] supports circulating chiral edge states, or a dissipative quantum wire in class BDI ((T, C, S) = (+1, +1, 1)), has Majorana end states [25], etc.

D. Edge state formation
While the classification of bulk topological phases requires the presence of a purity gap, at the boundary this gap closes. In fact, one may take the closing of the purity gap as a definition of the phase boundary. This interpretation follows from the identification of the covariance matrix with a band 'Hamiltonian', and its ground state as the carrier of a topological index. Changes in the ground state require the closure of a gap, presently the purity gap in the spectrum of Γ. Edge states are the low lying eigenstates of Γ spatially confined to that boundary region.
While the existence of edge states is a robust feature granted by topology, the accessibility and manipulability of these states -a feature required by, e.g., quantum information applications -is another matter. In fact, we here run into the seemingly paradoxical situation, where the closure of the gap around zero eigenvalue in Γ, i.e. the definition of the edge appears to contradict the accessibility of these states. To see how, consider the density matrixρ = exp(− 1 2 A †T ΘA) projected to the edge state subspace. Assuming Θ to be diagonalized, Θ = diag({ γ }) individual states |γ ≡ a † γ |0 are occupied with probability f ( γ ). For γ 0, the Fermi function approaches the value 1/2, which means that the state of the edge is given by a mixed state defined by equally occupied and empty edge states. For example, a Majorana wire whose two-Majorana edge space would be in an equal weight mixture of its two states, and hence useless for 'qubit' applications.
However, there is a loophole in the argument. It presumes that the spectral gap remains open at the boundary. To see why this matters, consider the solution of the differential equation (19) governing the evolution of the covariance matrix, with O(t) ≡ e −iKt Oe iK † t , and initial state Γ 0 . If the spectral gap, defined as in Eq. (21), is finite, the stationary state Γ of Eq. (20) is approached exponentially fast. However, now consider the different situation, where the approach of the boundary goes along with a simultaneous vanishing of all three generators of the dynamics, H, D, P within the subspace in which the purity gap closes. In this case, the solution of the evolution equation projected to that space will retain information on the initial state, and we have an edge capable of storing information. The condition of simultaneous vanishing of all generators is not quite as stringent as it may seem. For example, it is realized in dissipative systems possessing a dark state or multiple of these forming a dark space. A dark state is a state ρ = |Ψ Ψ| stationary under the full generator of the Lindbladian dynamics Eq. (12). This condition requires the simultaneous vanishing of all partial generators within the dark space. For an example and the discussion of the ensuing edge space, see section IV D.

IV. CASE STUDY
In this section, we consider a one-dimensional lattice model in the non-equilibrium symmetry class BDI (T, C, S) = (+1, +1, 1) for H, D, and P to illustrate the general concepts introduced above, and connect to various other concepts currently discussed in the literature. The model is defined by a chain of L sites i containing two orbitals 1, 2 indicated by red dots in Fig. 2. For the moment, we assume periodic boundary conditions, a system with edges will be considered below. The irreversible dynamics is governed by number-conserving jump operators, so that the Lindbladian is quartic in fermion operators. We consider a half filled system in which a Hartree-Fock style linearization of the evolution around a macroscopically filled state is possible, following the construction principles of [25,61]. The ensuing mean field dynamics is controlled by a parameter ϑ such that for ϑ = 0 the stationary state is a product state of decoupled equal weight superposition (spin-x) states defined along the solid links in the figure. In the opposite extreme, ϑ = π/2, it is a spin-x state defined via the dashed lines. For generic values, the eigenstates of the effective Hamiltonian Θ are spatially extended, and the full ground state is characterized by a winding number which changes in a topological phase transition at ϑ = π/4. In the following, we discuss both purely dissipative protocols stabilizing this state, and generalizations including an added Hamiltonian.

A. Interacting model
Let a 1,i and a 2,i be annihilation operators for orbitals 1 and 2 at site i, and define the two-component operator a i = (a 1,i , a 2,i ) T . Now consider the set of transformed two-component operators defined aŝ acting on neighboring lattice sites.l i andr i are related to the original operators by unitary transformationsl i = V L a i andr i = V R a i , with V L = e i π 4 Σy (E 11τ + E 22 ) and V R = e i π 4 Σy (E 11 + E 22τ ), such that {l a,i } and likewise {r a,i } (a = 1, 2 denoting the orbital index) generate a fermion algebra as well. Here (E ij ) kl = δ i k δ j l are matrices in orbital space, Σ i are Pauli matrices in the same orbital space, andτ is a lattice translation operatorτ a i = a i+1 .
To get some intuition for the operatorsl, consider the N particle configuration defined by occupation of alll 2 states: |L ≡ il † 2,i |0 . Since we have 2N sites in total, this is a half filled state. Withl † 2,i = 1 √ 2 (−a † 1,i+1 + a † 2,i ) we observe that it is a product state defined by an equal weight superposition (or spin x-state) defined along the solid bonds in the figure. Conversely, |R ≡ ir † 2,i |0 is a product state of spin-x hybridizations along the dashed bonds. Second the states |L , |R afford an interpretation of free fermion ground states of distinct topological order. To see how, consider the operatorsĤ L ≡ il † T i Σ zl i . The state |L is the gapped ground state of this Hamiltonian at half filling. On the other hand, substitution of the definition (24) shows thatĤ L = q a †T q w † w a q , with w = e iq in momentum space. This comparison identifies |L as the ground state of a gapped chiral parent Hamiltonian, where the winding of the phase z = exp(iq) as a function of q ∈ [0, 2π) defines a topological invariant W = 1 2πi dqz −1 ∂ q z = 1. Likewise, the state |R is the ground state of a HamiltonianĤ R with winding W = −1.
In the following, we construct dissipative protocols driving into the states |L , |R , or more generally a competition between them. Depending on the balance at which these reference states enter we expect different topological phases, with a transition between them. Specifically, it will be instructive to consider two different Lindblad protocols (c = cos ϑ, s = sin ϑ), whereL x is the Lindblad operator defined relative to the fermion operator choicex, (i.e.l,r, or cl + sr respectively). TheL x are bilinears in jump operatorsL as in Eq. (10). Suppressing the site index we define the latter asL where = 1, 2,x =l,r, orx = cl + sr for the second protocol, and κ setting the coupling strength. Note that these operators are number conserving and act locally in real space. Specifically, the jump operatorL a1 depletes the upper x-band andL a2 fills the lower, so that the joint action of the two fills the ground state associated to thê x-operators. The jump operatorsL an , a, n = 1, 2, define a Lindblad equation Eq. (10) withĤ = 0 and the sum extending over the four configurations and lattice sites, α ≡ (a, n, i).
The 'incoherent' operatorL inc puts the two operators cooling into |L and |R , respectively, into competition, while the operatorL coh cools into the ground or 'dark' state defined by a coherent superposition ofl andr. Either way, we expect a phase transition at ϑ = π/4, where c = s = 1 √ 2 and the left and the right operator algebra couple at equal strength.

B. Symmetries
Turning to symmetries, we first note that the Lindblad operators in theL x model are T = 1 symmetric with auxiliary unitary U T = 1 in the sense of our earlier discussion. The reason is thatV L,R = V L,R are real, and hence TL an T −1 =L an for real κ, c, s.
Second, we have a S = 1 symmetry under the chiral symmetry operation with auxiliary matrix U S = Σ z . Noting that the matrices V L,R in Eq. (24) satisfy the symmetry it is straightforward to verify that SL an S −1 = (−1) a+1L an , wheren = 2, 1 for n = 1, 2. Since theL's appear pairwise and are summed over, the sign factor drops out, and L x has a chiral symmetry. Finally, the composition T • S induces a C = +1 symmetry under the unitary map Ca i C −1 ≡ Σ z a † i , with U C = Σ z . We thus have a model with BDI symmetry (T, C, S) = (+1, +1, 1) on the second quantized level. C. Linearized model L x is quartic in fermion operators and hence defines a nonlinear problem. However, this nonlinearity can be Hartree-Fock decoupled, thanks to the macroscopically large number of particles in the problem. To see how, consider the bilinears (no sum convention here) 1x 1 at half filling a ai a † ai = 1/2, used in the final step. Likewise,L † a2L a2 → κ 2x2x † 2 . Withx = V a, and V the linear transformation defining the site operatorsx via a (e.g. V = V R in Eq. (24) forx =r), we thus find where A = (a, a † ) T as before. Comparison with Eq. (13) shows that the principal building blocks of the dynamical generator are given by where we rescaled κ → 4κ for notational simplicity. Since there are no terms coupling the particle and hole sector in Nambu space, we focus on the particle blocks throughout, and redefine D → κV † V , and P → iκV † Σ z V for simplicity. Using Eq. (27), we note that in this first quantized representation, the chiral symmetry acts as D → D S = −Σ z DΣ z = −D, and P → P S = +Σ z P Σ z = −P . This is in accordance with table II and implies that D and P are diagonal and off-diagonal in orbital space, respectively. Time reversal manifests itself in the reality of the matrices D and iP . In passing, we note that this reality can be used to represent the system in the Majorana representation of the BDI chain [79]. For example, with η = 1 √ 2 (a † + a), ν = 1 i √ 2 (a † − a), it is straightforward to verify that a †T Da = −iη T Dν. However, the representation change will not play an essential role throughout.
We now have everything together to explore the stationary states stabilized by this dynamics. As a warmup exercise, consider the casex =l. The corresponding matrices V L are unitary, and a quick calculation using Eq. (24) yields the momentum representation with matrix structure in orbital space. Substitution of these expressions into Eq. (20) defines the momentum representation of the covariance matrix We thus conclude that the system cools at uniform rate (D = const.1) into the |L state with its winding number W = 1. The covariance matrix (and, equivalently, the effective Hamiltonian, Θ) are chiral in the sense of table II, such that we have a dissipatively realized class BDI system with stable topological phase. Next consider the incoherent driving protocol. Since the Lindblad generators are linear in the matrices P and D, we now have a linear superposition, D = κ(c + s)1, and P = iκ(E 12 (ct † + st) + h.c.). Substitution of these expressions into Eq. (20) gives The straightforward computation of the winding number shows that W = sgn(c − 1 Now, compare this to the driveL x defined by the linear superpositionx = cl + sr. Up to a normalization factor, the algebra ofx-operators satisfies canonical commutation relations, {x a,q ,x † a ,q } = δ aa δ qq (1 + 2sc cos q), implying that we cool into a pure state. Proceeding as above, we define a matrix V = cV L + sV R , which is 'almost' unitary, V † V = (1 + sc(t + t † ))1, reflecting the scaled commutation relations. Using Eq. (30), we find D = κ(1+2sc cos q), and P = iκ(E 12 (ce −iq +s) 2 e iq +h.c.), which gives Γ coh ≡z coh E 12 + z coh E 21 , z coh = c 2 e iq + s 2 e −iq + 2sc 1 + 2sc cos q .
Again, we have a topological phase transition at c = s = 1 √ 2 . However, its critical signatures are markedly different from those of the incoherent case. The covariance matrix now has unit-norm eigenvalues, |z coh,q | = 1, reflecting the purity of the stationary state. Interpreted as a topological band insulator Hamiltonian, it thus resembles a Hamiltonian with 'flattened' non-dispersive spectrum. The corresponding purity gap equals unity, except at the phase transition point, where the absence of a well defined limit in z coh,q as q → ±π reflects the singularity required to have a unit norm curve z coh,q change its winding number around the origin. However, the damping gap, generally open in the incoherent case, now does close at criticality.

D. Edge states
Both, in the coherent and the incoherent protocol, an edge is defined by the closure of the purity gap. In the incoherent protocol, the spectral gap remains open, and in the light of the discussion of section III D this implies a non-manipulable edge in the sense that the edge covariance matrix will approach Γ = 0 exponentially fast, reflecting a fully mixed state. However, in the coherent protocol the situation is more interesting. Here, too, a (now singular) vanishing of the purity gap defines the edge. However, this singularity goes along with a vanishing of both D and P . The situation can be described both, in a continuum representation via a smooth variation of the parameters (c, s) through the critical point, or directly on the lattice. We here choose the latter option, and identify a simultaneous null-space of the operators D and iP in real space, for a system cut open as indicated in Fig. 2. To this end, we first identify a vector |Ψ 0 annihilated by V . With V = cV L + sV R and V L,R given by Eq.  (30) implies that |Ψ 0 is a simultaneous eigenstate of P and D.
For unitary systems whose Hamiltonian has the same topology as Γ coh , these states are the Majorana edge states of the BDI chain (the limit c = 1 corresponding to its 'sweet spot', where the decoupling of an edge Majorana is manifest in the lattice representation). Presently, the same topological setting manifests in a dark space, spanned by |Ψ 0 Ψ 0 |, |Ψ 1 Ψ 1 |, where the many body states |Ψ 0,1 differ in the occupation of the complex fermion corresponding to the two Majoranas centered at the left and right edge, respectively. The decoupling of |Ψ 0 from the dynamics implies that arbitrary mixed statesρ = i=0,1 p i |Ψ i Ψ i |, p 0 + p 1 = 1 are solutions of the dynamical equations and we have a freely manipulable edge space.

E. Adding a reversible contribution
While our so far discussion focused on purely dissipative dynamics, it is straightforward to include a Hamiltonian contribution H q = hĥ q · Σ, where the unit vectorĥ q describes the action of the Hamiltonian in orbital space, and h is its strength (here assumed uniform for simplicity.) For definiteness, we study the influence of this Hamiltonian on a purely left winding dissipative background,x =l. Parameterizing the dissipative generators (31) associated to the Lindblad operators in Eq. (24) in the same way, D q = κ1, and P q = iκp q · Σ, wherê p q = (cos q, sin q, 0) T , it is straightforward to verify that the stationary covariance matrix assumes the form This solution demonstrates the influence of the Hamiltonian on the stationary state. In general, only Hamiltonians obeying the symmetry condition of table II, H = +Σ z HΣ z , orĥ x =ĥ y = 0 define a Γ-matrix chiral in the sense Σ z ΓΣ z = −Γ. In other words, it takes a Hamiltonian diagonal in the orbital space to preserve the off-diagonality of Γ. (One may trace the origin of this perhaps unexpected finding back to the fact that the anti-unitary irreversible chiral symmetry operation S = T • C does not involve a sign change of physical time.) If in contrast H would obey the equilibrium relation H = −Σ z HΣ z , chiral symmetry of the steady state would in general be lost. We also note that a symmetry preserving H does not commute with P (the Pauli matrix structure again) and hence degrades the purity of the state: for non-vanishing coupling h, the eigenvalues of Γ no longer have unit modulus. While the purity gap gets affected, the spectral gap, κ, remains untouched, as long as [D, H] = 0, which the setup above assumes. Finally, there is one interesting exception to the general rule above: in cases, where the dynamics stabilizes a dark state, the symmetry condition on H gets lifted. A dark state is a zero eigenstate of the full Lindblad generator. Translated to the language of the matrices, H, P, D, this requires commutativity of all three of them, and P 2 = D 2 . While the latter has been a feature built into our model from the outset, the commutativity [H, P ] = 0 impliesĥ q ×p q = 0, and hence the vanishing of the second term in Eq. (35). In this way, the Hamiltonian decouples, and we are left with the pure and chiral configuration described by the first term, i.e. the projector onto the dark state.

F. Non-Hermitian SSH model and exceptional points
Against this background, let us address a few generalizations of Hermitian one-dimensional systems with chiral symmetry [80], which have been discussed in the recent literature. Specifically, the non-Hermitian SSH model [39,42,49] is defined by the matrix where z ± = z ± ∆z are arbitrary complex numbers. Can this model feature as a generator of dissipative dynamics, be associated to a symmetry class, or define a fermionic topological state? Turning to the first part of the question, we need to understand how H nh relates to the general triptych H, D, P of unitary, dissipation and fluctuation generators, respectively. Since the imaginary parts of the eigenvalues λ ± = ± √ z + z − have indefinite sign, the direct identification H nh ? = H − iD does not define a legitimate dissipation generator. As a cure, one could add an overall unit matrix −iκ 1, with κ ≥ |Imλ ± | to make it sign definite, such case is considered below. Here, we discuss the alternative interpretation H nh = H 1 +i H 2 = H+P as the sum of a Hermitian part H 1 ≡ H and a fluctuation generator H 2 ≡ −iP . Choosing ∆z ≡ iκ for definiteness, a quick calculation shows that the fluctuation generator iH 2 = P becomes identical to the matrix P of the winding protocol Eq. (31). However, the definition of a consistent Lindbladian dynamics requires the balancing presence of a dissipation generator D of matching strength such as D = κ1. We are led to the conclusion that H nh by itself does not define a consistent dynamical protocol, while the generator defined by the matrices H = H 1 , D = κ1 and P = iH 2 does. Referring back to the discussion of the previous section, this generator stabilizes the covariance matrix (35). However, this matrix, and the corresponding effective Hamiltonian, Θ, do not have chiral symmetry since H 1 does not fulfill the non-equilibrium chirality condition of table II. Hence Θ does not possess topological ground states. Although H nh 'looks chiral', it does not define a fermionic topological state of Lindbladian dynamics.
Role of exceptional points -Exceptional points are terminal points of branch cut singularities forming in the complex spectra of non-Hermitian linear operators [81]. At these points eigenvalues merge, and topological numbers counting the multivaluedness of the corresponding eigenvalue can be determined. A prototypical setup of this sort is defined by a non-Hermitian matrix N = N 1 +iN 2 with Pauli matrix structure and Hermitian and anti-Hermitian contributions, N 1 = n 1,µ Σ µ , N 2 = n 2,µ Σ µ , where µ = 0, 1, 2, 3, Σ 0 = 1, and real coefficients n 1,µ , n 2,µ . The semi-positivity constraint reads n 2,0 ≥ | n 2 | ≥ 0. Its eigenvalues λ ± = n 1,0 + in 2,0 ± ( n 1 + i n 2 ) 2 , merge at the terminal point of the branch cut defined by the vanishing of the argument of the square root. The classification of exceptional points according to the symmetries of their parent operators, N , has attracted a lot of attention in the recent literature [38][39][40][41][42][43][44][45][46][47][48][49]. Do exceptional points affect the phases forming as stationary limits of driving protocols? To answer this question, we need to take a look at non-Hermitian linear combinations of the three constituent operators (H, D, P ). Specifically, we have seen that K = H − iD governs the approach towards the stationary limit. It is straightforward to realize exceptional points by adaption of the chiral model dynamics studied above. As an example, consider the matrices V L in (24) modified as V L → V L (E 11 + rE 22 ) with real r. Inspection of Eq. (30) shows that this changes the matrix Adding to this the Hamiltonian contribution H = H 1 of the SSH model Eq. (36), we obtain a model with dissipation generator K = H − iD, whose eigenvalues are given by Both eigenvalues have negative imaginary part, a necessary condition to define a valid dissipation generator. For |z| = (1−r 2 ) 1/2 κ/2 an exceptional point with degenerate eigenvalues is realized. Does this non-analyticity affect the nature of the system's stationary phases? To understand what is happening, consider the general formula for the covariance matrix Eq. (20), represented as a frequency integral over the retarded and advanced propagators G R ≡ (ω − K) −1 and G A ≡ (ω − K † ) −1 . With Imλ ± < 0 in the lower complex plane, G R/A is analytic in the upper/lower half of the complex plane, with simple poles at ω = λ ± and ω =λ ± , respectively. At first sight, it looks like the residue of the pole integration might depend on the values ω = λ ± with their non-analytic dependence on system parameters. However, actually doing the integrals one finds that the resulting expression for the covariance matrix is a rational function of these parameters.
In other words, the presence of exceptional points has no effect on the stationary long-time phase. While we have no proof demonstrating this feature in the most general terms, it is straightforward to verify for translationally invariant dynamical generators with two internal bands, and we suspect it to be of general nature. Non-analyticities in the complex eigenvalue spectra of dissipation generators may affect the transient dynamical stages on the way towards the stationary limit or the dynamical response of the stationary state [8][9][10][11][12][13][14][15][16]66]. The physical significance of exceptional points hosted in the non-Hermitian matrix K for the dynamical evolution is clearly read off Eq. (23), in the perhaps clearest way for a system that approaches an infinite temperature state (P = 0). However, they do not appear to affect the ensuing stationary phases and their topological classification in themselves.

V. BEYOND THE MARKOVIAN LIMIT
Our discussion so far focused on the Markovian case of a memoryless environment. This excludes important settings, notably those stabilizing quantum thermal distributions of fermion systems. The extension to non-Markovian situations is achieved via the Keldysh path integral formalism [82,83]. In this section, we will discuss how the apparatus of symmetries manifests itself in the path integral, and then apply it to situations outside the Markovian limit. The main subject of this section is the extension of the time reversal transformation of quantum mechanics, Fock space T combined with time inversion t → −t, to irreversible equilibrium dynamics. The ensuing thermal time reversal transformation, Eq. (3) is essential to the identification of the symmetries ofĤ,D,P in equilibrium settings. Naturally, the symmetries ofĤ coincide with those of the standard zero temperature setting familiar from the literature. However, they differ from those found for non-equilibrium dynamics, table II.

A. Symmetries in the Keldysh path integral
Within the fermion Keldysh framework, physical observables are computed as expectation values of a coherent sate path integral (see Refs. [84][85][86] for review). Specifically, the covariance matrix now assumes the form where the kernel may have time, M = M(t−t ), or equivalently frequency dependence. The functional integration ≡ D(η, ν) is over two Grassmann variables η = (η c , η q ) T , ν = (ν c , ν q ) T where each component η c,q = (η c,q 1 , η c,q 2 ) T is subject to a Nambu doubling, and carries a Hilbert space index. The Pauli matrices τ i act in c/q space. η c/q are often referred to as 'classical' (c) and 'quantum' (q) components of the fields. We note that complex conjugation has no meaning for these integration variables. However, with the Fourier convention ((dω) = dω 2π ) we have ξ ω = ξ −ω , as for real variables. The Keldysh path integral contains the information previously expressed in terms of the covariance matrix. To see how, we use the rules of Gaussian integration to obtain For frequency independent P ω = P , K ω = K, this reduces to the covariance matrix Eq. (20) of the Markovian framework. However, the advantage of the Keldysh path integral is that it allows us to go beyond that, and include processes with memory. Referring for a more detailed discussion to section V B, the most important representative of this category is the thermal Fermi-Dirac distribution. In this case, the damping and fluctuation kernels are related by the fluctuation dissipation relation Eq. (47), which implies that the process remains non-Markovian (P cannot be approximated by a constant) at all frequencies, ω.
The practical identification of symmetries within the path integral formalism differs somewhat from our previous strategy. There, we had asked what symmetry transformations leave equations of motion invariant. Presently, it is more natural to ask what transformations leave the path integral action unchanged. Since equations of motion are implied by the path integral (even though they may assume an effectively intractable form in non-Markovian situations), the two strategies are equivalent. Within the present approach, we will use various freedoms in manipulating the integral, notably the option to exchange the order of variables. We here illustrate these features on a simple consistency check, namely the path integral verification of the Hermiticity of Γ. This introduces the manipulations required in the subsequent discussion of anti-unitary symmetries. Expressed in path integral language, In the second line, we took the transpose in the action (catching a minus sign due to the anticommutativity of the Grassmann fields), and in the third used M † ω = −τ 3 M ω τ 3 and a variable transform η → τ 3 η, ν → −τ 3 ν , to arrive at an expression identical to the original integral, except differently named integration variables.
Similar manipulations applied to the covariance matrix transformed as indicated in the last column of table II (for details, cf. Appendix D) lead to the condition (U X act as unit matrices in Keldysh space) Substitution of Eq. (39) into these relations yields conditions for the blocks, P , K = H − iD identical to those listed in table II. However, recall that for a matrix kernel with time dependence, complex conjugation M ω = M −ω implies a change in the frequency variable.

B. Symmetries in systems with detailed balance
Above we saw that the stabilization of a T-symmetric effective Hamiltonian, Θ = Θ T from a dynamical process containing a Hamiltonian contribution, H, requires the Hamiltonian to transform as H → −H T . But how can this be reconciled with the familiar case of thermal equilibrium, where a T-symmetric effective Hamiltonian Θ = βH forms via thermalization (likewise an irreversible process) of a system with T-symmetric H? The resolution to this seeming paradox lies in a point mentioned in the introduction, namely that the T operation appropriate to the description of Markovian irreversible dynamics leaves the physical time parameter, t, untouched. By contrast, physical time reversal in unitarily evolving systems is described by an operation T 0 ≡ T • E 0 , where E 0 is the operation E 0 : t → −t in a 'theory space' containing t as an external parameter. Systems at thermal equilibrium are not time reversal symmetric in the strict sense of unitary evolution. However, they obey a principle of 'micro-reversibility', or detailed balance, which makes the reflection of time a meaningful operation.
To see how this comes about, consider the correlation function tr(ρÂ tBt ) with Heisenberg evolved oper-atorsÂ t = e iĤtÂ e −iĤt , and thermal, normalizedρ = Z −1 e −βĤ with partition sum Z. We work in second quantized representation, as indicated by the carets. As a straightforward consequence of the cyclic invariance of the trace, one obtains the Kubo-Martin-Schwinger (KMS) relation [23,24], tr(ρÂ tBt ) = tr(ρB t Â t+iβ ), i.e. an operator reordering relative toρ at the expense of a shift of the time parameter into the complex plane. For Heisenberg evolved operators, this equation is based on the specific form ofρ, and hence is a signature of the thermal state. Time reversal plays no role, up to now. Now use the relation to fuse the KMS relation with an anti-linear transformation in Fock space. (In a first quantized representation, the auxiliary relation follows from tr(Θ † T ) = tr(U TΘ † U −1 T ) = tr(Θ T ) = tr(Θ). For a verification in second quantized representation, see Appendix E1). A straightforward combination of Eqs. (42) and (43) then leads to where we assumed T-invarianceĤ T =Ĥ of the Hamiltonian.
This equation states the invariance of two-time correlation functions under a simultaneous application of the Fock space anti-linear transformation T and the transformation E β (Eq. (2)) acting on functions of time E β f (t) = f (−t−iβ). Provided thatÂ T =Â andB T =B are T-invariant, the combined operation T β = T • E β of Eq. (3) defines a symmetry of the correlation function. Notice that the operation T β is the product of two transformations T and E β acting in different spaces, namely Fock space and functions defined over Fock space, respectively. For β → 0 we obtain the standard operation of quantum mechanical time reversal, i.e. anti-linear T followed by an inversion of time. Here, the 'infinite temperature limit' simply means that in this case, the symmetry makes a statement for unrestricted traces of operators.
The above symmetry of correlation functions under thermal time reversal T β suggests that this operation might define a symmetry of the theory in general. To verify this expectation, we consider the Keldysh functional [87][88][89][90][91][92] and analyze how its action transforms under T β . We do not assume a thermal distribution just yet, but will obtain it as part of the criteria required for invariance. In Keldysh language, the inversion of time amounts to (a) an exchange of the forward (+) and the backward (−) time contour, (b) a contour dependent sign due to the reverse ordering of fermion/Grassmann fields along the contour, summarized jointly as ν ± → ±ν ∓ , where ν ± = (ν c ± ν q )/ √ 2 are the fields on the contours ±, and (c) an inversion of the time parameter ν(t) → ν(−t), or ν ω → ν −ω (which effectively changes the sign of all terms odd in frequency parameters such as (dω)η T −ω ων ω ). In the ν = (ν c , ν q ) T representation via 'classical' and 'quantum' fields, the combined effect of these operations assumes the form ξ ω → (iτ y )ξ −ω , ξ = ν, η. This needs to be supplemented by (d) a shift by iβ into the complex time domain. In the frequency representation, the shift operation on the individual contours assumes the form of a multiplicative factor exp(±ωβ/2). Turning to the (ν c , ν q ) form, we obtain the full representation of E β , as for both ξ = ν, η. Having established the representation of E β on Keldysh fields, we now need to figure out what effect T β = T • E β has on a Keldysh action with kernel M as in Eq. (39), with generally frequency dependent K ω = H ω − iD ω , and P ω . Referring to Appendix F for details, we find that the theory is invariant under thermal time reversal T β provided that (i) where the notation emphasizes, that the symmetry not only tolerates, but actually requires frequency dependence of the operators P and D. The latter is constrained by condition (ii) which is an implementation of the fluctuation-dissipation relation. This identity conditions fluctuations (P ) and dissipation (D) to each other, via a frequency dependent factor which one may consider the definition of the global equilibrium temperature. Notice how the conditions for T differ from those listed for the non-equilibrium case in table II, with sign changes owed to the active transformation of the time parameter. By contrast, the transformation C does not relate to the time parameter in either case and remains unchanged. However, the combined transformation S = T•C inherits the changes from T, The unmodified transformation laws under C, together with the modified ones in Eqs. (46,48), yield the equilibrium column of table II, generalizing the standard symmetry classification [5,22] to finite temperature settings.
Reflecting the fact that for Gaussian systems in equilibrium Θ = βH, these conditions are identical to those in the third column (Θ/Γ) of table II. While our analysis was performed for a Gaussian setting, it is natural to expect that the same symmetries characterize the self energies Σ(D, P ) forming in an interacting system relaxing into an equilibrium configuration at arbitrary temperature T . In this case, the symmetries are inherited from that of the microscopic parent theory under Eq. (3), where T acts on Fock space operators and β is set by the temperature of a background bath determining the system's temperature. (The role of this bath can be played by the system itself, in which case the temperature is determined by the energy of the initial state from which the thermalizing evolution departs.) We finally note that, as in the complementary out of equilibrium case, two gaps stabilize a topological phase: a spectral gap |eigenval(H)| > 0, and a purity gap in Γ eq = tanh βH/2, realized for temperatures smaller infinity, β > 0, as long as the spectral gap is open [72].

C. Scope of the equilibrium symmetry conditions
Our discussion above showed that their combined symmetry under Fock space operations, X, and generalized time reversal, E β , defines the symmetries of microreversible systems different from that of out of equilibrium systems. This makes one wonder just how general the scope of the micro-reversible framework is. Can it be extended beyond the category of equilibrium systems?
We first note that a symmetry class in the sense of our present discussion is defined by a set of operators (H, D, P ) sharing a certain set of conditions under application of X = T, C, S. (The realization of the symmetry in the dynamical evolution may or may not include an additional transformation of time under E β .) Individual deformations of these operators do not leave the symmetry class, provided they do not violate the symmetries of the constituent operators. This is an important disclaimer. For example, a lattice Hamiltonian invariant under the combined application of lattice inversion and time reversal does not define a class, because the addition of static disorder would violate the symmetry.
In this reading, the irreversible relaxation into an equilibrium configuration does not define a class, as it requires the 'fine tuning' D ∝ P , Eq. (47). Physically, a configuration D ∝ P builds up at long times when a system acts as its own thermalizing bath; before reaching that stage, traces of the initial state -possibly with different symmetry properties -may remain visible. In this sense, D ∝ P defines an attractive 'surface' in the 'space' H, D, P , provided the conditions for thermalization are met.
To make this point more concrete, we consider the long time limit of the covariance matrix, obtained from Eq. (40) as where in the second line we emphasize the dynamical interpretation of the covariance matrix by defining the retarded and advanced propagators For simplicity, we neglect optional frequency dependences in H and D, but in view of Eq. (47) not in P ω . On this basis, we ask under what circumstances does the covariance matrix have symmetry under T. Out of equilibrium, the answer is given by the first three columns of table II. To see this in explicit terms, compute Γ T = U T ΓU † T as Provided the symmetries hold as stated, a change of variables ω → −ω brings us back to the original expression, Γ T = Γ. We repeat that this symmetry requires oddity H T = −H of the system Hamilton operator.
In equilibrium, we have the 'fine tuning' Eq. (47) which implies a different option to establish T: In this case, the numerator may be written as −2iP ω = 2 tanh( βω , and the covariance matrix assumes the form i.e. an integral over the spectral function of the system (in general broadened by the coupling to the bath establishing equilibrium) over the Fermi-Dirac distribution function. In this case, Γ = Γ T holds if H = +H T , as required by the equilibrium column of table Eq. (II). However, notice that this symmetry crucially relies on the proportionality D ∝ P ω , Eq. (47). It gets broken by even mild departures away from equilibrium realized, e.g. by coupling to two baths kept at different temperatures, T 1 and T 2 . In this case, we have D = D 1 + D 2 which in general is no longer proportional to P 1 + P 2 = i tanh( β1ω 2 )D 1 + i tanh( β2ω 2 )D 2 . Of course, the out of equilibrium symmetry relation still holds, provided H, D, P satisfy the required criteria.
This also provides us with the opportunity to point out that both in and out of equilibrium, we never encounter the Hermitian adjoint, η, as a natural part of symmetry operations. For example, this operation would act on our Green functions as Considering the time representation of these propagators, Ref. [33] noted that (G + (t)) η = G − (−t) is a natural operation exchanging retarded and advanced propagators consistent with causality. Combined with T, this defines the operation T η ≡ T • η (cf. Eq. (6a) of Ref. [33]), there introduced as the "unique way to extend Hamiltonian symmetries to Lindbladian symmetries". Our conclusions are different. As evidenced by Eq. (49), the advanced and retarded propagators appear in Lindbladian evolution in combination. Application of the standard antilinear operation T (no η involved) to the state exchanges G + ↔ G − . Since these two operators appear in a paired configuration, cf. Eq. (49), the operation η has no place in the symmetry analysis of Lindbladian state evolution. However, it is of relevance in cases where the symmetries of the dissipation generators K = H − iD are considered in isolation. In this case, H = H T and D = D T implies a symmetry K Tη = H T −iD T = K, with physical consequences, e.g., in the statistical theory of the decay of resonances of open quantum systems (cf. Ref. [93]).
In view of the fact that a symmetry of K under T η may have observable consequences for open quantum systems, one may ask just how general such symmetries are. For example, it is natural to expect that if a system and its environment are time reversal invariant in the sense of unitary quantum mechanics, this symmetry is inherited by the dissipative operator K after the integration over environmental degrees of freedom. In Appendix G we show that this is not the case in general. Only in equilibrium does quantum mechanical time reversal invariance guarantee T η invariance of the reduced theory. This finding underpins our general statement that time inversion out of equilibrium is meaningless in general.
Finally, we note that the classifications of non-Hermitian matrices K = H − iD [17][18][19] all have in common that they preserve the transformation laws of the Hermitian contribution H under the anti-unitary symmetries familiar from the ground state classification. In light of the above discussion, such an extension to non-Hermitian matrices is only possible under conditions of global thermodynamic equilibrium.

VI. CONCLUSIONS
In this paper, we classified the symmetries governing the dynamics of open fermionic quantum matter. Symmetries were defined as linear or anti-linear transformations -represented in Fock space, or in the first quantized language of matrices for free systems -leaving the irreversible equations of motion invariant. While this rationale resembles the one applied in the identification of symmetries in unitary quantum time evolution, two principles make the out of equilibrium case different. First, unitary state evolution, or, somewhat more generally, the micro-reversible approach to a thermal equilibrium configuration, is governed by a single linear operator, the system Hamiltonian. By contrast, the out of equilibrium generators considered here comprise three linear operators, describing unitary evolution, dissipation, and fluctuations, respectively. All three must obey individually defined symmetry conditions for the full dynamics to be symmetric. The second difference concerns time itself. In unitary dynamics, the application of anti-linear symmetries (i → −i) is matched with an extraneously imposed inversion of time (t → −t) to leave the quantum time evolution operator, exp(−iĤt), invariant. Within the more general class of equilibrium processes, this operation is generalized to the shift inversion E β , Eq. (2), (E β f (t) = f (−t − iβ)), likewise designed to keep the dynamics invariant. Combined with the anti-unitary Fock space transformation it defines the thermal time reversal T β = T • E β , extending quantum mechanical time reversal to irreversible equilibrium dynamics. However, out of equilibrium time reversal becomes unphysical, meaning that anti-linear symmetries -featuring in six out of ten symmetry classes -have a fundamentally different representation.
In view of these differences, it is remarkable that the stationary states stabilized by equilibrium or nonequilibrium dynamics can be fully classified by identical symmetries. For example, when we say that a system is in symmetry class BDI, what we mean is that it obeys an anti-linear T-symmetry (squaring to unity) and a linear C-symmetry, likewise squaring to one. If the stationary state is Gaussian, its effective Hamiltonian can be represented as a real matrix possessing a block-off diagonal 'chiral' structure. The identical realization of symmetries relies on the stationarity of the asymptotic states -where the meaning of time is lost, or, more formally, time inversion E β acts as an identity and all symmetry operations reduce to their action in Fock space. It implies a strong principle of universality with obvious practical consequences. Notably, the information contained in the periodic table of topological insulators and superconductors universally applies to the classification of Gaussian stationary states both in and out of equilibrium.
However, the equivalent representation of symmetries in the stationary limit does not extend to the dynamical processes stabilizing them. Here, the presence or absence of E β is key and leads to the identification of twenty 'dynamical symmetry classes', distinct by the symmetry representations in either case. Ten of them describe the asymptotic approach towards a stationary equilibrium configuration, the other ten the approach to stationary non-equilibrium. Where the former assume the form of the ten well-known symmetry conditions for fermionic Hamiltonians, the latter require 30 = 3 × 10 conditions for the generators of unitary evolution, dissipation and fluctuations, respectively. The hallmark of the different symmetry representations in-and out of equilibrium is the sign difference in the T-transformation of the Hamiltonian contribution to the generator of dynamics. For example, the approach to a BDI symmetric non-equilibrium stationary state, admits a Hamiltonian contribution which, however, must be odd under the BDI symmetry rule for Hermitian matrix generators, while in equilibrium evenness is required. Referring back to the different treatment of time, the equilibrium symmetry classes reflect invariance of the dynamical approach under the joint application of Fock space symmetries and E β , while the non-equilibrium classes do not engage the latter.
In this paper, we illustrated the above symmetry principles, and the consequences for state topologies on the simple case of the BDI chain. However, it is relatively straightforward to construct other realizations by reverse engineering. Starting from a model effective Hamiltonian Θ of specified symmetry and topological ground state, one thus asks which Lindblad operatorsL α (see Eq. (10)) stabilize this state. By design, the dissipation and fluctuation generators defining these operators via Eqs. (11) and (13) are then conditioned via the symmetry relations of table II.
While the emphasis in this work has been on the bulk classification of phases, the next stage will be a more thorough exploration of the physics at the edge. Ideally, one would like to use the dissipatively stabilized topological edge as a resource for the storage and manipulation of quantum information. This requires the decoupling of the edge from the very physical mechanisms stabilizing it. For example, the edge space will be decoupled if it is realized as the dark space of an effective Lindblad equation. However, it will be interesting to explore if the isolation of the edge space can be effected in different ways, building on combined principles of symmetries and conservation laws as in Ref. [61].
Another direction of research concerns the dynamical processes leading to a stationary limit. In this paper, stationarity was attained in a competition of dissipative damping (K = H − iD), and fluctuations (P ). However, there are alternative ways to describe the quantum stochastic process driving the relaxation: Inspection of the Lindblad equation shows that it contains the non-Hermitian combinations H − P and H + P as operators acting to the left and right of the density matrix, while D − iP acts from both sides. In this decomposition, the first term describes the short time relaxation of quantum trajectories, interspersed by 'quantum jumps' described by the last. The competition between the two can be accessed in dynamically resolved ways by post-selection or measurement protocols [94]. It will be interesting to explore topological signatures in full counting statistics [95,96] in the above language. This may define topological structures of the symmetry constrained matrix operators H, D, P different from those considered in this paper.
Finally, one may ask how bosonic systems fit into the general framework. While the definition of fundamental symmetries extends to bosonic Fock spaces, the manifestations of these symmetries in concrete states are strikingly different: individual states can be multiply, or even macroscopically occupied, in which case the dynamics becomes (semi)-classical and quantum noise less of an issue. The action of state transformations confined by symmetries may be non-compact (think of the bosonic Bogoliubov transformation), and the stability of macroscopic stationary states becomes an issue. The latter compromises topology of Gaussian states of bosons [97] and makes the presence of interactions necessary [98]. At the same time, macroscopic bosonic quantum states likely define a more natural application field for the physics of non-Hermitian matrices than the strongly fluctuating fermion matter discussed here.
From a yet more general stance one may notice that, as with the physics of unitarily evolving quantum matter, the short range entangled symmetry protected fermionic phases considered here represent a relatively simple form of topological matter. With promising first steps taken in concerning the fate of fractional Quantum Hall states in open systems [99], the fascinating problem of extending the framework to dissipative variants of fractional or long range entangled matter is still out there and awaiting exploration.
research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.
We dedicate this work to the memory of our friend and colleague Federico Tonielli, whose brilliant promise in Theoretical Physics was cut short at age 27.

Appendix A: Symmetry classes
For the convenience of the reader, we here summarize the 10 symmetry classes, along with the dimensions where they admit topological ground states with Z or Z 2 classification. The labels Z, Z 2 , 0 in the table denote the possible topological invariants. For example, the system considered in section IV is defined in d = 1 and symmetry class (T, C, S) = (+1, +1, 1), or BDI.

Appendix B: Action of unitary symmetries
Here we elaborate on the condition, that anti-unitary symmetries T, C, S be realized within the irreducible representation spaces of the system's unitary symmetries U. This is best explained on the basis of examples. For instance, the plain operator exchange Ca i C −1 = a † i defines a symmetry structure of every free fermion operator, O, entirely on the basis of Fermi statistics: the operation acts on the Nambu operators as A i → A † i = σ x A i , and using the equivalence of these representations, where the trace comes from the Kronecker δ in {a i , a † j } = δ ij . Momentarily ignoring the trace, we have the relation O = −σ x O T σ x , which for Hermitian O defines an operator of class D. Now consider a situation with particle number conservation. In this case,Ô commutes with the unitary operator U = exp(iαN ), whereN = a † i a i and α ∈ u(1). However, C does not, CN C −1 = N −N . The operation C couples different sectors of conserved particle number such that the above principle is violated. To understand the consequences, note that number conservation implies block diagonality in Nambu space, O ≡ bdiag(o, −o T ). The operation C connects the two blocks, however, absent a physical coupling remains meaningless. The example illustrates how a second quantized operator of relatively higher symmetry (number conservation) can have a lesser symmetry on the matrix level (just Hermiticity, class A, rather than D) -in other words, symmetries or 'structures' on the first quantized level need not be rooted in actual symmetries of the many-body context. However, if particle number conservation is violated by, e.g., an 'order parameter' a † i ∆ ij a † j , the previously isolated sectors of definite number get combined to an enlarged representation space. C now acts within this space and does define a meaningful 'BCS' matrix structure Similarly, consider the example of spin rotation symmetry from Sec. II A, U s , where the plain C : a i → a † i does not commute and violates the unitarity principle, while C s : a σ → (σ y ) σσ a † σ does not. The operation C s thus defines a symmetry class (C) as O = −σ y O T σ y .

Appendix C: Symmetries in Gaussian state evolution
We here derive table II by subjecting the Lindblad equation (16) to the symmetry operations in their second quantized incarnation, Eq. (4). T-invariance -Application of T in second quantized incarnation leads to This equation becomes identical to the untransformed one, providedĤ T = −Ĥ,P T = −P , andD T =D, P T = −P . Comparison with table II shows that these conditions are met if the matrices H, P, D satisfy the conditions listed in the first row. C-invariance -In a similar manner, the application of C yields We change the representation of the jump term in the second line as 14,15) = A T U † C (D + iP )U CρC A † (9) = A T (−D C + iP C )ρ C A † , and substitution back into the equation yields ∂ tρC = − i(Ĥ C −P C )ρ C + iρ C (Ĥ C +P C ) By the same rationale as in the previous case, the invariance of the equation requires the matrix transformations listed in the second row of table II. S-invariance -Testing for S is not necessary, as it is a consequence of the combined presence of C and T. However, it is instructive to see how this symmetry manifests itself in the Lindblad equation without reference to the composition. Application of S leads to Once again, the term in the second line requires special attention:

Appendix D: Symmetries in Keldysh representation
We here discuss how the symmetry of Gaussian states as expressed by the last column of table II leads to Eqs. (41) for the generally non-Markovian matrix kernels M generating the dynamics. T-invariance -According to table II, T-invariance means the existence of a unitary matrix U T such that Γ = Γ T = U TΓ U † T . Starting from the representation (38) this becomes In the third line, we transformed variables U T ν → ν and η T U † T → η T (note that U T ωτ 1 U † T = ωτ 1 ), and changed ω → −ω in the frequency integral. A sufficient condition for the invariance of the Γ matrix then is (M T ) ω = U T M −ω U † T = −M ω , which is the first line of Eq. (41). C-invariance -In the path integral formalism, Cinvariance, Γ = Γ C = −U T C Γ TŪ C , is probed as where in the crucial fourth line we swapped the order of variables both in the action and the pre-exponential variables (picking up a sign in the process). Except for a different naming of the dummy variables, η ↔ ν, the final expression equals the original one, provided M satisfies the C entry in Eq. (41). S-invariance -Finally, S-invariance, Γ = Γ S = −U T S Γ †Ū S is established as which leads to the final entry in Eq. (41).
Microscopic model -We consider a model with two fermionic modes a = (a 1 , a 2 ), each coupled to a bath with modes b µ = (b µ,1 , b µ,2 ). We choose bath temperatures T 1,2 , which for T 1 = T 2 breaks equilibrium. The system-bath Hamiltonian readŝ Here, the 1, Σ x,y,z act in band space andĤ s ≡Ĥ +Ĥ int describes the system Hamiltonian, featuring a free part and a density-density self-interaction, which we will see plays an important role. The termsĤ c andĤ b model the coupling to a harmonic bath. With the real couplings, g µ , the above system-bath Hamiltonian is time reversal invariant, with transformations in band space as The microscopic Keldysh action for the non-interacting contribution to this setting reads, in the frequency domain (from now on, a = (a c,1 , a c,2 , a q,1 , a q,2 ) T , b µ = (b c,µ1 , b c,µ2 , b q,µ1 , b q,µ2 ) T ) where H is the above system Hamiltonian matrix, the Pauli matrix τ x acts in Keldysh space, and t ± ω = 1 2 (tanh β1ω 2 ± tanh β2ω 2 ). The infinitesimal parameter κ > 0 defines the causality of the Green functions, and it fixes the state of the bath: inverting the above matrix we find G cc,µ = 2πiδ(ω − µ )(t + ω 1 + t − ω Σ z ), where the departure from equilibrium is measured by a non-zero coefficient t − ω , and the appearance of a matrix Σ z in band space. Notice that the distribution mismatch breaks the T-invariance of the action, Σ z,T = −Σ z . Although the distribution functions couple to the action only infinitesimally via κ, they do feed back into the system dynamics on an O(1) level, as the following discussion shows.
Effective Lindblad model -Integrating out the bath, we obtain the effective system action S eff = (dω)ā T G −1 eff,ω a − λ dt(n 1c n 2q + n 1c n 2q ), with , where d ω = π µ g 2 µ δ( µ − ω) and n ic =ā ic a ic + a iq a iq , n iq =ā ic a iq +ā iq a ic , i = 1, 2, and we neglected a Lamb shift renormalizing the system Hamiltonian, which is unimportant for the present discussion.
At this level, it looks like the induced matrix generator K ≡ H − id ω 1 is T-symmetric. However, this changes once the interaction Hamiltonian is taken into account. To first order in perturbation theory, this generates a selfenergy correction ∼ λ (dω)G + eff,ω P eff,ω G − eff,ω , where the Σ z matrix contained in P eff,ω reflects the absence of equilibrium. Substituting this expression into Eq. (G4), we induce a term ∝ Σ z in K. This term is a consequence of the sensitivity of the self energy to the bath distribution functions, which for T 1 = T 2 break time reversal. In this case, K T = K and the symmetry under non-Hermitian time reversal T η [33] is lost.
The upshot of the above discussion is that out of equilibrium, the microscopic T-symmetry of a sys-tem+environment Hamiltonian does not stabilize an induced T η symmetry of the matrix generator K. Only in equilibrium, the full symmetry of the theory (including the bath distribution functions) under T β descends to a T η symmetry of this operator, as outlined in section V. Out of equilibrium, the only way to stabilize an antilinear T-symmetry of a dynamical evolution and its stationary phases (then without reference to time reversal) is via