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

Andrew Lingenfelter<EMAIL_ADDRESS>Pritzker School of Molecular
Engineering, University of Chicago, Chicago, IL 60637, USA  Department of
Physics, University of Chicago, Chicago, IL 60637, USA  Mingxing Yao  Pritzker
School of Molecular Engineering, University of Chicago, Chicago, IL 60637, USA
Andrew Pocklington  Pritzker School of Molecular Engineering, University of
Chicago, Chicago, IL 60637, USA  Department of Physics, University of Chicago,
Chicago, IL 60637, USA  Yu-Xin Wang (王语馨)  Pritzker School of Molecular
Engineering, University of Chicago, Chicago, IL 60637, USA  Abdullah Irfan
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL
61801, USA  Wolfgang Pfaff  Department of Physics, University of Illinois at
Urbana-Champaign, Urbana, IL 61801, USA  Aashish A. Clerk  Pritzker School of
Molecular Engineering, University of Chicago, Chicago, IL 60637, USA

###### Abstract

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.

##  I Introduction

The recent intense interest in driven-dissipative quantum systems has many
distinct motivations. One key goal is to understand how tailored dissipative
processes [1, 2] could be used to stabilize entangled quantum states in both
the few and many-body regime, with possible applications to quantum
information processing (see e.g. [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22]). A second, seemingly distinct line of work seeks
to understand the unique properties of non-equilibrium steady states (NESS)
that arise in driven quantum spin chains from the interplay of driving,
lattice dynamics and dissipation [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35]. Here, certain exact solutions have been especially valuable [23, 24,
25, 26, 27].

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

Figure 1:  (a) Two XX-coupled $N$ qubit chains, with boundary collective loss
(rate $\gamma$, mediated by a waveguide) and Rabi drives (strength $\Omega$,
detuning $\pm\Delta$). (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 $\tilde{\Omega}$ (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
$\Omega<\gamma$.

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 real-space 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 double-chain system has potential
advantages over other proposals [20, 16, 21, 19, 17, 18] 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 [5, 4, 6,
7, 8]. For large Rabi drive amplitude $\Omega$, 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, 19, 17, 18], the underlying mechanism had not been
elucidated. For finite $\Omega$, the two qubit systems need $\Omega>\gamma$
for significant steady state entanglement, where $\gamma$ is the waveguide-
induced dissipation. We find that this is no longer the case for larger
systems: for $N\geq 2$, strong steady state entanglement only requires the
potentially much weaker condition $\Omega>\sqrt{J\gamma}$, 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 $\Omega$. 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-Bell-pair entanglement
stabilization.

##  II XX-coupled qubit chains with boundary dissipation and driving

###  II.1 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

$\displaystyle\partial_{t}\hat{\rho}$
$\displaystyle=-i[\hat{H},\hat{\rho}]+\gamma\mathcal{D}[\hat{c}]\hat{\rho}$
(1) $\displaystyle\hat{H}$ $\displaystyle=\hat{H}_{\rm drive}+\hat{H}_{\rm
XX}+\hat{H}_{\rm diss},$ (2) $\displaystyle\hat{c}$
$\displaystyle=\hat{\sigma}_{A,1}^{-}+\hat{\sigma}_{B,1}^{-},$ (3)

where the Hamiltonian terms are given by

$\displaystyle\hat{H}_{\rm drive}$
$\displaystyle=\frac{\Omega}{2}\left(\hat{\sigma}_{A,1}^{x}+\hat{\sigma}_{B,1}^{x}\right)+\frac{\Delta}{2}\left(\hat{\sigma}_{A,1}^{z}-\hat{\sigma}_{B,1}^{z}\right),$
(4) $\displaystyle\hat{H}_{\rm XX}$
$\displaystyle=\frac{1}{2}\sum_{j=1}^{N-1}\sum_{s=A,B}J_{j}\left(\hat{\sigma}_{s,j}^{+}\hat{\sigma}_{s,j+1}^{-}+{\rm
h.c.}\right).$ (5) $\displaystyle\hat{H}_{\rm diss}$
$\displaystyle=\frac{1}{2}i\nu\gamma\left(\hat{\sigma}_{A,1}^{+}\hat{\sigma}_{B,1}^{-}-{\rm
h.c.}\right).$ (6)

The Lindblad dissipator
$\mathcal{D}[\hat{c}]\,\cdot=\hat{c}\,\cdot\,\hat{c}^{\dagger}-\\{\hat{c}^{\dagger}\hat{c},\,\cdot\,\\}/2$
describes collective loss on the site 1 qubits (with
$\hat{\sigma}^{-}=|0\rangle\langle 1|$, $\hat{\sigma}^{z}=|1\rangle\langle
1|-|0\rangle\langle 0|$).

Consider first the driving of our system. Qubits $A1$ and $B1$ are Rabi-driven
at the same frequency and amplitude $\Omega$. We however take qubit $A1$
($B1$) to be detuned by $+\Delta$ ($-\Delta$) from the drive frequency.
Treating these drives within the rotating wave approximation, we obtain the
rotating frame Hamiltonian $\hat{H}_{\rm 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: $\hat{H}_{\rm 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 $\gamma$) 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
$\hat{H}_{\rm diss}$ given by Eq. (6) [41], where $\nu$ is a directionality
factor, $-1\leq\nu\leq 1$. When $\nu=0$ the waveguide is perfectly
bidirectional and when $\nu=+1$ ($\nu=-1$) the waveguide is perfectly
unidirectional with system $A$ (system $B$) upstream of $B$ ($A$).

###  II.2 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 $\hat{c}_{\rm
pair}=u\hat{\sigma}_{A,1}^{-}+v\hat{\sigma}_{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 qubit-waveguide couplings [7, 16]; this is also
challenging in many setups.

A simpler approach for stabilizing remote two qubit entanglement (using only
Rabi drives and passive loss) is provided by the $N=1$ version of Eq. (1).
This both unifies and generalizes the entanglement stabilization schemes
studied in Refs. [5, 4, 6, 7]. The bidirectional waveguide case ($\nu=0$)
yields the schemes of Refs. [6, 7] and the unidirectional waveguide case
($\nu=\pm 1$) yields the scheme introduced in Ref. [4].

The general $N=1$ system has a pure steady state
$\hat{\rho}_{1}(t\to\infty)=|\psi_{1}\rangle\langle\psi_{1}|$ given by (up to
normalization)

$\displaystyle|\psi_{1}\rangle$
$\displaystyle=\sqrt{2}\Gamma|00\rangle+\Omega|S\rangle,$ (7)
$\displaystyle\Gamma$ $\displaystyle\equiv\Delta-\frac{1}{2}i\nu\gamma,$ (8)

where $\Gamma\in\mathbb{C}$ is a generalized complex detuning. Here we
introduce the singlet and triplet entangled states

$\displaystyle|S\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(|01\rangle-|10\rangle),\quad|T\rangle=\frac{1}{\sqrt{2}}(|01\rangle+|10\rangle),$
(9)

in addition to the unentangled vacuum $|00\rangle$ and the “doublon”
$|11\rangle$. Eq. (7) is the _unique_ two-qubit steady state for any
$|\Omega/\Gamma|<\infty$, and as $|\Omega/\Gamma|\to\infty$, it approaches a
perfect Bell state $|\psi_{1}\rangle\to|S\rangle$. Note, however, that in this
limit the dissipative gap closes and a second _impure_ steady state emerges
[9, 8, 7].

Going forward, one key goal will be to extend this entanglement stabilization
to the case where each chain has $N\gg 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\gg 1$ chain is complicated and crucially
depends on the interplay between the boundary driving and lattice dynamics.
However, this interplay also gives the $N\gg 1$ chain advantages over
squeezing based schemes by enabling large amounts of entanglement to be
stabilized without requiring $\Omega\gg\gamma$. To describe the stabilized
state for $N\gg 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}\rangle=|0_{A,j}1_{B,j}\rangle$ to denote the state of the qubit
pair on dimer site $j$. For example, $|S_{j}\rangle$ denotes a Bell pair
spanning qubits $Aj$ and $Bj$.

###  II.3 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

$\displaystyle\partial_{t}\hat{\rho}_{A}$
$\displaystyle=-i[\hat{H}_{A},\hat{\rho}_{A}]+\gamma\mathcal{D}[\hat{\sigma}_{A,1}^{-}]\hat{\rho}_{A},$
(10) $\displaystyle\hat{H}_{A}$
$\displaystyle=\frac{\Omega}{2}\hat{\sigma}_{A,1}^{x}+\frac{\Delta}{2}\hat{\sigma}_{A,1}^{z}+\frac{1}{2}\sum_{j}J_{j}\big{(}\hat{\sigma}_{A,j}^{+}\hat{\sigma}_{A,j+1}^{-}+{\rm
h.c.}\big{)}.$

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\rangle$), 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 $\hat{H}_{A}$ alone. As we
show in App. A, upon making a Jordan-Wigner transformation, the Rabi drive in
$\hat{H}_{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: $\partial_{t}\hat{\rho}=-i[\hat{H}_{{\rm
fermi},\eta},\hat{\rho}]+\gamma\mathcal{D}[(-1)^{\hat{\eta}^{\dagger}\hat{\eta}}\hat{c}_{1}]\hat{\rho}$,
where $\hat{H}_{{\rm fermi},\eta}$ is a quadratic fermionic Hamiltonian,
$\hat{c}_{1}$ is the fermion lowering operator on the dissipative site, and
$\hat{\eta}$ 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 $|\psi\rangle$ 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.
$\nu=1$), simply tracing out the $B$ chain yields the steady state of the
single-chain system in Eq. (10).

$\displaystyle\hat{\rho}_{A}={\rm tr}_{B}|\psi\rangle\langle\psi|,$ (11)

This corresponds to a new many-body application of the coherent quantum
absorber technique introduced in Ref. [4] and extended in Refs. [37, 36, 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 two-time 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 $\Omega$

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
$\Omega$ 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, 19,
17, 18]). For finite drives $\Omega$, 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}\rangle$ state.

Figure 2:  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 $J_{B}=J_{A}+\delta$
(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:
$\hat{H}_{\Delta}=(\delta/2)(\hat{\sigma}_{A,2}^{z}\pm\hat{\sigma}_{B,2}^{z})$.
The purity of the numerically computed steady state $\hat{\rho}_{\rm
ss}=\hat{\rho}(t\to\infty)$ is shown in each case vs. the perturbation
strength $\delta$. The steady state is only pure when $\delta=0$. We use
unperturbed system parameters $\Omega=0.5\gamma$, $\Delta=0.1\gamma$, and
$J=0.25\gamma$.

###  III.1 $\Omega\to\infty$ 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 $\Omega$:

$\displaystyle|\psi\rangle\rightarrow|\psi_{\infty}\rangle$
$\displaystyle=|S_{1}T_{2}S_{3}T_{4}S_{5}T_{6}\cdots(S/T)_{N}\rangle.$ (12)

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 phase of these pairs
alternates from $|S\rangle$ or $|T\rangle$ as one moves down the chain as
indicated. Up to this local phase variation, the state is translationally
invariant. This result is in fact the consequence of a much 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, 19, 17, 18]. 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 $\hat{\mathcal{L}}$
acting on a 2-qubit system that stabilizes an arbitrary state $|\psi\rangle$
with fixed excitation number parity, i.e.
$|\psi\rangle=a|00\rangle+b|11\rangle$ or
$|\psi\rangle=a|01\rangle+b|10\rangle$. 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
$\Omega\rightarrow\infty$, 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). For convenience, we use a different gauge choice
for the $B$ chain (i.e. $|1\rangle\rightarrow-|1\rangle$ for the $B$ qubits),
such that the XX couplings now have opposite signs in the $B$ chain versus the
$A$ chain. Specifically, if $\hat{\mathcal{L}}$ stabilizes
$|\psi_{1}\rangle=a|(00)_{1}\rangle+b|(11)_{1}\rangle$, then we are interested
in the Hamiltonian

$\displaystyle\hat{H}_{\mathrm{XX},2}$
$\displaystyle=J\sum_{s=A,B}(-1)^{s}\left[\hat{\sigma}^{x}_{s,1}\hat{\sigma}^{x}_{s,2}+\hat{\sigma}^{y}_{s,1}\hat{\sigma}^{y}_{s,2}\right].$
(13)

This Hamiltonian now has an extremely convenient feature: no matter what the
parameters $a,b$, the “replicated”, dimerized state
$|\psi\rangle=|\psi_{1}\psi_{2}\rangle$ is a zero eigenstate, as can be
confirmed by a simple direct computation

$\displaystyle\hat{H}_{\mathrm{XX},2}|\psi_{1}\psi_{2}\rangle$
$\displaystyle=0.$ (14)

This tells us that $|\psi\rangle=|\psi\rangle_{1}\otimes|\psi\rangle_{2}$ is a
steady state of the 4-qubit dynamics defined by

$\displaystyle\partial_{t}\hat{\rho}$
$\displaystyle=\left(\hat{\mathcal{L}}\otimes\mathds{1}\right)\hat{\rho}-i[\hat{H}_{\mathrm{XX},2},\hat{\rho}].$
(15)

Hence, without changing the dissipative stabilization mechanism at all, we can
use the passive Hamiltonian interaction to propagate entanglement to a second
qubit pair. Even more strikingly, we can now add a third pair of qubits (again
via mirrored XX couplings). The same arguments tell us that the steady state
will be a replicated entangled state, i.e. a product of cross-chain dimers,
where each entangled dimer is in the state $|\psi\rangle$. Iterating this
argument, we can show that for arbitrary $N$, if we let

$\displaystyle\hat{H}_{\mathrm{XX},N}$
$\displaystyle=\sum_{i=1}^{N-1}\sum_{s=A,B}J_{i}(-1)^{s}\left[\hat{\sigma}^{x}_{s,i}\hat{\sigma}^{x}_{s,i+1}+\hat{\sigma}^{y}_{s,i}\hat{\sigma}^{y}_{s,i+1}\right],$
(16) $\displaystyle|\psi\rangle$
$\displaystyle=\bigotimes_{i=1}^{N}\bigg{[}a|(00)_{i}\rangle+b|(11)_{i}\rangle\bigg{]},$
(17)

then $|\psi\rangle$ is still a steady state of the dynamics

$\displaystyle\left(\hat{\mathcal{L}}\otimes\mathds{1}^{\otimes
N-1}\right)|\psi\rangle\langle\psi|-i\left[\hat{H}_{\mathrm{XX},N},|\psi\rangle\langle\psi|\right]=0.$
(18)

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\to(-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}\rangle$) 111The argument here does not prove uniqueness of the
replicated steady state. In many cases, including that studied here, the lack
of additional symmetries precludes additional steady states.. The change of
gauge thus explains why Eq. (12) is a product of staggered $|S\rangle$ and
$|T\rangle$ instead of uniform $|S\rangle$ states.

Figure 3:  Entanglement replication via passive exchange couplings. Suppose
there is some dissipative dynamics $\hat{\mathcal{L}}$ that stabilizes a two
qubit entangled state $|\psi\rangle$ on the first pair of qubits in two
exchange-coupled chains. If that state has definite parity, then the product
state $|\psi_{1}\psi_{2}\rangle$ is a zero energy eigenstate of passive
exchange couplings $\hat{H}_{\rm 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
$|\psi_{1}\psi_{2}\cdots\psi_{N}\rangle$, replicating the entanglement on site
1 down a pair of arbitrarily long chains.

The fact that the Hamiltonian perfectly replicates fixed parity states can be
understood intuitively from the fact that it commutes with $\hat{S}^{2}_{s}$
the total spin operator and
$\hat{S}^{z}_{s}=\hat{\sigma}_{s,1}^{z}+\hat{\sigma}_{s,2}^{z}$ 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.

###  III.2 $\Omega<\infty$: pure steady-state condensate of paired holes

We now turn to the more general (and experimentally relevant) case of a non-
infinite drive amplitude $\Omega$. 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\geq 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}\rangle|(T)_{j+1}\rangle\to|(00)_{j}\rangle|(00)_{j+1}\rangle$.

Certain features of our state can be understood intuitively. For finite
$\Omega$, 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 $\hat{H}_{\rm 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” $|\bullet_{j}\rangle$ and
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{j}\rangle$
and “hole” $|\circ_{j}\rangle$ states via

$\displaystyle|\circ_{j}\rangle$ $\displaystyle\equiv|(00)_{j}\rangle,$ (19)
$\displaystyle|\bullet_{j}\rangle$
$\displaystyle\equiv\frac{1}{\sqrt{2}}\hat{\tau}_{j}^{\dagger}|\circ_{j}\rangle=\begin{cases}|S_{j}\rangle&j~{}{\rm
odd}\\\ |T_{j}\rangle&j~{}{\rm even}\end{cases},$ (20)
$\displaystyle|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{j}\rangle$
$\displaystyle\equiv\frac{1}{\sqrt{2}}\hat{\lambda}_{j}^{\dagger}|\circ_{j}\rangle=\begin{cases}|T_{j}\rangle&j~{}{\rm
odd}\\\ |S_{j}\rangle&j~{}{\rm even}\end{cases},$ (21)

where we implicitly define the dimer ladder operators
$\hat{\tau}_{j}^{\dagger},\hat{\tau}_{j}$ that create and destroy the dimer
particles $|\bullet_{j}\rangle=|(S/T)_{j}\rangle$ (for $j$ odd/even) when
acting on $|\circ_{j}\rangle$ or $|\bullet_{j}\rangle$, respectively.
Similarly $\hat{\lambda}_{j}^{\dagger},\hat{\lambda}_{j}$ create and destroy
the particles
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{j}\rangle=|(T/S)_{j}\rangle$
(for $j$ odd/even) when acting on $|\circ_{j}\rangle$ or
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{j}\rangle$,
respectively. Note that there is a hard-core constraint that prevents a
$|{\bullet}\rangle$ and a
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
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
$|\psi_{\infty}\rangle=|\bullet\bullet\bullet\cdots\rangle$.

Using the dimer particle representation defined in Eq. (20), we introduce an
operator that creates a delocalized hole-pair:

$\displaystyle\hat{Q}=\frac{1}{2\bar{J}\sqrt{N}}\sum_{j=1}^{N-1}J_{j}(-1)^{j}\hat{\tau}_{j}\hat{\tau}_{j+1}.$
(22)

Here

$\displaystyle\bar{J}=\sqrt{\frac{1}{N-1}\sum_{j}J_{j}^{2}}$ (23)

is the RMS hopping rate. Acting on the reference “filled” state
$|\psi_{\infty}\rangle$, the operator $\hat{Q}$ creates a superposition state
where each term corresponds to a pair of adjacent holes at a different
location 222For convenience, we normalize the operator by $1/\bar{J}\sqrt{N}$
so it produces an approximately normalized state when acting on
$|\psi_{\infty}\rangle$ and other nearly-filled states
(${||}\hat{Q}|\psi_{\infty}\rangle{||}^{2}\approx{||}\hat{Q}^{2}|\psi_{\infty}\rangle{||}^{2}\approx
1$)..

At a heuristic level, the phases in $\hat{Q}$ will give each hole pair a net
momentum, allowing them to become zero-energy eigenstates of the “kinetic
energy” $\hat{H}_{\rm XX}$. More formally, we show in App. D that $\hat{Q}$
commutes with $\hat{H}_{\rm XX}$. Its action on an $\hat{H}_{\rm XX}$
eigenstate will thus produce another eigenstate with the same energy. In
particular:

$\displaystyle\hat{H}_{\rm
XX}(\hat{Q}|\psi_{\infty}\rangle)=\hat{Q}\hat{H}_{\rm
XX}|\psi_{\infty}\rangle=0.$ (24)

An entire tower of zero-energy $\hat{H}_{\rm XX}$ eigenstates can thus be
generated by repeated application of $\hat{Q}$ to $|\psi_{\infty}\rangle$,
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 delocalized particle $|\bullet_{j}\rangle$ (if $N$ is odd).

Figure 4:  First few terms of the steady state. The expansion of the steady
state in the number of holes added to $|\psi_{\infty}\rangle$ (cf. Eq. (12)).
The expansion coefficients $a_{n}\sim 1/\tilde{\Omega}^{n}$ are the weights of
the $n$th-hole components of the state; they can be found analytically from
Eq. (25). Here we take uniform $J_{j}=\bar{J}$. For nonuniform $J_{j}$ the
components of each $a_{n}$ term are weighted by factors of $J_{j}/\bar{J}$.

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 $\hat{Q}^{m}|\psi_{\infty}\rangle$ provide
us with a large class of states 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 $|\psi_{Q}\rangle$ given by the “pair
condensate” state

$\displaystyle|\psi_{Q}\rangle$
$\displaystyle=\left(1+\frac{\Gamma}{\Omega}\hat{\tau}_{1}\right)\exp\left[\frac{\sqrt{N}}{\tilde{\Omega}^{2}}\hat{Q}\right]|\psi_{\infty}\rangle$
(25)

where $\hat{\tau}_{1}$ is the dimer lowering operator that removes the
particle on site 1 (cf. Eq. (20)) and the reference state
$|\psi_{\infty}\rangle$ is given by Eq. (12). The dimensionless drive strength
$\tilde{\Omega}$ appearing in the exponential is

$\displaystyle\tilde{\Omega}$
$\displaystyle\equiv\frac{\Omega}{\sqrt{\Gamma\bar{J}}},$ (26)

with $\Gamma$ given by Eq. (8) and $\bar{J}$ by Eq. (23). The first few terms
of $|\psi_{Q}\rangle$ 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}=\bar{J}$) $a_{0}=1$,
$a_{1}=\Gamma/\Omega=\sqrt{\Gamma/\bar{J}}/\tilde{\Omega}$,
$a_{2}=1/\tilde{\Omega}^{2}$, etc; this is a power series in
$1/\tilde{\Omega}$, with $a_{j}\sim 1/\tilde{\Omega}^{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
$|\psi_{Q}\rangle$ 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
$\langle\hat{n}_{j}\rangle=\frac{1}{2}\langle\hat{\tau}_{j}^{\dagger}\hat{\tau}_{j}\rangle$.
Details are provided in App. F.

##  IV Real-space hole pairing and consequences for the steady state

###  IV.1 Hole pairing as a kinetic constraint

Figure 5:  Flavor-change hopping forces hole-pairing in zero-energy
eigenstates. (a) $\hat{H}_{\rm 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
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
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
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
always ends up sandwiched between two holes.

The hole-pairing in $|\psi_{Q}\rangle$ 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) $\hat{H}_{\rm
XX}$ causes $|{\bullet}\rangle$ and
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
particles to change flavor when they hop, and (ii) the hard-core constraint
forbids a $|{\bullet}\rangle$ and an
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
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 $\hat{H}_{\rm XX}|\\!\bullet\bullet\rangle=0$. If we now try
to form zero kinetic energy states (i.e., zero-energy eigenstates of
$\hat{H}_{\rm 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 flavor-changing hopping, see Fig. 5(b), but
delocalizing a pair is successful (see Fig. 5(c)).

The creation of zero kinetic energy eigenstates of $\hat{H}_{\rm XX}$ via the
hole-pairing operator $\hat{Q}$ (cf. Eq. (22)) is somewhat analogous to
$\eta$-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 $\eta$-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 $\eta$-pairing ($\hat{\eta},\hat{\eta}^{\dagger}$
form a closed representation of SU(2)) is not found in $Q$-pairing: $\hat{Q}$
and $\hat{Q}^{\dagger}$ 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 $|\psi_{\infty}\rangle$ [39].
Furthermore, certain terms can be added to the XX chain to make the model
nonintegrable while preserving this RSGA structure, thus guaranteeing that
$\hat{Q}^{n}|\psi_{\infty}\rangle$ 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 $\hat{Q}$ has the same properties, generating
eigenstates of the 1D chain when acting on a filled state. These hole-paired
eigenstates may be accessible in e.g. the ultracold atoms platform proposed in
Ref. [38].

###  IV.2 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. $|\tilde{\Omega}|\gtrsim 2$. This correlation is directly observable
in the single-chain system of Fig. 1(a) as $Z$ magnetization correlations
between adjacent sites,
$\langle\hat{\sigma}^{z}_{A,j}\hat{\sigma}^{z}_{A,j+1}\rangle\neq 0$. The
correspondence between magnetization and hole density follows from the fact
that the dimer holes $|\circ\rangle$ are polarized but the particles
$|\bullet\rangle$ are depolarized (cf. Eqs. (19) and (20)):

$\displaystyle\langle\bullet_{j}|\hat{\sigma}^{z}_{A,j}|\bullet_{j}\rangle=0;\qquad\langle\circ_{j}|\hat{\sigma}^{z}_{A,j}|\circ_{j}\rangle=-1.$
(27)

Thus, $(-\hat{\sigma}^{z}_{A,j})$ acts as a local hole number operator when
acting on the steady state. We define the $z$-magnetization correlation
function for the $A$ chain as

$\displaystyle
C_{zz}(j,k)=\frac{\langle\hat{\sigma}^{z}_{A,j}\hat{\sigma}^{z}_{A,k}\rangle-\langle\hat{\sigma}^{z}_{A,j}\rangle\langle\hat{\sigma}^{z}_{A,k}\rangle}{\sqrt{\big{(}\langle\hat{\sigma}_{A,j}^{z}\rangle+\langle\hat{\sigma}_{A,j}^{z}\rangle^{2}\big{)}\big{(}\langle\hat{\sigma}_{A,k}^{z}\rangle+\langle\hat{\sigma}_{A,k}^{z}\rangle^{2}\big{)}}},$
(28)

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 averaged over the whole chain. For strong driving, the average
correlations between adjacent sites saturates to $\overline{C_{zz}(j,j\pm
1)}\to 0.5$ because in this regime, a hole on site $j$ is always paired with a
hole on either $j\pm 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\pm 1)\approx 0.5$ for
$\tilde{\Omega}>1$ is a clear indication of hole pairing.

Figure 6:  Density correlations on adjacent sites are a direct observable
consequences of hole pairing. The hole density correlation function
$C_{zz}(j,k)$ (cf. Eq. (28)) is averaged over averaged over the entire $N=40$
chain, while holding $|k-j|$ fixed. The averaged $\overline{C_{zz}(j,k)}$ is
plotted vs. drive strength $|\tilde{\Omega}|$. The nearest neighbor
($|k-j|=1$) correlations saturate to $0.5$ with increasing $|\tilde{\Omega}|$,
indicating that as the filled state $|\psi_{\infty}\rangle$ is approached,
deviations from $|\psi_{\infty}\rangle$ in the bulk are due to holes paired on
adjacent sites. For this plot we hold $|\Gamma|=\bar{J}$ fixed.

###  IV.3 Universal density scaling

There are two dimensionless parameters in our model: the dimensionless drive
amplitude $\tilde{\Omega}$ (c.f. Eq. (26)) and the ratio
$\zeta\equiv\sqrt{\Gamma/\bar{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

$\displaystyle|\psi_{Q}\rangle$
$\displaystyle=\left(1+\frac{\zeta}{\tilde{\Omega}}\hat{\tau}_{1}\right)\exp\left[\frac{\sqrt{N}}{\tilde{\Omega}^{2}}\hat{Q}\right]|\psi_{\infty}\rangle,$
(29)

which suggests that the excitation density in the bulk is controlled only by
$\tilde{\Omega}$. We now show explicitly that the excitation density (i.e. $Z$
magnetization density) is indeed intensive, and scales universally with the
single parameter $|\tilde{\Omega}|$ in the regime $|\tilde{\Omega}|\gtrsim 2$.

Consider first the limit where $|\zeta|\to 0$ while $|\tilde{\Omega}|$ remains
fixed, such that we can ignore the boundary term in Eq. (29). We thus have
$|\psi_{Q}\rangle=e^{\alpha\hat{Q}}|\psi_{\infty}\rangle$, where
$\alpha=\sqrt{N}/|\tilde{\Omega}|^{2}$. This expression mimics a bosonic
coherent state $|\alpha\rangle=e^{\alpha\hat{a}^{\dagger}}|0\rangle$, where
hole pairs are the bosons, $|\psi_{\infty}\rangle$ is the hole-pair vacuum and
$\hat{Q}$ is the hole-pair creation operator. As we show in App. I, this
analogy to bosonic coherent states can be made precise when $N\gg 1$ and the
hole density $\bar{m}$ satisfies

$\displaystyle\bar{m}$
$\displaystyle\equiv\frac{1}{N}\sum_{j}\langle-\hat{\sigma}^{z}_{A,j}\rangle\ll
1.$ (30)

In this low hole density regime, we may neglect the hard-core constraint that
prevents two hole-pairs from occupying the same sites, and can thus accurately
estimate the total hole number as

$\displaystyle M=2|\alpha|^{2}=\frac{2N}{|\tilde{\Omega}|^{4}}\equiv
N\bar{m}.$ (31)

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 $|\tilde{\Omega}|\gtrsim 2$ to approach the infinite
temperature state.

When $|\zeta|>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 $|\tilde{\Omega}|$ but also $|\zeta|$. The hole density $\bar{m}$ thus
obtains a $\mathcal{O}(1/N)$ correction to the intensive universal scaling

$\displaystyle\bar{m}=\frac{2}{|\tilde{\Omega}|^{4}}+\frac{2|\zeta|^{2}}{N(|\tilde{\Omega}|^{2}+2|\zeta|^{2})}.$
(32)

There is thus a transition from universal $\sim|\tilde{\Omega}|^{-4}$ scaling
to an $N$-dependent $\sim|\tilde{\Omega}|^{-2}$ scaling at
$|\tilde{\Omega}|^{2}\approx N/|\zeta|^{2}$, or equivalently
$\Omega/\bar{J}\approx\sqrt{N}$. The universal scaling behavior is exact for
any $|\zeta|$ when the dissipative site density is excluded, as is shown in
Fig. 7(a), and the dissipative site corrections for $|\zeta|>0$ are shown in
Fig. 7(b).

Figure 7:  Universal scaling of steady state $Z$ magnetization (hole density)
The $Z$ magnetization / hole density $\bar{m}$ (cf. Eq. (30)) is shown for
chains of widely varying length as a function of inverse effective drive
strength $1/|\tilde{\Omega}|$. (a) Here, we plot the hole density, excluding
the dissipative site, which shows universal scaling
$\bar{m}=2/|\tilde{\Omega}|^{4}$ for $1/|\tilde{\Omega}|\gtrsim 1$. (b) Here,
we include the dissipative site and show that $\bar{m}$ now has deviations
$\sim 1/N|\tilde{\Omega}|^{2}$ appearing when $|\tilde{\Omega}|^{2}\approx
N/|\zeta|^{2}$. For both plots, $|\zeta|^{2}\equiv|\Gamma/\bar{J}|=0.05$
except where indicated otherwise, and we vary only $\Omega$.

###  IV.4 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 $\tilde{\Omega}$. 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}=\bar{J}$) is given by the CDW form:

$\displaystyle|\Psi_{\rm cdw}\rangle=\frac{1}{\sqrt{\lceil
N/2\rceil}}\sum_{j}(-1)^{j}|\bullet_{N-2j}\rangle,$ (33)

where $\lceil X\rceil$ is the ceiling function, and $|\bullet_{j}\rangle$ 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
$|\bullet_{j}\rangle$ is weighted by an additional factor $(\cdots
J_{j-4}J_{j-2})(J_{j+1}J_{j+3}\cdots)/\bar{J}^{\lceil(N-1)/2\rceil}$ (and the
requisite correction to the normalization).

For the single dissipative spin chain, Fig. 1(b), the non-equilibrium steady
state $\hat{\rho}_{A,{\rm cdw}}={\rm tr}_{B}|\Psi_{\rm
cdw}\rangle\langle\Psi_{\rm cdw}|$ is an equal mixture of vacuum and the
single-particle CDW:

$\displaystyle\hat{\rho}_{A,{\rm cdw}}$
$\displaystyle=\frac{1}{2}\Big{[}|0\rangle\langle 0|+|\Phi_{\rm
cdw}\rangle\langle\Phi_{\rm cdw}|\Big{]},$ (34) $\displaystyle|\Phi_{\rm
cdw}\rangle$ $\displaystyle=\frac{1}{\sqrt{\lceil
N/2\rceil}}\sum_{j}(-1)^{j}\hat{\sigma}_{A,N-2j}^{+}|0\rangle.$ (35)

Despite not being a pure state, the single chain CDW is a very low entropy
state, $S(\hat{\rho}_{A,{\rm cdw}})=\ln 2$, for which any particle density
necessarily arises from the single coherently delocalized excitation. Here,
$|\Phi_{\rm cdw}\rangle$ is the single chain equivalent of $|\Psi_{\rm
cdw}\rangle$ given above, and its components have the same weighting factors
when the $J_{j}$ are disordered.

Figure 8:  The emergence of single particle “CDWs” in finite length spin
chains. (a) Local particle density $\langle\hat{n}_{j}\rangle$ is plotted for
each site of 40 and 41 site chains, scaled by length $N$. For the given
parameters, $\tilde{\Omega}=0.1$ and $|\zeta|=0.005$, and taking uniform
$J_{j}$, the odd-length chain has a single particle delocalized across all odd
sites, hence the average particle density is $\bar{n}=1/2N\approx 0.012$. The
even length chains do not have charge density waves, and have much smaller
average density $\bar{n}=N|\tilde{\Omega}|^{4}/8\approx 0.0005$. The
modulation of local density across the chain is due to the highly-correlated
two-particle state. (b) Particle density $\bar{n}$ is plotted vs. drive
strength $|\tilde{\Omega}|$ for two odd-length chains ($N=11$ and $N=10001$)
in the odd-length CDW driving strength regime (cf. Eq. (38)) for two different
$|\zeta|^{2}=10^{-4}$ and $|\zeta|^{2}=10^{-2}$. In the limit $|\zeta|\to 0$
(dashed curves), the single particle CDW persists for any arbitrarily small
$|\tilde{\Omega}|>0$. In both plots, the particle density is
$\bar{n}=\frac{1}{N}\sum_{j}\langle\hat{\sigma}_{A,j}^{+}\hat{\sigma}_{A,j}^{-}\rangle$.

To see why these single particle states emerge, consider the steady state, Eq.
(29), in the limit $|\zeta|\to 0$:
$|\psi_{Q}\rangle=e^{\sqrt{N}\hat{Q}/\tilde{\Omega}^{2}}|\psi_{\infty}\rangle$.
In this limit, all holes in the state are created by some power of $\hat{Q}$
acting on $|\psi_{\infty}\rangle$, with the component
$\hat{Q}^{m}|\psi_{\infty}\rangle$ having $2m$ holes. For an $N$-length chain,
$\hat{Q}$ can thus act up to $m=\lfloor N/2\rfloor$ times to produce new
states, after which $\hat{Q}^{m>\lfloor N/2\rfloor}|\psi_{\infty}\rangle=0$.
For an even chain, all $2m=N$ particles can be removed, taking the filled
state to vacuum:

$\displaystyle\hat{Q}^{\lfloor
N/2\rfloor}|\psi_{\infty}\rangle\propto|0\rangle,~{}~{}N~{}{\rm even}.$ (36)

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 $|\bullet\rangle$ 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:

$\displaystyle\hat{Q}^{\lfloor N/2\rfloor}|\psi_{\infty}\rangle$
$\displaystyle\propto|\bullet\circ\circ\circ\cdots\rangle$ (37)
$\displaystyle\,-|\circ\circ\bullet\circ\cdots\rangle+\cdots,~{}~{}N~{}{\rm
odd}.$

Thus for $|\zeta|\ll 1$, there is a regime of sufficiently weak
$|\tilde{\Omega}|$ 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
$\bar{n}=N|\tilde{\Omega}|^{4}/8$ of the even chain is found using Eq. (143)
in App. (J).

The analysis follows analogously in the limit $|\zeta|\to\infty$, where now
$|\psi_{Q}\rangle=e^{\sqrt{N}\hat{Q}/\tilde{\Omega}^{2}}\hat{\tau}_{1}|\psi_{\infty}\rangle$
(as $\hat{\tau}_{j}$ commutes with $\hat{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
$|\psi_{Q}\rangle$ to the lowest few orders in $|\tilde{\Omega}|$. Leaving the
details to App. J, we can show that two distinct scales emerge: a
$N$-dependent upper limit on drive strength, $|\tilde{\Omega}|^{2}\ll 1/N$,
and a $|\zeta|$\- and $N$-dependent lower limit that differs for even and odd
length chains:

$\displaystyle\frac{|\zeta|^{2}}{N}\ll|\tilde{\Omega}|^{2}\ll\frac{1}{N},$
$\displaystyle\quad N~{}{\rm odd},$ (38)
$\displaystyle\frac{1}{N|\zeta|^{2}}\ll|\tilde{\Omega}|^{2}\ll\frac{1}{N},$
$\displaystyle\quad N~{}{\rm even}.$ (39)

When $|\tilde{\Omega}|^{2}>1/N$, many particles can populate the chain,
destroying the CDW ordering, and when $|\tilde{\Omega}|^{2}<|\zeta|^{\pm 2}/N$
(for $N$ odd/even), the average particle density vanishes as
$\sim|\tilde{\Omega}|^{2}$. Here, Eq. (34) is no longer an equal mixture of
$|0\rangle$ and $|\Phi_{\rm cdw}\rangle$, but increasingly weighted toward
vacuum with $|\tilde{\Omega}|\to 0$. The emergence of CDWs, and the dependence
with both $N$ and $|\zeta|$, 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 $\Omega$, 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
$\Omega\rightarrow\infty$ 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 $|\tilde{\Omega}|\gtrsim 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
$\Omega\ll\gamma$ (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.

###  V.1 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 $|\Omega/\Gamma|\to\infty$, the steady state of
the $N=1$ double chain is not unique [7, 9, 8], but can be any state of the
form $\hat{\rho}_{1}=\nu|S\rangle\langle
S|+\frac{1}{4}(1-\nu)\hat{\mathds{I}}$ for any $-1/3\leq\nu\leq 1$, see App. M
for details. One can readily show that the maximally mixed state
$\hat{\mathbb{I}}/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.

###  V.2 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.

Figure 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 $\gamma\tau_{\rm rel}$ (red) and the
infidelity of the steady state to the maximally entangled state
$|\psi_{\infty}\rangle$ (cf. Eq. (12)),
$1-\langle\psi_{\infty}|\hat{\rho}|\psi_{\infty}\rangle$ (black), are shown as
functions of the hopping rate $\bar{J}/\gamma$ for $N=3$, $\Omega/\gamma=10$,
and $\Delta=0$. For the qubit-qutrit scheme, the relaxation time is optimized
over $\eta$. We also plot the relaxation time for the single chain system for
comparison. For a fixed state fidelity of 0.999 (achieved at
$\bar{J}=7\gamma$, dashed line), the relaxation times are $\gamma\tau_{\rm
2qb}=2.1\times 10^{4}$, $\gamma\tau_{\rm qutrit}=120$ and $\gamma\tau_{\rm
sc}=14$, for the qubit-qubit, qubit-qutrit, and single chain, respectively.

More concretely, starting from the double chain master equation (cf. Eq. (1))
for $N=1$ in the directional limit $\nu=1$, we promote qubit $B1$ to a qutrit
and modify its coupling to the waveguide via

$\displaystyle\hat{H}_{\mathrm{drive}}$
$\displaystyle=\frac{\Omega}{2}\left(\sigma_{A,1}^{x}+|0\rangle\langle
1|_{B,1}+|1\rangle\langle 0|_{B,1}\right)$ (40)
$\displaystyle+\frac{\Delta}{2}\left(\hat{\sigma}_{A,1}^{z}-\big{[}|1\rangle\langle
1|_{B,1}-|0\rangle\langle 0|_{B,1}\big{]}\right),$
$\displaystyle\hat{H}_{\mathrm{diss}}$
$\displaystyle=\frac{i\gamma}{2}\left(\sigma_{A,1}^{+}[|0\rangle\langle
1|_{B,1}+\eta|1\rangle\langle 2|_{B,1}]-{\rm h.c.}\right),$ (41)
$\displaystyle\hat{L}$
$\displaystyle=\sqrt{\gamma}\left(\sigma_{A,1}^{-}+|0\rangle\langle
1|_{B,1}+\eta|1\rangle\langle 2|_{B,1}\right),$ (42)

for which the master equation now reads
$\partial_{t}\hat{\rho}=-i[\hat{H}_{\mathrm{drive}}+\hat{H}_{\mathrm{diss}},\hat{\rho}]+\mathcal{D}[\hat{L}]\rho$.
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 $\eta$). The result is a dissipative
interaction that allows the state $|11\rangle$ to pass a single photon through
the waveguide at a rate $\eta\gamma$ and become $|02\rangle$, which can in
turn decay into the state $|01\rangle$. The effective interaction (no jump
Hamiltonian) of such a process is $-i\eta\gamma\sigma_{A}^{-}(|2\rangle\langle
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
$\hat{H}_{\mathrm{XX}}$. We stress that all remaining qubits are just qubits:
it is only $B1$ where we need to make use of the higher $|2\rangle$ level. A
direct calculation shows that the steady state found in Eq. (25) is still a
zero-energy eigenstate of the new $\hat{H}_{\rm
drive}+\hat{H}_{\mathrm{diss}}$, as well as the new jump operator $\hat{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
$\tau_{\rm 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 $\eta$. These results are
as shown in Fig. 9(b). Here, we fix all other parameters except the uniform
hopping rate $\bar{J}$, and show how both the fidelity with the ideal
dimerized entangled state and $\tau_{\mathrm{rel}}$ vary with $\bar{J}$. For
the all-qubit chain, small $\bar{J}$ makes $|\tilde{\Omega}|$ large, thus the
fidelity to the maximally entangled state $|\psi_{\infty}\rangle$ (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 $\eta$ as a function of
$\bar{J}/\gamma$, 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”
definite-parity 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 $1$D
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 [29, 16]. 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.

_Acknowledgements._ We thank P. Rabl and J. Agustí for useful discussions. We
also thank the anonymous referee who pointed out the connection of our hole
pairing construction to quantum scars, and the existence of strong symmetries
in a modified model with Creutz ladder style couplings. This work was
supported by the National Science Foundation QLCI HQAN (NSF Award No.
2016136). AL, AP, MY, YXW and AC acknowledge support from the Army Research
Office under Grant No. W911NF-23-1-0077, and from the Simons Foundation
through a Simons Investigator Award (Grant No. 669487). This work was
completed in part with resources provided by the University of Chicago’s
Research Computing Center.

_Note added:_ While completing our manuscript, we became aware of a related
but independent work on autonomously stabilizing many-qubit entanglement;
unlike our study, the setup in this work used squeezed light and explicitly
directional qubit-qubit couplings [51].

## References

  * Poyatos  _et al._ [1996] J. F. Poyatos, J. I. Cirac, and P. Zoller, Quantum Reservoir Engineering with Laser Cooled Trapped Ions, Physical Review Letters 77, 4728 (1996).
  * Plenio and Huelga [2002] M. B. Plenio and S. F. Huelga, Entangled Light from White Noise, Physical Review Letters 88, 197901 (2002).
  * Kraus and Cirac [2004] B. Kraus and J. I. Cirac, Discrete Entanglement Distribution with Squeezed Light, Physical Review Letters 92, 013602 (2004).
  * Stannigel  _et al._ [2012] K. Stannigel, P. Rabl, and P. Zoller, Driven-dissipative preparation of entangled states in cascaded quantum-optical networks, New Journal of Physics 14, 063014 (2012).
  * Schirmer and Wang [2010] S. G. Schirmer and X. Wang, Stabilizing open quantum systems by Markovian reservoir engineering, Physical Review A 81, 062306 (2010).
  * Motzoi  _et al._ [2016] F. Motzoi, E. Halperin, X. Wang, K. B. Whaley, and S. Schirmer, Backaction-driven, robust, steady-state long-distance qubit entanglement over lossy channels, Physical Review A 94, 032313 (2016).
  * Govia  _et al._ [2022] L. C. G. Govia, A. Lingenfelter, and A. A. Clerk, Stabilizing two-qubit entanglement by mimicking a squeezed environment, Physical Review Research 4, 023010 (2022).
  * Brown  _et al._ [2022] T. Brown, E. Doucet, D. Ristè, G. Ribeill, K. Cicak, J. Aumentado, R. Simmonds, L. Govia, A. Kamal, and L. Ranzani, Trade off-free entanglement stabilization in a superconducting qutrit-qubit system, Nature Communications 13, 3994 (2022).
  * Doucet  _et al._ [2020] E. Doucet, F. Reiter, L. Ranzani, and A. Kamal, High fidelity dissipation engineering using parametric interactions, Physical Review Research 2, 023370 (2020).
  * Vollbrecht  _et al._ [2011] KarlGerdH. Vollbrecht, C. A. Muschik, and J. I. Cirac, Entanglement Distillation by Dissipation and Continuous Quantum Repeaters, Physical Review Letters 107, 120502 (2011).
  * Ma  _et al._ [2019] S.-l. Ma, X.-k. Li, X.-y. Liu, J.-k. Xie, and F.-l. Li, Stabilizing Bell states of two separated superconducting qubits via quantum reservoir engineering, Physical Review A 99, 042336 (2019).
  * Ma  _et al._ [2021] S.-l. Ma, J. Zhang, X.-k. Li, Y.-l. Ren, J.-k. Xie, M.-t. Cao, and F.-l. Li, Coupling-modulation–mediated generation of stable entanglement of superconducting qubits via dissipation, Europhysics Letters 135, 63001 (2021).
  * Agustí  _et al._ [2022] J. Agustí, Y. Minoguchi, J. M. Fink, and P. Rabl, Long-distance distribution of qubit-qubit entanglement using Gaussian-correlated photonic beams, Physical Review A 105, 062454 (2022).
  * Horn  _et al._ [2018] K. P. Horn, F. Reiter, Y. Lin, D. Leibfried, and C. P. Koch, Quantum optimal control of the dissipative production of a maximally entangled state, New Journal of Physics 20, 123010 (2018).
  * Diehl  _et al._ [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, Quantum states and phases in driven open quantum systems with cold atoms, Nature Physics 4, 878 (2008).
  * Pocklington  _et al._ [2022] A. Pocklington, Y.-X. Wang, Y. Yanay, and A. A. Clerk, Stabilizing volume-law entangled states of fermions and qubits using local dissipation, Physical Review B 105, L140301 (2022).
  * Zippilli  _et al._ [2015] S. Zippilli, J. Li, and D. Vitali, Steady-state nested entanglement structures in harmonic chains with single-site squeezing manipulation, Physical Review A 92, 032319 (2015).
  * Angeletti  _et al._ [2023] J. Angeletti, S. Zippilli, and D. Vitali, Dissipative stabilization of entangled qubit pairs in quantum arrays with a single localized dissipative channel, Quantum Science and Technology 8, 035020 (2023).
  * Zippilli  _et al._ [2013] S. Zippilli, M. Paternostro, G. Adesso, and F. Illuminati, Entanglement Replication in Driven Dissipative Many-Body systems, Physical Review Letters 110, 040503 (2013).
  * Yanay and Clerk [2018] Y. Yanay and A. A. Clerk, Reservoir engineering of bosonic lattices using chiral symmetry and localized dissipation, Physical Review A 98, 043615 (2018).
  * Zippilli and Vitali [2021] S. Zippilli and D. Vitali, Dissipative Engineering of Gaussian Entangled States in Harmonic Lattices with a Single-Site Squeezed Reservoir, Physical Review Letters 126, 020402 (2021).
  * Ma  _et al._ [2017] S. Ma, M. J. Woolley, I. R. Petersen, and N. Yamamoto, Pure Gaussian states from quantum harmonic oscillator chains with a single local dissipative process, Journal of Physics A: Mathematical and Theoretical 50, 135301 (2017).
  * Žnidarič [2010] M. Žnidarič, Exact solution for a diffusive nonequilibrium steady state of an open quantum chain, Journal of Statistical Mechanics: Theory and Experiment 2010, L05002 (2010).
  * Prosen [2008] T. Prosen, Third quantization: A general method to solve master equations for quadratic open Fermi systems, New Journal of Physics 10, 043026 (2008).
  * Sigurdsson  _et al._ [2017] H. Sigurdsson, A. J. Ramsay, H. Ohadi, Y. G. Rubo, T. C. H. Liew, J. J. Baumberg, and I. A. Shelykh, Driven-dissipative spin chain model based on exciton-polariton condensates, Physical Review B 96, 155403 (2017).
  * Popkov  _et al._ [2019] V. Yu. Popkov, D. Karevski, and G. M. Schütz, Exact Results for the Isotropic Spin-1/2 Heisenberg Chain With Dissipative Boundary Driving, Theoretical and Mathematical Physics 198, 296 (2019).
  * Yamanaka and Sasamoto [2023] K. Yamanaka and T. Sasamoto, Exact solution for the Lindbladian dynamics for the open XX spin chain with boundary dissipation, SciPost Physics 14, 112 (2023).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. Prosen, and P. Prelovšek, Many-body localization in the Heisenberg XXZ magnet in a random field, Physical Review B 77, 064426 (2008).
  * Žnidarič [2015] M. Žnidarič, Relaxation times of dissipative many-body quantum systems, Physical Review E 92, 042143 (2015).
  * Medvedyeva and Kehrein [2014] M. V. Medvedyeva and S. Kehrein, Power-law approach to steady state in open lattices of noninteracting electrons, Physical Review B 90, 205410 (2014).
  * Prosen and Pižorn [2008] T. Prosen and I. Pižorn, Quantum Phase Transition in a Far-from-Equilibrium Steady State of an XY Spin Chain, Physical Review Letters 101, 105701 (2008).
  * Prosen and Žnidarič [2009] T. Prosen and M. Žnidarič, Matrix product simulations of non-equilibrium steady states of quantum spin chains, Journal of Statistical Mechanics: Theory and Experiment 2009, P02035 (2009).
  * Prosen [2011] T. Prosen, Open XXZ Spin Chain: Nonequilibrium Steady State and a Strict Bound on Ballistic Transport, Physical Review Letters 106, 217206 (2011).
  * Parmee and Cooper [2020] C. D. Parmee and N. R. Cooper, Steady states of a driven dissipative dipolar XXZ chain, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 135302 (2020).
  * Vitagliano  _et al._ [2010] G. Vitagliano, A. Riera, and J. I. Latorre, Volume-law scaling for the entanglement entropy in spin-1/2 chains, New Journal of Physics 12, 113049 (2010).
  * Roberts and Clerk [2020] D. Roberts and A. A. Clerk, Driven-Dissipative Quantum Kerr Resonators: New Exact Solutions, Photon Blockade and Quantum Bistability, Physical Review X 10, 021022 (2020).
  * Roberts  _et al._ [2021] D. Roberts, A. Lingenfelter, and A. A. Clerk, Hidden Time-Reversal Symmetry, Quantum Detailed Balance and Exact Solutions of Driven-Dissipative Quantum Systems, PRX Quantum 2, 020336 (2021).
  * Mamaev  _et al._ [2019] M. Mamaev, I. Kimchi, M. A. Perlin, R. M. Nandkishore, and A. M. Rey, Quantum Entropic Self-Localization with Ultracold Fermions, Physical Review Letters 123, 130402 (2019).
  * Moudgalya  _et al._ [2020] S. Moudgalya, N. Regnault, and B. A. Bernevig, Eta-pairing in Hubbard models: From spectrum generating algebras to quantum many-body scars, Physical Review B 102, 085140 (2020).
  * Žnidarič [2013] M. Žnidarič, Coexistence of Diffusive and Ballistic Transport in a Simple Spin Ladder, Physical Review Letters 110, 070602 (2013).
  * Metelmann and Clerk [2015] A. Metelmann and A. A. Clerk, Nonreciprocal Photon Transmission and Amplification via Reservoir Engineering, Physical Review X 5, 021025 (2015).
  * Colpa [1979] J. H. P. Colpa, Diagonalisation of the quadratic fermion Hamiltonian with a linear part, Journal of Physics A: Mathematical and General 12, 469 (1979).
  * Roberts and Clerk [2023a] D. Roberts and A. A. Clerk, Competition between Two-Photon Driving, Dissipation, and Interactions in Bosonic Lattice Models: An Exact Solution, Physical Review Letters 130, 063601 (2023a).
  * Roberts and Clerk [2023b] D. Roberts and A. A. Clerk, Exact Solution of the Infinite-Range Dissipative Transverse-Field Ising Model, Physical Review Letters 131, 190403 (2023b).
  * Note [1] The argument here does not prove uniqueness of the replicated steady state. In many cases, including that studied here, the lack of additional symmetries precludes additional steady states.
  * Note [2] For convenience, we normalize the operator by $1/\bar{J}\sqrt{N}$ so it produces an approximately normalized state when acting on $|\psi_{\infty}\rangle$ and other nearly-filled states (${||}\hat{Q}|\psi_{\infty}\rangle{||}^{2}\approx{||}\hat{Q}^{2}|\psi_{\infty}\rangle{||}^{2}\approx 1$).
  * Yang [1989] C. N. Yang, Eta pairing and off-diagonal long-range order in a Hubbard model, Physical Review Letters 63, 2144 (1989).
  * Zhang [1990] S. Zhang, Pseudospin symmetry and new collective modes of the Hubbard model, Physical Review Letters 65, 120 (1990).
  * Irfan  _et al._ [2024] A. Irfan, M. Yao, A. Lingenfelter, X. Cao, A. A. Clerk, and W. Pfaff, Loss resilience of driven-dissipative remote entanglement in chiral waveguide quantum electrodynamics (2024), arxiv:2404.00142 [quant-ph] .
  * Alaeian  _et al._ [2019] H. Alaeian, C. W. S. Chang, M. V. Moghaddam, C. M. Wilson, E. Solano, and E. Rico, Creating lattice gauge potentials in circuit QED: The bosonic Creutz ladder, Physical Review A 99, 053834 (2019).
  * Agustí  _et al._ [2023] J. Agustí, X. H. H. Zhang, Y. Minoguchi, and P. Rabl, Autonomous Distribution of Programmable Multiqubit Entanglement in a Dual-Rail Quantum Network, Physical Review Letters 131, 250801 (2023).
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Über das Paulische Äquivalenzverbot, Zeitschrift für Physik 47, 631 (1928).
  * Buča and Prosen [2012] B. Buča and T. Prosen, A note on symmetry reductions of the Lindblad equation: Transport in constrained open spin chains, New Journal of Physics 14, 073007 (2012).
  * Kraus  _et al._ [2008] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Preparation of entangled states by quantum Markov processes, Physical Review A 78, 042307 (2008).
  * Moudgalya  _et al._ [2022] S. Moudgalya, B. A. Bernevig, and N. Regnault, Quantum many-body scars and Hilbert space fragmentation: A review of exact results, Reports on Progress in Physics 85, 086501 (2022).
  * Weinberg and Bukov [2017] P. Weinberg and M. Bukov, QuSpin: A Python package for dynamics and exact diagonalisation of quantum many body systems part I: Spin chains, SciPost Physics 2, 003 (2017).
  * Karamlou  _et al._ [2022] A. H. Karamlou, J. Braumüller, Y. Yanay, A. Di Paolo, P. M. Harrington, B. Kannan, D. Kim, M. Kjaergaard, A. Melville, S. Muschinske, B. M. Niedzielski, A. Vepsäläinen, R. Winik, J. L. Yoder, M. Schwartz, C. Tahan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Quantum transport and localization in 1d and 2d tight-binding lattices, npj Quantum Information 8, 1 (2022).
  * Fang  _et al._ [2017] K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Generalized non-reciprocity in an optomechanical circuit via synthetic magnetism and reservoir engineering, Nature Physics 13, 465 (2017).
  * Kannan  _et al._ [2023] B. Kannan, A. Almanakly, Y. Sung, A. Di Paolo, D. A. Rower, J. Braumüller, A. Melville, B. M. Niedzielski, A. Karamlou, K. Serniak, A. Vepsäläinen, M. E. Schwartz, J. L. Yoder, R. Winik, J. I.-J. Wang, T. P. Orlando, S. Gustavsson, J. A. Grover, and W. D. Oliver, On-demand directional microwave photon emission using waveguide quantum electrodynamics, Nature Physics 19, 394 (2023).
  * Vool  _et al._ [2016] U. Vool, S. Shankar, S. O. Mundhada, N. Ofek, A. Narla, K. Sliwa, E. Zalys-Geller, Y. Liu, L. Frunzio, R. J. Schoelkopf, S. M. Girvin, and M. H. Devoret, Continuous Quantum Nondemolition Measurement of the Transverse Component of a Qubit, Physical Review Letters 117, 133601 (2016).
  * Mirhosseini  _et al._ [2019] M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B. Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang, and O. Painter, Cavity quantum electrodynamics with atom-like mirrors, Nature 569, 692 (2019).
  * Somoroff  _et al._ [2023] A. Somoroff, Q. Ficheux, R. A. Mencia, H. Xiong, R. Kuzmin, and V. E. Manucharyan, Millisecond Coherence in a Superconducting Qubit, Physical Review Letters 130, 267001 (2023).

##  Appendix A Interacting fermion model of the dissipative spin chain

Given a boundary driven/dissipative spin chain, often the first step towards
obtaining a solution is performing a Jordan-Wigner transform into free
fermions [52]. Define the canonical (Dirac) fermions:

$\displaystyle\hat{c}_{j}$
$\displaystyle=\left(\prod_{i=1}^{j-1}\hat{\sigma}_{i}^{z}\right)\hat{\sigma}_{j}^{-}.$
(43)

Then we can rewrite the spin Hamiltonian [c.f. $\hat{H}_{A}$ in Eq. (10)] as

$\displaystyle\hat{H}_{\mathrm{fermi}}$
$\displaystyle=\frac{\Omega}{2}(\hat{c}_{1}+\hat{c}_{1}^{\dagger})+\Delta\hat{c}_{1}^{\dagger}\hat{c}_{1}-\frac{1}{2}\sum_{j=1}^{N-1}J_{j}\left(\hat{c}_{j}^{\dagger}\hat{c}_{j+1}+{\rm
h.c.}\right).$ (44)

The Hamiltonian is a sum of quadratic and linear fermion terms, and can be
exactly diagonalized (using the procedure outlined in [42]).

We are of course interested in the dissipative dynamics. The full master
equation in terms of fermions is then

$\displaystyle\dot{\hat{\rho}}$
$\displaystyle=-i[\hat{H}_{\mathrm{fermi}},\hat{\rho}]+\gamma\mathcal{D}[\hat{c}_{1}]\hat{\rho}.$
(45)

It also has both quadratic and linear fermion terms. While one might assume
that this master equation is exactly solvable, this is not the case. If the
fermionic Lindbladian contains both a Hamiltonian term linear in fermion
operators, along with linear dissipation, then this will generically
correspond to an interacting problem. The easiest way to see this is simply to
compute the equation of motion for the linear expectation value of a fermionic
operator under just the dissipative dynamics. Because the Hamiltonian has a
linear term, the even and odd moments are dynamically connected and so we must
consider these linear expectation values.

$\displaystyle\partial_{t}\langle\hat{c}_{j}\rangle$
$\displaystyle=\gamma\left\langle\mathcal{D}[\hat{c}_{1}]^{\dagger}\hat{c}_{j}\right\rangle=\frac{\gamma}{2}\left\langle\hat{c}_{1}^{\dagger}[\hat{c}_{j},\hat{c}_{1}]+[\hat{c}_{1}^{\dagger},\hat{c}_{j}]\hat{c}_{1}\right\rangle,$
(46) $\displaystyle=-\frac{\gamma}{2}\left\langle
4\hat{c}_{1}^{\dagger}\hat{c}_{1}\hat{c}_{j}+\delta_{1,j}\hat{c}_{j}\right\rangle.$
(47)

We see that these first moments are coupled to third moments. In a similar
manner, one can observe that all odd moments of degree $n$ couple to odd
moments of degree $n+2$, and so the equations of motions do not close on
themselves, signifying that this is not a simple free fermion model.

Figure 10:  Lindblad spectrum of a single chain with 2 qubits. We plot the
eigenvalues $\lambda$ of the single chain master equation Eq. (10) for
$N=2,\Delta=0$ and $J=\gamma$. Left: In the absence of any drive $\Omega=0$,
the spectrum has the normal mode form expected for a quadratic fermionic
Lindblad master equation. Right: For a non-zero drive $\Omega=0.6\gamma$, 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.

More formally, one could try and use the standard diagonalization technique
for a linear fermionic Hamiltonian [42] and introduce a fictitious fermion
$\hat{\eta}$ to homogenize the Hamiltonian and make everything quadratic. This
is equivalent to rewriting the Hamiltonian Eq. (44) as

$\displaystyle\hat{H}_{\mathrm{fermi,\eta}}$
$\displaystyle=\frac{\Omega}{2}(\hat{c}_{1}+\hat{c}_{1}^{\dagger})(\hat{\eta}-\hat{\eta}^{\dagger})+\Delta\hat{c}^{\dagger}\hat{c}$
$\displaystyle\ \ \ \
-\frac{1}{2}\sum_{j=1}^{N-1}J_{j}\left(\hat{c}_{j}^{\dagger}\hat{c}_{j+1}+{\rm
h.c.}\right).$ (48)

Now, the Majorana $(\hat{\eta}+\hat{\eta}^{\dagger})$ is conserved by the
Hamiltonian, and so if this were a closed system we would be done as the two
Hamiltonians $\hat{H}_{\mathrm{fermi,\eta}}$ and $\hat{H}_{\mathrm{fermi}}$
would be isospectral (with $\hat{H}_{\mathrm{fermi,\eta}}$ doubly degenerate).
However, this is an open system, and we need to consider the dissipation. In
this case, the Majorana $(\hat{\eta}+\hat{\eta}^{\dagger})$ only constitutes a
weak symmetry [53] of the Lindbladian as it anticommutes with the jump term:

$\displaystyle\\{\hat{\eta}+\hat{\eta}^{\dagger},\hat{c}_{1}\\}$
$\displaystyle=0.$ (49)

This is problematic, as without further modification, the dissipation will
cause unphysical jumps between the two conserved sectors of
$\hat{H}_{\mathrm{fermi},\eta}$.

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 $(\hat{\eta}+\hat{\eta}^{\dagger})$ a strong symmetry
[53]. This is achieved by the master equation

$\displaystyle\dot{\hat{\rho}}$
$\displaystyle=-i[\hat{H}_{\mathrm{fermi,\eta}},\hat{\rho}]+\gamma\mathcal{D}[(-1)^{\hat{\eta}^{\dagger}\hat{\eta}}\hat{c}_{1}]\hat{\rho}.$
(50)

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 $\eta$, 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 $\hat{c}_{1}$,
as the action of the superoperator $\mathcal{D}[\hat{c}_{1}]$ cannot be
written solely in terms of the quadratic operator
$\hat{c}^{\dagger}_{1}\hat{c}_{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 $\lambda$) for the $N=2$ version of the master equation Eq. (10)
For $\Omega=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

###  B.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

$\displaystyle\hat{H}$
$\displaystyle=\hat{\sigma}^{x}_{A,1}\hat{\sigma}^{x}_{A,2}+\hat{\sigma}^{y}_{A,1}\hat{\sigma}^{y}_{A,2}-\left(\hat{\sigma}^{x}_{B,1}\hat{\sigma}^{x}_{B,2}+\hat{\sigma}^{y}_{B,1}\hat{\sigma}^{y}_{B,2}\right)$
(51)

We can immediately observe that
$[\hat{H},\hat{S}^{2}_{A,B}]=[\hat{H},\hat{S}^{z}_{A,B}]=0$, and so we can
diagonalize it in terms of the singlet and triplet states
$|s_{A},m_{A}\rangle\otimes|s_{B},m_{B}\rangle$ where $s,m$ are the standard
quantum numbers of total spin and z-angular momentum. It is quite easy to
observe that

$\displaystyle(\hat{\sigma}^{x}_{1}\hat{\sigma}^{x}_{2}+\hat{\sigma}^{y}_{1}\hat{\sigma}^{y}_{2})|s,m\rangle$
$\displaystyle=\delta_{m,0}(-1)^{s+1}|s,m\rangle$ (52)
$\displaystyle\implies\hat{H}|s_{A},m_{A}\rangle\otimes|s_{B},m_{B}\rangle$
$\displaystyle=\left(\delta_{m_{A},0}(-1)^{s_{A}+1}-\delta_{m_{B},0}(-1)^{s_{B}+1}\right)|s_{A},m_{A}\rangle\otimes|s_{B},m_{B}\rangle$
(53)
$\displaystyle\implies\hat{H}|s_{A},m_{A}\rangle\otimes|s_{B},m_{B}\rangle$
$\displaystyle=0$ (54)

The only remaining piece of the puzzle is the observation that tensoring
together two copies of a fixed (even) parity state $|\psi\rangle$ is diagonal
in this spin basis.

$\displaystyle|\psi\rangle$
$\displaystyle=\sum_{i\in\\{0,1\\}}\psi_{i}|i\rangle_{A}\otimes|i\rangle_{B}$
(55) $\displaystyle\implies|\psi\rangle_{1}\otimes|\psi\rangle_{2}$
$\displaystyle=\sum_{i,j\in\\{0,1\\}}\psi_{i}\psi_{j}|ij\rangle_{A}\otimes|ij\rangle_{B}$
(56) $\displaystyle=\sum_{s,m}c_{s,m}|s,m\rangle\otimes|s,m\rangle,$ (57)

where in the last line we go from the computational basis in Eq. (56) to the
total spin basis in Eq. (57). Since the coefficient $\psi_{i}\psi_{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. $|00\rangle$ in the
computational basis is $|1,-1\rangle$ in the total spin basis already, and
vice versa for $|11\rangle\leftrightarrow|1,1\rangle$. Thus, it only requires
checking the cross terms. Denote
$|S\rangle=|0,0\rangle=(|01\rangle-|10\rangle)/\sqrt{2}$ the singlet state,
and $|T\rangle=|1,0\rangle=(|01\rangle+|10\rangle)/\sqrt{2}$ the triplet. Then

$\displaystyle|01\rangle\otimes|01\rangle+|10\rangle\otimes|10\rangle$
$\displaystyle=\frac{1}{2}(|S\rangle+|T\rangle)\otimes(|S\rangle+|T\rangle)+\frac{1}{2}(|S\rangle-|T\rangle)\otimes(|S\rangle-|T\rangle)$
$\displaystyle=|S\rangle\otimes|S\rangle+|T\rangle\otimes|T\rangle$ (58)

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
$|\overline{i}\rangle\equiv|1-i\rangle=\hat{\sigma}^{x}|i\rangle$)

$\displaystyle|\psi\rangle$
$\displaystyle=\sum_{i}\psi_{i}|i\rangle_{A}\otimes|\overline{i}\rangle_{B}$
(59) $\displaystyle\implies|\psi\rangle_{1}\otimes|\psi\rangle_{2}$
$\displaystyle=\sum_{i,j}\psi_{i}\psi_{j}|ij\rangle_{A}\otimes|\overline{i}\overline{j}\rangle_{B}$
(60)
$\displaystyle=\sum_{s,m}c_{s,m}|s,m\rangle\otimes|\overline{s,m}\rangle.$
(61)

However, bit flip on the entire $B$ chain leaves the Hamiltonian invariant:

$\displaystyle\hat{\sigma}^{x}_{B,1}\hat{\sigma}^{x}_{B,2}\hat{H}\hat{\sigma}^{x}_{B,1}\hat{\sigma}^{x}_{B,2}=\hat{H}$
(62)

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
$\hat{H}_{{\rm XX},i}$ as the XX Hamiltonian acting on the $i$ and $i+1$ pair
of spins, then the total Hamiltonian would be

$\displaystyle\hat{H}_{\mathrm{total}}|\Psi\rangle_{\mathrm{total}}=\left(\sum_{i=1}^{N-1}J_{i}\hat{H}_{{\rm
XX},i}\right)|\psi\rangle^{\otimes N}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N-1}|\psi\rangle_{1}\dots\left(J_{i}\hat{H}_{{\rm
XX},i}|\psi\rangle_{i}|\psi\rangle_{i+1}\right)\dots|\psi\rangle_{N}$ (63)
$\displaystyle=$ $\displaystyle 0,$ (64)

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.

###  B.2 Heisenberg Couplings

To observe that Heisenberg couplings can replicate any state, it is first
important to observe that, given a fixed (even) parity state
$|\psi\rangle=u|00\rangle+v|11\rangle$, then

$\displaystyle\left(\hat{\sigma}_{A,1}^{z}\hat{\sigma}_{A,2}^{z}-\hat{\sigma}_{B,1}^{z}\hat{\sigma}_{B,2}^{z}\right)|\psi\rangle\otimes|\psi\rangle=0$
(65)

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 $|\phi\rangle$, then by
the Schmidt decomposition, there exist local unitaries
$\hat{U}_{1},\hat{U}_{2}$ such that

$\displaystyle\left(\hat{U}_{1}\otimes\hat{U}_{2}\right)|\phi\rangle$
$\displaystyle=u|00\rangle+v|11\rangle\equiv|\psi\rangle$ (66)

for some $u,v$. Let’s define $\hat{H}_{H}$ to be the isotropic Heisenberg
Hamiltonian. Then we have that

$\displaystyle\hat{H}_{H}|\phi\rangle\otimes|\phi\rangle$
$\displaystyle=\hat{H}_{H}\left(U_{1}^{\dagger}\otimes
U_{1}^{\dagger}\right)_{A}\otimes\left(U_{2}^{\dagger}\otimes
U_{2}^{\dagger}\right)_{B}|\psi\rangle\otimes|\psi\rangle$ (67)

Since the isotropic Heisenberg Hamiltonian is invariant under uniform local
unitary rotations, it can be commuted through to annihilate the fixed parity
state:

$\displaystyle\hat{H}_{H}|\phi\rangle\otimes|\phi\rangle$ $\displaystyle=0$
(68)

as desired.

###  B.3 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
$\mathcal{L}_{0}$ acting on qubits at site $A_{0},B_{0}$ that stabilizes a
fixed parity state $|\psi\rangle$,
$\mathcal{L}_{0}(|\psi\rangle\langle\psi|)=0$. Next, we will extend this to
the next layer in the graph by defining the next set of qubits on sites
$A_{1},\dots,A_{n}$ and $B_{1},\dots,B_{n}$, so that the full system
Liouvillian is now

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{0}\otimes\mathds{1}^{\otimes n}-i[\hat{H}_{{\rm
XX},n},\cdot]$ (69) $\displaystyle\hat{H}_{{\rm XX},n}$
$\displaystyle=\sum_{s=A,B}\sum_{i=1}^{n}\mathrm{sgn}(s)J_{i}\left(\sigma_{s,0}^{x}\sigma_{s,i}^{x}+\sigma_{s,0}^{y}\sigma_{s,i}^{y}\right)$
(70)

where we define $\mathrm{sgn}(A)=1,\mathrm{sgn}(B)=-1$. Now at this point,
defining $|\Psi\rangle=|\psi\rangle^{\otimes n+1}$, it is simple to observe
that since the Hamiltonian is simply a sum of terms acting independently on
the different dimers, $\hat{H}_{{\rm XX},n}|\Psi\rangle=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-$\frac{1}{2}$$\times$spin-$\frac{1}{2}$ Hilbert space
which share $m=\pm 1$ states but have orthogonal $m=0$ states. The $m=-1$
state is the two-qubit vacuum $|00\rangle$, the $m=+1$ state is the double-
excited state $|11\rangle$, and the orthogonal pair of $m=0$ states are
$|S\rangle=(|01\rangle-|10\rangle)/\sqrt{2}$ and
$|T\rangle=(|01\rangle+|10\rangle)/\sqrt{2}$.

We define the lowering operators $\hat{\tau}_{j}$ and $\hat{\lambda}_{j}$ on
each site $j$ that each destroy one of $|S_{j}\rangle$ or $|T_{j}\rangle$. The
$\hat{\tau}_{j}$ destroy $|S_{j=1,3,5,\dots}\rangle$ and
$|T_{j=2,4,6,\dots}\rangle$ – the $|\bullet\rangle$ particles – and the
$\hat{\lambda}_{j}$ destroy the opposite $m=0$ states – the
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
particles. Explicitly, in terms of the qubit operators, the lowering operators
are:

$\displaystyle\hat{\tau}_{j}$
$\displaystyle\equiv\hat{\sigma}_{B,j}^{-}+(-1)^{j}\hat{\sigma}_{A,j}^{-}$
(71) $\displaystyle=\sqrt{2}\begin{cases}|00_{j}\rangle\langle
S_{j}|-|S_{j}\rangle\langle 11_{j}|&j~{}{\rm odd}\\\ |00_{j}\rangle\langle
T_{j}|-|T_{j}\rangle\langle 11_{j}|&j~{}{\rm even}\end{cases},$
$\displaystyle\hat{\lambda}_{j}$
$\displaystyle\equiv\hat{\sigma}_{B,j}^{-}-(-1)^{j}\hat{\sigma}_{A,j}^{-}$
(72) $\displaystyle=\sqrt{2}\begin{cases}|00_{j}\rangle\langle
S_{j}|+|S_{j}\rangle\langle 11_{j}|&j~{}{\rm even}\\\ |00_{j}\rangle\langle
T_{j}|+|T_{j}\rangle\langle 11_{j}|&j~{}{\rm odd}\end{cases}.$

Each ladder operator and its adjoint forms a spin-1 representation of SU(2),
with the same $\hat{S}^{z}_{j}$ completing the algebra for each. That is, the
commutation relations
$[\hat{\tau}_{j}^{\dagger},\hat{\tau}_{k}]=[\hat{\lambda}_{j}^{\dagger},\hat{\lambda}_{k}]=2\hat{S}^{z}_{j}\delta_{jk}$
are simultaneously satisfied for the operator

$\displaystyle\hat{S}^{z}_{j}=\frac{1}{2}\left(\hat{\sigma}_{A,j}^{z}+\hat{\sigma}_{B,j}^{z}\right)=|11_{j}\rangle\langle
11_{j}|-|00_{j}\rangle\langle 00_{j}|.$ (73)

The two flavors of ladder operators on each site $j$ do not commute with each
other:

$\displaystyle[\hat{\tau}^{\dagger}_{j},\hat{\lambda}_{j}]=[\hat{\lambda}^{\dagger}_{j},\hat{\tau}_{j}]$
$\displaystyle=\hat{\sigma}_{B,j}^{z}-\hat{\sigma}_{A,j}^{z}$ (74)
$\displaystyle=2\left(|S_{j}\rangle\langle T_{j}|+{\rm h.c.}\right).$

Notice that these commutators are Hermitian and act to change flavor on a
given dimer site. Finally products of same-site ladder operators are:

$\displaystyle\hat{\tau}_{j}^{2}$
$\displaystyle=-\hat{\lambda}_{j}^{2}=-2(-1)^{j}|00_{j}\rangle\langle
11_{j}|,$ (75) $\displaystyle\hat{\tau}_{j}\hat{\lambda}_{j}$
$\displaystyle=0,$ (76)

This completes the dimer algebra.

We use the following notation of states in the dimer representation:
$|\circ_{j}\rangle$ and $|(11)_{j}\rangle$ are the vacuum and double-excited
state, respectively. The single-excited states are
$|\bullet_{j}\rangle=\frac{1}{\sqrt{2}}\hat{\tau}^{\dagger}_{j}|\circ_{j}\rangle$
and
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{j}\rangle=\frac{1}{\sqrt{2}}\hat{\lambda}^{\dagger}_{j}|\circ_{j}\rangle$.
Finally, we define the steady-state subspace

$\displaystyle\mathcal{H}_{s}={\rm
span}\,\Big{\\{}\sum_{j}|\circ_{j}\rangle,|\bullet_{j}\rangle\Big{\\}}$ (77)

(so-called because $|\psi_{Q}\rangle\in\mathcal{H}_{s}$) which contains all
states with only $|{\bullet}_{j}\rangle$ and $|{\circ}_{j}\rangle$. The
projector into $\mathcal{H}_{s}$ is given by

$\displaystyle\hat{\mathcal{P}}_{s}=\bigotimes_{j}\left(|\circ_{j}\rangle\langle\circ_{j}|+|\bullet_{j}\rangle\langle\bullet_{j}|\right).$
(78)

Now we rewrite Eq. (1) in the dimer representation. The collective loss
operator (cf. Eq. (3)) is simply

$\displaystyle\hat{c}=\hat{\lambda}_{1},$ (79)

which makes immediately obvious that its dark subspace on site 1 is spanned by
$|\circ_{1}\rangle$ and $|\bullet_{1}\rangle$ since
$\hat{\lambda}_{1}|(11)_{1}\rangle=\sqrt{2}|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{1}\rangle$
and
$\hat{\lambda}_{1}|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{1}\rangle=\sqrt{2}|\circ_{1}\rangle$.
The Hamiltonian terms (cf. Eqs. (4),(5), and (6)) are

$\displaystyle\hat{H}_{\rm XX}$
$\displaystyle=\frac{1}{4}\sum_{j=1}^{N-1}J_{j}\left(\hat{\tau}^{\dagger}_{j}\hat{\lambda}_{j+1}+\hat{\lambda}^{\dagger}_{j}\hat{\tau}_{j+1}+{\rm
h.c.}\right),$ (80) $\displaystyle\hat{H}_{\rm drive}$
$\displaystyle=\frac{\Omega}{2}\left(\hat{\lambda}^{\dagger}_{1}+{\rm
h.c.}\right)-\frac{\Delta}{2}\big{[}\hat{\lambda}^{\dagger}_{1},\hat{\tau}_{1}\big{]},$
(81) $\displaystyle\hat{H}_{\rm diss}$
$\displaystyle=\frac{1}{4}i\nu\gamma\left(\hat{\lambda}^{\dagger}_{1}\hat{\tau}_{1}-{\rm
h.c.}\right).$ (82)

$\hat{H}_{\rm XX}$ becomes a flavor-changing exchange interaction and the Rabi
drives become a single drive acting on the
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
flavor. The drive detuning in $\hat{H}_{\rm drive}$ and the nonreciprocity-
induced exchange in $\hat{H}_{\rm diss}$ appear in different ways to achieve
the same effect: a flavor change that swaps $|\bullet_{1}\rangle$ and
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{1}\rangle$.
Note that the commutator in $\hat{H}_{\rm drive}$ is Hermitian, which can be
seen using Eq. (74).

##  Appendix D Hole-pairing and XX eigenstates in qubit chains and Fermi-
Hubbard chains

###  D.1 Hole-pairing operator

The hole-pairing operator $\hat{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 $\hat{H}_{\rm XX}$ (cf. Eq. (5)) within the
steady-state subspace $\mathcal{H}_{s}$ (cf. Eq. (77)) to produce new
eigenstates of $\hat{H}_{\rm XX}$ with the same energy and also in
$\mathcal{H}_{s}$. We prove this at the operator level by showing that when
$\hat{Q}$ acts on states within the steady-state subspace $\mathcal{H}_{s}$,
it commutes with both $\hat{H}_{\rm XX}$ and the projector into the subspace
$\hat{\mathcal{P}}_{s}$ (cf. Eq. (78)),

$\displaystyle[\hat{Q},\hat{H}_{\rm XX}]\hat{\mathcal{P}}_{s}$
$\displaystyle=0,$ (83)
$\displaystyle[\hat{Q},\hat{\mathcal{P}}_{s}]\hat{\mathcal{P}}_{s}$
$\displaystyle=0.$ (84)

Together, these imply that for a given $\hat{H}_{\rm XX}$ eigenstate
$|\phi\rangle\in\mathcal{H}_{s}$ with energy $E$, the action of $\hat{Q}$ on
$|\phi\rangle$ produces another eigenstate with energy $E$:

$\displaystyle(\hat{H}_{\rm XX}-E)(\hat{Q}|\phi\rangle)$
$\displaystyle=0,\quad\hat{Q}|\phi\rangle\in\mathcal{H}_{s}.$ (85)

In this way, we find that by repeated application of $\hat{Q}$ on an
$\hat{H}_{\rm XX}$ eigenstate, a tower of some finite number $n_{0}<\infty$ of
degenerate eigenstates is returned.

In general the commutators $[\hat{Q},\hat{H}_{\rm XX}]\neq 0$ and
$[\hat{Q},\hat{\mathcal{P}}_{s}]\neq 0$. A direct computation using Eqs. (78)
and (80) yields the commutators

$\displaystyle[\hat{Q},\hat{H}_{\rm XX}]$
$\displaystyle=\frac{1}{8\sqrt{N}}\sum_{j}(-1)^{j}\left[\frac{J_{j}^{2}}{\bar{J}}[\hat{\tau}_{j},\hat{\lambda}_{j}^{\dagger}]\left(\hat{\tau}_{j+1}^{2}-\hat{\tau}_{j-1}^{2}\right)-2\frac{J_{j}J_{j-1}}{\bar{J}}\hat{S}_{j}^{z}\left(\hat{\tau}_{j+1}\hat{\lambda}_{j-1}-\hat{\tau}_{j-1}\hat{\lambda}_{j+1}\right)\right],$
(86) $\displaystyle[\hat{Q},\hat{\mathcal{P}}_{s}]$
$\displaystyle=-\frac{1}{\sqrt{N}}\sum_{j}(J_{j}/\bar{J})(-1)^{j}\bigotimes_{k\neq
j,j+1}\left(|\circ_{k}\rangle\langle\circ_{k}|+|\bullet_{k}\rangle\langle\bullet_{k}|\right)\times$
(87)
$\displaystyle\qquad\qquad\qquad\Big{[}|{\bullet}_{j}\circ_{j+1}\rangle\langle(11)_{j}{\bullet}_{j+1}|+|{\circ}_{j}{\bullet}_{j+1}\rangle\langle{\bullet}_{j}(11)_{j+1}|+|{\bullet}_{j}{\bullet}_{j+1}\rangle\langle(11)_{j}(11)_{j+1}|\Big{]}.$

Noting that
$\hat{\tau}_{j}^{2}|\bullet_{j}\rangle=\hat{\lambda}_{j}|\bullet_{j}\rangle=0$
and
$\hat{\tau}_{j}^{2}|\circ_{j}\rangle=\hat{\lambda}_{j}|\circ_{j}\rangle=0$, we
immediately see that within the steady-state subspace, $[\hat{Q},\hat{H}_{\rm
XX}]\hat{\mathcal{P}}_{s}=0$. By inspection,
$[\hat{Q},\hat{\mathcal{P}}_{s}]\hat{\mathcal{P}}_{s}=0$.

###  D.2 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 tight-binding 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\to\infty$, 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, $\hat{H}=\hat{H}_{0}+\hat{H}_{\rm
SOC}$, where

$\displaystyle\hat{H}_{0}$
$\displaystyle=J\sum_{j,\sigma}\left(\hat{c}^{\dagger}_{j,\sigma}\hat{c}_{j+1,\sigma}+{\rm
h.c.}\right)+\frac{U}{2}\sum_{j}\hat{n}_{j,\uparrow}\hat{n}_{j,\downarrow},$
(88) $\displaystyle\hat{H}_{\rm SOC}$
$\displaystyle=\frac{\Omega}{2}\sum_{j}(-1)^{j}(\hat{c}_{j,\uparrow}^{\dagger}\hat{c}_{j,\downarrow}+{\rm
h.c.}).$ (89)

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, $\Omega$ is the laser Rabi drive strength, and
$\hat{n}_{j,\sigma}=\hat{c}_{j,\sigma}^{\dagger}\hat{c}_{j,\sigma}$. As in
Ref. [38], we define a new set of fermion operators

$\displaystyle\hat{a}_{j,\uparrow}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\hat{c}_{j,\uparrow}+(-1)^{j}\hat{c}_{j,\downarrow}\right),$
(90) $\displaystyle\hat{a}_{j,\downarrow}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\hat{c}_{j,\uparrow}-(-1)^{j}\hat{c}_{j,\downarrow}\right),$
(91)

which obey the canonical anticommutation relations
$\\{\hat{a}_{j,\sigma},\hat{a}_{k,\tau}\\}=\delta_{jk}\delta_{\sigma\tau}$. In
this basis, the Hamiltonian is

$\displaystyle\hat{H}$
$\displaystyle=J\sum_{j,\sigma}\left(\hat{a}_{j,\sigma}^{\dagger}\hat{a}_{j+1,\bar{\sigma}}+\mathrm{h.c.}\right)+\frac{U}{2}\sum_{j}\hat{n}_{j,\uparrow}\hat{n}_{j,\downarrow}$
(92)
$\displaystyle\qquad+\frac{\Omega}{2}\sum_{j}\Big{(}\hat{n}_{j,\uparrow}-\hat{n}_{j,\downarrow}\Big{)},$

where $\bar{\sigma}$ denotes the spin-flip of $\sigma$ and now
$\hat{n}_{j,\sigma}=\hat{a}_{j,\sigma}^{\dagger}\hat{a}_{j,\sigma}$. We thus
have the desired spin-flip hopping. Also note the energy splitting $\Omega$
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 $\Omega=U\to\infty$, here
we wish to consider a slightly different limit: $U\to\infty$ with
$\Omega,J<\infty$. This effects the on-site hard-core repulsion between
opposite spins, thus we have the effective model

$\displaystyle\hat{H}_{\infty}$
$\displaystyle=\frac{\Omega}{2}\sum_{j}\Big{(}\hat{n}_{j,\uparrow}-\hat{n}_{j,\downarrow}\Big{)}$
(93)
$\displaystyle+J\prod_{k}\Big{(}1-\hat{n}_{k,\uparrow}\hat{n}_{k,\downarrow}\Big{)}\sum_{j\sigma}\left(\hat{a}_{j,\sigma}^{\dagger}\hat{a}_{j+1,\overline{\sigma}}+\mathrm{h.c.}\right).$

With the hard-core repulsion, the $N$-particle ferromagnetic states
$|\Xi_{\uparrow}\rangle=|{\uparrow}{\uparrow}{\uparrow}\cdots\rangle$ and
$|\Xi_{\downarrow}\rangle=|{\downarrow}{\downarrow}{\downarrow}\cdots\rangle$
are eigenstates of $\hat{H}_{\infty}$, with energies $+N\frac{\Omega}{2}$ and
$-N\frac{\Omega}{2}$, respectively. We now introduce the fermionic hole-
pairing operators (one for each spin)

$\displaystyle\hat{Q}_{\sigma}=\frac{1}{\sqrt{N}}\sum_{j}(-1)^{j}\hat{a}_{j,\sigma}\hat{a}_{j+1,\sigma},$
(94)

that create pairs of adjacent holes with staggered phases, in perfect analogy
with Eq. (22). We also denote by $\hat{\mathcal{P}}_{\sigma}$ the projector
into each spin subspace, in analogy with Eq. (78).

Just as with the dimer spin chain case, we find that the hole-pairing
operators here commute with the tight-binding Hamiltonian when restricted to
the appropriate subspace, but the total energy changes due to the removal of
two particles:

$\displaystyle[\hat{Q}_{\uparrow},\hat{H}_{\infty}]\hat{\mathcal{P}}_{\uparrow}=-\Omega\hat{Q}_{\uparrow},\quad[\hat{Q}_{\downarrow},\hat{H}_{\infty}]\hat{\mathcal{P}}_{\downarrow}=+\Omega\hat{Q}_{\downarrow}.$
(95)

Thus, when acting on their respective $|\Xi_{\sigma}\rangle$, powers of
$\hat{Q}_{\sigma}$ produce towers of Hamiltonian eigenstates with energies
ranging from $\pm N\frac{\Omega}{2}$ to either $0$ or $\pm\frac{\Omega}{2}$
for either even-length or odd-length chains, respectively. Notice that due to
the SOC energy splitting, $\hat{Q}_{\sigma}$ raises or lowers the eigenstate
energy by a multiple of the SOC drive strength $\Omega/2$, furthering the
analogy with $\eta$-pairing [47, 48], for which $\eta$-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 $|\psi_{Q}\rangle$, Eq. (25), is indeed a
pure steady state of Eq. (1). For $|\psi_{Q}\rangle$ 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
$\hat{c}|\bullet_{1}\rangle=\hat{c}|\circ_{1}\rangle=0$ and $|\psi_{Q}\rangle$
has only those components on site 1. It remains to show that
$|\psi_{Q}\rangle$ is an eigenstate of $\hat{H}$, which we do by way of a
variational ansatz.

We consider a two-parameter variational ansatz
$|\psi^{\prime}[\alpha,\beta]\rangle$ and show that it is an exact eigenstate
for a set of uniquely determined variational parameters $\alpha$, $\beta$.
First note that for $N=1$, the steady state (cf. Eq. (7)) is

$\displaystyle|\psi_{1}\rangle=\left(1+\frac{\Gamma}{\Omega}\hat{\tau}_{1}\right)|\bullet\rangle.$
(96)

This is a zero-energy eigenstate of the boundary Hamiltonian

$\displaystyle\hat{H}_{1}\equiv\hat{H}_{\rm drive}+\hat{H}_{\rm diss}$ (97)

(cf. Eqs. (81) and (82)). Thus we construct a variational ansatz which
includes this wavefunction in the $N=1$ case, but we replace $\Gamma/\Omega$
with a variational parameter $\beta$ for generality. This correction to
$|\psi_{\infty}\rangle$ 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
$\hat{H}_{\rm XX}$. One might expect that a linear combination of zero energy
$\hat{H}_{\rm XX}$ eigenstates with all possible number of hole pairs is
needed. A “hole-pair condensate” i.e., the exponential of $\hat{Q}$ acting on
$|\psi_{\infty}\rangle$, is a simple way to achieve that, thus we make the
ansatz

$\displaystyle|\psi^{\prime}[\alpha,\beta]\rangle=(1+\beta\hat{\tau}_{1})e^{\alpha\hat{Q}}|\psi_{\infty}\rangle.$
(98)

Evaluating $\hat{H}|\psi^{\prime}[\alpha,\beta]\rangle$, there are three
nonzero terms

$\displaystyle\hat{H}|\psi^{\prime}[\alpha,\beta]\rangle=\left[\hat{H}_{1}+\beta\hat{H}_{1}\hat{\tau}_{1}+\beta[\hat{H}_{\rm
XX},\hat{\tau}_{1}]\right]e^{\alpha\hat{Q}}|\psi_{\infty}\rangle,$ (99)

as $\hat{H}_{\rm XX}e^{\alpha\hat{Q}}|\psi_{\infty}\rangle=0$. In what
follows, we use the fact that
$e^{\alpha\hat{Q}}|\psi_{\infty}\rangle\in\mathcal{H}_{s}$ (cf. Eq. (77))
contains no
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}\rangle$
or doublon states $|11\rangle$. Thus we may always write
$e^{\alpha\hat{Q}}|\psi_{\infty}\rangle=\hat{\mathcal{P}}_{s}e^{\alpha\hat{Q}}|\psi_{\infty}\rangle$
(cf. Eq. (78)). Evaluating the boundary Hamiltonian restricted to the the
steady state subspace yields

$\displaystyle\hat{H}_{\rm drive}\hat{\mathcal{P}}_{s}$
$\displaystyle=\frac{\Omega}{2}\hat{\lambda}_{1}^{\dagger}-\frac{\Delta}{2}\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1},$
(100) $\displaystyle\hat{H}_{\rm diss}\hat{\mathcal{P}}_{s}$
$\displaystyle=\frac{1}{4}i\nu\gamma\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1}.$
(101)

So we have
$\hat{H}_{1}\hat{\mathcal{P}}_{s}=\frac{1}{2}\Omega\hat{\lambda}_{1}^{\dagger}-\frac{1}{2}\Gamma\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1}$,
where $\Gamma$ is defined in Eq. (8). Likewise, the commutator is given by
$[\hat{H}_{\rm
XX},\hat{\tau}_{1}]\hat{\mathcal{P}}_{s}=\frac{1}{4}J_{1}\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1}\hat{\tau}_{2}$.
Thus we have

$\displaystyle\hat{H}$ $\displaystyle|\psi^{\prime}[\alpha,\beta]\rangle=$
(102)
$\displaystyle\left[\frac{1}{2}\Omega\hat{\lambda}_{1}^{\dagger}+\beta\frac{1}{4}J_{1}\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1}\hat{\tau}_{2}+\frac{1}{2}\big{(}\beta\Omega-\Gamma\big{)}\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1}\right]e^{\alpha\hat{Q}}|\psi_{\infty}\rangle.$

Commuting everything past the exponential $e^{\alpha\hat{Q}}$, one can readily
show that
$[\hat{\lambda}^{\dagger}_{1}\hat{\tau}_{1},\hat{Q}]\hat{\mathcal{P}}_{s}=[\hat{\lambda}^{\dagger}_{1}\hat{\tau}_{1}\hat{\tau}_{2},\hat{Q}]\hat{\mathcal{P}}_{s}=0$.
The only non-zero commutator is

$\displaystyle[\hat{\lambda}^{\dagger}_{1},e^{\alpha\hat{Q}}]\hat{\mathcal{P}}_{s}=-\alpha(J_{1}/2\bar{J}\sqrt{N})e^{\alpha\hat{Q}}\hat{\lambda}^{\dagger}_{1}\hat{\tau}_{1}\hat{\tau}_{2}.$
(103)

This commutator is evaluated by first noting that
$[\hat{\lambda}^{\dagger}_{1},\hat{Q}]\hat{\mathcal{P}}_{s}=-(J_{1}/2\bar{J}\sqrt{N})\hat{\lambda}^{\dagger}_{1}\hat{\tau}_{1}\hat{\tau}_{2}$
and
$[\hat{\lambda}^{\dagger}_{1}\hat{\tau}_{1}\hat{\tau}_{2},\hat{Q}]\hat{\mathcal{P}}_{s}=0$.
Thus $[[\hat{\lambda}^{\dagger}_{1},\hat{Q}],\hat{Q}]\hat{\mathcal{P}}_{s}=0$.
Then we can use the general results that for two operators $\hat{A}$,
$\hat{B}$ satisfying $[[\hat{A},\hat{B}],\hat{B}]=0$, the commutator
$[\hat{A},e^{\alpha\hat{B}}]=\alpha[\hat{A},\hat{B}]e^{\alpha\hat{B}}=\alpha
e^{\alpha\hat{B}}[\hat{A},\hat{B}]$. Using the hard-core constraint, we also
have $\hat{\lambda}^{\dagger}_{1}|\psi_{\infty}\rangle=0$, hence the action of
$\hat{H}$ on the variational ansatz reduces to the following two terms

$\displaystyle\hat{H}|\psi^{\prime}[\alpha,\beta]\rangle=$ (104)
$\displaystyle
e^{\alpha\hat{Q}}\left[\frac{J_{1}}{4}\left(\beta-\alpha\frac{\Omega}{\bar{J}\sqrt{N}}\right)\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1}\hat{\tau}_{2}+\frac{1}{2}\big{(}\beta\Omega-\Gamma\big{)}\hat{\lambda}_{1}^{\dagger}\hat{\tau}_{1}\right]|\psi_{\infty}\rangle.$

Since both terms involve the creation of a
$|{\mathrel{\raisebox{0.38748pt}{${\scriptscriptstyle\blacksquare}$}}}_{1}\rangle$,
the eigenstate must have eigenvalue 0, and each term must vanish separately as
they remove differing numbers of $|\bullet_{j}\rangle$. Immediately, we see
that letting $\beta=\Gamma/\Omega$ (as predicted by Eq. (96)) and
$\alpha=\sqrt{N}\Gamma\bar{J}/\Omega^{2}$ yields $|\psi_{Q}\rangle$ as
desired.

We conclude this appendix with an important comment on uniqueness. While we do
not have a proof that this is the _unique_ steady state of the master equation
for any $N>1$, we have reason to expect that it is unique so long as
$|\tilde{\Omega}|<\infty$. As noted above, the original $N=1$ problem does
have a unique steady state $|\psi_{1}\rangle$ (cf. Eq. (7)) for any finite
driving strength $|\Omega/\Gamma|<\infty$. Numerical exact diagonalization of
the Liouvillian for up to $N=4$ finds Eq. (25) to be the unique steady state.
Moreover, there is no obvious symmetry in the problem that would allow for a
steady state degeneracy (except in the limit $\Omega\to\infty$, see Sec. V),
thus we expect this steady state to be generically unique.

##  Appendix F Recursive steady state construction

###  F.1 Steady state recursion relation

The steady state $|\psi\rangle$ given by Eq. (25) can alternatively be
constructed recursively in the length of the chains $N$. While this
construction does not lead to further analytic insights into the structure and
properties of the steady state beyond what can be obtained from Eq. (25), it
does allow for the numerically amenable calculation of expectation values and
correlation functions. In what follows, it is convenient to define the
dimensionless parameters

$\displaystyle\zeta\equiv\sqrt{\frac{\Gamma}{\bar{J}}},\quad\eta_{j}\equiv\frac{J_{j}}{\bar{J}}$
(105)

The two $\zeta$ and $\tilde{\Omega}$ (cf. Eq. (26)) are sufficient to describe
overall properties of the steady state (i.e., $|\psi\rangle$ can be written
using only these two parameters when the $\\{J_{j}\\}$ are hidden away in
$\hat{Q}$), and the set of $\\{\eta_{j}\\}$ encode all intra-chain hopping
disorder.

The starting point for the recursive construction is the pair of steady states
$|\psi_{1}\rangle$ and $|\psi_{2}\rangle$ for the $N=1$ and $N=2$ chains,
respectively. The $N=1$ steady state $|\psi_{1}\rangle$ (cf. Eq. (7)) is given
in terms of $\tilde{\Omega}$ and $\zeta$ by

$\displaystyle|\psi_{1}\rangle$
$\displaystyle=\frac{1}{\mathcal{N}_{1}}\left[\sqrt{2}\frac{\zeta}{\tilde{\Omega}}|\circ_{1}\rangle+|\bullet_{1}\rangle\right],$
(106) $\displaystyle\mathcal{N}_{1}^{2}$
$\displaystyle=\frac{2|\zeta|^{2}}{|\tilde{\Omega}|^{2}}+1,$ (107)

where $\mathcal{N}_{1}$ is the normalization. For the later use in computing
correlation functions, it is crucial that the states be normalized. The $N=2$
steady state is

$\displaystyle|\psi_{2}\rangle$
$\displaystyle=\frac{1}{\mathcal{N}_{2}}\left[\frac{\eta_{1}}{\mathcal{N}_{1}\tilde{\Omega}^{2}}|{\circ}_{1}{\circ}_{2}\rangle-|\psi_{1}\rangle\otimes|{\bullet}_{2}\rangle\right],$
(108) $\displaystyle\mathcal{N}_{2}^{2}$
$\displaystyle=\frac{\eta_{1}^{2}}{\mathcal{N}_{1}^{2}|\tilde{\Omega}|^{4}}+1.$
(109)

By construction this is a dark state of the dissipation, and one can readily
verify that it is a zero-energy eigenstate of the Hamiltonian
$\hat{H}_{2}=\hat{H}_{1}+\hat{H}_{{\rm XX},2}$.

For all $n\geq 3$, the steady state is given recursively by

$\displaystyle|\psi_{n}\rangle$
$\displaystyle=\frac{1}{\mathcal{N}_{n}}\left[\frac{\eta_{n-1}}{\mathcal{N}_{n-1}\tilde{\Omega}^{2}}|\psi_{n-2}\rangle\otimes|{\circ}_{n-1}{\circ}_{n}\rangle\right.$
(110) $\displaystyle\qquad\qquad~{}~{}\left.+(-1)^{\lfloor
n/2\rfloor}|\psi_{n-1}\rangle\otimes|\bullet_{n}\rangle\right],$
$\displaystyle\mathcal{N}_{n}^{2}$
$\displaystyle=\frac{\eta_{n-1}^{2}}{\mathcal{N}_{n-1}^{2}|\tilde{\Omega}|^{4}}+1.$
(111)

One readily verifies that this is an eigenstate of the Hamiltonian
$\hat{H}_{n}=\hat{H}_{1}+\hat{H}_{{\rm XX},n}$ using the induction hypotheses
$\hat{H}_{n-1}|\psi_{n-1}\rangle=\hat{H}_{n-2}|\psi_{n-2}\rangle=0$. The
hypotheses hold for $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$, thus completing
the inductive proof that $\hat{H}_{n}|\psi_{n}\rangle=0$. Notice that the
pairing of holes in the steady state appears clearly in the recursion
relation. The strong $\tilde{\Omega}$ limit
$|\psi\rangle\to|\bullet\bullet\bullet\cdots\rangle$ is also evident by
neglecting terms at least $\mathcal{O}(1/\tilde{\Omega})$ in small
$1/\tilde{\Omega}\ll 1$.

###  F.2 Correlation functions

The recursive construction of $|\psi\rangle$ 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
$\hat{n}_{j}\equiv\frac{1}{2}\hat{\tau}_{j}^{\dagger}\hat{\tau}_{j}$. Since
the recursion relation for $|\psi_{n}\rangle$ is in the chain length $n$, it
is necessary to denote the chain length for which expectation values are
taken, e.g.

$\displaystyle\langle\hat{n}_{j}\rangle_{n}=\frac{1}{2}\langle\psi_{n}|\hat{\tau}_{j}^{\dagger}\hat{\tau}_{j}|\psi_{n}\rangle.$
(112)

Here the state is normalized by construction,
$\langle\psi_{n}|\psi_{n}\rangle=1$. Using the recursion relation Eq. (110),
we expand $\langle\hat{n}_{j}\rangle_{n}$ in terms of the expectation values
evaluated for $n-1$ and $n-2$ length chains:

$\displaystyle\langle\hat{n}_{j}\rangle_{n}=\frac{\eta_{n-1}^{2}}{\mathcal{N}_{n}^{2}\mathcal{N}_{n-1}^{2}|\tilde{\Omega}|^{4}}\langle\hat{n}_{j}\rangle_{n-2}+\frac{1}{\mathcal{N}_{n}^{2}}\langle\hat{n}_{j}\rangle_{n-1},$
(113)

where $\mathcal{N}_{k}$ are the normalization factors given by Eq. (111).

This expression tells us that the expectation value
$\langle\hat{n}_{j}\rangle_{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. (113). This
recursion terminates at the expectation value $\langle\hat{n}_{j}\rangle_{j}$,
i.e. for a length $n=j$ chain. The termination $\langle\hat{n}_{j}\rangle_{j}$
is readily evaluated using Eq. (110) directly and is

$\displaystyle\langle\hat{n}_{j}\rangle_{j}=\langle\psi_{j}|\hat{n}_{j}|\psi_{j}\rangle=\frac{1}{\mathcal{N}_{j}^{2}}.$
(114)

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 $\hat{Q}$ and $|\psi_{\infty}\rangle$ (cf.
Eq. (12)) constitute a restricted spectrum generating algebra (RSGA), and show
that these states have area-law entanglement entropy across a particular
bipartition, and thus are scar states of a simple ladder model.

###  G.1 RSGA in a nonintegrable model

In Ref. [39], a restricted spectrum generating algebra of order $M$ for a
system with Hamiltonian $\hat{H}$ is defined by an operator
$\hat{\eta}^{\dagger}$ and a state $|\psi_{0}\rangle$ satisfying: (i)
$\hat{H}|\psi_{0}\rangle=E_{0}|\psi_{0}\rangle$, (ii)
$\hat{H}_{1}|\psi\rangle=\mathcal{E}\hat{\eta}^{\dagger}|\psi_{0}\rangle$,
$\hat{H}_{n}|\psi_{0}\rangle=0~{}\forall~{}n,2\leq n\leq M$, and (iv)
$\hat{H}_{n}\neq 0,~{}n\leq M$, $\hat{H}_{n}=0,~{}n=M+1$. Here
$\hat{H}_{n}=[\hat{H}_{n-1},\hat{\eta}^{\dagger}]$ are successive commutators
with $\hat{\eta}^{\dagger}$, and $\hat{H}_{0}=\hat{H}$ is the system
Hamiltonian. Then the states $(\hat{\eta}^{\dagger})^{n}|\psi_{0}\rangle$ are
eigenstates of $\hat{H}$ with equally spaced energies $E_{0}+n\mathcal{E}$,
for all $n$ for which $(\hat{\eta}^{\dagger})^{n}|\psi_{0}\rangle\neq 0$.

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 $\hat{H}_{\delta}$):

$\displaystyle\hat{H}$ $\displaystyle=\hat{H}_{\rm XX}+\hat{H}_{\rm
rung}+\hat{H}_{\delta},$ (115) $\displaystyle\hat{H}_{\rm rung}$
$\displaystyle=\sum_{j=1}^{L}\left[g\left(\hat{\sigma}^{x}_{A,j}\hat{\sigma}^{x}_{B,j}+\hat{\sigma}^{y}_{A,j}\hat{\sigma}^{y}_{B,j}\right)+\mu\hat{\sigma}^{z}_{A,j}\hat{\sigma}^{z}_{B,j}\right],$
$\displaystyle\hat{H}_{\delta}$
$\displaystyle=\delta\sum_{i=1}^{L/2}\left(\hat{\sigma}_{A,2i}^{+}\hat{\sigma}_{B,2i}^{-}+{\rm
h.c.}\right).$ (116)

This model is nonintegrable (even for $\mu=0$), and it conserves total
$Z$-axis magnetization [40],

$\displaystyle\hat{M}=\sum_{j}\hat{\sigma}_{A,j}^{z}+\hat{\sigma}_{B,j}^{z}.$
(117)

The rung terms split the $|S\rangle$ and $|T\rangle$ states: $\hat{H}_{\rm
rung}|S\rangle=-2g|S\rangle$, $\hat{H}_{\rm rung}|T\rangle=+2g|T\rangle$, and
acts as a chemical potential by giving an energy cost to adding either holes
or particles: $\hat{H}_{\rm rung}|\circ\rangle=\mu|\circ\rangle$,
$\hat{H}_{\rm rung}|(11)\rangle=\mu|(11)\rangle$. The added $\hat{H}_{\delta}$
splits the two reference states $|\psi_{\infty}\rangle=|STST\cdots\rangle$ and
$|\tilde{\psi}_{\infty}\rangle=|TSTS\cdots\rangle$ via
$\hat{H}_{\delta}|S_{2i}\rangle=-2\delta|S_{2i}\rangle$ and
$\hat{H}_{\delta}|T_{2i}\rangle=2\delta|T_{2i}\rangle$, to lift some
degeneracies between the hole-paired states and other eigenstates of the
system.

It is straightforward to see that the $\hat{Q}$ operator constructed in App. D
and the reference state $|\psi_{\infty}\rangle$ (cf. Eq. (12)) form a RSGA of
order 2 under this Hamiltonian. In particular, we have

$\displaystyle\hat{H}|\psi_{\infty}\rangle$
$\displaystyle=E_{0}|\psi_{\infty}\rangle,~{}~{}E_{0}=\begin{cases}0&L~{}{\rm
even}\\\ -g&L~{}{\rm odd}\end{cases},$ (118)
$\displaystyle\hat{H}_{1}|\psi_{\infty}\rangle$
$\displaystyle=2\mu\hat{Q}|\psi_{\infty}\rangle,$
$\displaystyle\hat{H}_{2}|\psi_{\infty}\rangle$
$\displaystyle=0,~{}~{}\hat{H}_{3}=0,$

where $\hat{H}_{n}=[\hat{H}_{n-1},\hat{Q}]$ with $\hat{H}_{0}=\hat{H}$. Eq.
(118) implies that the set of hole-paired states

$\displaystyle|\psi_{n}\rangle$
$\displaystyle=\hat{Q}^{n}|\psi_{\infty}\rangle,$ (119)
$\displaystyle\hat{H}|\psi_{n}\rangle$
$\displaystyle=(E_{0}+2n\mu)|\psi_{n}\rangle,$ (120)

are a tower of exact eigenstates of the nonintegrable ladder, with equal
energy spacing $E_{n}-E_{n-1}=2n\mu$, for all $n<L/2$.

Similar calculations show that $\hat{Q}^{\dagger}$ acting on
$|\psi_{\infty}\rangle$ is also a RSGA of order 2, but this time creating
particle (i.e., $|(11)_{j}\rangle$) pairs. Similarly, if we define a
$\tilde{Q}$ which is identical to $\hat{Q}$ but with the replacement
$\hat{\tau}_{j}\hat{\tau}_{j+1}\mapsto\hat{\lambda}_{j}\hat{\lambda}_{j+1}$,
and a new reference state
$|\tilde{\psi}_{\infty}\rangle=|TSTSTS\cdots\rangle$, we also find that
$\tilde{Q}$ and $\tilde{Q}^{\dagger}$ with $|\tilde{\psi}_{\infty}\rangle$
form order 2 RSGAs, creating hole and particle pairs, respectively.

###  G.2 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\sim 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 $\alpha,~{}\beta$ be the
left, right halves of the chain, respectively). It is clear that the reference
states $|\psi_{\infty}\rangle,~{}|\tilde{\psi}_{\infty}\rangle$ 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\to\infty$: $S_{\alpha}\sim L^{0}$,
thus these eigenstates are quantum scars.

To demonstrate the anomolously low EE of the hole-paired states, we perform
numerical exact diagonalization of Eq. (115) 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 $\sim 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.

Figure 11:  Hole-paired states as scars in the double XX chain ladder model.
The entanglement entropy of the ladder model eigenstates (cf. Eq. (115)) of
the $\alpha$ 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 hole-paired states are denoted with stars; in each
magnetization sector there are two states corresponding to hole-pairs created
in each reference state $|\psi_{\infty}\rangle$ and
$|\tilde{\psi}_{\infty}\rangle$. For both plots the parameters are
$\mu/J=0.2$, $g/J=0.5$, and $\delta/g=0.3$.

###  G.3 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(\hat{\rho}_{\alpha,n})\sim L^{0}$, where
$\hat{\rho}_{\alpha,n}={\rm
tr}_{\beta}[\hat{Q}^{n}|\psi_{\infty}\rangle\langle\psi_{\infty}|\hat{Q}^{\dagger
n}]$. Here we denote the partitions by $\alpha,~{}\beta$; the $\alpha$
partition contains the first $L/2$ qubits of each chain $A$ and $B$. First, we
split the $\hat{Q}$ operator into three terms:

$\displaystyle\hat{Q}=\hat{Q}_{\alpha}+\hat{Q}_{\beta}+\hat{Q}_{L/2},$ (121)

where the first term acts only on the $\alpha$ partition (sites 1 through
$L/2$) and the second on the $\beta$ partition and the last term
$\hat{Q}_{L/2}\propto\hat{\tau}_{L/2}\hat{\tau}_{L/2+1}$ places a hole pair