Multidimensional dark space and its underlying symmetries: towards dissipation-protected qubits

Quantum systems are always subject to interactions with an environment, typically resulting in decoherence and distortion of quantum correlations. It has been recently shown that a controlled interaction with the environment may actually help to create a state, dubbed as {\em ``dark''}, which is immune to decoherence. To encode quantum information in the dark states, they need to span a space with a dimensionality larger than one, so different orthogonal states act as a computational basis. We devise a symmetry-based conceptual framework to engineer such degenerate dark spaces (DDS), protected from decoherence by the environment. We illustrate this construction with a model protocol, inspired by the fractional quantum Hall effect, where the DDS basis is isomorphic to a set of degenerate Laughlin states. The long-time steady state of our driven-dissipative model exhibits thus all the characteristics of degenerate vacua of a unitary topological system. This approach offers new possibilities for storing, protecting and manipulating quantum information in open systems.

Quantum systems are always subject to interactions with an environment, typically resulting in decoherence and distortion of quantum correlations. It has been recently shown that a controlled interaction with the environment may actually help to create a state, dubbed as "dark", which is immune to decoherence. To encode quantum information in the dark states, they need to span a space with a dimensionality larger than one, so different orthogonal states act as a computational basis. We devise a symmetrybased conceptual framework to engineer such degenerate dark spaces (DDS), protected from decoherence by the environment. We illustrate this construction with a model protocol, inspired by the fractional quantum Hall effect, where the DDS basis is isomorphic to a set of degenerate Laughlin states. The long-time steady state of our driven-dissipative model exhibits thus all the characteristics of degenerate vacua of a unitary topological system. This approach offers new possibilities for storing, protecting and manipulating quantum information in open systems.
It is believed that dissipation conspires against coherence of quantum states, rendering them to be close to a classical ensemble. This belief was recently challenged by approaches aimed at incorporating both drive and dissipation to reach a correlated coherent steady state [1][2][3][4][5][6][7][8][9][10][11] . One remarkable example has been the idea of harnessing dissipation to purify non-trivial topological states [12][13][14] . This is achieved by a careful interplay between radiationinduced drive, and coupling to an external bath, that provides a desired relaxation channel. A sequence of excitations and relaxations generates, in the long time limit, a non-equilibrium steady state that decouples from the external drive creating a decoherence free subspace 15,16 dubbed a dark state. This idea opens a way to engineer a rich variety of non-trivial stationary states, going well beyond thermal states of equilibrium systems.
In order to use this approach to design (and ultimately manipulate) qubits, it is necessary to engineer a non-equilibrium steady space, which is at least twodimensional 14,17 . Here we develop a framework to construct driven dissipative schemes with degenerate dark spaces (DDS). We achieve this goal by analyzing the role of symmetries in dissipative dynamics. Specifically, we claim that the dimensionality of DDS is given by the period of the projective symmetry representation 18 , inherent to the system's evolution (which is considered to be Lindbladian 19,20 ). To this end we extend the discussion of Lindbladian symmetries [21][22][23] to include those that are realized projectively, providing a link between the projective representations and the dimensionality of the DDS density operator.
We illustrate this framework by studying driven dissipative evolution of a correlated one-dimensional (1D) system, inspired by Laughlin quantum Hall states with ν = 1/m filling fractions (m is an odd integer) in a quasi-1D strip (the so-called "thin torus limit"). This evolution is described by a Lindbladian master equation that possesses a DDS. The latter is spanned by m orthogonal vectors, isomorphic to the set of many-body Laughlin ground states on the torus. This correlated DDS has an extra advantage of being exactly described by computationally convenient matrix product states (MPS). We design a systematic protocol, based on adiabatically varied Lindbladians, that maximizes the purity and fidelity (overlap with the dark space) of its ultimate steady states.
As a warm up for the symmetry discussion, let's consider the Hamiltonian case. If a Hamiltonian is invariant under the action of a symmetry group G, then action of the group elements, g ∈ G, on a state is implemented by a unitary representation 18 . In particular, the eigenstates of the Hamiltonian can be labeled by eigenvalues of the symmetry operator. As quantum states form rays in the Hilbert space, such that states |ψ and e iφ |ψ are equivalent, it is natural to consider representations that satisfy the group multiplication rule up to a phase, i.e. projective representations 18 , defined as D(g 1 )D(g 2 ) = e iφ(g1,g2) D(g 1 g 2 ), (1) where D(g) is a representation of a group element g ∈ G.
Every projective representation is characterized by the set of phases ω 2 (g 1 , g 2 ) = e iφ(g1,g2) , known as a 2-cocycle, which are strongly constrained by associativity 18 . For a quantum system invariant under the projective representation (1), the period of the 2-cocycle (i.e the minimum number m such that [ω 2 (g 1 , g 2 )] m = 1 for all g 1 , g 2 ) determines the dimension of the degenerate space. Given a non-trivial 2-cocycle, representations of at least some group elements do not commute, even for an Abelian group G, i.e. for some g, h ∈ G, [D(g), D(h)] = 0. One can thus label the eigenstates by eigenvalues of say D(g) and generate a different state, with the same energy by acting on it with D(h). Notice that this argument implies degeneracy of an entire spectrum. Degenerate dark states in Lindbladian evolution.-The way symmetry operators act on the Lindbladian operators is rather different from the Hamiltonian case. In a system with a combined unitary and dissipative dynamics, the most general Markovian evolution of the density matrix operator, ρ, is described by the quantum master equation,ρ = L(ρ), with Here H is an effective Hamiltonian that represents the unitary evolution. Generically non-Hermitian, quantum jump operators, i , describe environment-induced dissipation effects 19,20,24 . A Lindbladian is invariant under an irreducible unitary representation D(g) with g an element of G, if the Hamiltonian and the quantum jump operators satisfy 21,22 (also known as weak symmetry 21 ) where U (g) is a unitary matrix that depends on g. In particular, if σ is an eigenmatrix of L with an eigenvalue λ σ , i.e. L(σ) = λ σ σ, then D(g)σD † (g) is an eigenmatrix with the same eigenvalue. These eigenmatrices obtained from σ by conjugation can represent either the same, or different states. For a projective representation D(g), satisfying (1), D(g)σD † (g) and σ are necessarily different for some element g. This is because if for a particular element h ∈ G, D(h)σD † (h) = σ, then σ and D(h) share eigenvectors, so one can take an element g such that D(g) does not commute with D(h). Given that [D(g), D(h)] = 0, D(g)σD † (g) and σ do not share eigenvectors, meaning that they are different. By virtue of the Schur's lemma 18 , the only case where this logic fails is for σ a fully mixed state, proportional to the identity matrix. Focusing on the case of projective representations of Abelian groups, where D(h)D(g) = e iφ(h,g) D(g)D(h) (hereφ(g 1 , g 2 ) = φ(g 1 , g 2 ) − φ(g 2 , g 1 ) is again a cocycle) one can determine the dimension of the degenerate subspace in terms of the factor e iφ(h,g) . Focusing on an eigenvector |e of D(h) with eigenvalue e iα , the operator D(g) acts on this state as a cyclic raising operator, as D(h)D(g)|e = e i(α+φ(h,g)) D(g)|e . Forφ(h, g) = 2πa m , with a, m co-prime integers, one can raise a state m times before it returns to itself. This implies that all eigenspaces of L are m-fold degenerate. In particular, the stationary subspace, defined by the eigenvalue λ = 0, is m-dimensional. Note that the degeneracies of the eigenspaces of L do not translate into degeneracies of the density matrix in general, as the the eigenmatrices of L need not to be self-adjoint, but can come in pairs {σ, σ † } with eigenvalues {λ σ , λ * σ }. For the DDS, this is not a problem as λ = 0.
Although conditions presented above are sufficient for the existence of DDS, they are actually excessive and impractical. Indeed, the symmetry operators D(g) split the entire Hilbert space into sectors of different quantum numbers that do not mix during the evolution. In other words, there is a certain number of conservation laws, which confine the long time evolution of any initial state to only a limited fraction of DDS. To access the entire DDS, this needs to be avoided. In order to achieve this, we consider states ρ in the DDS that satisfy We call such DDS frustration free, as ρ is annihilated individually by each quantum jump operator, i . For systems that satisfies Eq. (4), one can deform quantum jump operators in a way that the symmetry is broken in all the decaying subspaces, while it is maintained within the DDS. With this in mind, we define dressed quantum jump operators, that do not satisfy Eq.
i are for now arbitrary local operators. Dynamics, generated by the dressed operators, does not, in general, obey any conservation laws. Yet,˜ i still satisfy Eq. (4), which preserves the DDS and its degenerate multidimensional nature. For properly dressed quantum jump operators, a generic initial state evolves into a state within DDS, which has projections on all of its basis eigenmatrices.
Below we demonstrate these considerations on a 1D model, borrowing intuition from the well-studied physics of the Laughlin states on a torus 25,26 . In particular, we demonstrate that the system is driven to DDS regardless of the nature of an initial state (pure or mixed). We also devise adiabatic time-dependent Lindbladians that guarantee that the initial state is fully driven into the DDS, resulting in a state with a maximized purity.
Laughlin states in a narrow torus geometry.-A quantum Hall droplet of N electrons subject to a magnetic flux N Φ = mN (in units of the flux quanta) and filling fraction N/N Φ ≡ ν = 1/m (with odd integer m) develops nontrivial correlations that are reproduced by Laughlin wavefunctions 27 . In a torus with periods L x and L y (with distances measured in units of the magnetic length), the area is related to the flux as L x L y = 2πN Φ . The Laughlin states at filling ν = 1/m correspond to exact zero energy states of a local Hamiltonian 28 , which, after projecting onto the lowest Landau level (LLL), takes the form H = n † 0,n 0,n + † 1,n 1,n , where the operators s,n (s = 0, 1) are 29 (see supplemental material) Here c n destroys an electron at orbital n; η(l) ∝ e −(κl) 2 is a fast decaying function in the narrow torus limit, κ 2 = 2π NΦ Lx Ly 1, (see supplemental material). The crucial property of Laughlin states |Ψ , that makes them useful for our discussion of the Lindbladian evolution, is that they satisfies s,n |Ψ = 0, for s = 0, 1 and all n.
In the quantum Hall context, the operators D(g) of Eq. (3) correspond to inserting fluxes through the two periods of the torus, Fig. 1. In the 1D representation they are the translation operator T and the operator U , that measures the center of mass of the particles in orbital space. They are given by wheren l andn k are the number operators at position l and with momentum k, correspondingly. Note that T and U satisfy T U T † = e −2πiν U . They thus provide a projective representation of the group G = Z m ×Z m with m = 1/ν. This group is represented by One can choose a basis |Ψ a , such that e.g. operator U is diagonal, with its eigenvalues, ω a , along the diagonal. In this basis, the operator T acts as a raising operator, because, if U |Ψ a = e i 2πa m |Ψ a , then U T |Ψ a = e i 2π(a+1) m T |Ψ a . In the context of the quantum Hall effect, states |Ψ a are the m-fold degenerate Laughlin ground states on the torus 25,30 .
Note that, although the construction of the Lindbladian using the projective symmetry ensures that the eigenvalue λ = 0 of L is m-fold degenerate, the frustration-free condition Eq (4) enlarges the degeneracy of λ = 0 to m 2 . This is seen as follows: The symmetry operators U and T ensure that the density matrices |Ψ a Ψ a | for a = 0 . . . m − 1 share the same eigenvalue (they are related by conjugation with T ). In general the matrices |Ψ a Ψ a+p | are related, for a fixed p, by conjugation with T . Now these m different families (each labeled by p = 0 . . . m − 1) share the same eigenvalue between them due to the frustration-free condition which ensures that if |Ψ a is annihilated by the quantum jump operators, then all the matrices |Ψ a Ψ b | are as well.
Lindblad operators from Quantum Hall physics.-In the narrow torus limit, κ 1, one can truncate expressions for operators s,n in Eq. (5), which become short-range in n. In this limit we have 31 where β = η( 3 2 )/η( 1 2 ) = 3e −2κ 2 . These operators transform as U s, s,j and T s,j T † = s,j+1 . Hereafter we regard β as an arbitrary parameter. We now employ these results to construct quantum jump operators that drive the system into the frustration free DDS. In contrast with the s,i operators in the previous section, here these operators represent processes in a real lattice of N Φ sites, where c i destroys a fermion at the lattice site i. We assume m = 3 and the fermion density ν = 1/3.
We also assume a purely dissipative evolutionρ = L(ρ) with H = 0, Eq. (2), and the quantum jump operators Here operators Q i are linear combinations of the operators that annihilate the Laughlin states. Parameters t, B, A and β are determined by a realization of the dissipative dynamics and are non-universal. The DDS only depends on β. For the purposes of this work we take them as free parameters. Realization of Lindblad operators˜ i .-The dissipative evolution dictated by the Q i operators can be obtained by coupling two systems: a target fermion chain (where c i destroys a fermion at site i) with no intrinsic dynamics (H 0 = 0) and a fermionic chain with Hubbard interactions described by the Hamiltonian

fermion in this at position i in this chain and
We assume that the chemical potential E 1 is the largest energy scale in the system. These two chains interact through an external classical radiation, with cou- where Ω is the intensity of the radiation, α parameterizes the spatial laser envelope, and ω is the frequency of the monochromatic light. The role of the driving is to excite particles from the target chain to the interacting chain. The particles then relax to the lower energy state interacting with the bath, that provides dissipation. These components are shown in Fig. 2a As a motivation for this construction, let's consider how a charge-density-wave (CDW) state decouples from the dynamics in this context, in the limit J 1 << U and α = 0, becoming a dark state. At first order in Ω, exciting a single particle from the c to the f chain is strongly suppressed as E 1 is large. But in second order ( , the radiation Hamiltonian can create two states in the f chain, which can bound in a doublon, consisting on a tightly bound pair of fermions (due to the Hubbard attraction). The wavefunction of the doublon decays exponentially with the distance d between its two constituents as t d with t ∼ J 1 /U 1. This means that, for the laser to create a doublon, the particles in the c chain should be near each other, as the laser acts locally. In particular, the state which locally contains nearby fermions c † i c † i+1 |0 will be affected by radiation, as well as c † i c † i+2 |0 , where |0 is the state with no particles. The first local configuration that is not affected by radiation, becoming dark is , as the fermions are too far away to be excited into a configuration with a non-vanishing matrix element with the doublon. The whole system will decouple when reaches one of the CDW states. Increasing the range of the laser (by letting α = 0), creates superpositions.
The Lindblad operators (8) are obtained considering transitions between the doublon band and the low energy band in the system (more details in the supplemental information). After performing the rotating wave approximation to account for the time dependence of the radiation Hamiltonian, the dynamics of the system occurs between the lower c, and the doublon bands (see Fig 2b upper panel). Using second order perturbation theory in Ω, we obtain the matrix elements of the transitions between the doublon states |d i = f † i f † i+1 |0 and the lower band states |i, j = c † i c † j |0 , The transition process from lower band to doublon reads (after adiabatic elimination of the while A, B, and β entering the definition of Q satisfy A = α, B = αA Ω and β = α 2 . The fermion operators in Q i all act on the chain c. Finally, taking into account the transitions from the doublon back to the c chain, which is mainly mediated by the dissipation with the bath, and integrating out the doublon states, we arrive at the Lindblad operators (8).
Structure of the DDS.-Heuristically, one can understand the roles of R † i and Q i as follows. Operator Q i checks if at the site i the state matches the local configuration of one of the Laughlin-like states. If true, it gives zero and the system stops evolving locally; if false, the operator R † i Q i scrambles the particles. As long as this process can efficiently mix the particles locally, all states in the Hilbert space eventually evolve into DDS, spanned by the three Laughlin states. Crucially, the decay into these states is a consequence of the projective symmetry, enforcing existence of the degenerate space with Q i ρ = 0 and thusρ = L(ρ) = 0.
The basis of such DDS is formed by the Laughlin states where N is a normalization factor, a = 0, 1, 2 and The state |• • • i represents the three consecutive empty sites at positions (3i − 2, 3i − 1, 3i), while a full dot represents an occupied site, e.g.
The dark space basis vectors, |Ψ a , are related by a translation T by one site, as T |Ψ a = |Ψ a⊕1 , where ⊕ is an addition modulo 3. In the basis |Ψ a , the operators T and U are represented by the 3 × 3 matrices Within DDS the density matrix is ρ = 2 a,b=0 with ab a 3 × 3 positive semidefinite, hermitian matrix with unit trace. In general, the structure of the density matrix within DDS depends on initial conditions, which determine the parameters ab . The basis vectors that define the DDS depend explicitly on the parameter β. For β = 0 the DDS is spanned by a mixture of three different classical CDW 30 , e.g. | • • • • • • . . . , with periodic density and sharp structure factor in each basis vector. Changing this parameter modifies the average local density. The latter is given by indicating that, if β is known, the local density measurements are enough to determine the probabilities ρ aa . Starting from a CDW configuration, the system is evolved using a protocol with a given β. This generates a DDS state that depends on β. The top row represents the pure state 00 = 1, while the second row represents the mixed state 00 = 0.5, 11 = 0.3, 22 = 0.2. Different colors are used to help track the changes in average occupation at each site and to highlight the 3-site periodicity of the density. b.-Static structure factor S1(k) for different values of β, for a pure state with ρ11 = 1. At β = 0 the system is in the CDW state, with a definite spatial periodicity, indicated by the peaks in S1(k) at k = 0, 2π 3 , 4π 3 . Increasing |β|, the system becomes more homogeneous. c.-Entanglement entropy of the DDS as a function of β.
An alternative way of characterizing the DDS is through the correlation function of local observables. To highlight the relation with the CDW states, we focus on the static structure factor shown in Fig. 3b. As β increases, the system transitions from a crystal-like state with the well defined translational symmetry of three sites, into a more homogeneous state, where the density is uniform across the system.
Finally, to describe the quantum nature of the DDS basis states, we compute their entanglement entropy. We separate the degrees of freedom of the system in two complementary large regions A and A c and define the partial density matrix ρ A = Tr A c (|Ψ a Ψ a |). The entanglement entropy is then S(β) = −Tr(ρ A ln(ρ A )). The result is shown in Fig. 3c. We observe that the entanglement entropy is monotonic with β. It reaches its maximum value of 2 ln(2) for MPS of bond dimension 2 at |β| → ∞.
Time evolution and global diagnostics.-Now that we have constructed a dissipative evolution that drives the system into the DDS, we discuss how the system approaches the DDS. We analyze the Lindbladian evolution with the quantum jump operators (8) numerically. The decay into the DDS is evaluated using a quench protocol: the system is initiated in the CDW state | • • • • • • . . . , which is one of the dark states of the Lindbladian at β = 0. This state is evolved then using the Lindblad operators (8) with A = B = t = 1 and β ∈ [0, 1] for simplicity (the results are qualitatively similar for slightly varying these parameters). In order to obtain the evolution of the system we perform exact diagonalization using Runge-Kutta (RK) integration 35 of the master equation, for systems of sizes up to L = 15.
To characterize the steady state mixture, we compute the purity of the state, defined as γ(t) = Tr{ρ 2 (t)}. From  Fig. 4a, we find that the purity approaches 1/3 with larger system sizes. This is indeed the case for a sudden quench from β(t ≤ 0) = 0 to β(t > 0) = 1. In this scenario, the system explores an extensive portion of the Hilbert space, becoming highly mixed, as seen in the intermediate region of Fig. 4a, where the purity plateaus at a minimum. Only after the system is sufficiently mixed, it starts approaching DDS and its purity increases. The information about the initial state is practically lost in the intermediate mixing process, and the eventual steady state is a highly mixed state within DDS.
Convergence to DDS, spanned by Laughlin-like dark states {|Ψ a } a=0,1,2 , may be visualized by where P DDS = 2 a=0 |Ψ a Ψ a | is a projector onto the DDS. Figure 4b shows that the system indeed evolves towards the Laughlin-like DDS, proving that this is the only non-decaying subspace of the Lindbladian evolution. At large times, one finds 1 − D DDS (t) ∝ e −λ0t , where the rate λ 0 is given by the lowest non-zero eigenvalue of the Lindblad operator.
One notices that λ 0 is slowly decreasing upon increasing the system size. We note that local observables, like the particle density, fast approach their steady state values in a way which is independent of the system size, Fig.4a. This separation of scales indicates that, while locally the system reaches a configuration that is close to the dark states that span the DDS, globally it takes much longer to fully reach the DDS.
If instead of quenching the system into β = 1, we quench it into β 1, we observe a very different be- havior. Here the system does not have to explore an extensive part of the Hilbert space before it reaches the DDS. As a result, the purity remains close to 1 at all times, as can be seen in Fig. 4c. Adiabatic evolution.-Although the previous analysis shows that the system does not generically end up in a pure state, it is possible to increase the purity of the final state by performing an adiabatic evolution from a pure state [37][38][39] . To illustrate this, we evolve the system from an initial state given by a superposition of the three CDW configurations. This allows us to characterize the coherences in the MPS basis throughout the adiabatic evolution (sec. III in SI). Individual CDW states can be created using existing experimental techniques 40 . We then evolve this state with the Lindblad operators (8) using a time-dependent β parameter: β(t) = ∆ · t for 0 < t ≤ 1/∆ and β(t) = 1 for t > 1/∆, where ∆ is the ramp velocity. The purity of the final state depends on ∆ as shown in Fig 5a. For small enough ∆, the system does not explore the whole Hilbert space, but instead remains almost pure throughout its entire evolution. This mechanism can be used to achieve a purity arbitrarily close to unity. For larger ramp velocities, the system rapidly departs from the initial state, exploring the many-body Hilbert space, as shown for intermediate times in Fig.5b, before leaking back to the DDS, which remains the only attractor of the dynamics. This increases the departure of the steady state from a pure state and erases the information about the initial state. The steady state purity as a function of the ramp velocity is shown in Fig.5c. Conclusions.-We have shown that to achieve a DDS the Lindbladian evolution should have an underlying symmetry, admitting a projective representation. The period of its 2-cocycle determines the dimensionality of the dark space. Reaching a DDS, protected against environmental influence, offers a way of maintaining quantum information. To manipulate this information, it is necessary to have an access to high purity states within the DDS. We found that an adiabatically varying Lindblad operator allows to reach such nearly-pure, entangled states. We have demonstrated these ideas by studying the thin torus limit of the ν = 1/3 fractional quantum Hall state of matter. Being able to generate and manipulate states within a DDS may be utilized for quantum information processing platforms. The many-body nature of the state renders it less fragile against local disturbances.
Acknowledgments.-We are indebted to R. Fazio for valuable discussions. R.S. acknowledges funding from by EPSRC grant EP/M02444X/1, and the ERC Start- In this section we provide supplemental information not been inserted in the main text. We discuss the analysis of the Lindbladian gap of the dissipative model with and without perturbations, as well as an expanded analysis of the adiabatic evolution. We revisit the map of the Laughlin state from two dimensions to the orbital basis, where it becomes effectively one dimensional. Finally we include a longer discussion about the experimental realization of the Lindblad operators.
The dissipative evolutionρ = L(ρ) withĤ = 0 studied in the main text is defined by the quantum jump operators In Sec. (I) we show our results for the Lindbladian gap in finite system sizes with no perturbations, while the case of the dissipative evolution with imperfections are presented in Sec. (II). In Sec.(III) we show details of the adiabatic evolution analysis.

I. LINDBLADIAN GAP IN FINITE SYSTEM SIZES
In this Section we show our analysis for the dissipative gap of the Lindbladian for finite system sizes. We obtain the gap in two different forms: (i) directly by exact diagonalization of the Lindbladian superoperator, or (ii) indirectly by the asymptotic decay rate (ADR) of the quantum state dynamics. While exact diagonalization allows us to study the Lindbladian gap for sizes up to L ∼ 12, the asymptotic decay rate analysis allows the study of larger L ∼ 15 system sizes.
The exact diagonalization analysis is performed by first describing the Lindbladian superoperator as a linear operator in the extended Hilbert space, and then diagonalizing the vectorized Lindbladian |L . The eigenvalues λ j of the Lindbladian have non-positive real values, with the DDS described by those with zero eigenvalue. The gap of the Lindbladian corresponds in this way to the eigenvalue with largest nonzero real part. Ordering the eigenvalues according to their real part (λ j ) ≥ (λ j+1 ), for j = 1 to the extended Hilbert space dimension, the dissipative model of this work have in this way λ j = 0 for j = 1, ..., 9 and the gap is described by the eigenvalue λ 10 . On the other hand, one can also extract the gap by studying the dynamics of the quantum state. In the long time limit only the slower decaying modes of the Lindbladian are relevant to the dynamics and the expectation value of any observable O(t) = T r(Ôρ(t)) is approximated by where λ ADR < 0 corresponds to the slower decay mode (or asymptotic decay mode) for the observableÔ. In Fig.(S6) we show our results for the Lindbladian gap and ADR analysis of the DDS(t) dynamics. We observe that the gap from exact diagonalization matches the one obtained from ADR. Furthermore we see that the dissipative gap increases going from 9 to 12 sites, and decreases again for 15 sites (to slightly the same value as 9 sites), which is the limit size that we can achieve. Although instructive, these numerical results for small systems do not allow us to unequivocally determine the nature of the dissipative gap in the thermodynamic limit. It is suggestive, however, for the presence of either (i) a gapped Lindbladian in the thermodynamics limit, or (ii) a gapless Lindbladian with a slow decaying of the gap with system size, excluding e.g. the possibility of an exponential closure of the gap with system size which would be detrimental for the preparation and manipulation of DDS in quantum information tasks.
In the main text, we showed that the existence of the symmetry operators U and T ensures that all of the eigenvalues of the Lindbladian are at least m-fold degenerate. It is good to stress that the eigenmatrices of the Lindbladian do not translate directly into density matrices. In particular, this means that the symmetry operators U and T do not generate directly a symmetry on the density matrices. The Lindbladian is a linear operator acting on the extended Hilbert space V = H ⊗ H, Eq. S16 which for concreteness we can take to be spanned by the vector states |i |j , each defined within a copy of H . The vectorized density matrices form a subspace S ⊂ V which is not a vector space, as the sum of two elements of S does not generically belong to S (e.g. given ρ 1,2 two positive definite operators with unit trace, the linear combination αρ 1 + βρ 2 with α, β complex numbers is not necessarily positive definite).

II. LINDBLADIAN PERTURBATIONS
In this Section we study the effects of Lindbladian imperfections on the DDS. We consider the following different imperfections: • Additional set of dephasing dissipative channels: describing fluctuations of on-site energies, which tend to suppress coherences between classical particle number basis of the fermionic system.
• Imperfections in the current Lindblad operators of the form of extra hopping terms: • Coherent Hamiltonian competing with the dissipative dynamics: describing the case in which despite the dissipation being the leading term, there are still some non-negligible coherent local dynamics within the fermionic system. We consider the simplest form of hopping terms, but small local coherent perturbations lead to the same conclusions.
• Additional set of decay dissipative channels: We consider the Lindbladian of the manuscript with parameters A = B = β = t = 1 for a system with L = 9 sites and analyse the effects of the imperfections described in Eqs.(S18),(S19),(S20) on the spectral properties of the Lindbladian. We see that in the presence of an imperfection the DDS is not completely degenerated anymore, with a splitting of the degeneracy proportional to the imperfection strength. We obtain that λ1 = 0 (by definition the Lindbladian has always a zero eigenvalue) while λj=2,3,...,9 < 0 and are approximately equal to each other. We also show in the figure for clarity the eigenvalue λj=10 for = 0, corresponding to the Lindbladian gap in the unperturbed case, which is also approximately equal to the cases with imperfections in the considered range of strengths.
corresponding to losses of particles in the optical lattice due to an extra channel, driving the fermionic system towards an empty vacuum state Our results for the first three imperfections above are shown in Fig.(S7). We see that perturbation theory provides a qualitative picture: in the regime of small perturbations (in units of the rate of the original Lindbladian, 1) while imperfections in the jump operators lead to a linear splitting of the DDS, λ 2 ∼ , a Hamiltonian perturbation leads to a quadratic dependence λ 2 ∼ 2 . Thus, as long as the perturbation is small compared to the unperturbed gap there is a time window between the system entering the DDS and the system characteristics of the state being destroyed by the imperfections, which gives a possibility to effectively use these states for quantum information tasks.
The case of an additional set of decay dissipative channels follows in a similar form. In this case the steady state of the evolution is the vacuum state for decay > 0. However, as above, if the perturbation is small there is a time window over which the effects on the DDS are negligible. One may obtain the characteristic time of the dissipative decay effects by the dynamics of the total number of particlesN in the system, which in the Heisenberg picture is described by L † [N ] = − decayN , i.e., N (t) ∼ e − decay t . Thus the effects of particle losses in the system are relevant for times of the order t ∼ 1/ decay , similarly to the other imperfections in the quantum jump operators considered above.

III. ADIABATIC EVOLUTION
In this Section we expand the discussion of the adiabatic evolution in the dissipative dynamics. Using the same protocol as in the main manuscript, we evolve the system from an initial state given by a superposition of the three charge density wave configurations, precisely, where |Ψ a are the Laughlin states for β = 0 (i.e., product charge density wave configurations). We then evolve this state with the Lindblad operators (Eq.(S15)) using a time-dependent β parameter: β(t) = ∆ · t for 0 < t ≤ 1/∆ and β(t) = 1 for t > 1/∆, where ∆ is the ramp velocity.
We analyse the purity γ(t) of the quantum state during the dynamics, as well as its overlap onto the DDS. For the later we first describe the 3 × 3 density matrixρ DDS representing the expectation value of the quantum state in the Laughlin dark states, (ρ DDS (t)) a,a ≡ Ψ a (t)|ρ(t)|Ψ a (t) for a, a = 0, 1, 2, where |Ψ a (t) are the Laughlin dark states for the Lindblad parameters at time t. We study (i) the projection of the quantum state over the DDS, given by the diagonal terms of the matrix, D DDS (t) = a ρ DDS (t) a,a ; (ii) how the coherence of the initial state evolve under the adiabatic evolution, quantified by and (iii) the distinguishability of the full matrix with the initial state, quantified by the trace norm as follows, Both coherence C(t) as the distinguishability D(t) should be small if the evolved quantum state do not differ significantly from the initial state. We show our results in Fig.(S8). We see that for slow ramp rates the initial quantum state characteristics (purity and coherences) are preserved. We find in particular that for very slow rates ∆ 1 the steady state properties have polynomial corrections compared to their initial conditions, i.e., γ(t → ∞) − γ(0), C(t → ∞) and D(t → ∞) ∼ ∆ cte . Interestingly also to notice that the coherence is approximately constant after an initial time t ∼ O(1/∆) (i.e., β(t) ∼ 1), while it is the diagonal terms of the DDS subspace (D DDS (t)) that still shows non trivial dynamics after this initial transient time.

IV. MAPPING FRACTIONAL QUANTUM HALL GROUNDSTATE TO ONE DIMENSIONAL MODEL.
In this section we revisit the exact mapping of the Laughlin state 27 into a one dimensional state 41,42 . We will be interested in filling fractions ν < 1. Recalling that a 2DEG in a strong magnetic field displays Landau levels, we assume that the relevant physics occurs in the lowest Landau level (LLL). We place the system into a 2D torus with linear sizes L x and L y and area A = L x L y sin θ, defined by the region in the upper half complex plane enclosed by the points w = (0, L x , L y τ, L x + L y τ ). This torus is characterized by the modular parameter τ = L y /L x e iθ = τ 1 + iτ 2 , (Im(τ ) > 0, θ ∈ [0, π]). Following 25 we introduce the translation operators t(L) = exp(L · (∇ − ieA) − iL x y + iL y x) (here L = (L x , L y ) are measured in units of the magnetic length) which correspond to the usual translation operators in terms of the canonical momentum, and an extra space dependent phase. The single particle wavefunction satisfies the boundary conditions t(L a )Ψ = e iφa Ψ, with L a a translation over the lattice vectors L 1 = (L x , 0) and L 2 = L y (cos θ, sin θ). Both conditions can be satisfied if the flux over the torus LxLy sin θ 2π = N Φ is integral. We parameterize the coordinates on the torus by z =z/L x withz = L x (x + yτ ), where x ∈ [0, 1] and y ∈ [0, 1].
The relation with the usual Cartesian coordinates is x 1 = L x (x + τ 1 y) and x 2 = L y τ 2 y The single particle wavefunction has the form Ψ = e − 1 2 (Im(τ )L1y) 2 f (z), where f (z) is an entire (holomorphic) function in the complex plane. The we use units where /eB = 1. In the Landau gauge A = −Byx, the boundary conditions read where the phases φ a correspond to the fluxes piercing the torus in the two orthogonal directions a = 1, 2. From these relations it follows that´dz d dz ln(f (z)) = 2πiN Φ , which implies that the function f (z) has N Φ zeroes. The single particle wavefunction that satisfies the boundary conditions (S25) and has N Φ zeroes is given by the generalized theta function and z n = 2πn+φ2−τ φ1

2πNs
. This corresponds to a normalizable wavefunction for Im(τ ) > 0. The zeroes of ϕ n (z; τ, φ 1 , φ 2 ) are located at z = z n + 1 2 + m + 1 2 + n τ NΦ . As shown in Ref 28, the Laughlin state at filling ν = 1/3 is the zero energy exact ground state of the Landau problem with the interaction H = V 0´d r|∇ρ(r)| 2 , where ρ(r) = ψ † (r)ψ(r) and r = (x 1 , x 2 ). The projection of the electron operator into the first Landau level is ψ = n ϕ n (r)c n , where c n destroys a state at occupation n. The interaction Hamiltonian projected onto the first Landau level becomes In this sum the pair of numbers (j, k) are all integers or all half integers and satisfy 0 < j < N Φ , 0 < k, l, < N Φ /2. Separating both cases, and defining κ 2 = 2π NΦ Lx Ly gives the operators s,n (Eq. 5) in the main text.

V. PHYSICAL REALIZATION
We consider a one dimensional optical lattice (system) immersed in condensate that acts as a bath to the system, providing dissipation. Each site in the optical lattice consists of a potential-well accommodating two single particle levels denoted by c and f Here a i,σ is a fermionic annihilation operator at site i = 1 . . . N Φ and level σ = {0, 1}. The relation with the main text operators is a i,0 = c i and a i,1 = f i . The Hamiltonian includes an attractive (U > 0) interaction between neighbour particles in the level σ = 1. The operator n iσ = a † iσ a iσ measures the occupation at the site i and level σ. We assume that the number of particles in the system is given by N e = N Φ /3.
To capture the essential physics generated by the interaction, we first study the two particle problem. Defining the two particle state |σσ = i,j χ σσ ij a † i,σ a † j,σ |0 , with |0 the state with no particles (vacuum), the Schrödinger equation for the wavefunction χ s ij = 1 j , the discrete Laplace operator. Introducing the central and relative coordinates R = a(j + l)/2 and r = a(j − l) the wavefunction can be written as χ jl = e iRK χ(K) r = e iRK q e irqχs q where we have introduced the total and relative momentum K = k 1 + k 2 and 2q = k 1 − k 2 . The Schrödinger equation for χ s whereχ = q 2 sin qχ s q . For fixed center of mass momentum, the bound state energy E d (K) is found by solving self-consistently Eq. (S29). We consider the regime J 0 ∼ 0, along with |U | E 1 and J1 E1 J1 |U | 1. In this case, the bound state energy E d (k ∼ 0) is far below the bottom of the (1, 0)-pair band, but still above the (0, 0) two-particle band. The amplitude for tunneling between two 0 states is taken to be negligible compared to all other energy scales (J 0 ∼ 0). This implies a flat band for the (0, 0)-pair. The two-particle energies E σσ of the continuous Bloch bands are A doublon state of definite momentum is created by the combination d † The state |d K ≡ d † K |0 is normalized as d K |d K = δ K K , with |0 a state with no particles.

VI. LASER DRIVING AND COUPLING TO A BATH
We are interested in the dynamics generated between the low lying doublon states and the (0, 0) lower energy band. We can induce Raman transitions between states in these bands using an external driving laser with Raman detuning ∆ = 2E 1 − U − ω. The interaction between the (classical) radiation of amplitude Ω and frequency ω, and the system is H rad = Ω cos(ωt) i f † i (c i + α(c i−1 + c i+1 )) + h.c with α = A1 A0 1. The amplitudes A m decay fast with m as they represent the matrix element of the different Wannier functions and the laser radiation profile. We couple our system to a 3-dimensional (3D) Bose Einstein Condensate (BEC). This coupling is realized by immersing the system in the BEC which acts as a dissipative, memoryless bath (valid when the spectral function of the bath is almost constant around the frequency of the laser ω). The bath is very efficient to de-excite the system and does not induce excitations to higher energy bands as its temperature T satisfy T ω so thermal excitations in the bath cannot induce excitations in the system. This type of bath has been used to obtain superconducting states realizing Majorana zero modes 12 . The quasiparticle excitations of this bath are described by the Hamiltonian H bath = k E k b † k b k , where the bogoliubov quasiparticles b k have mass m b and propagate with velocity c b . These bosonic quasiparticles have energy E k = k(c 2 b + k 2 /4m 2 b ) 1/2 ,with k = √ k 2 = (k 2 x + k 2 y + k 2 z ) 1/2 the magnitude of the three dimensional momentum of the Bogoliubov quasiparticles. Here b k (b † k ) destroys (creates) a bosonic bogoliubov excitation of momentum k. We work in the regime where the excitation induced by the laser is relaxed by the bath immediately, such that states with multiple excitations are not created. This approximation corresponds to the limit of weakly far detuned laser E d |∆| |g|, |Ω|. The interaction between the system and the condensate is described by the density interaction H bath/sys = g´drδρ s (r)δρ bath (r) 2 , where g is the strength of the system/bath coupling, which is assumed to be small. The density δρ s (r) is the density of fermions at the point r in 3D, so it involves the Wannier wavefunctions of the fermions. The density δρ bath corresponds to the phonon waves around the equilibrium BEC density induced by the interaction with the fermions in the optical lattice. 2 We move to a rotating frame of reference (the rotating wave approximation, RWA), where the problem, in the reduced Hilbert space which consists of the two lowest bands and the bath, looks static. This is achieved by introducing the time dependent unitary transformation U t = exp (iH sys t) . The penalty is that terms involving higher bands of the extended Hilbert space are fast oscillating. By coarse graining in time we are allowed to get rid of these terms, which amounts to projecting out the higher bands. The radiation Hamiltonian is modified by the rotating wave approximation accordingly. In this rotating basis, the interaction Hamiltonian with the bath H bath/sys is also modified as H int = U t H int U † t . Retaining only the slow oscillating terms (with frequency ∼ ω) in the system-bath interaction, amounts to consider just the transition between the doublon and the (0, 0) band. Each of these transitions is accompanied by the creation (or destruction) of a phonon excitation in the bath.
Tracing over the bath, using the Born-Markov approximation, results in a quantum master equation for the effective density matrix in the reduced Hilbert space involving the bottom and the doublon bands. For small amplitude of the external drive Ω, such that |∆| |g|, |Ω|, no more than one doublon is excited at any given time. We may then trace out the doublon band, and obtain a closed dissipative equation of motion within the Hilbert space of the lowest band. After this adiabatic elimination of the doublon band, we find the quantum master equatioṅ ρ = L(ρ) = NΦ i=1 γ ˜ i ρ˜ † i − 1 2 {˜ † i˜ i , ρ} , with γ = Ω 2 g 2 ∆ 2 Γ 0 and Γ 0 parameterizing details of the bath/system coupling. The quantum jump operator is in turn˜ , and Q i = 1,i + A( 0,i−1 + 0,i+1 ) + B( 1,i+1 + 1,i−1 ).