Exact Results for a Boundary-Driven Double Spin Chain and Resource-Efficient Remote Entanglement Stabilization

We derive an exact solution for the steady state of a setup where two $XX$-coupled $N$-qubit spin chains (with possibly non-uniform couplings) are subject to boundary Rabi drives, and common boundary loss generated by a waveguide (either bidirectional or unidirectional). For a wide range of parameters, this system has a pure entangled steady state, providing a means for stabilizing remote multi-qubit entanglement without the use of squeezed light. Our solution also provides insights into a single boundary-driven dissipative $XX$ spin chain that maps to an interacting fermionic model. The non-equilibrium steady state exhibits surprising correlation effects, including an emergent pairing of hole excitations that arises from dynamically constrained hopping. Our system could be implemented in a number of experimental platforms, including circuit QED.

In this work, we present new exact results for two boundary-driven spin models that are directly relevant to both of the above motivations.The first (Fig. 1(a)) consists of two passively-coupled N -qubit chains that hang off the same waveguide.We show that for arbitrary N , this system has a pure, highly-entangled steady, even for weak driving and with certain kinds of disorder.The second (Fig. 1(b)) is a single qubit chain with boundary Rabi driving and loss, which somewhat surprisingly corresponds to an interacting fermionic model.We nonetheless obtain an exact result for the NESS by considering the directional waveguide limit of the double-chain system: the double-chain pure state is the purification of the desired NESS.This represents the first use of the hidden time-reversal / quantum absorber exact solution method [4,36,37] to a non-trivial system where interactions are not long-ranged.This solution is also a rare example of an exactly solved coherently-driven 1D spin chain model.In contrast, the few existing examples of exactly solvable boundary-driven spin chains have involved purely incoherent drives [23][24][25][26][27]. * lingenfelter@uchicago.eduFIG. 1.(a) Two XX-coupled N qubit chains, with boundary collective loss (rate γ, mediated by a waveguide) and Rabi drives (strength Ω, detuning ±∆).(b) A single XX chain with a Rabi drive and loss on the boundary, which corresponds to an interacting fermionic model.The steady state here is obtained from the pure steady state of the two-chain system.(c) For strong dimensionless driving Ω (cf.Eq. ( 26)), the steady state approaches a product of dimer Bell pairs.For finite driving, the exact steady state is obtained by doping this state with delocalized, paired "holes".Strong inter-chain entanglement can be achieved even if Ω < γ.
Our exact solutions provide a wealth of insights relevant to understanding correlations in the NESS, and to the application of remote entanglement stabilization.Despite the lack of any explicit attractive interactions in our systems, their steady states exhibit strong realspace pairing correlations.In fact, we show that the pure steady state of the double-chain system can be exactly written as a condensate of paired holes, where a hole here corresponds to an interchain dimer of qubits that are both in the vacuum state (see Fig. 1(c)).We discuss how this pairing has clear observable consequences, and how it ultimately arises from a kinetically constrained hopping process (something that could also be studied in certain non-dissipative cold atom systems [38]).This pairing mechanism also has a direct connection to the emergence of quantum scar states: it represents a restricted spectrum generating algebra of the two-chain Hamiltonian [39], and can be used to construct a tower of many-body scar states in related nonintegrable ladder models such as that studied in Ref. [40].We also discuss regimes where this pairing structure directly results in an NESS with features reminiscent of charge density wave order.
In terms of entanglement stabilization, our doublechain system has potential advantages over other proposals [16][17][18][19][20][21] as it does not require the preparation and transport of high-fidelity squeezed light, but only simple Rabi drives and passive waveguide couplings.It unifies and extends to the multi-qubit regime previously known two-qubit schemes [4][5][6][7][8].For large Rabi drive amplitude Ω, our system is just one example of a more general mechanism for replicating definite-parity two qubit states using simple XX (hopping) couplings.While this was seen previously in different systems [16][17][18][19], the underlying mechanism had not been elucidated.For finite Ω, the two qubit systems need Ω > γ for significant steady state entanglement, where γ is the waveguide-induced dissipation.We find that this is no longer the case for larger systems: for N ≥ 2, strong steady state entanglement only requires the potentially much weaker condition Ω > √ Jγ, where J is the hopping (XX coupling) in each chain.We also show that the method introduced in Ref. [8] for speeding up the dissipative stabilization of a two-qubit entangled state can be extended to situations with multiple qubits, leading to a dramatic acceleration of our protocol.The rest of this paper is organized as follows.In Sec.II we introduce our two basic models, while in Sec.III we give the exact pure steady state of the double chain, discuss its structure, and explain the general entanglement replication mechanism that applies for large Ω.In Sec.IV we explain three surprising features of the steady state: effective hole-pairing, universal single-parameter scaling of the excitation density, and the emergence of single particle states with charge density wave order.In Sec.V we discuss the application of Fig. 1(a) to remote many-Bellpair entanglement stabilization.

II. XX-COUPLED QUBIT CHAINS WITH BOUNDARY DISSIPATION AND DRIVING
A. The two-chain model Consider the setup in Fig. 1(a): two passively coupled N -qubit chains (A and B), with boundary driving and correlated dissipation.The dynamics is described by the Lindblad master equation where the Hamiltonian terms are given by The Lindblad dissipator • }/2 describes collective loss on the site 1 qubits (with σ− = |0⟩⟨1|, σz = |1⟩⟨1| − |0⟩⟨0|).
Consider first the driving of our system.Qubits A1 and B1 are Rabi-driven at the same frequency and amplitude Ω.We however take qubit A1 (B1) to be detuned by +∆ (−∆) from the drive frequency.Treating these drives within the rotating wave approximation, we obtain the rotating frame Hamiltonian Ĥdrive , given by Eq. ( 4).We will take the remaining qubits in each chain to be resonant with the drive frequency.
Within each chain, excitations can hop between adjacent qubits.This is described by simple nearest-neighbor XX couplings: ĤXX , given by Eq. (5).While the hopping amplitudes J j in each chain can vary from bond to bond, we require that the hopping across a particular bond j is the same for chain A and B; as we will see in Sec.III, this mirror symmetry is necessary to obtain a pure steady state.The passive exchange couplings we use here are natural in many experimental settings.For example, in superconducting circuits they could be realized straightforwardly with capacitive couplings.
Finally, we turn to the dissipation-mediated coupling between the two chains.Qubits A1 and B1 experience collective loss (at rate γ) due to a common coupling to a Markovian reservoir.We focus on the case where this bath is an open waveguide structure, and the two chains are spatially separated (in the limit where non-Markovian effects associated with a finite propagation time can be neglected).We consider two types of waveguides: (i) a bidirectional waveguide that supports both left-and right-propagating waves, or (ii) a unidirectional waveguide that supports only, e.g., right-propagating waves.Note that a bidirectional waveguide requires precise spacing of the qubits to engineer collective loss (see e.g.Ref. [7]); such control is not needed for the directional setup.When the waveguide is not fully bidirectional (e.g.qubits couple preferentially to right-propagating modes versus left-propagating modes), it induces the effective exchange Hamiltonian Ĥdiss given by Eq. ( 6) [41], where ν is a directionality factor, −1 ≤ ν ≤ 1.When ν = 0 the waveguide is perfectly bidirectional and when ν = +1 (ν = −1) the waveguide is perfectly unidirectional with system A (system B) upstream of B (A).

B. Remote two-qubit entanglement stabilization
A key motivation for our two chain setup is the ability to stabilize large amounts of remote entanglement.To understand the challenge here, we first review the simpler problem of dissipatively stabilizing entanglement between two remote qubits.One generic approach is to use squeezed light (as first introduced by Kraus et al [3], and further studied in [10,13,18]).Such schemes ultimately rely on generating correlated "pairing" dissipation, with a Lindblad jump operator ĉpair = uσ − A,1 + vσ + B,1 .While conceptually appealing, such pairing-based protocols are experimentally challenging, given the difficulty of preparing and propagating high quality squeezed light.Recent work shows that pairing dissipation can be realized without squeezed light, instead using modulated qubitwaveguide couplings [7,16]; this is also challenging in many setups.
The general N = 1 system has a pure steady state ρ1 (t → ∞) = |ψ 1 ⟩⟨ψ 1 | given by (up to normalization) where Γ ∈ C is a generalized complex detuning.Here we introduce the singlet and triplet entangled states in addition to the unentangled vacuum |00⟩ and the "doublon" |11⟩.Eq. ( 7) is the unique two-qubit steady state for any |Ω/Γ| < ∞, and as |Ω/Γ| → ∞, it approaches a perfect Bell state |ψ 1 ⟩ → |S⟩.Note, however, that in this limit the dissipative gap closes and a second impure steady state emerges [7][8][9].Going forward, one key goal will be to extend this entanglement stabilization to the case where each chain has N ≫ 1 qubits.As we show, this is a priori a non-trivial exercise.Unlike schemes that use squeezing dissipation, the entanglement structure of the N ≫ 1 chain is complicated and crucially depends on the interplay between the boundary driving and lattice dynamics.However, this interplay also gives the N ≫ 1 chain advantages over squeezing based schemes by enabling large amounts of entanglement to be stabilized without requiring Ω ≫ γ.To describe the stabilized state for N ≫ 1, it will be useful to use a basis where we specify the state of each cross chain dimer j.We will use the notation like |(01) j ⟩ = |0 A,j 1 B,j ⟩ to denote the state of the qubit pair on dimer site j.For example, |S j ⟩ denotes a Bell pair spanning qubits Aj and Bj.

C. Boundary driven dissipative quantum spin chain
A second motivation for our work comes from the seemingly simpler single-chain system depicted in Fig. 1(b): a chain of XX-coupled qubits with local loss and Rabi driving on one boundary site of the chain.The master equation for this system is As with other boundary-driven spin chain systems, we are interested in understanding the NESS of this setup.
For weak drives one expects the NESS approaches a product state (all qubits in |0⟩), whereas for strong driving, one instead expects an infinite temperature state.How one interpolates between these limits (and how the corresponding NESS depends on master equation parameters and the disordered hopping) is at first glance unclear.Surprisingly, Eq. ( 10) cannot be mapped to a system of free fermions, even though this is possible for the Hamiltonian ĤA alone.As we show in App.A, upon making a Jordan-Wigner transformation, the Rabi drive in ĤA yields a linear-in-fermions term, something that can be treated using the method of Ref. [42].However, when applied to the full master equation, this necessarily yields an interacting fermionic problem.Specifically, the loss dissipator becomes nonlinear in the fermionic master equation: where Ĥfermi,η is a quadratic fermionic Hamiltonian, ĉ1 is the fermion lowering operator on the dissipative site, and η is the fermionic lowering operator of an auxiliary mode introduced per the method of Ref. [42] (see App.A for details).
Eq. ( 10) thus corresponds to an interacting fermionic system without translational invariance (the hoppings can be disordered).Surprisingly, we are able to find an exact analytic description of its NESS.We do this by first solving for the pure steady state |ψ⟩ of the double-chain system in Eq. ( 1), something that can be be done analytically.If we then focus on the case where the waveguide is directional from A to B (i.e.ν = 1), simply tracing out the B chain yields the steady state of the single-chain system in Eq. (10).
This corresponds to a new many-body application of the coherent quantum absorber technique introduced in Ref. [4] and extended in Refs.[36,37,43,44].The existence of this exact solution implies that the single-chain system has a "hidden time reversal symmetry" which enforces Onsager time-symmetry of a certain class of twotime correlation functions [37].
As presented in more detail in Sec.IV, our exact solution for this boundary-driven spin chain reveals a number of surprising features in the NESS, including regimes of strong long-range correlations and even structures reminiscent of charge density wave order.

III. PURE ENTANGLED STEADY STATE FOR ARBITRARY N AND Ω
We now introduce a key result of this work: the boundary-driven double spin chain in Eq. ( 1) has a pure steady state for arbitrary N , drive strength Ω and hoppings J j .Even though the 1+1-qubit system has generically a unique pure steady state, a priori there is no reason to expect that this will also be true when N > 1.Indeed, Fig. 2 shows that for a generalized version of our 2 + 2 qubit model, the steady state will be impure if the two hopping amplitudes differ, or if we detune the second pair of qubits from the drive.
We find surprisingly that these two conditions (mirror symmetry of hoppings, no detunings of additional qubits) are enough to guarantee a pure steady state for arbitrary N , and for arbitrary choices of the parameters in Eq. ( 1).In the infinite-drive limit, the steady state has a simple translationally-invariant dimerized form that can be understood from a general replication argument that we present below (and that applies to many other setups [16][17][18][19]).For finite drives Ω, the steady state has a far more complicated form that is neither dimerized nor translationaly invariant.Our exact analytic expression nonetheless provides a simple picture for the state: it is a condensate of paired "hole" excitations, where holes correspond to cross-chain dimers that are in the vacuum |(00) j ⟩ state.

A. Ω → ∞ limit: generic entanglement replication via XX couplings
The form of this pure state of Eq. ( 1) becomes extremely simple in the limit of strong driving Ω: This is a highly entangled state of the two chains, that factors into a product of Bell pairs on each cross-chain dimer j.
The existence of a pure steady state is special.We consider two kinds of Hamiltonian perturbations to the N = 2 version of the system in Fig. 1(a), demonstrating that the emergence of a pure steady state is not generic.We break the mirror symmetry in the hopping rates via JB = JA + δ (cf. the discussion below Eq. ( 5)).We also consider the addition of both equal detunings and opposite detunings to the second site of each chain: Ĥ∆ = (δ/2)(σ z A,2 ± σz B,2 ).The purity of the numerically computed steady state ρss = ρ(t → ∞) is shown in each case vs. the perturbation strength δ.The steady state is only pure when δ = 0. We use unperturbed system parameters Ω = 0.5γ, ∆ = 0.1γ, and J = 0.25γ.more general "entanglement replication" phenomenon associated XX couplings and definite-parity dimer states; we explain this in what follows.This mechanism also explains and unifies the replication phenomena seen in several previous works [16][17][18][19].We note that the general nature of the replication mechanism we present here was not discussed in earlier works.
Imagine, as in Fig. 3, that we have dissipative dynamics L acting on a 2-qubit system that stabilizes an arbitrary state |ψ⟩ with fixed excitation number parity, i.e. |ψ⟩ = a|00⟩ + b|11⟩ or |ψ⟩ = a|01⟩ + b|10⟩.The actual method for stabilization is unimportant, only the state matters.For now, we will focus on even parity states for concreteness, but the analysis is identical for odd parity ones.Note for our specific system, for the N = 1 case and Ω → ∞, the stabilized state has a definite parity, hence the following arguments apply.
Next, imagine passively coupling a second pair of qubits to the first with an XX Hamiltonian, c.f. Eq. (5) then |ψ⟩ is still a steady state of the dynamics The analysis follows in the exact same manner if one uses odd parity states instead of even parity ones.Elsewhere in this manuscript, we work in a gauge with uniform coupling signs (so that there is no factor of (−1) s in Eq. ( 5), contrast with Eq. ( 16)).The replication here goes through exactly the same, where we can make a local sign flip on every other dimer, moving the relative phase onto the b coefficient in Eq. ( 17), b → (−1) i b.Therefore, this replication argument gives a general proof that Eq. ( 12) is a pure steady state of Eq. (1) in the infinite driving limit (as in this limit, the N = 1 problem has a definite (odd) parity pure steady state |S 1 ⟩) [45].The change of gauge thus explains why Eq. ( 12) is a product of staggered |S⟩ and |T ⟩ instead of uniform |S⟩ states.
The fact that the Hamiltonian perfectly replicates fixed parity states can be understood intuitively from the fact that it commutes with Ŝ2 s the total spin operator and Ŝz s = σz s,1 + σz s,2 the collective Z operator of each chain s = A, B. More details are in App.B 1.More generally, one can demonstrate that given any arbitrary two-qubit entangled state, Heisenberg couplings can be used to perfectly replicate this state down the chain, alleviating the parity constraint.Details are in App.B 2. Moreover, for both the Heisenberg coupling or XX couplings, significantly more complex geometries than chains can be used, generalizing [18].For more details, see App.B 3.

FIG. 3.
Entanglement replication via passive exchange couplings.Suppose there is some dissipative dynamics L that stabilizes a two qubit entangled state |ψ⟩ on the first pair of qubits in two exchange-coupled chains.If that state has definite parity, then the product state |ψ1ψ2⟩ is a zero energy eigenstate of passive exchange couplings ĤXX (cf.Eq. ( 13)).Thus the product state is a steady state of the stabilization dynamics and exchange couplings.Therefore, any number of qubit pairs may be added via exchange couplings, resulting in the tensor product steady state |ψ1ψ2 • • • ψN ⟩, replicating the entanglement on site 1 down a pair of arbitrarily long chains.

B. Ω < ∞: pure steady-state condensate of paired holes
We now turn to the more general (and experimentally relevant) case of a non-infinite drive amplitude Ω.For finite strength driving, the steady state for N = 1, Eq. ( 7), no longer has definite parity.As such, the replication argument of the previous section does not apply when we now consider larger systems, and there are no general arguments that would guarantee the existence of a pure steady state for N ≥ 2. Remarkably, we find that for arbitrary parameters, Eq. (1) has a pure steady state, albeit one that is far more complicated than the dimerized, translationally invariant state in Eq. (12).As we now show, this state can be exactly understood as a condensate of paired "holes".Here, a "hole pair" excitation is defined as starting with the "filled" state in Eq. ( 12), and then replacing the state of two adjacent dimers with the vacuum state, i.e. |(S) j ⟩|(T ) j+1 ⟩ → |(00) j ⟩|(00) j+1 ⟩.
Certain features of our state can be understood intuitively.For finite Ω, it is reasonable to expect the presence of holes, i.e. fewer qubit excitations than in the infinite driving limit.Further, these holes should be delocalized throughout the chain in order to have an eigenstate of the kinetic energy term ĤXX .This motivates looking for a pure steady state having delocalized hole excitations (with the density of holes scaling inversely with drive amplitude).More unexpected is our finding that in the steady state, these holes must be paired on adjacent sites.
To present our solution, it is convenient to map the double chain system to a 1D "dimer chain", where each site of the new 1D chain has local Hilbert space dimension 4, and corresponds to a dimer of the original system.With this mapping, we introduce two flavors of "particle" where we implicitly define the dimer ladder operators τ † j , τj that create and destroy the dimer particles |• j ⟩ = |(S/T ) j ⟩ (for j odd/even) when acting on |• j ⟩ or |• j ⟩, respectively.Similarly λ † j , λj create and destroy the particles |■ j ⟩ = |(T /S) j ⟩ (for j odd/even) when acting on |• j ⟩ or |■ j ⟩, respectively.Note that there is a hard-core constraint that prevents a |•⟩ and a |■⟩ from simultaneously occupying a site.We can neglect the remaining basis state for each dimer for now, as this state does not appear in the pure steady state of interest (see App. C for more details).With our new representation, the filled state Eq. ( 12) is thus Using the dimer particle representation defined in Eq. ( 20), we introduce an operator that creates a delocalized hole-pair: Here is the RMS hopping rate.Acting on the reference "filled" state |ψ ∞ ⟩, the operator Q creates a superposition state where each term corresponds to a pair of adjacent holes at a different location [46].
At a heuristic level, the phases in Q will give each hole pair a net momentum, allowing them to become zeroenergy eigenstates of the "kinetic energy" ĤXX .More formally, we show in App.D that Q commutes with ĤXX .Its action on an ĤXX eigenstate will thus produce another eigenstate with the same energy.In particular: An entire tower of zero-energy ĤXX eigenstates can thus be generated by repeated application of Q to |ψ ∞ ⟩, each state having an increasing number of hole pairs.The maximal-hole state in this tower corresponds to the empty state (if N is even), or a state with a single delo- Recall that in general, to obtain a pure steady state we require a state that is both an eigenstate of the Hamiltonian, and annihilated by relevant dissipators.The family of states Qm |ψ ∞ ⟩ provide us with a large class of states First few terms of the steady state.The expansion of the steady state in the number of holes added to |ψ∞⟩ (cf.Eq. ( 12)).The expansion coefficients an ∼ 1/ Ωn are the weights of the nth-hole components of the state; they can be found analytically from Eq. ( 25).Here we take uniform Jj = J.For nonuniform Jj the components of each an term are weighted by factors of Jj/ J. that are compatible with the Hamiltonian and connected to the pure steady state in the infinite-drive limit.One might expect that they can be used to construct the steady state for the finite-drive amplitude system.We find that apart from a boundary correction, this is indeed the case.As we rigorously show in App.E, for any set of parameters, Eq. ( 1) has a pure steady state |ψ Q ⟩ given by the "pair condensate" state where τ1 is the dimer lowering operator that removes the particle on site 1 (cf.Eq. ( 20)) and the reference state |ψ ∞ ⟩ is given by Eq. ( 12).The dimensionless drive strength Ω appearing in the exponential is with Γ given by Eq. ( 8) and J by Eq. ( 23).The first few terms of |ψ Q ⟩ are shown in Fig. 4. Up to overall normalization, the coefficients a j can be read off from Eq. ( 25), e.g., for a uniform chain (J j = J) a 0 = 1, a 1 = Γ/Ω = Γ/ J/ Ω, a 2 = 1/ Ω2 , etc; this is a power series in 1/ Ω, with a j ∼ 1/ Ωj .As we will discuss in more detail, this exact solution also immediately lets us understand the NESS of the non-trivial single chain system in Fig. 1(b).
Finally, we note that there is an equivalent construction of |ψ Q ⟩ that is recursive in the length N of the chains.The recursive construction enables the efficient numerical evaluation of expectation values and correlation functions, e.g.particle density ⟨n j ⟩ = 1 2 ⟨τ † j τj ⟩.Details are provided in App.F.

IV. REAL-SPACE HOLE PAIRING AND CONSEQUENCES FOR THE STEADY STATE A. Hole pairing as a kinetic constraint
The hole-pairing in |ψ Q ⟩ is surprising at first glance, as there is no attractive interaction or other explicit pairing mechanism in our model.Hole pairing turns out to be a consequence of two facts: (i) ĤXX causes |•⟩ and |■⟩ particles to change flavor when they hop, and (ii) the hard-core constraint forbids a |•⟩ and an |■⟩ to simultaneously occupy a site.As shown in Fig. 5(a), with details in App.C, a particle can swap positions with a hole and change flavor in the process.Two adjacent particles of the same flavor cannot hop due to the hard-core constraint, hence ĤXX |• •⟩ = 0.If we now try to form zero kinetic energy states (i.e., zero-energy eigenstates of ĤXX ) by delocalizing hole excitations, we find that they must be paired on adjacent sites.Delocalizing a single hole to get a zero-energy eigenstate fails due to the flavorchanging hopping, see Fig. 5(b), but delocalizing a pair is successful (see Fig. 5(c)).
The creation of zero kinetic energy eigenstates of ĤXX via the hole-pairing operator Q (cf.Eq. ( 22)) is somewhat analogous to η-pairing found in Fermi-Hubbard lattices [47,48].In both cases, the pairing operator generates exact Hamiltonian eigenstates with zero kinetic energy by delocalizing a pair of excitations throughout the system.A key distinction however is that in η-pairing, each paired excitation is spatially local, i.e. a pair of fermions on one lattice site, whereas in Q-pairing, each hole-pair occupies adjacent lattice sites.We also note that the algebraic structure of η-pairing (η, η † form a closed representation of SU(2)) is not found in Q-pairing: Q and Q † do not form a closed SU(2) group.
While not a true symmetry of the Hamiltonian, the hole-pairing operator has a special relation to the Hamiltonian of the double-chain system, via a structure that was introduced in the context of many-body scar states.Specifically, it constitutes a restricted spectrum generating algebra (RSGA) of double-chain Hamiltonian when acting on |ψ ∞ ⟩ [39].Furthermore, certain terms can be added to the XX chain to make the model nonintegrable while preserving this RSGA structure, thus guaranteeing that Qn |ψ ∞ ⟩ remain exact eigenstates.In App.G we discuss how this makes the corresponding hole-pairing states true many-body scar states in a nonintegrable ladder system (a model related to that studied in Ref. [40]).
Having understood the route to hole pairing in our model, we can also postulate other Hamiltonian models where this will occur.For example, a 1D Fermi-Hubbard chain with a spin orbit interaction can exhibit effective flavor-changing hopping.For strong interactions, it can thus also exhibit hole-pairing in a subset of its eigenstates (see App. D).The fermionic analog of the hole-pairing operator Q has the same properties, generating eigenstates of the 1D chain when acting on a filled state.These holepaired eigenstates may be accessible in e.g. the ultracold atoms platform proposed in Ref. [38].

B. Density correlations due to hole pairing
The structured hole pairing in Eq. ( 25) immediately gives rise to spatial density correlations, something that is most apparent when the hole density is low, i.e. | Ω| ≳ 2. This correlation is directly observable in the single-chain system of Fig. 1(a) as Z magnetization correlations between adjacent sites, ⟨σ z A,j σz A,j+1 ⟩ ̸ = 0.The correspondence between magnetization and hole density follows from the fact that the dimer holes |•⟩ are polarized but the particles |•⟩ are depolarized (cf.Eqs. ( 19) and ( 20)): Thus, (−σ z A,j ) acts as a local hole number operator when acting on the steady state.We define the zmagnetization correlation function for the A chain as normalized such that C zz (j, j) = 1.See App.H for a discussion of the nonstandard normalization.In Fig. 6, we show C zz (j, k) for fixed distance |k − j| and aver- aged over the whole chain.For strong driving, the average correlations between adjacent sites saturates to C zz (j, j ± 1) → 0.5 because in this regime, a hole on site j is always paired with a hole on either j ± 1, but is unlikely to be correlated with any other site.In contrast, there are no appreciable correlations in this regime for larger distances.We note that the correlation functions C zz (j, k) can be measured in either the double-chain system or the single-chain system.In either case, the saturation of C zz (j, j ±1) ≈ 0.5 for Ω > 1 is a clear indication of hole pairing.

C. Universal density scaling
There are two dimensionless parameters in our model: the dimensionless drive amplitude Ω (c.f.Eq. ( 26)) and the ratio ζ ≡ Γ/ J.One might naturally expect that bulk steady state properties would depend on both these parameters.However, the exact result Eq. (25) shows that this is not the case.We can write this as which suggests that the excitation density in the bulk is controlled only by Ω.We now show explicitly that the excitation density (i.e.Z magnetization density) is indeed intensive, and scales universally with the single parameter | Ω| in the regime | Ω| ≳ 2. Consider first the limit where |ζ| → 0 while | Ω| remains fixed, such that we can ignore the boundary term in Eq. ( 29).We thus have This expression mimics a bosonic coherent state |α⟩ = e αâ † |0⟩, where hole pairs are the bosons, |ψ ∞ ⟩ is the hole-pair vacuum and Q is the hole-pair cre-ation operator.As we show in App.I, this analogy to bosonic coherent states can be made precise when N ≫ 1 and the hole density m satisfies In this low hole density regime, we may neglect the hardcore constraint that prevents two hole-pairs from occupying the same sites, and can thus accurately estimate the total hole number as We thus find for large drives an intensive scaling of hole density.Note that this result immediately implies that for an arbitrarily long single chain system, one only needs | Ω| ≳ 2 to approach the infinite temperature state.When |ζ| > 0, the coherent state analogy is still valid, but there is a boundary correction to Eq. (31).Because site 1 directly sees the dissipation, it can be occupied by an isolated hole and thus its occupation depends not only on | Ω| but also |ζ|.The hole density m thus obtains a O(1/N ) correction to the intensive universal scaling There is thus a transition from universal ∼ | Ω| −4 scaling to an N -dependent ∼ | Ω| −2 scaling at | Ω| 2 ≈ N/|ζ| 2 , or equivalently Ω/ J ≈ √ N .The universal scaling behavior is exact for any |ζ| when the dissipative site density is excluded, as is shown in Fig. 7(a), and the dissipative site corrections for |ζ| > 0 are shown in Fig. 7(b).

D. Single particle "charge density waves"
For J j = J, the spin chain systems in Fig. 1 are translationally invariant except for the boundary j = 1 site.For long chains, one might thus expect that the steady state density is also translationally invariant, except for edge effects near j = 1.Surprisingly, this expected translational invariance is strongly broken in the steady state for weak drives Ω.As a consequence of the hole pair condensate of Eq. ( 25), there is a regime where the steady state corresponds to a single excitation that is localized on either the even-j or odd-j sublattice, i.e. a kind of single particle charge density wave (CDW).
As we will show, for weak drives, the double spin chain system's steady state (for uniform J j = J) is given by the CDW form: where ⌈X⌉ is the ceiling function, and |• j ⟩ is given by Eq. ( 20).This describes a single particle that is delocalized either across all odd-j sites in odd-length chains or across all even-j sites in even-length chains.For disordered J j , each component |• j ⟩ is weighted by an additional factor (• • • J j−4 J j−2 )(J j+1 J j+3 • • • )/ J⌈(N−1)/2⌉ (and the requisite correction to the normalization).
For the single dissipative spin chain, Fig. 1(b), the nonequilibrium steady state ρA,cdw = tr B |Ψ cdw ⟩⟨Ψ cdw | is an equal mixture of vacuum and the single-particle CDW: Despite not being a pure state, the single chain CDW is a very low entropy state, S(ρ A,cdw ) = ln 2, for which any particle density necessarily arises from the single coherently delocalized excitation.Here, |Φ cdw ⟩ is the single chain equivalent of |Ψ cdw ⟩ given above, and its components have the same weighting factors when the J j are disordered.
To see why these single particle states emerge, consider the steady state, Eq. (29), in the limit |ζ| → 0: In this limit, all holes in the state are created by some power of Q acting on |ψ ∞ ⟩, with the component Qm |ψ ∞ ⟩ having 2m holes.For an N -length chain, Q can thus act up to m = ⌊N/2⌋ times to produce new states, after which Qm>⌊N/2⌋ |ψ ∞ ⟩ = 0.For an even chain, all 2m = N particles can be removed, taking the Eq. ( 38)) for two different |ζ| 2 = 10 −4 and |ζ| 2 = 10 −2 .In the limit |ζ| → 0 (dashed curves), the single particle CDW persists for any arbitrarily small | Ω| > 0. In both plots, the particle density is n = 1 filled state to vacuum: However, for an odd chain, all but 1 of the N = 2m + 1 particles can be removed.The real-space pairing of holes on neighbouring sites requires that the remaining single |•⟩ particle is confined to the odd-j sites, resulting in a CDW-like structure where a single particle is delocalized over the odd-j sublattice only: Thus for |ζ| ≪ 1, there is a regime of sufficiently weak | Ω| for which odd length chains exhibit a charge density wave consisting of a single delocalized particle.The emergence of a CDW in an odd-length chain, and the lack of a CDW in an even-length chain for the same parameters, is shown in Fig. 8(a).The particle density n = N | Ω| 4 /8 of the even chain is found using Eq.(J6) in App.(J).
The analysis follows analogously in the limit |ζ| → ∞, τ1 |ψ ∞ ⟩ (as τj commutes with Q).Site 1 thus always has a hole, and we repeat the above analysis on the remaining N − 1 sites.Thus, even chains have a CDW-like structure where a single particle is delocalized over the even-j sites.
The parameter regime in which CDWs emerge can be found by expanding |ψ Q ⟩ to the lowest few orders in | Ω|.Leaving the details to App.J, we can show that two distinct scales emerge: a N -dependent upper limit on drive strength, | Ω| 2 ≪ 1/N , and a |ζ|and N -dependent lower limit that differs for even and odd length chains: When | Ω| 2 > 1/N , many particles can populate the chain, destroying the CDW ordering, and when | Ω| 2 < |ζ| ±2 /N (for N odd/even), the average particle density vanishes as ∼ | Ω| 2 .Here, Eq. ( 34) is no longer an equal mixture of |0⟩ and |Φ cdw ⟩, but increasingly weighted toward vacuum with | Ω| → 0. The emergence of CDWs, and the dependence with both N and |ζ|, is shown in Fig. 8(b).

V. RESOURCE FOR CROSS-CHAIN (REMOTE) ENTANGLEMENT STABILIZATION
As discussed in Sec.II, the double qubit chain system in Fig. 1(a) (where collective loss is provided by passive couplings to a waveguide) is a potentially powerful setup for stabilizing large amounts of steady-state remote entanglement.The scheme is also very resource efficient: it foregoes the complication and resource overhead of using squeezed light in favor of local driving, and it only requires passive hopping between qubit in each chain.It also does not require precise fine-tuning of parameters, nor does it require extremely strong driving to approach maximal entanglement of the chains.We now use insights obtained from our exact solution Eq. ( 25) to better understand this potential application.
The exact solution tells us that for any drive strength Ω, we will have a pure steady state with some degree of entanglement between the remote chain-A and chain-B qubits.This entanglement is maximal in the Ω → ∞ limit, where the steady state becomes a dimerized product of cross-chain maximally entangled Bell pairs.A natural question is how strong must the Rabi drive be to achieve this level of entanglement.The exact solution provides a succinct and surprising answer here: one only needs that the effective drive amplitude | Ω| ≳ 2, as in this regime the density of "holes" is very small, implying the steady state is very close to the ideal dimerized state.We stress that this condition is independent of N (even though drives are only applied to the first qubit in each change), and further, that one can achieve this condition even if Ω ≪ γ (the drive does not need to overwhelm dissipation if hopping is sufficiently weak).
Of course, these considerations neglect a crucial second issue: one cares both about the amount of entanglement in the dissipative steady state, as well as the time needed to prepare this state (i.e. the characteristic system relaxation time, or the inverse dissipative gap).This timescale will also directly determine the susceptibility of our scheme to additional unwanted dissipative processes (e.g.waveguide loss, qubit dephasing and relaxation).
A full study of the effects of waveguide loss and qubit dissipation on entanglement stabilization in a circuit QED realization of Fig. 1(a) is presented in a complementary work [49], but we briefly discuss the basic requirements to realize the scheme in circuit QED in App.K.
Here, we instead focus on a fundamental aspect of the relaxation time physics in our double chain scheme.While the qubit-only version suffers from a fundamental tradeoff between speed and entanglement, we show below that by generalizing the local 2-qubit scheme introduced in Ref. [8] to a multi-qubit setup with directional dissipation, we can dramatically improve this seemingly unavoidable tradeoff.

A. Slowdown in large drive limit
Recent works [6][7][8][9] on dissipatively preparing entanglement between two remote qubits have observed that the dissipative gap closes as the target entangled state approaches a perfect Bell pair (e.g., as the vacuum component of Eq. ( 7) vanishes).Refs.[8,9] showed that this is a generic property of two-qubit systems, and that it arises due to an approximate conservation of total angular momentum that becomes an exact symmetry in the infinite driving limit.In this limit, a second impure steady state emerges in the subspace orthogonal to the target Bell state.As one approaches the point of added symmetry, the transition rate out of the orthogonal subspace and into the target Bell state becomes extremely small, leading to a vanishing dissipative gap.We briefly review this argument in App.L.
In the infinite driving limit |Ω/Γ| → ∞, the steady state of the N = 1 double chain is not unique [7][8][9], but can be any state of the form ρ1 = ν|S⟩⟨S| + 1 4 (1 − ν) Î for any −1/3 ≤ ν ≤ 1, see App. M for details.One can readily show that the maximally mixed state Î/4 is replicated via the XX Hamiltonian.Therefore, the near-symmetry that causes a slowdown in the N = 1 system persists for N > 1 because the near-steady infinite temperature state is replicated down the chain, and thus the chain cannot relax out of that state except by the very slow dissipative population transfer at the boundary.12)), 1 − ⟨ψ∞|ρ|ψ∞⟩ (black), are shown as functions of the hopping rate J/γ for N = 3, Ω/γ = 10, and ∆ = 0.For the qubit-qutrit scheme, the relaxation time is optimized over η.We also plot the relaxation time for the single chain system for comparison.For a fixed state fidelity of 0.999 (achieved at J = 7γ, dashed line), the relaxation times are γτ 2qb = 2.1 × 10 4 , γτqutrit = 120 and γτsc = 14, for the qubit-qubit, qubit-qutrit, and single chain, respectively.

B. Speeding up stabilization with a qutrit
A recent work by Brown et al. [8] theoretically proposed and experimentally demonstrated that the slowdown associated with dissipatively stabilizing 2-qubit Bell pairs can be circumvented by promoting one of the qubits to a qutrit, in a system for which the dissipation is both local and reciprocal (i.e.mediated by common coupling to a damped cavity mode).They demonstrated that the near-symmetry that conserves total angular momentum is no longer present in a qubit-qutrit system.This makes the degenerate dark state vanish in the large drive limit.
Here, we show that this scheme can now be extended to the directional version of our double chain system by promoting the down-steam qubit B1 to a qutrit, see Fig. 9(a).This leaves Eq. ( 25) as the pure steady state while avoiding the symmetry-induced slowdown, and allows a dramatic stabilization speed up for arbitrary N without sacrificing the fidelity with the perfect dimerized entangled state.
More concretely, starting from the double chain master equation (cf.Eq. ( 1)) for N = 1 in the directional limit ν = 1, we promote qubit B1 to a qutrit and modify its coupling to the waveguide via for which the master equation now reads Physically, this means that now the qutrit B1 can produce a photon in the chiral waveguide either via a 1 − 0 relaxation event or a 2 − 1 relaxation event (with relative matrix elements η).The result is a dissipative interaction that allows the state |11⟩ to pass a single photon through the waveguide at a rate ηγ and become |02⟩, which can in turn decay into the state |01⟩.The effective interaction (no jump Hamiltonian) of such a process is −iηγσ − A (|2⟩⟨1| B ).Because this process explicitly breaks the conservation of angular momentum in the 2 qubit subspace, it circumvents the slow down previously observed.
From here, one can take our qubit-qutrit system and add back the remaining N − 1 qubits in each chain, and the hopping Hamiltonian ĤXX .We stress that all remaining qubits are just qubits: it is only B1 where we need to make use of the higher |2⟩ level.A direct calculation shows that the steady state found in Eq. ( 25) is still a zero-energy eigenstate of the new Ĥdrive + Ĥdiss , as well as the new jump operator L, and so it remains a dissipative steady state; the only thing that has changed in adding the third level is the dynamics, which should now be significantly faster.This is demonstrated in Fig. 9(b) for an N = 3 system.Numerically, we observe a significant improvement in the relaxation time scale τ rel (as determined by the inverse dissipative gap of the full Lindbladian) of a N = 3 system, when we promote site B1 to a qutrit and optimize over the qutrit 2-1 transition coefficient η.These results are as shown in Fig. 9(b).Here, we fix all other parameters except the uniform hopping rate J, and show how both the fidelity with the ideal dimerized entangled state and τ rel vary with J.For the all-qubit chain, small J makes | Ω| large, thus the fidelity to the maximally entangled state |ψ ∞ ⟩ (cf.Eq. ( 12)) is high, but the relaxation slows down.As Fig. 9(b) shows, there is a dramatic improvement of the relaxation time: over two orders of magnitude at a state fidelity of 0.999.The optimized value of η as a function of J/γ, as well as the speed-up of an N = 2 system, is shown in App.N. We expect that the qutrit scheme speeds up the stabilization time for larger N systems as well.

VI. CONCLUSION
Our work presents an exact analytic solution for the steady state of two different spin chain models with boundary dissipation and driving.As discussed, these solutions reveal a number of surprising correlation effects (e.g. the effective pairing of holes), and lay the groundwork for a potentially powerful route to dissipative stabilization of remote multi-qubit entanglement.We also elucidated a general mechanism for "replicating" definiteparity two-qubit entangled states using passive XX couplings in a double qubit chain, and demonstrated that the approach of Ref. [8] for avoiding slowdowns in dissipative entanglement stabilization could be extended from a 2 qubit situation to a setup with many qubits and directional dissipation.
In future work, it will be extremely interesting to explore whether the ideas introduced here could be extended to more complex systems, where multiple 1D XX qubit chains are attached to the same common waveguide.This could potentially be a source of stabilized multi-partite, multi-qubit remote entanglement.It would also be interesting to explore further the dynamics of our solvable dissipative spin chain models.As discussed, the solvability of the non-trivial single chain model in Fig. 1(a) can be ultimately traced to a surprising hidden time-reversal symmetry [37].Understanding how this symmetry constrains the dynamics and Liouvillian spectrum could be an extremely rich direction for future research.It would also be interesting to understand whether the scaling of the dissipative gap in the two chain model of Fig. 1(b) could be improved beyond the usual 1/N 3 scaling that is found in a variety of integrable spin chain models [16,29].Finally, it would be interesting to study the spectra and NESS of other two-chain models.We already demonstrated in App.G that our hole pairing states correspond to many-body scar states in a closed, non-intergrable ladder system.Adding dissipation here could be extremely interesting.Further, one could extend our two-chain model to a ladder system with Creutz ladder style couplings [50] along the full length of the chain.If one tunes the diagonal interchain couplings t d to be equal to the intra-chain XX couplings J, then this system possesses at least N − 1 strong symmetries.It thus has multiple dissipative steady states, making it another interesting system worthy of further study.10.Lindblad spectrum of a single chain with 2 qubits.We plot the eigenvalues λ of the single chain master equation Eq. ( 10) for N = 2, ∆ = 0 and J = γ.Left: In the absence of any drive Ω = 0, the spectrum has the normal mode form expected for a quadratic fermionic Lindblad master equation.Right: For a non-zero drive Ω = 0.6γ, the normal mode structure is lost, e.g.summing single-excitation eigenvalues does not correctly predict higher lying eigenvalues.This provides a direct confirmation that with driving and loss, the single chain problem is not equivalent to a quadratic fermionic Lindbladian.
Now, the Majorana (η + η † ) is conserved by the Hamiltonian, and so if this were a closed system we would be done as the two Hamiltonians Ĥfermi,η and Ĥfermi would be isospectral (with Ĥfermi,η doubly degenerate).However, this is an open system, and we need to consider the dissipation.In this case, the Majorana (η + η † ) only constitutes a weak symmetry [53] of the Lindbladian as it anticommutes with the jump term: This is problematic, as without further modification, the dissipation will cause unphysical jumps between the two conserved sectors of Ĥfermi,η .
To correct this problem, we must also modify the linear-in-fermion jump operator so that these unphysical jumps do not occur.Formally, we must make the conservation of (η + η † ) a strong symmetry [53].This is achieved by the master equation Note that this is ultimately equivalent to first introducing an auxilliary spin in Eq. ( 10) (preceding the first lattice site), and then performing the Jordan-Wigner transform.Note that the jump operator is now cubic, and therefore the system is explicitly interacting.We also note that this procedure is consistent with the general rules outlined in Ref. [42]: when introducing the auxiliary fermion η, all linear-in-fermion operator terms must be modified to ensure that they have the correct matrix elements in the expanded space.This rule must be applied to the jump operator ĉ1 , as the action of the superoperator D[ĉ 1 ] cannot be written solely in terms of the quadratic operator ĉ † 1 ĉ1 .As a final confirmation that the single-chain qubit system is not equivalent to free fermions, in Fig. 10, we plot the full Liouvillian spectrum (eigenvalues λ) for the N = 2 version of the master equation Eq. (10) For Ω = 0 (left panel), the eigenvalues have the normal mode structure expected for a free fermion Lindblad master equation (as can be found using third quantization [24]).This implies, e.g., that eigenvalues corresponding to two particle excitations are formed by summing eigenvalues associated with single particle excitations.With a non-zero drive (right panel), this structure is clearly lost.The Lindblad spectrum with both drive and loss no longer has the form expected for a free fermion Lindbladian, consistent with our conclusions above.

Appendix B: Replication 1. Intuition for fixed parity states
The fact that the passive XX couplings can replicate fixed parity states is simple enough to check, but is not immediately intuitive.To get some more intuition, let's write the Hamiltonian for a pair of dimers as We can immediately observe that [ Ĥ, Ŝ2 A,B ] = [ Ĥ, Ŝz A,B ] = 0, and so we can diagonalize it in terms of the singlet and triplet states |s A , m A ⟩ ⊗ |s B , m B ⟩ where s, m are the standard quantum numbers of total spin and z-angular momentum.It is quite easy to observe that (σ x The only remaining piece of the puzzle is the observation that tensoring together two copies of a fixed (even) parity state |ψ⟩ is diagonal in this spin basis.
where in the last line we go from the computational basis in Eq. (B6) to the total spin basis in Eq. (B7).Since the coefficient ψ i ψ j is a symmetric tensor, it is diagonal in the spin basis; intuitively, this is because the spin states are all either symmetric (s = 1) or antisymmetric (s = 0).Thus, the tensor product of s = 0 with s = 1 is antisymmetric, which multiplied by a symmetric tensor is identically zero.However, the product of two symmetric or two antisymmetric tensors won't be.This can also be checked very straightforwardly by direct computation.
as expected.
Going to a fixed (odd) parity state can be deduced in the same manner by observing that (defining a bit flip operation via |i⟩ ≡ |1 − i⟩ = σx |i⟩) However, bit flip on the entire B chain leaves the Hamil- and hence the argument still holds.This shows that the XX Hamiltonian annihilates every tensor product of identical fixed parity states.Thus, we can repeat this argument for an arbitrarily long chain of identical, fixed parity states.If we denote ĤXX,i as the XX Hamiltonian acting on the i and i + 1 pair of spins, then the total Hamiltonian would be so the total 2N -site XX Hamiltonian has as a many-body zero-energy eigenstate composed of N copies of an arbitrary fixed parity 2-qubit state.

Heisenberg Couplings
To observe that Heisenberg couplings can replicate any state, it is first important to observe that, given a fixed (even) parity state |ψ⟩ = u|00⟩ + v|11⟩, then σz Combining this with the fact that, as shown in the previous section, this state is annihilated by the XX couplings, we see that it is annihilated by any XXZ Hamiltonian.
From here, we point out that given an arbitrary state |ϕ⟩, then by the Schmidt decomposition, there exist local unitaries Û1 , Û2 such that for some u, v. Let's define ĤH to be the isotropic Heisenberg Hamiltonian.Then we have that Since the isotropic Heisenberg Hamiltonian is invariant under uniform local unitary rotations, it can be commuted through to annihilate the fixed parity state: as desired.

Replication in more complicated graphs
We will now demonstrate the claim made in the main text that the replication mechanism can work on more complicated graphs than a 1D chain.In fact, states can be replicated down any tree-like structure (i.e., with no closed loops) that has exactly one symmetry axis.
The proof is a simple generalization of what we have already shown: we will show that, given two dissipatively stabilized qubits, one can attach an arbitrary number of qubit pairs off of these using symmetric XX couplings.This generates one level of the tree graph, and then simple bootstrapping shows that arbitrary trees are possible.
Let's assume once again that there exists a Liouvillian operator L 0 acting on qubits at site A 0 , B 0 that stabilizes a fixed parity state |ψ⟩, L 0 (|ψ⟩⟨ψ|) = 0. Next, we will extend this to the next layer in the graph by defining the next set of qubits on sites A 1 , . . ., A n and B 1 , . . ., B n , so that the full system Liouvillian is now where we define sgn(A) = 1, sgn(B) = −1.Now at this point, defining |Ψ⟩ = |ψ⟩ ⊗n+1 , it is simple to observe that since the Hamiltonian is simply a sum of terms acting independently on the different dimers, ĤXX,n |Ψ⟩ = 0 by the exact same logic as presented before.Thus, a single qubit can sustain any number of pairs branching off of it.However, this means that each of those are now dissipatively stabilized fixed parity states, and so we can repeat the argument to branch more qubits off, generating trees.
At this point, it is crucial to note that, since there are multiple qubits branching off a single pair, now it is important that the J i terms are all distinct if you want a unique steady state.If there are degeneracies in the J i parameters, then one generates a permutation symmetry where there are multiple degenerate pairings between A qubits and B qubits, and so the steady state will necessarily be degenerate.

Appendix C: Dimer representation
Each dimer chain site can be described as a pair of non-commuting spin-1's embedded in the spin-1 2 ×spin- We define the lowering operators τj and λj on each site j that each destroy one of |S j ⟩ or |T j ⟩.The τj destroy |S j=1,3,5,... ⟩ and |T j=2,4,6,... ⟩ -the |•⟩ particles -and the λj destroy the opposite m = 0 states -the |■⟩ particles.Explicitly, in terms of the qubit operators, the lowering operators are: Each ladder operator and its adjoint forms a spin-1 representation of SU(2), with the same Ŝz j completing the algebra for each.That is, the commutation relations [τ † j , τk ] = [ λ † j , λk ] = 2 Ŝz j δ jk are simultaneously satisfied for the operator The two flavors of ladder operators on each site j do not commute with each other: = 2 (|S j ⟩⟨T j | + h.c.) .
Notice that these commutators are Hermitian and act to change flavor on a given dimer site.Finally products of same-site ladder operators are: This completes the dimer algebra.
We use the following notation of states in the dimer representation: |• j ⟩ and |(11) j ⟩ are the vacuum and double-excited state, respectively.The single-excited states are Finally, we define the steady-state subspace Now we rewrite Eq. ( 1) in the dimer representation.The collective loss operator (cf.Eq. ( 3)) is simply which makes immediately obvious that its dark subspace on site 1 is spanned by The Hamiltonian terms (cf.Eqs. ( 4), (5), and ( 6)) are ĤXX = 1 4 ĤXX becomes a flavor-changing exchange interaction and the Rabi drives become a single drive acting on the |■⟩ flavor.The drive detuning in Ĥdrive and the nonreciprocityinduced exchange in Ĥdiss appear in different ways to achieve the same effect: a flavor change that swaps |• 1 ⟩ and |■ 1 ⟩.Note that the commutator in Ĥdrive is Hermitian, which can be seen using Eq.(C4).
Appendix D: Hole-pairing and XX eigenstates in qubit chains and Fermi-Hubbard chains

Hole-pairing operator
The hole-pairing operator Q defined by Eq. ( 22) is a central character in the analytical description of the steady state of Eq. ( 1).Here, we prove that it acts on eigenstates of ĤXX (cf.Eq. ( 5)) within the steady-state subspace H s (cf.Eq. (C7)) to produce new eigenstates of ĤXX with the same energy and also in H s .We prove this at the operator level by showing that when Q acts on states within the steady-state subspace H s , it commutes with both ĤXX and the projector into the subspace Ps (cf.Eq. (C8)), Together, these imply that for a given ĤXX eigenstate |ϕ⟩ ∈ H s with energy E, the action of Q on |ϕ⟩ produces another eigenstate with energy E: (D3) In this way, we find that by repeated application of Q on an ĤXX eigenstate, a tower of some finite number n 0 < ∞ of degenerate eigenstates is returned.

Hole-pairing in a Fermi-Hubbard model
The hole-pairing observed in the dimer chains is not unique to that system.The requirements for a 1D tightbinding chain to have eigenstates of paired holes (or some suitable paired excitation) are that (i) there exists 2 flavors or species of excitation, (ii) the exchange couplings change the flavor when the particles hop, and (iii) there is a hard-core constraint preventing double occupation on a site.The Fermi-Hubbard tight-binding chain with onsite repulsion U can, in the hard-core interaction limit U → ∞, fulfill these requirements, albeit in a nonstandard basis.A recent work by Mamaev et al. [38] introduced a proposal of a Fermi-Hubbard model with a staggered laser drive that flips spin on each site.This laser drive induces an effective spin-orbit coupling that causes spin-flip (flavor-change) hopping in an appropriate excitation basis.In the hard-core repulsion limit, hole-paired states emerge as eigenstates of the Hamiltonian.
The model introduced in Ref. [38] is of a 1D Fermi-Hubbard lattice with an added spin-orbit coupling (SOC) laser drive, Ĥ = Ĥ0 + ĤSOC , where Here J is the hopping rate (taken to be uniform for simplicity, but can be generalized to non-uniform J, like the spin chains), U is the Hubbard potential, Ω is the laser Rabi drive strength, and nj,σ = ĉ † j,σ ĉj,σ .As in Ref. [38], we define a new set of fermion operators which obey the canonical anticommutation relations {â j,σ , âk,τ } = δ jk δ στ .In this basis, the Hamiltonian is where σ denotes the spin-flip of σ and now nj,σ = â † j,σ âj,σ .We thus have the desired spin-flip hopping.Also note the energy splitting Ω between spin-up and spin-down particles, thus making this a good basis with respect to the SOC laser drive.Whereas Ref. [38] explores the physics in the limit Ω = U → ∞, here we wish to consider a slightly different limit: U → ∞ with Ω, J < ∞.This effects the on-site hard-core repulsion between opposite spins, thus we have the effective model With the hard-core repulsion, the N -particle ferromagnetic states eigenstates of Ĥ∞ , with energies +N Ω 2 and −N Ω 2 , respectively.We now introduce the fermionic hole-pairing operators (one for each spin) that create pairs of adjacent holes with staggered phases, in perfect analogy with Eq. ( 22).We also denote by Pσ the projector into each spin subspace, in analogy with Eq. (C8).Just as with the dimer spin chain case, we find that the hole-pairing operators here commute with the tightbinding Hamiltonian when restricted to the appropriate subspace, but the total energy changes due to the removal of two particles: Thus, when acting on their respective |Ξ σ ⟩, powers of Qσ produce towers of Hamiltonian eigenstates with energies ranging from ±N Ω 2 to either 0 or ± Ω 2 for either evenlength or odd-length chains, respectively.Notice that due to the SOC energy splitting, Qσ raises or lowers the eigenstate energy by a multiple of the SOC drive strength Ω/2, furthering the analogy with η-pairing [47,48], for which η-paired states have energy splitting that are a multiple of the Hubbard potential U .
Appendix E: Existence proof for the pure steady state Here, we will rigorously prove that |ψ Q ⟩, Eq. ( 25), is indeed a pure steady state of Eq. ( 1).For |ψ Q ⟩ to be a pure steady state, it is sufficient for it to be an eigenstate of the Hamiltonian and a dark state of the dissipation [54].The latter property is satisfied as ĉ|• 1 ⟩ = ĉ|• 1 ⟩ = 0 and |ψ Q ⟩ has only those components on site 1.It remains to show that |ψ Q ⟩ is an eigenstate of Ĥ, which we do by way of a variational ansatz.
We consider a two-parameter variational ansatz |ψ ′ [α, β]⟩ and show that it is an exact eigenstate for a set of uniquely determined variational parameters α, β.First note that for N = 1, the steady state (cf.Eq. ( 7)) is This is a zero-energy eigenstate of the boundary Hamiltonian Ĥ1 ≡ Ĥdrive + Ĥdiss (E2) (cf.Eqs.(C11) and (C12)).Thus we construct a variational ansatz which includes this wavefunction in the N = 1 case, but we replace Γ/Ω with a variational parameter β for generality.This correction to |ψ ∞ ⟩ cannot be enough for N > 1 because the hole on site 1 will be delocalized throughout the chain in order to have an eigenstate of ĤXX .One might expect that a linear combination of zero energy ĤXX eigenstates with all possible number of hole pairs is needed.A "hole-pair condensate" i.e., the exponential of Q acting on |ψ ∞ ⟩, is a simple way to achieve that, thus we make the ansatz For all n ≥ 3, the steady state is given recursively by

Correlation functions
The recursive construction of |ψ⟩ provides a convenient way to numerically compute correlation functions.Using Eq. ( 25) directly presents some analytic challenges that are avoided in the recursion.For the sake of clarity, we focus here on the dimer particle number nj ≡ 1 2 τ † j τj .Since the recursion relation for |ψ n ⟩ is in the chain length n, it is necessary to denote the chain length for which expectation values are taken, e.g.
Here the state is normalized by construction, ⟨ψ n |ψ n ⟩ = 1.Using the recursion relation Eq. (F6), we expand ⟨n j ⟩ n in terms of the expectation values evaluated for n−1 and n − 2 length chains: where N k are the normalization factors given by Eq. (F7).This expression tells us that the expectation value ⟨n j ⟩ n , evaluated for a chain of length n, is given in terms of the expectation value evaluated on chains of length n − 1 and n − 2, which are similarly given by the recursive expression Eq. (F9).This recursion terminates at the expectation value ⟨n j ⟩ j , i.e. for a length n = j chain.The termination ⟨n j ⟩ j is readily evaluated using Eq.(F6) directly and is Any expectation value or correlation function can be evaluated in a similar way.
Appendix G: Restricted spectrum generating algebra and relation to scarring Here we briefly demonstrate that Q and |ψ ∞ ⟩ (cf.Eq. ( 12)) constitute a restricted spectrum generating algebra (RSGA), and show that these states have arealaw entanglement entropy across a particular bipartition, and thus are scar states of a simple ladder model.
We consider a simple nonintegrable ladder model, specifically the model discussed in Ref. [40] of two XX chains with an XXZ interaction on each rung (to which we add a term Ĥδ ): It is straightforward to see that the Q operator constructed in App.D and the reference state |ψ ∞ ⟩ (cf.Eq. ( 12)) form a RSGA of order 2 under this Hamil-tonian.In particular, we have where Ĥn = [ Ĥn−1 , Q] with Ĥ0 = Ĥ.Eq. (G4) implies that the set of hole-paired states are a tower of exact eigenstates of the nonintegrable ladder, with equal energy spacing E n − E n−1 = 2nµ, for all n < L/2.Similar calculations show that Q † acting on |ψ ∞ ⟩ is also a RSGA of order 2, but this time creating particle (i.e., |(11) j ⟩) pairs.Similarly, if we define a Q which is identical to Q but with the replacement τj τj+1 → λj λj+1 , and a new reference state we also find that Q and Q † with | ψ∞ ⟩ form order 2 RS-GAs, creating hole and particle pairs, respectively.

Scarring in a nonintegrable ladder model
For the hole-paired states to be true quantum scars in the ladder model, it is not enough that they are generated from a restricted spectrum generating algebra.A key hallmark of quantum scarring is the existence of eigenstates of a nonintegrable system which are highly anomalous in some observable quantities when compared to other eigenstates of similar energy and thus violate strong eigenstate thermalization hypothesis (ETH), which one would otherwise naively expect to hold [55].One common quantity to compare is the entanglement entropy (EE) across some biparition of the quantum system.Generically, eigenstates of nonintegrable systems have volume law EE across a given bipartition; in a 1D system we thus expect S ∼ L for a generic eigenstate.
In the ladder model, we consider the bipartition defined by cutting the system lengthwise at site L/2 (letting the partitions α, β be the left, right halves of the chain, respectively).It is clear that the reference states |ψ ∞ ⟩, | ψ∞ ⟩ have zero EE across the bipartition, as they are tensor products of a Bell state on each dimer site.Furthermore, as we show rigorously below, the EE of hole-paired states follows an area law in the limit L → ∞: S α ∼ L 0 , thus these eigenstates are quantum scars.
To demonstrate the anomolously low EE of the holepaired states, we perform numerical exact diagonalization of Eq. (G1) for a length L = 8 chain (i.e., 16 total spins) using QuSpin [56].For two fixed magnetization sectors, M = −2, −4, we plot the EE of each eigenstate across the central bond vs. eigenstate energy in Fig. 11.We find that, as is typical of nonintegrable systems, the typical energy eigenstate near the center of the band has high EE (which we expect to be ∼ L).The hole-paired states, identified by stars in the figure, are states in the middle of the eigenspectrum and clearly have anomalously low EE.

Area law entanglement entropy of hole-paired states
For completeness, we rigorously show that the entanglement entropy of the hole-paired states across the central cut at L/2 of the length L chain follows an area law, scaling as S(ρ α,n ) ∼ L 0 , where ρα,n = tr Here we denote the partitions by α, β; the α partition contains the first L/2 qubits of each chain A and B. First, we split the Q operator into three terms: where the first term acts only on the α partition (sites 1 through L/2) and the second on the β partition and the last term QL/2 ∝ τL/2 τL/2+1 places a hole pair across the central cut.Now using the fact that the Q states are orthogonal, specifically that this relation holds term-by-term, and that Q2 L/2 = 0, we can write a generic Q state as Because Qα/β only act on their respective halves of the chain, the terms Similarly for the ∝ QL/2 terms, the only entangling operator is QL/2 , with all powers of Qα/β only acting on their halves of the chain.Thus we can efficiently perform the partial trace over the β partition of the density matrix.In particular, terms with unique powers of Qβ and QL/2 are orthogonal, thus the partial trace over the β partition is Here |•ψ ∞ ⟩ β is the β partition of QL/2 |ψ ∞ ⟩, which has a hole on site L/2 + 1.
For a given fixed n, when we take the limit L → ∞ we have || Qn ) n , as each half of the chain has ≈ L/2 places where each hole-pair can be placed.
Thus tr[ρ α,n ] ∼ (L/2) n 2n n + O(L n−1 ).Denoting the probabilities of terms where a hole-pair does not span the central cut by p k , and denoting the probabilities of terms with a hole-pair spanning the cut by pk , we find We can neglect the pk for L → ∞ (keeping n fixed), so the entanglement entropy of a hole-paired state is arealaw: Note that this argument only holds for a 2n-hole state when taking L → ∞ with n fixed.
Appendix H: Hole-pair correlation function Here we derive the correlation function Eq. ( 28) and discuss its nonstandard normalization.It is natural to measure hole correlations using the dimer operators τj defined in Eq. (C1) as τ z j = [τ † j , τj ] has the matrix elements ⟨• j |τ z j |• j ⟩ = −2 and ⟨• j |τ z j |• j ⟩ = 0 (and zero offdiagonal matrix elements) within the steady state subspace {|• j ⟩, |• j ⟩}.Thus, we define the hole-pair correlation function as which is a standard connected correlation function normalized by the onsite fluctuations.To relate this to measurable quantities in the single A chain alone, we first note that τ z j = σz A,j + σz B,j , and that within the steady state subspace, ⟨• j |σ z s,j |• j ⟩ = −1 and ⟨• j |σ z s,j |• j ⟩ = 0, for s ∈ {A, B}.Thus, when evaluating the correlation functions for an arbitrary state within the steady state subspace, we have ⟨σ z A,j ⟩ = ⟨σ z B,j ⟩, and for j ̸ = k, it is readily shown that all ⟨σ z s,j σz s ′ ,k ⟩ are equivalent for any s, s ′ ∈ {A, B}.Therefore, the connected correlation function for j ̸ = k can be expressed in A-chain observables as For j = k, the expectation values ⟨σ z s,j σz s ′ ,j ⟩ are not all equivalent.The connected correlation functions can be reduced to As a final step, we find that within the steady state subspace, the nonzero matrix elements of σz Thus when restricted to this subspace, the expression 2 + 2⟨σ z A,j σz B,j ⟩ has the same matrix elements as −4⟨σ z A,j ⟩.While these are not equivalent as operators, we thus have the expectation value equivalence in the steady state subspace.Therefore, the correlation function defined above in terms of τ z j is equivalent to Eq. ( 28) when evaluated on the steady state (and any state in the steady state subspace, generally).

ζ-independent upper bound
We first consider the limit |ζ| → 0 in which odd-length chains have CDWs.We estimate the particle density of even and odd length chains by expanding |ψ⟩ to second order in | Ω|: For even chains, the particle numbers of the two terms are 2 and 0, respectively, while for odd chains the particle numbers are 3 and 1.Thus the total chain particle numbers are (using subscripts e and o for even-and oddlength chains, respectively): The onset of the CDW scale occurs when these particle numbers start to diverge from each other -⟨ N ⟩ e vanishes with decreasing Ω whereas ⟨ N ⟩ o saturates to 1 particle.We thus need to estimate the ratio of state norms.
We do not need to estimate the absolute magnitude of the state norms, but only their ratios.Ignoring any boundary effects and assuming that all hole configurations of each state are equally likely, the state norm of Qk |ψ ∞ ⟩ can be given as A(k) × (# configs of k hole−pairs) where A(k) is an amplitude weighting function.Dropping the floor notation for simplicity, for the nearly empty chains, each configuration of QN/2−1 |ψ ∞ ⟩ produces up to one configuration of QN/2 |ψ ∞ ⟩ (neglecting the O(1/N ) rare configurations of three adjacent particles in odd chains) with a permutation factor N/2. Thus, A(N/2) = 1 N (N/2)A(N/2 − 1), where the extra 1/N comes from the normalization of Q (cf.Eq. ( 22)).
For even chains, there are (N/2)(N/2−1)/2 configurations of N/2−1 hole-pairs (i.e., 2 particles), of which N/2 produce the vacuum state (i.e., N/2 hole-pairs).Thus, the state norm ratio for even chains is (J5) From these state norms, we immediately obtain the particle numbers which are valid when we assume the drive strength satisfies Thus we have the upper bound on | Ω| below which CWDs emerge in odd length chains in the |ζ| = 0 limit.This analysis is essentially identical in the |ζ| → ∞ limit via the same reasoning used to explain the emergence of single particle states in even-length chains.Because site 1 always has a hole, we may basically neglect it and consider chains of length N −1 in an effective |ζ| = 0 limit, thus obtaining exactly the same CDW emergence scale.The only difference here is that even chains have CDWs and odd chains do not.

ζ-dependent lower bound
With the upper bound of the CDW regime established, we now turn to the |ζ|-dependent lower bound.This bound is characterized by the point at which the particle number of chains with CDWs starts to fall below 1.Note that the CDWs persist below the lower bound, but with a particle number that vanishes as ∼ | Ω| 2 , i.e., there is less than one particle in the chain on average.
Here we focus on the case of odd chains for |ζ| ≪ 1 -as above the even chain in the |ζ| ≫ 1 regime follows analogously.We expand the odd chain state to first order in | Ω| as These two terms are the CDW and vacuum, respective.
Note that only one configuration of the CDW can be acted upon with τ1 to produce vacuum, hence we immediately obtain the state norm ratio || QN/2 |ψ ∞ ⟩|| 2 /|| QN/2 τ1 |ψ ∞ ⟩|| 2 = N/2.Thus, computing the particle number yields from which we immediately derive the lower bound on the single particle CDW for odd chains Note that for the even chain, we have the same condition but with the inverse |ζ| −1 .Finally we note that if |ζ| ≈ 1, the CDWs of either parity chain will not be particularly pronounced because the vacuum will be reached before the 2 and 3-particle states are fully suppressed.
Appendix K: Notes on experimental realization in circuit QED Here, we briefly describe the necessary components and setup to realize on a circuit QED platform both the twochain model (for remote entanglement stabilization) and the single chain model.For the remote entanglement scheme, a detailed study of waveguide and qubit losses will be given in a complementary work [49].

Remote entanglement realization
We consider two chains of N qubits each, with nearest neighbor capacitive coupling (one may use tunable couplings [57]), and Rabi drives applied to qubit 1 of each chain: Ĥ = (K1) Here, ω q s,j are the qubit frequencies, J is the hopping rate, Ω is the Rabi drive strength, and ω 0 is the common Rabi frequency.We take qubits 2 through N of each chain to be resonant with each other and with the Rabi drives, thus ω q s,j>1 = ω 0 .If the driven qubits of each chain are not resonant with the other qubits, they must be oppositely detuned from ω 0 : ∆ ≡ |ω q s,1 − ω 0 |.Note that the scheme is perfectly robust to any disorder in the hopping rates along the two chains so we can replace J with J j in Eq. (K1).Moving to a common rotating frame at ω 0 and making a rotating wave approximation, we arrive at the desired Hamiltonian  4) and ( 5)).Now it remains to engineer the collective dissipation.In a remote circuit QED realization, the collective dissipation is engineered using a waveguide that links the two remote chains together.In the main text we consider both bidirectional and unidirectional waveguides.The driven qubits are coupled to the waveguide with rate γ.To realize the collective dissipation using a bidirectional waveguide, the qubits must be properly positioned along the waveguide a distance ∆x = nλ 0 /2 apart, i.e., an integer number of half-wavelengths of the drive frequency ω 0 [7].Precise spacing control of qubits along a waveguide has been demonstrated in state-of-the-art waveguide QED experiments [58,59].A unidirectional waveguide can be constructed using microwave circulators that couple the output fields of the qubits to only one propagating direction of the waveguide.Here, there is no spacing requirement to get collective loss.Moreover, the nonreciprocity of the waveguide automatically induces the dissipative exchange Hamiltonian Ĥdiss (cf.Eq. ( 6)) [41].Combined with the Hamiltonian derived above, we arrive at the desired master equation, Eq. (1).

Qubit dephasing
Any experimental realization of a scheme must contend with unwanted dissipation, and the remote entanglement stabilization scheme is no exception.We consider a superconducting circuit realization using transmon qubits.The dominant source of dissipation affecting this scheme is unwanted qubit dephasing, as dephasing quickly degrades the coherence of the entangled states.Other sources of dissipation, namely qubit relaxation and losses in the waveguide, will be discussed in more detail in a complementary work [49].Here we provide an estimate for how much transmon dephasing can be tolerated using realistic driving, hopping and dissipation rates.A heuristic for determining how much qubit dephasing can be tolerated is to compare the typical qubit coherence time T 2 to the stabilization time τ rel of the ideal system (i.e., the relaxation time of the slowest Liouvillian eigenmode).When τ rel ≪ T 2 , the qubit dephasing will have relatively small effects on the steady state, but when τ rel ≈ T 2 , the performance of the scheme is rapidly degraded.
Based on previous waveguide QED experiments [59][60][61], engineered collective dissipation rates have been demonstrated ranging from γ/2π = 20 MHz to γ/2π = 100 MHz.Here we take a rather conservative set of parameters γ = Ω = J = 2π × 3 MHz.In the dephasingfree system with 3 + 3 qubits, these parameters lead to 90% fidelity to the ideal tensor product of Bell states |ST S⟩.We find the stabilization rate numerically to be τ rel ≈ 156/γ = 9 µs.This is considerably faster than the typical qubit coherence time, which have been shown to be typically at least T 2 ∼ 20 µs [59] up to T 2 ∼ 1 ms [62].Moreover, as discussed in Sec.V, utilizing the natural third level of transmon accelerates the dissipation process by orders of magnitude.For a 3 + 3 system with optimal qutrit coupling, here instead we aim for higher Bell state fidelity of 99% with γ = 2π × 3 MHz, and + h.c. .

FIG. 5 .
FIG. 5. Flavor-change hopping forces hole-pairing in zero-energy eigenstates.(a) ĤXX (cf.Eq. (5)) swaps a particle and a hole in the dimer chain and changes the particle flavor in the process.(b) Due to the flavor change in hopping, a single delocalized hole is not an eigenstate of the chain because the final states of a hole hopping to the left and hole hopping to the right are distinguishable: the |■⟩ either ends up to the right of the hole or to the left of the hole, respectively.(c) Adjacent paired holes can form zero-energy eigenstates via destructive interference because the |■⟩ always ends up sandwiched between two holes.

6 C 3 |j − k| = 4 FIG. 6 .
FIG.6.Density correlations on adjacent sites are a direct observable consequences of hole pairing.The hole density correlation function Czz(j, k) (cf.Eq. (28)) is averaged over averaged over the entire N = 40 chain, while holding |k−j| fixed.The averaged Czz(j, k) is plotted vs. drive strength | Ω|.The nearest neighbor (|k − j| = 1) correlations saturate to 0.5 with increasing | Ω|, indicating that as the filled state |ψ∞⟩ is approached, deviations from |ψ∞⟩ in the bulk are due to holes paired on adjacent sites.For this plot we hold |Γ| = J fixed.

FIG. 8 .
FIG. 8.The emergence of single particle "CDWs" in finite length spin chains.(a) Local particle density ⟨nj⟩ is plotted for each site of 40 and 41 site chains, scaled by length N .For the given parameters, Ω = 0.1 and |ζ| = 0.005, and taking uniform Jj, the odd-length chain has a single particle delocalized across all odd sites, hence the average particle density is n = 1/2N ≈ 0.012.The even length chains do not have charge density waves, and have much smaller average density n = N | Ω| 4 /8 ≈ 0.0005.The modulation of local density across the chain is due to the highly-correlated twoparticle state.(b) Particle density n is plotted vs. drive strength | Ω| for two odd-length chains (N = 11 and N = 10001) in the odd-length CDW driving strength regime (cf.Eq.(38)) for two different |ζ| 2 = 10 −4 and |ζ| 2 = 10 −2 .In the limit |ζ| → 0 (dashed curves), the single particle CDW persists for any arbitrarily small | Ω| > 0. In both plots, the particle density is n = 1

FIG. 9 .
FIG. 9. Speeding up entanglement stabilization time of the nonreciprocal double chain using a qutrit.(a)The nonreciprocal double chain is modified by replacing the downstream B1 qubit with a qutrit and engineering the nonreciprocal coupling to include the 2-1 transition of the qutrit.(b) The numerically computed relaxation time γτ rel (red) and the infidelity of the steady state to the maximally entangled state |ψ∞⟩ (cf.Eq. (12)), 1 − ⟨ψ∞|ρ|ψ∞⟩ (black), are shown as functions of the hopping rate J/γ for N = 3, Ω/γ = 10, and ∆ = 0.For the qubit-qutrit scheme, the relaxation time is optimized over η.We also plot the relaxation time for the single chain system for comparison.For a fixed state fidelity of 0.999 (achieved at J = 7γ, dashed line), the relaxation times are γτ 2qb = 2.1 × 10 4 , γτqutrit = 120 and γτsc = 14, for the qubit-qubit, qubit-qutrit, and single chain, respectively.

FIG. 11 .
FIG. 11.Hole-paired states as scars in the double XX chain ladder model.The entanglement entropy of the ladder model eigenstates (cf.Eq. (G1)) of the α partition (i.e., sites 1 to L/2) is plotted vs. the eigenstate energy.Here we consider a length L = 8 ladder (16 total spins) and plot the eigenstates in the M = −2 magnetization sector (i.e., 2-hole sector) and the M = −4 (4-hole) sector.The holepaired states are denoted with stars; in each magnetization sector there are two states corresponding to hole-pairs created in each reference state |ψ∞⟩ and | ψ∞⟩.For both plots the parameters are µ/J = 0.2, g/J = 0.5, and δ/g = 0.3.