Supplemental Material for “Hidden time-reversal in driven XXZ spin chains:
exact solutions and new dissipative phase transitions”

###### Contents

  1. I Coherent quantum absorber and exact solution
    1. I.1 CQA solution to the coherent driven-dissipative spin chain
    2. I.2 Proof of steady state uniqueness
    3. I.3 Observables
  2. II Phase transition
    1. II.1 Quantum to classical mapping
    2. II.2 Weak dissipation limit
    3. II.3 Free Energy and Critical Exponent
    4. II.4 Finite size scaling
  3. III Special points of anisotropy /J
    1. III.1 No phase transition at special points
    2. III.2 Resonances in magnetization at special points
    3. III.3 Survival of the phase transition in the thermodynamic limit
  4. IV Transport
  5. V Hidden Time Reversal Symmetry
  6. VI CQA solution for the incoherent model

##  I Coherent quantum absorber and exact solution

###  I.1 CQA solution to the coherent driven-dissipative spin chain

In this section, we derive the exact solution for the coherent boundary
driven-dissipative model. The bulk coherent dynamics is described by XXZ
Hamiltonian:

$\hat{H}_{\rm
XXZ}=\sum_{j=1}^{N-1}(J\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+J\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{y}+\Delta\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z}).$
(S1)

The system is driven by a coherent drive
$\frac{\Omega}{2}\hat{\sigma}_{N}^{x}$ at site $N$ and coupled to zero
temperature reservoir with jump operator $\hat{L}=\sqrt{\gamma}\sigma_{1}^{-}$
on site $1$. The dynamics is governed by Lindblad master equation

$\displaystyle\partial_{t}\hat{\rho}=\hat{\mathcal{L}}\hat{\rho}=-i[\hat{H}_{\rm
XXZ}+\frac{\Omega}{2}\hat{\sigma}_{N}^{x},\hat{\rho}]+\gamma\mathcal{D}[\hat{\sigma}_{1}^{-}]\hat{\rho}.$
(S2)

We seek the exact steady state solution $\mathcal{L}\hat{\rho}_{\rm ss}=0$ of
this model for all parameters.

We employ the “coherent quantum absorber” (CQA) method to find the exact
solution [1, 2]. The basic idea is to channel the dissipation of the system
(call it system “A”) via a chiral waveguide to an absorber system (“B”). For a
specific choice of absorber, the combined system + absorber relaxes to a pure
steady state $|\psi_{\rm CQA}\rangle$ [1]. Moreover, because the system is
upstream the absorber, one may trace out the absorber to obtain the A system
steady state: $\hat{\rho}_{\rm ss}={\rm Tr}_{B}(\outerproduct{\psi_{\rm
CQA}}{\psi_{\rm CQA}})$. A special class of dissipative systems possess a
hidden time reversal symmetry (hTRS) and as a consequence have a particularly
simple form of absorber: an mirrored copy of the system with the Hamiltonian
$\hat{H}^{(B)}=-\hat{H}^{(A)}$ [2]. As we show next, we can solve this model
using the CQA construction with $\hat{H}^{(B)}=-\hat{H}^{(A)}$, where
$\hat{H}^{(A)}=\hat{H}_{\rm XXZ}^{(A)}+(\Omega/2)\hat{\sigma}^{x}_{A,N}$ [See
Fig. 1 in the main text], hence the model possesses hTRS. The system A and
absorber B are coupled to a common chiral waveguide reservoir that mediates
collective dissipation
$\hat{L}^{(AB)}=\hat{\sigma}_{A,1}^{-}-\hat{\sigma}_{B,1}^{-}$. The chiral
waveguide induces an exchange Hamiltonian $\hat{H}_{\rm
c}=-\frac{i\gamma}{2}(\hat{\sigma}_{A,1}^{+}\hat{\sigma}_{B,1}^{-}-{\rm
h.c.})$ which is required to ensure nonreciprocal dynamics. The full system
and absorber master equation is

$\displaystyle\partial_{t}\hat{\rho}_{\rm AB}$ $\displaystyle=\mathcal{L}_{\rm
AB}\hat{\rho}_{\rm AB}=-i[\hat{H}^{(A)}_{\rm XXZ}-\hat{H}^{(B)}_{\rm
XXZ}+\frac{\Omega}{2}(\hat{\sigma}_{A,N}^{x}-\hat{\sigma}_{B,N}^{x})+\hat{H}_{\rm
c},\hat{\rho}_{\rm
AB}]+\gamma\mathcal{D}[\hat{\sigma}_{A,1}^{-}-\hat{\sigma}_{B,1}^{-}]\hat{\rho}_{\rm
AB}$ $\displaystyle=-i[\hat{H}^{(AB)},\hat{\rho}_{\rm
AB}]+\gamma\mathcal{D}[\hat{L}^{(AB)}]\hat{\rho}_{\rm AB}.$ (S3)

Note that one may trace out the absorber system from this master equation and
recover the system A master equation [cf. Eq. (S2)].

To find the pure steady state $\ket{\psi_{\rm CQA}}$ of Eq. (S3), it is
sufficient to construct a state that satisfies the following dark state
conditions [3],

$\hat{H}^{(AB)}\ket{\psi_{\rm CQA}}=0,\quad\hat{L}^{(AB)}\ket{\psi_{\rm
CQA}}=0.$ (S4)

(In general it is only required that $|\psi_{\rm CQA}\rangle$ is an
eigenstate, but for models with hTRS, the eigenvalue is zero [2].) To
construct the steady state satisfying these dark state conditions, we make a
matrix product state (MPS) ansatz for the solution following the previous work
[4]:

$|\psi_{\rm
CQA}\rangle=\sum_{s_{1},\cdots,s_{N}}\mathbf{v_{L}^{\dagger}}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}\mathbf{v_{R}}|s_{1}\cdots
s_{N}\rangle.$ (S5)

Here, each state $|s_{i}\rangle$ is a two-spin dimer state of the A spin and B
spin of site $j$ and is represented in the triplet-singlet basis
${|0_{i}\rangle}\equiv{|{\downarrow_{A,i}\downarrow_{B,i}}\rangle}$,
${|1_{i}\rangle}\equiv{|{\uparrow_{A,i}\uparrow_{B,i}}\rangle}$,
${|S_{i}\rangle}\equiv({|{\uparrow_{A,i}\downarrow_{B,i}}\rangle}-{|{\downarrow_{A,i}\uparrow_{B,i}}\rangle})/\sqrt{2}$,
and
${|T_{i}\rangle}\equiv({|{\uparrow_{A,i}\downarrow_{B,i}}\rangle}+{|{\downarrow_{A,i}\uparrow_{B,i}}\rangle})/\sqrt{2}$,
where ${|{\uparrow}\rangle}$, ${|{\downarrow}\rangle}$ are the eigenstates of
$\hat{\sigma}^{z}$ with eigenvalues $\pm 1$, respectively. The $\hat{A}_{s}$
matrices are taken to be translationally invariant (i.e., identical on every
site) but the full state $|\psi_{\rm CQA}\rangle$ is not translationally
invariant due to the boundary vectors $\mathbf{v_{L}}$ and $\mathbf{v_{R}}$.
Inspired by explicit solutions found for very small system sizes, we make the
ansatz that the only allowed spin states are the triplet subspace
$s_{i}\in\\{0_{i},1_{i},T_{i}\\}$. Then, we use the dark state conditions Eq.
(S4) to determine the MPS coefficients of $\hat{A}_{s}$, $\mathbf{v_{L}}$ and
$\mathbf{v_{R}}$.

The first dark state condition $\hat{L}_{\rm AB}|\psi_{\rm CQA}\rangle=0$
imposed by the jump term requires
$(\hat{\sigma}_{A,1}^{-}-\hat{\sigma}_{B,1}^{-})\ket{s_{1}}=0$. This restricts
the first site to $s_{1}\in\\{0_{1},T_{1}\\}$. This imposes the constraint

$\mathbf{v_{L}^{\dagger}}\hat{A}_{1}=0.$ (S6)

We are free to choose the left vector $\mathbf{v_{L}}$ and we take it to be
$\mathbf{v_{L}}=\ket{0}$.

The second dark state condition is $\hat{H}^{(AB)}|\psi_{\rm CQA}\rangle=0$.
This requires that the state be a zero-energy eigenstate of the Hamiltonian
$\hat{H}^{(AB)}$. Using our ansatz that each site $i$ of the double-system
chain contains only $s_{i}\in\\{0_{i},1_{i},T_{i}\\}$ (but no $T_{1}$ on site
1), we can see, mechanistically, how a dark state of $\hat{H}^{(AB)}$ can be
constructed: each action of $\hat{H}^{(AB)}$ on any state containing only the
allowed dimer states $s_{i}$ always produces a single “defect” dimer state
$|S_{i}\rangle$. Since this state is excluded from our ansatz, $|\psi_{\rm
CQA}\rangle$ must be a zero-energy eigenstate constructed so that all defects
produced by the action of $\hat{H}^{(AB)}$ destructively interfere. As an
illustrative example, consider the boundary coherent drive term $\hat{H}_{\rm
d}^{(AB)}=\frac{\Omega}{2}(\hat{\sigma}^{x}_{A,N}-\hat{\sigma}^{x}_{B,N})$. It
annihilates $|T_{N}\rangle$ and creates defects from $|0_{N}\rangle$ and
$|1_{N}\rangle$:

$\hat{H}_{\rm
d}^{(AB)}\ket{0_{N}}=\frac{\Omega}{\sqrt{2}}\ket{S_{N}},\quad\hat{H}_{\rm
d}^{(AB)}\ket{1_{N}}=-\frac{\Omega}{\sqrt{2}}\ket{S_{N}},\quad\hat{H}_{\rm
d}^{(AB)}\ket{T_{N}}=0.$ (S7)

One finds similar results for the dissipative boundary exchange Hamiltonian
$\hat{H}_{\rm
c}=-\frac{i\gamma}{2}(\hat{\sigma}_{A,1}^{+}\hat{\sigma}_{B,1}^{-}-\hat{\sigma}_{A,1}^{-}\hat{\sigma}_{B,1}^{+})$,
and one finds that the action of the bulk Hamiltonian $\hat{H}_{\rm
XXZ}^{(AB)}=\hat{H}_{\rm XXZ}^{(A)}-\hat{H}_{\rm XXZ}^{(B)}$ on each adjacent
pair $|s_{i}s^{\prime}_{i+1}\rangle$ is either zero or it generates a single
defect.

The creation of these defects and the zero-energy eigenstate condition
requiring that all defects destructively interfere, taken together, can be
shown to constrain the form of the MPS matrices $\hat{A}_{s}$ in _exactly_
exactly the same way as the so-called “isolated defect operator” method
employed in Refs. [4, 5]. The CQA dark state condition
$\hat{L}^{(AB)}|\psi_{\rm CQA}\rangle=0$ is identical to the dissipative
boundary condition found in Ref. [4], but it provides a more physically
transparent reason for the boundary condition. Following Ref. [4], we choose a
matrix representation $\mathcal{V}$ in the semi-infinite space
$\\{|k\rangle:k=0,1,\dots\\}$ with the three (translationally-invariant)
matrices

$\hat{A}_{T}=\sum_{k=0}^{\infty}a_{k}\ket{k}\bra{k},\quad\hat{A}_{1}=\sum_{k=0}^{\infty}b_{k}\ket{k+1}\bra{k},\quad\hat{A}_{0}=\sum_{k=0}^{\infty}c_{k}\ket{k}\bra{k+1},$
(S8)

with coefficients $a_{k}$, $b_{k}$, and $c_{k}$ determined by the recursion
relations coming from the Hamiltonian eigenstate condition in the bulk. Note
that while these matrices are formally defined in an infinite-dimensional
space, for any finite length chain $N<\infty$ we can cut the space off at
level $k=N$ because the left vector $\mathbf{v_{L}}=|0\rangle$ is localized to
$k=0$, hence $\bra{0}(\hat{A}_{0})^{N}\ket{k>N}=0$. Working in units of the
hopping rate $J\equiv 1$, we obtain two recursion relations that determine
$a_{k}$ and the product $b_{k}c_{k}$:

$\displaystyle a_{k+1}-2\Delta a_{k}+a_{k-1}=0,$ (S9) $\displaystyle
b_{k}c_{k}-b_{k-1}c_{k-1}=a_{k}(\Delta a_{k}-a_{k-1}).$ (S10)

The initial conditions for $a_{0}$, $a_{1}$, and $b_{0}c_{0}$ are determined
by the dissipative boundary. The dark state condition Eq. (S6) is
automatically satisfied by picking the left vector $\mathbf{v_{L}}=\ket{0}$
and $\hat{A}_{1}$ a raising operator, but the exchange Hamiltonian
$\hat{H}_{\rm c}$ generates extra defect terms in the Hamiltonian dark state
condition for site 1 (compared to sites only acted upon by the bulk
Hamiltonian) which, in turn, impose the initial conditions

$\displaystyle a_{0}=1,\quad a_{1}=\Delta+\frac{i\gamma}{2},\quad
c_{0}b_{0}=\frac{i\gamma}{2}.$ (S11)

In the easy-axis regime $\Delta/J<1$ and with the above initial conditions,
the recursion relations for $a_{k}$ and the product $b_{k}c_{k}$ have the
exact analytic solution

$\displaystyle
a_{k}=\cos(k\eta)+\frac{i\gamma}{2}\frac{\sin(k\eta)}{\sin\eta},$ (S12)
$\displaystyle
b_{k}c_{k}=\frac{1}{2}\sin[(k+1)\eta]\left[\frac{i\gamma}{\sin\eta}\cos(k\eta)-(1+\frac{\gamma^{2}}{4\sin^{2}\eta})\sin(k\eta)\right],$
(S13)

where $\cos\eta\equiv\Delta/J$. The dissipation dark state condition
$\hat{L}^{(AB)}|\psi_{\rm CQA}\rangle=0$ is satisfied, and the Hamiltonian
zero-energy eigenstate condition is satisfied on all sites except the
coherently driven site $N$.

To fully satisfy the energy eigenstate condition, we must find the
coefficients $\alpha_{k}$ of the right vector
$\mathbf{v_{R}}=\sum_{k=0}^{N}\alpha_{k}\ket{k}$. The interference between the
defects generated by the bulk Hamiltonian $\hat{H}_{\rm XXZ}^{(AB)}$ and the
coherent drive $\hat{H}_{\rm d}^{(AB)}$ yields the recursion relation

$c_{k}\alpha_{k+1}+\frac{\sqrt{2}}{\Omega}(\Delta
a_{k}-a_{k+1})\alpha_{k}-b_{k-1}\alpha_{k-1}=0,\quad{\rm for}\;k>0.$ (S14)

We take the initial conditions for the recursion relation to be $\alpha_{0}=1$
and $\alpha_{-1}=0$. Then we have a well-defined recursion relation, but one
that depends on the coefficients $b_{k}$ and $c_{k}$ separately. Notice that
the bulk recursion relation Eq. (S10) only determines the product
$b_{k}c_{k}$, which reflects a gauge freedom in the MPS representation,
whereas Eq. (S14) requires $b_{k}$ and $c_{k}$ separately. In fact, we can
always perform a similarity transformation $\hat{V}$ so that
$\tilde{\hat{A}}_{s}=\hat{V}\hat{A}_{s}\hat{V}^{-1}$, and
$\tilde{\mathbf{v}}=\hat{V}\mathbf{v}$. While in general, $\hat{V}$ can be any
invertible matrix, here we consider the restricted similarity transformation
$\hat{V}=\sum_{k=0}^{N}v_{k}\ket{k}\bra{k}$ . This specific transformation
leaves the MPS matrices still in tridiagonal form Eq. S8. The transformation
$\hat{V}$ leaves the physical state $\ket{\psi_{\rm CQA}}$ invariant. Yet it
transforms the auxiliary operator and right vector coefficients to

$\tilde{c}_{k}=c_{k}\frac{v_{k+1}}{v_{k}},\quad\tilde{b}_{k}=b_{k}\frac{v_{k}}{v_{k+1}},\quad\tilde{a}_{k}=a_{k},\quad\tilde{\alpha}_{k}=\frac{1}{v_{k}}\alpha_{k},$
(S15)

while keeping the recursion relations Eqs. (S9), (S10) and (S14) valid. As we
will see, this gauge freedom allows us to find a particularly simple solution
to the recursion relation Eq. (S14) in the weak dissipation limit.

###  I.2 Proof of steady state uniqueness

In this section we briefly present the proof for the uniqueness of the NESS
using Frigerio’s first theorem [6, 7]. The theorem states that if the master
equation has a full rank steady state $\hat{\rho}_{\rm ss}$, and if the set
$\\{\hat{L}_{j},\hat{L}_{j}^{\dagger},\hat{H}\\}$ of jump operators
$\\{\hat{L}_{j}\\}$ and Hamiltonian $\hat{H}$ generate all operators in the
Hilbert space, then the steady state is unique.

First we must compute the exact steady state density matrix $\hat{\rho}_{\rm
ss}$ of the original master equation (S2). It is found from $|\psi_{\rm
CQA}\rangle$ by tracing out the absorber

$\hat{\rho}_{\rm ss}=\Tr_{B}(\outerproduct{\psi_{\rm CQA}}{\psi_{\rm
CQA}})=\Tr_{B}(\hat{U}_{B}\outerproduct{\psi_{\rm CQA}}{\psi_{\rm
CQA}}\hat{U}_{B}^{\dagger}).$ (S16)

Here it is convenient to perform a spin-flip of all spins in the absorber
system by acting with the spin-flip unitary
$\hat{U}_{B}=\prod_{i}\hat{\sigma}_{B,i}^{x}$. As we show next, this spin flip
of the absorber allows us to write the steady state $\hat{\rho}_{\rm ss}$ as a
Cholesky decomposition, $\hat{\rho}_{\rm ss}=\hat{\Psi}\hat{\Psi}^{\dagger}$
where $\hat{\Psi}$ is a lower triangular matrix. Using this decomposition, we
then prove that the steady state is full rank.

We define the matrix
$\Psi_{s_{a},s_{b}}=\bra{s_{a},s_{b}}\hat{U}_{B}\ket{\psi_{\rm CQA}}$ with
matrix elements indexed by $s_{a}=(s_{a,1},s_{a,2},\dots,s_{a,N})$ and
$s_{b}=(s_{b,1},s_{b,2},\dots,s_{b,N})$ where
$s_{a/b,i}\in\\{{\uparrow},{\downarrow}\\}$ are the spin states of the system
A and absorber B, respectively. With the definition we find the matrix
elements of $\hat{\rho}_{\rm ss}$ to be

$\displaystyle\bra{s_{a}}\hat{\rho}_{\rm
ss}\ket{s_{a}^{\prime}}=\sum_{s_{b}}\bra{s_{a},s_{b}}\hat{U}_{B}\ket{\psi_{\rm
CQA}}$ $\displaystyle\bra{\psi_{\rm
CQA}}\hat{U}^{\dagger}_{B}\ket{s_{a}^{\prime},s_{b}}=\sum_{s_{b}}\Psi_{s_{a},s_{b}}\Psi_{s_{b},s_{a}^{\prime}}^{*}$
(S17) $\displaystyle\implies\hat{\rho}_{\rm ss}$
$\displaystyle=\hat{\Psi}\hat{\Psi}^{\dagger}.$ (S18)

The matrix elements of $\hat{\Psi}$ can be found explicitly by substituting
Eq. (S5) into $\Psi_{s_{a},s_{b}}=\bra{s_{a},s_{b}}\hat{U}_{B}\ket{\psi_{\rm
CQA}}$ to get

$\hat{\Psi}=\sum_{s_{i}\in\\{\pm,0\\}}(\mathbf{v_{L}^{\dagger}}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}\mathbf{v_{R}})\times\hat{\sigma}_{1}^{s_{1}}\otimes\cdots\otimes\hat{\sigma}_{N}^{s_{N}},$
(S19)

where the $2$ by $2$ matrices are
$\hat{\sigma}^{+}=\outerproduct{\uparrow}{\downarrow}$,
$\hat{\sigma}^{-}=\outerproduct{\downarrow}{\uparrow}$, and
$\hat{\sigma}^{0}=\frac{1}{\sqrt{2}}(\outerproduct{\uparrow}{\uparrow}+\outerproduct{\downarrow}{\downarrow})$.
These directly correspond to the double chain dimer state $|1_{i}\rangle$,
$|0_{i}\rangle$, and $|T_{i}\rangle$, respectively. Likewise, the relabeling
of the spin indices applies to the $\hat{A}_{s}$ matrices: $s\in\\{0,1,T\\}\to
s\in\\{-,+,0\\}$ [cf. Eq. (S8)]. To show that $\hat{\Psi}$ is a lower
triangular matrix, recall that $\mathbf{v_{L}^{\dagger}}\hat{A}_{+}=0$ and
$\mathbf{v_{L}^{\dagger}}\hat{A}_{0}=\mathbf{v_{L}^{\dagger}}$. Then, we see
that Eq. (S19) _cannot_ have terms where a $\hat{\sigma}_{i}^{+}$ precedes the
first $\sigma_{j>i}^{-}$. Hence, $\Psi$ is a lower triangular matrix;
therefore, it is the Cholesky decomposition of $\hat{\rho}_{\rm ss}$.

Equipped with the NESS $\hat{\rho}_{\rm ss}=\hat{\Psi}\hat{\Psi}^{\dagger}$,
we can now prove that it is full rank (i.e., positive definite). The matrix
elements of $\hat{\Psi}$ are directly connected to the MPS solution via

$\bra{\underline{s}^{\prime}}\hat{\Psi}\ket{\underline{s}}=\mathbf{v_{L}^{\dagger}}\hat{A}_{s_{1}^{\prime}-s_{1}}\hat{A}_{s_{2}^{\prime}-s_{2}}\cdots\hat{A}_{s_{N}^{\prime}-s_{N}}\mathbf{v_{R}}.$
(S20)

Here $\underline{s}=(s_{1},\dots,s_{N})$ is a configuration vector in the spin
basis, with each $s_{i}\in\\{0,1\\}$. We see by direct calculation that the
diagonal elements of $\hat{\Psi}$ are
$\bra{\underline{\nu}}\hat{\Psi}\ket{\underline{\nu}}=\mathbf{v_{L}^{\dagger}}(\hat{A}_{0})^{N}\mathbf{v_{R}}=a_{0}^{N}\alpha_{0}=1$.
Since $\hat{\Psi}$ is lower triangular, its eigenvalues are all 1, hence it is
full rank. Therefore, for an arbitrary state $|\phi\rangle\neq 0$, we have

$\displaystyle\langle\phi|\hat{\rho}_{\rm
ss}|\phi\rangle=\langle\phi|\hat{\Psi}\hat{\Psi}^{\dagger}|\phi\rangle=||\hat{\Psi}^{\dagger}|\phi\rangle||^{2}>0,$
(S21)

where the strict inequality follows from $\hat{\Psi}^{\dagger}$ being full
rank. Therefore, we have shown $\hat{\rho}_{\rm ss}$ to be full rank.

Next, we show that $\\{\hat{H}_{\rm
XXZ},\hat{\sigma}_{1}^{-},\hat{\sigma}_{1}^{+}\\}$ generates the whole algebra
of the Hilbert space. Note that while the full Hamiltonian includes the
coherent driving $\hat{H}=\hat{H}_{\rm XXZ}+\hat{H}_{\rm d}$, the driving term
is redundant for this proof. First, observe that

$\hat{\sigma}_{2}^{-}=\frac{1}{4}\hat{\sigma}_{1}^{z}[\hat{\sigma}_{1}^{-},[\hat{H}_{\rm
XXZ},\hat{\sigma}_{1}^{z}]],$ (S22)

where we assumed $J=1$ for simplicity. This also generates
$\hat{\sigma}_{2}^{+}$ by conjugation. In general, for all sites $j>1$, we can
inductively apply

$\hat{\sigma}_{j+1}^{-}=\frac{1}{4}\hat{\sigma}_{j}^{z}[\hat{\sigma}_{j}^{-},[\hat{H}_{\rm
XXZ},\hat{\sigma}_{j}^{z}]]-\hat{\sigma}_{j-1}^{-}.$ (S23)

By conjugation this also generates $\hat{\sigma}_{j+1}^{+}$. Since
$\sigma_{j}^{-}$, $\sigma_{j}^{+}$, and their products generate the full
algebra on each site, this procedure generates the full Hilbert space algebra.

With the two conditions of Frigerio’s first theorem satisfied, we have proved
that $\hat{\rho}_{\rm ss}=\hat{\Psi}\hat{\Psi}^{\dagger}/{\rm
Tr}(\hat{\Psi}\hat{\Psi}^{\dagger})$ is the unique steady state of Eq. (S2).
There is one subtlety here: in the absence of a coherent drive, $\Omega=0$,
the steady state is not full rank, but rather the completely polarized pure
state $|{\downarrow\cdots\downarrow}\rangle$. Generically, for any $\Omega\neq
0$, the result holds. Finally, we point out that this proof of uniqueness
extends to our previous work of the boundary driven XX chain [8].

As a final note, it has been found that the “isolated defect operator” method
yields exact solutions to other dissipative XXZ chain models in the form of
Cholesky decompositions of the steady states [4, 9, 5]. Due to the
similarities between the CQA approach and the previous exact solutions, both
in the form of the matrix product ansatz and in the form of the solution as a
Cholesky decomposition, we suspect that these previously solved models – and
others amenable to the isolated defect operator method – could be solved using
the CQA approach. Indeed, below in Section VI we construct the CQA system and
find the doubled-system pure state solution for the model first solved in Ref.
[4].

###  I.3 Observables

In this section, we describe how to compute the relevant observables
(magnetization and current) from the MPS solution Eq. S5.

The main quantity of interest, the average magnetization, is defined by

$\displaystyle\langle{m}\rangle\equiv\frac{1}{N}\sum_{j=1}^{N}\Tr(\hat{\sigma}_{j}^{z}\hat{\rho}_{\rm
ss})=\frac{1}{N}\sum_{j=1}^{N}\frac{\bra{\psi_{\rm
CQA}}\hat{\sigma}_{j}^{z}\ket{\psi_{\rm CQA}}}{\bra{\psi_{\rm
CQA}}\ket{\psi_{\rm CQA}}},$ (S24)

where $\hat{\sigma}_{j}^{z}$ in the final expression can be either a system A
operator or an absorber B operator. Here we compute it in terms of the exact
solution MPS state Eq. (S5).

The local spin operators $\hat{\sigma}_{j}^{z}$ have matrix elements
$\bra{0_{j}}\hat{\sigma}_{j}^{z}\ket{0_{j}}=-1$,
$\bra{T_{j}}\hat{\sigma}_{j}^{z}\ket{T_{j}}=0$, and
$\bra{1_{j}}\hat{\sigma}_{j}^{z}\ket{1_{j}}=1$ in the triplet spin basis. All
of the other terms vanish. The total magnetization is, thus, diagonal in this
basis:

$m_{\\{s\\}}=\sum_{j=1}^{N}\bra{\\{s^{\prime}\\}}\hat{\sigma}_{j}^{z}\ket{\\{s\\}}=\delta_{\\{s\\},\\{s^{\prime}\\}}(m_{+}-m_{-}).$
(S25)

Here $\ket{\\{s\\}}\equiv\ket{s_{1},s_{2},\cdots,s_{N}}$ is the spin
configuration and we denote the number of $1$s and $0$s by
$m_{+}=\sum_{j=1}^{N}|\bra{1_{j}}\ket{\\{s\\}}|^{2}$ and
$m_{-}=\sum_{j=1}^{N}|\bra{0_{j}}\ket{\\{s\\}}|^{2}$, respectively. We relate
these quantities to the coefficients of the MPS matrices and right vector as
follows. As $\hat{A}_{0(1)}$ always move the level in the auxiliary space by
$\pm 1$, the total magnetization is mapped to the distance from 0 in the
auxiliary space after taking $N$ steps (e.g., the completely polarized
configuration $\\{{\downarrow\cdots\downarrow}\\}$ with total magnetization
$-N$ connects $\mathbf{v_{L}^{\dagger}}=\langle 0|$ on the left to the
$\alpha_{N}|N\rangle$ component of the right vector). The coefficient of
configuration $\\{s\\}$ with $m_{\\{s\\}}=k$ is therefore
$\mathbf{v_{L}^{\dagger}}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}|k\rangle\langle
k|\mathbf{v_{R}}=\bra{0}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}\ket{k}\alpha_{k}$.
The magnetization density can be given compactly as

$\displaystyle\langle\hat{m}\rangle$
$\displaystyle=-\frac{1}{N}\frac{\sum_{k=0}^{N}\sum_{\\{s\\}}k|\bra{0}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}\ket{k}|^{2}|\alpha_{k}|^{2}}{\sum_{k=0}^{N}\sum_{\\{s\\}}|\bra{0}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}\ket{k}|^{2}|\alpha_{k}|^{2}}=-\frac{1}{N}\frac{\sum_{k=0}^{N}k\bra{0,0}\hat{T}^{N}\ket{k,k}|\alpha_{k}|^{2}}{\sum_{k=0}^{N}\bra{0,0}\hat{T}^{N}\ket{k,k}|\alpha_{k}|^{2}}.$
(S26)

Here we defined the vectorized auxiliary space $\mathcal{V}\otimes\mathcal{V}$
by mapping $\ket{k}\bra{k}\to\bra{k,k}$, and we introduced the transfer matrix
$\hat{T}=\sum_{s}(\hat{A}_{s}\otimes\hat{A}_{s}^{*})$.

Because $\hat{T}$ always keeps the state on the diagonal $\ket{j,j}$, we can
further reduce the vectorized auxiliary space to its diagonal subspace, with
states denoted by $\ket{\mathbf{j}}\equiv\ket{j,j}$. By doing so, we find the
magnetization density to be

$\langle\hat{m}\rangle=-\frac{1}{N}\frac{\sum_{k=0}^{N}k\bra{\mathbf{0}}\hat{T}^{N}\ket{\mathbf{k}}|\alpha_{k}|^{2}}{\sum_{k=0}^{N}\bra{\mathbf{0}}\hat{T}^{N}\ket{\mathbf{k}}|\alpha_{k}|^{2}}=-\sum_{k=0}^{N}\frac{k}{N}p_{k},$
(S27)

where the transfer matrix $\hat{T}$ is defined by

$\displaystyle\hat{T}\equiv\sum_{k=0}^{N}|a_{k}|^{2}\ket{\mathbf{k}}\bra{\mathbf{k}}+|b_{k}|^{2}\ket{\mathbf{k+1}}\bra{\mathbf{k}}+|c_{k}|^{2}\ket{\mathbf{k}}\bra{\mathbf{k+1}}.$
(S28)

Eq. S27 expresses the magnetization as a weighted average over the probability
distribution
$p_{k}=\bra{\mathbf{0}}\hat{T}^{N}\ket{\mathbf{k}}|\alpha_{k}|^{2}/(\sum_{k=0}^{N}\bra{\mathbf{0}}\hat{T}^{N}\ket{\mathbf{k}}|\alpha_{k}|^{2})$.
Given the exact form of $a_{k}$, $b_{k}$, and $c_{k}$ found from the solutions
Eqs. (S12) and (S13), and by solving Eq. (S14), we see that the statistical
weight comes from two parts: the boundary contribution $|\alpha_{k}|^{2}$,
which represents the interplay between drive and bulk properties, and the
transfer matrix $\bra{\mathbf{0}}\hat{T}^{N}\ket{\mathbf{k}}$, which reflects
the bulk dynamics and the dissipation. As we show in Section II.2, in the weak
dissipation limit $\gamma/J\ll 1$, there is an MPS gauge in which the
distribution $\\{p_{k}\\}$ takes a relatively simple form.

Similarly, there are other quantities can be defined in the vectorized
auxiliary space $\mathcal{V}\otimes\mathcal{V}$. Here we omit the detailed
calculation and directly state the results. By defining the partition function
$Z_{N}=\sum_{k=0}^{N}\bra{\mathbf{0}}\hat{T}^{N}\ket{\mathbf{k}}|\alpha_{k}|^{2}$,
the spin-spin correlation is found to be
$\langle\hat{\sigma}^{z}_{a}\hat{\sigma}^{z}_{b}\rangle=\frac{1}{Z_{N}}\sum_{k=0}^{N}\bra{\mathbf{0}}\hat{T}^{a-1}\hat{Z}\hat{T}^{b-a}\hat{Z}\hat{T}^{N-b}\ket{\mathbf{k}}|\alpha_{k}|^{2}$
for $b>a$ and where
$\hat{Z}=\sum_{k}|b_{k}|^{2}|\mathbf{k+1}\rangle\langle\mathbf{k}|-|c_{k}|^{2}|\mathbf{k}\rangle\langle\mathbf{k+1}|$.
Similarly, the spin current is $\langle
j\rangle=\frac{\gamma}{J}\frac{Z_{N-1}}{Z_{N}}$. We see that the complexity of
calculating observable quantities is equivalent to diagonalizing the $\hat{T}$
matrix, with a total time complexity of $\mathcal{O}(N^{4})$. This allows us
to readily calculate steady state properties for systems with up to thousands
of spins.

##  II Phase transition

###  II.1 Quantum to classical mapping

In this section we show that the transfer matrix of Eq. (S28) can be
interpreted as the transition matrix of a classical random walk model on a
one-dimensional semi-infinite lattice. Specifically, we can interpret the
matrix elements $\bra{\mathbf{k}}(\hat{T}^{\dagger}){}^{N}\ket{\mathbf{0}}$ as
being the probability that a particle initialized on site 0 ends on site $k$
after $N$ time steps (with each time step corresponding to a random hop
described by $T^{\dagger}$). Because the coefficients
$\bra{\mathbf{k}}(\hat{T}^{\dagger}){}^{N}\ket{\mathbf{0}}$ directly help
determine the average magnetization via Eq. (S27), we can then reinterpret the
magnetization as the mean distance that a particle initialized on site 0
travels in $N$ time steps, weighted by a post-selection distribution
$|\alpha_{k}|^{2}$.

To map the problem of calculating
$\bra{\mathbf{k}}(\hat{T}^{\dagger}){}^{N}\ket{\mathbf{0}}$ to a classical
random walk on the lattice one-dimensional semi-infinite lattice comprising
$\\{\ket{\mathbf{k}}\\!,\,\mathbf{k}=0,1,\cdots\\}$, we can define a classical
Markov process with the master equation:

$p(k,t+1)=|a_{k}|^{2}p(k,t)+|c_{k-1}|^{2}p(k-1,t)+|b_{k}|^{2}p(k+1,t).$ (S29)

Heuristically, at each discrete time step $t$ a particle on site $k$ of the
auxiliary lattice can either hop to the left with with probability
$|b_{k-1}|^{2}$, hop right with probability $|c_{k}|^{2}$, or remain on the
same site with probability $|a_{k}|^{2}$. The probability distribution
$p(k,N)$, after taking $N$ steps from the initial condition
$p(k,0)=\delta_{k,0}$, encodes the transfer matrix elements
$\bra{\mathbf{k}}(\hat{T}^{\dagger}){}^{N}\ket{\mathbf{0}}=p(k,N)$.

For Eq. (S29) to meaningfully represent a stochastic process, we must enforce
probability conservation of the coefficients [cf. Eq. (S8)]:
$|a_{k}|^{2}+|b_{k-1}|^{2}+|c_{k}|^{2}=1$. This can always be enforced using
the gauge freedom of the solution Eq. S15. First notice from the recursion
relations Eqs. (S9) and (S10) that we can freely scale all coefficients by a
constant factor $C$:

$\displaystyle\tilde{a}_{k}=Ca_{k},\quad\tilde{b}_{k}=Cb_{k},\quad\tilde{c}_{k}=Cc_{k}.$
(S30)

We have the gauge freedom in Eq. (S15) to fix the ratio
$|b_{k}|^{2}/|c_{k}|^{2}$ through a judicious choice of $\\{v_{k}\\}$.
Therefore, probability conservation can be enforced by solving

$\displaystyle
C^{2}\left[|a_{k}|^{2}+\frac{\absolutevalue{v_{k-1}}^{2}}{\absolutevalue{v_{k}}^{2}}|b_{k-1}|^{2}+\frac{\absolutevalue{v_{k+1}}^{2}}{\absolutevalue{v_{k}}^{2}}|c_{k}|^{2}\right]=1,\quad\forall
k=0,1,\cdots N-1.$ (S31)

Note that we only need to satisfy Eq. S31 for $k=0,1,\cdots N-1$. This is
because we are only interested in the master equation Eq. S29 subject to the
initial condition $p(k,0)=\delta_{k,0}$ and for times $t\leq N$.

We can solve Eq. S31 by starting from $k=0$ (note that $b_{-1}=0$) and
iteratively solving every equation for the ratio
$\frac{\absolutevalue{v_{k}}^{2}}{\absolutevalue{v_{k+1}}^{2}}$ until $k=N-1$
with

$\frac{\absolutevalue{v_{k+1}}^{2}}{\absolutevalue{v_{k}}^{2}}=(|c_{k}|^{2})^{-1}\quantity(1/C^{2}-|a_{k}|^{2}-\frac{\absolutevalue{v_{k-1}}^{2}}{\absolutevalue{v_{k}}^{2}}|b_{k-1}|^{2}),$
(S32)

with the initial condition
$|v_{1}|^{2}/|v_{0}|^{2}=(|c_{0}|^{2})^{-1}(1/C^{2}-|a_{0}|^{2})$. It is
important of course that the solutions satisfy
$\frac{\absolutevalue{v_{k}}^{2}}{\absolutevalue{v_{k+1}}^{2}}>0$. This can
guaranteed by properly choosing $C$. We start from a specific gauge $c_{k}=1$.
For bounded sequences $|a_{k}|^{2}<a^{2}$ and $|b_{k}|^{2}<b^{2}$, we pick a
constant $C$ to be

$1/C^{2}=1+a^{2}+b^{2}.$ (S33)

Such choice of $C$ results in
$|v_{1}|^{2}/|v_{0}|^{2}=(1/C^{2}-|a_{0}|^{2})>(1+b^{2})>1$. Next we show that
$|v_{k+1}|^{2}/|v_{k}|^{2}>1$ using induction. Given
$|v_{k}|^{2}/|v_{k-1}|^{2}>1$, we have

$\frac{\absolutevalue{v_{k+1}}^{2}}{\absolutevalue{v_{k}}^{2}}=\quantity(1/C^{2}-|a_{k}|^{2}-\frac{\absolutevalue{v_{k-1}}^{2}}{\absolutevalue{v_{k}}^{2}}|b_{k-1}|^{2})>\quantity(1+(1-\frac{\absolutevalue{v_{k-1}}^{2}}{\absolutevalue{v_{k}}^{2}})b^{2})>1,$
(S34)

By doing so, we show $|v_{k+1}|^{2}/|v_{k}|^{2}>1$ by the specific
construction of $C$. The net result is that we can always pick a gauge such
that the relevant elements of the $T^{\dagger}$ matrix correspond to a valid
transition matrix of a classical 1D hopping model. While the procedure we have
described is applicable for all parameter choices of our model, there is no
guarantee that this can be done in a manner that yields closed-form analytic
solutions. In the next section, we show that in the special case of weak
dissipation ($\gamma\rightarrow 0$), this mapping to a classical model can be
done analytically.

###  II.2 Weak dissipation limit

In this section we analyze the exact solution in the weak dissipation limit
$\gamma/J\ll 1$ to obtain analytic insight into the phase transition in
$\Omega$. Recall that we have exact analytic expressions for the MPS
coefficients $a_{k}$ and $b_{k}c_{k}$ in the $\Delta/J<1$ regime given by Eqs.
(S12) and (S13), the recursion relation Eq. (S14) for $\alpha_{k}$, and a
gauge freedom in choosing $b_{k}$ and $c_{k}$. In the weak dissipation limit,
a convenient gauge is

$\displaystyle b_{k}=\frac{1}{\sqrt{2}}\sin[(k+1)\eta],$ (S35) $\displaystyle
c_{k}=\frac{1}{\sqrt{2}}\left[\frac{i\gamma}{J\sin\eta}\cos(k\eta)-(1+\frac{\gamma^{2}}{4J^{2}\sin^{2}\eta})\sin(k\eta)\right],$
(S36)

where $\cos\eta=\Delta/J$ and we have restored $J$. Using these solutions and
Eq. (S12) for $a_{k}$, and keeping only the $\gamma$-independent terms, the
recursion relation Eq. (S14) simplifies to

$\alpha_{k+1}-\frac{2J\sin\eta}{\Omega}\alpha_{k}+\alpha_{k-1}=0$ (S37)

with initial conditions $\alpha_{0}=1$ and $\alpha_{1}=\sin\eta/\Omega$. Here
we assumed $\sin k\eta\neq 0$. This is generically true except at the set of
special points discussed in Sec. III.

The solution to Eq. (S37) is a Chebyshev polynomial of the first kind
$\alpha_{k}=T_{k}(\frac{J}{\Omega}\sin\eta)$. Across the critical driving
strength $\Omega_{c}\equiv J\sin\eta=\sqrt{J^{2}-\Delta^{2}}$, there is a
change in the functional form of the solution:

$|\alpha_{k}|^{2}=\begin{cases}\cos^{2}k\theta,\leavevmode\nobreak\
\leavevmode\nobreak\
\,\theta=\arccos(\Omega_{c}/\Omega)&\frac{\Omega}{\Omega_{c}}>1,\\\
\cosh^{2}k\theta,\leavevmode\nobreak\ \theta={\rm
arccosh}(\Omega_{c}/\Omega)&\frac{\Omega}{\Omega_{c}}<1.\end{cases}$ (S38)

Above the critical drive $\Omega/\Omega_{c}>1$, we have $|\alpha_{k}|^{2}\leq
1$ and oscillatory behaviour as a function of $k$. Below the critical drive
$\Omega/\Omega_{c}<1$, the $|\alpha_{k}|^{2}$ grow exponentially in $k$. This
transition in the behavior of $|\alpha_{k}|^{2}$ leads to non-analytic
behavior in the spin chain magnetization and the free energy of the classical
model later shown in Section II.3.

The transfer matrix has a particularly simple form in the weak dissipation
limit; this allows us to derive insights about the phase transition from an
effective classical random walk model. Keeping only the zeroth order in
$\gamma$ terms in $a_{k}$, $b_{k}$, and $c_{k}$ [except for $c_{0}$ whose
leading order is $\mathcal{O}(\gamma^{2})$, see discussion at the end of Sec.
II.2], we see that the classical master equation (S29) in this limit yields
simple reciprocal hopping rates that vary quasiperiodically along the
auxiliary lattice:
$|b_{k-1}|^{2}=|c_{k}|^{2}=w_{k}\equiv\frac{1}{2}\sin^{2}k\eta$ and
$|a_{k}|^{2}=1-2w_{k}$. Hence the classical master equation is

$p(k,t+1)=(1-2w_{k})p(k,t)+w_{k-1}p(k-1,t)+w_{k+1}p(k+1,t).$ (S39)

In general the effective wavelength $2\pi/\eta$ of the modulation of the
$w_{k}$ will not be commensurate with the auxiliary lattice, so the classical
master equation corresponds to a model of quenched disorder. We focus on the
final distribution after $N$ time steps, $f_{N}(k)\equiv
p(k,N)=\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}$. To get
insight into the physics of our classical model, consider first a related,
simplified model where we ignore disorder, and assume
$w_{k}=\overline{w_{k}}=\frac{1}{4}$. We find for this case
$f_{N}(k)=\frac{1}{2^{N}}\exp(-Ng(k/N))$, where $g(\xi)=2\xi\,{\rm
artanh}(\xi)+\ln(1-\xi^{2})$ for the reduced coordinate $\xi=k/N$. As we see
in Fig. S1(a), this result does not agree with the numerically-computed
distribution using the $w_{k}=\frac{1}{2}\sin^{2}k\eta$ of the exact classical
model. Thus, the disorder-free approximation is not sufficient to obtain the
correct functional form of $f_{N}(k)$.

Noting that the weights $w_{k}=\frac{1}{2}\sin^{2}k\eta$ are quasi-periodic in
$k$, we expect that for $\eta\neq\pi\frac{l}{m}$, the arguments of the sin
function, $(k\eta\,{\rm mod}\,\pi)$ are distributed uniformly in the interval
$[0,\pi)$. Thus, we make the approximation that the weights are distributed in
the interval $[0,1/2)$ with probability distribution

$\displaystyle P(w)=\frac{2}{\pi\sqrt{2w(1-2w)}}.$ (S40)

This classical single-particle model is similar to the Sinai random walk [10],
but with a crucial difference: there is a non-zero probability $1-2w_{k}$ to
remain on site $k$. As we see in Fig. S1(a), this model estimates the true
transfer matrix elements $f_{N}(\xi)$ extremely well. We seek for a specific
form taking out the system size contribution

$g(\xi)\equiv-\frac{1}{N}\ln\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}.$
(S41)

We numerically find sub-diffusive scaling $\langle k^{2}\rangle\sim N^{2/3}$
and $g(\xi)\sim\xi^{3/2}$, as shown in Fig. S1(b). In contrast, if we took the
weights to be uniform, $w_{k}=\frac{1}{4}$, then we would find a diffusive
behaviour, i.e. $g(\xi)\sim\xi^{2}$.

Finally we turn to the behavior of the transfer matrix elements
$\bra{\mathbf{0}}\hat{T}^{N}\ket{\mathbf{k}}=\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}$
[cf. Eq. (S28)] in thermodynamic limit $N\to\infty$. These are the matrix
elements relevant for analytically describing the phase transition in the
coherent drive $\Omega$. In the weak dissipation limit, we keep only zeroth
order in $\gamma$ terms for all coefficients except $c_{0}$ which would vanish
to $\mathcal{O}(\gamma^{2})$. The transfer matrix is

$\hat{T}^{\dagger}=\sum_{k=0}^{\infty}\cos^{2}k\eta\outerproduct{\mathbf{k}}{\mathbf{k}}+\frac{1}{2}\sin^{2}\eta(k+1)\outerproduct{\mathbf{k}}{\mathbf{k+1}}+\frac{1}{2}\sin^{2}\eta
k\outerproduct{\mathbf{k+1}}{\mathbf{k}}+\frac{\gamma^{2}}{2\sin^{2}\eta}\outerproduct{\mathbf{1}}{\mathbf{0}}$
(S42)

We see that
$\bra{\mathbf{0}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}=\mathcal{O}(1)$ in
$\gamma$ but all other matrix elements are
$\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}=\mathcal{O}(\gamma^{2})$.
When taking the thermodynamic and weak dissipation limits, it is important to
take $N\to\infty$ while keeping $\gamma>0$. Taking this order of limits, each
$k>0$ element of $\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}$ is
at least $\mathcal{O}(N)$ whereas the $k=0$ element is $\mathcal{O}(1)$.
Therefore, in the thermodynamic limit we can safely ignore the single
contribution $\langle\mathbf{0}|(\hat{T}^{\dagger})^{N}|\mathbf{0}\rangle$.
Then we can take $\gamma\to 0$, keeping leading $\mathcal{O}(\gamma^{2})$
terms only.

Figure S1: (a) $g(\xi)$, the log of the $N$-step transfer matrix element from
site $0$ to $k$, c.f. Eq. (S41), is plotted versus $\xi=k/N$. This function
plays the role of a potential in our effective free energy. The quasi-periodic
distribution coming from the exact solution (orange line) is well approximated
by the random environment with the same weight distribution $P(w)$ (blue
line). The results from a simplified model with uniform hopping $w_{k}=1/4$
(green dashed line) yields diffusive behaviour, and does not agree with the
full model. (b) The scaling of the potential $g(\xi)$ with $\xi$ is shown near
the origin. For $\gamma\rightarrow 0$, it follows $\xi^{3/2}$ (dashed line).
As the dissipation $\gamma$ is increased, we find that the power law changes
continuously. (c) Comparing field term $h(\xi)$ for different dissipation.
Finite dissipation leads to different slope. (d) Field term $h(\xi)$ in the
free energy $F(\xi)=g(\xi)+h(\xi)$ are plotted vs. $\xi$ using the exact
solution ($\gamma/J=10^{-3}$). Here $h(\xi)$ changes across the critical drive
from $0$ above $\Omega_{c}$ to a linear field below.

###  II.3 Free Energy and Critical Exponent

In this subsection, we explain how to understand the phase transition from an
effective free energy and derive the critical exponent in the weak dissipation
limit. We can rewrite the average magnetization Eq. (S27) in terms of a
dimensionless “free energy” $F(\xi)$

$\langle\hat{m}\rangle=-\frac{\sum_{\xi}\xi
e^{-NF(\xi)}}{\sum_{\xi}e^{-NF(\xi)}},$ (S43)

where $\xi=k/N\in[0,1]$ is the reduced coordinate. The free energy
$F(\xi)=g(\xi)+h(\xi)$ decomposes into the $\Omega$-independent “potential”
term $g(\xi)$ [Eq. S41] and the $\Omega$-dependent “external field” term

$\displaystyle h(\xi)\equiv-\frac{1}{N}\ln|\alpha(\xi)|^{2}.$ (S44)

In the thermodynamic limit, the mean magnetization is found by minimizing the
free energy $F(\xi)$ over $\xi$. Using the weak dissipation expression Eq. S38
for $|\alpha(\xi)|^{2}$, we find that the field term is continuous but non-
analytical in the thermodynamic limit due to the discontinuous derivative
$\partial_{\Omega}h(\xi)$ at $\Omega_{c}$:

$\lim_{N\to\infty}h(\xi)=-\lim_{N\to\infty}\frac{1}{N}\ln|\alpha_{\xi}|^{2}=\begin{cases}-\lim_{N\to\infty}\frac{1}{N}\ln\cos^{2}N\theta\xi=0,&\Omega>\Omega_{c},\\\
-\lim_{N\to\infty}\frac{1}{N}\ln\cosh^{2}N\theta\xi=-2\theta\xi,&\Omega<\Omega_{c}.\end{cases}$
(S45)

The non-analytical behavior of the field $h(\xi)$ explains the mechanism
behind the phase transition. In the ordered phase $\Omega>\Omega_{c}$, the
field term is turned off $h(\xi)=0$. The free energy $F(\xi)=g(\xi)$ has a
global minimum at $\xi^{*}=0$. Whereas when $\Omega<\Omega_{c}$, the free
energy takes the form of $F(\xi)=g(\xi)-2\theta\xi$, with a global minimum
$\xi^{*}>0$.

Note that in our chosen gauge, the mapping to a classical stochastic model is
only valid for $\gamma\rightarrow 0$. We stress however that this is just a
question of interpretation. The expressions for $\langle m\rangle$ and our
effective free energy remain valid no matter the value of $\gamma$. The
behavior of the external field $h(\xi)$ is numerically computed for
$\gamma/J=10^{-2}$, using the exact solution in Fig. S1(d) for $N=1000$. We
see that $h(\xi)$ suddenly switches on for $\Omega<\Omega_{c}$.

We can obtain a critical exponent $\beta$ for the scaling of magnetization
with the drive below the critical point,
$\langle\hat{m}\rangle\sim\delta^{\beta}$, where
$\delta=\frac{\Omega_{c}-\Omega}{\Omega_{c}}\ll 1$. Near the transition, we
assume the potential is a power law $g(\xi)=g_{0}+g_{1}\xi^{m}$ for some
$m>0$. The global minimum of the free energy is, thus,

$\displaystyle\delta
F(\xi^{*})=0\implies\xi^{*}=\xi_{0}\theta^{\frac{1}{(m-1)}}.$ (S46)

Recall that $\cos\theta=\Omega_{c}/\Omega$ for $\Omega>\Omega_{c}$, and
$\cosh\theta=\Omega_{c}/\Omega$ for $\Omega<\Omega_{c}$. Further, as discussed
in Section II.2, we numerically found that $g(\xi)\sim\xi^{3/2}$ in the weak
dissipation limit, implying that $m=3/2$. Therefore, we predict a critical
exponent of $\beta=1$ in the $\gamma/J\to 0$ limit.

Away from the weak dissipation limit, we numerically find that $h(\xi)$
undergoes a similar transition as in Eq. S45. However, both $h(\xi)$ and
$g(\xi)$ show continues dependence on $\gamma$ for $\Omega<\Omega_{c}$ [see
Fig. S1(b) and Fig. S1(c)]. As we will see in the next section, the phase
transition still occurs at the same point even away from weak dissipation
limit, but the critical exponent $\beta$ varies continuously with $\Omega$.

###  II.4 Finite size scaling

In this section, we perform numerical analysis of the phase transition away
from the weak dissipation limit. First, we show that the critical drive
$\Omega_{c}=\sqrt{J^{2}-\Delta^{2}}$, derived analytically in the weak
dissipation limit, remains valid even away from that limit. Numerically, we
diagnose the critical drive by the discontinuity in the magnetization
susceptibility $\partial_{\Omega}\langle\hat{m}\rangle$ (in the weak
dissipation limit, this corresponds to the sharp change in the free energy
external field $h(\xi)$ in Eq. (S45)). Numerically, we observe in Fig. S2 that
the susceptibility for different system sizes cross each other all at one
point. Moreover, this crossing agrees extremely well with the analytically
predicted $\Omega_{c}=\sqrt{J^{2}-\Delta^{2}}$, even a way from the weak
dissipation limit.

Next, we perform a finite size scaling analysis and observe scaling collapse
near the critical drive $\delta=\frac{\Omega_{c}-\Omega}{\Omega_{c}}\ll 1$ to
a single universal scaling function:

$\langle m\rangle=N^{-a}\mathcal{M}(\delta N^{b})$ (S47)

for scaling exponents $a$ and $b$. In particular, near the transition, we find
the scaling function $\langle\hat{m}\rangle\sim\delta^{\beta}N^{\beta b-a}$
where the fitted scaling exponents $a$ and $b$ are related to the observed
critical exponent $\beta=\frac{a}{b}$. In Fig. S3(a) we reproduce Fig. 2(c) of
the main text _without_ rescaling, which verifies that the scaling collapse to
the above universal function is a real effect. Furthermore, we see that there
is clear scaling collapse to a single scaling function even as the dissipation
$\gamma/J$ is varied. Despite the $\gamma$-dependence of the critical exponent
$\beta$, we still find universal scaling collapse, as shown in Fig. S3(b) and
(c), and that $\beta=\frac{a}{b}$ holds.

In Fig. S3 we focused on dissipation rates $\gamma/J\sim 1$ instead of in the
weak dissipation limit. There is a subtlety in analyzing the weak dissipation
limit related to the weak bond in the transfer matrix
$\bra{\mathbf{1}}\hat{T}^{\dagger}\ket{\mathbf{0}}=\frac{\gamma^{2}}{\sin^{2}\eta}$
discussed in Eq. S42. We argued that in the thermodynamic limit, taking care
order the limits $N\to\infty$ before $\gamma\to 0$, this weak bond does not
matter as the contribution to the magnetization of
$\bra{\mathbf{0}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}/\sum_{k}^{N}\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}\sim\mathcal{O}(N^{-1})$
vanishes. For finite $N$ and sufficiently weak $\gamma/J\ll 1$, this assertion
is no longer valid. In fact, due to the $\mathcal{O}(\gamma^{2})$ suppression
of $\bra{\mathbf{1}}\hat{T}^{\dagger}\ket{\mathbf{0}}$ relative to
$\bra{\mathbf{0}}\hat{T}^{\dagger}\ket{\mathbf{0}}$, we observe a kind of
“condensation” physics, wherein there is an unusually large occupation of the
fully depolarized $\langle\hat{m}\rangle=0$ state. In this regime, we
approximate
$\bra{\mathbf{0}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}=1+\mathcal{O}(\gamma^{2})\approx
1$. Thus, the magnetization is

$\langle\hat{m}\rangle=-\frac{1}{N}\frac{\sum_{k>0}^{N}k\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}|\alpha_{k}|^{2}}{1+\sum_{k>0}^{N}\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}|\alpha_{k}|^{2}}=-\frac{\gamma^{2}\int\differential\xi\exp(-NF(\xi))\xi}{1+\gamma^{2}\int\differential\xi\exp(-NF(\xi))},$
(S48)

where the condensation effects appear as the $\mathcal{O}(1)$ term in the
normalization. This additional term is significant when
$\gamma^{2}\int\differential\xi\exp(-NF(\xi))\approx\gamma^{2}N^{2/3}\ll 1$
and vanishes by taking the thermodynamic limit first. However, we see that for
small $\gamma$, we require very large system sizes to eliminate this effect.
In Fig. S4 we see strong finite-size effects for $\delta\ll 1$ due to this
behavior. Fig. S4(a) is a plot of the approximate magnetization Eq. (S48) for
various system sizes and Fig. S4(b) and (c) show the magnetization calculated
numerically from the exact solution. Nevertheless, our approximation captures
the essential physics and we see the finite size effect vanishes for large
system size in Fig. S4(c). Moreover, this allows us to extract the critical
exponent $\beta$ by fitting the scaling of
$\langle\hat{m}\rangle\sim\delta^{\beta}$ in the range of $\delta$
sufficiently large to avoid the effects of the condensation physics [see Fig.
S4(b) and (c)]. We confirm that as $\gamma/J\to 0$, we recover the predicted
critical exponent $\beta=1$.

Figure S2: The magnetic susceptibility
$\partial_{\Omega}\langle\hat{m}\rangle$ is plotted as a function of $\Omega$
near the critical drive $\Omega_{c}$. Here we use the analytic prediction for
$\Omega_{c}=\sqrt{J^{2}-\Delta^{2}}$. Each plot is for a different (labeled)
value of dissipation $\gamma/J$. We find a small relative error on the order
$10^{-3}$ between the analytically predicted $\Omega_{c}$ and the crossing
point of the numerically calculated finite-size susceptibilities. We suspect
this error could itself be a finite-size effect that slowly vanishes in the
thermodynamic limit. Figure S3: (a) The finite size scaling plot in the main
text, Fig. 2(c), without the scaling with system size with
$\delta=(\Omega_{c}-\Omega)/\Omega_{c}$. (b) and (c) Finite size scaling for
weak dissipation $\gamma/J=0.2$ and strong dissipation $\gamma/J=10$. For
panel (b), although the dissipation is weak, it is still sufficiently strong
to avoid the condensation physics. Figure S4: (a) Critical behavior for weak
dissipation using the approximate expression Eq. S48 for
$\langle\hat{m}\rangle$ with $\gamma/J=10^{-2}$. We use the analytical
expression for free energy and weak dissipation limit $g(\xi)\sim\xi^{3/2}$.
The dashed line is the infinite size limit. This suggests we should fit to the
power law scaling part to overcome the finite size effect. (b) and (c) Full
numerical results of critical behavior for weak dissipation using the exact
solution. Comparing with the approximation in panel (a), we see that the
approximation captures the essential effects of the condensation physics for
sufficiently small $\gamma/J$. By fitting to the power law part, we obtain the
critical exponent $\beta$. As $\gamma/J$ decreases, the exponent approaches
the analytic prediction $\beta=1$.

##  III Special points of anisotropy /J

###  III.1 No phase transition at special points

In this section, we discuss properties of the exact solution when the
anisotropy (or interaction) parameter $\Delta$ is tuned to a special point
$\Delta_{l,m}$ defined by $\Delta_{l,m}/J=\cos\frac{l}{m}\pi\equiv\cos\eta$
for $m,l\in\mathbb{Z}$. We find that at these special values of $\Delta$,
there is no critical behaviour as a function of $\Omega$ in the vicinity of
$\Omega_{c}$. We can establish this lack of a phase transition rigorously in
the limit of weak dissipation $\gamma\to 0$. Numerically, we find that the
lack of critical behavior persists even for $\gamma\sim J$.

First, using the exact solution constraint equation Eq. S13, we immediately
see that $\Delta=\Delta_{l,m}$ are special because the product of the MPS
coefficients $b_{nm-1}c_{nm-1}=0$ for every positive integer $n\in\mathbb{N}$.
In general, we can always choose an MPS gauge for which $b_{nm-1}=0$; however,
for finite $\gamma>0$ we must have all $c_{k}\neq 0$ in order to have a
boundary state with $\alpha_{k+1}<\infty$ [cf. Eq. (S14)]. For the classical
random walk model, the vanishing $b_{nm-1}=0$ coefficients imply that a
particle on site $nm$ can never hop further to the left, i.e. to sites $k$
with $k\leq nm-1$.

In the weak dissipation limit, a suitable gauge for which $b_{nm-1}=0$ is the
same gauge given in Eqs. (S35) and (S36). Notice that in this gauge,
$c_{nm}\propto\gamma$ and so becomes arbitrarily small for $\gamma\to 0$. In
the zero dissipation limit, we thus have that each site $nm$ in the auxiliary
lattice ($n\in\mathbb{N}$) becomes an absorbing site: there is no hopping out
of this site, either to the right or the left. Recall that our exact solution
MPS corresponds to a particle starting on site $0$ of the auxiliary lattice
and attempting to hop to other sites. For $\gamma\to 0$ and for
$\Delta=\Delta_{l,m}$, we thus have that a particle that starts on site 0 can
never move past site $m$, and if it reaches this site, will become stuck
there. Therefore, at a special point $\eta=\frac{l}{m}\pi$ (and for weak
$\gamma$), the spin chain has an effective chain length $N_{*}=m$, even in the
thermodynamic limit. This implies that there can be no phase transition at the
special values of anistropy $\Delta_{l,m}$ in the zero dissipation limit.

Away from the weak dissipation limit, only the cutoff $b_{nm-1}=0$ remains.
This implies that starting in site $nm$, hopping to the left is still
impossible, but there is an amplitude for hopping to the right. Despite there
no longer being a strict localization of the dynamics in the auxiliary
lattice, we nevertheless have numerical evidence that there still is no
critical behavior near the critical drive $\Omega_{c}$ at the special points
even for $\gamma\sim J$. As we see in Fig. S5(a) and (b), there is not a sharp
transition in the magnetization at $\Omega_{c}$. In Fig. S5(a), we see that
the magnetization smoothly decreases and is highly $\gamma$-dependent for the
special point $\eta=\frac{2}{5}\pi$. In Fig. S5(b), we see that for fixed
$\gamma$ and $\eta=\frac{2}{5}\pi$, the magnetization susceptibility does not
continue to sharpen with increasing system size. This is consistent with the
notion of a effective chain length $N_{*}<N$. We also find in Fig. 4(a) of the
main text that the entanglement entropy shows an area law at the special point
$\eta=\frac{\pi}{3}$ vs. logarithmic scaling $S_{N/2}\sim\log N$ away from
special points.

As a final remark, we also note that at a special point for any $\gamma$, the
matrix elements $|a_{k}|^{2}=|a_{k+m}|^{2}$, $|b_{k}|^{2}=|b_{k+m}|^{2}$, and
$|c_{k}|^{2}=|c_{k+m}|^{2}$ becomes periodic with period $m$ in the auxiliary
space $\mathcal{V}$. Recall that there is always a gauge in which
$b_{nm-1}=0$. Together with periodicity, this reduces the transfer matrix
$\hat{T}$ to an block upper triangular form:

$T=\begin{pmatrix}T_{m}&O_{m}&0&\cdots\\\ 0&T_{m}&O_{m}&\\\ 0&0&T_{m}&\\\
\vdots&&&\ddots&O_{m}\\\ &&&&T_{m}\end{pmatrix}.$ (S49)

Here $T^{\dagger}_{m}$ and
$(O_{m})_{jk}=\delta_{j,m}\delta_{k,1}|c_{m-1}|^{2}$ are $m\times m$ matrices.
This transfer matrix structure leads to finite number of (degenerate)
eigenvalues and exponentially decaying spatial correlations. We leave
determining if and how exactly this structure prevents critical behavior for
future investigation.

Figure S5: Special points of anisotropy $\Delta/J=\cos\frac{l}{m}\pi$ do not
show critical behavior. (a) Magnetization is plotted vs. drive strength
$\Omega/J$ for the special point $\eta=\frac{2}{5}\pi$ for various dissipation
rates $\gamma/J$. No sharp transition in the magnetization is observed, even
for a $N=500$ site chain. (b) For the same $\eta$ as (a), the magnetization
susceptibility is plotted vs drive strength for fixed $\gamma/J=1$ but for
various system sizes. Contrasting with the non-special anisotropy point
$\Delta/J=0.2$ [cf. Fig. S2], here we see no evidence for a sharpening of the
change in susceptibility around $\Omega_{c}\approx 0.95$.

###  III.2 Resonances in magnetization at special points

In this section, we explain briefly how the resonance features, shown in Fig.
3 of the main text, emerge from the classical random walk model. For
$\gamma/J\ll 1$, the transfer matrix in the convenient gauge defined by Eqs.
(S35) and (S36), for the $p$th unit cell
$\hat{T}^{\dagger}_{m;p}=\sum_{k=pm}^{pm+m-1}\cos^{2}k\eta\ket{\mathbf{k}}\bra{\mathbf{k}}+\frac{1}{2}\sin^{2}[(k+1)\eta]\ket{\mathbf{k}}\bra{\mathbf{k+1}}+\frac{1}{2}\sin^{2}k\eta\ket{\mathbf{k+1}}\bra{\mathbf{k}}+\frac{\gamma^{2}}{2\sin^{2}\eta}\ket{\mathbf{1}}\bra{\mathbf{0}}$.
Here $T_{m;p}$ is the $p$th diagonal block $T_{m}$ of Eq. (S49). To this order
in $\gamma$, we have no amplitude for hopping from site $k=(pm-1)\to k=pm$
(though the amplitude for the reverse process is non-zero). To next order in
$\gamma$, rightwards hopping to this bond is possible but with a very weak
amplitude $\sim\gamma^{2}/(\sin^{2}\eta)$. For an auxiliary-space particle
initialized on site 0 to reach a final coordinate $k=nm+s$, the particle has
to pass a weak bond $\gamma^{2}/(\sin^{2}\eta)$ at least $n$ times. The
transfer matrix is therefore
$f_{N}(k)=\bra{\mathbf{k}}(\hat{T}^{\dagger})^{N}\ket{\mathbf{0}}\approx
C(n,N)f(s)(\frac{\gamma^{2}}{\sin^{2}\eta})^{n}$ to leading order in $\gamma$.
Here $C(n,N)$ counts the number configurations that cross $n$ weak bonds to
reach $k$, and $f(s)$ is a function that depends only on the position inside
of the unit cell. This leads to a potential that gets deeper as $\gamma\to 0$.
Therefore, to move away from the $\langle\hat{m}\rangle=0$ configuration, a
larger field $\theta$ is required (thus, a correspondingly _weaker_ drive
$\Omega$, cf. Eq. (S38)). For $\gamma\ll\Omega$, we expect the magnetization
to be close to zero, and we see this numerically: decreasing the dissipation
results in higher resonance peaks around the special anisotropy points which
agrees with our crude argument.

###  III.3 Survival of the phase transition in the thermodynamic limit

Lastly, we argue that although there are infinitely many special points
$\eta=\frac{l}{m}\pi$ in the thermodynamic limit (i.e., all rationals $l/m$),
the phase transition for generic points persists. The crucial observation is
that the resonances get sharper with increasing $N$. We assume that this
sharpening persists so that they become arbitrarily sharp in the thermodynamic
limit. At any fixed $N$, even though any generic non-special $\eta$ will fall
within the resonance width of some special point, only two things can be:
Either that special point will have a transfer matrix cutoff dimension
$m\gtrsim N$ (in the weak dissipation limit), in which case there is no
resonance. Or there is a resonance, in which case we can always increase $N$
until the generic point no longer falls within the resonance width of the
problem special point.

We can be a little bit more formal, again only assuming the sharpening of
resonance peaks with $N$. Similar to the argument for the fractal points of of
the Drude weight [11], we assume the generic point $\Delta/J=\cos\lambda\pi$
where $\lambda$ is a irrational number. By Dirichlet’s approximation theorem,
we can always find a rational $\lambda_{q}=\frac{p}{q}$ that approximates
$\lambda$ with an error bounded by $|\lambda_{q}-\lambda|\lesssim 1/q^{2}$.
Thus, better approximations of $\lambda$ necessarily require larger $q$ for
generic $\lambda$, and so have larger effective chain lengths $N_{*}\sim q$.
One may worry about cases where $\lambda$ is very close to a special point
with a very small $q$ (e.g., $\eta=\frac{\pi}{3}$). We can always find a
rational number $\lambda_{s}$ satisfying $\lambda<\lambda_{s}<\lambda_{q}$
(taking $\lambda<\lambda_{q}$) that better approximates $\lambda$, but with a
greater denominator $s$. So for sufficiently long chains, any non-special
anisotropy point will exhibit critical behavior around $\Omega_{c}$.

##  IV Transport

In this section, we briefly introduce the transport property of the coherent
boundary driven-dissipative spin chain Eq. S2. The spin current is defined by
$\langle
j\rangle\equiv{i\langle(\hat{\sigma}_{A,m}^{+}\hat{\sigma}_{A,m+1}^{-}-{\rm
H.c.})\rangle_{\rm ss}}$. We numerically find that the current scaling with
system size to be the same as in 1D spin chains driven by incoherent pumping
[4, 5, 12]: ballistic $\langle j\rangle\sim N^{0}$ in the easy-axis regime
$\Delta/J<1$, insulating $\langle j\rangle\sim\exp(-\alpha N)$ for
$\Delta/J>1$, and sub-diffusive $j\sim N^{-2}$ at the Heisenberg point
$\Delta=1$. These scaling laws are shown in Fig. S6. The current is ballistic
for $\Delta/J<1$ in both phases $(1)$ and $(2)$ [cf. Fig. 2 of the main text].
It is $\Omega$-dependent, and its non-analytic behavior similar to that of the
average magnetization is numerically observed at the critical drive
$\Omega_{c}$. Finally, we note that although phases $(1)$ and $(3)$ are
asymptotically connected at $\Omega=0$, with each approaching the
ferromagnetic state, they remain distinct in their transport properties.

Figure S6: Transport properties of the coherent boundary driven spin chain in
the three different phases [cf. Fig. 2 of the main text]. The current $\langle
j\rangle$ is plotted as a function of system size $N$. In the easy-axis regime
$\Delta/J<1$ the current is ballistic, and differs only by a prefactor between
the low density phase $(1)$ and the high density phase $(2)$. At the
Heisenberg point $\Delta/J=1$, the current scaling is sub-diffusive, $\langle
j\rangle\sim N^{-2}$, as shown by the green curve closely following the black
dashed line $\sim N^{-2}$. In the inset we plot the current in the
$\Delta/J>1$ regime with $\Delta/J=1.2$. The current falls exponentially with
$N$, indicative of an insulating phase. For all curves we pick $\gamma/J=1$.

##  V Hidden Time Reversal Symmetry

In this section, derive the Onsager symmetry mentioned in the main text. As a
consequence of hidden time reversal symmetry, there are constraints on
observables, including that a particular subset of two-time correlation
functions obey Onsager time symmetry [2]:

${\rm Tr}(\hat{X}(t)\mathcal{J}[\hat{Y}]\hat{\rho}_{\rm ss})={\rm
Tr}(\hat{Y}(t)\mathcal{J}[\hat{X}]\hat{\rho}_{\rm ss}).$ (S50)

Here $\mathcal{J}[\hat{O}]=\hat{\rho}^{1/2}_{\rm ss}\hat{T}_{\rm
h}\hat{O}^{\dagger}\hat{T}_{\rm h}^{-1}\hat{\rho}^{-1/2}_{\rm ss}$ is the time
reversal superoperator with the hidden time reversal antiunitary operator
$\hat{T}_{\rm h}$. This symmetry is nonlinear in the steady state density
matrix, so the time-reversal of generic observables is not necessarily an
easily experimentally measurable quantity. However, there is a class of
operators that are invariant under time reversal,
$\mathcal{J}[\hat{X}]=\hat{X}$, hence obey standard Onsager symmetry

$\langle\hat{X}(t)\hat{Y}(0)\rangle_{\rm
ss}=\langle\hat{Y}(t)\hat{X}(0)\rangle_{\rm ss}.$ (S51)

This subset always include the effective non-Hermitian Hamiltonian
$\hat{H}_{\rm eff}=\hat{H}-\frac{i}{2}\sum_{k}\hat{L}_{k}$ and the jump
operators $\hat{L}_{k}$. One may generate several other time-reversal
symmetric operators by taking linear combinations and products of
$\hat{H}_{\rm eff}$ and $\hat{L}_{k}$ and using the property
$\mathcal{J}[\hat{X}\hat{Y}]=\mathcal{J}[\hat{Y}]\mathcal{J}[\hat{X}]$. This
implies, for example, that the commutator is anti-symmetric,
$\mathcal{J}\\{[\hat{X},\hat{Y}]\\}=-[\hat{X},\hat{Y}]$, and the nested
commutator is symmetric,
$\mathcal{J}\\{[[\hat{X},\hat{Y}],\hat{Y}]]\\}=[[\hat{X},\hat{Y}],\hat{Y}]]$.

In the case of the coherently driven model, the non-Hermitian Hamiltonian is
$\hat{H}_{\rm eff}=\hat{H}_{\rm
XXZ}+\frac{\Omega}{2}\hat{\sigma}^{x}_{N}+\frac{i\gamma}{2}\hat{\sigma}_{1}^{+}\hat{\sigma}_{1}^{-}$,
and the single jump operator is $\hat{\sigma}_{1}^{-}$. The latter can be
readily recovered through measurements of $\hat{\sigma}^{x}_{1}$ and
$\hat{\sigma}^{y}_{1}$. Taking the nested commutator of $\hat{H}_{\rm eff}$
with $\hat{\sigma}_{1}^{-}$, which is time-reversal symmetric, we find

$\mathcal{J}\left[-\frac{1}{2}[[\hat{H}_{\rm
eff},\hat{\sigma}_{1}^{-}],\hat{\sigma}_{1}^{-}]\right]=-\frac{1}{2}[[\hat{H}_{\rm
eff},\hat{\sigma}_{1}^{-}],\hat{\sigma}_{1}^{-}]=\hat{\sigma}_{1}^{-}\hat{H}_{\rm
eff}\hat{\sigma}_{1}^{-}=J\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{-}.$ (S52)

The expectation value of this operator can be recovered from joint $X$ and $Y$
measurements of the two spins. Thus, the correlation function
$\langle\hat{\sigma}_{1}^{-}(t)(\hat{\sigma}_{1}^{-}\sigma_{2}^{-})(0)\rangle$
obeys Onsager symmetry, as noted in the main text:

$\langle\hat{\sigma}_{1}^{-}(t)(\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{-})(0)\rangle_{\rm
c}=\langle(\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{-})(t)\hat{\sigma}_{1}^{-}(0)\rangle_{\rm
c}.$ (S53)

Hidden time reversal symmetry can be frequently broken by either Hamiltonian
or dissipative perturbations. For example, consider when the dissipative bath
is at finite temperature, with thermal occupation $\bar{n}_{\rm th}\neq 0$.
The non-zero thermal occupation modifies Eq. (S2) to
$\partial_{t}\hat{\rho}=\hat{\mathcal{L}}\hat{\rho}=-i[\hat{H}_{\rm
XXZ}+\frac{\Omega}{2}\hat{\sigma}_{N}^{x},\hat{\rho}]+\bar{n}_{\rm
th}\gamma\mathcal{D}[\hat{\sigma}_{1}^{+}]\hat{\rho}.+(1+\bar{n}_{\rm
th})\gamma\mathcal{D}[\hat{\sigma}_{1}^{-}]\hat{\rho}$. In this case, the
Onsager symmetry of the correlation functions in Eq. S53 is broken [See Fig.
S7].

Figure S7: The Onsager symmetry for two-time correlation functions
$\hat{X}=\hat{\sigma}_{1}^{-}$ and $\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{-}$
[cf. Fig. 4(b) in the main text] is broken when the dissipative bath is at
finite temperature, here $\bar{n}_{\rm th}=0.1$. This is shown for $N=3$,
$\Omega/J=0.4$, $\Delta/J=0.2$ and $\gamma/J=1$.

##  VI CQA solution for the incoherent model

As mentioned in the main text, the CQA method is also applicable to other
previously solved models [9, 4, 13]. Here we demonstrate the CQA solution to a
well known exact solution: the incoherent pumping and loss XXZ chain solved in
Ref. [4]. The system is described by the Lindblad master equation

$\partial_{t}\hat{\rho}=\mathcal{L}\hat{\rho}=-i[\hat{H}_{\rm
XXZ},\hat{\rho}]+\gamma\mathcal{D}[\hat{\sigma}_{1}^{-}]\hat{\rho}+\gamma\mathcal{D}[\hat{\sigma}_{N}^{+}]\hat{\rho},$
(S54)

where the Hamiltonian is the usual $\hat{H}_{\rm
XXZ}=\sum_{j=1}^{N-1}(J\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+J\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{y}+\Delta\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z})$
with loss on site 1 and incoherent pumping on site $N$. We take the pumping
and loss rates to be equal. we make the CQA ansatz by constructing a mirrored
absorber, $\hat{H}^{(B)}=-\hat{H}^{(A)}$, and coupling the system and absorber
via two chiral waveguides, yielding the CQA master equation

$\displaystyle\partial_{t}\hat{\rho}_{\rm AB}$ $\displaystyle=\mathcal{L}_{\rm
AB}\hat{\rho}_{\rm AB}=-i[\hat{H}^{(AB)},\hat{\rho}_{\rm
AB}]+\gamma\mathcal{D}[\hat{L}_{\rm AB}^{(1)}]\hat{\rho}_{\rm
AB}+\gamma\mathcal{D}[\hat{L}_{\rm AB}^{(2)}]\hat{\rho}_{\rm AB}.$ (S55)
$\displaystyle\hat{H}^{(AB)}$ $\displaystyle=\hat{H}^{(A)}_{\rm
XXZ}-\hat{H}^{(B)}_{\rm XXZ}+\hat{H}_{\rm c,1}+\hat{H}_{\rm c,N},$ (S56)
$\displaystyle\hat{H}_{\rm c,1}$
$\displaystyle=-\frac{i\gamma}{2}(\hat{\sigma}_{A,1}^{+}\hat{\sigma}_{B,1}^{-}-{\rm
h.c.}),\quad\hat{H}_{\rm
c,N}=-\frac{i\gamma}{2}(\hat{\sigma}_{A,N}^{-}\hat{\sigma}_{B,N}^{+}-{\rm
h.c.}).$ (S57)

Here $\hat{L}_{AB}^{(1)}=\hat{\sigma}_{A,1}^{-}-\hat{\sigma}_{B,1}^{-}$ and
$\hat{L}_{AB}^{(2)}=\hat{\sigma}_{A,N}^{+}-\hat{\sigma}_{B,N}^{+}$ are the
collective dissipators mediated by the chiral waveguides.

A pure steady state $\ket{\psi_{\rm CQA}^{I}}$ must satisfy the dark state
conditions

$\hat{H}_{\rm AB}\ket{\psi_{\rm CQA}^{I}}=0,\quad\hat{L}^{(1)}_{\rm
AB}\ket{\psi_{\rm CQA}^{I}}=0,\quad\hat{L}^{(2)}_{\rm AB}\ket{\psi_{\rm
CQA}^{I}}=0.$ (S58)

Notice that Eq. S54 and Eq. S3 are identical except for the driving term on
site $N$: $\hat{H}_{\rm d}$ in Eq. S54 is replaced by $\hat{H}_{{\rm c},N}$.
We make the same MPS ansatz $|\psi^{I}_{\rm
CQA}\rangle=\sum_{s_{1},\cdots,s_{N}}\mathbf{v_{L}^{\dagger}}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}\mathbf{v_{R}}|s_{1}\cdots
s_{N}\rangle$ as before. We do not assume these to be the same matrices and
vectors as for the coherently driven model, but the Hamiltonian dark state
condition in the bulk and on the loss site $1$ can be satisfied by
constructing the same matrices as above [cf. Eq. (S8)], and by taking
$\mathbf{v_{L}^{\dagger}}=\langle 0|$. Moreover, the dark state condition for
the pumping on site $N$ imposes the constraint $\hat{A}_{0}\mathbf{v_{R}}=0$,
which is satisfied by $\mathbf{v_{R}}=|0\rangle$. Note that the boundary
vectors being identical preserves the $\mathbb{Z}_{2}$ symmetry of the master
equation.

We, therefore, find the CQA solution of the incoherently driven model to be

$|\psi^{I}_{\rm
CQA}\rangle=\sum_{s_{1},\cdots,s_{N}}\bra{0}\hat{A}_{s_{1}}\cdots\hat{A}_{s_{N}}\ket{0}|s_{1}\cdots
s_{N}\rangle,$ (S59)

where $s_{i}\in\\{0_{i},1_{i},T_{i}\\}$, and the MPS matrix coefficients are
given by Eqs. (S12) and (S13). One may make the gauge choice $c_{k}=1$, hence
Eq. (S13) gives the coefficients $b_{k}$. We note that our CQA solution is a
specific purification of the exact solution density matrix found in Ref. [4];
moreover, we may follow the same steps as above and trace out the absorber
system to write $\hat{\rho}_{\rm
ss}^{I}=\hat{\Psi}_{I}\hat{\Psi}_{I}^{\dagger}$ as a Cholesky decomposition.
Finally we point out that $|\psi_{\rm CQA}\rangle$ is a pure _current-
carrying_ steady state of Eq. S55.

## References

  * Stannigel  _et al._ [2012] K. Stannigel, P. Rabl, and P. Zoller, Driven-dissipative preparation of entangled states in cascaded quantum-optical networks, New J. Phys. 14, 063014 (2012).
  * 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).
  * 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, Phys. Rev. A 78, 042307 (2008).
  * Prosen [2011] T. Prosen, Exact nonequilibrium steady state of a strongly driven open XXZ chain, Phys. Rev. Lett. 107, 137201 (2011).
  * Prosen [2015] T. Prosen, Matrix product solutions of boundary driven quantum chains, J. Phys. A Math. Theor. 48, 373001 (2015).
  * Zhang and Barthel [2024] Y. Zhang and T. Barthel, Criteria for davies irreducibility of markovian quantum dynamics, J. Phys. A Math. Theor. 57, 115301 (2024).
  * Frigerio [1977] A. Frigerio, Quantum dynamical semigroups and approach to equilibrium, Letters in Mathematical Physics 2, 79 (1977).
  * Lingenfelter  _et al._ [2024] A. Lingenfelter, M. Yao, A. Pocklington, Y.-X. Wang, A. Irfan, W. Pfaff, and A. A. Clerk, Exact results for a boundary-driven double spin chain and resource-efficient remote entanglement stabilization, Phys. Rev. X. 14, 021028 (2024).
  * Karevski  _et al._ [2013] D. Karevski, V. Popkov, and G. M. Schütz, Exact matrix product solution for the boundary-driven lindblad XXZ chain, Phys. Rev. Lett. 110, 047201 (2013).
  * Sinai [1983] Y. G. Sinai, The limiting behavior of a one-dimensional random walk in a random medium, Theory of Probability & Its Applications 27, 256 (1983).
  * Agrawal  _et al._ [2020] U. Agrawal, S. Gopalakrishnan, R. Vasseur, and B. Ware, Anomalous low-frequency conductivity in easy-plane XXZ spin chains, Phys. Rev. B. 101, 224415 (2020).
  * Bertini  _et al._ [2021] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, Finite-temperature transport in one-dimensional quantum lattice models, Rev. Mod. Phys. 93, 025003 (2021).
  * Prosen [2014] T. Prosen, Exact nonequilibrium steady state of an open hubbard chain, Phys. Rev. Lett. 112, 030603 (2014).