Dissipation-induced topological insulators: A no-go theorem and a recipe

Nonequilibrium conditions are traditionally seen as detrimental to the appearance of quantumcoherent many-body phenomena, and much effort is often devoted to their elimination. Recently this approach has changed: It has been realized that driven-dissipative could be used as a resource. By proper engineering of the reservoirs and their couplings to a system, one may drive the system towards desired quantum-correlated steady states, even in the absence of internal Hamiltonian dynamics. An intriguing category of equilibrium many-particle phases are those which are distinguished by topology rather than by symmetry. A natural question thus arises: which of these topological states can be achieved as the result of purely dissipative Lindblad-type (Markovian) evolution? Beside its fundamental importance, it may offer novel routes to the realization of topologicallynontrivial states in quantum simulators, especially ultracold atomic gases. In this work I go beyond previous studies in giving a general answer to this question in the context of Gaussian (“integer”) topological states, concentrating on 2D Chern insulators as the main example. On the one hand I prove a no-go theorem which shows that such dynamics cannot lead to a topological robust (gapped Liouvillian) unique pure steady state as the result of finite-range Liouvillian. On the other hand, I construct a recipe showing that, under the above conditions, a pure topological steady state may result if exponetially-local Liouvillian is allowed. If strictly local evolution is insisted upon, a mixed steady state arbitrarily close to the desired pure state may be obtained. I will also show how such dynamics could be realized with ultracold atoms and similar systems, and how the resulting states may be detected and topologically-classified. Extension to other types of topological insulators and superconductors is also discussed.

the Lindblad type [2], which are driven towards quantum correlated steady states [3][4][5]11]. Thus, instead of having to tweak the Hamiltonian of the system as well as to couple it to a bath to drive it towards the ground state, one could arrive at the same state by the coupling to the reservoir only, without any internal Hamiltonian dynamics 1 .
Can the states achieved in this way be topologically nontrivial? Previous work concentrated on chiral p-wave superconductors in 1D and 2D, and provided intriguing results: In the former case it was found that a pure topologically-nontrivial dark state can be realized in a concrete model of ultracold fermions coupled to a bosonic environment, and that this state is accompanied by the appearance of decoupled Majorana modes at the ends of the system, as in equilibrium [7,15]. In the latter case it was shown that the pure dark state so generated is actually topologically trivial (though decoupled Majorana modes still exist at the centers of vortices) [8,9]. Modification of the procedure allowed to arrive at topologically-nontrivial state, which is however of low-purity and thus far from the corresponding equilibrium state [13,16]. The limitations on attaining topological order dissipatively were also studied [6,10,14,17]. In the context of superconducting nanostructures it has been shown how augmenting Hamiltonian dynamics by bath engineering may aid in stabilizing conventional (non-dissipative) fractional quantum Hall states [12]. The scope and generality of these results, as well as the situation for other types of systems, remains unclear. The issue of topological classification of mixed nonequilibrium steady states is also far from resolved [48][49][50][51][52][53][54][55][56].
In this work I give a general answer to these questions in the context of fermions with quadratic driven-dissipative Lindblad dynamics, concentrating on Chern-insulator/lattice integer quantum Hall states, that is, 2D systems in symmetry class A [25][26][27]. I formulate and prove a no-go theorem, which shows that if the Liouvillian superoperator is composed of terms with finite range in real space, exponential decay (at a finite rate in the thermodynamic limit) in time towards a pure steady state with nonzero Chern number (or any other topological state above 1D) cannot be attained. However, a mixed steady state as close as desired to such a pure state is possible. Alternatively, a pure steady state could be reached if the Liouvillian is exponentially (or even faster)-local in real space. I present a concrete recipe for the implementation of the scenarios which are allowed by this theorem. A possible physical system for the implementation of the recipe is schematically depicted in Fig. 1: Ultracold atoms occupy an optical lattice, which may be modified so as to completely suppress Hamiltonian tunneling dynamics, in manners to be specified below. Rather, the lattice is coupled in a local fashion (e.g., optically) to reservoirs which can supply or take atoms out of the system 2 . I show below how to properly tune these couplings according to the mentioned recipe, so that the induced dissipative dynamics may drive the system to a steady state with the special properties specified above. I also explain the topological classification of the resulting (possibly mixed) state vis-à-vis the previous discussion in the literature, and go beyond it by proposing an observable allowing to detect the quantized topological index in a cold-atom experiment.
The paper is organized as follows: After setting the notation and posing the main results as a theorem in Sec. 2, I present the general recipe in Sec. 3, and then show in Sec. 4 that the recipe cannot be improved via the no-go part of the theorem. I then present an example and a concrete cold-atom setup for the realization of the recipe in Sec. 5, and explain how a topological index could be defined and measured in Sec. 6. Finally I conclude Figure 1: Schematic depiction of a possible realization of the system considered in this work: Fermionic atoms (red spheres) occupy an optical lattice. The lattice is engineered to suppress ordinary tunneling. Instead, the system is connected (e.g., optically) to reservoirs (elliptic blobs) by local couplings (black arrows). The reservoirs exchange atoms with the system and drive it towards a steady state with specific properties, such as being topologically nontrivial.
in Sec. 7 with a discussion of the extension of the results to other topological classes, as well as of future research directions.

Preliminaries and statement of the main results
Let me start with some preliminaries. The state of a system coupled to an environment is described by a density matrix, ρ. When the coupling to the reservoir is weak, the reservoir is much larger than the system (and therefore is negligibly-affected by the system), and has much faster dynamics than the system (conditions which are widely applicable in practice in optically-driven systems, especially ultracold gases), one can integrate out the environment and arrive at a Markovian (memoryless) Lindblad master equation for the density matrix of the system [2]: where the Liouville superoperatorL contains two contributions. The first one describes the Hamiltonian part of the evolution, which may include "Lamb shift" corrections due to the coupling to the bath. The second one is the dissipative part. L i are known as the "Lindblad operators", with a corresponding nonnegative Hermitian rate matrix γ ij , whose eigenvalues are the various decay and decoherence rates. This is the most general Markovian equation that ρ could obey, which would preserve its Hermiticity, positivity, and overall normalization, Tr(ρ) = 1.
In this work we are going to concentrate on quadratic Lindblad dynamics. For an arbitrary number of baths, if the total number of fermions in the system and the baths is conserved (but fermions may be exchanged between the system to the baths), the most general Lindblad equation assumes the form where a m is an annihilation operator of a fermion in a single particle state m,h S mm the Hermitian single-particle system Hamiltonian (which may include corrections due to the coupling to the baths), and γ out mm and γ in mm are nonnegative Hermitian rates of particles leaving and entering the system, respectively. For completeness the derivation is presented in Appendix A. Although I will not dwell much on this case, let me mention that when pairing is present, one should add to the Hamiltonian terms of the form ∆ mm a m a m +h.c., and similarly add Lindblad terms of the form γ pair mm [2a m ρ S a m − a m a m ρ S − ρ S a m a m ]/2 + h.c., with antisymmetric matrices ∆ mm and γ pair mm . Let me now pose the main question of this study: Can one engineer quadratic Lindblad dynamics for a many-body fermionic system with the above general form which obeys the following requirements: • Unique steady state ρ ss ; • Finite gap in the spectrum of the Liouville superoperatorL, ensuring that the steady state is approached (in terms of few-particle observables) exponentially fast in time, at a finite rate which independent of the system size 3 . The finite gap also serves to suppress the effects of small perturbations, as in equilibrium; • Finite range, namely, in a real space representation the matricesh, γ out γ in , as well as ∆ and γ pair , have elements that only connect real space orbitals whose distance is below some upper bound. This is to be contrasted with "local", where elements connecting orbitals at any distance are allowed, provided they decay at least exponentially fast with the distance; where the resulting steady state ρ ss is: • Pure, that is, can be written as ρ ss = |Ψ d Ψ d | in terms of a pure dark state |Ψ d ; • Topologically nontrivial, that is, the dark state |Ψ d is characterized by a nonzero topological index, e.g., the Chern number in 2D.
As I will show, one has to give up the strict realization of at least one of these desiderata. This can be stated more formally as follows: Theorem. A finite range quadratic 2D fermionic Lindbladian cannot have a pure unique steady state which is approached exponentially fast in time at a finite rate independent of the system size, if the pure steady state has nonzero Chern number. However, such a steady state is possible for a local Lindbladian. Alternatively, with a finite range Lindbladian one may obtain a mixed steady state which is approached exponentially fast in time at a finite rate, and which is arbitrarily close (exponentially fast in the range) to a pure steady state that has nonzero Chern number.
I will start by giving in Sec. 3 a recipe for the realization of either of the allowed cases, which is physically motivated by the possibilities offered in cold-atom systems, of the type depicted schematically in Fig. 1. Then I prove in Sec. 4 the impossibility ("no-go") of the excluded case, showing that the recipe is the best one can indeed achieve. For concreteness, I will concentrate on 2D systems in symmetry class A [25][26][27], that is, lattice quantum Hall/Chern insulators, but will comment on other symmetry classes along the way and in the Conclusions, Sec. 7.

Recipe
I will now describe a recipe for realizing a topologically-nontrivial state of noninteracting fermions as a unique steady state of local dynamics obeying the constraints allowed by the Theorem just stated. As we will see, purely-dissipative dynamics [no Hamiltonian part, H = 0 in Eq. (1)], will be sufficient. I will assume that the topologically-nontrivial state in question may arise as the ground state (with the lowest band filled) of some finite-range (to be relaxed below to just local) tight-binding reference Hamiltonian of fermions on a clean d-dimensional lattice (concrete examples will be given below) 4 : where i, i label lattice unit cells, and µ, µ = 1, . . . , n s is a basis of states within a unit cell (possibly including spin), and where a iµ obey standard fermionic anticommutation relations 5 . One may then Fourier transform to operators a kµ with a given crystal momentum k, and diagonlize the resulting k-dependent n s × n s Hamiltonian matrix h ref µµ (k) to find the eigenstates and eigenenergies (λ = 1, . . . , n s is the band index), Let us first assume that the lowest band of the reference Hamiltonian is flat [dispersionless, cf. Fig. 2(a)] in which case we can set ε 1 (k) = 0, without loss of generality; later on we will relax this assumption. Suppose now that the lattice in question has been constructed, e.g., optically in an ultracold atomic system, but that the Hamiltonian tunneling dynamics has been suppressed (by means to be described in Sec. 5 below), cf. Fig. 1. A many-body state describing fermions filling up the lowest band (λ = 1) completely and leaving the other bands (λ ≥ 2) empty can then be reached from a trivial initial state via finite-range purely-dissipative dynamics. The key idea is to implement not the reference Hamiltonian itself, but rather purely-dissipative dynamics induced by modification of reference Hamiltonian. This modification involves introducing two internal states of the fermions (e.g., hyperfine states in cold atoms), one which is fully trapped by a confining potential and one which is not trapped in at least one direction, with annihilation operators a iµ and b iµ , respectively. The modified Hamiltonian is then similar to the reference Hamiltonian, but with each tunneling term being replaced by an optically-assisted hopping term which also modifies the internal state of the atom from a trapped to untrapped or vice versa. In the rotating frame this amounts to: This Hamiltonian must be supplemented by another term, which includes the confinement of the a fermions, e.g., in the perpendicular direction for a 2D optical lattice, and the possibility of motion of the b atoms along that direction. The latter can be modeled, for example, by assuming that the 2D lattice just introduced is one layer in a 3D lattice. Denoting the crystal momentum in that direction by q z (k is still the in-plane momentum) we define the annihilation operators b iµqz , such that b iµ = qz b iµqz / √ N z , N z being the size of the system (number of sites) in the z direction. Correspondingly we add to the Hamiltonian the term where ε z (q z ) is the dispersion due to the motion in the z direction, and ε 0 > 0 is the detuning between the atomic a-b transition frequency (including zero point motion energies) and the frequency of the light inducing the optically-assisted hopping, that is, the residual 1D kinetic energy of the atom after an a → b transition. I will assume that b atoms that escape the system run away to infinity without ever returning, i.e., that the b fermions effectively constitute a bath with a negative infinite chemical potential. Hence, the bath density matrix ρ E obeys Tr(ρ E b † iµqz b i µ q z ) = 0 and Tr(ρ E b iµqz b † i µ q z ) = δ ii δ µµ δq z q z . If the escape rate of the b fermions is much faster than the optical transition rates, every optically-assisted hopping event leads to a loss of an atom. The conditions mentioned above for the validity of the Lindblad Markov equation are thus fulfilled. Then, we may integrating out the b fermions [using the second form of Eq. (5)] in the manner described in Appendix A and derive the following Lindbladian for the trapped a atoms 6 : where, by the Fermi golden rule, γ out Here ν 0 is the 1D (z direction) local density of states per unit cell of the untrapped atomic internal state b emitted with kinetic energy ε 0 ; a concrete expression will be given in Sec. 5 below.
By our previous assumptions [ε 1 (k) = 0 and separated by a finite gap from the other ε λ≥2 (k) = 0], the escape rate is finite for all the higher bands (λ ≥ 2) but is exactly zero for the lowest band (λ = 1). Thus, if one starts at the trivial, completely filled state for the a atoms, it will exponentially decay towards the desired pure many body dark state |Ψ d , in which the lowest band is completely filled and all the others are empty. As a matter of fact, the nature of the final state does not depend on the Markovian approximation used to derive Eq. (7), and stems already from the Hamiltonian (5). This Hamiltonian coupled all the a atoms to the reservoir (and thus dooms them to leave the system sooner or later), except those that occupy the states of the lowest band of the reference Hamiltonian 7 . Let me note that, for the purely evaporative dynamics presented so far, the final state is not unique and depends on the initial conditions, since if one of the states in the λ = 1 band is initially empty, it will not be refilled. To remedy this I would later on introduce a "refilling" reservoir, which could similarly be engineered to fill in only states in the λ = 1 band. But before going into that, let us first continue the discussion of the evaporative dynamics, Eq. (7).
The success of the above evaporative procedure crucially depends on the flatness of the desired λ = 1 band of the reference Hamiltonian. However, as has been appreciated for some time [58][59][60][61][62][63] but proven only quite recently [64][65][66], a topologically-nontrivial (having nonzero Chern number in the absence of symmetries, or arising from a reference Hamiltonian with some symmetry and having a nontrivial topological index appropriate for that symmetry class) flat band cannot arise in 2D and above from a finite-range gapped reference Hamiltonian 8 , that is, with matrix elements h ref µµ (i − i ) which vanish for i − i larger than some fixed value. This implies similar range restriction for the Lindbladian (7) built out of it. In the next subsection I will show that this is not a peculiarity of the current construction, but rather a result of a general no-go theorem. Here I will instead discuss what one may still achieve and how.
One possibility is to forgo finite range and content ourselves with locality, that is, , and hence γ out µµ (i − i ), which decay rapidly with i − i . An exponential decay is easy seen to be sufficient: One may replace the Hamiltonian by a projector onto the bands λ ≥ 2, which, due to the gap in the original Hamiltonian, will have exponentiallydecaying matrix elements in real space [67][68][69]. In fact, even faster decay is possible, as is exemplified by the Kapit-Mueller Hamiltonian [70], which is characterized by matrix elements decaying exponentially with the square of the distance, and leads to a lowest flat band with nonzero Chern number. In fact, I am not aware on any bound on the rate of decay, except that h ref µµ (i − i ) cannot have finite support (finite range). This is to be contrasted with the corresponding Wannier functions, which can only decay algebraically with distance for a band with nonzero Chern number [71][72][73].
Another possibility is to keep the finite range of the terms in the Lindbladian, and instead forgo the requirement of purity of the steady state. A finite-range reference Hamiltonian may have a very narrow topological lowest band, in the sense that its width is much smaller than the gap separating it from the next band, λ = 2 [58][59][60][61][62][63]. For convenience we set the zero of energy at the minimum of the lowest band [see Fig. 2 I now add a competing term, which tries to populate all the lattice states at a constant rate γ in . This could be realized by using yet another atomic species/internal state c in a setup similar to b except that the c bath is filled with a large number of atoms (effectively making it an infinite chemical potential reservoir), and that the matrix elements coupling it to a in the analogue of Eq. (5) are proportional to the unit matrix, δ ii δ µµ . The combined LindbladianL =L out +L in is quadratic, and thus has a Gaussian steady state (a unique one, sinceL in has no zero modes [57]), which is fully characterized by the positive and Hermitian single-particle reduced density matrix, G ss µµ (i − i ) ≡ Tr(ρ ss a † iµ a i µ ). The latter is diagonal in the basis of the eigenmodes of the reference Hamiltonian, with the steady-state mode occupancies n ss Choosing γ in in the range γ out 1,max γ in γ out 2,min , single-particle states belonging to the lowest band will be almost filled, while the other bands would remain almost empty, as desired. The optimal value of γ in will be of order (γ out 1,max γ out 2,min ) 1/2 , for which all the deviations of the occupancies of the mixed steady state from those of the pure ground state of the reference Hamiltonian (with chemical potential in the gap between the bands λ = 1 and λ = 2) are of the order of the flatness ratio ε ref 1. The latter decreases exponentially with the range of the terms in the reference Hamiltonian [63,[67][68][69]. This immediately translates into an exponential decrease (with the range of the Lindbladian) of the deviation between few particle observables (observables depending on modes whose number does not scale with the system size) of the mixed steady state and the corresponding pure state values. Therefore, the required range to obtain a fixed deviation scales logarithmically with the desired deviation size, independently of the system size. As for more general observables, they follow the 1 distance of the mixed steady state from the pure desired one, which also decays exponentially with the range, but with a prefactor proportional to the system size. Hence, the required range to achieve a fixed deviation of these quantities varies logarithmically with both the desired deviation and the system size.
Furthermore, any deviation of few-particle observables from their steady state value will decay exponentially in time at a finite rate which is independent of the system size, γ out λ (k) + γ in ≥ γ in . For example, the mode occupancies would decay as which is a particular case of Eq. (12) below. Hence, the time required to achieve a certain fixed deviation varies logarithmically with the desired deviation size, independently of the system size. Again, general observables will follow the 1 distance of the density matrix from the steady state ρ ss , which will feature a prefactor proportional to the system size. As a result, the time needed for the 1 distance to become smaller than some fixed deviation will scale logarithmically with both the desired deviation and the system size. Let me also mention that the finite rate ≥ γ in amounts to a finite gap in the gap of the Liouvillian, which would act to protect the system from perturbations (in analogy to the excitation gap in the equilibrium case [25][26][27]). It should also be noted that one may create a similar state without invoking the incoming Lindbladian (8). Instead one may start from a completely filled initial state and apply the outgoing Lindbladian (7) for a finite time ∼ (γ out 1,max γ out 2,min ) −1/2 . The resulting simplification of the setup should be weighted, however, against the non-uniqueness of the steady state, and thus the lack of stabilization of the state after its creation.
Finally, let me come back to the implementation a local (not short-ranged) Lindbladian with a pure unique steady state. For this purpose one would need to improve the construction of the refilling Lindbladian, Eq. (7), and choose a space-dependent γ in µµ (i − i ) built out of a flat-band reference Hamiltonian which is a projector into the λ = 1 band [similar but opposite to the construction of γ out µµ (i − i )], which will refill only states belonging to that band. With this one would have γ out λ=1 (k) = γ in λ =1 (k) = 0, and hence n ss λ (k) = δ λ,1 , corresponding to the desired pure steady state.

No-go theorem
I will now proceed to prove the negative side of the theorem stated in Sec. 2, namely that one cannot do better than the above recipe. This implies that one cannot achieve exponential decay in time at a finite (system size independent) rate towards a pure dark state with nonzero Chern number in 2D if the Lindbladian has a finite range. I will do so by showing that this would amount to having a finite-range Hermitian operator in 2D featuring eigenbands with nonzero Chern number separated by gaps from other bands, which is known to be impossible [64][65][66]. Similar results hold for other topological classes [66].
To set the stage for the proof, let us examine general quadratic Lindblad dynamics, described by Eq. (2). This gives rise to the following evolution equation for the singleparticle density matrix G (defined as in the previous Section but now not necessarily in the steady state): The steady state single-particle density matrix is the solution of ∂ t G ss = 0. The decay towards it is given by The positivity of γ out and γ in implies that the real parts of the eigenvalues of (γ out + γ in )/2 ± ih S are nonnegative. To ensure a finite decay rate towards the steady state in the thermodynamic limit, the lower bound of these real parts must be positive and finite as the system size goes to infinity, i.e., there should be a finite gap between the real part of the spectrum of (γ out + γ in )/2 ± ih S and zero. With this we can return to the no-go theorem. Its proof consists of two steps: 1. The main one is to show that a pure steady state implies that the positive Hermitian rate matrices γ out and γ in have flat bands with eigenvalue zero, corresponding to the filled and empty single-particle states, respectively. Since these bands should be topologically nontrivial, if γ out and γ in have finite range, they are gapless.
2. The second step is to demonstrate that the gaplessness of γ out and γ in implies that the real part of the spectrum of (γ out + γ in )/2 ± ih S is not separated by a finite gap from zero, hence the decay rate towards the steady state vanishes in the thermodynamic limit, in contradiction to our assumptions.
The first step results from considering the steady state. Let us work in the singleparticle basis where G ss is diagonal, which is a basis of Bloch states if the system is translationally invariant. Since ρ ss is Gaussian, if it is also pure, the eigenvalues of G ss , namely the eigen-occupancies n ss λ (k), are all either 0 or 1, corresponding to empty and filled states, respectively. We will write all the matrices appearing in Eq. (11) in terms of blocks in the empty-filled basis, e.g., . In this basis G ss = 0 0 0 I , with I being a unit matrix. Substituting into Eq. (11) we find that γ out 11 = 0 and γ in 00 = 0, which, due to the positivity of the rate matrices, also implies that γ out 01 = γ in 01 = 0 and γ out 10 = γ in 10 = 0. In addition we geth S 01 = 0 andh S 10 = 0, which will become important in the next step. Thush . This is a natural result: Completely filled (empty) modes are modes which do not couple at all to the outgoing (incoming) bath, and thus form a flat band with eigenvalue zero of the Hermitian matrix γ out (γ in ), which we assume to have a finite range. However, since the respective modes have a nonzero Chern number [64][65][66], this is impossible unless γ out (γ in ) is not gapped, i.e., the positive Hermitian block γ out 00 (γ in 11 ) has eigenvalues which are either zero or vanish in the thermodynamic limit.
I now turn to the second step. By Eq. (12), the last conclusion (gaplessness of γ out and γ in with respect to zero) implies for purely dissipative dynamics (h S = 0), or wheñ h S = 0 but there is only one empty or one filled band, the absence of a finite (systemsize-independent) exponential decay rate towards the steady state. The general case of h S = 0 and arbitrary number of bands requires a slightly more involved argument. We have observed above that the filled (empty) bands are flat bands with zero eigenvalue of the Hermitian matrix γ out (γ in ), whose elements are short-range in real space, and thus are trigonometric polynomials in the k basis. The filled (empty) bands wavefunctions thus obey linear equations with coefficients which are trigonometric polynomials in k, i.e., these bands form a polynomial bundle in the terminology of Refs. [65,66]. On the other hand, as mentioned above, a finite decay rate of G towards its steady state value G ss in the thermodynamic limit implies a finite gap between zero and the real part of the spectrum of (γ out + γ in )/2 ± ih S = γ out 00 /2 ± ih S 00 0 0 γ in 11 /2 ± ih S

11
, which holds separately for each of its two nonzero blocks. There will then be a finite gap between the spectra of the two blocks of, e.g., γ out /2 + ih S = γ out 00 /2 + ih S 00 0 0 ih S

11
, since the lower block has purely imaginary eigenvalues. By the arguments of Refs. [65,66] this would imply that the filled and empty bands form analytic bundles 10 . But since these bundles were argued to be polynomial, they cannot be topologically nontrivial [65,66], in contradiction to our assumptions. This completes the proof. Similar results would hold for other topological classes above 1D, if one enforces that γ out or γ in obey the respective symmetry [66].
Let us note that one can further extend the model to include cases of "superconducting" reservoirs, where in Eq. (1) the Hamiltonian may contain pairing terms with either two creation operators or two annihilation operators, and in the Lindbladian there could be "pairing" terms where L i and L † j are both creation or both annihilation operators. In that case it would be useful to define Majorana fermions, . One may then write down an analogous equation to Eq. (11) and use similar arguments to reach the same conclusion, namely that a topologically-nontrivial pure state with superconducting correlations, such as a chiral p-wave (class D [25][26][27]) state in 2D, cannot be obtained at a finite rate in time as the unique steady state of a gapped finite-range quadratic Lindbladian, but could be created if only locality is imposed. This in retrospect explains the difference found in previous studies between 1D and 2D dissipative pairing systems: In 1D the Kitaev chain [74] has a well-known flat-band limit, allowing a finite-range Lindbladian with a pure topological steady state, as found in Refs. [7][8][9]. In 2D this is not possible [64][65][66], so only a mixed steady state can be reached, as surmised in Ref. [13]. The results presented here bear some relation to the tradeoff between purity and short range in tensor network representations of topological states [69].

Example: Dissipative Hofstadter model with cold atoms
As an example for the usage of the above technique, I will apply it to realize a lattice integer quantum Hall state. Here the reference Hamiltonian is the Hofstadter model [75], a tight binding model with nearest-neighbor hopping on a square lattice pierced by a magnetic flux, nx,ny a nx,ny + Je i2παny a † nx+1,ny a nx,ny + Ja † nx,ny+1 a nx,ny + h.c., where α, the magnetic flux per plaquette in units of the flux quantum, is taken as rational, α = p/q (p, q are integers with no common factor). Fig. 3(a) shows its spectrum on a cylinder with periodic boundary conditions in the x direction for α = 1/5. It features n s = q = 5 bands as well as chiral midgap modes localized at the two edges of the system. The flatness ratio for the lowest band is below 4.2 · 10 −2 1. For the recipe one needs to implement the two-flavor analogue of this Hamiltonian, cf. Eq. (5). There are numerous proposals for engineering the necessary gauge fields artificially in cold-atom systems [35][36][37][38][39][40][41][42][43][44], many of which could be adapted to create the purely dissipative dynamics described here. As an example I will concentrate on the  (13) with periodic boundary conditions in the x direction (that is, on a cylinder whose axis is in the y direction) for α = 1/5 and J 0 = 2.93J (set to bring the minimal energy to zero). In between the q = 5 bands there are chiral midgap edge modes localized at the two ends of the system. (b) The spectrum of the steady state single-particle reduced density matrix, Eq. (9) [1 − n ss λ is plotted so that the order of states is as in (a)] resulting from the Lindblad dynamics, Eq. (7)-(8), with γ in = 0.2πν 0 J 2 . The four upper bands appear together at the top, but are actually separated by small gaps. (c) The same with the reference Hamiltonian (13) containing the next-nearest-neighbor hopping terms of the Kapit-Mueller Hamiltonian [70].
construction suggested in [76,77]. In the current context it would involve fermionic twolevel atoms on a 3D tetragonal optical lattice, with lattice constants d and d z in the xy plane and z direction, respectively. The b atoms are free to propagate in 3D, but the a atoms are confined to a single 2D plane with a steep enough confinement potential that the confinement level spacing is much larger than all the other energy scales introduced below (except ω 0 ). The lattice is modulated in the x and y directions in the way depicted in Fig. 4. Employing a tight-binding single-band description, this corresponds to a Hamiltonian H V = nx,ny V nx,ny a † nx,ny a nx,ny + nx,ny,qz ε z (q z ) + V nx,ny + ω 0 b † nx,ny,qz b nx,ny,qz , where a nx,ny annihilates an a atom in the Wannier orbital w a (x − n x d, y − n y d, z), including the confinement in the z direction, while b nx,ny,qz annihilates a b atom in the state nz w b (x − n x d, y − n y d, z − n z d z )e inzqzdz / √ N z , localized in the 2D plane but propagating with crystal momentum q along the z direction, which is N z -layers long. ω 0 is the a → b transition frequency (including differences in zero point motion energies), ε z (q) = 4J z sin 2 (qd/2), and the modulated potential is V nx,ny = [δ nx≡1 (mod 2) + δ ny≡1 (mod 2) ]ε y + [δ nx≡2 (mod 4) − δ nx≡1 (mod 4) ]ε x . The resulting mismatch in energies between neighboring sites prevents usual hopping for both the a and b states, provided the bare hoppings in the 2D plane as well as J z are all much smaller than ε x and ε y .
One may now turn on lasers with frequencies chosen so as to allow photon-assisted nearest-neighbor tunneling events which change the internal states of the atoms (a and b, or trapped and un-trapped, depicted as blue and orange, respectively, in Fig. 4). A laser beam with 3D wavevector k L = (k Lx , k Ly , k Lz ), frequency ω L , and Rabi frequency (amplitude) Ω L allows for transitions between a nxny and b n x n y ≡ qz b n x n y qz / √ N z provided that the frequency mismatch obeys 0 < ω L − (ω 0 + V n x n y − V nxny ) < 4J z . One may pick the following frequencies, as indicated in Fig. 4: • ε 0 + ω 0 ± ε y (dotted arrows), allowing nearest-neighbor state flip hopping along y, as well as along x between n x ≡ 3 (mod 4) and n x + 1; Figure 4: Lattice modulation scheme along the y (left) and x (right) directions corresponding to Eq. (14). Each lattice site may contain fermionic atoms with two internal states [blue and red, the latter not being trapped, corresponding to a and b in Eq. (16)] whose energies differ by ω 0 . A period-2 modulation with amplitude ε y is applied in both directions. Further period-4 modulation with amplitude ε x is applied in the x direction only. Optically assisted hopping is created by 6 laser beams with appropriate frequencies (green/pink solid/dashed/dotted arrows), whose directions are chosen so as to reproduce the phase factors in the Hofstadter Hamiltonian, Eq. (13) [76,77]; an additional beam with frequency ω 0 (not shown) allows for an on-site term. This gives rise to the Hofstadter system-bath coupling Hamiltonian (16), and thus to the evaporative Lindbladian (7).
Here 0 < ε 0 ε x , ε y is a small detuning. For each laser, the resulting optically-induced transition amplitude is J L = Ω L 2 e ik Lx nxd+ik Ly nyd dx dy dzw * a (x−n x d, y−n y d, z)w b (x−n x d, y−n y d, z)e ik Lx x+ik Ly y+ik Lz z . (15) One may therefore use the Rabi frequencies to control the magnitudes of the various J L and make them equal to J for all the beams except the onsite one, for which it should equal J 0 . In addition, the directions of the laser wavevectors should be set so that their wavevectors are in the yz plane with the appropriate projection onto the y axis to give the tunneling matrix elements along the x axis the phase specified in the reference Hamiltonian (13), exactly as in [76,77] 11 . In the rotating frame with respect to the V nxny and ω 0 terms in Eq. (14) one thus obtains the following Hamiltonian: nx,ny b nx,ny,qz + Je i2παny a † nx+1,ny b nx,ny,qz + Ja † nx,ny+1 b nx,ny,qz + h.c.
This is exactly the sum of Eqs. (5) and (6) for the specific reference Hamiltonian (13).
Integrating out the b fermions thus gives the desired outgoing Lindbladian, Eq. (7) with , the local density of states per unit cell due to the z direction motion.
As advocated in Ref. [77], this construction is best-suited to two-electron atoms, such as 171 Yb. One may use the clock transition states, that is, the ground state 1 S 0 , and the long-lived first excited state 3 P 0 , as the a and b states, respectively. Most of the parameters could be set as suggested there, with a couple of important differences: (i) In the system (xy) plane both the a and b states should occupy the same square lattice, hence one may use the magic wavelength to produce it (as assumed for simplicity in Fig. 4). However, one may employ any wavelength for which the polarizability of the two states has the same sign, which would allow to go further away from atomic resonances and reduce heating; (b) For the confinement in the z direction one should use the wavelength for which the b-state polarizability vanishes (which was the assumption in the previous discussion), or any wavelength for which the b-state polarizability is opposite to that of the a state, thus confining the latter but not the former atoms. The lattice modulation in the z direction could be created by any wavelength with non-vanishing b polarizability.
One may then set the 3D square optical lattice depth to be of the order of 10E R , where E R ∼ 2π · 1 kHz is the free space recoil energy of an atom due to the absorption or emission of a single photon at the relevant wavelength. The corresponding ordinary tunneling amplitudes would then be of order ∼ 5 · 10 −2 E R . Choosing the modulations ε x , ε y of order of a few E R would then strongly suppress ordinary tunneling in the xy plane. On the other hand, choosing the laser amplitudes about 3 times smaller than the lattice band spacing (to prevent exciting higher bands) could give optically assisted tunneling amplitudes of order J ∼ 10 −2 E R . As follows from Eq. (7) and the subsequent discussion, by reducing ε 0 , that is, choosing frequencies such that the residual kinetic energy of the b atoms is small, one may increase ν 0 = 1/[π ε 0 (4J z − ε 0 )] and thus the rates γ out . The only limitation is that γ out should be small enough with respect to E R so as not to mix in unwanted states such as higher bands of the optical lattice. Thus one may achieve γ out ∼ 2πν 0 J 2 ∼ 10 −1 E R . 12 The corresponding timescale would be a couple orders of magnitude smaller than the spontaneous emission time of the b state as well as the typical time of heating loss, which are above 1 sec in this parameter regime [77].
With this we saw how to realize the outgoing Lindbladian (7). If desired one may also similarly realize the incoming Lindbladian (8) utilizing a reservoir of atoms in a third state c, as described in Sec. 3. The added complexity with respect to the equilibrium case is offset by the efficiency of the resulting relaxation process: While cooling fermionic atoms to low temperatures so as to approach an equilibrium ground state is generally hard [33,34], here the relaxation is engineered to bring the system to the desired state at a finite rate ≥ γ in , independent of the system size.
To demonstrate the utility of this approach I plot in Fig. 3(b) the resulting spectrum of the single-particle steady-state density matrix (9) [the occupancies n λ (k)], which deviate from the ideal values of 0 and 1 by less than 4%. The only exceptions are the edge modes, whose occupancies vary continuously between 0 to 1. This is a consequence of the steady state being topologically nontrivial, as discussed in the next Section. The behavior is qualitatively similar to the equilibrium case at small but finite temperature, though with a non-Fermi-Dirac distribution (with respect to the reference Hamiltonian), cf. Eq. (9). Employing the Kapit-Mueller reference Hamiltonian [70] truncated to a finite range may allow for arbitrarily low values of the flatness ratio. For example, adding to the Hofstadter reference Hamiltonian (13) the small next-nearest-neighbor terms of the Kapit-Mueller Hamiltonian results in steady state populations deviating by less than 0.7% from their ideal values, as depicted in Fig. 3(c).
One may now introduce different perturbations, and study the stability of the state against them. The effects of two of them are presented in Fig. 5: (i) A deviation, ∆J 0 , of the onsite parameter J 0 [cf. Eq. (13)] from its optimal value; (ii) Finite Hamiltonian in Eq. (1), corresponding to residual nearest-neighbor tunneling with matrix element J . As expected, the finite minimal decay rate (≥ γ in ) of all the modes under the combined Lindbladian (7)-(8) stabilizes the topological state against perturbations, so long as they do not reach the order of magnitude of the unperturbed couplings.
6 Topological classification and detection out of equilibrium 6

.1 Topological classification
In what sense should the mixed state produced by the above procedure be considered topological? To answer this question I will mostly follow the approach put forward in Refs. [7-9, 13, 52], though using a somewhat different argument. Topological classification of the ground states of quadratic fermionic Hamiltonians [28][29][30][31] is based on spectrally flattening the first quantized Hamiltonian, i.e., on replacing the eigenenergies by their signs [with respect to the chemical potential µ 0 = (ε 1,max + ε 2,min )/2; see Fig. 2

(c)] and defining
This is a (first-quantization) operator that assigns the value −1 to all the occupied states, and +1 to all the unoccupied ones, and thus obeys Q 2 (k) = 1. Topological classification then proceeds as the study of the structure of the space of all possible such operators Q(k).
Here comes the crucial observation: In equilibrium at zero temperature the matrix representing Q(k) is simply related to the equilibrium single-particle reduced density matrix, Q µµ (k) = 1 − 2Tr(ρ ss a † kµ a kµ ) = 1 − 2G ss µµ (k). Out of equilibrium, for a quadratic Lindbladian of the type considered here, the steady state is still characterized by the steady state single-particle reduced density matrix G ss µµ (k), but its eigenvalues (the occupancies) are not strictly 0 or 1 [cf. Fig. 3(b)-(c)]. However, as long as there is a gap between the occupancies which are closer to 0 and those which are closer to 1 (akin to the "purity gap" discussed in Refs. [7-9, 13, 52] 13 ), one can extend the definition of Q(k) to this case as where n 0 can take any value in the purity gap (similarly to the chemical potential in equilibrium zero-temperature systems), e.g., n 0 = 1/2. so that it is again amenable to topological classification by the procedures of Refs. [28][29][30][31]. Another way to formulate this approach is by defining the "entanglement Hamiltonian" between the system and the baths, H ref ≡ − log ρ ss . In the current case it is quadratic in the fermionic operators, with the same eigenmodes |kλ and with eigenvalues log{[1 − n ss λ (k)]/n ss λ (k)}. As long as it is gapped, one may employ it to define the usual topological indices through its flattened version, which is exactly Q(k) given by Eq. (18). For paired states one would need to include anomalous averages (of two creation or two annihilation operators) in the definition of Q µµ (k) as well. Let us note that in the presence of disorder the quasi-momentum k is not well defined, but one may introduce instead periodic boundary conditions with phases φ x , φ y , . . . in the different spatial directions, and classify the dependence of Q on these phases [26,78].
This approach has the pleasant property of assigning to equilibrium noninteracting systems at finite temperature the same sharply-quantized topological index as at zero temperature, reflecting the fact that the protected edge modes of such systems are not modified as function of temperature (although response functions do change, of course). This should be contrasted with other recent proposals, such as treating a non-pure steady state as a mixture of pure states and averaging over the topological index of these states, which does not lead to a quantized value [48]. Another recent approach advocates characterizing the density matrix in terms of the Uhlmann phase, which is quantized but shows spurious abrupt changes at a particular finite temperature in equilibrium [49,50]. The topological indexes defined in this work are not only free from these problems, but, as I will show in the next subsection, are experimentally measurable in cold atom systems.
Continuing with our dissipative Hofstadter example, one may now calculate the integer topological index, the first Chern number [79], where the integration is over the first Brillouin zone. This gives the value 1 for the states depicted in Fig. 3(b)-(c), which implies that there are edge eigenmodes of G ss whose occupancies must straddle the gap in the spectrum of G ss , as indeed seen in Fig. 3(b)-(c). Let me also note that Eq. (10) shows that starting from, e.g., the trivial completely occupied [n λ (k) = 1 for all λ and k] or completely empty [n λ (k) = 0 for all λ and k] states, a gap in the spectrum of the occupancies will open up at arbitrarily short times, due to the gap in the spectrum of γ out . Thus, a nonzero Chern number and edge modes in the occupancy gap will immediately appear.
It should be noted that these edge modes are not decoupled modes of the dynamics. On the contrary, by Eq. (12) each mode in the system is damped by a rate γ out λ (k) + γ in ≥ γ in , 13 As noted in these works, in contrast with equilibrium zero-temperature systems, the purity gap might close by varying the Lindbldian continuously without closing the gap of of the Lindbldian. For example, in our system one may increase γ in to arbitrarily large values, and thus actually increase the gap of the Lindbldian, while driving the system towards the trivial completely occupied state, n ss λ (k) = 1.
as explained above, so any deviation from the steady state, either in the bulk or at the edge, would decay at a finite rate, independent of the system size. As a matter of fact, since by Eq. (7) γ out λ (k) = 2πν 0 [ε ref λ (k)] 2 , the decay rates of the edge modes span the gap between the decay rates of the almost-filled (λ = 1) and almost-empty (λ ≥ 2) bands.
One can shed further light on these results from a different perspective. There are different ways to define topological order. Some of them depend on the existence of anyons, manifested by topological degeneracy and entanglement entropy, and thus exclude integer quantum Hall or Chern insulators [67]. Others define a topologically-ordered state as one that cannot be created, starting from some trivial product state, by a finite depth quantum circuit, or, equivalently, finite time (not scaling with the system size) evolution with a local Hamiltonian [80]. The last definition is usually taken to include Chern insulators, due to their gapless edge modes, which thus have long-range (power-law) correlations that cannot be created at a finite time by a local Hamiltonian. In our system the dynamics is dissipative but still local, hence one would expect similar scaling of the relaxation time to reach a topologically-nontrivial state. Indeed, such an approach has recently been suggested for the definition of topological order out of equilibrium [56]. Yet we found that pure Chern insulator states can be dissipatively created exponentially fast in time at a finite rate (finite gap of the Lindbladian) by a local (though not finite-range) Lindbaldian, which seems to be at odds with this argument. The resolution comes from the topological edge modes of G ss . Their existence implies that one can only achieve the desired pure state in the bulk of the system. On the edge there must be edge modes of the density matrix whose average occupancies span the gap between 0 and 1 [ Fig. 3(b)-(c)]. In a qualitatively similar manner to a finite temperature equilibrium state, this causes the edge correlations in the steady state of our system to decay exponentially with distance, so indeed there is no issue with creating these correlations at a finite rate by local Lindblad dynamics. In particular, this shows that while local Hamiltonian evolution of a pure state cannot to change its Chern number [81][82][83][84], local Lindblad evolution can, implying that Chern insulators, and actually all noninteracting topological states in any dimension, are not topological by the definition of Ref. [56].

Detection of nonequilibrium topology
How could the topologically-nontrivial states be detected experimentally? Here too one may employ many of the techniques suggested in equilibrium [35][36][37][38][39][40][41][42][43][44]. I will examine some of them concentrating on the dissipative Hofstadter state as an example. It should be noted that, differently from equilibrium, here the edge modes do not propagate but rather display pure decay, due to the purely-dissipative dynamics. Since there is no Hamiltonian term in the master equation (1), the conductivity cannot be straightforwardly defined. Hence, other indicators are needed.
One approach is to probe the real-space distribution of the atoms [85]. As shown in Fig. 6(a), it attains the constant value of 1/q per lattice site (corresponding to the lowest band, λ = 1, of the reference Hamiltonian (13) being almost filled, and all the others, λ ≥ 2, being almost empty) in bulk of the system, in analogy with the incompressibility of the corresponding equilibrium state. The partially-filled edge states [see Fig. 3(b)-(c)] could be revealed by taking the derivative of the density map with respect to the filling rate γ in , which is significant only near the edge [ Fig. 6(b)]. However, this is not an unambiguous indicator of a topologically-nontrivial state.
Different information can be revealed by examining the quasi-momentum distribution function [ Fig. 7(a)], which can be inferred from a time-of-flight measurement, i.e., releasing the atoms and taking images of the expanding cloud [33,34]. The q-fold structure in the k x direction reflects the periodicity of the Hofstadter reference Hamiltonian (13), which is Figure 6: (a) A color map of the real-space particle distribution of the dissipative Hofstadter model on a 50 × 50 lattice with in a circular trap [simulated by adding a large state-independent evaporation rate outside the trap to L out , Eq. (7)]. Inside the trap the density ≈ 1/q = 1/5 atoms per unit site, since only the lowest band of the reference Hamiltonian (13) is filled. The parameter values are as in Fig. 3. (b) A color map of the derivative of the local particle number with respect to γ in /(2πν 0 J 2 ). It is significant only near the edge, thus revealing the edge modes from Fig. 3.
hidden in the real-space distribution, Fig. 6. To find the Chern number itself, Eq. (19), one may take a hybrid time-of-flight image [86] (see also [87]), in which the atoms are released in the x direction but are kept confined in the y direction, leading to a hybrid real space-momentum space distribution [ Fig. 7(b)]: A nonzero Chern number, C 1 = 0, implies that the spatial density profile in the y direction winds by C 1 unit cells of the reference Hamiltonian (13) (C 1 q lattice sites) as k x winds through the Brillouin zone. Thus, the mixed state Chern number defined above is sharply-defined and measurable in an experiment. Let me note (following the discussion in the previous subsection of mixedstate topological classification) that for an equilibrium system at finite temperature this mixed-state Chern number will retain its value all the way up to infinite temperature, as only then would the contrast between the maxima and minima in Fig. 7(b) disappear.

Conclusions
To conclude, this work exposes the capabilities and limitations of dissipative preparation of topological states using quadratic Lindblad dynamics. We have seen that while finiterange gapped quadratic Lindbladians cannot lead to an exponential decay in time at a finite (system size independent) rate towards a unique pure steady state with nonzero Chern number, they can lead to a mixed state arbitrarily close to the desired pure one, which could be realized with cold atoms using currently available experimental techniques. The pure limit may be achieved if exponentially-local Lindbladians are allowed. I have also discussed the topological classification of such states, using the single-particle reduced density matrix, and the implications of topological nontriviality, such as the edge modes of the single-particle density matrix and the winding number of mixed position-momentum density maps.
While this study concentrated on Chern insulators, that is, 2D systems in symmetry class A, the extension to other topological classes [25][26][27] is clear. In particular, as mention Figure 7: (a) A contour plot of the "quasi-momentum" (Fourier transform defined, here only, with respect to a 1 × 1, rather than 1 × q = 1 × 5, unit cell) distribution of the dissipative Hofstadter model on the first Brillouin zone, which can be extracted from a time-of-flight measurement [33,34]. The q = 5 fold periodicity in the k x direction reflects the periodicity of the Hofstadter reference Hamiltonian (13). The parameter values are as in Fig. 3. (b) A color map of the mixed quasi momentum-real space distribution (along the x and y directions, respectively), which can be inferred from a hybrid time-of-flight measurement [86]. The spatial density profile in the y direction winds by one unit cell (of size q = 5) as k x winds through the Brillouin zone, showing that C 1 = 1 [cf. Eq. (19)]. The system is a 50-site-wide strip whose long axis is x; similar results hold in the presence of a trap [86]. k x and k y are here dimensionless, or, equivalently, measured in units of 1/d, the inverse 2D lattice spacing.
in Sec. 4, the no-go theorem holds for these as well, following the no-go theorem for finiterange Hamiltonians with flat bands in topologically nontrivial states in two and higher dimensions [66]. One may of course argue that in the dissipative context, the limitations imposed by, e.g., time-reversal symmetry, lose their physical significance. If these are relaxed, pure steady states can result from finite-range gapped Lindbladian. This work was centered at cold atom systems, but the ideas introduced here should find applications to other types of "quantum simulators", including, among others, superconducting nanocircuits, quantum dots, trapped ions, and nuclear and impurity spins [32].
These results give rise to many important questions for future study. Among them is the effects of including the Hamiltonian in the Lindblad dynamics (1). While this neither modifies the results of the no-go theorem nor is necessary for our recipe as we have seen above, it may allow one to define currents and hence, e.g., the Hall conductivity. Our recent results [88] indicate that the latter is not generally quantized in the drivendissipative case. A similar conclusion holds for Laughlin-type charge pumping [89]. Thus, the equilibrium relation between topology, edge modes, and quantized responses does not generalize unmodified to nonequilibrium situations, an interesting topic for further research.
Even more exciting directions open up when going beyond quadratic dynamics, thus allowing for the creation of strongly-correlated dissipative states. In the topological context, states which were distinct in the noninteracting case may become topologically-equivalent (first realized in 1D [90,91]), while new phases could arise [19,[58][59][60][61][62]67,80,92,93]. In particular, above 1D the inclusion of interactions may lead to "topologically-ordered states" (e.g., fractional quantum Hall), whose excitations display fractional charge and statistics (anyons). The statistics could even be nonabelian, raising the possibility of topological quantum computation [94]. It would be interesting to understand the dissipative analogues of these, and the even more exciting possibility of novel topological phenomena which occur only out of equilibrium.
representation, ρ I S = Tr E ρ I SE : ∂ t ρ I S (t) = −iTr E H I SE (t), ρ I SE (t 0 ) − t t 0 dt Tr E H I SE (t), H I SE (t ), ρ I SE (t ) .
Assuming that at t 0 the system and the environment are uncorrelated, ρ I SE (t 0 ) = ρ I S (t 0 ) ⊗ ρ I E (t 0 ), the first term on the right hand side of the last equation vanishes, since it reduces to a sum over product of averages of single fermion operators over the system and environment initial states. As for the second term, it is second order in the system-environment coupling. Thus, to the lowest Born approximation one may take a factorized form for the system+environment density matrix, and further assume the environment state is approximately constant in time (i.e., the baths are "large", hence not affected by the system), ρ I SE (t) ≈ ρ I S (t) ⊗ ρ I E . Moreover, one may make the Markov approximation, namely assume that the bath dynamics is much faster than that of the system, so for the values of t which significantly contribute in the last equation one may approximate ρ I S (t ) ≈ ρ I S (t). This also allows to take the lower integral limit to −∞. We can finally write The right hand side then breaks into a sum of averages of two environment operators over the environment state, times a term containing two system operators.
Let us now specialize to our quadratic fermionic system: where η → 0 + is a convergence factor, and P denotes Cauchy's principal value. Twice the real parts of these coefficients, γ out p ≡ 2Re(Γ out p ) and γ in p ≡ 2Re(Γ in p ), are the dissipation rates (rates of particles leaving the system into the bath or entering the system from the bath, respectively) as given by the Fermi golden rule. The imaginary parts, ∆ε out p ≡ Im(Γ out p ) and ∆ε in p ≡ −Im(Γ in p ), are corrections to the system eigenenergies ("Lamb shifts"). Returning to the Schrodinger picture and to the original basis, we obtain the fermionic Lindblad master equation for the system density matrix ρ S (which is denoted by ρ in the rest of the paper), Eq. (2), withh S mm = p (ε S p + ∆ε out p + ∆ε in p )(u S pm ) * u S pm being the Hermitian system Hamiltonian modified by the Lamb shift corrections, while γ out mm = p γ out p u S pm (u S pm ) * , γ in mm = p γ in p (u S pm ) * u S pm are Hermitian positive semidefinite decay rate matrices.
As noted in Sec. 2 above, this derivation can be extended to the case when fermion number in the system+environment is not conserved due to the presence of pairing terms (a m a m , a m b n , b n b n and their conjugates) in the total Hamiltonian (20). The end result would be (beyond modifications to the expressions for the previously-defined coefficients) that the Hamiltonian part of the Lindblad equation would be supplemented by new terms of the form ∆ mm a m a m + h.c.. Similarly, the Lindbladian would include new terms of the form γ pair mm [2a m ρ S a m − a m a m ρ S − ρ S a m a m ]/2 + h.c. with antisymmetric matrices ∆ mm and γ pair mm .