# Liouvillian skin effect in an exactly solvable model

Fan Yang  Department of Physics, Stockholm University, AlbaNova University
Center, 106 91 Stockholm, Sweden  Qing-Dong Jiang  Department of Physics,
Stockholm University, AlbaNova University Center, 106 91 Stockholm, Sweden
Tsung-Dao Lee Institute and School of Physics and Astronomy, Shanghai Jiao
Tong University, 200240, China  Emil J. Bergholtz  Department of Physics,
Stockholm University, AlbaNova University Center, 106 91 Stockholm, Sweden

(August 27, 2024)

###### Abstract

The interplay between dissipation, topology and sensitivity to boundary
conditions has recently attracted tremendous amounts of attention at the level
of effective non-Hermitian descriptions. Here we exactly solve a quantum
mechanical Lindblad master equation describing a dissipative topological Su-
Schrieffer-Heeger (SSH) chain of fermions for both open boundary condition
(OBC) and periodic boundary condition (PBC). We find that the extreme
sensitivity on the boundary conditions associated with the non-Hermitian skin
effect is directly reflected in the rapidities governing the time evolution of
the density matrix giving rise to a Liouvillian skin effect. This leads to
several intriguing phenomena including boundary sensitive damping behavior,
steady state currents in finite periodic systems, and diverging relaxation
times in the limit of large systems. We illuminate how the role of topology in
these systems differs in the effective non-Hermitian Hamiltonian limit and the
full master equation framework.

##  I Introduction

Topological phenomena in the non-Hermitian (NH) realm has attracted ample
interest during the past few years [1, 2, 3]. Compared to their conventional
Hermitan counterparts [4, 5, 6], NH effective Hamiltonians exhibit an entirely
different catalog of gapped and gapless topological phases [1, 2, 3, 7, 8, 9,
10, 11, 12, 13, 14]. The arguably most dramatic effect is caused by a
macroscopic piling up of eigenstates at the boundaries of the system [15, 16,
17], leading to a spectral sensitivity that grows exponentially with system
size when coupling the boundaries [18]. This phenomenology has been dubbed the
NH skin effect [19] whereby the celebrated bulk-boundary correspondence in
Hermitian systems is replaced by a dichotomy described by either non-Bloch
band invariants [19] or in terms of a biorthogonal bulk boundary
correspondence [18]. This has led to a blossoming field of research [19, 18,
15, 3, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 16, 17,
27, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45].

Experiments displaying the NH bulk-boundary correspondence have so far been
limited to classical systems including mechanical [40, 41], electrical [42,
43], and photonic [44, 45] platforms. Moreover, the remarkable sensitivity to
boundary conditions has recently been suggested to be harnessed in
applications such as NH topological sensors [38], and a quantum input-output
theory of such systems has been developed [46, 47]. A fully consistent quantum
mechanical description in terms of Lindbland master equations [48] appropriate
for Markovian dissipative systems [49, 50, 51, 52, 53, 54] has earlier been
studied and fruitfully employed in the context of preparing or stabilizing
Hermitian topological phases [55, 56, 57, 58, 59]. Recent pioneering work has
highlighted that also in such fully quantum mechanically consistent
descriptions, the sensitivity to boundary conditions remains [60, 61, 62, 63,
64, 65, 66, 67], but it is fair to say that a comprehensive understanding is
still lacking.

FIG. 1: Illustration of a dissipative SSH chain under periodic and open
boundary conditions. Lindblad jump operators (solid arrows) are introduced to
a total of $2N$ bonds with an explicit form in Eqs. (2.6), representing the
gain and loss from the environment. To switch from PBC to OBC, we remove one
$B$ site in the last unit cell. As a consequence, the bond dissipators at the
boundary (dashed arrows) only act on single $A$ sites marking two ends of the
open chain, and share the form of Eqs. (2.7).

Here we advance the understanding of this problem by providing a complete
analytical solution of a class of dissipative fermionic chains for both open
and periodic boundary conditions (Fig. 1). Remarkably, we find that there is a
NH SSH Hamiltonian $H_{\text{S}}$ (2.8) that fully diagonalizes the
Liouvillian: All normal modes can be expressed in terms of the eigenstates
$H_{\text{S}}$, and the rapidities, $\beta_{m}$, of the Liouvillian are simply
related to the energy eigenvalues $E_{m}$ of $H_{\text{S}}$ as

$\displaystyle\beta_{m}=\text{const}+iE_{m}.$ (1.1)

This directly implies that the topological properties and the skin effect of
$H_{\text{S}}$ carry over to the quantum context, including fluctuations and
quantum jumps, although the interpretation, as we illustrate, is somewhat
altered and there is no skin effect in the steady state reached at
sufficiently long times.

The effective Hamiltonian, $H_{\text{eff}}$, of our dissipative model ignoring
quantum jumps is also of the form of a NH SSH model in Eq. 2.57 though the
effective parameters reflect different aspects of the underlying dissipation
as illustrated in Fig. 2. Using recent insights into exactly solving the full
set of eigenstates and eigenvalues of NH SSH models [34], this allows us to
carry out a detailed study of dynamical phenomena in the Lindblad setting,
comparing to the much more studied NH phenomenology and revealing several
interesting features, including anomalous damping behavior and diverging
relaxation times (in large systems).

FIG. 2: Illustration of the full dissipative model in (a). (b) shows the
Hamiltonian $H_{\text{S}}$ in Eq. 2.8 whose eigenstates diagonalizes the full
dissipative model (a), and (c) illustrates the effective short-time
Hamiltonian $H_{\text{eff}}$ in Eq. 2.57 describing the short-time dynamics
before any quantum jump has taken place. Intriguingly, both (b) and (c) are NH
SSH models although their asymmetric hopping parameters $\gamma_{i}$ and
$\eta_{i}$ in Eqs. (2.12) reflect different aspects of the microscopic
dissipative processes in (a). $\gamma_{i}$ evaluates the total strength of
loss and gain dissipations on a given bond, while the imbalance between the
two gives rise to $\eta_{i}$. All models exhibit skin effects as discussed in
the text.

The structure of the paper is organized as follows. In Section II, we
introduce the fermionic bond-dissipative SSH model in a Majorana fermion
representation. The exact solutions to the quadratic Lindbladian, which can be
mapped to a generalized NH SSH chain up to a total rapidity shift, are
constructed according to the periodic and open boundary conditions. We point
out the differences and links between the full Liouvillian spectrum and the
energy spectrum of the truncated effective Hamiltonian that neglects quantum
jumps. In Section III, we analyze the configurations of the non-equilibrium
steady state (NESS). When the boundary is closed, the Liouvillian gap vanishes
in an intermediate window of hopping amplitudes and a persistent current flow
is identified there. In Section IV, we discuss the Liouvillian skin effects in
an open quantum system. The anomalous quantum dynamics are manifested in the
relaxation of the current, the modulated damping behavior, and the lifetime
dependence of an edge mode on the system size.

##  II The Model

Our model is built on a bond-dissipative SSH chain of spinless fermions, shown
in Fig. 1. There are hopping terms along the chain with alternating amplitudes
$t_{1}$ and $t_{2}$:
$\mathcal{H}=\sum_{j}t_{1}a_{j,A}^{\dagger}a_{j,B}+t_{2}a_{j+1,A}^{\dagger}a_{j,B}+{\text{H.c}}$.
In the presence of single-particle loss and gain bond dissipations, the set of
Lindblad operators takes the general form:
$L^{l}=\sum_{(j,\alpha)}f^{l}_{(j,\alpha)}a_{j,\alpha}$ and
$L^{g}=\sum_{(j,\alpha)}f^{g}_{(j,\alpha)}a^{\dagger}_{j,\alpha}$.
$a_{j,\alpha}^{\dagger}$ ($a_{j,\alpha}$) creates (annihilates) a fermion in
the $j$th unit cell that belongs to the sublattice $\alpha=A$ or $B$ and they
satisfy the fermionic anti-commutation relations:
$\\{a_{j,\alpha},a_{j^{\prime},\alpha^{\prime}}^{\dagger}\\}=\delta_{j,j^{\prime}}\delta_{\alpha,\alpha^{\prime}}$.
Under different boundary conditions, we choose a total of $n=2N(2N-1)$ sites
for PBC (OBC) with $N$ an integer. Subjected to open boundary, the last unit
cell is broken with an empty $B$ site. The odd number of sites helps to
stabilize a zero-energy boundary mode and eventually leads to an exact
Liouvillian spectrum.

The full dynamics of the generic model is captured by the Lindblad equation
[54, 48],

$\displaystyle\frac{d\rho}{dt}=\hat{\mathcal{L}}\rho\coloneqq-i[\mathcal{H},\rho]+\sum_{\mu}(L_{\mu}\rho
L_{\mu}^{\dagger}-\frac{1}{2}\\{L_{\mu}^{\dagger}L_{\mu},\rho\\}),$ (2.1)

where $\mu$ denotes the summation over all types of Lindblad dissipators.

It is important to note that with linear dissipators, the master Eq. 2.1
generates quadratic Lindbladians, from which the relaxation process of any
observable can be studied by solving the equation of motion either in the
Majorana fermion representation [68, 69, 70, 71, 55, 56, 72, 73, 67] or in the
original complex fermion basis [60, 63]. While the two approaches are
equivalent, we employ the first method due to its advantage of obtaining the
exact Liouvillian spectrum instantaneously.

###  II.1 Liouvillians of non-interacting fermions

Before addressing the specifics of our model, we give a brief review of this
general third quantization approach through which the quadratic Lindbladians
are diagonalized in the canonical basis of Majorana fermions [68, 69, 70].

Starting with $n$ complex fermions, the reduced density matrix of the system
$\rho=\sum_{jk}\rho_{jk}|a_{j}\rangle\langle a_{k}|$ lives in a Hilbert space
of dimension $2^{n}\times 2^{n}$. One can always construct a pair of Majorana
fermions out of one complex fermion, for instance,
$\omega_{j}=a_{l}^{\dagger}+a_{l}$ and $\omega_{k}=-i(a_{l}^{\dagger}-a_{l})$
such that they become their own anti-particles
$\omega_{j}^{\dagger}=\omega_{j}$ and satisfy anti-commutation relations
$\\{\omega_{j},\omega_{k}\\}=2\delta_{j,k}$. The Hamiltonian and the
dissipators take the matrix form depending on the details of the mapping:
$\mathcal{H}=\sum_{j,k}w_{j}H_{j,k}w_{k}$, $L^{\nu}=\sum_{j}l_{j}^{\nu}w_{j}$.

The Hilbert space is now represented in the $2^{2n}$-dimensional Liouville
space $\mathcal{K}$ expanded by the new set of Majorana operators:
$P_{\underline{\alpha}}=w_{1}^{\alpha_{1}}w_{2}^{\alpha_{2}}\cdots
w_{2n}^{\alpha_{2n}}$ with $\alpha_{j}\in\\{0,1\\}$. We apply the notation
$\underline{x}=(x_{1},x_{2},\dots,x_{k})^{T}$ to represent a vector of scalars
or operators. Over the space $\mathcal{K}$, it is convenient to define the
adjoint creation and annihilation linear maps through $\varphi$ fermions:
$\varphi_{j}|P_{\underline{\alpha}}\rangle=\delta_{\alpha_{j},1}|w_{j}P_{\underline{\alpha}}\rangle$,
$\varphi_{j}^{\dagger}|P_{\underline{\alpha}}\rangle=\delta_{\alpha_{j},0}|w_{j}P_{\underline{\alpha}}\rangle$,
which obey $\\{\varphi_{j},\varphi_{k}^{\dagger}\\}=\delta_{j,k}$. One
verifies that
$|w_{j}P_{\underline{\alpha}}\rangle=(\varphi_{j}^{\dagger}+\varphi_{j})|P_{\underline{\alpha}}\rangle$,
$|P_{\underline{\alpha}}w_{j}\rangle=\mathcal{P}_{F}(\varphi_{j}^{\dagger}-\varphi_{j})|P_{\underline{\alpha}}\rangle$.
The operator $\mathcal{P}_{F}=(-1)^{\hat{\mathcal{N}}}$ denotes the Fermi
parity associated with the fermion number
$\hat{\mathcal{N}}=\sum_{j=1}^{2n}\varphi_{j}^{\dagger}\varphi_{j}=\underline{\varphi}^{\dagger}\cdot\underline{\varphi}$.
Since $[\hat{\mathcal{L}},\mathcal{P}_{F}]=0$ and $(\mathcal{P}_{F})^{2}=1$,
the parity is conserved and takes the value $\pm 1$.

In the canonical basis $|P_{\underline{\alpha}}\rangle$, the Liouvillian in
Eq. 2.1 contains the even ($+$) and odd ($-$) parity sectors:

$\displaystyle\hat{\mathcal{L}}=$
$\displaystyle-4i\underline{\varphi}^{\dagger}\cdot H\underline{\varphi}$
$\displaystyle+\frac{(1+\mathcal{P}_{F})}{2}\left[-\underline{\varphi}^{\dagger}\cdot(M+M^{T})\underline{\varphi}+\underline{\varphi}^{\dagger}\cdot(M-M^{T})\underline{\varphi}^{\dagger}\right]$
$\displaystyle+\frac{(1-\mathcal{P}_{F})}{2}\left[-\underline{\varphi}\cdot(M+M^{T})\underline{\varphi}^{\dagger}+\underline{\varphi}\cdot(M-M^{T})\underline{\varphi}\right].$
(2.2)

Here, the matrix $M=\sum_{\nu}M^{\nu}$ arises from the loss and gain
dissipators: $M^{\nu}=(l^{\nu})^{T}(l^{\nu})^{*}$ and $\nu=g,l$.

Since a physical observable must contain an even number of fermionic
operators, one may focus on only a definite parity sector, e.g., the even
parity sector ($\mathcal{P}_{F}=1$) in this paper. Correspondingly, the
dynamics of the system is governed by

$\displaystyle\hat{\mathcal{L}}_{+}=\frac{1}{2}\begin{pmatrix}\underline{\varphi}^{\dagger}\cdot&\underline{\varphi}\cdot\end{pmatrix}\begin{pmatrix}-X^{\dagger}&iY\\\
0&X\end{pmatrix}\begin{pmatrix}\underline{\varphi}\\\
\underline{\varphi}^{\dagger}\end{pmatrix}-A_{0},$ (2.3)

with $X=-4iH+M+M^{T}$, $Y=-2i(M-M^{T})$, and $A_{0}=\frac{1}{2}\text{Tr}[X]$.
After proper diagonalization, one is able to express the Liouvillian in terms
of rapidities $\beta_{m}$ and normal master modes (NMMs) $b^{\prime}_{m}$,
$b_{m}$:

$\displaystyle\hat{\mathcal{L}}_{+}=-\sum_{m=1}^{2n}\beta_{m}b^{\prime}_{m}b_{m},$
(2.4)

with the band index $m$ and NMMs satisfying the anti-commutation relations
$\\{b^{\prime}_{m},b_{l}\\}=\delta_{m,l}$.

Notably, the complex rapidity spectrum contains rich physics. An initial state
$\rho_{0}$ of positive fermion parity approaches the NESS after a long-time
evolution:
$\rho_{\text{ss}}=\left.e^{\hat{\mathcal{L}}_{+}t}\rho_{0}\right|_{t\to\infty}$.
While the imaginary part of the rapidity $\beta_{m}$ encodes the phase
oscillation frequency, the real part of $\beta_{m}$ reveals the relaxation
speed of the system to the steady state and $\text{Re}({\beta_{m}})\geq 0$ is
required naturally. To describe the asymptotic decay rate quantitively, a
spectral Liouvillian gap [74, 66] can be defined as

$\displaystyle\Delta=2\min\\{\text{Re}[\beta_{m}]\\}.$ (2.5)

Extra insight comes from the upper triangular structure of
$\hat{\mathcal{L}}_{+}$. It infers that the rapidity spectrum must coincide
with the eigenvalues of the matrix $X$ in the diagonal block, which is NH and
also called the damping matrix. The off-diagonal block $Y$, on the other hand,
shapes the configurations of NMMs and NESS [71, 72].

###  II.2 Exactly solvable models

Next, we apply the third quantization approach to the bond-dissipative SSH
chain in Fig. 1 and derive exact PBC and OBC Lindblad spectra. Compared with
previous studies, our model is more general and include the special cases of
Refs. [73, 60]. We start from a generic set of linear bond dissipators acting
on both $t_{1}$ and $t_{2}$ bonds,

$\displaystyle\begin{cases}L^{l}_{1,j}&=\sqrt{\gamma^{l}_{1}}(a_{j,A}-ia_{j,B})\\\
L^{g}_{1,j}&=\sqrt{\gamma^{g}_{1}}(a_{j,A}^{\dagger}+ia_{j,B}^{\dagger}),\end{cases}$
$\displaystyle\begin{cases}L^{l}_{2,j}&=\sqrt{\gamma^{l}_{2}}(a_{j,B}-ia_{j+1,A})\\\
L^{g}_{2,j}&=\sqrt{\gamma^{g}_{2}}(a_{j,B}^{\dagger}+ia_{j+1,A}^{\dagger}).\end{cases}$
(2.6)

Under PBC, the index $j$ runs over $N$ unit cells. When the boundary opens up
with the last $B$ site taken away, the associated dissipators are curtailed
simultaneously:

$\displaystyle\begin{cases}L^{l}_{1,N}&=\sqrt{\gamma^{l}_{1}}a_{N,A}\\\
L^{g}_{1,N}&=\sqrt{\gamma^{g}_{1}}a_{N,A}^{\dagger},\end{cases}\quad\begin{cases}L^{l}_{2,N}&=(-i)\sqrt{\gamma^{l}_{2}}a_{1,A}\\\
L^{g}_{2,N}&=i\sqrt{\gamma^{g}_{2}}a_{1,A}^{\dagger}.\end{cases}$ (2.7)

It turns out that the damping matrix $X$ in the Majorana representation can be
transformed to the NH SSH Hamiltonian in Eq. 2.15:

$\displaystyle\mathcal{H}_{\text{S}}$
$\displaystyle=\sum_{j}(t_{1}+\gamma_{1})\varphi_{j,A}^{\dagger}\varphi_{j,B}+(t_{1}-\gamma_{1})\varphi_{j,B}^{\dagger}\varphi_{j,A}$
$\displaystyle\phantom{=}+(t_{2}+\gamma_{2})\varphi_{j,B}^{\dagger}\varphi_{j+1,A}+(t_{2}-\gamma_{2})\varphi_{j+1,A}^{\dagger}\varphi_{j,B}.$
(2.8)

The strengths of asymmetric hopping terms $\gamma_{1}$, $\gamma_{2}$ take the
value $2\gamma_{i}=|\gamma_{i}^{l}|+|\gamma^{g}_{i}|$. In the Appendix, we
construct the exact spectrum of $\mathcal{H}_{\text{S}}$ as a direct
generalization of Ref. [18, 34] to the new limit $\gamma_{2}\neq 0$. It
enables us to build an exact solution to $\hat{\mathcal{L}}_{+}$ under both
PBC and OBC (with a total number of sites $n=2N$ and $n=2N-1$, respectively),
using the eigenvectors of $\mathcal{H}_{\text{S}}$.

####  II.2.1 Majorana representation

As a first step, let us define the Liouville space $\mathcal{K}$ by a mapping
from $n$ spinless fermions to $2n$ Majorana particles:

$\displaystyle\begin{cases}a_{j,A}&=\frac{1}{2}(c_{j,A}-id_{j,A}),\\\
a_{j,B}&=\frac{1}{2}(d_{j,B}+ic_{j,B}).\end{cases}$ (2.9)

For later convenience, we regroup $c$ and $d$ Majorana fermions into a whole
set $\\{w\\}$ under the vector notation:
$\underline{w}=(w_{1},w_{2},\dots,w_{2n})^{T}=(c_{1},\dots,c_{n},d_{1},\dots,d_{n})^{T}$.
Accordingly, in the Lindblad Eq. 2.1, the operators
$\mathcal{H}=\sum_{j,k}w_{j}H_{j,k}w_{k}$ and
$M_{ij}=\sum_{\nu=g,l}(l^{\nu}_{i,\mu})^{T}(l^{\nu}_{\mu,j})^{*}$ with
$L_{\mu}^{\nu}=\sum_{j}l_{\mu,j}^{\nu}w_{j}$ take the following matrix forms:

$\displaystyle\underline{w}^{T}H\underline{w}$
$\displaystyle=\begin{pmatrix}\underline{c}^{T}&\underline{d}^{T}\end{pmatrix}\begin{pmatrix}H_{0}&0\\\
0&H_{0}\end{pmatrix}\begin{pmatrix}\underline{c}\\\
\underline{d}\end{pmatrix},$ $\displaystyle\underline{w}^{T}M\underline{w}$
$\displaystyle=\begin{pmatrix}\underline{c}^{T}&\underline{d}^{T}\end{pmatrix}\begin{pmatrix}M_{1}&iM_{2}\\\
-iM_{2}&M_{1}\end{pmatrix}\begin{pmatrix}\underline{c}\\\
\underline{d}\end{pmatrix}.$ (2.10)

$H_{0}$ and $M_{1,2}$ are $n\times n$ matrices holding entries, for instance,
under OBC:

$\displaystyle H_{0}$
$\displaystyle=\frac{i}{4}\begin{pmatrix}0&t_{1}&&&&&&&\\\
-t_{1}&0&-t_{2}&&&&&&\\\ &t_{2}&0&t_{1}&&&&&\\\ &&-t_{1}&0&&&&&\\\
&&&&&\ddots&&&\\\ &&&&&&0&t_{1}&\\\ &&&&&&-t_{1}&0&-t_{2}\\\
&&&&&&&t_{2}&0\end{pmatrix},$ $\displaystyle M_{1}$
$\displaystyle=\frac{\gamma}{2}\cdot\mathbb{1}_{n\times n}+\frac{1}{2}$
$\displaystyle\times\begin{pmatrix}0&\gamma_{1}&&&&&&&\\\
\gamma_{1}&0&-\gamma_{2}&&&&&&\\\ &-\gamma_{2}&0&\gamma_{1}&&&&&\\\
&&\gamma_{1}&0&&&&&\\\ &&&&&\ddots&&&\\\ &&&&&&0&\gamma_{1}&\\\
&&&&&&\gamma_{1}&0&-\gamma_{2}\\\ &&&&&&&-\gamma_{2}&0\end{pmatrix},$
$\displaystyle M_{2}$ $\displaystyle=M_{1}(\gamma_{i}\rightarrow\eta_{i}).$
(2.11)

$\gamma_{i}$’s and $\eta_{i}$’s stand for the sum and the imbalance of loss
and gain dissipations:

$\displaystyle\gamma=\gamma_{1}+\gamma_{2},\quad
2\gamma_{i}=|\gamma_{i}^{l}|+|\gamma^{g}_{i}|,$
$\displaystyle\eta=\eta_{1}+\eta_{2},\quad
2\eta_{i}=|\gamma_{i}^{l}|-|\gamma^{g}_{i}|.$ (2.12)

It can be seen immediately that the Majorana representation in Eq. 2.9 is
better adapted to diagonalize the Liouvillian in Eq. 2.3: in the adjoint
fermion basis
$\underline{\varphi}^{T}=(\underline{\varphi}_{c}^{T},\underline{\varphi}_{d}^{T})$,
it incorporates the matrix blocks

$\displaystyle X=\begin{pmatrix}X_{c}&0\\\ 0&X_{d}\end{pmatrix},\quad
Y=4\begin{pmatrix}0&M_{2}\\\ -M_{2}&0\end{pmatrix},$ (2.13)

with $X_{c}=X_{d}=-4iH_{0}+2M_{1}$. Notice that the damping matrix $X$ is
diagonal and depends only on the total strength of dissipations $\gamma_{i}$.
By contrast, the matrix $Y$ is off-diagonal (thus couples $\varphi_{c}$ and
$\varphi_{d}$-fermions) and depends only on the imbalance of gain and loss
$\eta_{i}$. Taking into account the identical structure shared by $X_{c}$ and
$X_{d}$, one concludes that the rapidity spectrum determined by the full
damping matrix is at least doubly degenerate,

$\displaystyle\hat{\mathcal{L}}_{+}=-\sum_{m=1}^{n}\beta_{m}(b^{\prime}_{c,m}b_{c,m}+b^{\prime}_{d,m}b_{d,m}),$
(2.14)

where $\beta_{m}=\beta_{c,m}=\beta_{d,m}$.

It is not difficult to find a unitary transformation $U_{n\times
n}=\text{diag}\\{1,i,1,i,\dots,1,i,1\\}$ under which $X_{c(d)}$ is mapped to
the generalized NH SSH Hamiltonian in Eq. 2.8:

$\displaystyle X_{c}=X_{d}=\gamma\cdot\mathbb{1}+iUH_{\text{S}}U^{-1}.$ (2.15)

The matrix form $H_{\text{S}}$ is given in Eq. A8. We thus reveal one
important relation aforementioned in Eq. 1.1 for the rapidity spectrum of our
model:

$\displaystyle\beta_{m}=\gamma+iE_{m},$ (2.16)

where $E_{m}$ represents the eigenvalues of $H_{\text{S}}$.

It can be checked directly that the equality in Eq. 2.16 holds true for the
PBC spectrum as well. Going to the momentum space, we set the lattice spacing
to unity and adopt the Fourier transform
$\varphi_{c,{(j,\alpha)}}=\frac{1}{\sqrt{N}}\sum_{q}e^{iqj}\varphi_{c,\alpha}(q)$
with $q=2\pi m^{\prime}/N$, $m^{\prime}=-N/2,-N/2+1,\dots,0,\dots,N/2-1$, such
that the anti-commutation relations are satisfied:
$\\{\varphi_{c,\alpha}(q),\varphi_{c,\alpha^{\prime}}^{\dagger}(q^{\prime})\\}=\delta_{q,q^{\prime}}\delta_{\alpha,\alpha^{\prime}}$.
The same transform is applied to $\varphi_{d}$-fermions. In the basis of
$\underline{\varphi}(q)=(\varphi_{c,A}(q),\varphi_{c,B}(q),\varphi_{d,A}(q),\varphi_{d,B}(q))^{T}$,
the Liouvillian in Eq. 2.3 turns into

$\displaystyle\hat{\mathcal{L}}^{\text{PBC}}_{+}=$
$\displaystyle\frac{1}{2}\sum_{q}\begin{pmatrix}\underline{{\varphi}}^{\dagger}(q)\cdot&\underline{{\varphi}}(-q)\cdot\end{pmatrix}$
$\displaystyle\times\begin{pmatrix}-X^{\dagger}(q)&iY(q)\\\
0&X(q)\end{pmatrix}\begin{pmatrix}\underline{{\varphi}}(q)\\\
\underline{{\varphi}}^{\dagger}(-q)\end{pmatrix}$ $\displaystyle-A_{0}.$
(2.17)

$X(q)$ and $Y(q)$ become $4\times 4$ matrices inheriting the same structures
as before in Eq. 2.13 with ingredients expressed in terms of Pauli matrices:

$\displaystyle H_{0}(q)$ $\displaystyle=(-it_{2}\sin
q)\cdot\sigma^{x}+i(t_{1}+t_{2}\cos q)\cdot\sigma^{y},$ $\displaystyle
M_{1}(q)$
$\displaystyle=\frac{1}{2}[\gamma\cdot\mathbb{1}+(\gamma_{1}-\gamma_{2}\cos
q)\cdot\sigma^{x}-\gamma_{2}\sin q\cdot\sigma^{y}],$ $\displaystyle M_{2}(q)$
$\displaystyle=\left.M_{1}(q)\right|_{\gamma_{i}\to\eta_{i}}.$ (2.18)

The mapping in Eq. 2.15 from the damping matrix to the NH SSH Hamiltonian
follows naturally under the unitary transformation $U=\text{diag}\\{1,i\\}$
with $H_{\text{S}}(q)$ shown in Eqs. (A2). Therefore, we confirm the validity
of the equality relation in Eq. 2.16 for the PBC rapidity spectrum.

####  II.2.2 Changing boundaries from PBC to OBC

We proceed to construct the complete set of eigenvectors of the damping matrix
based on the mapping in Eq. 2.15. Let us write the generic eigenvalue
equations of the NH SSH Hamiltonian:

$\displaystyle H_{\text{S}}\
\tilde{\underline{\psi}}_{Rm}=E_{m}\tilde{\underline{\psi}}_{Rm},\quad
H_{\text{S}}^{\dagger}\
\tilde{\underline{\psi}}_{Lm}=E_{m}^{*}\tilde{\underline{\psi}}_{Lm},$ (2.19)

of which the exact solutions under different boundary conditions are derived
in the Appendix. It renders that the pair of eigenvectors of $X_{c(d)}$ can be
constructed as

$\displaystyle\underline{\psi}_{Rm}=U\
\tilde{\underline{\psi}}_{Rm},\quad\underline{\psi}_{Lm}=U\
\tilde{\underline{\psi}}_{Lm},$ (2.20)

with corresponding eigenvalues $\beta_{m}$ in consistency with the relation
(2.16):

$\displaystyle X_{c(d)}\
{\underline{\psi}}_{Rm}=\beta_{m}{\underline{\psi}}_{Rm},\quad
X_{c(d)}^{\dagger}\
{\underline{\psi}}_{Lm}=\beta_{m}^{*}{\underline{\psi}}_{Lm}.$ (2.21)

Under this construction, the biorthogonal normalization of the left and right
eigenstates [75, 18, 34] are respected:

$\displaystyle{\underline{\psi}}_{Lm}^{*}\cdot{\underline{\psi}}_{Rl}=\tilde{\underline{\psi}}_{Lm}^{*}\cdot\tilde{\underline{\psi}}_{Rl}=\delta_{m,l}.$
(2.22)

When the boundary is switched from PBC to OBC as depicted in Fig. 1, we can
extract the rapidity spectrum directly from $E_{m}$ solved in Eqs. A3, A9 and
A14:

$\displaystyle\beta_{\pm}^{\text{PBC}}(q)$ $\displaystyle=\gamma\pm
i[t_{1}^{2}+t_{2}^{2}-(\gamma_{1}^{2}+\gamma_{2}^{2})$
$\displaystyle\phantom{=}+2(t_{1}t_{2}+\gamma_{1}\gamma_{2})\cos
q+2i(t_{1}\gamma_{2}+t_{2}\gamma_{1})\sin q]^{1/2},$
$\displaystyle\beta_{\pm}^{\text{OBC}}(q)$ $\displaystyle=\gamma\pm
i[t_{1}^{2}+t_{2}^{2}-(\gamma_{1}^{2}+\gamma_{2}^{2})$
$\displaystyle\phantom{=}+2\sqrt{(t_{1}^{2}-\gamma_{1}^{2})(t_{2}^{2}-\gamma_{2}^{2})}\cos
q]^{1/2},$ $\displaystyle\beta_{m=0}^{\text{OBC}}$ $\displaystyle=\gamma.$
(2.23)

In the OBC spectrum, the band index $m\in\\{(\pm,q),0\\}$ is assigned to
$n=2N-1$ bands with discrete modes $q=\pi m^{\prime}/N$,
$m^{\prime}=1,2,\dots,N-1$. Given an odd number of sites, there emerge two
right and left boundary modes at zero energy $E_{0}=0$ with exponential
localization factors: $r_{R}=-(t_{1}-\gamma_{1})/(t_{2}+\gamma_{2})$,
$r_{L}=-(t_{1}+\gamma_{1})/(t_{2}-\gamma_{2})$ [see also Eq. A9]. As for the
bulk spectrum, analogous to the simplified NH SSH model with only one
asymmetric hopping term $\gamma_{1}$ [19, 20, 34], one finds up to a shift in
$q$ a general relation of rapidities between two boundary conditions:

$\displaystyle\beta_{\pm}^{\text{OBC}}(q)=\beta_{\pm}^{\text{PBC}}(q-i\ln(r)),$
(2.24)

where
$r=\sqrt{(t_{1}-\gamma_{1})(t_{2}-\gamma_{2})/[(t_{1}+\gamma_{1})(t_{2}+\gamma_{2})]}$.

It is important to note that under OBC, in the region $|r_{L}^{*}r_{R}|<1$ or,
equivalently,

$\displaystyle\left|t_{1}^{2}-\gamma_{1}^{2}\right|<\left|t_{2}^{2}-\gamma_{2}^{2}\right|,$
(2.25)

the complete set of eigenvectors exhibits a non-zero biorthogonal polarization
[18] and holds a non-trivial non-Bloch topological invariant [19]. It further
indicates at $|r_{L}^{*}r_{R}|=1$, the gap of $E_{m}$ closes [an alternative
argument is given above Eq. A13]. As a result, when $\mathcal{H}_{\text{S}}$
in Eq. 2.8 reduces to the Hatano-Nelson model at $t_{1}=t_{2}$,
$\gamma_{1}=\gamma_{2}$ [76, 77, 78, 2], a collapse of the exact bulk states
is expected. One is nevertheless able to study the behavior of the system
around these gap closing points via approximate variational states [34].
Meanwhile, in a dissipative quantum system composed of tight-binding bosons
and non-linear (quadratic) Lindblad operators, a Liouvillian can be
constructed in such a way that its diagonal subspace resembles the Hatano-
Nelson model, thus giving rise to a similar Liouvillian skin effect [66].

####  II.2.3 Exactly solvable NMMs

We are now prepared to get an analytical set of NMMs for the Liouvillian in
Eq. 2.3, in particular, under the open boundary condition. The essence is to
remove in the upper triangular structure the off-diagonal block $Y$ that
entangles adjoint fermions $\varphi_{c}$ and $\varphi_{d}$. More precisely,

$\displaystyle\begin{pmatrix}-X^{\dagger}&iY\\\
0&X\end{pmatrix}=W\begin{pmatrix}-X^{\dagger}&0\\\ 0&X\end{pmatrix}W^{-1}.$
(2.26)

We find a solution for the transformation above

$\displaystyle W=\begin{pmatrix}\mathbb{1}_{2n\times 2n}&C\\\
0&\mathbb{1}_{2n\times 2n}\end{pmatrix},$ (2.27)

where the covariance matrix $C$ satisfies

$\displaystyle X^{\dagger}C+CX=iY.$ (2.28)

It is easy to discern that if ${\eta}/\gamma={\eta_{1}}/{\gamma_{1}}$ or
${\eta_{2}}/{\gamma_{2}}$, the covariance matrix holds a simple structure

$\displaystyle C=\frac{i\eta}{\gamma}\begin{pmatrix}0&\mathbb{1}_{n\times
n}\\\ -\mathbb{1}_{n\times n}&0\end{pmatrix}.$ (2.29)

By definition in Eqs. (2.12), the solvable limit encompasses the following
possibilities:

$\displaystyle\gamma_{1}\gamma_{2}=0,\quad\frac{\gamma_{1}}{\eta_{1}}=\frac{\gamma_{2}}{\eta_{2}}\neq
0\quad\text{and}\quad\eta_{1}=\eta_{2}=0.$ (2.30)

A deeper understanding of the covariance matrix comes from the pairing
function of Majorana fermions:
${C}_{jk}(t)=-\text{Tr}[w_{j}w_{k}\rho(t)]+\delta_{jk}$. In its time
evolution, by applying the Lindblad Eq. 2.1 in combination with the anti-
commutation relations of Majorana fermions, one arrives at

$\displaystyle\partial_{t}{C}(t)=-{C}(t)X-X^{\dagger}{C}(t)+iY.$ (2.31)

For a steady state, $\partial_{t}{C}_{s}=0$. Therefore,

$\displaystyle C_{s}\equiv C.$ (2.32)

It can be interpreted that in the solvable limit of the Lindbladian in Eqs.
(2.30), the covariance matrix encodes a stationary pairing pattern favoured by
Majorana fermions at large times:

$\displaystyle\langle
c_{j,\alpha}d_{j^{\prime},\alpha^{\prime}}\rangle_{\text{ss}}$
$\displaystyle=-\frac{i\eta}{\gamma}\delta_{j,j^{\prime}}\delta_{\alpha,\alpha^{\prime}},$
$\displaystyle\langle
c_{j,\alpha}c_{j^{\prime},\alpha^{\prime}}\rangle_{\text{ss}}$
$\displaystyle=\langle
d_{j,\alpha}d_{j^{\prime},\alpha^{\prime}}\rangle_{\text{ss}}=\delta_{j,j^{\prime}}\delta_{\alpha,\alpha^{\prime}}.$
(2.33)

A non-vanishing covariance matrix plays a role in counteracting the effect of
the imbalanced gain and loss dissipations ($\eta$). Back to the original
Hilbert space expanded by spinless fermions $a$, a more physical picture can
be drawn as follows. Since
$n_{j,\alpha}=a_{j,\alpha}^{\dagger}a_{j,\alpha}=(1-ic_{j,\alpha}d_{j,\alpha})/2$,
the dissipative fermionic SSH chain has a tendency to evolve towards a
uniformly distributed configuration with an occupation number:

$\displaystyle n_{\text{ss},(j,\alpha)}=\frac{\gamma-\eta}{2\gamma}\in[0,1].$
(2.34)

There is no correlation between $a$ fermions. In general, one can express the
steady state in the form of a product state,

$\displaystyle\left|\Psi_{\text{ss}}\right>=\bigotimes_{j=1}^{n_{\text{tot}}}\left(\sqrt{\frac{\gamma+\eta}{2\gamma}}\left|0\right>_{j}+e^{i\theta_{j}}\sqrt{\frac{\gamma-\eta}{2\gamma}}\left|1\right>_{j}\right),$
(2.35)

where $\theta_{j}$ denotes an arbitrary phase difference.

In the end, with a covariance matrix fulfilling the transformation in Eq.
2.26, we can relate the NMMs in Eq. 2.14 to the left and right eigenvectors of
the damping matrix in Eqs. (2.21):

$\displaystyle\begin{cases}\
b^{\prime}_{c,m}=\underline{\psi}_{Lm}^{*}\cdot\underline{\varphi}_{c}^{\dagger}\\\
\
b_{c,m}=\underline{\psi}_{Rm}\cdot(\underline{\varphi}_{c}-\frac{i\eta}{\gamma}\underline{\varphi}_{d}^{\dagger}),\end{cases}$
$\displaystyle\begin{cases}\
b^{\prime}_{d,m}=\underline{\psi}_{Lm}^{*}\cdot\underline{\varphi}_{d}^{\dagger}\\\
\
b_{d,m}=\underline{\psi}_{Rm}\cdot(\underline{\varphi}_{d}+\frac{i\eta}{\gamma}\underline{\varphi}_{c}^{\dagger}).\end{cases}$
(2.36)

The anticommutation relations are guaranteed by the biorthonormality in Eqs.
(2.22):
$\\{b^{\prime}_{c,m},b_{c,l}\\}=\\{b^{\prime}_{d,m},b_{d,l}\\}=\delta_{m,l}$
and all others zero.

In the same manner, under PBC, the covariance matrix takes a simple analytical
form in the solvable limit in Eqs. (2.30),

$\displaystyle C(q)=\frac{i\eta}{\gamma}\begin{pmatrix}0&\mathbb{1}_{2\times
2}\\\ -\mathbb{1}_{2\times 2}&0\end{pmatrix},$ (2.37)

and the set of NMMs associated with rapidities
$\beta_{\nu=\pm}^{\text{PBC}}(q)$ is given by

$\displaystyle\begin{cases}\
b^{\prime}_{c,\nu}(q)=\underline{\psi}^{*}_{L,\nu}(q)\cdot\underline{\varphi}_{c}^{\dagger}(-q)\\\
\
b_{c,\nu}(q)=\underline{\psi}_{R,\nu}(q)\cdot\left[\underline{\varphi}_{c}(-q)-\frac{i\eta}{\gamma}\underline{\varphi}_{d}^{\dagger}(q)\right],\end{cases}$
$\displaystyle\begin{cases}\
b^{\prime}_{d,\nu}(q)=\underline{\psi}^{*}_{L,\nu}(q)\cdot\underline{\varphi}_{d}^{\dagger}(-q)\\\
\
b_{d,\nu}(q)=\underline{\psi}_{R,\nu}(q)\cdot\left[\underline{\varphi}_{d}(-q)+\frac{i\eta}{\gamma}\underline{\varphi}_{c}^{\dagger}(q)\right].\end{cases}$
(2.38)

###  II.3 Spectrum and topology

In this section, we reveal the topology of the Liouvillian starting from the
spectral winding of the NH damping matrix. Known for a NH system, the
conventional bulk-boundary correspondence is broken. Yet the prevalence of the
exceptional points (EPs) in the OBC rapidity spectrum implies the Liouvillian
skin effect, of which more physical consequences will be discussed in Section
IV. We further analyze the formation of the Liouvillian gap that is found to
be highly sensitive to the boundary conditions.

####  II.3.1 Spectral winding number and exceptional topology

Based on different classification schemes [73, 79], open fermion matter falls
into one of the ten NH Bernard-LeClair symmetry classes. Given a quadratic
Liouvillian [73], the classification can be defined through the damping matrix
$X$, or, in our case, $X_{c(d)}$:

$\displaystyle Z$ $\displaystyle=-iX_{c(d)},$
$\displaystyle=-i(\gamma_{1}+\gamma_{2})\cdot\mathbb{1}-[i\gamma_{1}+t_{2}\sin
q-i\gamma_{2}\cos q]\cdot\sigma^{x}$
$\displaystyle\phantom{=}+[(t_{1}+t_{2}\cos q)+i\gamma_{2}\sin
q]\cdot\sigma^{y}.$ (2.39)

The matrix $Z$ resembling the Hamiltonian in the closed limit preserves the
time-reversal, particle-hole, and pseudo-anti-Hermiticity (PAH, or generalized
chiral) symmetries: $Z=\sigma^{x}Z^{T}\sigma^{x}$, $Z=-Z^{*}$, and
$Z=-\sigma^{x}Z^{\dagger}\sigma^{x}$. Hence, each subspace of the damping
matrix belongs to class BDI with a $\mathbb{Z}$ classification in 1D. It
should be noted that at the edge, two real Majorana fermions $c$ and $d$
recombine into one complex fermion $a$ as indicated by Eq. 2.9. When the
dissipations are turned on, this edge mode shares a finite lifetime with a
contribution coming from the effective Liouvillian gap (non-vanishing as
indicated by the purely imaginary total energy shift in $Z$) and another from
the non-Hermiticity of the damping matrix [see Eq. 4.9].

In the presence of the PAH symmetry, a $\mathbb{Z}$ classification is captured
by the topological invariant, spectral winding number. By shifting the
reference point to $(0,-i\gamma)$ in the complex plane, it is equivalent to
evaluate the winding of the Bloch Hamiltonian $H_{\text{S}}(q)$ [2, 3]:

$\displaystyle\nu=\frac{1}{2\pi i}\int_{-\pi}^{\pi}dq\
\partial_{q}\ln\\{\det[H_{\text{S}}(q)]\\}.$ (2.40)

Remarkably, when the spectral winding number becomes non-trivial, the left and
right boundary modes in the OBC spectrum in Eq. A9 start to localize at
different ends of the chain:

$\displaystyle|\nu|=1,\quad\text{sgn}[\ln(|r_{L}|)]\neq\text{sgn}[\ln(|r_{R}|)].$
(2.41)

For $\gamma_{1}\gamma_{2}\geq 0$, the above topological regime resides in

$\displaystyle
t_{1}t_{2}>0,\quad\begin{cases}\left||t_{1}|-|t_{2}|\right|<|\gamma_{1}|+|\gamma_{2}|\\\
||\gamma_{1}|-|\gamma_{2}||<|t_{1}|+|t_{2}|,\end{cases}$ $\displaystyle
t_{1}t_{2}<0,\quad\begin{cases}\left||t_{1}|-|t_{2}|\right|<||\gamma_{1}|-|\gamma_{2}||\\\
|\gamma_{1}|+|\gamma_{2}|<|t_{1}|+|t_{2}|,\end{cases}$
$\displaystyle\phantom{t_{1}t_{2}<0\ }\text{or}\
\begin{cases}\left||t_{1}|-|t_{2}|\right|>||\gamma_{1}|-|\gamma_{2}||\\\
|\gamma_{1}|+|\gamma_{2}|>|t_{1}|+|t_{2}|.\end{cases}$ (2.42)

In the framework of our model, the bond dissipators entail $\gamma_{1}\geq
0,\gamma_{2}\geq 0$. Figure 3 shows the dependence of the spectral winding
number on one of the dissipation strengths with the signs of two symmetric
hopping terms being either the same ($t_{1}=t_{2}$) or the opposite
($t_{1}=-t_{2}$). As soon as $|\nu|\neq 0$, the left and right boundary modes
exponentially pile up at different ends. This unique feature of
$\mathcal{H}_{\text{S}}$ can be applied to the design of the NH topological
sensors exhibiting anomalous sensitivity that grows exponentially with the
system size [38].

FIG. 3: Spectral winding number (gray dots) and localization determinants for
the left and right boundary modes (green and blue lines) and for the bulk
modes (orange curve) as a function of the bond dissipation strength
$\gamma_{2}$ at (a) $t_{1}=1,t_{2}=1,\gamma_{1}=0.2$; (b)
$t_{1}=-1,t_{2}=1,\gamma_{1}=1.5$. For the evaluation of the spectral winding
number $\nu$, we turn the integral in Eq. 2.40 into a discrete summation over
$N=500$ unit cells.

Inside the bulk spectrum, however, it is clear to see that the topological
invariant $\nu$ obtained from the Bloch Hamiltonian fails to locate the
boundary zero modes at $E_{m\neq 0}^{\text{OBC}}=0$ [if compared with Eq.
A13]. In a closed NH system, the concept of conventional bulk-boundary
correspondence has thus been generalized to allow the reconstruction of
topological quantities in the biorthogonal basis [18, 19] such that the
occurrence of the boundary modes are accurately predicted [see also Eq. 2.25].
What happens to an open quantum system? Similarly, its relaxation dynamics
have the Liouvillian skin effect once the NMMs of the Liouvillian pile up
exponentially close to the boundary [60, 66]. It is then crucial to look at
the behavior around the EPs, arising naturally from the Lindblad master
equation [80]:

$\displaystyle t_{i}=\pm\gamma_{i},\qquad i=1,2.$ (2.43)

At EPs, the geometric multiplicity of the Liouvillian is smaller than the
algebraic multiplicity which holds an order that scales with the system size
[16, 20, 1]. We check that for $t_{i}=\pm\gamma_{i}$, the rapidity spectrum
$\beta_{m}^{\text{OBC}}$ in Eqs. (2.23) indeed has one or three eigenvalues.
Approaching one of the EPs, the set of NMMs merges into one or three linearly
independent eigenstates. Moreover, after the mapping of the damping matrix to
$\mathcal{H}_{\text{S}}$ in Eq. 2.8, the adjoint fermions are only permitted
to hop in one direction when $t_{i}=\pm\gamma_{i}$, so all NMMs become exactly
localized at that one end. Consequently, we envision the most drastic
Liouvillian skin effect in close proximity to EPs, which, in terms of the
momentum shift parameter $r$ defined for linking two rapidities in Eq. 2.24,
is manifested as

$\displaystyle
t_{i}\to\pm\gamma_{i}\quad\Leftrightarrow\quad\left|\ln|r^{2}|\right|\gg 0.$
(2.44)

As expected, with a sufficiently large system size, the Liouvillian skin
effect is fully determined by the localization behavior of the bulk modes [see
also Eqs. (A16)$-$(A18)] and the influence of the boundary modes is
negligible. By varying the dissipation strength, alongside the spectral
winding number, Fig. 3 compares the responses in different localization
determinants: $r_{R}^{j}(r_{L}^{j})$, $r^{j}(r^{-j})$ for the piling of the
right (left) boundary and bulk modes at unit cell $j$. To conclude, while not
characterized by the Bloch topological invariant, the Liouvillian skin effect
is embodied in the exceptional topology of the Liouvillian, or, more
precisely, the damping matrix.

####  II.3.2 Liouvillian gap

We go on to study the development of the Liouvillian gap in relation to
various hopping amplitudes and dissipation strengths together with its
response to different boundary conditions. Apparently, the imaginary part of
complex energy $E$ in Eqs. (2.23) is bounded by $\pm(\gamma_{1}+\gamma_{2})$.
Hence, $\Delta=2\min\\{\text{Re}[\beta_{m}]\\}\geq 0$. Figure 4 shows the real
part of the rapidity spectrum as a function of $t_{1}$ under PBC and OBC.

FIG. 4: Real part of the rapidity spectra for $\gamma_{1}=1.5,t_{2}=1$ and
$N=46$. The grey lines show the structure of the periodic system while the
blue and red ones indicate the gap formed by the bulk and edge modes under
OBC. We vary the hopping amplitude $t_{1}$ and set different bond dissipation
strength between unit cells: (a) $\gamma_{2}=0$, (b) $\gamma_{2}=0.5$.

Subjected to a periodic boundary, with finite bond dissipation inside unit
cells $\gamma_{1}\neq 0$, one observes drastically different gap closing
behaviors with and without the second inter-unit-cell dissipation:

$\displaystyle\gamma_{2}$ $\displaystyle=0,\quad\Delta^{\text{PBC}}=0\text{
for }|t_{1}|\leq|t_{2}|;$ $\displaystyle\gamma_{2}$ $\displaystyle\neq
0,\quad\Delta^{\text{PBC}}=0\text{ at }t_{1}=t_{2}.$ (2.45)

We find that after the introduction of the second bond dissipation, the gap
closing line of the Liouvillian discovered in Ref. [60] becomes unstable and
shrinks to a point in the phase diagram. It would lead to a collapse of the
non-trivial (quasi-)NESS and a termination of the steady-state current in a
wide parameter range (see more details in Section III.2).

Meanwhile, the gap-closing points vanish completely when the boundary opens
up. In this circumstance, we resolve the Liouvillian gap from the exact
spectrum in Eqs. (2.23):

$\displaystyle\Delta^{\text{OBC}}=2\gamma-\sum_{i=1,2}2\sqrt{\gamma_{i}^{2}-t_{i}^{2}}\cdot\theta(\gamma_{i}-|t_{i}|).$
(2.46)

Here, $\theta(x)$ denotes the Heaviside step function which is defined as:
$\theta(x)=0$ if $x<0$; $\theta(x)=1$ if $x>0$. Meeting the EPs at
$t_{i}=\pm\gamma_{i}$, a discontinuity appears in the derivatives of the gap
with respect to the symmetric hopping term
$\partial_{t_{i}}\Delta^{\text{OBC}}$ (see Fig. 4) as well as the asymmetric
hopping term $\partial_{\gamma_{i}}\Delta^{\text{OBC}}$ (see Fig. 7).

Lastly, we mention briefly the effects of single-site loss and gain on the
Liouvillian. Given a general set of loss and gain dissipators acting on
individual sites,

$\displaystyle L^{l}_{0,(j,\alpha)}=\sqrt{\gamma^{l}_{0}}a_{j,\alpha},\quad
L^{g}_{0,(j,\alpha)}=\sqrt{\gamma^{g}_{0}}a_{j,\alpha}^{\dagger},$ (2.47)

the only change to the Liouvillian becomes a modified constant term in the
damping matrix that ultimately lifts the minimum of the Liouvillian gap,

$\displaystyle\gamma=\gamma_{0}+\gamma_{1}+\gamma_{2},\quad\Delta\geq
2\gamma_{0},$ (2.48)

with $2\gamma_{0}=|\gamma_{0}^{l}|+|\gamma_{0}^{g}|$. The gap now can no
longer be removed by closing the boundary. Fortunately, the on-site
dissipations will not add any NH term to the damping matrix, thus not altering
the Liouvillian skin effect. Therefore, we implicitly assume
$\gamma_{0}^{l}=\gamma_{0}^{g}=0$.

###  II.4 Correlation function and comparison with the effective Hamiltonian

In this section, we derive a closed form of the single-particle correlation
function from the exact eigenmodes of the damping matrix. The rapidity
spectrum determines the time-dependent part of the two-point correlator. Any
observable consisting of even number fermionic operators can then be
constructed by Wick’s theorem. From the perspective of the correlation
function, we compare the physics of the effective Hamiltonian that neglects
Lindblad quantum jump operators with the picture of the full Liouvillian.
Notably, when the bonds are subjected to purely loss dissipations, the two
mechanisms become identical.

####  II.4.1 Single-particle correlator

Let us start by resolving the time-dependent pairing function for Majorana
fermions. In the equation of motion in Eq. 2.31, the constant matrix $Y$ can
be replaced by the covariance (or steady state) matrix from Eq. 2.28. Through
a change of variable $\tilde{C}(t)=C(t)-C_{\text{ss}}$, one reaches

$\displaystyle\partial_{t}\tilde{C}(t)=-\tilde{C}(t)X-X^{\dagger}\tilde{C}(t).$
(2.49)

Starting from an arbitrary initial configuration that is not trivial
$\tilde{C}(0)\neq 0$, we can integrate the above equation and implement the
diagonalized damping matrix in the exponential:

$\displaystyle
X=\sum_{m}\sum_{\mu=c,d}\beta_{m}|\Theta^{\mu}_{Rm}\rangle\langle\Theta^{\mu}_{Lm}|,$
(2.50)
$\displaystyle|\Theta^{c}_{R(L)m}\rangle=\begin{pmatrix}\underline{\psi}_{R(L)m}\\\
0\end{pmatrix},\quad|\Theta^{d}_{R(L)m}\rangle=\begin{pmatrix}0\\\
\underline{\psi}_{R(L)m}\end{pmatrix}.$

Taking into account the biorthogonality of the basis and the fact that the
damping matrix is real, $X^{*}=X$, we arrive at

$\displaystyle\tilde{C}(t)=\sum_{m,m^{\prime}}\sum_{\mu,\mu^{\prime}}e^{-(\beta_{m}+\beta_{m^{\prime}})t}|\Theta^{\mu^{\prime}*}_{Lm^{\prime}}\rangle\langle\Theta^{\mu^{\prime}*}_{Rm^{\prime}}|\tilde{C}(0)|\Theta^{\mu}_{Rm}\rangle\langle\Theta^{\mu}_{Lm}|.$
(2.51)

At $t=0$, without loss of generality, throughout the text we choose the system
to be in a static configuration with each site completely filled:
$|\Psi_{0}\rangle=\bigotimes_{j=1}^{n_{\text{tot}}}|1\rangle_{j}$, which
corresponds to

$\displaystyle\tilde{C}(0)=-\frac{i(\gamma+\eta)}{\gamma}\begin{pmatrix}0&\mathbb{1}_{n\times
n}\\\ -\mathbb{1}_{n\times n}&0\end{pmatrix}.$ (2.52)

Apparently, $\tilde{C}(0)$ selects $\mu\neq\mu^{\prime}$. With no pairing
between Majorana fermions of the same species in the initial state, $\langle
c_{j}c_{k}\rangle_{t}=\langle d_{j}d_{k}\rangle_{t}=0$.

Next, we go back to the physical space and define the single-particle
correlation function in the spinless fermion language:
$Q_{jk}(t)=\text{Tr}[a_{j}^{\dagger}a_{k}\rho(t)]$. After the mapping to
Majorana fermions in Eq. 2.9, one rewrites it in terms of the pairing
function,

$\displaystyle
Q_{jk}(t)=\frac{i}{4}\sigma(j,k)\left[C_{j,k+n}(t)+C_{k,j+n}(t)\right],$
(2.53)

with $n$ the total number of sites. The phase factor depends on whether the
correlation resides on the same sublattice or not:

$\displaystyle\sigma(j,k)=\begin{cases}1,&j+k=\text{even}\\\
(-1)^{j}\cdot(-i),&j+k=\text{odd}.\end{cases}$ (2.54)

Combined with Eqs. 2.51 and 2.52, the single-particle correlator takes the
explicit form in terms of the exact solutions of the damping matrix in Eqs.
(2.21):

$\displaystyle
Q_{jk}(t)=\left(\frac{\gamma+\eta}{2\gamma}\right)\sigma(j,k)\sum_{m,m^{\prime}}\sum_{l=1}^{n_{\text{tot}}}e^{-(\beta_{m}+\beta_{m^{\prime}})t}$
$\displaystyle\phantom{===}\psi_{Lm}^{*}(j)\psi_{Lm^{\prime}}^{*}(k)\cdot\psi_{Rm}(l)\psi_{Rm^{\prime}}(l).$
(2.55)

####  II.4.2 Effective Hamiltonian in the absence of quantum jumps

For non-Gaussian Lindbladians, on the other hand, it is useful to study the
short-time dynamics by ignoring quantum jumps in the Lindblad equation [55]:
$\sum_{\mu}L_{\mu}\rho L_{\mu}^{\dagger}$. Here, we compare the effective
Hamiltonian description with the full Lindblad master equation framework.

Without quantum jumps, the time evolution of the density can be described by

$\displaystyle\partial_{t}\rho=-i(\mathcal{H}_{\text{eff}}\rho-\rho\mathcal{H}_{\text{eff}}^{\dagger}),$
(2.56)

where $H_{\text{eff}}=H-\frac{i}{2}\sum_{\mu}L_{\mu}^{\dagger}L_{\mu}$. It
turns out that the structure of the effective Hamiltonian becomes drastically
different from the damping matrix [60]. We verify that rather than the total
strength $\gamma$, the non-Hermiticity of $\mathcal{H}_{\text{eff}}$ is
related to the imbalance $\eta$ between loss and gain dissipations:

$\displaystyle\mathcal{H}_{\text{eff}}=$
$\displaystyle\sum_{j=1}^{N-1}(t_{1}-\eta_{1})a_{j,A}^{\dagger}a_{j,B}+(t_{1}+\eta_{1})a_{j,B}^{\dagger}a_{j,A}$
$\displaystyle+(t_{2}-\eta_{2})a_{j,B}^{\dagger}a_{j+1,A}+(t_{2}+\eta_{2})a_{j+1,A}^{\dagger}a_{j,B}$
$\displaystyle-i\eta\sum_{(j,\alpha)}a_{(j,\alpha)}^{\dagger}a_{(j,\alpha)}-is_{0}.$
(2.57)

The purely imaginary energy shift scales with the size of the system:
$s_{0}=[(\gamma-\eta)/2]\cdot
n_{\text{tot}}=\sum_{i=1,2}(|\gamma_{i}^{g}|/2)\cdot n_{\text{tot}}$. It is
convenient to resolve the single-particle correlator directly from
$\mathcal{H}_{\text{eff}}$ according to the equation of motion in Eq. 2.56,

$\displaystyle
Q_{jk,\text{eff}}(t)=\sigma(j,k)\sum_{m,m^{\prime}}\sum_{l=1}^{n_{\text{tot}}}e^{-(\beta_{m,\text{eff}}+\beta_{m^{\prime},\text{eff}})t}$
$\displaystyle\phantom{===}\psi_{Lm}^{\eta*}(j)\psi_{Lm^{\prime}}^{\eta*}(k)\cdot\psi_{Rm}^{\eta}(l)\psi_{Rm^{\prime}}^{\eta}(l),$
(2.58)

with an effective rapidity spectrum:

$\displaystyle\beta_{m,\text{eff}}=\eta+iE_{m}(\eta)+\left(\frac{\gamma-\eta}{2}\right)\cdot
n_{\text{tot}}.$ (2.59)

To work with the same basis as the damping matrix, we have applied the
transformation in Eqs. (2.20). In the rapidity spectrum, $s_{0}$ prevents the
Liouvillian gap from turning negative when
$\eta=\sum_{i=1,2}(|\gamma_{i}^{l}|-|\gamma_{i}^{g}|)/2<0$.

By neglecting the quantum jumps in the Lindblad equation, we find that the
dynamics of the dissipative system are not properly captured by the effective
Hamiltonian at all times. For the short-time interval, regardless of the total
dissipation strength $\gamma$, the imbalance $\eta$ determines the Lindblad
spectrum and makes the Liouvillian gap increase with the system size rather
than stay a constant value as suggested by Eq. 2.16. In the long-time limit,
the EPs of the exact solutions move to $t_{i}=\pm\eta_{i}$, leaving the
Liouvillian skin effect unpredictable compared with Eq. 2.44. Ultimately, the
system always decays to an empty chain as there is no residual matrix $Y$ in
the equation of motion that can add up to a finite stationary occupation of
fermions according to Eqs. (2.31)$-$(2.35).

In spite of all the discrepancies, however, once gain dissipators are
suppressed on the bonds $\gamma_{1}^{g}=\gamma_{2}^{g}=0$, we reach one
special point where

$\displaystyle\gamma=\eta,\quad Q_{jk,\text{eff}}(t)=Q_{jk}(t).$ (2.60)

It infers that for the SSH chain with bond dissipations that only lead to
losses, though truncated, the effective Hamiltonian encapsulates the full
dynamics at arbitrary times.

##  III Non-equilibrium steady states

In Section II, we have already revealed the configuration of the trivial
steady state in Eq. 2.35. It is a unique NESS under OBC and persists in PBC as
long as the Liouvillian gap is not vanishing. With a gapless rapidity
spectrum, NESS becomes degenerate and carries a stationary current that is
independent of the dissipation strengths when the gain and loss contributions
are in balance.

FIG. 5: (a) Absolute value and (b) the imaginary part of the rapidity spectra
as a function of $t_{1}$ for $\gamma_{1}=1.5,\gamma_{2}=0,t_{2}=1$, and
$N=46$. Comparisons are made between the open (blue, red) and periodic (gray)
boundary conditions. For the periodic system, the green dashed line indicates
the emergence of a degenerate NESS at $|t_{1}|=|t_{2}|$ while the intermediate
region $|t_{1}|<|t_{2}|$ supports a quasi-NESS mode with a closed Liouvillian
gap [see Fig. 4(a)] and a non-zero imaginary rapidity.

###  III.1 Open boundary: Uniqueness of NESS

First, we look at the steady state from the OBC rapidity spectrum. As shown in
Fig. 5(a), every NMM has a non-zero rapidity: $\beta_{m}\neq 0$, $\forall m$.
It implies that the density matrix of the steady state
$\rho_{\text{ss}}=\left.e^{\hat{\mathcal{L}}_{+}t}\rho_{0}\right|_{t\to\infty}$
is uniquely determined by the left and right vacua of the Liouvillian in Eq.
2.14: $\langle\text{NESS}^{\prime}|b^{\prime}_{\mu,m}=0$,
$b_{\mu,m}|\text{NESS}\rangle=0$ for $\forall m$ and $\mu\in\\{c,d\\}$. A
proper normalization can be chosen as
$\langle\text{NESS}^{\prime}|\text{NESS}\rangle=\text{Tr}\rho_{\text{ss}}=1$.
Moreover, taking into account $\text{Re}({\beta_{m}})>0$ (see Fig. 4), any
initial state decays to the trivial steady state in Eq. 2.35 with a relaxation
rate proportional to the strictly positive Liouvillian gap in Eq. 2.46. In
particular, at $\gamma_{2}=0$,

$\displaystyle\Delta^{\text{OBC}}=\begin{cases}\
2\gamma_{1}-2\sqrt{\gamma_{1}^{2}-t_{1}^{2}},&|t_{1}|<\gamma_{1}\\\ \
2\gamma_{1},&|t_{1}|\geq\gamma_{1}.\end{cases}$ (3.1)

When $|t_{1}|\geq\gamma_{1}$, all the bulk and edge modes stabilized by the
open boundary share the same Liouvillian gap $2\gamma_{1}$. Whereas for
$|t_{1}|<\gamma_{1}$, the Liouvillian gap decreases but stays positive as long
as $t_{1}\neq 0$. In this regime, the modes with the slowest decay rate appear
at $m=(\nu,q)\in\\{(+,0),(-,\pi)\\}$.

###  III.2 Periodic boundary: Degeneracy, quasi-NESS, and stationary current

From Section II.3, the Liouvillian gap can be closed by switching the boundary
condition to PBC, thus lifting the degeneracy of NESS. For non-zero
$\gamma_{1}$ and $\gamma_{2}$, Fig. 4 shows the real part of the PBC rapidity
spectrum holds a gapless point at $t_{1}=t_{2}$. It is easy to check that the
rapidity of the bulk mode $m^{*}=(+,-\pi)$ vanishes completely:
$\beta_{m^{*}}=0$. Therefore, the right set of the steady states becomes
three-fold degenerate: $s_{R,0}=|\text{NESS}\rangle$,
$s_{R,c}=b^{\prime}_{c,m^{*}}|\text{NESS}\rangle$,
$s_{R,d}=b^{\prime}_{d,m^{*}}|\text{NESS}\rangle$ (accordingly, the left set
is expanded by
$\langle\text{NESS}^{\prime}|,\langle\text{NESS}^{\prime}|b_{c,m^{*}}$ and
$\langle\text{NESS}^{\prime}|b_{d,m^{*}}$). In contrast to OBC, after a long-
time evolution, the final state of the periodic system now depends on the
initial configuration and may appear as a superposition among different NESSs:
$\rho_{\text{ss}}=\rho_{\text{ss}}(\rho_{0})$.

When $\gamma_{2}=0$, from Fig. 5(a), the parameter regime allowing the three-
fold degenerate NESS can be extended to $|t_{1}|=|t_{2}|$ where the zero-
rapidity bulk mode appears at
$m^{*}=(+,-\frac{\pi}{2}[\text{sgn}(t_{1}t_{2})+1])$. Meanwhile, with
$\gamma_{2}$ suppressed, one observes interesting features in the relaxation
behavior reflected by the Liouvillian gap [see also Fig. 4(a)]:

$\displaystyle\Delta^{\text{PBC}}=\begin{cases}\ 0,&|t_{1}|\leq|t_{2}|\\\ \
2\gamma_{1}-2\sqrt{\gamma_{1}^{2}-(|t_{1}|-|t_{2}|)^{2}},&|t_{2}|<|t_{1}|\leq
t_{c}\\\ \
2\gamma_{1}-\frac{2\gamma_{1}|t_{2}|}{\sqrt{t_{1}^{2}-\gamma_{1}^{2}}},&|t_{1}|>t_{c},\end{cases}$
(3.2)

where the critical value is identified as
$t_{c}=(|t_{2}|+\sqrt{t_{2}^{2}+4\gamma_{1}^{2}})/2$. In the region
$|t_{1}|<|t_{2}|$, the Liouvillian gap closes at the bulk modes
$m^{*}=(\pm,\pm\arccos[-t_{1}/t_{2}])$. We call them quasi-steady states,
which are stationary states characteristic of a vanishing decay rate and a
finite phase oscillation frequency as shown in Fig. 5(b):
$|\text{Re}(\beta_{m^{*}})|=0,|\text{Im}(\beta_{m^{*}})|=\sqrt{t_{2}^{2}-t_{1}^{2}}$.
Similar to the degenerate NESS, the final state can also select the quasi-
steady states without any decay in the probability density. Once
$|t_{1}|>|t_{2}|$, the Liouvillian gap of the PBC spectrum opens up but
remains smaller than OBC. The bulk mode that dominates the relaxation process
should be the one with a minimal decay rate. It changes from
$(+,-\frac{\pi}{2}[\text{sgn}(t_{1}t_{2})+1])$ to
$(\pm,\pm\arccos[t_{1}t_{2}/(\gamma_{1}^{2}-t_{1}^{2})])$ when the growing
amplitude $|t_{1}|$ goes past $t_{c}$. It should be noted that the gap
solution in Eq. 3.2 is valid for both the strong ($\gamma_{1}>|t_{2}|$) and
weak ($\gamma_{1}\leq|t_{2}|$) dissipations.

Compared to the trivial steady state in Eq. 2.35, the degenerate and quasi-
NESSs can be viewed as plane waves of fermions with fixed momenta $m^{*}$ on
top of a static uniform occupation. It is rather important to distinguish the
degenerate and quasi-NESSs from the trivial one, especially on account of the
former two being a direct consequence of the closing of the Liouvillian gap.
We find the current flow [72, 81] is such an ideal observable. Defined as
$j_{c}(t)=(\frac{i}{n_{\text{tot}}})\sum_{j}[\langle
a_{j}^{\dagger}a_{j+1}\rangle_{t}-\langle a_{j+1}^{\dagger}a_{j}\rangle_{t}]$,
the time-dependent current flow can be conveniently obtained from the single-
particle correlator in Eq. 2.55:

$\displaystyle j_{c}(t)$
$\displaystyle=\frac{1}{2N}\left(\frac{\gamma+\eta}{\gamma}\right)\sum_{j=1}^{N}\sum_{l=1}^{2N}\sum_{m,m^{\prime}\in\text{all
bands}}e^{-(\beta_{m}+\beta_{m^{\prime}})t}$
$\displaystyle\times\left[-\psi_{Lm}^{*}(2j)\psi_{Lm^{\prime}}^{*}(2j-1)+\psi_{Lm}^{*}(2j)\psi_{Lm^{\prime}}^{*}(2j+1)\right]$
$\displaystyle\times\psi_{Rm}(l)\psi_{Rm^{\prime}}(l).$ (3.3)

Our focal point is to study the behavior of the current in different NESSs. At
larger times, only the Liouvillian gapless modes $m^{*}$ satisfying
$\text{Re}(\beta_{m^{*}})=0$ survive. It leads to

$\displaystyle j_{\text{ss}}$ $\displaystyle\stackrel{{\scriptstyle
t\to\infty}}{{=}}\frac{1}{2N}\left(\frac{\gamma+\eta}{\gamma}\right)\sum_{j=1}^{N}\sum_{l=1}^{2N}\sum_{m,m^{\prime}\in
m^{*}}e^{-i\text{Im}({\beta_{m}}+\beta_{m^{\prime}})t}$
$\displaystyle\times\left[-\psi_{Lm}^{*}(2j)\psi_{Lm^{\prime}}^{*}(2j-1)+\psi_{Lm}^{*}(2j)\psi_{Lm^{\prime}}^{*}(2j+1)\right]$
$\displaystyle\times\psi_{Rm}(l)\psi_{Rm^{\prime}}(l).$ (3.4)

Let us begin with the special limit $\gamma_{2}=0$ where all three types of
NESSs coexist and assume $t_{2}>0$. For the trivial steady state
$m^{*}=\varnothing$, thus the current vanishes in the end

$\displaystyle\left.j_{\text{ss}}\right|_{\gamma_{2}=0}=0,\quad|t_{1}|>t_{2}.$
(3.5)

At phase transition points, the degenerate NESS generated by two bulk modes
$m^{*}=(+,-\frac{\pi}{2}[\text{sgn}(t_{1}t_{2})+1])$ supports a current flow
with an amplitude

$\displaystyle\left.j_{\text{ss}}\right|_{\gamma_{2}=0}=\begin{cases}\
0,&t_{1}=-t_{2};\\\ \
\frac{1}{2N}\left(\frac{\gamma+\eta}{\gamma}\right),&t_{1}=t_{2}.\end{cases}$
(3.6)

In the intermediate gapless region, there emerges quasi-NESS from the modes
$m^{*}=(\pm,\pm\arccos[-t_{1}/t_{2}])$. Analytically, the stationary current
can be expanded in the orders of $1/N$:

$\displaystyle\left.j_{\text{ss}}\right|_{\gamma_{2}=0}\stackrel{{\scriptstyle
t\to\infty}}{{=}}\left(\frac{\gamma+\eta}{\gamma}\right)\left\\{\frac{1}{2N}\left(\frac{t_{1}+t_{2}}{t_{2}}\right)\right.$
$\displaystyle+\frac{1}{2N^{2}}\left.\left[\left(\frac{t_{1}+t_{2}}{t_{2}}\right)\cos(2t_{2}\alpha
t)+\alpha\sin(2t_{2}\alpha t)\right]\right\\}+\cdots,$ (3.7)

with $\alpha=\sqrt{1-(t_{1}/t_{2})^{2}}$. The time dependence in the steady
steady current arises from the finite phase oscillation frequency, a unique
feature possessed by quasi-NESS. Given sufficiently large system size, to the
leading order $O(N^{-1})$, we get

$\displaystyle\left.j_{\text{ss}}\right|_{\gamma_{2}=0}=\frac{1}{2N}\left(\frac{\gamma+\eta}{\gamma}\right)\left(\frac{t_{1}+t_{2}}{t_{2}}\right),\quad|t_{1}|<t_{2}.$
(3.8)

For non-zero $\gamma_{1}$ and $\gamma_{2}$, on the other hand, the steady-
state current is carried by the gapless bulk mode $m^{*}=(+,-\pi)$ at
$t_{1}=t_{2}$, while it vanishes elsewhere. Therefore,

$\displaystyle\left.j_{\text{ss}}\right|_{\gamma_{1}\cdot\gamma_{2}\neq
0}=\frac{1}{2N}\left(\frac{\gamma+\eta}{\gamma}\right)\delta_{t_{1},t_{2}}.$
(3.9)

One immediately notices that when the gain and loss dissipations are in
balance, namely $\eta=0$, the prefactor $(\gamma+\eta)/\gamma\to 1$. The
steady state current is then independent of the dissipation strengths
$\gamma_{i}$. In Fig. 6, we confirm the analytical predictions on the steady
state current for $\gamma_{2}=0$ in Eqs. (3.5)$-$(3.8) by a measurement of
$j_{c}(t)$ at time $\gamma t=10^{5}$. With balanced gain and loss bond
dissipations, we verify the stationary current remains the same under
different values of $\gamma_{1}$. When $\gamma_{2}\neq 0$, in the same manner,
we still find a persistent current at the Liouvillian gap closing point
regardless of the choices of $\gamma_{2}$ [see also Fig. 8(a)].

By contrast, in the description of the effective Hamiltonian, once
$\gamma\neq\eta$ the real part of the rapidity spectrum for a periodic chain
in Eq. 2.59 is always gapped and a persistent current will not be observed in
any allowed parameter regime.

FIG. 6: Steady-state current in the periodic chain with balanced gain and loss
($\eta=0$) as a function of $t_{1}$ for $\gamma_{1}>0,\gamma_{2}=0,t_{2}=1$,
and $N=32$. The dark dots correspond to the stationary current flow
$j_{\text{ss}}=j_{c}(t)$ measured at the time $\gamma_{1}t=10^{5}$. At
$|t_{1}|=t_{2}$, the numerical current matches well with the analytical
prediction [light dots, given by Eq. 3.6]. For $|t_{1}|<t_{2}$, the solid line
indicates the leading order approximation in Eq. 3.7. Outside this region, the
current vanishes due to the opening of the Liouvillian gap.

##  IV Anomalous quantum dynamics

In this section, we search for dynamical signatures of the Liouvillian skin
effect in dissipative quantum systems, originating from the piling up of the
NMMs exponentially close to an open boundary. Compared with previous studies
[60, 63, 67], we show the relaxation behaviors directly obtained from our
exact solutions for an odd number of sites $n=2N-1$ and, at the same time,
include the impact of the second bond dissipators. Apart from a diverging
lifetime without gap closing [66], we find that other global observables such
as a tail of a dynamical current flow and a chiral damping wavefront center
can also serve as good probes of the Liouvillian skin effect, the nature of
which will be related to the exceptional topology of the damping matrix.

###  IV.1 Relaxation of current flow

FIG. 7: The real part of the Lindblad spectrum as a function of the
dissipation strength $\gamma_{2}$ with varied hopping amplitudes: (a)
$t_{1}=t_{2}=1$, (b) $t_{1}=-0.8,t_{2}=1$. We choose $\gamma_{1}=1.5$ and take
$N=12$ unit cells. Fixing $\eta=0$ or $|\gamma_{1}^{l}|=|\gamma_{1}^{g}|$,
$|\gamma_{2}^{l}|=|\gamma_{2}^{g}|$, the lower panel shows the time evolution
of the current circulating along the chain with $\gamma_{1}=1.5$,
$\gamma_{2}=0.5$, $\gamma=\gamma_{1}+\gamma_{2}=2$, and the hopping regimes
are chosen in accordance with the upper panel: (c) $t_{1}=t_{2}=1$, (d)
$t_{1}=-0.8,t_{2}=1$. The green line indicates a steady-state current
$j_{\text{ss}}=0.0417(\simeq 1/24)$ measured at longer times $\gamma
t=10^{4},10^{5}$.

As revealed in Section III.2, for a general set of bond dissipators
($\gamma_{1}\neq 0,\gamma_{2}\neq 0$), a non-vanishing current in the steady
state in Eq. 3.9 helps us to distinguish the degenerate NESS from the trivial
one on a periodic chain. Figure 7 depicts the relaxation process of the
current before reaching the steady state under different boundary conditions.
For PBC, the closing of the Liouvillian gap at $t_{1}=t_{2}$ shown by Fig.
7(a) sustains a stationary current in Fig. 7(c) which, after averaging over
all sites, scales with the inverse of the system size $\sim 1/N$. Deviating
from the gap closing point, the current flow vanishes at relatively short
times. Changing the boundary condition to OBC, the gapless mode in the
rapidity spectrum disappears. All bulk and edge NMMs immediately pile up
exponentially at one of the boundaries, thus terminating the current flow in
Fig. 7(c). If we zoom in to look at the region where the Liouvillian gaps of
two spectra are comparable [for instance, $\gamma_{2}=0.5$ in Fig. 7(b)], the
behaviors of the current flow seen from Fig. 7(d) turn out to be less
sensitive to the boundary conditions.

It infers that from the perspective of the current, the Liouvillian skin
effect is better captured when the gap of the PBC rapidity spectrum is closed.
In Fig. 8, we thus fix $t_{1}=t_{2}$ and study the dynamics of the current
flow occurring at the Liouvillian gap-closing point in a wide range of
dissipation strengths. Consistent with our earlier prediction in Eq. 3.9 for
balanced gain and loss dissipators ($\eta=0$), Fig. 8(a) shows that even in
the presence of very weak dissipations ($\gamma_{1}=0.02,\gamma_{2}=0.01$),
the current flow along a periodic chain saturates to a finite value identical
to the limit of strong dissipations. However, once the boundary opens up,
driven by weak dissipations, the current flow decays with oscillations and has
a much shorter relaxation time as indicated in Fig. 8(b). It also represents
the behaviors of those points far away from EPs ($t_{i}=\pm\gamma_{i}$),
bringing about weak Liouvillian skin effect. For stronger dissipations or
$|\ln(|r^{2}|)|\gg 0$, the amplitude of the current is enlarged as the NMMs
continue to pile up in the same direction. The current flow also shares a
longer relaxation time before the system evolves to the trivial NESS hosting a
static uniform distribution of fermions. The relaxation time or the tail of
the current flow is determined by two factors, the effective Liouvillian gap
$\Delta_{\text{eff}}$ and the correlation length $\xi(r^{2})$ of the system.
It is equal to the lifetime of the particle at one boundary in favour of the
Liouvillian skin effect, of which more detailed analysis can be found in
Section IV.3.

FIG. 8: Relaxation of the electronic current under different boundary
conditions: (a) PBC, (b) OBC. We vary the bond dissipation amplitudes
$\gamma_{1},\gamma_{2}$ while fixing $t_{1}=t_{2}=1$ in a system composed of
$N=12$ unit cells. Given a long-time evolution, the current in a periodic
chain saturates to the steady-state value $j_{\text{ss}}=0.0417(\simeq 1/24)$
regardless of the strengths of $\gamma_{i}$. When the open boundary builds up,
the current vanishes immediately after reaching the peak. Its relaxation time,
equal to the lifetime of the particle at the edge under the Liouvillian skin
effect, depends on the system correlation length $\xi(r^{2})$ and the
effective Liouvillian gap $\Delta_{\text{eff}}$ (4.9).

###  IV.2 Density evolution under damping

Next, we focus on the chiral damping phenomena in the particle-number
distribution of the dissipative chain [60]. We find the damping itself
displays the Liouvillian skin effect beyond the gapless point $t_{1}=t_{2}$.
By turning on the second bond dissipators $\gamma_{2}$, we show the center of
the chiral damping wavefront can be tuned via the localization parameter
$r^{2}$ of the bulk NMMs. It enables us to associate the Liouvillian skin
effect with the existence of the EPs in the damping matrix at which
$|r^{2}|\to 0$ or $\infty$. We further clarify the sensitivity of dynamical
observables to the boundary conditions at longer times. Our results are
consistent with one of the earlier studies [63], demonstrating that probes in
the short-time domain are not sufficient to distinguish the relaxation
processes of open quantum matter under different boundary conditions.

####  IV.2.1 Link with exceptional topology

First, we associate the chiral damping of the particle-number operator with
the exceptional topology of the damping matrix. In the full Lindblad master
equation framework, the evolution of the particle number to NESS,
$\tilde{n}_{j}(t)=n_{j}(t)-n_{j,\text{ss}}$ can be built on our exact single
particle correlator $\langle a_{j}^{\dagger}a_{j}\rangle_{t}$ from Eq. 2.55:

$\displaystyle\tilde{n}_{j}(t)$
$\displaystyle=\left(\frac{\gamma+\eta}{2\gamma}\right)\sum_{m,m^{\prime}}\sum_{l=1}^{n_{\text{tot}}}e^{-(\beta_{m}+\beta_{m^{\prime}})t}$
$\displaystyle\phantom{===}\psi_{Lm}^{*}(j)\psi_{Lm^{\prime}}^{*}(j)\cdot\psi_{Rm}(l)\psi_{Rm^{\prime}}(l).$
(4.1)

Under OBC, a large system size enables us to safely neglect the contributions
from the boundary eigenmodes in Eq. A9. Given exact solutions of the left and
right bulk eigenmodes in Eqs. 2.21, A17 and A18, one derives an asymptotic
scaling of particle number evolution for sites in the odd sublattice $A$:

$\displaystyle\tilde{n}_{2j-1}(t)\sim\frac{1}{r^{2j}}e^{-\Delta_{\text{eff}}t}\cdot\left(\frac{\gamma+\eta}{2\gamma}\right).$
(4.2)

A similar expression can be obtained for the even sublattice $B$. Here, the
damping factor $r^{2}$ comes from the localization of the left bulk
eigenmodes. When $|r^{2}|\to 0(\infty)$, for instance, the left eigenmodes are
localized on the right (left) end of the chain, pushing the chiral wavefront
towards the same boundary while the right eigenmodes center on the opposite.
In the time-dependent part of the particle number operator, we replace the
original gap $\Delta^{\text{OBC}}$ with an effective Liouvillian gap
$\Delta_{\text{eff}}$. The comparison of two quantities is made in Fig. 13.
Since we are interested in the full dynamics of the relaxation process,
instead of taking the long-time limit and keeping track of the slowest
decaying mode, we should include all bulk modes with various decaying rates.
The summation then leads to an effective Liouvillian gap that is also useful
in the evaluation of the particle lifetime in Eq. 4.9.

The two extremities of the Liouvillian skin effect $|r^{2}|\to 0$ or $\infty$
occur at EPs: $t_{i}=\pm\gamma_{i}$, when the adjoint fermions in the damping
matrix are only allowed to hop in one direction after the mapping to
$\mathcal{H}_{\text{S}}$ in Eq. 2.8. It is noteworthy that with a large system
size, the chiral damping condition can be further relaxed to

$\displaystyle|r^{2}|\ll 1\quad\text{or}\quad|r^{2}|\gg 1,$ (4.3)

irrespective of whether the PBC rapidity spectrum is gapless or not.

So far, one discerns that within the framework of the full Lindbladian, the
Liouvillian skin effect is closely related to the exceptional topology of the
non-Bloch damping matrix, whereas it does not depend on the topology of the
Bloch NH Hamiltonian $H_{\text{S}}(q)$ as illustrated in Fig. 3. Meanwhile,
considering $\gamma_{i}$’s and $\eta_{i}$’s belong to two sets of independent
parameters in Eqs. (2.12), the effective Hamiltonian in Eq. 2.57 does not play
a role here neither. The EPs of $\mathcal{H}_{\text{eff}}$ are shifted to
$t_{i}=\pm\eta_{i}$ and the topological regime of ${H}_{\text{eff}}(q)$ for
$\eta_{1}\eta_{2}\geq 0$ lies in Eqs. (2.42) with a replacement:
$\gamma_{i}\to\eta_{i}$.

####  IV.2.2 Modulation of chiral wavefront center

FIG. 9: Polarization of the chiral damping wavefront under OBC. We choose a
system of $N=14$ unit cells. The initial state of the system is completely
filled: $n_{j}(0)=1,\forall j$. At $\eta=0$, the reference steady state
(trivial one) becomes $n_{j,\text{ss}}=0.5,\forall j$. When $|r^{2}|\ll 1$ or
$|r^{2}|\gg 1$, the polarization of $\tilde{n}(t)=n(t)-n_{\text{ss}}$
approaches $1$ or $0$ at longer times. These are cases where the wavefront
terminates at the right or left boundary, leaving the chiral damping behavior
the most distinguished (see also Fig. 10).

We embark on the characterization of the chiral damping behavior emerging in
the dissipative SSH chain when subjected to an open boundary. Let us define
the polarization of the damping process according to

$\displaystyle\Delta
P(t)=\frac{\sum_{j=1}^{n_{\text{tot}}}j\cdot\tilde{n}_{j}(t)}{{n_{\text{tot}}}\cdot\sum_{j=1}^{n_{\text{tot}}}\tilde{n}_{j}(t)},$
(4.4)

where, as depicted in Fig. 1, the length of chain is chosen to be odd
${n_{\text{tot}}}=2N-1$. Without loss of generality, we restrict ourselves to
the case when the gain and loss dissipations are in balance:
$\eta_{1}=\eta_{2}=0$. Starting from a completely filled chain, the NESS
should then be half filling at each site in Eq. 2.34. Plugging the analytical
expression for particle number in Eq. 4.1, the tendencies of polarization
under different hopping amplitudes and dissipation strengths are shown in Fig.
9. Figure 10 further illustrates the motion of the chiral damping wavefront
along the chain as time goes by.

FIG. 10: Dynamics of the particle-number operator
$\tilde{n}_{j}(t)={n}_{j}(t)-n_{j,\text{ss}}$ under OBC. We keep the same set
of Lindblad bond dissipators and the same initial condition as Fig. 9. The
left column corresponds to the upper three curves in Fig. 9 associated with
$|r^{2}|$: (a) $0.074$, (c) $0.60$, (e) $0.94$, and the right column produces
the damping of the lower three curves linked with $|r^{2}|$: (b) $0.99$, (d)
$1.8$, (f) $4.5$. The wavefront is visible in the square root of the relative
density.

At the initial time, $\Delta P=1/2$, it refers to a uniform particle number
distribution consistent with our initial condition. Around $|r^{2}|=1$, when
the Liouvillian skin effect is quite weak, the polarization stays close to
$1/2$ and the chiral damping wavefront is absent in Figs. 10(e) and 10(b). As
$|r^{2}|\ll 1$ ($\gg 1$), the NMMs of the Liouvillian start to pile up towards
the right (left) end such that the particles closer to that boundary are
granted a longer lifetime, thus pushing the polarization to $1$ $(0)$. From
Figs. 10(a) and 10(f), the chiral wavefront also appears most distinguished in
these two limits. In the intermediate parameter regime of $|r^{2}|$, the
polarization evolves to a finite value in between $[0,1/2]$ or $[1/2,1]$ in
accompaniment with a damping wavefront growing obscure, as shown by Figs.
10(d) and 10(c). Therefore, by tuning the parameter $r^{2}$, we are able to
modulate the polarization or the center of the chiral damping wavefront
regardless of the topology of the PBC rapidity spectrum before the boundary
opens up. The left and right columns of Fig. 10 correspond to the gap closing
($t_{1}=t_{2}$) and gap opening ($t_{1}\neq t_{2}$) points in the PBC
spectrum, respectively.

On a side note, in Fig. 9 even at relatively large times, $\Delta P$ will not
converge to $1/2$ for its measurement is targeted on the excess of the
particle number over a uniformly half-filled NESS, rather than the real
occupation number.

####  IV.2.3 Boundary sensitivity

From the chiral damping phenomena, the Liouvillian skin effect takes place
upon changing the boundary from PBC to OBC. Since the system evolves from a
completely filled initial state to a uniform steady state, the damping serves
as a global effect. Nevertheless, we can also resolve the motion of a single
particle and study the sensitivity of the open system to boundary conditions
on broader grounds.

It is argued in Ref. [63] that in the thermodynamical limit, the single-
particle Green’s function of NH systems is independent of the boundary
conditions. We show in the following within the Lindblad master equation
framework, the argument is indeed true for the evolution of the density matrix
at short times. However, as soon as the motion of a particle involves the
edge, due to the gap closing and the Liouvillian skin effect, there emerge
drastic differences in the relaxation process under PBC or OBC.

Let us start by putting one particle in the middle of the chain:
$n_{j}(0)=\delta_{j,N}$. The initial condition corresponds to a Majorana
pairing configuration,

$\displaystyle\tilde{C}(0)=-\frac{i(\eta-\gamma)}{\gamma}\begin{pmatrix}0&\mathbb{1}_{n\times
n}\\\ -\mathbb{1}_{n\times n}&0\end{pmatrix}-i\begin{pmatrix}0&D\\\
-D&0\end{pmatrix},$ (4.5)

where the $n\times n$ matrix $D$ only holds one non-zero element
$D_{(N,N)}=1$. The time evolution of the particle number at the $j$th site can
be described by

$\displaystyle\tilde{n}_{j}(t)$
$\displaystyle=\sum_{m,m^{\prime}}e^{-(\beta_{m}+\beta_{m^{\prime}})t}\psi_{Lm}^{*}(j)\psi_{Lm^{\prime}}^{*}(j)$
$\displaystyle\left[\psi_{Rm}(N)\psi_{Rm^{\prime}}(N)+\left(\frac{\eta-\gamma}{2\gamma}\right)\sum_{l=1}^{n_{\text{tot}}}\psi_{Rm}(l)\psi_{Rm^{\prime}}(l)\right].$
(4.6)

To make a direct comparison with Ref. [63], we suppress the gain dissipations
on both bonds: $\gamma_{1}^{g}=\gamma_{2}^{g}=0$ without violating the
solvable constraint in Eqs. (2.30):
${\gamma_{1}}/{\eta_{1}}={\gamma_{2}}/{\eta_{2}}=1$. Since $\gamma=\eta$, the
second term in the particle-number operator vanishes, and the trivial NESS
refers to an empty chain indicated by Eq. 2.35.

Also, we focus on the relaxation process at the same point $t_{1}=t_{2}$ as
Ref. [63], where the gap of the PBC rapidity spectrum closes for arbitrary
dissipation strengths $\gamma_{1}$ and $\gamma_{2}$. The inclusion of
$\gamma_{2}\neq 0$ offers more tunability on the Liouvillian skin effect
parameter $r^{2}$. As shown in Fig. 11, the motion of the particle does not
depend on the boundary condition until its trajectory hits the edge. Once that
long-time evolution is permitted, contrary to Ref. [63], the motion differs a
lot under PBC or OBC. We observe a persistent current $j_{\text{ss}}=1/N$
circulating along the periodic chain while the current quickly terminates
under the influence of an open boundary.

Meanwhile, we compare the responses of the particle motion to the presence of
a strong or weak Liouvillian skin effect. When the skin effect dominates
[$|r^{2}|\to 0$ in Figs. 11(a) and 11(b)], the wave packet becomes less
dispersive and oriented towards the right end under both PBC and OBC. The
current also vanishes completely at the open boundary without any reflection.
Both features are in marked contrast to the case of the weak Liouvillian skin
effect [$|r^{2}|\to 1$ shown by Figs. 11(c) and 11(d)].

One concludes that regardless of the initial conditions in the Lindblad setup,
the Liouvillian skin effect is visible and manifested in the evolution of the
density matrix given a probe of open quantum matter at longer times.

FIG. 11: Boundary effects on the motion of a fermionic particle. At the
beginning, we put one particle in the middle of the chain: $n_{N}(0)=1$ with
$N=14$ and $n_{j}(0)=0$, elsewhere. With gain dissipators suppressed,
$\gamma_{1}^{g}=\gamma_{2}^{g}=0$,
$\gamma=(|\gamma_{1}^{l}|+|\gamma_{2}^{l}|)/2=\eta$. The reference steady
state (trivial one) corresponds to an empty chain: $n_{j,\text{ss}}=0,\forall
j$. The top panel shows the case of the strong Liouvillian skin effect at
$t_{1}=t_{2}=1,\gamma_{1}=0.2,\gamma_{2}=0.8,|r^{2}|=0.074$ under boundary
conditions: (a) PBC and (b) OBC. The lower panel compares the region where the
Liouvillian skin effect becomes weaker:
$t_{1}=t_{2}=1,\gamma_{1}=0.1,\gamma_{2}=0.05,|r^{2}|=0.74$ under boundary
conditions: (c) PBC and (d) OBC.

###  IV.3 Lifetime of non-equilibrium particle at the edge

In the last section, we resolved the relaxation time of dynamical observables
under the Liouvillian skin effect. For a measurement, any local observable of
the operator $\hat{Q}$ has reached the steady state if
$|\tilde{Q}(t)|=|Q(t)-Q_{\text{ss}}|\ll||\hat{Q}||=(\text{Tr}[\hat{Q}^{\dagger}\hat{Q}])^{1/2}$.
Along the lines of Ref. [66], one can thus define the maximal relaxation time
$\tau$ of the system as that of the slowest decaying mode:
$|c|e^{-\tau\Delta}=e^{-1}\ll 1$ with $\Delta$ the Liouvillian gap. The
presence of the Liouvillian skin effect enhances the amplitude of this mode,
showing an exponential localization tendency near the boundary: $|c|\sim
e^{\mathcal{O}(L/\xi)}$. Here, $\xi$ denotes the correlation length of the
system. It leads to

$\displaystyle\tau\sim\frac{1}{\Delta}+\frac{1}{\xi}\cdot\frac{L}{\Delta}.$
(4.7)

In our dissipative SSH model, we look at the Liouvillian skin effect region
$|r^{2}|<1$ and take the dynamical observable as the occupation of a particle
residing at the right boundary. More importantly, we make fewer assumptions
than Ref. [66] by defining the lifetime of the non-equilibrium particle
according to

$\displaystyle|\tilde{n}_{2N-1}(\tau)|=e^{-l}|\tilde{n}_{2N-1}(0)|,$ (4.8)

where $l$ is taken as a positive integer. The precision of the lifetime can be
improved by increasing the value of $l$. The definition in Eq. 4.8 has the
advantage of including the contributions from all decaying modes, staying
closer to a real measurement.

For simplicity, we start with a completely filled chain and balanced gain and
loss dissipators $\eta_{1}=\eta_{2}=0$. Then,
$|\tilde{n}_{2N-1}(0)|=1-1/2=1/2$. From the asymptotic scaling of the particle
number evolution in Eq. 4.2, we establish that

$\displaystyle\tau_{2N-1}\sim\frac{1}{\Delta_{\text{eff}}}+\ln(r^{-2})\cdot\frac{N}{\Delta_{\text{eff}}}.$
(4.9)

The correlation length of the system in Eq. 4.7 thus satisfies

$\displaystyle\xi(r^{2})\simeq\frac{2}{|\ln(|r^{2}|)|},$ (4.10)

where the factor of $2$ comes from the identification $L=2N-1$. Consistent
with Ref. [66], we find the lifetime of the particle at the edge grows
linearly with the system size without a closing of the effective Liouvillian
gap.

Next, based on our exact solutions, we verify numerically the above relations
by varying the system size and the skin effect parameter $r^{2}$. It is also
important to explore two limits, the weak and strong dissipations where the
formation of the effective Liouvillian gap differs.

At weak dissipations $\gamma_{1}<|t_{1}|$ and $\gamma_{2}<|t_{2}|$, all bulk
and boundary modes share the same Liouvillian gap [see Eqs. (2.23) and Fig.
4]. Hence, $\Delta_{\text{eff}}=\Delta^{\text{OBC}}=2(\gamma_{1}+\gamma_{2})$.
Figure 12(a) plots the numerical scaling of the particle lifetime at the edge
$\tau_{2N-1}=\tau$ under different system sizes at $\gamma_{2}=0$. The linear
dependence on $N$ in Eq. 4.9 becomes more visible as the length of the chain
increases. With stronger Liouvillian skin effect $|r^{2}|\ll 1$, the linear
scaling holds true for a relatively small system size $N\sim 8$. The system
correlation length can be extracted from the slope and Fig. 12(b) confirms our
analytical prediction in terms of the skin effect parameter in Eq. 4.10.
Approaching the EP, $\gamma_{1}\to t_{1}=1$, the small discrepancy results
from the fact that the strong polarization of the bulk eigenmodes challenges
the numerical precision of the damping matrix decomposition in Eqs. (2.50).

FIG. 12: (a) Lifetime of the out-of-equilibrium mode on the right boundary as
a function of the system length $L=2N-1$. We take
$|\tilde{n}_{2N-1}(\tau)|=e^{-2}|\tilde{n}_{2N-1}(0)|=0.068$. (b) Correlation
length as a function of the bond dissipation strength $\gamma_{1}$. To extract
the correlation length $\xi$ in the linear relation in Eq. 4.7, a finite
system size in the range of $8$ to $16$ unit cells is employed. Meanwhile,
$\tau$ is chosen at the value
$|\tilde{n}_{2N-1}(\tau)|=e^{-3}|\tilde{n}_{2N-1}(0)|=0.025$. For both plots,
the initial conditions are kept the same as Fig. 9: the filling of each site
evolves from $1$ to $1/2$. The remaining parameters are fixed accordingly:
$t_{1}=t_{2}=1$ and $\gamma_{2}=0$.

At strong dissipations, on the other hand, from Fig. 4 different bulk and
boundary modes form distinct real rapidity spectra. As a consequence, the
effective Liouvillian gap incorporating effects of all decaying modes differs
from $\Delta^{\text{OBC}}$ in Eq. 2.46. By definition,

$\displaystyle\Delta_{\text{eff}}$
$\displaystyle=\Delta^{\text{OBC}},\quad\gamma_{1}<|t_{1}|,\
\gamma_{2}<|t_{2}|;$ $\displaystyle\Delta_{\text{eff}}$
$\displaystyle>\Delta^{\text{OBC}},\quad\text{ otherwise }.$ (4.11)

Still, we can estimate its value from the established linear relation for the
particle lifetime in Eq. 4.9 and at the same time, fix the correlation length
by the skin effect parameter $r^{2}$ in Eq. 4.10. In Fig. 13, we show the
formation of the effective Liouvillian gap beyond weak dissipations:
$\gamma_{1}>|t_{1}|$, $\gamma_{2}<|t_{2}|$. Indeed, $\Delta^{\text{OBC}}$
provides a lower bound for the effective Liouvillian gap.

On the contrary, when $\gamma\neq\eta$, the lifetime of a particle at the
right boundary under the truncated $\mathcal{H}_{\text{eff}}$ would scale as

$\displaystyle\left.\tau_{2N-1}\right|_{\mathcal{H}_{\text{eff}}}\sim\frac{1}{\Delta}\simeq\frac{1}{2(\gamma-\eta)}\cdot\frac{1}{N}$
(4.12)

in the large $N$ limit. Without the quantum jumps, the lifetime decreases with
the system size. As expected, the effective Hamiltonian becomes problematic in
describing the long-time dynamics.

FIG. 13: Effective Liouvillian gap (green dots) as a function of $\gamma_{1}$
beyond weak dissipations: $\gamma_{1}>|t_{1}|$, $\gamma_{2}<|t_{2}|$. We take
$t_{1}=t_{2}=1,\gamma_{2}=0.2$, and
$\tau\to|\tilde{n}_{2N-1}(\tau)|=e^{-5}|\tilde{n}_{2N-1}(0)|=0.0034$. The
fitting of the lifetime $\tau$ in Eq. 4.9 is performed with a finite system
size $N=8\sim 14$. The inset shows its linear dependence at
$\gamma_{1}=2.5,|r^{2}|=0.286$. The background corresponds to the real part of
the rapidity spectrum generated by the bulk modes (blue curves) and boundary
modes (red curve) in a system of $N=14$ unit cells. $\Delta_{\text{OBC}}$
(black curve) provides a lower bound for the effective Liouvillian gap.

##  V Discussion

We have studied the Liouvillian skin effect in a bond-dissipative SSH model,
of which the rapidity spectrum and NMMs are exactly solvable. This has
illuminated the relation between the NH skin effect, two different effective
Hamiltonians, $H_{\text{eff}}$, $H_{\text{S}}$, and the full quantum master
equation description. A number of dynamical phenomena, such as diverging
relaxation times, inherited in the quantum setup originate from the anomalous
boundary sensitivity of the NH Hamiltonians and have been investigated in
detail. This paper also provides a solid platform to resolve the entanglement
spectrum and identify transitions between different NH topological phases [82]
whose experimental realization is feasible, e.g., in ultra-cold atoms with a
momentum lattice [59].

Dynamical probes of the quantum Fisher information [47] and the application as
NH topological sensors [46] provide intriguing questions for future
investigation. Meanwhile, the current model reveals the Liouvillian skin
effect in the relaxation process while posing the question on the search for
steady states that inherit the exceptional topology of the NH damping matrix.
It is found that using monitored quantum circuits, periodic measurement allows
access to biorthogonal steady state observables [29]. Future works could also
be centered around the interaction effects on generic Lindbladians. For
instance, the dynamical mean-field theory can be implemented to uncover the
interplay of dissipation and environmental fluctuations in the open quantum
matter that consists of interacting fermions [79] or bosons [83]. A mixture of
particle species is a good path to the discovery of further variations of the
Liouvillian skin effect. From a broader theoretical point of view, effective
field theories with bosonization prove powerful in studying dissipative
Luttinger liquids [84]. In presence of weak symmetries, one may tackle even
strongly interacting Liouvillians [85]. The inclusion of disorder and chaos
may be possible provided the solvable limit in a disorder-averaging SYK
Lindbladian [86].

## Acknowledgements

This work has benefited from useful discussions with Elisabet Edvardsson,
Yuchi He, Lukas König, Daniel Varjas, and Kang Yang. The authors were
supported by the Swedish Research Council (VR) and the Wallenberg Academy
Fellows program as well as the project Dynamic Quantum Matter of the Knut and
Alice Wallenberg Foundation. Q.D.J was also supported by Pujiang Talent
Program No. 21PJ1405400.

## APPENDIX: EXACT SOLUTION OF THE GENERALIZED NON-HERMITIAN SSH CHAIN

In the Appendix, we present a detailed derivation on the exact solution of the
NH SSH chain in Eq. 2.8. The generic results for the simplified model [18, 34]
are extended here to allow two asymmetric hopping terms.

###  .1 Periodic boundary condition

First, we address the case of a periodic boundary:
$\varphi_{N+1,A(B)}=\varphi_{1,A(B)}$. For convenience, we set the lattice
spacing to unity. Through the Fourier transform
$\varphi_{{j,\alpha}}=\frac{1}{\sqrt{N}}\sum_{q}e^{iqj}\varphi_{\alpha}(q)$,
with $q=2\pi m^{\prime}/N$, $m^{\prime}=-N/2,-N/2+1,\dots,0,\dots,N/2-1$, the
anti-commutation relations for fermionic operators are respected:
$\\{\varphi_{\alpha}(q),\varphi_{\alpha^{\prime}}^{\dagger}(q^{\prime})\\}=\delta_{q,q^{\prime}}\delta_{\alpha,\alpha^{\prime}}$.
In the new basis $\underline{\varphi}(q)=(\varphi_{A}(q),\varphi_{B}(q))^{T}$,
one rewrites the Hamiltonian into

$\displaystyle\mathcal{H}_{\text{S}}=\sum_{q}\underline{\varphi}^{\dagger}(q)\cdot
H_{\text{S}}(q)\underline{\varphi}(q).$ (A1)

In terms of Pauli matrices, the matrix elements read

$\displaystyle H_{\text{S}}(q)$ $\displaystyle=\vec{h}(q)\cdot\vec{\sigma},$
$\displaystyle h_{x}(q)$
$\displaystyle=t_{1}+t_{2}\cos(q)+i\gamma_{2}\sin(q),$ $\displaystyle
h_{y}(q)$ $\displaystyle=i\gamma_{1}-i\gamma_{2}\cos(q)+t_{2}\sin(q),$
$\displaystyle h_{z}(q)$ $\displaystyle=0.$ (A2)

From the eigenvalue Eqs. (2.19), we solve the energy spectrum,
$E_{m}^{\text{PBC}}=E^{\text{PBC}}_{\nu=\pm}(q)$,

$\displaystyle E^{\text{PBC}}_{\pm}(q)$
$\displaystyle=\pm\left[t_{1}^{2}+t_{2}^{2}-(\gamma_{1}^{2}+\gamma_{2}^{2})+2(t_{1}t_{2}+\gamma_{1}\gamma_{2})\cos
q\right.$
$\displaystyle\phantom{=}\left.+2i(t_{1}\gamma_{2}+t_{2}\gamma_{1})\sin
q\right]^{\frac{1}{2}},$ (A3)

together with the right and left eigenvectors:

$\displaystyle\tilde{\underline{\psi}}_{R,\nu}(q,\gamma_{i})$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}(t_{1}+\gamma_{1}+(t_{2}-\gamma_{2})e^{-iq})/E_{\nu}(q)\\\
1\end{pmatrix},$ $\displaystyle\tilde{\underline{\psi}}_{L,\nu}(q,\gamma_{i})$
$\displaystyle=\tilde{\underline{\psi}}_{R,\nu}(q,-\gamma_{i}).$ (A4)

The second mutual relation can be proved by taking into account
$H_{\text{S}}^{\dagger}(q,\gamma_{i})=H_{\text{S}}(q,-\gamma_{i})$ and
$E_{\nu}^{*}(q,\gamma_{i})=E_{\nu}(q,-\gamma_{i})|_{\text{PBC}}$. The
normalization of eigenvectors is chosen according to the biorthogonality:

$\displaystyle\tilde{\underline{\psi}}_{L,\nu}^{*}(q)\cdot\tilde{\underline{\psi}}_{R,\nu^{\prime}}(q)=\delta_{\nu,\nu^{\prime}}.$
(A5)

The Hamiltonian now shares a diagonalized structure:

$\displaystyle
H_{\text{S}}=\sum_{(\nu,q)}\underline{\varphi}^{\dagger}(q)\cdot\left[E_{\nu}(q)\tilde{\underline{\psi}}_{R,\nu}(q)\cdot\tilde{\underline{\psi}}^{*}_{L,\nu}(q)\right]\underline{\varphi}(q).$
(A6)

For the evaluation of the single-particle correlation acting on the real space
in Eq. 2.55, the components $\tilde{\psi}_{R,(\nu,q)}(j,\alpha)$,
$\tilde{\psi}^{*}_{L,(\nu,q)}(j,\alpha)$ are obtained by the Fourier
transform:

$\displaystyle\begin{pmatrix}\tilde{\psi}_{R,(\nu,q)}(j,A)\\\
\tilde{\psi}_{R,(\nu,q)}(j,B)\end{pmatrix}$
$\displaystyle=\frac{e^{iqj}}{\sqrt{N}}\ \tilde{\underline{\psi}}_{R,\nu}(q),$
$\displaystyle\begin{pmatrix}\tilde{\psi}^{*}_{L,(\nu,q)}(j,A)\\\
\tilde{\psi}^{*}_{L,(\nu,q)}(j,B)\end{pmatrix}$
$\displaystyle=\frac{e^{-iqj}}{\sqrt{N}}\
\tilde{\underline{\psi}}^{*}_{L,\nu}(q),$ (A7)

where we have used $\langle
0|\varphi_{\alpha^{\prime}}(q^{\prime})\varphi^{\dagger}_{\alpha}(q)|0\rangle=\delta_{\alpha,\alpha^{\prime}}\delta_{q,q^{\prime}}$.

###  .2 Open boundary condition

Under OBC, we break the $N$th unit cell by taking away the last $B$ site. The
NH SSH Hamiltonian is then expressed in an explicit $n\times n$ matrix form:
$\mathcal{H}_{\text{S}}=\underline{\varphi}^{\dagger}\cdot
H_{\text{S}}\underline{\varphi}$ with
$\varphi=(\varphi_{(1,A)},\varphi_{(1,B)},\dots,\varphi_{(N,A)})^{T}$ and

$\displaystyle\ \ H_{\text{S}}=$
$\displaystyle\begin{pmatrix}0&t_{1}+\gamma_{1}&&&&&\\\
t_{1}-\gamma_{1}&0&t_{2}+\gamma_{2}&&&&&&\\\ &t_{2}-\gamma_{2}&0&&&&&\\\
&&&&\ddots&&\\\ &&&&&0&t_{2}+\gamma_{2}\\\
&&&&&t_{2}-\gamma_{2}&0\end{pmatrix}.$ (A8)

For $n=2N-1$ an odd number of sites, $H_{\text{S}}$ is shown to be exactly
solvable at $\gamma_{1}\neq 0$, $\gamma_{2}=0$ [18, 34]. Below, we generalize
the exact OBC solution to the region: $\gamma_{1}\neq 0$, $\gamma_{2}\neq 0$.
The spectrum recovers the known solution by setting $\gamma_{2}\to 0$ and
retains the main features in the new limit.

On one hand, there exist two zero-energy boundary states that are fully
suppressed on the sublattice $B$,

$\displaystyle
E_{0}=0:\quad\tilde{\underline{\psi}}_{R0}=\mathcal{N}_{R}\begin{pmatrix}r_{R}\\\
0\\\ r_{R}^{2}\\\ 0\\\ \dots\\\ 0\\\ r_{R}^{N}\end{pmatrix},\
\tilde{\underline{\psi}}_{L0}=\mathcal{N}_{L}\begin{pmatrix}r_{L}\\\ 0\\\
r_{L}^{2}\\\ 0\\\ \dots\\\ 0\\\ r_{L}^{N}\end{pmatrix},$ (A9)

with associated parameters $r_{R}=-(t_{1}-\gamma_{1})/(t_{2}+\gamma_{2})$,
$r_{L}=-(t_{1}+\gamma_{1})/(t_{2}-\gamma_{2})$. The two eigenstates are
biorthogonal to each other,

$\displaystyle\tilde{\underline{\psi}}_{L0}^{*}\cdot\tilde{\underline{\psi}}_{R0}=1,$
(A10)

which brings about a constraint on the normalization factors:
$\mathcal{N}_{L}^{*}\mathcal{N}_{R}={(1-r_{L}^{*}r_{R})}/[{r_{L}^{*}r_{R}(1-(r_{L}^{*}r_{R})^{N})]}$.
It is easy to discern that, depending on the absolute values of $r_{L}$ and
$r_{R}$, the left and right boundary modes can be localized at different ends:
$\text{sgn}[\ln(|r_{L}|)]\neq\text{sgn}[\ln(|r_{R}|)]$. We find the parameters
satisfying this condition all stay in the topological regime in Eqs. (2.42)
where the spectral winding number is non-trivial.

On the other hand, to derive the bulk spectrum, we present an intuitive
approach by analogy to Refs. [19, 34]. Our first step is the identification of
the gap-closing points when the boundary is opened. One way is to find a
transformation matrix $R$ such that the NH SSH Hamiltonian is mapped to its
Hermitian counterpart $\bar{H}_{\text{S}}$:

$\displaystyle\bar{H}_{\text{S}}=R^{-1}H_{\text{S}}R.$ (A11)

A proper construction leads to $R=R_{1}R_{2}$, where
$R_{1}=\text{diag}\\{1,r_{1},r_{1},r_{1}^{2},r_{1}^{2},\dots,r_{1}^{N-1},r_{1}^{N-1}\\}$,
$R_{2}=\text{diag}\\{1,1,r_{2},r_{2},r_{2}^{2},r_{2}^{2},\dots,r_{2}^{N-2},r_{2}^{N-1}\\}$,
and $r_{1}=\sqrt{(t_{1}-\gamma_{1})/(t_{1}+\gamma_{1})}$,
$r_{2}=\sqrt{(t_{2}-\gamma_{2})/(t_{2}+\gamma_{2})}$. The Hermitian SSH chain
is embedded with the anisotropic hopping strengths

$\displaystyle\bar{t}_{1}=\sqrt{t_{1}^{2}-\gamma_{1}^{2}},\quad\bar{t}_{2}=\sqrt{t_{2}^{2}-\gamma_{2}^{2}},$
(A12)

with the gapless phase transition occurring at $|\bar{t}_{1}|=|\bar{t}_{2}|$.
$H_{\text{S}}$ inherits from the transformation these gap closing points:

$\displaystyle|t_{1}^{2}-\gamma_{1}^{2}|=|t_{2}^{2}-\gamma_{2}^{2}|.$ (A13)

An alternative method of reproducing the gapless points in the OBC spectrum is
to study the biorthogonal polarization of the NH system [18, 34], which
changes its integer value at $|r_{L}^{*}r_{R}|=1$, in consistency with our
result in Eq. A13.

Under the transformation in Eq. A11, one can also link the bulk eigenstates of
$H_{\text{S}}$ to the ones of $\bar{H}_{\text{S}}$ by
$\tilde{\underline{\psi}}_{R}=R\bar{\underline{\psi}}_{R}$ and
$\tilde{\underline{\psi}}_{L}=\bar{\underline{\psi}}_{L}R^{-1}$, which
indicates the piling up of the right (left) states at one end with an
exponential localization factor $r^{j}$ ($r^{-j}$), where
$r=r_{1}r_{2}=\sqrt{r_{R}/r_{L}}$. It is equivalent to a change of momentum in
the Bloch phase factor $e^{iqj}$ from $q$ to $q-i\ln(r)$, thus allowing us to
build the OBC bulk spectrum from the PBC one [19, 20]. More precisely, at a
fixed momentum $q$, like the PBC spectrum the bulk energies always come in
pairs ($\nu=\pm$): $E_{m\neq 0}=E_{\nu}^{\text{OBC}}(q)$. The two spectra are
related by

$\displaystyle E_{\pm}^{\text{OBC}}(q)=E_{\pm}^{\text{PBC}}(q-i\ln(r))$
$\displaystyle=\pm\sqrt{t_{1}^{2}+t_{2}^{2}-(\gamma_{1}^{2}+\gamma_{2}^{2})+2\sqrt{(t_{1}^{2}-\gamma_{1}^{2})(t_{2}^{2}-\gamma_{2}^{2})}\cos(q)}.$
(A14)

Given an even number of sites $n=2N$, the OBC spectrum obtained in this way
applies only to the large $N$ limit [19]. For an odd number of sites $n=2N-1$,
with additional mirror symmetry [87, 34], it can be shown that the bulk
spectrum in Eq. A14 becomes exact for any finite $N$.

Let us verify by building the exact left and right eigenvectors. We choose
$q=\pi m^{\prime}/N$, $m^{\prime}=1,\dots,N-1$, so that in total the index
$m\in\\{0,(\pm,q)\\}$ reproduces a complete set of $2N-1$ bands. A trial
function for the right eigenstates can be written as

$\displaystyle\tilde{\underline{\psi}}_{R\nu}(q,\gamma_{i})=\frac{1}{\sqrt{2N}}\begin{pmatrix}\tilde{\underline{\psi}}_{R\nu}(q,1)\\\
\tilde{\underline{\psi}}_{R\nu}(q,2)\\\ \cdots\\\
\tilde{\underline{\psi}}_{R\nu}(q,N)\end{pmatrix},$ (A15)

where the component in the $j$th unit cell should be linked with the PBC
eigenvector in Eq. A4 by the same momentum shift:

$\displaystyle\tilde{\underline{\psi}}_{R\nu}(q,j)\sim
r^{j}e^{iqj}\tilde{\underline{\psi}}_{R\nu}^{\text{PBC}}(q-i\ln(r)).$ (A16)

On top of that, mirror symmetry in the spectrum
$E_{\pm}^{\text{OBC}}(q)=E_{\pm}^{\text{OBC}}(-q)$ enforces a superposition of
two wave functions with opposite momenta. The relative amplitude of this
superposition is determined by the boundary condition of a broken last unit
cell:
$|\tilde{\underline{\psi}}_{R\nu}(q,0)\rangle_{B}=|\tilde{\underline{\psi}}_{R\nu}(q,N)\rangle_{B}=0$.
One arrives at

$\displaystyle\tilde{\underline{\psi}}_{R\nu}(q,j)$
$\displaystyle=r^{j}\left(e^{iqj}\underline{u}_{R\nu}(q)-e^{-iqj}\underline{u}_{R\nu}(-q)\right),$
(A17) $\displaystyle\underline{u}_{R\nu}(q)$
$\displaystyle=\frac{i}{\sqrt{2}}\begin{pmatrix}(t_{1}+\gamma_{1}+(t_{2}-\gamma_{2})r^{-1}e^{-iq})/E_{\nu}(q)\\\
1\end{pmatrix}.$

Taking into account
$H_{\text{S}}^{\dagger}(\gamma)=H_{\text{S}}^{T}(\gamma)=H_{\text{S}}(-\gamma)$
and $E_{\nu}(q,\gamma_{i})=E_{\nu}(q,-\gamma_{i})|_{\text{OBC}}$, the left and
right eigenvectors are related in a way slightly different from PBC:

$\displaystyle\tilde{\underline{\psi}}_{L\nu}(q,\gamma_{i})=\tilde{\underline{\psi}}_{R\nu}^{*}(q,-\gamma_{i}).$
(A18)

It is straightforward to check that the energies in Eq. A14 and eigenstates in
Eqs. (A15)$-$(A18) satisfy the eigenvalue Eqs. (2.19) for any finite system
size $N$. The biorthogonal normalization conditions among the bulk and the
boundary states are respected as well:

$\displaystyle\tilde{\underline{\psi}}^{*}_{L\nu}(q)\cdot\tilde{\underline{\psi}}_{R\nu^{\prime}}(q^{\prime})=\delta_{\nu,\nu^{\prime}}\delta_{q,q^{\prime}},$
$\displaystyle\tilde{\underline{\psi}}^{*}_{L\nu}(q)\cdot\tilde{\underline{\psi}}_{R0}=\tilde{\underline{\psi}}^{*}_{L0}\cdot\tilde{\underline{\psi}}_{R\nu}(q)=0.$
(A19)

## References

  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-Hermitian systems, Rev. Mod. Phys. 93, 015005 (2021).
  * Gong  _et al._ [2018] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Topological phases of non-Hermitian systems, Phys. Rev. X 8, 031079 (2018).
  * Shen  _et al._ [2018] H. Shen, B. Zhen, and L. Fu, Topological band theory for non-Hermitian Hamiltonians, Phys. Rev. Lett. 120, 146402 (2018).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Kawabata  _et al._ [2019a] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-Hermitian physics, Phys. Rev. X 9, 041015 (2019a).
  * Zhou and Lee [2019] H. Zhou and J. Y. Lee, Periodic table for topological bands with non-Hermitian symmetries, Phys. Rev. B 99, 235112 (2019).
  * Budich  _et al._ [2019] J. C. Budich, J. Carlström, F. K. Kunst, and E. J. Bergholtz, Symmetry-protected nodal phases in non-Hermitian systems, Phys. Rev. B 99, 041406(R) (2019).
  * Yoshida  _et al._ [2019] T. Yoshida, R. Peters, N. Kawakami, and Y. Hatsugai, Symmetry-protected exceptional rings in two-dimensional correlated systems with chiral symmetry, Phys. Rev. B 99, 121101(R) (2019).
  * Kawabata  _et al._ [2019b] K. Kawabata, T. Bessho, and M. Sato, Classification of exceptional points and non-Hermitian topological semimetals, Phys. Rev. Lett. 123, 066405 (2019b).
  * Delplace  _et al._ [2021] P. Delplace, T. Yoshida, and Y. Hatsugai, Symmetry-protected multifold exceptional points and their topological characterization, Phys. Rev. Lett. 127, 186602 (2021).
  * Mandal and Bergholtz [2021] I. Mandal and E. J. Bergholtz, Symmetry and higher-order exceptional points, Phys. Rev. Lett. 127, 186601 (2021).
  * Stålhammar and Bergholtz [2021] M. Stålhammar and E. J. Bergholtz, Classification of exceptional nodal topologies protected by $\mathcal{PT}$ symmetry, Phys. Rev. B 104, L201104 (2021).
  * Lee [2016] T. E. Lee, Anomalous edge state in a non-Hermitian lattice, Phys. Rev. Lett. 116, 133903 (2016).
  * Martinez Alvarez  _et al._ [2018] V. M. Martinez Alvarez, J. E. Barrios Vargas, and L. E. F. Foa Torres, Non-Hermitian robust edge states in one dimension: Anomalous localization and eigenspace condensation at exceptional points, Phys. Rev. B 97, 121401(R) (2018).
  * Xiong [2018] Y. Xiong, Why does bulk boundary correspondence fail in some non-Hermitian topological models, J. Phys. Commun. 2, 035043 (2018).
  * Kunst  _et al._ [2018] F. K. Kunst, E. Edvardsson, J. C. Budich, and E. J. Bergholtz, Biorthogonal bulk-boundary correspondence in non-Hermitian systems, Phys. Rev. Lett. 121, 026808 (2018).
  * Yao and Wang [2018] S. Yao and Z. Wang, Edge states and topological invariants of non-Hermitian systems, Phys. Rev. Lett. 121, 086803 (2018).
  * Kunst and Dwivedi [2019] F. K. Kunst and V. Dwivedi, Non-Hermitian systems and topology: A transfer-matrix perspective, Phys. Rev. B 99, 245116 (2019).
  * Herviou  _et al._ [2019] L. Herviou, J. H. Bardarson, and N. Regnault, Defining a bulk-edge correspondence for non-Hermitian Hamiltonians via singular-value decomposition, Phys. Rev. A 99, 052118 (2019).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Topological origin of non-Hermitian skin effects, Phys. Rev. Lett. 124, 086801 (2020).
  * Edvardsson  _et al._ [2019] E. Edvardsson, F. K. Kunst, and E. J. Bergholtz, Non-Hermitian extensions of higher-order topological phases and their biorthogonal bulk-boundary correspondence, Phys. Rev. B 99, 081302(R) (2019).
  * Longhi [2019] S. Longhi, Probing non-Hermitian skin effect and non-Bloch phase transitions, Phys. Rev. Research 1, 023013 (2019).
  * Leykam  _et al._ [2017] D. Leykam, K. Y. Bliokh, C. Huang, Y. D. Chong, and F. Nori, Edge modes, degeneracies, and topological numbers in non-Hermitian systems, Phys. Rev. Lett. 118, 040401 (2017).
  * Lee and Thomale [2019] C. H. Lee and R. Thomale, Anatomy of skin modes and topology in non-Hermitian systems, Phys. Rev. B 99, 201103(R) (2019).
  * Lee  _et al._ [2019] C. H. Lee, L. Li, and J. Gong, Hybrid higher-order skin-topological modes in nonreciprocal systems, Phys. Rev. Lett. 123, 016805 (2019).
  * Liu  _et al._ [2019] T. Liu, Y.-R. Zhang, Q. Ai, Z. Gong, K. Kawabata, M. Ueda, and F. Nori, Second-order topological phases in non-Hermitian systems, Phys. Rev. Lett. 122, 076801 (2019).
  * [29] C. Fleckenstein, A. Zorzato, D. Varjas, E. J. Bergholtz, J. H. Bardarson, and A. Tiwari, Non-Hermitian topology in monitored quantum circuits, arXiv:2201.05341 .
  * Zirnstein  _et al._ [2021] H. G. Zirnstein, G. Refael, and B. Rosenow, Bulk-boundary correspondence for non-Hermitian Hamiltonians via Green functions, Phys. Rev. Lett. 126, 216407 (2021).
  * Borgnia  _et al._ [2020] D. S. Borgnia, A. J. Kruchkov, and R.-J. Slager, Non-Hermitian boundary modes and topology, Phys. Rev. Lett. 124, 056802 (2020).
  * Yokomizo and Murakami [2019] K. Yokomizo and S. Murakami, Non-Bloch band theory of non-Hermitian systems, Phys. Rev. Lett. 123, 066404 (2019).
  * Yang  _et al._ [2021] K. Yang, S. C. Morampudi, and E. J. Bergholtz, Exceptional spin liquids from couplings to the environment, Phys. Rev. Lett. 126, 077201 (2021).
  * Edvardsson  _et al._ [2020] E. Edvardsson, F. K. Kunst, T. Yoshida, and E. J. Bergholtz, Phase transitions and generalized biorthogonal polarization in non-Hermitian systems, Phys. Rev. Research 2, 043046 (2020).
  * Koch and Budich [2020] R. Koch and J. C. Budich, Bulk-boundary correspondence in non-Hermitian systems: Stability analysis for generalized boundary conditions, Eur. Phys. J. D 74, 70 (2020).
  * Schomerus [2020] H. Schomerus, Nonreciprocal response theory of non-Hermitian mechanical metamaterials: Response phase transition from the skin effect of zero modes, Phys. Rev. Research 2, 013058 (2020).
  * Brzezicki and Hyart [2019] W. Brzezicki and T. Hyart, Hidden Chern number in one-dimensional non-Hermitian chiral-symmetric systems, Phys. Rev. B 100, 161105(R) (2019).
  * Budich and Bergholtz [2020] J. C. Budich and E. J. Bergholtz, Non-Hermitian topological sensors, Phys. Rev. Lett. 125, 180403 (2020).
  * Zhang  _et al._ [2020] K. Zhang, Z. Yang, and C. Fang, Correspondence between winding numbers and skin modes in non-Hermitian systems, Phys. Rev. Lett. 125, 126402 (2020).
  * Brandenbourger  _et al._ [2019] M. Brandenbourger, X. Locsin, E. Lerner, and C. Coulais, Non-reciprocal robotic metamaterials, Nat. Commun. 10, 4608 (2019).
  * Ghatak  _et al._ [2020] A. Ghatak, M. Brandenbourger, J. van Wezel, and C. Coulais, Observation of non-Hermitian topology and its bulk–edge correspondence in an active mechanical metamaterial, PNAS 117, 29561 (2020).
  * Helbig  _et al._ [2020] T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. Molenkamp, C. Lee, A. Szameit, M. Greiter, and R. Thomale, Generalized bulk–boundary correspondence in non-Hermitian topolectrical circuits, Nat. Phys. 16, 747 (2020).
  * Hofmann  _et al._ [2020] T. Hofmann, T. Helbig, F. Schindler, N. Salgo, M. Brzezińska, M. Greiter, T. Kiessling, D. Wolf, A. Vollhardt, A. Kabaši, C. H. Lee, A. Bilušić, R. Thomale, and T. Neupert, Reciprocal skin effect and its realization in a topolectrical circuit, Phys. Rev. Research 2, 023265 (2020).
  * Xiao  _et al._ [2020] L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Non-Hermitian bulk–boundary correspondence in quantum dynamics, Nat. Phys. 16, 761 (2020).
  * Weidemann  _et al._ [2020] S. Weidemann, M. Kremer, T. Helbig, T. Hofmann, A. Stegmaier, M. Greiter, R. Thomale, and A. Szameit, Topological funneling of light, Science 368, 311 (2020).
  * Koch and Budich [2022] F. Koch and J. C. Budich, Quantum non-Hermitian topological sensors, Phys. Rev. Research 4, 013113 (2022).
  * McDonald and Clerk [2020] A. McDonald and A. A. Clerk, Exponentially-enhanced quantum sensing with non-Hermitian lattice dynamics, Nat. Commun. 11, 5382 (2020).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * Langen  _et al._ [2015] T. Langen, R. Geiger, and J. Schmiedmayer, Ultracold atoms out of equilibrium, Annu. Rev. Condens. Matter Phys. 6, 201 (2015).
  * Diehl  _et al._ [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. Büchler, and P. Zoller, Quantum states and phases in driven open quantum systems with cold atoms, Nat. Phys. 4, 878 (2008).
  * 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).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantum computation and quantum-state engineering driven by dissipation, Nat. Phys. 5, 633 (2009).
  * Krauter  _et al._ [2011] H. Krauter, C. A. Muschik, K. Jensen, W. Wasilewski, J. M. Petersen, J. I. Cirac, and E. S. Polzik, Entanglement generated by dissipation and steady state entanglement of two macroscopic objects, Phys. Rev. Lett. 107, 080503 (2011).
  * Breuer and Petruccione [2007] H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Oxford, 2007).
  * Diehl  _et al._ [2011] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Topology by dissipation in atomic quantum wires, Nat. Phys. 7, 971 (2011).
  * Bardyn  _et al._ [2013] C. E. Bardyn, M. A. Baranov, C. V. Kraus, E. Rico, A. İmamoglu, P. Zoller, and S. Diehl, Topology by dissipation, New J. Phys. 15, 085001 (2013).
  * Budich  _et al._ [2015] J. C. Budich, P. Zoller, and S. Diehl, Dissipative preparation of Chern insulators, Phys. Rev. A 91, 042117 (2015).
  * Liu  _et al._ [2021] Z. Liu, E. J. Bergholtz, and J. C. Budich, Dissipative preparation of fractional Chern insulators, Phys. Rev. Research 3, 043119 (2021).
  * He  _et al._ [2022] P. He, Y.-G. Liu, J.-T. Wang, and S.-L. Zhu, Damping transition in an open generalized Aubry-André-Harper model, Phys. Rev. A 105, 023311 (2022).
  * Song  _et al._ [2019] F. Song, S. Yao, and Z. Wang, Non-Hermitian skin effect and chiral damping in open quantum systems, Phys. Rev. Lett. 123, 170401 (2019).
  * Wanjura  _et al._ [2020] C. C. Wanjura, M. Brunelli, and A. Nunnenkamp, Topological framework for directional amplification in driven-dissipative cavity arrays, Nat. Commun. 11, 3149 (2020).
  * Liu  _et al._ [2020] C.-H. Liu, K. Zhang, Z. Yang, and S. Chen, Helical damping and dynamical critical skin effect in open quantum systems, Phys. Rev. Research 2, 043167 (2020).
  * Mao  _et al._ [2021] L. Mao, T. Deng, and P. Zhang, Boundary condition independence of non-Hermitian Hamiltonian dynamics, Phys. Rev. B 104, 125435 (2021).
  * Wanjura  _et al._ [2021] C. C. Wanjura, M. Brunelli, and A. Nunnenkamp, Correspondence between non-Hermitian topology and directional amplification in the presence of disorder, Phys. Rev. Lett. 127, 213601 (2021).
  * McDonald  _et al._ [2022] A. McDonald, R. Hanai, and A. A. Clerk, Nonequilibrium stationary states of quantum non-Hermitian lattice models, Phys. Rev. B 105, 064302 (2022).
  * Haga  _et al._ [2021] T. Haga, M. Nakagawa, R. Hamazaki, and M. Ueda, Liouvillian skin effect: Slowing down of relaxation processes without gap closing, Phys. Rev. Lett. 127, 070402 (2021).
  * [67] Z. Zhou and Z. Yu, Skin Effect in Quadratic Lindbladian Systems: An Adjoint Fermion Approach, arXiv:2110.09874 .
  * Prosen [2008] T. Prosen, Third quantization: A general method to solve master equations for quadratic open Fermi systems, New J. Phys. 10, 043026 (2008).
  * Prosen and Žunkovič [2010] T. Prosen and B. Žunkovič, Exact solution of Markovian master equations for quadratic Fermi systems: Thermal baths, open XY spin chains and non-equilibrium phase transition, New J. Phys. 12, 025016 (2010).
  * Prosen [2010] T. Prosen, Spectral theorem for the Lindblad equation for quadratic open fermionic systems, J. Stat. Mech. 2010, P07020 (2010).
  * Prosen and Ilievski [2011] T. Prosen and E. Ilievski, Nonequilibrium phase transition in a periodically driven XY spin chain, Phys. Rev. Lett. 107, 060403 (2011).
  * van Caspel  _et al._ [2019] M. van Caspel, S. E. Tapias Arze, and I. Pérez Castillo, Dynamical signatures of topological order in the driven-dissipative Kitaev chain, SciPost Phys. 6, 026 (2019).
  * Lieu  _et al._ [2020] S. Lieu, M. McGinley, and N. R. Cooper, Tenfold way for quadratic Lindbladians, Phys. Rev. Lett. 124, 040401 (2020).
  * Minganti  _et al._ [2018] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Spectral theory of Liouvillians for dissipative phase transitions, Phys. Rev. A 98, 042118 (2018).
  * Brody [2014] D. C. Brody, Biorthogonal quantum mechanics, J. Phys. A 47, 035305 (2014).
  * Hatano and Nelson [1996] N. Hatano and D. R. Nelson, Localization transitions in non-Hermitian quantum mechanics, Phys. Rev. Lett. 77, 570 (1996).
  * Hatano and Nelson [1997] N. Hatano and D. R. Nelson, Vortex pinning and non-Hermitian quantum mechanics, Phys. Rev. B 56, 8651 (1997).
  * Hatano and Nelson [1998] N. Hatano and D. R. Nelson, Non-Hermitian delocalization and eigenfunctions, Phys. Rev. B 58, 8384 (1998).
  * Altland  _et al._ [2021] A. Altland, M. Fleischhauer, and S. Diehl, Symmetry classes of open fermionic quantum matter, Phys. Rev. X 11, 021037 (2021).
  * Hatano [2019] N. Hatano, Exceptional points of the Lindblad operator of a two-level system, Mol. Phys. 117, 2121 (2019).
  * Benatti  _et al._ [2021] F. Benatti, R. Floreanini, and L. Memarzadeh, Exact steady state of the open XX-spin chain: Entanglement and transport properties, PRX Quantum 2, 030344 (2021).
  * Sayyad  _et al._ [2021] S. Sayyad, J. Yu, A. G. Grushin, and L. M. Sieberer, Entanglement spectrum crossings reveal non-Hermitian dynamical topology, Phys. Rev. Research 3, 033022 (2021).
  * Scarlatella  _et al._ [2021] O. Scarlatella, A. A. Clerk, R. Fazio, and M. Schiró, Dynamical mean-field theory for Markovian open quantum many-body systems, Phys. Rev. X 11, 031018 (2021).
  * Yamamoto  _et al._ [2022] K. Yamamoto, M. Nakagawa, M. Tezuka, M. Ueda, and N. Kawakami, Universal properties of dissipative Tomonaga-Luttinger liquids: Case study of a non-Hermitian XXZ spin chain, Phys. Rev. B 105, 205125 (2022).
  * McDonald and Clerk [2022] A. McDonald and A. A. Clerk, Exact solutions of interacting dissipative systems via weak symmetries, Phys. Rev. Lett. 128, 033602 (2022).
  * [86] L. Sá, P. Ribeiro, and T. Prosen, Lindbladian Dissipation of Strongly-Correlated Quantum Matter, arXiv:2112.12109 .
  * Kunst  _et al._ [2019] F. K. Kunst, G. van Miert, and E. J. Bergholtz, Extended Bloch theorem for topological lattice models with open boundaries, Phys. Rev. B 99, 085427 (2019).