Non-Unitary Quantum Many-Body Dynamics using the Faber Polynomial Method

Rafael D. Soares1,2*, M. Schirò2

1 Université Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France

2 JEIP, UAR 3573 CNRS, Collège de France, PSL Research University, 11 Place
Marcelin Berthelot, 75321 Paris Cedex 05, France

<EMAIL_ADDRESS>
## Abstract

Efficient numerical methods are still lacking to probe the unconventional
dynamics of quantum many-body systems under non-unitary evolution. In this
work, we use Faber polynomials to numerically simulate both the dynamics of
non-Hermitian systems and the quantum jumps unravelling of the Lindblad
dynamics. We apply the method to the non-interacting and interacting Hatano-
Nelson models evolving from two different setups: i) a Néel state, and ii) a
domain wall. In the first case, we study how interactions preserve the initial
magnetic order against the skin effect. In the second example, we present
numerical evidence of the existence of an effective hydrodynamic description
for the domain-wall melting problem in the non-interacting limit.
Additionally, we investigate both the conditional and unconditional dynamics
of the quantum jump unravelling in two quantum spin chains, which exhibit
either the non-Hermitian or the Liouvillian skin effect. This numerical method
inherently generalises the well-established method based on Chebyshev
polynomials to accommodate non-Hermitian scenarios.

###### Contents

  1. 1 Introduction
  2. 2 Non-Unitary Dynamics of Open Quantum Systems
  3. 3 Faber Polynomial Method
    1. 3.1 Warm-Up: Unitary Evolution
    2. 3.2 Non-Unitary Evolution
      1. 3.2.1 Convergence
  4. 4 Application to Non-Hermitian Gaussian Systems
    1. 4.1 Benchmark: Dynamics of the Hatano-Nelson model
    2. 4.2 Domain Wall Melting for Hatano-Nelson
  5. 5 Application to Non-Hermitian Many-Body Systems
    1. 5.1 Magnetisation Dynamics in the Interacting Hatano-Nelson Model
    2. 5.2 Effect of interaction in the Domain-Wall melting for Hatano-Nelson
  6. 6 Application to Quantum Jumps Unravelling
    1. 6.1 Magnetisation and Entanglement Dynamics in Monitored Spin Chains
  7. 7 Conclusions
  8. A Further Details on Faber Polynomials
  9. B General Features on the Hatano-Nelson Model
  10. C Comparison with MPS based Methods

##  1 Introduction

In recent years, the scientific community has shown a growing interest in
elucidating the distinctive characteristics of many-body quantum systems
subjected to effective non-unitary dynamics. In quantum mechanics, non-unitary
dynamics typically arises when a closed quantum system interacts with an
external environment, leading to dissipation, decoherence, or wave function
collapse that disrupts the usual unitary Schrödinger evolution. While a fully
microscopic description of both environment and system is a daunting task, in
many cases of experimental and theoretical interest one can assume the
dynamics of the environment to be sufficiently fast, allowing for the
derivation of a local-in-time, Markovian, non-unitary evolution for the system
of interest [1, 2].

Different types of Markovian open quantum system dynamics have been considered
in the literature. A first relevant example is provided by systems that evolve
according to the Lindblad master equation [3, 4, 5]. Here, the evolution of
the system density matrix is generated by a non-unitary (super)operator, the
Lindbladian. Many-body versions of Lindblad master equations have been studied
in a number of contexts and with different objectives, from dissipative phase
transitions [6] to quantum transport [7]. Another class of non-unitary
dynamics arises for continuously monitored quantum systems [2, 8, 9], whose
stochastic evolution - a so-called quantum trajectory [10] \- is described by
a non-unitary unravelling of the Lindblad master equation [11]. In particular,
under the quantum jump unravelling, a deterministic non-unitary evolution is
driven by a non-Hermitian Hamiltonian between stochastic measurements [12, 13,
14]. Recent works have raised interest in possible phase transitions in the
entanglement structure of those quantum many-body trajectories [15, 16, 17,
18, 19, 20, 21].

Finally, as a last example of non-unitary Markovian dynamics, we can consider
the one generated by a purely non-Hermitian Hamiltonian. Non-Hermitian physics
emerges intrinsically in various domains of physics beyond quantum physics,
encompassing photonics [22], hydrodynamics [23] and active matter [24, 25]. In
the context of open quantum systems, a non-Hermitian evolution can be obtained
by post-selecting the quantum trajectories corresponding to no quantum jumps,
i.e. the no-click limit [26]. Non-Hermitian quantum systems show anomalous
static and dynamical properties, which are attracting widespread interest.
Among these, we mention the unconventional propagation of quantum correlations
[27, 28, 29], the distinct entanglement transitions generated by time
evolution [30, 31, 32, 33, 34] or their extraordinary sensitivity to boundary
conditions, also known as the skin effect [35, 36, 37, 38, 39, 40]. The latter
manifests itself through the unusual localisation of all single-particle
eigenstates at the system’s edges under open boundary conditions [41, 42, 43,
44, 45]. Furthermore, it has unique signatures in the dynamics, leading to
non-reciprocal transport[46].

Besides these theoretical developments, unlike their Hermitian counterparts,
the toolbox of computational many-body physics is more limited when it comes
to non-unitary dynamics. In this work, we introduce a new method to tackle the
quantum dynamics of a non-unitary system. Our approach is based on expanding
the evolution operator in Faber polynomials [47, 48]. This numerical approach
is a natural generalisation of its Hermitian counterpart, the Chebyshev
polynomial method for time evolution [49], which has been the primary choice
for efficiently simulating nonequilibrium transport phenomena in both
interacting [50] and non-interacting systems [51, 52, 53]. Although time
evolution integrators based on Faber polynomials have already been proposed in
some works [54, 55, 56], for example, in the simulation of electromagnetic
wave propagation through passive media [57, 58], it appears that their full
potential remains largely unexplored. Namely, in the simulation of non-
Hermitian quantum systems, they could have a significant impact owing to their
numerical stability and adjustable accuracy compared to other methods such as
integrators based on Runge-Kutta [59] or Trotterization techniques [60].

In this manuscript, we test, benchmark and apply the Faber polynomial method
to investigate the time evolution of particle density, charge current, and
entanglement in several setups involving the interacting and non-interacting
Hatano-Nelson model [61, 62], a paradigmatic non-Hermitian model showing non-
reciprocity [63] and the skin effect at the single particle level, and non-
Hermitian quantum spin chains. In addition, we merge the Faber polynomial
method with quantum jumps in order to simulate the dynamics of the full
Lindblad master equation through a suitable unravelling and the conditional
dynamics encoded in the entanglement of quantum trajectories.

The manuscript is structured as follows. In Sec. 2 we set the stage and define
the classes of non-unitary dynamics that we will focus on throughout this
manuscript. In Sec. 3 we describe the Faber polynomial method for non-unitary
dynamics and discuss its convergence. In Sec. 4 we present our first
application to quadratic (Gaussian) non-Hermitian systems. In particular, we
explore the melting of a domain-wall state under the Hatano-Nelson
Hamiltonian. Furthermore, in Sec. 5, we focus on the spin version of the many-
body interacting Hatano-Nelson chain, examining the evolution of both an
initial Néel state and a domain wall under the influence of non-reciprocal
hopping and interactions. Finally, in Sec. 6, we apply the method to the
stochastic quantum jump dynamics obtained by the unravelling of a Lindblad
master equation. Sec. 7 summarises our conclusions and discusses potential
future research directions.

##  2 Non-Unitary Dynamics of Open Quantum Systems

In this work, we focus on the dynamics of open quantum many-body systems
described by a Hamiltonian $\mathcal{H}$ and a set of independent
environments. A typical example will be a quantum spin chain connected on each
site to an external bath. In practice, we will always assume a Markovian
description of the environment, which can also be identified as a measurement
apparatus that monitors a certain physical property of the system, for example
the particle density [17, 19, 20, 21]. However, our primary focus is on the
non-unitary dynamics of the system, obtained by tracing out the environment.
This is naturally modelled by the Lindblad master equation and its unravelling
[9]. In this setting, we consider two types of quantum dynamics: (i) the
stochastic evolution of the system conditioned to a given set of measurement
outcomes, and (ii) the dynamics of the averaged state. In the first case, the
system evolves according to the stochastic Schrödinger equation,

$\displaystyle d\ket{\psi(\xi_{t},t)}$
$\displaystyle=-idt\left[\mathcal{H}-\frac{i}{2}\sum_{\mu}\left(L^{\dagger}_{\mu}L_{\mu}-\langle
L^{\dagger}_{\mu}L_{\mu}\rangle_{t}\right)\right]\ket{\psi(\xi_{t},t)}$ (1)
$\displaystyle\quad+\sum_{\mu}\left(\frac{L_{\mu}}{\sqrt{\langle
L^{\dagger}_{\mu}L_{\mu}\rangle}}-1\right)d\xi_{\mu,t}\ket{\psi(\xi_{t},t)},$

where
$\langle\circ\rangle_{t}\equiv\langle\psi(\xi_{t},t)|\circ|\psi(\xi_{t},t)\rangle$,
and $\xi_{t}=\left\\{\xi_{\mu,t}\right\\}$ are a set of statistically
independent Poisson processes ${d\xi_{\mu,t}\in\\{0,1\\}}$ with average value
$\overline{d\xi_{\mu,t}}=dt\langle L^{\dagger}_{\mu}L_{\mu}\rangle_{t}$. The
above dynamics breaks down into two steps: a deterministic non-unitary
evolution driven by a non-Hermitian Hamiltonian

$\displaystyle\mathcal{H}_{\rm
nH}=\mathcal{H}-\frac{i}{2}\sum_{\mu}L^{\dagger}_{\mu}L_{\mu},$ (2)

and a series of stochastic quantum jumps at random times, at which the wave
function changes discontinuously (see second line of Eq. 1). We note that the
non-Hermitian evolution is normalised and state dependent. This is encoded in
the last term in the first line of Eq. 1. If one post-selects the quantum
trajectories over the records of no click, the dynamics is deterministic and
driven by $\mathcal{H}_{\rm nH}$ [33, 32, 31]. Otherwise, if one considers all
the trajectories and averages over the measurement outcomes, the conditional
density matrix
$\rho_{c}(\xi_{t},t)=|\psi(\xi_{t},t)\rangle\langle\psi(\xi_{t},t)|$, i.e.

$\rho(t)=\overline{\rho_{c}(\xi_{t},t)},$ (3)

evolves according to the Lindblad master equation with jump operators
$L_{\mu}$, i.e.

$d\rho(t)=-idt\left[\mathcal{H},\rho\right]+dt\sum_{\mu}\left(L_{\mu}\rho(t)L^{\dagger}_{\mu}-\frac{1}{2}\left\\{L^{\dagger}_{\mu}L_{\mu},\rho\right\\}\right).$
(4)

In both cases the basic building block of the non-unitary dynamics is the
evolution driven by a non-Hermitian Hamiltonian. In the next section, we
introduce the Faber polynomial method to accurately solve the time evolution
governed by a non-unitary Schrödinger equation,

##  3 Faber Polynomial Method

The knowledge of the time evolution operator, $\mathcal{U}(t)$, allows a
comprehensive description of the physical properties of a system when it is
far from equilibrium. This operator is necessary to propagate a given initial
state, $\ket{\Psi(t)}=\mathcal{U}(t)\ket{\Psi_{0}}$, allowing the calculation
of observables that characterise the nonequilibrium state. In principle, an
exact expression for the state is necessary as soon as one moves beyond the
scope of linear response theory. Nevertheless, this problem is equivalent to
solving for the spectrum and eigenstates of the Hamiltonian, as in the case of
a time-independent Hamiltonian,

$\mathcal{U}(t)=\exp\left(-i\mathcal{H}t\right).$ (5)

The idea behind both the Chebyshev (unitary evolution) and Faber (non-unitary)
polynomial methods is to perform an expansion of the time evolution in the
respective polynomial basis,

$\mathcal{U}(t)=\sum_{n=0}^{+\infty}c_{n}\left(t\right)\mathcal{P}_{n}\left(\mathcal{H}\right),$
(6)

where $c_{n}\left(t\right)$ is the $n^{\rm th}$ coefficient of the series
expansion and $\mathcal{P}_{n}$ is the $n^{\rm th}$ polynomial, which
corresponds to a Chebyshev polynomial of the first kind or to a Faber
polynomial, depending on the situation. Then the state after the time step,
$\delta t$, can be approximated by truncating the series expansion to the
order $N_{p}$,

$\ket{\Psi(t_{0}+\delta t)}\simeq\sum_{n=0}^{N_{p}-1}c_{n}(\delta
t)\ket{\Psi_{n}},$ (7)

where we define
$\ket{\Psi_{n}}=\mathcal{P}_{n}\left(\mathcal{H}\right)\ket{\Psi\left(t_{0}\right)}$.
As will be demonstrated, the coefficients $c_{n}(\delta t)$ decrease as the
order $n$ increases. Moreover, the states $\ket{\Psi_{n}}$ are efficiently
computed through the recurrence relations that the polynomials satisfy. To
compute the subsequent level of the expansion, the main computational task is
the application of the system’s Hamiltonian onto a particular state.
Consequently, the most demanding operation involves only multiplying the
Hamiltonian by a limited group of vectors, leading to a resource usage that
increases linearly with $\dim\left(\mathcal{H}\right)$ for sparse matrices or
quadratically with $\dim\left(\mathcal{H}\right)$ for dense matrices. Linear
scaling is expected for Hamiltonians characterising a system with short-range
interactions or hoppings. These principles are exactly those underpinning
Kernel Polynomial Methods [64], which has become an essential computational
resource in condensed matter physics, particularly for calculating various
spectral quantities[65, 66, 67, 68, 69].

###  3.1 Warm-Up: Unitary Evolution

When the time evolution is generated by a Hermitian Hamiltonian, one can
expand Eq. 5 using Chebyshev polynomials of the first kind (for further
details check Tal-Ezer and Kosloff [49]),

$\mathcal{U}\left(t\right)=\sum_{n=0}^{\infty}c_{n}(t)T_{n}\left(\tilde{\mathcal{H}}\right),\quad
c_{n}(t)=\frac{2}{1+\delta_{n,0}}(-i)^{n}J_{n}(\lambda t),$ (8)

where $\tilde{\mathcal{H}}=\mathcal{H}/\lambda$ is the rescaled Hamiltonian,
$J_{n}(x)$ is the $n^{\text{th}}$ Bessel Function of the first kind, and
$T_{n}$ is the $n^{\text{th}}$ Chebyshev polynomial of the first kind. It is
necessary to rescale the Hamiltonian so that its eigenvalues fall within the
domain of the definition of polynomials, the open interval
$\left(-1,1\right)$. It is always possible to do this given that the
Hamiltonian is always bounded in finite-size lattice models. Using the
recursion relation of the Chebyshev polynomials, the states $\ket{\Psi_{n}}$
in Eq. 7 are computed on the run using,

$\displaystyle\ket{\Psi_{0}}$ $\displaystyle=\ket{\Psi(t_{0})},$ (9)
$\displaystyle\ket{\Psi_{1}}$
$\displaystyle=\tilde{\mathcal{H}}\ket{\Psi_{0}},$
$\displaystyle\ket{\Psi_{n+1}}$
$\displaystyle=2\tilde{\mathcal{H}}\ket{\Psi_{n}}-\ket{\Psi_{n-1}},\quad n\geq
2,$

with $|\Psi_{n}\rangle=T_{n}(\tilde{\mathcal{H}})|\Psi\rangle$.

This expansion is feasible solely because the Hamiltonian is Hermitian.
Therefore, in the context of open quantum systems where a non-Hermitian
operator dictates the dynamics, a different set of polynomials is necessary.
Under these circumstances, the spectrum is defined in the complex plane and
the operator has right and left eigenvectors, which can be distinct.

###  3.2 Non-Unitary Evolution

For a non-Hermitian Hamiltonian, the propagator Eq. 5 is expanded using Faber
polynomials [70, 71] instead. These are a familiar tool in complex analysis,
used as a polynomial basis to represent a complex-valued function within the
domain $\mathcal{D}$ in which it is analytical. The Faber polynomials are
generated by conformal mapping, which maps the complement of a closed disk of
radius $\rho$ to the complement of the region containing all the spectra of
the Hamiltonian (the domain $\mathcal{D}$). In our work, we assume
$\mathcal{D}$ to be an elliptic region containing all the eigenvalues of
$\mathcal{H}$. This fits our purposes, as the conformal mapping associated
with this shape generates a class of Faber polynomials with a minimum
recurrence relation (see the Appendix A for further details). The expansion of
Eq. 5 for a non-Hermitian Hamiltonian with this choice reduces to

$\mathcal{U}\left(t\right)=\sum_{n=0}^{+\infty}c_{n}\left(t\right)F_{n}\left(\tilde{\mathcal{H}}\right),\quad
c_{n}(t)=e^{-i\lambda
t\gamma_{0}}\left(\frac{-i}{\sqrt{\gamma_{1}}}\right)^{n}\;J_{n}\left(2\sqrt{\gamma_{1}}\lambda
t\right).$ (10)

where $F_{n}$ is the $\rm n^{\text{th}}$ Faber polynomial. The parameter
$\lambda$ is used to rescale the Hamiltonian so that the norm of
$F_{n}\left(z\right)$ is bounded [57], and it is obtained from the bounds of
the real and imaginary part of the spectra. $\gamma_{0}$ and, $\gamma_{1}$ are
associated with the details of the elliptic contour chosen, $\gamma_{0}$ is
the centre of the ellipse and $\gamma_{1}=1-b$, where $b$ is the semi major-
axis111Recall that the equation of an ellipse is:
$\frac{(x-x_{0})^{2}}{a^{2}}+\frac{(y-y_{0})^{2}}{b^{2}}=1$.. The ellipse must
be constructed to be close as possible to eigenvalues of $\mathcal{H}$,
thereby reducing the magnitude of $\lambda$. One can show [58] that this is
achieved using

$\begin{aligned} \lambda&=\frac{\left(\ell^{2/3}+p^{2/3}\right)^{3/2}}{2},\\\
\gamma_{1}&=\frac{\left(\tilde{p}^{2/3}+\tilde{\ell}^{2/3}\right)\left(\tilde{p}^{4/3}-\tilde{\ell}^{4/3}\right)}{4\lambda}\end{aligned},$
(11)

with
$\ell=\left[\max\left[\text{Im}\left(E\right)\right]-\min\left[\text{Im}\left(E\right)\right]\right]/2$,
$p=\left[\max\left[\text{Re}\left(E\right)\right]-\min\left[\text{Re}\left(E\right)\right]\right]/2$
, $\tilde{\ell}=\ell/\lambda$ and $\tilde{p}=p/\lambda$. Using the recursion
relation of the Faber Polynomials (consult Appendix A) the states
$\ket{\Psi_{n}}=F_{n}\left(\tilde{\mathcal{H}}\right)\ket{\Psi}$ are given by,

$\displaystyle\ket{\Psi_{0}}$ $\displaystyle=\ket{\Psi\left(t_{0}\right)}$
(12) $\displaystyle\ket{\Psi_{1}}$
$\displaystyle=\left(\tilde{\mathcal{H}}-\gamma_{0}\right)\ket{\Psi_{0}}$
$\displaystyle\ket{\Psi_{2}}$
$\displaystyle=\left(\tilde{\mathcal{H}}-\gamma_{0}\right)\ket{\Psi_{1}}-2\gamma_{1}\ket{\Psi_{0}}$
$\displaystyle\ket{\Psi_{n+1}}$
$\displaystyle=\left(\tilde{\mathcal{H}}-\gamma_{0}\right)\ket{\Psi_{n}}-\gamma_{1}\ket{\Psi_{n-1}},\quad
n>2$

In order to perform one time step of evolution, one simply has to truncate Eq.
10 up to the desired order $N_{p}$, and calculate the associated states
$\ket{\Psi_{n}}$. Through the relations in Eq. 12, one never has to store the
Hamiltonian matrix, only needing to store in memory at most the previous two
states, $\ket{\Psi_{n-1}}$ and $\ket{\Psi_{n}}$, to compute the following term
of the expansion $\ket{\Psi_{n+1}}$. Furthermore, the expansion coefficients
can be computed once at the beginning of the algorithm, as they depend only on
the chosen time step. Given this, the computation of the state of the system
after a time step scales linearly in the number of polynomials and linearly in
the Hilbert space dimension (assuming that the Hamiltonian has a sparse
representation on the used basis). This scaling can be improved by using
parallelisation techniques [68] and making use of the underlying symmetries of
the Hamiltonian [72].

An additional procedure is required when addressing purely non-Hermitian
dynamics, specifically, ensuring the normalisation of the quantum state
throughout the time evolution. Theoretically, this normalisation can only be
done prior to the computation of an observable. However, the algorithm may
exhibit instability if the coefficients fluctuate within the limits of machine
precision. Consequently, it is prudent to normalise the state following each
time step,

$\ket{\Psi\left(t+\delta t\right)}=\dfrac{\mathcal{U}\left(\delta
t\right)\ket{\Psi\left(t\right)}}{\left\|\mathcal{U}\left(\delta
t\right)\ket{\Psi\left(t\right)}\right\|}$ (13)

In the following section, we revise how to do this when dealing with fermionic
Gaussian states.

####  3.2.1 Convergence

In this section, we discuss the convergence properties of our algorithm,
illustrating in Fig. 1 how the absolute value associated with the $n^{th}$
Faber polynomial varies for different rescaled time steps, $\lambda\delta t$
noting that a greater time step naturally requires more polynomials to achieve
convergence. For numerical purposes, the Faber series is an exact
representation of the time-evolution operator if the coefficient of the last
polynomial used is within the machine precision. This argument holds because
of Hamiltonian rescaling, as demonstrated in [73], which guarantees that
$\forall_{m}\max_{z\in G}\left|F_{m}(z)\right|\leq 2$. Using the asymptotic
properties of the Bessel functions, we see that the weight of the coefficient
decreases with $n$ according to [74],

$\left|c_{n}\right|\sim\dfrac{\left(\lambda\delta t\right)^{n}}{n!}$ (14)

whenever $\lambda\delta t\ll n$. In the remaining of this work, the number of
polynomials is chosen such that the last coefficient of the absolute value of
the last coefficient of the expansion is of the order $\left|c_{n}\right|\sim
10^{-16}$.

Figure 1:  Absolute value of the coefficients associated with the Faber
expansion of the time-evolution operator. We represent the absolute of the
$n^{th}$ Faber polynomial for different rescaled time steps.

##  4 Application to Non-Hermitian Gaussian Systems

In this Section, we use the Faber polynomial method to study the dynamics of a
fermionic Gaussian non-Hermitian system. Here, the many-body wave function can
be expressed using a single-particle basis, and the Faber polynomials can be
employed to represent the single-particle propagator. In the case of
Hamiltonians with $U(1)$ symmetry, associated with particle number
conservation, for an initial Gaussian state with a well-defined particle
number $M$, the many-body state can always be represented in the form

$\ket{\Psi\left(t\right)}=\prod_{n=0}^{\text{M}-1}\left[\sum_{\ell=0}^{\text{L}-1}\text{U}_{\ell
n}\left(t\right)c_{\ell}^{\dagger}\right]\ket{\text{vac}},$ (15)

with $M$ the total number of particles and $L$ the total number of sites. The
time evolution is given by the following equation.

$i\frac{d}{dt}\textbf{\text{U}}_{n}=\textbf{\text{h}}\;\textbf{\text{U}}_{n},$
(16)

where $\textbf{\text{U}}_{n}$ is the $n^{th}$ column vector and h is the
single-particle Hamiltonian, a $\rm L\times L$ matrix. This translates the
evolution of the many-body state into the evolution of the $\rm M$ single-
particle states, each represented by a column in the matrix U. This is
expected, as the dynamics and characteristics of a non-interacting system can
be simplified to those of the single-particle Hamiltonian, assuming the
appropriate quantum statistics. For a typical tight-binding Hamiltonian with a
finite number of hopping terms, the complexity of our algorithm for a single
time step of evolution scales as $\mathcal{O}(N_{\text{p}}\cdot\rm L\cdot M)$.

Following the time step, it is essential to restore proper normalisation and
particle statistics. In the case of Gaussian states, this is achieved through
a QR decomposition [17], which guarantees that the U matrix is an isometry,
specifically $\text{U}^{\dagger}\text{U}=\mathcal{I}_{M\times M}$.

$\text{U}(t+\delta t)=\text{QR},$ (17)

with Q a unitary $\rm L\times M$ matrix. So the proper-normalised many-body
state is obtained by assigning $\text{U}(t+\delta t)=Q$. Even tough, thought
in this document we only study particle-number conserving models, this method
is easily also applied to non-particle conserving Hamiltonians using similar
techniques [75, 20]. Typically, this normalisation step is the most
computationally intensive part of the time-evolution process. However, unlike
other methods, this step can be executed less frequently since the time-step
does not need to be small. Thus, this computationally intensive procedure can
be minimised while still ensuring high accuracy in the time integration.

The method of directly evolving the state presents a more cost-effective
alternative to approaches that compute the equations of motion for all two-
point functions. The latter typically employs conventional ordinary
differential equation solvers, such as the fourth-order Runge-Kutta method.
Firstly, it necessitates the selection of a short time-step to ensure accurate
integration, a requirement not imposed by the Faber algorithm. Secondly, this
method involves the evolution of the full correlation matrix, which in
practice involves the evolution of $\rm L\cdot\left(L-1\right)/2$ elements. In
contrast, the evolution of U requires the evolution of $\rm L\times M$
elements with $\rm M\leq L$. Lastly, the Faber algorithm requires the storage
of only two vectors of size $\rm L$ in memory to evolve a given column, while
the integration of the correlation matrix demands the storage of four vectors
of size $\rm L$ to evolve each of the $\rm L$ columns.

###  4.1 Benchmark: Dynamics of the Hatano-Nelson model

The Hatano-Nelson (HN) model [61, 62] is a paradigmatic lattice model for non-
Hermitian phenomena. It corresponds to a chain of spinless fermions with an
asymmetric hopping,

$\mathcal{H}_{\text{HN}}=-\dfrac{1}{2}\sum_{\ell=0}^{L-1}\left(\left(J+\gamma\right)c_{\ell}^{\dagger}c_{\ell+1}+\left(J-\gamma\right)c_{\ell+1}^{\dagger}c_{\ell}\right),$
(18)

where $J$ is a hopping term to the first-neighbour, $\gamma\in\mathbb{R}$
parameterised the left-right imbalance in charge hopping, also called non-
reciprocity, and $c^{\dagger}_{\ell}\left(c_{\ell}\right)$ is the creation
(annihilation) operator which creates (destroys) a fermion on site $\ell$.
Under open boundary conditions (OBC), this non-reciprocal hopping gives rise
to a unique phenomenon of non-Hermitian systems, the skin effect [35, 36, 37,
38, 39, 40]. This corresponds to the localization of the single-particle
eigenstates at the edges of the system. In addition, the model has a huge
sensitivity to the boundary conditions: under periodic boundary conditions
(PBC), the single-particle spectrum encircles $E=0$ in the complex plane, with
the eigenstates manifesting as delocalised plane waves. In contrast, with OBC,
the spectrum is real for $\left|\gamma\right|<J$ and purely imaginary for
$\left|\gamma\right|>J$. Additionally, all right single-particle eigenstates
show an exponential localisation at the left (right) boundary for positive
(negative) $\gamma$ (check Appendix. B for further details).

Figure 2: Left panel: Comparison between the results obtained with the Faber
polynomial method and those reported in [34] for the entanglement entropy of
half of the chain for a total chain of size $L=100$. Right panel: Entanglement
entropy for half of the chain, with $\gamma=-0.05J$ using a different number
of polynomials and the time step of $0.1J^{-1}$. We note that we use a
symmetric definition for $\gamma$ with regard to [34]

In the following, we benchmark our results using the Faber polynomial method,
with those reported by Kawabata et al. [34]. That is, we investigate the
dynamics of the entanglement [76] associated with a segment of the chain,
denoted $\ell$. This is rigorously derived from the von Neumann entropy of the
reduced density matrix($\rho_{\ell}$) [77],

$S_{\ell}(t)=-\text{Tr}\left(\rho_{\ell}\ln\rho_{\ell}\right).$ (19)

$\rho_{\ell}$ is determined by tracing out the complementary degrees of
freedom of the subregion $\ell$. However, the Gaussianity of the state allows
us to use the standard techniques [78, 79, 80] to perform this computation
using the one-particle density matrix restricted to the lattice sites
belonging to the region $\ell$, $\mathcal{C}_{n,m\in\ell}=\left\langle
c^{\dagger}_{n}c_{m}\right\rangle$. Thus, Eq. 19 for a free fermionic system
is simplified to

$S_{\ell}(t)=-\text{Tr}\left(\left.\mathcal{C}\right|_{\ell}\ln\left.\mathcal{C}\right|_{\ell}+\left(\mathcal{I}_{\ell\times\ell}-\left.\mathcal{C}\right|_{\ell}\right)\ln\left(\mathcal{I}_{\ell\times\ell}-\left.\mathcal{C}\right|_{\ell}\right)\right),$
(20)

where $\mathcal{I}_{\ell\times\ell}$ is the $\ell\times\ell$ identity matrix.
Similarly to the authors of Ref. [34], we prepare our system in a charge
density wave state in the system with open boundaries,

$\ket{\Psi_{0}}=\left(\prod_{l=1}^{L/2}c_{2l}^{\dagger}\right)\ket{\text{vac}}.$
(21)

In Fig. 2 (left) we benchmark the dynamics of the half-chain entanglement
entropy $S_{L/2}(t)$ with the results of Ref. [34], for two values of $\gamma$
finding perfect agreement. In the right panel, we demonstrate convergence with
respect to the number of polynomials $N_{p}$, for a given time-step $\delta
t=0.1J^{-1}$. We validate the decreasing of the entanglement due to the
presence of non-Hermitian Skin effect [34]. In addition to the entanglement
entropy, we have also validated our results against other metrics presented in
Ref. [34]. For instance, we equally observe the initial charge-density wave
state rapidly evolving into a state with charge accumulation at one boundary
due to non-reciprocal hopping. In Fig. 3 (top panels) we plot the space-time
dynamics of particle density as well as a cut at long-times, describing the
steady-state density profile along the chain for different system sizes. We
see that for short systems, as compared to the single particle wave function
localisation length, the accumulation takes the form of a domain wall, while
upon increasing system size a finite slope emerges, which we have checked to
vanish exponentially with $L$. In Fig. 3 (bottom panels) we plot the dynamics
of the local current [34] defined as

$I_{\ell}=\dfrac{Ji}{2}\left\langle
c_{\ell+1}^{\dagger}c_{\ell}-c_{\ell}^{\dagger}c_{\ell+1}\right\rangle.$ (22)

We see that, consistently with the density plot, a finite (negative) current
flows in the bulk of the chain, for sufficiently long systems and long times,
while at the boundaries the current vanishes due to the localised charge in
the domain walls. This current is a feature of this non-equilibrium steady,
given that it is not present in the ground-state of the Hatano-Nelson model.
Besides, in contrast to the Hermitian scenario, the density profiles exhibit a
spatial gradient within the bulk of the chain. This non-trivial spatial
distribution of density arises from the single-particle skin effect.
Additionally, in contrast to the Hermitian scenario, non-Hermiticity permits
spatial variations in the local charge current while the density profile
remains time-independent. This phenomenon is facilitated by the fact that the
continuity equation associated to charge conservation in non-Hermitian systems
takes the form

$\partial_{t}\left\langle
c^{\dagger}_{\ell}c_{\ell}\right\rangle+\left(I_{\ell}-I_{\ell-1}\right)=\mathcal{T}_{\ell},$
(23)

where the additional term $\mathcal{T}_{\ell}$ corresponds to the sink/source
of particles due to coupling with the environment equal to,

$\mathcal{T}_{\ell}=-\dfrac{\gamma}{J}\sum_{n=0}^{L-2}\left(\left\langle\left\\{c^{\dagger}_{\ell}c_{\ell},I_{n}\right\\}\right\rangle-2\left\langle
I_{n}\right\rangle\left\langle
c^{\dagger}_{\ell}c_{\ell}\right\rangle\right).$ (24)

This term is unique to systems evolving under non-unitary dynamics generated
by a non-Hermitian Hamiltonian and has important consequences for the
transport properties of these systems [34, 81], as we will further discuss
below.

Figure 3:  Left - Time and spatial dependence of the particle density (Top
plot) and charge current (bottom) profile for a system size of 100 sites.
Right - Spatial dependence of the particle density (top) and charge current
(bottom) in the steady for different system sizes. Other parameters:
$\gamma=-0.8J$ Figure 4: Time evolution of the magnetisation profile of the
Hatano-Nelson model for different values of the non-Hermitian parameter
$\gamma$. The total system length corresponds to $L=256$. The dashed white
lines correspond to the effective penetration length.

###  4.2 Domain Wall Melting for Hatano-Nelson

In this Section, our focus is on exploring the impact of non-Hermiticity on
the temporal evolution of particle density, current, and entanglement profiles
in an HN system initialised in a domain-wall state, $\ket{\rm
DW}=\ket{111\cdots 10\cdots 000}$222Using the Jordan-Wigner [82, 83]
transformation the Hatano-Nelson model can be viewed in the spin-$1/2$
language as an XX chain with a non-reciprocal XX exchange term.. This
configuration has been the subject of extensive investigations in the unitary
case [84, 85, 86, 87], as it exemplifies the distinctive characteristics of
non-equilibrium dynamics. In addition, it has inspired the development of
generalised hydrodynamics (GHD) [88, 89, 90], which facilitates precise
calculations of charge and current profiles using a hydrodynamic description.
Quantum correlations and entanglement entropy can also be calculated under an
extension of this framework, quantum GHD [91]. In the non-Hermitian case, the
dynamics of an initial domain wall has been less studied, even in the simple
non-interacting HN model.

We start considering the time-evolution of the particle density profile under
the HN dynamics, that we plot in Fig. 4 for increasing values of the non-
reciprocal coupling $\gamma$. In the unitary case ($\gamma=0$) a clear light
cone is visible, corresponding to ballistically propagating quasiparticles. In
the non-Hermitian case a light cone is still visible at short times, at least
for $\left|\gamma\right|<J$, when the particle density satisfies a scaling
function $n\left(\ell/t\right)=f\left(\ell/(v_{\text{eff}}\;t\right)$, where
$v_{\rm eff}$ is an effective velocity. As the non-reciprocal coupling
increases, the light cone shrinks more and more, up to $\gamma=J$,
corresponding to the exceptional point of the HN, at which the domain-wall
state remains stable. We now focus on the short-time dynamics and discuss the
origin of the velocity renormalisation. From the numerical data we obtain (see
Fig. 5 (top left))

Figure 5: Top - Quasi-particle effective velocity as a function of non-
reciprocal coupling (left plot). Particle density profile as a function of
space for a fixed time (right plot). Bottom – Spatial profile of charge
current at different times for $\gamma=0.2J$ (right plot). The plot in the
bottom left shows the current profile as a function of space and time for
$\gamma=0.8J$. Out of the ballistic region, the renormalised GHD equations are
no longer valid. Other parameters: $L=256$.

$v_{\text{eff}}=J-\gamma.$ (25)

According to this, the domain wall spreads less for higher values of the non-
Hermiticity parameter, and the renormalised velocity vanishes at $\gamma=J$.
We note that this simple formula is distinct and is not related to those
derived for the propagation of wave packets through a non-Hermitian medium
[92, 93], or to the one obtained from the energy dispersion relation of the
model, which would be proportional to $\sqrt{J^{2}-\gamma^{2}}$ (see Fig. 5
(top left) and consult Appendix. B for the explicit derivation). Moreover, it
depends on the initial conditions chosen, as if the system was initialised in
the state $\ket{\cdots 000111\cdots}$ instead, the formula should be replaced
by $v_{\text{eff}}=J+\gamma$. The reduction in propagation velocity results in
a suppression of correlations relative to the Hermitian scenario, a phenomenon
previously observed in other non-Hermitian systems [33, 34, 94].

We can provide a simple argument to justify the renormalisation of the
velocity due to non-hermiticity. The renormalised velocity can be obtained by
using a local continuity equation for the particle density (Eq. 23) by working
the term $\mathcal{T}_{\ell}$. In the non-interacting limit, this term can be
expanded using Wick’s theorem in a local and non-local term,

$\begin{aligned}
\mathcal{T}_{\ell}=&-\dfrac{\gamma}{J}\left(I_{\ell}+I_{\ell-1}\right)\left(1-2\left\langle
c_{\ell}^{\dagger}c_{\ell}\right\rangle\right)+\\\
&-i\gamma\sum_{n\neq\\{\ell,\ell-1\\}}^{L-2}\left[\left\langle
c_{\ell}^{\dagger}c_{n+1}\right\rangle\left\langle
c_{n}^{\dagger}c_{\ell}\right\rangle-\left\langle
c_{n+1}^{\dagger}c_{\ell}\right\rangle\left\langle
c_{\ell}^{\dagger}c_{n}\right\rangle\right]\end{aligned}.$ (26)

The second term of this equation is highly non-local, and it should be zero
outside the light-cone region. We proceed by writing the local continuity
equation for a site $\ell$, outside the light-cone region, for an instant of
time before the arrival of the magnetisation wavefront. Under these conditions
$I_{\ell}=0$ and $\langle c^{\dagger}_{\ell}c_{\ell}=0\rangle$ for $\ell>0$ or
$I_{\ell-1}=0$ and $\langle c^{\dagger}_{\ell}c_{\ell}=1\rangle$ for $\ell<0$.
Furthermore, the non-local correlator vanishes, and so the continuity equation
is approximately given by

$\displaystyle\partial_{t}\left\langle n_{\ell}\right\rangle$
$\displaystyle-\left(1-\dfrac{\gamma}{J}\right)I_{\ell-1}$
$\displaystyle=0,\quad\ell>0,$ (27) $\displaystyle\partial_{t}\left\langle
n_{\ell}\right\rangle$
$\displaystyle+\left(1-\dfrac{\gamma}{J}\right)I_{\ell}$
$\displaystyle=0,\quad\ell<0.$

The $\pm$ sign reflects the different directions of propagation. This formula
holds only before the particle density wavefront reaches this site. As after,
there is the development of off-diagonal correlations within the light-cone
area.

Since at short times we can still identify a sharp light cone, it is tempting
to try to formulate a hydrodynamic description for the HN model. Although some
advances have been made within the framework of Linblad dynamics [95, 96],
formulating a hydrodynamic description for the non-Hermitian variant of the
XXZ model remains extremely challenging, even in the absence of interactions.
The nonlinearity of the equations of motion results in the loss of most local
conservation laws. Surprisingly, as we show for the Hatano-Nelson model in
Fig. 5, the hydrodynamic equations still hold in the initial time. To make a
concrete comparison between the numerics and the analytical predictions, we
proceed in a somewhat phenomenological way and incorporate the renormalised
velocity into the hydrodynamic expressions derived for the Hermitian case
[97].

With this in mind, the particle density in the spatial interval,
$-t\leq\ell\leq t$, is given by the following expression,

$n\left(\ell,t\right)=\dfrac{1}{\pi}\arccos\left(\dfrac{\ell}{v_{\text{eff}}\;t}\right).$
(28)

which perfectly matches the results of the full numerical calculations (see
Fig. 5 top right).

Figure 6: Entanglement entropy for the subsystem $\left[-L/2,\ell\right]$ in
different times. From right to left, the non-Hermiticity is $\gamma=0,0.2J$
and $0.4J$.

This is also applicable to the charge current, as we show in Fig. 5 (bottom
left). Finally, using the quantum GHD formalism [91, 97], it can likewise be
extended to the entanglement entropy of the region $\left[-L/2,\ell\right]$,

$I\left(\ell,t\right)=\dfrac{1}{\pi}\sqrt{1-\left(\dfrac{\ell}{v_{\text{eff}}\;t}\right)^{2}}.$
(29)
$S\left(\ell,t\right)=\dfrac{1}{6}\ln\left[v_{\text{eff}}\;t\left(1-\left(\dfrac{\ell}{v_{\text{eff}}\;t}\right)^{2}\right)^{3/2}\right]+c_{1},$
(30)

where $c_{1}\simeq 0.4785$ [98] as shown in the Fig. 6. Similarly to the
unitary case, the entanglement increases as a result of the development of
quantum correlations within a well-defined spatial region [99]. Ballistic
transport of the charge current prevails, leading to a region that expands
linearly over time without significant entropy production. This phenomenon is
characterised by growth $\ln(t)$, which is emblematic of a local quantum
quench protocol [100, 101].

In the unitary regime, ballistic propagation ceases within temporal scales
commensurate with the entire system size, which is applicable to systems of
finite dimensions. The non-reciprocal hopping stabilises the domain wall, thus
preventing it from melting. With a finite magnitude of non-reciprocal
coupling, the particle density wavefront is constrained to penetrate only to a
prescribed depth. This maximum penetration depth corresponds to the total
system size, $L$, for, $\gamma=0$ and tends to zero at the exceptional point
$\gamma=J$. It is as if the system length is renormalised to an effective one
given by

$L^{\ast}\sim\dfrac{J-\gamma}{J+\gamma}L.$ (31)

This length is marked by the white dashed lines in the Fig. 4. We just stress
that this reasoning is only valid for $\gamma>0$, since for $\gamma<0$ the
magnetisation wavefront necessarily reaches the system boundary. The
expression suggests that faster quasiparticles with velocity
$v^{\ast}=J+\gamma$ cease the ballistic propagation of the magnetisation
wavefront upon hitting the physical boundary. Thus, the maximum propagation
velocity allowed in the system increases compared to the Hermitian case, since
particles can hop from right to left with velocity $J+\gamma$ (assuming
$\gamma>0$). Following the initial propagation of the magnetisation wave
front, the system reaches a steady state characterised by the presence of a
charge current traversing the two domains. This current emerges due to a flux
of particles sourced from the environment, which are subsequently annihilated,
as depicted in Fig. 5. However, contrary to the setting described by Kawabata
[34], non-Hermiticity leads to the spatial suppression of this current.

##  5 Application to Non-Hermitian Many-Body Systems

We now consider an application of the Faber polynomial method to a full non-
Hermitian many-body problem. In this case, one can simplify the evaluation of
the evolution operator, while the cost of storing the state is still
exponential in system size since, as we stress again, in the Faber polynomial
method there is no approximation on the state which is fully represented in
the basis.

Figure 7: Temporal and spatial dependence of the magnetisation profile,
$\langle S_{\ell}^{z}\rangle$, for different values of $\Delta$ in the
interacting Hatano-Nelson model. Other parameters: $\gamma=0.8J$ and $L=24$.
Figure 8: Long-time behaviour of the magnetisation profile shown in Fig. 7,
$\langle S_{\ell}^{z}\rangle_{\infty}$, for different values of $\Delta$. We
see that at $\Delta=0$ an emergent domain wall is formed, which is stable for
small $\Delta$ (left panel). Increasing $\Delta$ the system develops a
potential drop in the middle of the chain, akin to diffusive dynamics, which
then further develops an oscillating patter as $\Delta$ increases (right
panel). Other parameters $\gamma=0.8J$.

###  5.1 Magnetisation Dynamics in the Interacting Hatano-Nelson Model

In this section, we study the effects of interactions in the dynamics of the
magnetisation profile and spin current on a non-Hermitian XXZ chain with a
non-reciprocal XX exchange term.

$\mathcal{H}=-\sum_{\ell=0}^{L-2}\left[\frac{\left(J+\gamma\right)}{2}S_{\ell}^{+}S_{\ell+1}^{-}+\frac{\left(J-\gamma\right)}{2}S_{\ell+1}^{+}S_{\ell}^{-}+\Delta
S_{\ell}^{z}S_{\ell+1}^{z}\right],$ (32)

where $J$ is a XX exchange term between neighbouring spins, $\gamma$ induces
an imbalance between the propagation of left/right magnetic excitations and
$\Delta$ is an Ising like exchange. This system can also be viewed as an
interacting Hatano-Nelson model by performing the Jordan-Wigner transformation
[82, 83]. The system is prepared at time zero in an unentangled Néel state,

$\ket{\Psi(t=0)}=\ket{\uparrow,\downarrow,\cdots,\downarrow,\uparrow},$ (33)

which is an eigenstate of the Ising part of the model, so, one can expect that
when $\Delta\gg J,\gamma$ the Néel order is preserved. For $\Delta=0$, when
the model reduces to the non-recriprocal XX chain (or Hatano-Nelson in
fermionic language), it is known that the initial Néel state gives rise to a
domain wall state at long times as all magnetic excitations are transported to
one of the edges of the system [34], just as in the Hatano-Nelson model which
exhibits charge accumulation at a boundary. The extent of this proximity is
governed by the degree of non-Hermiticity, which is regulated by the parameter
$\gamma$. Here, we are interested in understanding the role of interactions in
the magnetisation dynamics and the stability of this emergent domain wall. In
Fig. 7 we plot the spatial-temporal dynamics of the magnetisation $\langle
S^{z}_{\ell}(t)\rangle$ for an increasing value of $\Delta$. We see that at
short times there is a rapid reshuffling of magnetic excitations driven by the
non-reciprocity, towards a boundary accumulation. Interactions compete with
this process and tend to preserve the initial antiferromagnetic pattern at the
expense of boundary accumulation, as we see well in the right panel of Fig. 7.
To better characterise the long-time dynamics, we compute the average
magnetisation profile

$\left\langle
S^{z}_{\ell}\right\rangle_{\infty}=\mbox{lim}_{T\rightarrow{\infty}}\frac{1}{T}\int_{\tau^{\ast}}^{\tau^{\ast}+T}dt\;\langle
S^{z}_{\ell}(t)\rangle,$ (34)

where $\tau^{\ast}$ is the initial time corresponding to the reshuffling of
magnetic until there is a stable boundary accumulation. We plot in Fig. 8
$\left\langle S^{z}_{\ell}\right\rangle_{\infty}$ for different values of
$\Delta$. We see that the emergent domain-wall state generated by the non-
reciprocal exchange is stable at weak interactions (left panel). However, as
$\Delta$ increases, a novel region emerges that interpolates between the two
magnetic domains (right panel). In particular, we see that the system develops
a _magnetisation drop_ , similar to a potential drop in systems that show
diffusive transport. Upon further increasing $\Delta$ we see the emergence of
antiferromagnetic correlations on top of this magnetisation slope, which, as
expected, are frozen by the large interaction and do not decay away. The
result above is particularly intriguing as it suggests that non-reciprocal
coupling and interaction collaborate to establish a current-carrying steady
state in the centre of the system: the former driving the formation of a
domain wall that acts as source and drain, the latter providing the necessary
scattering term to dissipate and establish a finite average current. In Fig. 9
we plot the spatio-temporal dynamics of the local current profile as a
function of interaction and non-reciprocity. Comparable findings have been
reported [102] for the ground state of the interacting Hatano-Nelson model
with nearest-neighbour repulsion. In that study, they also noted that the
initial magnetisation profile (referred to as the real-space Fermi surface) is
disrupted by interactions that induce a Néel order. In this work, we
demonstrate that such a phenomenon can also be dynamically generated by the
non-unitary time-evolution.

Similarly to the non-interacting case, there is a current in this
interpolating region, which satisfies the same continuity equation as in Eq.
27. As clearly seen in Fig. 9, once again there is a competition between non-
reciprocity and the interaction parameter $\Delta$ in defining the size of
this intermediate region, where it is possible for the particle to enter from
the environment and give rise to this current. As one is looking at small
system sizes, it is not possible to create a stable steady current (like in
the non-interacting case). Due to this, we see oscillations in the direction
of the central current, and thus there is current flowing in both directions.
However, in the extreme case where $\gamma=J$, this oscillatory behaviour
ceases to exist, and the current has a fixed direction. When $\Delta$ is much
greater than $J$, the current disappears in the bulk region, where the Néel
order is maintained.

Figure 9: Temporal and spatial dependence of the local current profile,
$\langle I_{n}\rangle$. Other parameters: $\gamma=0.8J$ and $L=24$.

###  5.2 Effect of interaction in the Domain-Wall melting for Hatano-Nelson

Finally, we conclude this section by discussing the effect of interactions on
the non-Hermitian problem discussed in Sec. 4.2, namely an initial domain wall
state.

In the previous section, we have focused on the domain wall melting for the
non-interacting non-reciprocal Hatano-Nelson model, or in the spin analogue,
the XX chain as $\Delta=0$. For a conventional Hermitian XXZ spin chain, the
domain wall is only stabilised when the Ising exchange is greater than the XX
exchange, $\left|\Delta\right|>J$. In contrast, $\left|\Delta\right|<J$, the
domain wall melts, with a ballistic propagation of the magnetisation wavefront
[103]. The Heisenberg point, $\Delta=J$, is special given the existence of the
extra spin SU(2) symmetry, which allows for superdiffusive behaviour of the
spin current [86, 104]. Nevertheless, there is no Heisenberg point in the spin
version of the Hatano-Nelson model, as the non-Hermiticity explicitly breaks
the spin SU(2) symmetry. We observe that the interactions contribute to
prevent the domain wall from melting, as shown in Fig. 10. In a certain sense,
non-Hermiticity and interactions help to preserve the initial magnetic order,
which otherwise would be eroded by the dynamics. We have benchmarked this
result with matrix product of states (MPS) calculations presented in the
Appendix C.

Figure 10: Spatial and temporal dependence of the magnetisation profile for
different values of the Ising exchange parameter. Other parameters: spins
$\gamma=0.2J$ and $L=24$. Figure 11: Magnetisation profile in function of the
lattice site and time. Left - Dynamics with the jump operator,
$\sqrt{\left|\gamma\right|}\left(S^{-}_{\ell}-i\text{sgn}\left(\gamma\right)S^{-}_{\ell+1}\right)$.
Right - Dynamics with the jump operators,
$\sqrt{\gamma}S^{+}_{\ell}S^{-}_{\ell+1}$. The system was initially prepared
in a Néel State (Eq. 33). Other parameters: $\Delta=0$, $\gamma=0.8J$ and
$L=20$.

##  6 Application to Quantum Jumps Unravelling

In this Section, we combine the Faber polynomial technique with a high-order
Monte Carlo Wave Function algorithm [11, 9] to investigate the Quantum Jumps
unravelling of the Lindblad master equation, discussed in Sec. 2. In
particular, we address the stochastic Schrödinger equation by propagating the
initial state with the Faber polynomial method up to the time instance
($\tau$) when a quantum jump occurs. Therefore, within the time interval
$t\in\left[t_{0},t_{0}+\tau\right]$, the state evolves purely non-unitary
according to

$\ket{\psi(t)}=\dfrac{e^{-i\mathcal{H}t}\ket{\psi(t_{0})}}{\left\|e^{-i\mathcal{H}t}\ket{\psi(t_{0})}\right\|}.$
(35)

The time instance $\tau$ is obtained via the standard higher-order Monte Carlo
wave function technique. It corresponds to the specific time at which the norm
of the state equates to a random variable, $r$, drawn from a uniform
distribution over the interval $\left[0,1\right]$. Consequently, $\tau$ is
implicitly defined by the following equation,

$r=1-\bra{\psi(t_{0})}e^{i\mathcal{H}^{\dagger}\left(\tau-
t_{0}\right)}e^{-i\mathcal{H}\left(\tau-t_{0}\right)}\ket{\psi(t_{0})}.$ (36)

The quantum jump is applied by first selecting the quantum jump channel in
accordance to their probability,

$p_{\mu}=\dfrac{\left\langle
L^{\dagger}_{\mu}L_{\mu}\right\rangle}{\sum_{\mu}\left\langle
L^{\dagger}_{\mu}L_{\mu}\right\rangle},$ (37)

where the average is taken with the state, $\ket{\psi\left(t+\tau\right)}$.
Then, the post-jump state, $\ket{\psi\left(t+\tau^{+}\right)}$, is obtained by
applying the chosen jump operator, $L_{\alpha}$,

$\ket{\psi\left(t+\tau^{+}\right)}=\dfrac{L_{\alpha}\ket{\psi\left(t+\tau\right)}}{\sqrt{\bra{\psi\left(t+\tau\right)}L^{\dagger}_{\alpha}L_{\alpha}\ket{\psi\left(t+\tau\right)}}}$
(38)

This algorithm gives access to the full Lindbladian dynamics, by averaging
over the propagated quantum trajectories. Further, it also provides access to
the dynamics under continuous monitoring [20, 9]. This is achieved by
following the many-body quantum trajectory and computing, for example, non-
linear functions of the state. One such function is the entanglement entropy
of quantum trajectories, as we shall discuss below.

###  6.1 Magnetisation and Entanglement Dynamics in Monitored Spin Chains

As an application, we consider a quantum spin chain described by the Hermitian
XXZ model with Hamiltonian

$\mathcal{H}=-\sum_{\ell=0}^{L-2}\left[\frac{J}{2}S_{\ell}^{+}S_{\ell+1}^{-}+\frac{J}{2}S_{\ell+1}^{+}S_{\ell}^{-}+\Delta
S_{\ell}^{z}S_{\ell+1}^{z}\right],$ (39)

where $J$ is the XX exchange term between neighbouring spins and $\Delta$ an
Ising like exchange. We compare the dissipative dynamics generated by two
different types of jump operators. A first set of jumps that we consider
describe next-neighbour decoherence of spin excitations along the chain and
take the form

$\displaystyle L_{0}$ $\displaystyle=\sqrt{\left|\gamma\right|}S^{-}_{0},$
(40) $\displaystyle L_{1+\ell}$
$\displaystyle=\sqrt{\left|\gamma\right|}\left(S^{-}_{\ell}-i\text{sgn}\left(\gamma\right)S^{-}_{\ell+1}\right),\;\ell=\\{0,\cdots,L-2\\}$
$\displaystyle L_{L}$ $\displaystyle=\sqrt{\left|\gamma\right|}S^{-}_{L-1}.$

Interestingly, with these choice of jump operators the non-Hermitian
Hamiltonian associated to the no-click limit turns out to be given by a spin
version of the many-body Hatano-Nelson model. Indeed, the following
straightforward calculation yields the non-Hermitian Hamiltonian

$\begin{aligned}
\mathcal{H}_{\text{eff}}&=\mathcal{H}-\frac{i}{2}\sum_{\ell}L_{\ell}^{\dagger}L_{\ell}\\\
&=\mathcal{H}-\frac{\gamma}{2}\sum_{\ell=0}^{L-2}\left(S_{\ell}^{+}S^{-}_{\ell+1}-S_{\ell+1}^{+}S^{-}_{\ell}\right)-i\left|\gamma\right|\sum_{\ell=0}^{L-1}S_{\ell}^{+}S^{-}_{\ell}\end{aligned}.$
(41)

The last term does not affect the dynamics, as it just an overall background
decay. Although this model can be connected to the fermionic version of the
Hatano-Nelson model, the Linblad equation is not quadratic due to the Jordan-
Wigner strings in the terms of the form $S^{-}\rho S^{+}$.

We compare the dynamics generated by the quantum jumps above with the one
induced by a different set of jump operators that create a spin-flip
excitation from site $n+1$ to site $n$. These are read as follows,

$L_{\ell}=\sqrt{\gamma}S_{\ell}^{+}S_{\ell+1}^{-},$ (42)

with $\ell=\left\\{0,\cdots,L-2\right\\}$. Previous studies have shown that
these operators induce a phenomenon known as the Liouvillian skin effect [24,
105, 106, 107]. Specifically, there is an exponential localisation of the
Liouvillian modes at the boundaries of the system. This becomes evident when
the Linblad is projected onto the one-particle sector, where, at $J=0$ and
$\Delta=0$, it simplifies to an effective Hatano-Nelson Hamiltonian tuned to
its exceptional point [105],

$\mathcal{H}=\gamma\sum_{\ell=0}^{L-2}S^{+}_{\ell}S^{-}_{\ell+1}.$ (43)

The non-click Hamiltonian differs from the Hatano-Nelson model, representing
an XXZ chain with an imaginary Ising exchange term and boundary imaginary
magnetic fields,

$\begin{aligned}
\mathcal{H}_{\text{eff}}&=-\frac{J}{2}\sum_{\ell=0}^{L-2}\left[S_{\ell}^{+}S_{\ell+1}^{-}+\text{h.c}\right]-\left[\Delta-i\dfrac{\gamma}{2}\right]\sum_{\ell=0}^{L-2}S_{\ell}^{z}S_{\ell+1}^{z}\\\
&\quad+\frac{i\gamma}{4}\left(S_{0}^{z}-S_{L-1}^{z}\right)-\frac{i\gamma}{8}\left(L-1\right)\end{aligned}.$
(44)

In Fig. 11, we plot the magnetisation dynamics starting from a product state
and evolving under the two types of dissipative evolutions. In the left panel,
we plot the dynamics generated by Eq. 40. Under this set of jump operators,
manifestations of non-reciprocity and charge accumulation are still present in
the transient dynamics [108, 109]. We note that these are the same jump
operators considered in [109] with the phase, $\phi$, defined by them, tuned
to get the non-reciprocal regime. However, in the long time limit, a generic
state converges to a configuration with zero excitations,
$\ket{\downarrow\cdots\downarrow}$ [109] (see the left plot of Fig. 11). This
is attributed to the recycling terms of the form $c_{n}\rho c^{\dagger}_{n}$,
which effectively remove particles from the system. Additionally, the system
has another dark state within the one-magnon sector, which only plays a role
in decelerating the relaxation dynamics towards the fully polarised down state
[109]. The final state obtained under this Linbladian does not resemble at all
the magnetisation profile that one would have obtained in the no-click limit.

We then focus on the jump operators in Eq. 42 and plot in Fig. 11 (right
panel) the resulting magnetisation dynamics. For the sake of simplicity, we
consider the dynamics with $\Delta=0$. The dynamics driven by this model will
facilitate the accumulation of spins in an up state at the left extremity of
the chain. The imaginary Ising exchange term is responsible for diminishing
the state’s norm when adjacent spins are antialigned, precipitating a quantum
jump that propagates a spin towards the left edge. This phenomenon is depicted
in the spatial and temporal evolution of the magnetisation profile shown in
Fig. 11. In contrast to the Hatano-Nelson model, this model does not converge
to a steady state characterised by zero excitations. This distinction is
attributable to a difference in symmetry. The Hatano-Nelson model exhibits
only weak U(1) symmetry, whereas the present model possesses strong U(1)
symmetry, thereby ensuring that the system state is kept within the initial
magnetisation sector at all times.

Figure 12: Top - Time dependence of the conditional entanglement entropy for
different system sizes with $\gamma=0.8J$. On the left, the dynamics
corresponds to the jump operator
$\sqrt{\left|\gamma\right|}\left(S^{-}_{\ell}-i\text{sgn}\left(\gamma\right)S^{-}_{\ell+1}\right)$,
while on the right to $\sqrt{\gamma}S^{+}_{\ell}S^{-}_{\ell+1}$. Bottom Left:
Time evolution of the entanglement entropy for $\gamma=0.1J$. The inset
corresponds to the fit of the steady-state entanglement entropy to the law
$S_{L/2}\left(\infty\right)=a_{0}+a_{1}L$, with $a_{0}=\left(0.145\pm
0.002\right)$ and $a_{1}=\left(0.18\pm 0.03\right)$. Bottom Right - Steady-
state entanglement entropy for different system sizes as a function of the
non-hermitian parameter. These systems were studied under the jump operators
$\sqrt{\gamma}S^{+}_{\ell}S^{-}_{\ell+1}$. Other parameter: $\Delta=0.0J$.

This disparity in symmetry radically changes the temporal dynamics of the
conditionally averaged entanglement entropy, which is defined as the mean
entanglement entropy across all conceivable quantum trajectories,

$\bar{S}_{\ell}\left(t\right)=\int\mathcal{D}\xi_{t}\;\mathcal{P}\left(\xi_{t}\right)S_{\ell}\left(\xi_{t}\right).$
(45)

with $S\left(\xi_{t}\right)$ is the entanglement of a given quantum trajectory
that evolves according to Eq. 1.

For early times, the entanglement entropy in both models increases linearly.
Without non-Hermiticity, the system evolves under standard hermitian dynamics,
leading to the entanglement entropy saturating at a value proportional to the
system’s volume as it locally thermalises [110, 111]. However, in the presence
of non-Hermiticity and quantum jumps, this behaviour may drastically alter. In
the Hatano-Nelson model with jumps, the entanglement entropy drops to zero
after an initial period determined by the measurement rate $\gamma$. Each jump
causes a spin-flip, driving the system to a state devoid of magnetic
excitations. This occurs regardless of the system’s total size, resulting in a
trivial entanglement area law for any nonzero value of $\gamma$. The no-click
limit of the Hatano-Nelson model also has this area law scaling of the
entanglement entropy; however, this is driven by the single-particle skin
effect[34].

In contrast, the model described by Eq. 44 does not relax to a zero excitation
state, as seen in Fig. 12. The strong U(1) symmetry confines the dynamics to
the magnetisation sector of the initial state333The dynamics conserves the
total number of excitations in the initial state as the initial state is
eigenstate of $\sum_{n=0}^{L-1}S^{z}_{n}$ and the Hamiltonian in Eq. 44 as
U(1) symmetry.. Initially, the system has a linear growth of the entanglement
entropy, reminiscent of the unitary evolution. However, similarly to the
Hatano-Nelson model, the steady-state supports an area law entanglement for a
positive value of $\gamma$, as seen in the left plot of Fig.12. This can be
understood through the spectral properties of the no-click Hamiltonian (Eq.
44). The spectrum is always complex-valued, and so, in the no-click limit, the
steady-state corresponds to right-eigenstate with the slowest decaying mode in
the zero magnetisation sector. In particular, for certain values of
$\gamma/J$, the imaginary component of the spectrum is gapped. Thus, the
entanglement entropy inevitably follows an area law in the long-time limit,
similar to the ground state of a gapped Hamiltonian [112, 32]. The imaginary
gap ($\Delta_{\rm Im}$) in the zero magnetisation sector can be analytically
obtained in the limit of $\gamma\gg J$,

$\Delta_{\rm
Im}=-i\dfrac{\gamma}{2J}+\mathcal{O}\left(\dfrac{J}{\gamma}\right)$ (46)

The measurement apparatus effectively disentangles the system, pushing
excitations towards the left boundary, and preventing the formation of long-
range correlations. It is harder to completely confirm this for smaller values
of $\gamma$, where the unitary evolution dominates both the non-Hermitian and
stochastic terms. This mainly affects systems with a smaller dimension, where
finite-size effects are substantial, since the terms proportional to $L^{-n}$
with $n>0$ in the entanglement entropy cannot be ignored. For example, for
$\gamma=0.1J$, one cannot extrapolate the true entanglement scaling, as the
observed linear growth might be an artefact of small system sizes.
Nevertheless, the remaining data points clearly indicate the collapse of all
system sizes to the same value, thus revealing the area law nature of the
entanglement entropy.

##  7 Conclusions

Throughout this study, we have successfully used the Faber polynomial method
to characterise the non-unitary dynamics of both non-interacting and
interacting Hatano-Nelson models. Additionally, we have seamlessly integrated
this approach with a high-order Monte Carlo Wave Function algorithm to
rigorously examine both Lindblad and continuous monitoring dynamics.

In the non-interacting problem, we provided the first numerical evidence
supporting the existence of a valid hydrodynamic description for the melting
dynamics of a domain wall in the presence of non-Hermitian terms. This finding
encourages further developments to properly formalise a theory of generalised
hydrodynamics applicable to non-Hermitian systems.

Our study also reveals an intriguing competition between the Ising exchange
term, which tends to preserve the initial Néel order, and the non-reciprocal
XX coupling, which tends to form a domain-wall ordering. For comparable values
of $\Delta$ and $J$, the interaction allows the formation of an intermediate
region that interpolates between the two magnetic domains, allowing the flow
of current. However, we could not reach considerable system sizes to determine
if this region could support a non-equilibrium steady-state current flowing in
only one direction, as seen in the non-interacting case. It is clear, however,
that this cannot be the case for $\Delta\gg J$, as the dynamics preserve the
initial magnetic ordering, and the current can only exist in the interpolating
region between the Néel-ordered domain and the ferromagnetic one generated by
the non-reciprocal coupling. Conversely, our work shows that interactions and
non-reciprocity help preserve the initial magnetic order when the system is
initialised in a domain-wall setup, a result consistent with both the non-
interacting non-Hermitian problem and the interacting Hermitian case.

This study offers additional insights into the entanglement transition in
quantum spin chains exhibiting Non-Hermitian or Liouvillian skin effects. In
the Hatano-Nelson model with quantum jumps, we found that the area law
behaviour of entanglement entropy remains for any nonzero $\gamma$, similar to
the no-click limit. However, the sources are different: in the no-click limit,
the area law stems from the single-particle skin effect, while in the
monitored stochastic trajectory, it is due to the quantum jumps that relax the
system into the fully polarised down spin-state, a product state. On the other
hand, the dynamics with two-body jump operators, $L_{\ell}\propto
S^{+}_{\ell}S^{-}_{\ell+1}$, allow for a non-equilibrium steady state with a
magnetisation profile resembling the no-click limit of the Hatano-Nelson
model. Furthermore, the average entanglement entropy still follows an area law
for finite values of the ratio $\gamma/J$ in conditional dynamics. The
measurement apparatuses effectively disentangle the state, suppressing the
volume law that would otherwise be generated in unitary dynamics. This work
extends previous studies that focused on the one-particle sector [105, 34].

Overall, our results support the utility and applicability of Faber
polynomials in various research domains. This encompasses investigations into
measurement-induced phase transitions, open quantum systems, and the
exploration of purely non-unitary dynamics governed by non-Hermitian
Hamiltonians. Faber polynomials can also be used for computing general
spectral properties of non-Hermitian Hamiltonians, potentially replicating the
role of Chebyshev polynomials within the Kernel polynomial method [64]. This
would complement the existing Non-Hermitian Kernel polynomial method [113,
114], which still uses Chebyshev polynomials, but relies on hermitianisation
techniques, at the cost of working in a vector space with twice the original
dimension. Furthermore, the Faber polynomial approach could potentially be
combined with MPS, similar to the developments already made in the Hermitian
case [115, 116, 117, 118].

## Acknowledgements

##### Funding information

R.D.S acknowledges funding from Erasmus$+$ (Erasmus Mundus programme of the
European Union) and from the Institute Quantum-Saclay under the project
_QuanTEdu-France_. We acknowledge Collège de France IPH cluster where the
numerical calculations were carried out.

##  Appendix A Further Details on Faber Polynomials

As discussed previously, Faber polynomials serve as a polynomial basis to
represent complex-valued functions that are analytic within the domain
$\mathcal{D}$. These are generated by a conformal mapping $\xi(w)$ that maps
the complement of a closed disk of radius $\rho$ to the complement of
$\mathcal{D}$,

$\dfrac{\xi^{\prime}(w)}{\xi(w)-z}=\sum_{n=0}^{\infty}\dfrac{1}{w^{n+1}}F_{n}\left(z\right),$
(47)

with $F_{n}(z)$ being the $n^{th}$ Faber polynomial generated by the conformal
mapping $\xi(w)$, $z\in\mathcal{D}$ and $w$ is such that $|w|>\rho$. The
existence of such a map, which also satisfies the conditions
$\xi(w)/w\rightarrow 1$ in the limit $\left|w\right|\rightarrow\infty$, is
guaranteed by the Riemann mapping theorem [57, 119]. Furthermore, $\xi(w)$
admits a Laurent expansion at $w=\infty$ of the form,

$\xi(w)=w+\sum_{m\geq 0}\gamma_{m}w^{-m},$ (48)

where $\gamma_{m}\in\mathbb{C}$. Using the Laurent expansion and integrating
in the contour defined around the disk of radius $\rho$, it is straightforward
to check that the Faber polynomials satisfies the following recurrence
relation,

$F_{n+1}(z)=zF_{n}(z)-\sum_{j=0}^{n}\gamma_{j}F_{n-j}(z)-n\gamma_{n},\quad
n>0$ (49)

where $F_{0}(z)=1$. For our purposes, we are interested in using the Faber
polynomials to approximate a given function of our non-Hermitian Hamiltonian,
$f\left(\mathcal{H}\right)$. The domain $\mathcal{D}$, is defined by the
spectrum of the Hamiltonian. Using Eq. 47, the expression in a Faber series is
given by

$f(z):=\sum_{k=-\infty}^{+\infty}c_{k}F_{k}(z).$ (50)

with the coefficients given by the integral,

$c_{n}=\frac{1}{2\pi
i}\int_{|w|=\rho}\frac{f\left(\xi\left(w\right)\right)}{w^{n+1}}dw.$ (51)

As stated in the main text, we perform this integral with $\rho=1$ by properly
rescaling the Hamiltonian. Furthermore, we compute the Faber coefficients
related to an elliptic contour. For this contour, the conformal mapping
reduces to $\xi(w)=w+\gamma_{0}+\gamma_{1}w^{-1}$, where $\gamma_{0}$ is the
centre of the ellipse and $\gamma_{1}=1-b$, with $b$ the major semiaxis. This
maximises the memory efficiency of the algorithm, as the recurrence relation
in Eq. 49 only depends on the two previous polynomials,

$\begin{aligned} F_{0}(z)&=1\\\ F_{1}(z)&=z-\gamma_{0}\\\
F_{2}(z)&=\left(z-\gamma_{0}\right)F_{1}(z)-2\gamma_{1}F_{0}(z)\\\
F_{n+1}(z)&=\left(z-\gamma_{0}\right)F_{n}(z)-\gamma_{1}F_{n-1}(z),\;n\geq
2\end{aligned}.$ (52)

The coefficients $c_{n}$ of the Faber series presented in the main text (Eq.
10) are straightforwardly computed by performing the contour integral in Eq.
51 for the function $f(z)=e^{-i\delta t_{s}z}$. This is done through the use
of the identity [74],

$\exp\left(\dfrac{z}{2}\left[s+\dfrac{a}{s}\right]\right)=\sum_{n=-\infty}^{+\infty}\left(\dfrac{s}{i\sqrt{a}}\right)^{n}J_{n}\left(i\sqrt{a}z\right).$
(53)

The Faber Polynomials have two interesting limits; when $\gamma_{1}=0$ the
conformal mapping corresponds to a circle, thus the Faber Polynomials reduce
to the Taylor polynomials centred around $\gamma_{0}$,

$F_{n}(z)=\left(z-\gamma_{0}\right)^{n},$ (54)

Whereas, in the limit $\gamma_{1}=1$, the conformal mapping evolves the real
line, and so, the Faber Polynomials can be related with the Chebyshev
polynomials by

$\begin{aligned} F_{1}(z)&=T_{1}\left(\dfrac{z-\gamma_{0}}{2}\right),\\\
F_{2}(z)&=T_{2}\left(\dfrac{z-\gamma_{0}}{2}\right)-1,\\\
F_{n}(z)&=T_{n}\left(\dfrac{z-\gamma_{0}}{2}\right),n\geq 2\end{aligned}.$
(55)

##  Appendix B General Features on the Hatano-Nelson Model

In this Appendix, we review the essential characteristics of the non-
interacting Hatano-Nelson model, highlighting the importance of boundary
conditions on the resulting physical properties. Specifically, we examine the
diagonalisation of the non-interacting Hatano-Nelson model (Eq. 18) under both
open and periodic boundary conditions. For OBCs, the Hamiltonian can be
diagonalised through a similarity transformation. This is done using the
$\text{GL}\left(1\right)$ gauge transformation[34],

$\hat{p}_{\ell}^{\dagger}:=e^{\ell\theta}\hat{c}_{\ell}^{\dagger},\quad\hat{q}_{\ell}:=e^{-\ell\theta}\hat{c}_{\ell},$
(56)

with $\theta\in\mathbb{R}$ and both $p_{\ell}$ and $q_{\ell}$ two fermionic
operators satisfying the unusual anticommutation relations:
$\left\\{p^{\dagger}_{\ell},q_{n}\right\\}=\delta_{\ell,n}$,
$\left\\{p^{\dagger}_{\ell},p_{n}\right\\}=e^{2\ell\theta}\delta_{\ell,n}$,
$\left\\{q^{\dagger}_{\ell},q_{n}\right\\}=e^{-2\ell\theta}\delta_{\ell,n}$
and
$\left\\{p_{\ell},q_{n}\right\\}=\left\\{p^{\dagger}_{\ell},q^{\dagger}_{n}\right\\}=0$.
This indicates the biorthogonality of the Hamiltonian’s eigenbasis.

The Hamiltonian can subsequently be expressed in the form

$\mathcal{H}_{\text{HN}}=-\sum_{n=0}^{L-2}\left[\frac{e^{\theta}\left(J+\gamma\right)}{2}p_{n}^{\dagger}q_{n+1}+\frac{e^{-\theta}\left(J-\gamma\right)}{2}p_{n+1}^{\dagger}q_{n}\right]$
(57)

The $\theta$ parameter is chosen so that the final Hamiltonian becomes
hermitian,

$\theta=\frac{1}{2}\log\left(\frac{J-\gamma}{J+\gamma}\right).$ (58)

This reduces the Hamiltonian to

$\mathcal{H}_{\text{HN}}=-\frac{\sqrt{J^{2}-\gamma^{2}}}{2}\sum_{n=0}^{L-2}\left[p_{n}^{\dagger}q_{n+1}+p_{n+1}^{\dagger}q_{n}\right].$
(59)

The Hamiltonian can be diagonalised through a straightforward Fourier
transformation,

$\mathcal{H}_{\text{HN}}=-\sqrt{J^{2}-\gamma^{2}}\sum_{k}\cos\left(k\right)p_{k}^{\dagger}q_{k},$
(60)

where $k=\frac{\pi}{L+1}\;n,n\in\left\\{1,\cdots,L\right\\}$. The quasi-
particles generated by $p^{\dagger}_{k}$ and $q^{\dagger}_{k}$ are
nonorthogonal, as shown by their anticommutation relations. Furthermore, they
correspond to states exponentially localised at the left and right edges of
the chain,

$\displaystyle p_{k}^{\dagger}\ket{\text{vac}}$
$\displaystyle=\sqrt{\frac{2}{L+1}}\sum_{\ell=0}^{L-1}e^{\ell\theta}\sin\left(k\cdot\left(\ell+1\right)\right)c_{\ell}^{\dagger}\ket{\text{vac}},$
(61) $\displaystyle q_{k}^{\dagger}\ket{\text{vac}}$
$\displaystyle=\sqrt{\frac{2}{L+1}}\sum_{\ell=0}^{L-1}e^{-\ell\theta}\sin\left(k\cdot\left(\ell+1\right)\right)c_{\ell}^{\dagger}\ket{\text{vac}}.$

The parameter $\theta$ controls wave-function localisation, inducing a
characteristic length scale $l$,

$l\sim\left(\log\left(\frac{J-\gamma}{J+\gamma}\right)\right)^{-1}.$ (62)

This characteristic length scale emerges due to the skin effect and is
exclusive to the open-boundary condition scenario. Conversely, under PBCs, the
Hamiltonian can be diagonalised using a Fourier transform without employing
the GL gauge transformation,

$\mathcal{H}_{\text{HN}}=-\sum_{k\in\text{FbZ}}\left[J\cos\left(ka\right)+i\gamma\sin\left(ka\right)\right]c_{k}^{\dagger}c_{k}.$
(63)

Unlike the previous case where the eigenstates were confined to the chain
ends, under periodic boundary conditions both the right and left eigenstates
are indistinguishable and delocalised. Furthermore, the spectra is always
complex for all values of $\gamma$.

Figure 13: Comparison between the Faber polynomial expansion and the TEBD
(MPO) algorithm, shown on the left (right). The calculations on the left were
performed with $\gamma=0.0J$, $\Delta=0.5J$ and $L=18$, while on the right we
used $\gamma=0.2J$, $\Delta=0.4J$ and $L=18$.

##  Appendix C Comparison with MPS based Methods

In this Appendix, we compare the Faber Polynomial technique with MPS
calculations for the interacting domain-wall melting problem. First, using a
second-order time-evolution block decimation (TEBD) algorithm [120], and
second, a first-order matrix product operator (MPO) representation of the
time-evolution operator [121]. These methods were implemented using the
ITensor Library [122, 123]. The TEBD algorithm consists in performing a
Suzuki-trotter break-up of the time evolution operator,

$e^{-i\delta t\left(\mathcal{H}_{\rm even}+\mathcal{H}_{\rm odd}\right)}\simeq
e^{-i\delta t\mathcal{H}_{\rm even}}e^{-i\delta t\mathcal{H}_{\rm odd}},$ (64)

where the Hamiltonian in Eq. 32 is decomposed as a sum of two body terms
either acting on the left or right sites, $\mathcal{H}=\mathcal{H}_{\rm
odd}+\mathcal{H}_{\rm even}$. Whilem in the first-order MPO, we use an Euler
expansion of the time-evolution operator, $e^{-i\delta t\mathcal{H}}\simeq
1-i\delta t\mathcal{H}$, and represent the Hamiltonian through an MPO.

MPS calculations could be effective for this problem, as the entanglement
entropy increases logarithmically over time and the domain-wall order is
favoured by non-Hermiticity. However, the TEBD and MPO algorithms face errors
that the Faber algorithm avoids. Firstly, there is a restriction on the
maximum allowed time step. In the TEBD algorithm, the Hamiltonian is
decomposed into smaller two-body terms that can be exactly exponentiated, as
shown in Eq. 64. The approximation neglects the nonzero commutator term
$\left[\mathcal{H}_{\rm odd},\mathcal{H}_{\rm even}\right]$, causing an error
of order $\delta t^{2}$. To minimise this error, the time step must be small
relative to the problem’s energy scale. This is evident in the left panel of
Fig. 13, where longer integration times with a larger time step deviate from
results obtained with Faber polynomials. A comparable error is present in the
MPO algorithm due to the use of a first-order expansion. The Faber algorithm
circumvents this issue by precisely representing the time evolution operator
for any time step with any desired accuracy, though it requires the
computation of the appropriate number of polynomials. Another limitation of
the TEBD algorithm is the handling of long-range interactions with the Suzuki-
Trotter decomposition, which requires the use of advanced techniques such as
swap gates [124]. The TEBD and MPO techniques also face inaccuracies because
of the truncation of the MPS bond dimension, a problem that the Faber
algorithm avoids at the cost of working with a state encompassing the entire
Hilbert space dimension. Of course, using the Faber polynomial limits the
analysis to smaller system sizes. However, it can be advantageous in
situations where MPS calculations fail, such as when the entanglement entropy
scales proportionally to the system size.

## References

  * [1] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ ,  Oxford University Press (Oxford, England, 2002).
  * [2] H. M. Wiseman and G. J. Milburn,  _Quantum Measurement and Control_ ,  Cambridge University Press (Cambridge, England, 2009).
  * [3] G. Lindblad,  _On the generators of quantum dynamical semigroups_ ,  Communications in Mathematical Physics 48(2), 119 (1976),  10.1007/BF01608499.
  * [4] V. Gorini, A. Kossakowski and E. C. G. Sudarshan,  _Completely positive dynamical semigroups of N-level systems_ ,  Journal of Mathematical Physics 17(5), 821 (1976),  10.1063/1.522979,  https://pubs.aip.org/aip/jmp/article-pdf/17/5/821/19090720/821_1_online.pdf.
  * [5] D. Manzano,  _A short introduction to the Lindblad master equation_ ,  AIP Advances 10(2), 025106 (2020),  10.1063/1.5115323,  https://pubs.aip.org/aip/adv/article-pdf/doi/10.1063/1.5115323/12881278/025106_1_online.pdf.
  * [6] L. M. Sieberer, M. Buchhold, J. Marino and S. Diehl,  _Universality in driven open quantum matter_ (2023), 2312.03073.
  * [7] G. T. Landi, D. Poletti and G. Schaller,  _Nonequilibrium boundary-driven quantum systems: Models, methods, and properties_ ,  Rev. Mod. Phys. 94, 045006 (2022),  10.1103/RevModPhys.94.045006.
  * [8] K. Jacobs,  _Quantum Measurement Theory and its Applications_ ,  Cambridge University Press (Cambridge, England, 2014).
  * [9] G. T. Landi, M. J. Kewming, M. T. Mitchison and P. P. Potts,  _Current fluctuations in open quantum systems: Bridging the gap between quantum continuous measurements and full counting statistics_ ,  PRX Quantum 5(2), 020201 (2024),  10.1103/prxquantum.5.020201.
  * [10] H. Carmichael,  _Statistical Methods in Quantum Optics 1_ ,  Springer Science & Business Media (Berlin, Germany, 1999).
  * [11] A. J. Daley,  _Quantum trajectories and open many-body quantum systems_ ,  Advances in Physics 63(2), 77 (2014),  10.1080/00018732.2014.933502.
  * [12] J. Dalibard, Y. Castin and K. Mølmer,  _Wave-function approach to dissipative processes in quantum optics_ ,  Physical Review Letters 68(5), 580 (1992),  10.1103/physrevlett.68.580.
  * [13] M. B. Plenio and P. L. Knight,  _The quantum-jump approach to dissipative dynamics in quantum optics_ ,  Reviews of Modern Physics 70(1), 101 (1998),  10.1103/revmodphys.70.101.
  * [14] C. Gneiting, A. V. Rozhkov and F. Nori,  _Jump-time unraveling of markovian open quantum systems_ ,  Physical Review A 104(6), 062212 (2021),  10.1103/physreva.104.062212.
  * [15] Y. Li, X. Chen and M. P. A. Fisher,  _Quantum zeno effect and the many-body entanglement transition_ ,  Phys. Rev. B 98, 205136 (2018),  10.1103/PhysRevB.98.205136.
  * [16] B. Skinner, J. Ruhman and A. Nahum,  _Measurement-induced phase transitions in the dynamics of entanglement_ ,  Phys. Rev. X 9, 031009 (2019),  10.1103/PhysRevX.9.031009.
  * [17] X. Cao, A. Tilloy and A. De Luca,  _Entanglement in a fermion chain under continuous monitoring_ ,  SciPost Physics 7(2) (2019),  10.21468/scipostphys.7.2.024.
  * [18] Y. Fuji and Y. Ashida,  _Measurement-induced quantum criticality under continuous monitoring_ ,  Phys. Rev. B 102, 054302 (2020),  10.1103/PhysRevB.102.054302.
  * [19] O. Alberton, M. Buchhold and S. Diehl,  _Entanglement transition in a monitored free-fermion chain: From extended criticality to area law_ ,  Physical Review Letters 126(17), 170602 (2021),  10.1103/physrevlett.126.170602.
  * [20] X. Turkeshi, A. Biella, R. Fazio, M. Dalmonte and M. Schiró,  _Measurement-induced entanglement transitions in the quantum ising chain: From infinite to zero clicks_ ,  Phys. Rev. B 103, 224210 (2021),  10.1103/PhysRevB.103.224210.
  * [21] Y. L. Gal, X. Turkeshi and M. Schirò,  _Entanglement dynamics in monitored systems and the role of quantum jumps_ ,  10.48550/ARXIV.2312.13419 (2023).
  * [22] R. El-Ganainy, M. Khajavikhan, D. N. Christodoulides and S. K. Ozdemir,  _The dawn of non-hermitian optics_ ,  Communications Physics 2(1) (2019),  10.1038/s42005-019-0130-z.
  * [23] S. Chandramouli, N. Ossi, Z. H. Musslimani and K. G. Makris,  _Dispersive hydrodynamics in non-hermitian nonlinear schrödinger equation with complex external potential_ ,  Nonlinearity 36(12), 6798 (2023),  10.1088/1361-6544/ad065d.
  * [24] K. Sone and Y. Ashida,  _Anomalous topological active matter_ ,  Physical Review Letters 123(20), 205502 (2019),  10.1103/physrevlett.123.205502.
  * [25] K. Sone, Y. Ashida and T. Sagawa,  _Exceptional non-hermitian topological edge mode and its application to active matter_ ,  Nature Communications 11(1) (2020),  10.1038/s41467-020-19488-0.
  * [26] Y. Ashida, Z. Gong and M. Ueda,  _Non-hermitian physics_ ,  Advances in Physics 69(3), 249 (2020),  10.1080/00018732.2021.1876991.
  * [27] Y. Ashida and M. Ueda,  _Full-counting many-particle dynamics: Nonlocal and chiral propagation of correlations_ ,  Physical Review Letters 120(18), 185301 (2018),  10.1103/physrevlett.120.185301.
  * [28] A. Bácsi and B. Dóra,  _Dynamics of entanglement after exceptional quantum quench_ ,  Physical Review B 103(8), 085137 (2021),  10.1103/physrevb.103.085137.
  * [29] B. Dóra, M. A. Werner and C. u. u. u. u. P. m. c. Moca,  _Quantum quench dynamics in the luttinger liquid phase of the hatano-nelson model_ ,  Physical Review B 108(3), 035104 (2023),  10.1103/physrevb.108.035104.
  * [30] S. Gopalakrishnan and M. J. Gullans,  _Entanglement and purification transitions in non-hermitian quantum mechanics_ ,  Physical Review Letters 126(17), 170503 (2021),  10.1103/physrevlett.126.170503.
  * [31] Y. Le Gal, X. Turkeshi and M. Schirò,  _Volume-to-area law entanglement transition in a non-hermitian free fermionic chain_ ,  SciPost Physics 14(5) (2023),  10.21468/scipostphys.14.5.138.
  * [32] C. Zerba and A. Silva,  _Measurement phase transitions in the no-click limit as quantum phase transitions of a non-hermitean vacuum_ ,  SciPost Physics Core 6(3) (2023),  10.21468/scipostphyscore.6.3.051.
  * [33] X. Turkeshi and M. Schiró,  _Entanglement and correlation spreading in non-hermitian spin chains_ ,  Phys. Rev. B 107, L020403 (2023),  10.1103/PhysRevB.107.L020403.
  * [34] K. Kawabata, T. Numasawa and S. Ryu,  _Entanglement phase transition induced by the non-hermitian skin effect_ ,  Physical Review X 13(2), 021007 (2023),  10.1103/physrevx.13.021007.
  * [35] T. E. Lee,  _Anomalous edge state in a non-hermitian lattice_ ,  Physical Review Letters 116(13), 133903 (2016),  10.1103/physrevlett.116.133903.
  * [36] S. Yao and Z. Wang,  _Edge states and topological invariants of non-hermitian systems_ ,  Phys. Rev. Lett. 121, 086803 (2018),  10.1103/PhysRevLett.121.086803.
  * [37] F. K. Kunst, E. Edvardsson, J. C. Budich and E. J. Bergholtz,  _Biorthogonal bulk-boundary correspondence in non-hermitian systems_ ,  Physical Review Letters 121(2), 026808 (2018),  10.1103/physrevlett.121.026808.
  * [38] X. Zhang, T. Zhang, M.-H. Lu and Y.-F. Chen,  _A review on non-hermitian skin effect_ ,  Advances in Physics: X 7(1) (2022),  10.1080/23746149.2022.2109431.
  * [39] K. Zhang, Z. Yang and C. Fang,  _Universal non-hermitian skin effect in two and higher dimensions_ ,  Nature Communications 13(1) (2022),  10.1038/s41467-022-30161-6.
  * [40] N. Okuma and M. Sato,  _Non-hermitian topological phenomena: A review_ ,  Annual Review of Condensed Matter Physics 14(1), 83 (2023),  10.1146/annurev-conmatphys-040521-033133.
  * [41] J. Claes and T. L. Hughes,  _Skin effect and winding number in disordered non-hermitian systems_ ,  Physical Review B 103(14), l140201 (2021),  10.1103/physrevb.103.l140201.
  * [42] L. Li, C. H. Lee, S. Mu and J. Gong,  _Critical non-hermitian skin effect_ ,  Nature Communications 11(1) (2020),  10.1038/s41467-020-18917-4.
  * [43] K. Kawabata, M. Sato and K. Shiozaki,  _Higher-order non-hermitian skin effect_ ,  Physical Review B 102(20), 205118 (2020),  10.1103/physrevb.102.205118.
  * [44] X. Zhu, H. Wang, S. K. Gupta, H. Zhang, B. Xie, M. Lu and Y. Chen,  _Photonic non-hermitian skin effect and non-bloch bulk-boundary correspondence_ ,  Physical Review Research 2(1), 013280 (2020),  10.1103/physrevresearch.2.013280.
  * [45] Y. Ma and T. L. Hughes,  _Quantum skin hall effect_ ,  Physical Review B 108(10), l100301 (2023),  10.1103/physrevb.108.l100301.
  * [46] H. Geng, J. Y. Wei, M. H. Zou, L. Sheng, W. Chen and D. Y. Xing,  _Nonreciprocal charge and spin transport induced by non-hermitian skin effect in mesoscopic heterojunctions_ ,  Phys. Rev. B 107, 035306 (2023),  10.1103/PhysRevB.107.035306.
  * [47] P. K. Suetin,  _Fundamental properties of faber polynomials_ ,  Russian Mathematical Surveys 19(4), 121 (1964),  10.1070/rm1964v019n04abeh001155.
  * [48] J. H. Curtiss,  _Faber polynomials and the faber series_ ,  The American Mathematical Monthly 78(6), 577 (1971),  10.1080/00029890.1971.11992813.
  * [49] H. Tal-Ezer and R. Kosloff,  _An accurate and efficient scheme for propagating the time dependent schrödinger equation_ ,  The Journal of Chemical Physics 81(9), 3967 (1984),  10.1063/1.448136.
  * [50] A. Suresh, R. D. Soares, P. Mondal, J. P. S. Pires, J. M. V. P. Lopes, A. Ferreira, A. E. Feiguin, P. Plecháč and B. K. Nikolić,  _Electron-mediated entanglement of two distant macroscopic ferromagnets within a nonequilibrium spintronic device_ ,  Physical Review A 109(2), 022414 (2024),  10.1103/physreva.109.022414.
  * [51] J. P. Santos Pires, B. Amorim and J. M. Viana Parente Lopes,  _Landauer transport as a quasisteady state on finite chains under unitary quantum dynamics_ ,  Physical Review B 101(10), 104203 (2020),  10.1103/physrevb.101.104203.
  * [52] J. M. A. Pinho, J. P. S. Pires, S. M. João, B. Amorim and J. M. V. P. Lopes,  _From bloch oscillations to a steady-state current in strongly biased mesoscopic devices_ ,  Physical Review B 108(7), 075402 (2023),  10.1103/physrevb.108.075402.
  * [53] H. P. Veiga, S. M. João, J. M. A. Pinho, J. P. S. Pires and J. M. V. P. Lopes,  _Unambiguous simulation of diffusive charge transport in disordered nanoribbons_ ,  10.48550/ARXIV.2311.03983 (2023).
  * [54] Y. Huang, D. J. Kouri and D. K. Hoffman,  _General, energy-separable faber polynomial representation of operator functions: Theory and application in quantum scattering_ ,  The Journal of Chemical Physics 101(12), 10493 (1994),  10.1063/1.468481.
  * [55] W. Huisinga, L. Pesce, R. Kosloff and P. Saalfrank,  _Faber and newton polynomial integrators for open-system density matrix propagation_ ,  The Journal of Chemical Physics 110(12), 5538 (1999),  10.1063/1.478451.
  * [56] L. Schulz, B. Inci, M. Pech and D. Schulz,  _Subdomain-based exponential integrators for quantum liouville-type equations_ ,  Journal of Computational Electronics 20(6), 2070 (2021),  10.1007/s10825-021-01797-2.
  * [57] A. G. Borisov and S. V. Shabanov,  _Wave packet propagation by the faber polynomial approximation in electrodynamics of passive media_ ,  Journal of Computational Physics 216(1), 391 (2006),  10.1016/j.jcp.2005.12.011.
  * [58] H. Fahs,  _Investigation on polynomial integrators for time-domain electromagnetics using a high-order discontinuous galerkin method_ ,  Applied Mathematical Modelling 36(11), 5466 (2012),  10.1016/j.apm.2011.12.055.
  * [59] R. H. Landau, M. J. Páez and C. C. Bordeianu,  _Computational physics_ ,  Physics textbook. Wiley-VCH, Weinheim, 2., rev. and enl. ed., 1. reprint edn.,  ISBN 9783527406265,  Parallel als digitalisierte Ausg. erschienen (2011).
  * [60] N. Hatano and M. Suzuki,  _Finding Exponential Product Formulas of Higher Orders_ , pp. 37–68,  Springer Berlin Heidelberg,  ISBN 9783540315155,  10.1007/11526216_2 (2005).
  * [61] N. Hatano and D. R. Nelson,  _Localization transitions in non-hermitian quantum mechanics_ ,  Physical Review Letters 77(3), 570 (1996),  10.1103/physrevlett.77.570.
  * [62] N. Hatano and D. R. Nelson,  _Vortex pinning and non-hermitian quantum mechanics_ ,  Physical Review B 56(14), 8651 (1997),  10.1103/physrevb.56.8651.
  * [63] M. Fruchart, R. Hanai, P. B. Littlewood and V. Vitelli,  _Non-reciprocal phase transitions_ ,  Nature 592(7854), 363 (2021),  10.1038/s41586-021-03375-9.
  * [64] A. Weiße, G. Wellein, A. Alvermann and H. Fehske,  _The kernel polynomial method_ ,  Reviews of Modern Physics 78(1), 275 (2006),  10.1103/revmodphys.78.275.
  * [65] A. Braun and P. Schmitteckert,  _Numerical evaluation of green’s functions based on the chebyshev expansion_ ,  Physical Review B 90(16), 165112 (2014),  10.1103/physrevb.90.165112.
  * [66] A. Ferreira and E. R. Mucciolo,  _Critical delocalization of chiral zero energy modes in graphene_ ,  Phys. Rev. Lett. 115, 106601 (2015),  10.1103/PhysRevLett.115.106601.
  * [67] S. M. João and J. M. Viana Parente Lopes,  _Basis-independent spectral methods for non-linear optical response in arbitrary tight-binding models_ ,  Journal of Physics: Condensed Matter 32(12), 125901 (2019),  10.1088/1361-648x/ab59ec.
  * [68] S. a. M. João, M. Anđelkovió, L. Covaci, T. G. Rappoport, J. a. M. V. P. Lopes and A. Ferreira,  _Kite: high-performance accurate modelling of electronic structure and response functions of large molecules, disordered crystals and heterostructures_ ,  Royal Society Open Science 7(2), 191809 (2020),  10.1098/rsos.191809.
  * [69] S. M. João, J. M. Viana Parente Lopes and A. Ferreira,  _High-resolution real-space evaluation of the self-energy operator of disordered lattices: Gade singularity, spin–orbit effects and p-wave superconductivity_ ,  Journal of Physics: Materials 5(4), 045002 (2022),  10.1088/2515-7639/ac91f9.
  * [70] G. Faber,  _Ueber polynomische entwicklungen_ ,  Mathematische Annalen 57, 389 (1903).
  * [71] G. Faber,  _Uber polynomische entwicklungen ii_ ,  Mathematische Annalen 64(1), 116 (1907),  10.1007/bf01449884.
  * [72] A. W. Sandvik, A. Avella and F. Mancini,  _Computational studies of quantum spin systems_ ,  In _AIP Conference Proceedings_. AIP,  10.1063/1.3518900 (2010).
  * [73] S. W. Ellacott,  _Computation of faber series with application to numerical polynomial approximation in the complex plane_ ,  Mathematics of Computation 40(162), 575 (1983),  10.1090/s0025-5718-1983-0689474-7.
  * [74] I. S. Gradštejn, I. M. Ryzhik, D. Zwillinger and V. H. Moll, eds.,  _Table of integrals, series, and products_ ,  Academic Press, Waltham, MA, eighth edition edn.,  ISBN 0123849349 (2015).
  * [75] S. Bravyi,  _Lagrangian representation for fermionic linear optics_ ,  Quantum Info. Comput. 5(3), 216–238 (2005).
  * [76] N. Laflorencie,  _Quantum entanglement in condensed matter systems_ ,  Physics Reports 646, 1 (2016),  10.1016/j.physrep.2016.06.008.
  * [77] L. Amico, R. Fazio, A. Osterloh and V. Vedral,  _Entanglement in many-body systems_ ,  Rev. Mod. Phys. 80, 517 (2008),  10.1103/RevModPhys.80.517.
  * [78] I. Peschel,  _Calculation of reduced density matrices from correlation functions_ ,  Journal of Physics A: Mathematical and General 36(14), L205 (2003),  10.1088/0305-4470/36/14/101.
  * [79] S.-A. Cheong and C. L. Henley,  _Many-body density matrices for free fermions_ ,  Physical Review B 69(7), 075111 (2004),  10.1103/physrevb.69.075111.
  * [80] I. Peschel and V. Eisler,  _Reduced density matrices and entanglement entropy in free lattice models_ ,  Journal of Physics A: Mathematical and Theoretical 42(50), 504003 (2009),  10.1088/1751-8113/42/50/504003.
  * [81] X. Turkeshi, L. Piroli and M. Schirò,  _Density and current statistics in boundary-driven monitored fermionic chains_ ,  Phys. Rev. B 109, 144306 (2024),  10.1103/PhysRevB.109.144306.
  * [82] P. Jordan and E. Wigner,  _Über das paulische Äquivalenzverbot_ ,  Zeitschrift für Physik 47(9-10), 631 (1928),  10.1007/bf01331938.
  * [83] S. Sachdev,  _Quantum Phase Transitions_ ,  Cambridge University Press,  ISBN 9780521582544 (1999).
  * [84] T. Antal, Z. Rácz, A. Rákos and G. M. Schütz,  _Transport in the $\mathrm{XX}$ chain at zero temperature: Emergence of flat magnetization profiles_,  Phys. Rev. E 59, 4912 (1999),  10.1103/PhysRevE.59.4912.
  * [85] J. Lancaster, E. Gull and A. Mitra,  _Quenched dynamics in interacting one-dimensional systems: Appearance of current-carrying steady states from initial domain wall density profiles_ ,  Phys. Rev. B 82, 235124 (2010),  10.1103/PhysRevB.82.235124.
  * [86] G. Misguich, K. Mallick and P. L. Krapivsky,  _Dynamics of the spin- 12 heisenberg chain initialized in a domain-wall state_ ,  Physical Review B 96(19), 195151 (2017),  10.1103/physrevb.96.195151.
  * [87] V. Eisler and F. Maislinger,  _Hydrodynamical phase transition for domain-wall melting in the xy chain_ ,  Phys. Rev. B 98, 161117 (2018),  10.1103/PhysRevB.98.161117.
  * [88] B. Bertini, M. Collura, J. De Nardis and M. Fagotti,  _Transport in out-of-equilibriumxxzchains: Exact profiles of charges and currents_ ,  Physical Review Letters 117(20), 207201 (2016),  10.1103/physrevlett.117.207201.
  * [89] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura,  _Emergent hydrodynamics in integrable quantum systems out of equilibrium_ ,  Physical Review X 6(4), 041065 (2016),  10.1103/physrevx.6.041065.
  * [90] F. H. Essler,  _A short introduction to generalized hydrodynamics_ ,  Physica A: Statistical Mechanics and its Applications 631, 127572 (2023),  10.1016/j.physa.2022.127572.
  * [91] P. Ruggiero, P. Calabrese, B. Doyon and J. Dubail,  _Quantum generalized hydrodynamics_ ,  Physical Review Letters 124(14), 140603 (2020),  10.1103/physrevlett.124.140603.
  * [92] T. Orito and K.-I. Imura,  _Unusual wave-packet spreading and entanglement dynamics in non-hermitian disordered many-body systems_ ,  Physical Review B 105(2), 024303 (2022),  10.1103/physrevb.105.024303.
  * [93] H. Spring, V. Könye, F. A. Gerritsma, I. C. Fulga and A. R. Akhmerov,  _Phase transitions of wave packet dynamics in disordered non-hermitian systems_ ,  SciPost Physics 16(5) (2024),  10.21468/scipostphys.16.5.120.
  * [94] B. Barch,  _Locality, correlations, information, and non-hermitian quantum systems_ ,  10.48550/ARXIV.2405.16842 (2024).
  * [95] V. Alba and F. Carollo,  _Spreading of correlations in markovian open quantum systems_ ,  Physical Review B 103(2), l020302 (2021),  10.1103/physrevb.103.l020302.
  * [96] F. Carollo and V. Alba,  _Dissipative quasiparticle picture for quadratic markovian open quantum systems_ ,  Physical Review B 105(14), 144305 (2022),  10.1103/physrevb.105.144305.
  * [97] F. Rottoli, S. Scopa and P. Calabrese,  _Entanglement hamiltonian during a domain wall melting in the free fermi chain_ ,  Journal of Statistical Mechanics: Theory and Experiment 2022(6), 063103 (2022),  10.1088/1742-5468/ac72a1.
  * [98] B.-Q. Jin and V. E. Korepin,  _Quantum spin chain, toeplitz determinants and the fisher–hartwig conjecture_ ,  Journal of Statistical Physics 116(1–4), 79 (2004),  10.1023/b:joss.0000037230.37166.42.
  * [99] P. Calabrese and J. Cardy,  _Time dependence of correlation functions following a quantum quench_ ,  Phys. Rev. Lett. 96, 136801 (2006),  10.1103/PhysRevLett.96.136801.
  * [100] P. Calabrese and J. Cardy,  _Entanglement entropy and conformal field theory_ ,  Journal of Physics A: Mathematical and Theoretical 42(50), 504005 (2009),  10.1088/1751-8113/42/50/504005.
  * [101] J.-M. Stéphan and J. Dubail,  _Local quantum quenches in critical one-dimensional systems: entanglement, the loschmidt echo, and light-cone effects_ ,  Journal of Statistical Mechanics: Theory and Experiment 2011(08), P08019 (2011),  10.1088/1742-5468/2011/08/p08019.
  * [102] S. Mu, C. H. Lee, L. Li and J. Gong,  _Emergent fermi surface in a many-body non-hermitian fermionic chain_ ,  Physical Review B 102(8), 081115 (2020),  10.1103/physrevb.102.081115.
  * [103] J. Sirker,  _Transport in one-dimensional integrable quantum systems_ ,  SciPost Physics Lecture Notes (2020),  10.21468/scipostphyslectnotes.17.
  * [104] V. B. Bulchandani, S. Gopalakrishnan and E. Ilievski,  _Superdiffusion in spin chains_ ,  Journal of Statistical Mechanics: Theory and Experiment 2021(8), 084001 (2021),  10.1088/1742-5468/ac12c7.
  * [105] T. Haga, M. Nakagawa, R. Hamazaki and M. Ueda,  _Liouvillian skin effect: Slowing down of relaxation processes without gap closing_ ,  Physical Review Letters 127(7), 070402 (2021),  10.1103/physrevlett.127.070402.
  * [106] F. Yang, Q.-D. Jiang and E. J. Bergholtz,  _Liouvillian skin effect in an exactly solvable model_ ,  Physical Review Research 4(2), 023160 (2022),  10.1103/physrevresearch.4.023160.
  * [107] Z. Wang, Y. Lu, Y. Peng, R. Qi, Y. Wang and J. Jie,  _Accelerating relaxation dynamics in open quantum systems with liouvillian skin effect_ ,  Physical Review B 108(5), 054313 (2023),  10.1103/physrevb.108.054313.
  * [108] X. Li, M. A. Begaowe, S. Zhang and B. Flebus,  _Reciprocal reservoir induced non-hermitian skin effect_ ,  10.48550/ARXIV.2307.15792 (2023).
  * [109] S. E. Begg and R. Hanai,  _Quantum criticality in open quantum spin chains with nonreciprocity_ ,  Physical Review Letters 132(12), 120401 (2024),  10.1103/physrevlett.132.120401.
  * [110] P. Calabrese and J. Cardy,  _Evolution of entanglement entropy in one-dimensional systems_ ,  Journal of Statistical Mechanics: Theory and Experiment 2005(04), P04010 (2005),  10.1088/1742-5468/2005/04/p04010.
  * [111] M. Fagotti and P. Calabrese,  _Evolution of entanglement entropy following a quantum quench: Analytic results for thexychain in a transverse magnetic field_ ,  Physical Review A 78(1), 010306 (2008),  10.1103/physreva.78.010306.
  * [112] M. B. Hastings,  _An area law for one-dimensional quantum systems_ ,  Journal of Statistical Mechanics: Theory and Experiment 2007(08), P08024 (2007),  10.1088/1742-5468/2007/08/P08024.
  * [113] N. Hatano and J. Feinberg,  _Chebyshev-polynomial expansion of the localization length of hermitian and non-hermitian random chains_ ,  Physical Review E 94(6), 063305 (2016),  10.1103/physreve.94.063305.
  * [114] G. Chen, F. Song and J. L. Lado,  _Topological spin excitations in non-hermitian spin chains with a generalized kernel polynomial algorithm_ ,  Physical Review Letters 130(10), 100401 (2023),  10.1103/physrevlett.130.100401.
  * [115] A. Holzner, A. Weichselbaum, I. P. McCulloch, U. Schollwöck and J. von Delft,  _Chebyshev matrix product state approach for spectral functions_ ,  Physical Review B 83(19), 195115 (2011),  10.1103/physrevb.83.195115.
  * [116] F. A. Wolf, J. A. Justiniano, I. P. McCulloch and U. Schollwöck,  _Spectral functions and time evolution from the chebyshev recursion_ ,  Physical Review B 91(11), 115144 (2015),  10.1103/physrevb.91.115144.
  * [117] J. C. Halimeh, F. Kolley and I. P. McCulloch,  _Chebyshev matrix product state approach for time evolution_ ,  Physical Review B 92(11), 115130 (2015),  10.1103/physrevb.92.115130.
  * [118] H. D. Xie, R. Z. Huang, X. J. Han, X. Yan, H. H. Zhao, Z. Y. Xie, H. J. Liao and T. Xiang,  _Reorthonormalization of chebyshev matrix product states for dynamical correlation functions_ ,  Physical Review B 97(7), 075111 (2018),  10.1103/physrevb.97.075111.
  * [119] J. Bak and D. J. Newman,  _Complex Analysis_ ,  Springer New York,  ISBN 9781441972880,  10.1007/978-1-4419-7288-0 (2010).
  * [120] G. Vidal,  _Efficient simulation of one-dimensional quantum many-body systems_ ,  Physical Review Letters 93(4), 040502 (2004),  10.1103/physrevlett.93.040502.
  * [121] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck and C. Hubig,  _Time-evolution methods for matrix-product states_ ,  Annals of Physics 411, 167998 (2019),  https://doi.org/10.1016/j.aop.2019.167998.
  * [122] M. Fishman, S. White and E. Stoudenmire,  _The itensor software library for tensor network calculations_ ,  SciPost Physics Codebases (2022),  10.21468/scipostphyscodeb.4.
  * [123] M. Fishman, S. White and E. Stoudenmire,  _Codebase release 0.3 for itensor_ ,  SciPost Physics Codebases (2022),  10.21468/scipostphyscodeb.4-r0.3.
  * [124] E. M. Stoudenmire and S. R. White,  _Minimally entangled typical thermal state algorithms_ ,  New Journal of Physics 12(5), 055026 (2010),  10.1088/1367-2630/12/5/055026.