A possible route towards dissipation-protected qubits using a multidimensional dark space and its symmetries

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. Here, 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.

I t 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 nontrivial topological states [12][13][14] . This is achieved by a careful interplay between radiation-induced 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 nontrivial 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 two dimensions 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.

Results
Projective representation of symmetries in Hamiltonian dynamics. As a warm up for the symmetry discussion, let us consider the Hamiltonian case. If a Hamiltonian is invariant under the action of a symmetry group G, then action of the group elements, g 2 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 ψ j i and e iϕ ψ j 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 where D(g) is a representation of a group element g 2 G. Every projective representation is characterized by the set of phases ω 2 ðg 1 ; g 2 Þ ¼ e iϕ ðg 1 ;g 2 Þ , 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 2cocycle (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 nontrivial 2-cocycle, representations of at least some group elements do not commute, even for an Abelian group G, i.e.
for some g; h 2 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 DðgÞHD y ðgÞ ¼ H; DðgÞ' i D y ðgÞ ¼ (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 2 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 j i 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 j i ¼ e iðαþφðh;gÞÞ DðgÞ e j i. Forφðh; gÞ ¼ 2πa m , with a, m coprime 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 eigenmatrices of L need not to be selfadjoint, but can come in pairs {σ, σ † } with eigenvalues fλ σ ; λ Ã σ g. 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. (3), as' i ¼ R y i ' i , where R y 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 ¼ P n ð' y 0;n ' 0;n þ ' y 1;n ' 1;n Þ, where the operators ℓ s,n (s = 0, 1) are 29 Here c n destroys an electron at orbital n; ηðlÞ / e Àðκl Þ 2 is a fast decaying function in the narrow torus limit, . The crucial property of Laughlin states Ψ j i that makes them useful for our discussion of the Lindbladian evolution is that they satisfies ' s;n Ψ j i ¼ 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, which 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 TUT † = e −2πiν U. They thus provide a projective representation of One can choose a basis Ψ a j i, 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 In the context of the quantum Hall effect, states Ψ a j i 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 j i Ψ a h j for a = 0…m − 1 share the same eigenvalue (they are related by conjugation with T). In general, the matrices Ψ a j ihΨ aþp j 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 frustrationfree condition which ensures that if Ψ a j i is annihilated by the quantum jump operators, then all the matrices Ψ a j i Ψ b h j 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 (switching from the orbital guiding center index n to the real space site index i) s 2 Þ ' 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 where f i destroys a fermion in this at position i in this chain and n i ¼ f y i f i . 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 coupling Hamiltonian 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, which provides dissipation. These components are shown in Fig. 2a As a motivation for this construction, let us 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 ( Ω 2 2E 1 ÀUÀω ), 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 y i c y iþ1 0 j i will be affected by radiation, as well as c y i c y iþ2 0 j i, where 0 j i is the state with no particles. The first local configuration that is not affected by radiation, becoming dark is c y i c y iþ3 0 j i, 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 Supplementary Note 4). 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 j i ¼ f y i f y iþ1 0 j i and the lower band states i; j j i ¼ c y i c y j 0 j i. The transition process from lower band to doublon reads (after adiabatic elimination of the doublon) (8). The prefactor A Ω~Ω 2 /(2E 1 − U − ω), 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 y 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 y 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 Ψ a j i, where a = 0, 1, 2, which are annihilated by all composite operators, ' s;i Ψ a j i ¼ 0, for all s, i (and thus by the quantum jump operators). Assuming periodic boundary conditions, these states are given by 29 where N is a normalization factor, a = 0, 1, 2 and The state j i 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 j i, are related by a translation T by one site, as T Ψ a j i ¼ Ψ aÈ1 , where ⊕ is an addition modulo 3. In the basis Ψ a j i, the operators T and U are represented by the 3 × 3 matrices Within DDS the density matrix is ρ ¼ P 2 a;b¼0 ϱ ab Ψ a j i Ψ b h j, 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. j i , 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 . 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 j i Ψ a h jÞ. 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 2ln ð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 j i , 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 34 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Þ ¼ Trfρ 2 ðtÞg. 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 f Ψ a j ig a¼0;1;2 , may be visualized by where P DDS ¼ P 2 a¼0 Ψ a j i Ψ a h j 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 At β = 0 the system is in the CDW state, with a definite spatial periodicity, indicated by the peaks in S 1 (k) at k ¼ 0; 2π 3 ; 4π 3 . Increasing |β|, the system becomes more homogeneous. c Entanglement entropy of the DDS as a function of β. 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 behavior. 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 [35][36][37] . 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 38 . 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.

Discussion
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.

Methods
Mapping fractional quantum Hall ground state to one-dimensional model. In this section we revisit the exact mapping of the Laughlin state 27 into a one-dimensional state 39,40 . 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 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; θ 2 ½0; πÞ. Following ref. 25 we introduce the translation operators tðLÞ ¼ expðL Á ð∇ À ie AÞ À 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 spacedependent 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 L x L y 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 À 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 R dz d d z ln ðf ðzÞÞ ¼ 2πi N Φ , which implies that the function f(z) has N Φ zeroes. The single-particle wavefunction that satisfies boundary conditions (15) and has N Φ zeroes is given by the generalized theta function ϕ n ðz; τ; ϕ 1 ; ϕ 2 Þ ¼ e À 1 2 ImðτÞL x y ð and z n ¼ 2πn þϕ 2 Àτϕ 1 2πN s . 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 R drj∇ρðrÞj 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π gives the operators ℓ s,n (Eq. 5) in the main text.
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 À P N i;σ ðJ σ a y i;σ a iþ1;σ þ h:c:Þ and 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 neighbor particles in the level σ = 1. The operator n iσ ¼ a y 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 σσ 0 j i¼ P i;j χ σσ 0 ij a y i;σ a y j;σ 0 0 j i, with 0 j i the state with no particles (vacuum), the Schrödinger equation for the wavefunction χ s ij ¼ 1 ffiffi 2 p ðχ 11 ij À χ 11 ji Þ reads ÀJ 1 ðΔ i þ Δ j Þχ s ij À Uðδ i;jþ1 þ δ iþ1;j Þχ s ij ¼ ðE À 2E 1 þ 4J 1 Þχ s ij , with Δ i χ ij = χ i+1,j − 2χ ij + χ i−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 P 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 ij becomes where χ ¼ P 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. (19), leading to E d ðKÞ ¼ 2E 1 À U À 4 J 2 1 U cos 2 K 2 . We consider the regime J 0~0 , along with |U| ≫ E 1 and J 1 E 1 ( J 1 jUj ( 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 twoparticle energies E σσ 0 of the continuous Bloch bands are E σσ 0 ðK; qÞ ¼ ðσ þ σ 0 ÞE 1 À 2ðJ σ þ J σ 0 Þ cos Ka 2 cos qa þ 2ðJ σ À J σ 0 Þ sin Ka 2 sin qa;

Data availability
The data used to create the plots are available from the authors upon reasonable request. The figures were produced using Python and processed using Inkscape.

Code availability
The numerical codes and scripts to obtain the data are available from the authors upon reasonable request.