Out-of-equilibrium steady states of a locally driven lossy qubit array

We find a rich variety of counterintuitive features in the steady states of a qubit array coupled to a dissipative source and sink at two arbitrary sites, using a master equation approach. We show there are setups where increasing the pump and loss rates establishes long-range coherence. At sufficiently strong dissipation, the source or sink effectively generates correlation between its neighboring sites, leading to a striking density-wave order for a class of"resonant"geometries. This effect can be used more widely to engineer nonequilibrium phases. We show the steady states are generically distinct for hard-core bosons and free fermions, and differ significantly from the ones found before in special cases. They are explained by generally applicable ansatzes for the long-time dynamics at weak and strong dissipation. Our findings are relevant for existing photonic setups.


I. INTRODUCTION
Environmental decoherence has long been seen as an unavoidable roadblock to stabilizing quantum phases for long periods of time [1]. However, rapid advances in cooling and trapping techniques over the last decades have led to experimental platforms where the coupling to the environment can be controlled and even engineered to an unprecedented degree [2]. As several studies have shown, such tailored dissipation can be used to prepare novel quantum states [3][4][5]. The competition between Hamiltonian dynamics and incoherent dissipation can produce feature-rich steady states with no analogue in equilibrium condensed matter [6]. Understanding these nonequilibrium phases is of fundamental interest [7], with potential applications in quantum computing [8].
A prototypical experimental setup for exploring such states is a one-dimensional (1D) array of qubits coupled to local reservoirs. In particular, the qubits can be realized by hard-core bosons on a lattice [9], or equivalently, a spin-1/2 XY chain, and the reservoir(s) can be designed to inject or remove a particle (or flip spin) at a given site, as in Ref. [10]. Theoretical studies modeling the resulting dynamics have focused almost exclusively on the cases where the pump and loss occur at the ends of the chain. Then the system can be reduced to free fermions [11], enabling special analytical approaches that have been used to examine nonequilibrium transport [11][12][13][14][15][16][17] and phase transitions [18][19][20]. However, without additional Zeeman fields, the steady state for end drives is rather featureless, with no long-range order [18,21]. On the other hand, a recent work showed that for pump and loss at the center, there are multiple steady states with long-range coherence, that are distinct from free fermions and arise from a dynamical symmetry [22]. Such disparate results beg the question of what happens for generic pump-loss configurations, where pump and loss are neither both at the end nor both at the center. * E-mail: sd843@cam.ac.uk † E-mail: nrc25@cam.ac.uk Here we characterize the steady states for generic setups with a single pump and a single loss site, finding several counterintuitive features which can be probed in already existing platforms [10]. In Sec. III, we use perturbation theory supported by numerics to show there are dipole-like arrangements where long-range coherence is induced by increasing dissipation. In Sec. IV, we show the steady state is generically nonthermal even at weak dissipation, contrary to what is known for symmetric setups [23]. Further, hard-core bosons and free fermions can form qualitatively distinct steady-state correlations, although their density profiles are always reflection symmetric. These attributes are explained by a simple product ansatz of the single-particle modes. In Sec. V, we find that at strong dissipation, the chain is generally divided into a filled and an empty segment separated by a highentropy bulk. These segments are coupled by the source or sink which effectively produces correlation. Whenever two modes in neighboring segments come into resonance, this effect leads to striking density waves and long-range order. This is a geometric effect and can be generalized to multiple sources and sinks. We explain the oscillations by a modified ansatz, finding they are more robust in free fermions than in hard-core bosons.
These results highlight surprising phenomena that can arise in open many-body settings, elucidating differences between hard-core bosons and free fermions in 1D [24]. Our ansatzes apply to more general forms of dissipation, and reduce the numerical cost to linear in system size.

II. MODEL AND KNOWN SPECIAL CASES
We consider strongly interacting bosons on a 1D lattice in the hard-core limit [9], described by the Hamiltonian whereb † i is the boson creation operator, J is the hopping amplitude, and L is the number of sites. The hard-core constraint is encoded in the relationb †2 i = 0, which means no two bosons can occupy the same site. This leads to the commutation rules wheren i :=b † ib i is the local occupation. Such a system is equivalent to a spin-1/2 XX chain, and has been realized with cold atoms in optical lattices [25][26][27] and microwave photons in nonlinear resonators [10]. The Hamiltonian is reduced to free fermions by the Jordan-Wigner map wheref j are the free fermion operators. The transformed Hamiltonian simply readsĤ = − J i f † if i+1 + H.c. . The system of bosons is coupled to bosonic reservoirs that inject particles at a site p, if it is empty, and remove bosons from a site q, if it is occupied. Such sources and sinks have been engineered using transmon qubits in microwave circuits [10] and electron beams in optical traps [28]. We assume the reservoirs are Markovian, i.e., they relax to equilibrium much faster than the system dynamics and the coupling, which is standard for these setups [29]. Upon tracing out the environment, the density operatorρ is governed by the master equation [7,[29][30][31][32] where we have two Lindblad operators modeling the dissipation,L + := √ γ +b † p andL − := √ γ −bq , γ ± being the pump and loss rates, respectively. Note that Eq. (2) does not, in general, reduce the full dynamics to free fermions with (local) pump and loss. Instead, the dissipation mediates nonlocal interactions between the fermions.
The above system has been studied most widely when the source and sink are at opposite ends, i.e., p = 1 and q = L. Then the only term in Eq. (3) that differs from the free-fermion case isL −ρL † − = γ −fLPρPf † L , whereP is the total particle-number parity. SinceP is conserved by the Hamiltonian, one can show the dynamics decouple into sectors withPρP = ±ρ. Thus, the Liouvillian L is quadratic in the free fermions. In such cases, the full solution can be found from the spectral properties of a 4L × 4L matrix using quantization in the space of operators [11]. The steady state is identical to that of end-driven free fermions, characterized by a uniform bulk with short-range correlations [21], as shown in Fig. 1(a). We present a closed-form solution in the Supplement [33].
The end-driven case is in sharp contrast to the scenario where pump and loss both occur at the center site (for odd L), which we explored in a recent work [22]. Here, the system does not map onto free fermions. Instead, one has multiple steady states due to a symmetry operator C = −1/2+ kf † L+1−kf k , which splits the dynamics into (L − 1)/2 sectors with varying degrees of entanglement. The symmetry stabilizes particle-hole pairs at reflectionsymmetric sites k andk := L + 1 − k, leading to steady states with long-range coherence, as shown in Fig. 1(b) for the maximally entangled sector. In comparison, free fermions have an exponentially large set of steady states, as all odd single-particle states vanish at the center.
Outside the above two scenarios, it is known that driving with both pump and loss at a generic site (not center) ib j of hardcore bosons on a 1D lattice subject to incoherent pump with rate γ+ and loss with rate γ− (a) at opposite ends, with γ+ = γ− = 10J, and (b) at the center, with γ− = 4γ+. [34], which is an infinite-temperature state with chemical potential µ = ln(γ + /γ − ). We find free fermions also have the same steady state, except when a single-particle state vanishes at the drive site, producing degeneracies.

III. DISSIPATION INDUCED LONG-RANGE COHERENCE
As described above, long-range order is absent for end drives, and restored for center drives by a special symmetry, irrespective of the pump/loss rate [22]. Here we find examples where long-range coherence is established by increasing dissipation. In particular, consider a "dipole" setup where pump and loss occur at two neighboring sites in the middle, i.e., p = L/2 and q = L/2 + 1, for even L. This can be seen as a center-drive analogue, but there is no strong symmetry [35] and the steady state is unique. Using first-order perturbation theory at weak dissipation (γ ± J), we find the steady state (see Supplement [33]) whereQ := L/2 k=1f † kf k , andρ 0 is the product state with uniform occupation n 0 = γ + /(γ + + γ − ). The perturba-tionQ is reminiscent of the symmetry operatorĈ, and generates the antidiagonal correlations Thus, at weak dissipation, the coherences decay exponentially with distance, and are limited to nearest neighbors (k = L/2) for γ + = γ − , as shown in Fig. 2(a). Conversely, for strong dissipation, the steady state approachesρ step where all sites i ≤ L/2 are filled and i > L/2 are empty.
To first order in J/γ ± , we find (see Supplement [33]), As shown in Fig. 2(b), now the coherences span the entire system, similar to the center-driven case [ Fig. 1(b)]. These results imply that long-range coherence is induced by increasing the pump/loss relative to tunneling, which is confirmed by exact numerics [Figs. 2(c) and 2(d)]. The correlation length grows monotonically with dissipation, exceeding the system size for γ ± 5J. In contrast, for end drives [ Fig. 1(a)], coherences are limited to nearest neighbors at both weak and strong dissipation, yielding (see Supplement [33] Numerically, we find similar results in "dipole" setups with q = p + 1, whenever p divides L − p or vice versa (see Supplement [33] for examples). In Sec. V, we discuss more general scenarios where such geometric resonances stabilize long-range order in the Zeno limit.

IV. NONTHERMAL STEADY STATES AT WEAK DISSIPATION
For end drives as well as for the "dipole" setup considered in Sec. III, the steady state approaches the uniform infinite-temperature stateρ 0 in the weak-coupling limit. Here, the system has time to equilibrate between successive pump and loss events. Thus, one might expect the same steady state regardless of where those events occur. This is indeed true whenever the pump and loss act on reflection-symmetric sites [23] or at the same site (except center) [34]. However, we find those are the only setups where the conjecture holds. As shown in Figs. 3(a)-(b), the steady state is generically nonthermal, and different for free fermions and hard-core bosons, although they both have symmetric densities. These features can be understood by focusing instead on the single-particle modes which are unaffected by the Hamiltonian. They are given by unitary mapsF m = j c m,jfj where and have energy ε m = −2 J cos πm L+1 , m = 1, . . . , L. If the dissipation is small compared to the energy splitting, the modes become uncorrelated. Thus, we find the steady state is well approximated by a "product-of-modes" form ρ ≈ m (1−N m )|0 m 0 m |+N m |1 m 1 m |, set by the mode occupations N m . The lack of correlation explains why the density is symmetric. Using F † mFn ≈ δ mn N m yields n j := n j = m N m |c m,j | 2 , which gives nj = n j for any site j. (Recall thatj := L + 1 − j.) One can use the product-of-modes ansatz in Eq. (3) to construct approximate rate equations for N m . For free fermions, one finds (see Supplement [33]) i.e., the modes are uncoupled, with incoming and outgoing currents set by the weight of the respective mode at the pump and loss sites. To reach a uniform steady state such asρ 0 , one must have |c m,p | = |c m,q | ∀m, which is satisfied iff q = p or q =p. While this conclusion also holds for hard-core bosons, the rate equations are more complex as the string operator in Eq. (2) generates nonlinear coupling between the modes. Using F † mF † m N n N n (δ m,n δ m ,n −δ m,n δ m ,n ), we find (see derivation in the Supplement [33]) In numerical trials, Eq. (11) gives a unique steady state with 0 ≤ N m ≤ 1 for all m. Figure 3(a) shows the ansatz accurately describes the steady states in the weak-coupling limit, becoming exact for free fermions. The densities are peaked at the pump site and minimized at the loss site, with n q = 1 − n p for γ + = γ − . The same is mirrored at sitesp andq, which explains why choosing q =p [23] or q = p [34] gives half filling at all sites. In general, the density fluctuations are significantly smaller for hard-core bosons due to strong interactions. There are also qualitative differences which persist to large system sizes, as shown in Figs Here one has pump at one end and loss at the center. For free fermions, all the odd modes (even m) are immune to loss and thus fully filled, which is not the case for hardcore bosons [ Fig. 3(c)]. Thus, we find strikingly different densities and correlations. In particular, Fig. 3(d) shows that unlike bosons, fermions exhibit long-range order in the density correlations , which holds more generally at weak dissipation. Also, note the free fermions have degenerate steady states if any of the modes vanishes at both pump and loss sites, whereas for hard-core bosons the steady state is unique except for center drive [22]. Equations (10) and (11) apply for any quadratic Hamiltonian with nondegenerate spectrum, reducing the dynamics to L rate equations.

V. GEOMETRIC RESONANCE AND LONG-RANGE ORDER IN ZENO LIMIT
In Sec. III we found, for pump and loss at neighboring sites, q = p + 1, the steady state at strong dissipation approaches a step where sites 1 through p are filled and sites q through L are empty [Eq. (6)]. This can be understood as follows. For t 1/γ ± , sites p and q are pinned at occupation 1 and 0, respectively. By tracing over this subspace, one can show the remaining sites are governed by a master equation with an effective HamiltonianĤ eff and weak effective dissipation, as detailed in Ref. [36] for more general systems. In our setup,Ĥ eff simply describes hopping in the disjoint segments 1 to p − 1, p + 1 to q − 1, and q + 1 to L. For q = p + 1, the middle region is absent and the dissipation is given by Lindblad operatorŝ [33]). The former injects particles into the first segment until it is filled, and the latter removes all particles from the last segment. More generally, i.e., the source and sink induce correlated pump and loss at neighboring sites with rate Γ ± , coupling the segments. However, this is a second-order effect (Γ ± /J ∼ J 2 /γ 2 ± ), and the steady state remains uncorrelated whenever the energy splitting between modes in adjacent segments is large compared to Γ ± , as we explain below. Then the mid region behaves like an end-driven qubit array, reaching a uniform product state with density n (2) = Γ + /(Γ + + Γ − ). This is similar toρ 0 in Eq. (4) except n (2) decreases with γ + /γ − , a consequence of the Zeno effect [37]. Thus, the steady state generically consists of fully filled and empty sites separated by a region of high entropy.
This picture breaks down if any two modes in neighboring segments are resonant. Then they can remain coherent via the correlated pump and loss, producing characteristic density waves as in Fig. 4. To understand this feature, note the single-particle modes in segment ν are given byF mν ,jν as in Eq. (9), for m ν = 1, . . . , L ν , where L ν is the number of sites. They have energies ε mν := πm ν /(L ν + 1). Thus, off-resonant modes dephase at a rate ∆ε/ ∼ J, much faster than the dissipative coupling Γ ± . As a result, for t 1/Γ ± , adiabatic elimination gives is 1 for resonant modes and 0 otherwise. Since the energies are set by the wavenumber k (ν) mν , the resonance condition is equivalent to a single plane wave fitting into both segments ν and ν + 1, as in Fig. 4(a). The corresponding modes are seen to have nonzero steady-state correlation in Fig. 4(b). Such modes exist iff L ν + 1 divides L ν+1 + 1, or vice versa, making these arrangements special. Within each segment, one finds F (ν) † These attributes can be explained by approximate rate equations for the modes. For free fermions, one finds (see  derivation in Supplement [33]) where u (ν) mν ,Lν are the mode amplitudes at the boundary where pump or loss occurs, and m3 | 2 is an effective decay rate. The equations for the remaining segments follow by symmetry [33]. Note the occupation N For hard-core bosons, Eq. (13) gains an interaction terṁ T − , which can be approximated by pairwise contractions, as in Eq. (11). The main result is an increased decay rate, [33]), which weakens the correlations, making the resonant features less prominent.
To further elucidate the characteristic features of hardcore bosons and free fermions, we consider a set of pumploss configurations with p = 1, q = l+1, and L = (r +1)l, for positive integers l and r. Here the sink is positioned such that every r-th mode in the last segment is resonant with successive modes in the middle [ Fig. 5(a)]. In steady state, this leads to a train of density bumps in the former, with nodes at every l-th site [ Fig. 5(b)]. As expected from the preceding paragraph, the bumps are smaller for hard-core bosons, well reproduced by the resonant-modes ansatz. In addition, free fermions show recurrent longrange density correlations of period 2l [ Fig. 5(c)], similar to the center-driven case [22]. For hard-core bosons, these are less prominent, and cloaked in a uniform background due to interactions [ Fig. 5(d)]. The reduced density fluctuations for bosons originate from a faster decay of correlations. Using c m3 ∼ O(1/l) for hard-core bosons. Thus, we find a density bump n q+1 ≈ m3 |u m3 ∼ 1/L in the latter, whereas in the former, n q+1 ∼ O(1/r), regardless of l. This different scaling with system size is evident in Figs. 5(e) and 5(f). Since the physics is solely determined by geometry, the framework can readily incorporate multiple sources and sinks. Further, Eq. (13) and its bosonic counterpart remain valid in the presence of a trap as long as the spectrum in each segment is nondegenerate.

VI. SUMMARY AND OUTLOOK
We have characterized a rich class of steady states that arise in a prototypical setting of two-level systems driven by localized pump and loss. For different arrangements, we find dissipation induced long-range coherence (Fig. 2) and surprising geometric order at strong pump and loss (Fig. 5). These results could be observed by measuring local densities and density correlations in state-of-the-art photonic setups [38]. We have developed a general framework to approximate the long-time dynamics at weak and strong dissipation in terms of rate equations for the energy eigenmodes, which would be useful in other systems. Note the strongly dissipative limit is reached in practice for γ ± 10J. Quite generally, we find a local source or sink generates correlation between the neighboring sites [Eq. (12)], which can stabilize long-range order [22]. This effect is present whenever the lossy site equilibrates much faster (γ ± ) than the global relaxation rate (J/L) [36], and provides strong motivation for using such dissipation to engineer correlated states of matter [39]. Thus, it would be valuable to extend our study to interacting Hamiltonians, such as an XXZ chain [40] or a Hubbard model [7], and to higher dimensions [41] where there is no simple map between hard-core bosons and free fermions [42]. It would also be interesting to see how the effect is altered in the continuum. Note the local pump and loss drives current through an interacting medium. Thus, future work could harness such probes to extract useful information about the bulk [43] and investigate larger questions of nonequilibrium transport [44]. Supplement for "Out-of-equilibrium steady states of a locally driven lossy qubit array" As described in the main text, we consider hard-core bosons on a 1D lattice modeled by the Hamiltonian where the boson operators satisfy 1}. The bosons are coupled to Markovian reservoirs that inject particles at site p with rate γ + and removes particles from site q with rate γ − . The resulting dynamics are modeled by a master equation for the density matrixρ whereL + := √ γ +b † p andL − := √ γ −bq are Lindblad operators describing the incoherent pump and loss, respectively.
The system is equivalent to a spin-1/2 XX chain with local spin flips if one identifiesb i with the spin lowering operator. It can also be mapped onto fermions by a Jordan-Wigner transformation, where {f i ,f j } = 0 and {f i ,f † j } = δ ij . The transformed Hamiltonian describes free fermions, However, the Lindblad operatorsL ± are nonlocal in the fermions, mediating interactions.

A. CLOSED-FORM SOLUTION FOR END DRIVES
As shown in Ref. [S1] and discussed in the main text, for pump and loss at opposite ends (p = 1, q = L), the dynamics are identical to those of free fermions, i.e., withL + := √ γ +f † 1 andL − := √ γ −fL in Eq. (S2). The steady state is unique [S2] and has been solved exactly in terms of a matrix product ansatz [S3]. Here we present a closed-form solution by a more direct approach. We adopt a thermofield representation [S4, S5] where one defines a new set of operatorsf i andf † i that act on the density matrix by right multiplication, i.e.,f i ρ := ρf i andf † i ρ := ρf † i . Note we have omitted the hat for operators to reduce clutter. It is straightforward to verify the following relations for any two operators a and b: The Liouvillian L in Eq. (S2) can be expressed in this notation as where and To write down the steady state, we define the generators and A := We will show the steady state is given by ρ = e Gτ ρ 0 , where and ρ 0 is a uniform product state with occupation n 0 = γ + /(γ + + γ − ), To prove that Lρ = 0, we first use the identities in Eq. (S5) to find the commutators Using these results in Eqs. (S6) and (S9) yields The expression within parentheses can be evaluated by noting that ρ 0 = (γ + /γ − ) N up to normalization, where N is the total number operator. Hence, [H, ρ 0 ] = 0, or T ρ 0 = 0. One also finds, using Eqs. (S11) and (S13), Substituting the above results in Eq. (S14) gives Lρ = 0, thus showing ρ is indeed the steady state. The computation of ρ can be simplified by noting the generators in Eq. (S9) commute with one another, Thus, e Gτ factorizes into a product of exponentials. Further, W i,i acts locally and one can show W 2 i,i = 0 ∀i, which means the action of these local generators on ρ 0 can be written explicitly, yielding Furthermore, A L+1 = 0, so the above exponential reduces to a sum of the first L + 1 terms in its power series, which can be computed iteratively using A σ = . Such a simple algebraic structure of the steady state is related to the q-deformed SU(2) symmetry of the XXZ chain [S6].
For weak dissipation (γ ± J), Jτ . Thus, to linear order in γ ± /J, Eq. (S16) gives the steady state where we have used Eq. (S11) to simplify Aρ 0 . The perturbation to ρ 0 induces nearest-neighbor correlations, . (S18) , and the product state in Eq. (S16) approaches ρ Z where the first site is filled, the last site is empty, and the other sites are in a product state with occupation 1 − n 0 , To first order in J/γ ± , the steady state is given by which again yields nearest-neighbor correlations, Hence, the correlations are limited to nearest neighbors in both limits, and are purely imaginary, which corresponds to a probability current from the source at the first site to the sink at the last site.

B. PERTURBATIVE SOLUTIONS FOR A DIPOLE DRIVE
In Sec. III of the main text, we described a "dipole" arrangement of the pump and loss, where increasing dissipation establishes long-range coherence, in sharp contrast to the end-driven case studied above. In this "dipole" setup, the pump and loss occur at neighboring sites in the middle, p = L/2 and q = L/2 + 1 for even L. Here we derive the perturbation results for weak and strong dissipation quoted in the main text.
The Liouvillian L in Eq. (S2) can be restated as where D[x]ρ :=xρx † − {x †x ,ρ}/2. The steady state at weak dissipation approaches the product stateρ 0 ∝ (γ + /γ − )N , as in the end-driven geometry [see Eq. (S11)]. This is the zeroth order solution, which commutes with the Hamiltonian. The solution to first order in γ ± /J is given byρ w ≈ρ 0 +ρ 1 , such that We will show this is satisfied byρ First, we calculate the action of the dissipators on the unperturbed solution. As in Eqs. (S15a) and (S15b), we find Next, we find the commutator [Ĥ,ρ 1 ] using Eq. (S4) and the identity [ab, cd] = a{b, c}d + ca{b, d} − {a, c}bd − c{a, d}b, Combining Eqs. (S26) and (S27) readily gives the first-order condition in Eq. (S23). Note the perturbationρ 1 induces coherence between reflection-symmetric sites k and L + 1 − k, leading to the single-particle correlations (for k ≤ L/2) (S28) Thus, the correlations fall off exponentially with distance due to the string, and vanish for k < L/2 if γ + = γ − . At strong dissipation (γ ± J), the steady state approaches a step where all sites i ≤ L/2 are filled and all sites i > L/2 are empty (see Sec. V in the main text for a physical explanation). This pure state is expressed aŝ and annihilated by the dissipators D[b † L/2 ] and D[b L/2+1 ] in Eq. (S22). To first order in J/γ ± , the steady state is of the formρ s ≈ρ step +ρ 1 , such that We will show the perturbationρ 1 is again generated by the operatorQ in Eq. (S25), First, using the expression forĤ in Eq. (S4), one finds Next, acting the dissipator D[b † L/2 ] onρ 1 and usingρ stepbL/2 = 0 yields The same result is found for D[b L/2+1 ]ρ 1 . Thus, This exactly cancels −(i/ )[Ĥ,ρ step ] in Eq. (S32), satisfying the first-order condition in Eq. (S30). This single-particle correlations can be calculated by a similar procedure as in Eq. (S28), yielding (for k ≤ L/2) Now the string gives rise to constant-amplitude oscillations, stabilizing long-range coherence. Note the steady states in Eqs. (S24) and (S31) can be written as a compact matrix product operator, as in Refs. [S3, S6].

C. GEOMETRIC RESONANCE AND LONG-RANGE COHERENCE IN DIPOLE SETUPS
In the main text, we mentioned that the dissipation induced long-range coherence found above can be generalized to other "dipole" setups where the pump and loss act on neighboring sites (q = p + 1), provided p divides L − p, or vice versa. Here we present numerical examples. Figure S1 shows the single-particle density matrices in steady state for different dipole arrangements at strong dissipation. Here the pump and loss divide the system into weakly coupled segments, sites 1 through p and sites q through L, as explained in Sec. V of the main text. Long-range coherence is found when the two segments share a resonant mode, as sketched in the upper panels of Fig. S1. This is a higher-order analog of the geometric resonances discussed in the main text. In all of the examples, the coherences are limited to nearest neighbors at weak dissipation (γ ± J). Hence, the resonant geometries constitute a family of setups where long-range coherence is stabilized by increasing the pump and loss rates.
FIG. S1. Top panel shows schematic "dipole" setups where the pump and loss divide the system into two segments which may share a resonant mode, shown by black lines. Bottom panel shows steady-state correlations b † ib j for the corresponding setups at strong dissipation, γ+ = γ− = 20J. The numbers in parentheses give the size of the segments, p and L − p. The first four plots are obtained from exact diagonalization and the last plot is obtained by averaging over quantum trajectories [S7].

D. RATE EQUATIONS FOR WEAK DISSIPATION
In Sec. IV of the main text, we explained that when the dissipation is weak compared to the energy splitting of the single-particle modes, the modes become uncorrelated, leading to simplified rate equations for the mode occupations which characterize the steady state. Here we derive these approximate rate equations.
First, we consider the equation of motion for the expectation of a general observableÂ. Substituting Â = Tr(Âρ) into Eq. (S2) and using the cyclic property of trace, one finds The single-particle modes of the Hamiltonian in Eq. (S4) are of the formF m = j c m,jfj , such that m c * m,i c m,j = δ ij . One can invert this relation to findf The number operator for a mode is given byN m :=F † mFm which, by definition, commutes withĤ. The equation of motion for the mode occupation N m := N m is simpler if one has pump and loss of free fermions, i.e.,L + = √ γ +f † p The commutators can be found using Eq. (S37) and the relations {F m ,F n } = 0 and {F † m ,F n } = δ m,n , which gives Now we approximate the modes to be uncorrelated, i.e., F † nFm ≈ δ m,n N m , obtaining f † q [f q ,N m ] ≈ |c m,q | 2 N m . The pump term in Eq. (S38) is found by exchanging p ↔ q and particles with holes, yieldinġ Mapping the bosons onto fermions through Eq. (S3) and using the expansion in Eq. (S37), we find whereN q := i<qn i . The interaction can be simplified by writingN m = i,j c * m,i c m,jf † if j and noting that (−1)N q transformsf i to −f i only if i < q. Thus, where σ k = 1 for k ≥ 0 and −1 for k < 0. Rewriting thef j 's in terms of the modes in Eq. (S37), one finds where Thus, the string gives rise to quartic coupling among the modes. Note the above result is exact for hard-core bosons. We approximate the quartic terms using the product-of-modes ansatz,ρ ≈ N m |1 m 1 m | +N m |0 m 0 m |, presented in the main text (recall,x := 1 − x), finding of their individual eigenvalues and eigenvectors listed in Table I. The eigenvalues of L 0 are given by ξ ν,ν = γ + ξ p ν +γ − ξ q ν , and the corresponding eigenvectors areψ ν,ν =ψ p ν ⊗ψ q ν andφ ν,ν =φ p ν ⊗φ q ν . One can readily verify the eigenvectors are normalized such that Tr Hp (φ p µψ p ν ) = δ µ,ν and Tr Hq (φ q µψ q ν ) = δ µ,ν . As expected, the target stateψ 0,0 describes a filled pump site and an empty loss site in the subspace H 0 = H p ⊗ H q . Using Eq. (S56), we find Thus, the Hamiltonian projected onto the Zeno subspace simply describes hopping in three uncoupled segments. The only nonzero coefficients Y (µ,µ ),(ν,ν ) in Eq. (S57), with ξ µ,µ , ξ ν,ν = 0, are The corresponding dissipatorsĝ ν,ν in Eq. (S56) are obtained as (for q > p) with rates Γ ± := 4J 2 /γ ± . Therefore, the source and sink generates correlation between its neighboring sites through a second-order process, dissipatively coupling the segments in Eq. (S58). For pump and loss of free fermions instead of hard-core bosons, one finds the same expressions withb's replaced byf 's in Eqs. (S61).

F. RATE EQUATIONS FOR STRONG DISSIPATION
We showed above that the qubit array reduces to weakly coupled segments at strong dissipation. Since this coupling is small compared to tunneling, any off-resonant modes become uncorrelated at long times, to a good approximation, as discussed in Sec. V of the main article. This is similar to what happens at weak dissipation (see Sec. D). However, the modes that are resonant in neighboring segments can remain coherent, producing surprising steady states. Here we derive the approximate rate equations that explain these features.
The pump and loss divide the system into three segments as in Eq. (S58). The single-particle modes in segment ν, with L ν sites, have the formF Note, in this section we use ν to label the segments, and not eigenvalues as in the last section! The modes in neighboring segments are coupled through the dissipators in Eqs. (S61). For simplicity, we assume all three segments are present, i.e., 1 < p < q −1 < L−1, which giveŝ To find the equations of motion, it is useful to write the dissipators in terms of the modes. To this end, we first invert the unitary coefficients c (ν) mν ,jν to findf (ν) jν , then attach appropriate string operators according to Eq. (S3), yieldinĝ where u (ν) mν ,Lν , N 1 is the total occupation in the first segment, and N 1,2 is the total occupation in the first two segments. The extra minus signs arise from the filled pump site and do not affect the physics. Note the convention for N 1 is different from that used in Sec. D. As explained in the main text, the steady state is set by the occupations N m3 from first principles, and later work out those for the other segments by symmetry. The equation of motion for the occupations can be found from Eq. (S36), Using Eqs. (S64), one finds [L eff + ,N m3 ] = 0 and [b q+1 ,N . (S66) Next, we make the approximation that all off-resonant modes are uncorrelated and that the spectrum is nondegenerate in a given segment, i.e., F (ν) † The same rate equation is obtained if the hard-core bosons were replaced by free fermions in Eq. (S63). To find the rate equation for the correlations, we first generalize Eq. (S36) for a non-Hermitian operatorX. Using X = Tr(Xρ) in Eq. (S2), along with the effective Hamiltonians and dissipators, one finds We consider the correlation between two resonant modes m 2 and m 3 , i.e.,X =F (2) † m2F m3 + 2X f q−1 +f q+1 , Combining these results with Eqs. (S64) yields (2) † m2F where Γ f := Γ+ 2 |u m3 | 2 , and ζ b = 1 for the hard-core bosons we have been considering. For free fermion dissipators, one obtains the same equation without the quartic term, i.e., ζ b = 0.