# Breakdown of Linear Spin-Wave Theory in a Non-Hermitian Quantum Spin Chain

Julien Despres  Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
Leonardo Mazza  Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
Marco Schirò  JEIP, USR 3573 CNRS, Collège de France, PSL Research University,
11 Place Marcelin Berthelot, 75321 Paris Cedex 05, France

(August 28, 2024)

###### Abstract

We present the spin-wave theory of the excitation spectrum and quench dynamics
of the non-Hermitian transverse-field Ising model. The complex excitation
spectrum is obtained for a generic hypercubic lattice using the linear
approximation of the Holstein-Primakoff transformation together with the
complex bosonic Bogolyubov transformation. In the one-dimensional case, our
result compares very well with the exact quasiparticle dispersion relation
obtained via a fermionic representation of the problem, at least in the regime
of large dissipation and transverse field. When applied to the quench dynamics
we show however that the linear spin-wave approximation breaks down and the
bosonic theory is plagued by a divergence at finite times. We understand the
origin of this instability using a single mode approximation. While limited to
short times, we show that this approach allows us to characterize the dynamics
arising from the quench of the dissipative term and the structure of the Lieb-
Robinson light-cone of the propagation quantum correlations. Furthermore, for
the one-dimensional case, the linear spin-wave dynamics shows good agreement
with the exact fermionic solution, both for the local magnetization and the
spin-spin correlations.

##  I Introduction

In the last decades, the significant and major advances in the experimental
control of quantum matter have given dramatic impulse to the experimental and
theoretical investigation of the out-of-equilibrium dynamics of isolated
quantum lattice models [1; 2; 3; 4; 5; 6]. However, since a quantum system
cannot be perfectly isolated from its environment, the study of open quantum
systems has attracted a lot of interest in the last few years by taking into
account dissipative effects such as local dephasing noise [7; 8; 9], local
gain/loss processes [10; 11; 12; 13] and incoherent hopping [14]. These
quantum lattice models are governed by the Lindblad master equation [15]
describing the time evolution of the density matrix associated to the quantum
system. The latter contains a unitary part representing the standard time
evolution of the quantum system with its corresponding Hermitian Hamiltonian,
but also a non-unitary part describing the dissipation which is represented by
the so-called jump operators. In the no-click limit i.e. when no quantum jump
occurs, the Lindblad master equation can be rewritten as the time evolution of
a quantum system governed by a non-Hermitian Hamiltonian [16]. Its non-
Hermiticity is responsible for the non-unitary time evolution of the open
many-body quantum system [17].

Non-Hermitian quantum many-body systems host unexpected novel and rich
physical properties compared to their Hermitian counterparts. Their spectral
features have been the focus of several studies related to the existence of
exceptional points [18] and their associated topological properties [19; 20].
More recently their entanglement structure has drawn attention, both for what
concerns eigenstates in the biorthogonal formulation [21; 22; 23] and for
dynamics, where several entanglement transitions have been reported [24; 25;
26; 27; 28; 29; 30]. Similarly, the non-unitary dynamics after a quantum
quench or a Floquet driving has also been studied [31; 32]. While most of the
literature has focused on many-body _non-interacting systems_ , less is known
about the behavior of interacting non-Hermitian quantum many body systems. One
dimensional models, such as non-Hermitian spin chains or the interacting
version of Hatano-Nelson model [33], have been studied using exact
diagonalization [34; 35; 36] and extensions of Bethe Ansatz and bosonization
[37; 38; 39].

In this respect it is important to develop theoretical methods for non-
Hermitian systems, which are able to capture their spectral and dynamical
properties in particular in dimensions higher than one. For quantum spin
models a textbook approach is offered by spin-wave theory within the Holstein-
Primakoff transformation [40; 41], that provides a way to characterize the
quasiparticle spectrum of various interacting systems in the equilibrium low-
energy limit and which has been successfully extended to the unitary dynamics
[42; 43]. The goal of this work is to develop a spin-wave theory for a
prototypical non-Hermitian quantum spin chain, the Ising model in a complex
transverse field, in one, two and three dimensional lattices. We find that
while a straightforward generalization of linear spin-wave theory captures
well the spectrum of non-Hermitian quasiparticles, at least for strong values
of dissipation and transverse field, it fails dramatically in describing the
non-unitary dynamics at long times, particularly in cases where the initial
state and the long-time steady state are not smoothly connected by a small
number of quasiparticle excitations (spin flips). While one might have
expected in this case the non-Hermitian dynamics to lead to a growth
(amplification) of quasiparticle occupation numbers, thus quickly leaving the
regime of validity of linear spin wave theory, we show in fact that the
breakdown of linear spin-wave theory is much more severe and leads to a
finite-time divergence of the dynamics. We provide a simple analytical
understanding of this instability in terms of a single mode bosonic model.
Although limited to short-times, we show that our spin-wave analysis allows us
to discuss the dynamics of correlation spreading in the quantum Ising chain.
This topic has been extensively investigated for isolated quantum lattice
models [44; 45; 46; 47; 48; 49; 50; 51; 52; 53], but is still debated in the
specific case of non-Hermitian quantum lattice models and a large number of
points remain unanswered. For one dimension we benchmark this approach with
the exact solution obtained via Jordan-Wigner transformation into non-
interacting non-Hermitian fermions [54; 25], confirming the validity of the
proposed approach for the spectrum and for the short-time dynamics.

This article is organized as follows: in Sec. II we start by introducing the
non-Hermitian Ising model and we apply the linear spin-wave theory to deduce
its excitation spectrum. Then, in Sec. III we investigate the quench dynamics
of the resulting bosonic non-Hermitian model using the equations of motion for
the single particle correlation matrix. We highlight a dynamical instability
at finite time in the bosonic theory that we understand using a single mode
model. Using this theory we discuss the short-time spreading of spin-spin
correlation functions and the local magnetization in one and two dimensional
lattices, and compare with the exact fermionic solution in the former case.
Finally in Sec. IV we present our conclusions. In the Appendices we include
further technical details on our results.

##  II Model and spin-wave theory for the excitation spectrum

###  II.1 One-dimensional lattice

We focus on the spin-$1/2$ non-Hermitian transverse-field Ising model (TFIM)
on a one-dimensional chain of $N_{s}$ sites whose lattice spacing is fixed to
unity ($a=1$); for simplicity, we also set $\hbar=1$. The Hamiltonian
$\hat{H}$ reads:

$\hat{H}=J\sum_{R}\hat{S}^{x}_{R}\hat{S}^{x}_{R+1}-(h+i\gamma)\sum_{R}\hat{S}^{z}_{R},$
(1)

where $\hat{S}^{\alpha}_{R}$ refers to the spin operator acting on the lattice
site $R$ along the $\alpha\in\\{x,z\\}$ axis. $J>0$ denotes the
antiferromagnetic coupling between nearest neighbours spins, $h>0$ is the
constant and homogeneous transverse magnetic field and $\gamma$ is the
dissipative strength representing the local dephasing noise. From the
perspective of open quantum systems described by a Lindblad master equation
[15] the latter can be interpreted as a heating process arising from a jump
operators $\hat{L}_{R}=\sqrt{2\gamma}(1/2+\hat{S}^{z}_{R})$ when $\gamma>0$,
leading to a non-Hermitian Hamiltonian of the form in Eq. (1). In most of the
paper we will take $\gamma>0$, but in Sec. III.5 we will present a set of data
obtained for $\gamma<0$, which corresponds to a different dissipative process
with jump $\hat{L}_{R}=\sqrt{2|\gamma|}(1/2-\hat{S}^{z}_{R})$. In what
follows, we will mainly discuss the properties of the paramagnetic phase with
spins aligned with the magnetic field for $h\gg J$. In most of our study we
will also take $|\gamma|\gg J$. For $J=0$ the properties of this non-Hermitian
Hamiltonian are easy to discuss because the sites are decoupled. Let us focus
on a single site for simplicity. The $\ket{\uparrow}$ and $\ket{\downarrow}$
states are eigenvectors of $\hat{H}$ with eigenvalues
$\mp\frac{1}{2}(h+i\gamma)$. If we take $\ket{\uparrow}$ as a reference state,
we can thus introduce a complex excitation spectrum equal to $(h+i\gamma)$.
The goal of this section is to discuss how this is modified by the presence of
the antiferromagnetic coupling $J$.

In order to deduce the excitation spectrum of the non-Hermitian TFIM, we rely
on the Holstein-Primakoff transformation:

$\displaystyle\hat{S}^{+}_{R}=\left(1-\hat{a}^{{\dagger}}_{R}\hat{a}_{R}\right)^{1/2}\hat{a}_{R},$
(2a)
$\displaystyle\hat{S}^{-}_{R}=\hat{a}^{{\dagger}}_{R}\left(1-\hat{a}^{{\dagger}}_{R}\hat{a}_{R}\right)^{1/2},$
(2b) $\displaystyle\hat{S}^{z}_{R}=1/2-\hat{a}^{{\dagger}}_{R}\hat{a}_{R};$
(2c)

where the $\hat{a}_{R}$ ($\hat{a}^{{\dagger}}_{R}$) are annihilation
(creation) bosonic operators at the lattice site $R$ that obey canonical
commutation relations. While this is an exact mapping between spin $1/2$ and
bosons, the presence of the square root makes the transformation non-linear
and of limited practical use since the Hamiltonian of the TFIM written in the
bosonic language is highly non-linear. However, in the regime in which the
bosonic occupation number remains small,
$\langle\hat{a}^{{\dagger}}_{R}\hat{a}_{R}\rangle\ll 1$, which is the case for
example in the paramagnetic phase of the TFIM in equilibrium, one can expand
the square root entering the Holstein-Primakoff transformation and write to
leading order $\hat{S}^{+}_{R}=\hat{a}_{R}$ and
$\hat{S}^{-}_{R}=\hat{a}^{{\dagger}}_{R}$ so that

$\hat{S}^{x}_{R}=\frac{\hat{a}_{R}+\hat{a}^{{\dagger}}_{R}}{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{S}^{y}_{R}=\frac{\hat{a}_{R}-\hat{a}^{{\dagger}}_{R}}{2i},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{S}^{z}_{R}=\frac{1}{2}-\hat{a}^{{\dagger}}_{R}\hat{a}_{R}.$ (3)

This is the essence of the linear spin wave theory, which leads to a quadratic
bosonic Hamiltonian that is best expressed in momentum space (see Appendix B):

$\hat{H}=\frac{1}{2}\sum_{k}A_{k}\left(\hat{a}^{{\dagger}}_{k}\hat{a}_{k}+\hat{a}_{-k}\hat{a}^{{\dagger}}_{-k}\right)+B_{k}\left(\hat{a}^{{\dagger}}_{k}\hat{a}^{{\dagger}}_{-k}+\hat{a}_{k}\hat{a}_{-k}\right),$
(4)

with $A_{k}=h+(J/2)\cos(k)+i\gamma$ and $B_{k}=(J/2)\cos(k)$. The momentum-
dependent complex function $A_{k}$ is responsible for the non-hermiticity of
the Hamiltonian. Hamiltonian (4) is the spin-wave theory of the non-Hermitian
TFIM.

Figure 1: Analytical results concerning the complex excitation spectrum
$E_{k}$ of the non-Hermitian transverse Ising chain in the $z$ polarized phase
using the bosonic approach relying on the linear approximation to the
Holstein-Primakoff transformation, together with the complex bosonic
Bogolyubov transformation. The real and imaginary part of the excitation
spectrum $E_{k}$ denoted by $\Re(E_{k})$ and $\Im(E_{k})$, deduced from Eq.
(8), are represented as dashed magenta lines on the panels (a,c,e,g) and
(b,d,f,h) respectively. The excitation spectrum $E_{k}$ is compared to the
momentum-dependent function $A_{k}$ given by $A_{k}=h+(J/2)\cos(k)+i\gamma$,
see dashed-dotted blue lines on panels (a,c,e,g) and (b,d,f,h) for the real
and imaginary part of $A_{k}$ respectively. The complex quasiparticle
dispersion relation $E_{k}$ is also compared to the complex excitation
spectrum $\epsilon_{k}$ found using the fermionic approach relying on the
Jordan-Wigner transformation together with the complex fermionic Bogolyubov
theory. The real and imaginary part of the excitation spectrum $\epsilon_{k}$,
deduced from Eq. (15), are represented as solid green lines on the panels
(a,c,e,g) and (b,d,f,h) respectively. The different set of parameters are
$J=1$, $h=20$, $\gamma=10$ for the panels (a,b), $J=1$, $h=2$, $\gamma=1$ for
the panels (c,d), $J=1$, $h=0.5$, $\gamma=1$ for the panels (e,f) and $J=1$,
$h=0.1$, $\gamma=1$ for the panels (g,h).

A complex bosonic Bogolyubov transformation [55] diagonalizes the non-
Hermitian Hamiltonian at Eq. (4), which reads:

$\displaystyle\hat{a}_{k}^{{\dagger}}=$
$\displaystyle\bar{\hat{b}}_{k}\cosh\left(\frac{\theta_{k}}{2}\right)+\hat{b}_{-k}\sinh\left(\frac{\theta_{k}}{2}\right),$
(5a) $\displaystyle\hat{a}_{k}=$
$\displaystyle\bar{\hat{b}}_{-k}\sinh\left(\frac{\theta_{k}}{2}\right)+\hat{b}_{k}\cosh\left(\frac{\theta_{k}}{2}\right),$
(5b)

where $\theta_{k}\in\mathbb{C}$ is an even function of $k$ and $\hat{b}_{k}$
($\bar{\hat{b}}_{k}$) are annihilation (creation) bosonic operators for a
quasimomentum $k$ satisfying:

$[\hat{b}_{k},\bar{\hat{b}}_{k^{\prime}}]=\delta_{k,k^{\prime}},\quad[\hat{b}_{k},\hat{b}_{k^{\prime}}]=[\bar{\hat{b}}_{k},\bar{\hat{b}}_{k^{\prime}}]=0.$
(6)

They define a vacuum state $\ket{0}$ via $\hat{b}_{k}\ket{0}=0$ that, in the
Hermitian case, corresponds to the ground state; here, it is rather a
reference state for the bosonic expansion. Note that
$\bar{\hat{b}}_{k}\neq\hat{b}_{k}^{{\dagger}}$ because
$\theta_{k}=\mathrm{arctanh}\left(-B_{k}/A_{k}\right)\in\mathbb{C}$. The
diagonal form of $\hat{H}$ reads:

$\hat{H}=\frac{1}{2}\sum_{k}E_{k}\left(\bar{\hat{b}}_{k}\hat{b}_{k}+\hat{b}_{-k}\bar{\hat{b}}_{-k}\right),$
(7)

where the complex excitation spectrum $E_{k}$ is given by

$E_{k}=\mathrm{sgn}(\Re(A_{k}))\sqrt{A_{k}^{2}-B_{k}^{2}},$ (8)

and is plotted on Fig. 1.

When $J\ll h,\gamma$, the limit of our interest, it is possible to perform a
Taylor expansion of the spectrum:

$E_{k}\approx(h+i\gamma)+\frac{J}{2}\cos(k)-\frac{J^{2}}{8(h+i\gamma)}\cos^{2}(k).$
(9)

Fig. 1 highlights the fact that in general even an approximation at first
order in $J$ gives a good qualitative agreement for the real part of the
excitation spectrum, $\Re(E_{k})\approx h+\frac{J}{2}\cos(k)$. For $h>J/2$ the
spectrum is gapped whereas for $h\leq J/2$ it becomes gapless. In order to
describe the momentum dependence of the imaginary part of the excitation
spectrum, it is instead necessary to go at second order in $J$, since the
first order correction is real. We obtain
$\Im(E_{k})\approx\gamma\left(1+\frac{J^{2}}{8(h^{2}+\gamma^{2})}\cos^{2}(k)\right)$
and we observe that we obtain no stable modes in the limit $h,\gamma\gg J$.

Remarkably, the first order approximation yields a good approximation of the
data even when $J$ is of the order of $h$ and $\gamma$. Even if it is tempting
to use the bosonic model in this regime of parameters, one should always
remember that the approximation based on the linearization of the Holstein-
Primakoff transformation is valid only for $h,\gamma\gg J$.

###  II.2 One-dimensional lattice: Comparison with the exact fermionic
approach.

In order to assess the validity of the spin-wave theory developed in the
previous paragraph, we benchmark it with the exact results that can be
obtained in the one-dimensional case via the Jordan-Wigner transformation and
the mapping to fermions. In this case, the excitation spectrum of the non-
Hermitian TFIM can be obtained via a complex fermionic Bogolyubov theory.

We introduce the fermionic creation (annihilation) operators
$\hat{c}^{{\dagger}}_{R}$ $\left(\hat{c}_{R}\right)$ for the lattice site $R$
that satisfy the canonical anti-commutation relations. They can be related to
spin-1/2 operators via the Jordan-Wigner transformation:

$\displaystyle\hat{S}^{x}_{R}=\frac{1}{2}\hat{K}_{R}\left(\hat{c}^{{\dagger}}_{R}+\hat{c}_{R}\right),$
(10a)
$\displaystyle\hat{S}^{y}_{R}=\frac{i}{2}\hat{K}_{R}\left(\hat{c}^{{\dagger}}_{R}-\hat{c}_{R}\right),$
(10b)
$\displaystyle\hat{S}^{z}_{R}=\frac{1}{2}-\hat{c}^{{\dagger}}_{R}\hat{c}_{R},$
(10c)

with the non-local string operator
$\hat{K}_{R}=\prod_{R^{\prime}=0}^{R-1}\left(1-2\hat{c}^{{\dagger}}_{R^{\prime}}\hat{c}_{R^{\prime}}\right)$.

In this novel language, the one-dimensional non-Hermitian TFIM in Eq. (1) is a
quadratic function of the fermionic operators (here in momentum
representation, see Appendix B):

$\hat{H}=\sum_{k}\mathcal{A}_{k}\left(\hat{c}^{{\dagger}}_{k}\hat{c}_{k}-\hat{c}_{-k}\hat{c}^{{\dagger}}_{-k}\right)+\mathcal{B}_{k}\left(\hat{c}^{{\dagger}}_{k}\hat{c}^{{\dagger}}_{-k}+\hat{c}_{k}\hat{c}_{-k}\right),$
(11)

where the momentum-dependent functions $\mathcal{A}_{k}$ and
$\mathcal{B}_{k}$, responsible for the non-hermiticity of $\hat{H}$, are
defined as follows [56; 54; 25]: $\mathcal{A}_{k}=(J/4)\cos(k)+(h+i\gamma)/2$
and $\mathcal{B}_{k}=-i(J/4)\sin(k)$. We rely on a complex fermionic
Bogolyubov transformation [57] to diagonalize the quadratic Hamiltonian
$\hat{H}$ at Eq. (11) defined by

$\displaystyle\hat{c}^{{\dagger}}_{-k}=\cos\left(\frac{\theta_{k}}{2}\right)\bar{\hat{\eta}}_{-k}+i\sin\left(\frac{\theta_{k}}{2}\right)\hat{\eta}_{k},$
(12)
$\displaystyle\hat{c}_{k}=\cos\left(\frac{\theta_{k}}{2}\right)\hat{\eta}_{k}+i\sin\left(\frac{\theta_{k}}{2}\right)\bar{\hat{\eta}}_{-k},$
(13)

where $\theta_{k}\in\mathbb{C}$ is an odd function of $k$. $\hat{\eta}_{k}$
$\left(\bar{\hat{\eta}}_{k}\right)$ denotes the annihilation (creation)
fermionic Bogolyubov operator for a quasimomentum $k$. They fulfill the
canonical anti-commutation rules for fermions given by
$\\{\hat{\eta}_{k},\bar{\hat{\eta}}_{k^{\prime}}\\}=\delta_{k,k^{\prime}}$,
$\\{\hat{\eta}_{k},\hat{\eta}_{k^{\prime}}\\}=\\{\bar{\hat{\eta}}_{k},\bar{\hat{\eta}}_{k^{\prime}}\\}=0$.
Besides, they define the vacuum state $\hat{\eta}_{k}\ket{0}=0$. It is crucial
to underline that $\bar{\hat{\eta}}_{k}\neq\hat{\eta}^{{\dagger}}_{k}$ because
$\theta_{k}$ is complex.

After applying the complex fermionic Bogolyubov transformation, we find
$\theta_{k}=\mathrm{arctan}\left(i\mathcal{B}_{k}/\mathcal{A}_{k}\right)$, see
a plot in Appendix A where we make a comparison with the bosonic case, and the
Hamiltonian $\hat{H}$ becomes diagonal in terms of $\hat{\eta}_{k}$,
$\bar{\hat{\eta}}_{k}$ and has the following expression

$\hat{H}=\sum_{k}\epsilon_{k}\bar{\hat{\eta}}_{k}\hat{\eta}_{k},$ (14)

where the complex excitation spectrum $\epsilon_{k}$ has the following
expression

$\epsilon_{k}=2\mathrm{sgn}(\Re(\mathcal{A}_{k}))\sqrt{\mathcal{A}_{k}^{2}-\mathcal{B}_{k}^{2}}.$
(15)

In the paramagnetic phase $J\ll h,\gamma$, it can be expanded as:

$\epsilon_{k}\approx(h+i\gamma)+\frac{J}{2}\cos(k)+\frac{J^{2}}{8(h+i\gamma)}\sin^{2}(k).$
(16)

Hence, we get
$\Im(\epsilon_{k})\approx\gamma\left(1-\frac{J^{2}}{8(h^{2}+\gamma^{2})}\sin^{2}(k)\right)$
and we obtain no stable modes in the limit $h,\gamma\gg J$. We observe that
the results deduced from the bosonic spin wave theory agree at first order in
$J$ with the exact fermionic spectrum. This means that at the level of the
real part of the excitation spectrum the two theories are basically identical.
It is interesting to observe the qualitative resemblance of the imaginary part
of the spectrum in the fermionic and bosonic cases. The momentum dependence,
captured in both cases by the term proportional to $J^{2}$, in one case is
proportional to $+\cos^{2}(k)$ and in the other to $-\sin^{2}(k)$. Although
this leads to a quantitative disagreement, the qualitative momentum dependence
is the same. Fig. 1 presents a quantitative comparison of the bosonic spectrum
with the exact fermionic one.

###  II.3 Generalization to a square and hypercubic lattice

Having benchmarked the spin-wave theory for the non-Hermitian TFIM with the
exact result, we can now generalised it to higher-dimensional lattices, where
exact approaches are not known. We begin for simplicity with the case of a
square lattice:

$\hat{H}=J\sum_{\mathbb{R}}\left(\hat{S}^{x}_{\mathbb{R}}\hat{S}^{x}_{\mathbb{R}+\mathbb{x}}+\hat{S}^{x}_{\mathbb{R}}\hat{S}^{x}_{\mathbb{R}+\mathbb{y}}\right)-(h+i\gamma)\sum_{\mathbb{R}}\hat{S}^{z}_{\mathbb{R}},$
(17)

with $\mathbb{x}=(a_{x}\leavevmode\nobreak\ \leavevmode\nobreak\ 0)^{T}$ and
$\mathbb{y}=(0\leavevmode\nobreak\ \leavevmode\nobreak\ a_{y})^{T}$ and where
$a_{x}=a_{y}=a=1$.

We develop the spin-wave theory of the model for the paramagnetic phase $J\ll
h,\gamma$ via the Holstein-Primakoff transformation; we obtain a Hamiltonian
that is quadratic and is best written in momentum space:

$\hat{H}=\frac{1}{2}\sum_{\mathbf{k}}A_{\mathbf{k}}\left(\hat{a}^{{\dagger}}_{\mathbf{k}}\hat{a}_{\mathbf{k}}+\hat{a}_{-\mathbf{k}}\hat{a}^{{\dagger}}_{-\mathbf{k}}\right)+B_{\mathbf{k}}\left(\hat{a}^{{\dagger}}_{\mathbf{k}}\hat{a}^{{\dagger}}_{-\mathbf{k}}+\hat{a}_{\mathbf{k}}\hat{a}_{-\mathbf{k}}\right).$
(18)

Both momentum-dependent functions $A_{\mathbb{k}}$ and $B_{\mathbb{k}}$ are
simple generalisation of the one-dimensional results:

$\displaystyle
A_{\mathbb{k}}=\frac{J}{2}\left(\cos(\mathbb{k}\cdot\mathbb{x})+\cos(\mathbb{k}\cdot\mathbb{y})\right)+h+i\gamma,$
(19a) $\displaystyle
B_{\mathbb{k}}=\frac{J}{2}\left(\cos(\mathbb{k}\cdot\mathbb{x})+\cos(\mathbb{k}\cdot\mathbb{y})\right).$
(19b)

Using the complex bosonic Bogolyubov transformation defined in Eqs. (5a) and
(5b) and adapted to the two-dimensional case, the excitation spectrum
$E_{\mathbb{k}}$ of the non-Hermitian TFIM on a square lattice is given by

$E_{\mathbb{k}}=\mathrm{sgn}(\Re(A_{\mathbb{k}}))\sqrt{A_{\mathbb{k}}^{2}-B_{\mathbb{k}}^{2}},$
(20)

and is plotted on Fig. 2 for different sets of parameters.

As in the 1D case, we recover the symmetry $\mathbb{k}\to-\mathbb{k}$ as well
as the presence of a gap given by $h$ for the real part and $\gamma$ for the
imaginary part when the 2D non-Hermitian quantum system is deep in the $z$
polarized phase i.e. $h\gg J$, see first line on Fig. 2. For the regime where
$h\gtrsim J$, the gap in the imaginary part of $E_{\mathbb{k}}$ is still given
by the dissipation strength $\gamma$ contrary to the one in the real part
which is not given anymore by the transverse magnetic field $h$. Besides, on
the third line of Fig. 2, due to the specific value of the parameters, the
momentum-dependent function $\mathrm{sgn}(\Re(A_{\mathbb{k}}))$ in the
expression of $E_{\mathbb{k}}$ is not defined at the four points corresponding
to the edges of the first Brillouin zone and given by
$\mathbb{k}=(-\pi,-\pi),(-\pi,\pi),(\pi,-\pi),(\pi,\pi)$. Furthermore, for the
same set of parameters, we notice that there is a gap closing for the real
part of the excitation spectrum $\Re(E_{\mathbb{k}})$ meaning that the latter
is gapless. This leads to a breakdown of the bosonic approach based on the
linear spin-wave theory (truncation of Holstein-Primakoff transformation
together with the complex bosonic Bogolyubov transformation).

Figure 2:  Analytical results concerning the complex excitation spectrum of
the non-Hermitian transverse Ising model on a square lattice obtained via the
bosonic approach relying on the Holstein-Primakoff transformation together
with the complex bosonic Bogolyubov transformation. The real and imaginary
part of the excitation spectrum $E_{\mathbb{k}}$, denoted by
$\Re(E_{\mathbb{k}})$ (first column) and $\Im(E_{\mathbb{k}})$ (second column)
respectively, are plotted with respect to the momentum $k_{x}$ and $k_{y}$,
see Eq. (8). The considered numerical parameters are (first line) $J=1$,
$h=20$, $\gamma=10$, (second line) $J=1$, $h=2$, $\gamma=1$, (third line)
$J=h=\gamma=1$.

The generalization to a hypercubic lattice of dimension $D$ is
straightforward. Indeed, the excitation spectrum is still given by the
following equation

$E_{\mathbb{k}}=\mathrm{sgn}(\Re(A_{\mathbb{k}}))\sqrt{A_{\mathbb{k}}^{2}-B_{\mathbb{k}}^{2}},$
(21)

where $\mathbb{k}=(k_{1}\leavevmode\nobreak\ k_{2}\leavevmode\nobreak\
...\leavevmode\nobreak\ k_{D-1}\leavevmode\nobreak\ k_{D})^{T}$ is a
$D$-dimensional vector and where the momentum-dependent functions
$A_{\mathbb{k}}$ and $B_{\mathbb{k}}$ are given by

$\displaystyle
A_{\mathbb{k}}=h+i\gamma+\frac{J}{2}\sum_{i=1}^{D}\cos(\mathbb{k}\cdot\mathbb{d}_{i}),$
(22a) $\displaystyle
B_{\mathbb{k}}=\frac{J}{2}\sum_{i=1}^{D}\cos(\mathbb{k}\cdot\mathbb{d}_{i}),$
(22b)

with $\mathbb{d}_{1}=(a_{1}\leavevmode\nobreak\ 0\leavevmode\nobreak\
...\leavevmode\nobreak\ 0)^{T}$, $\mathbb{d}_{2}=(0\leavevmode\nobreak\
a_{2}\leavevmode\nobreak\ 0\leavevmode\nobreak\ ...\leavevmode\nobreak\
0)^{T}$, $\mathbb{d}_{j}=(0\leavevmode\nobreak\ ...\leavevmode\nobreak\
0\leavevmode\nobreak\ a_{j}\leavevmode\nobreak\ 0\leavevmode\nobreak\
...\leavevmode\nobreak\ 0)^{T}$, $\mathbb{d}_{D}=(0\leavevmode\nobreak\
...\leavevmode\nobreak\ 0\leavevmode\nobreak\ a_{D})^{T}$ and where
$a_{i}=a=1$ with $i\in[|1,D|]$.

##  III Spin-wave theory for the Quench Dynamics

We now investigate the out-of-equilibrium dynamics of the non-Hermitian TFIM
discussing a sudden global quench starting from an initial state
$\ket{\Psi_{0}}$ corresponding to the ground state of the Hermitian TFIM
($\gamma=0$) deep in the paramagnetic phase ($h\gg J$). At time $t=0$ we
switch on the dissipative coupling ($\gamma\neq 0$) and we study the time
evolution with the post-quench non-Hermitian Hamiltonian $\hat{H}$ of the
TFIM. Both pre- and post-quench Hamiltonians are confined relatively deep into
the paramagnetic phase where the spin-wave theory can be safely employed, at
least for the spectrum. As we will see this is not necessary the case for the
long-time dynamics. The initial state can be roughly thought of as a fully-
polarised product state $\ket{\uparrow\uparrow...\uparrow\uparrow}$, although
in the calculations we will use the correct expression, i.e. the ground-state
of the Hermitian TFIM ($\gamma=0$) for $h=5J$ where $J$ will be fixed to
unity. The time-evolved state reads:

$\ket{\Psi(t)}=\frac{e^{-i\hat{H}t}\ket{\Psi_{0}}}{||e^{-i\hat{H}t}\ket{\Psi_{0}}||},$
(23)

and the normalization factor is essential since $\hat{H}$ is non-Hermitian and
reduces the norm of the state, thus making the dynamics a non-linear function
of the state.

Let us begin our discussion considering the case $J=0$ where all lattice sites
are decoupled and focus on the single site dynamics which can be solved
exactly. For a generic initial state
$\ket{\Psi_{0}}=\alpha\ket{\uparrow}+\beta\ket{\downarrow}$, with
$|\alpha|^{2}+|\beta|^{2}=1$, the time-evolved state reads:

$\ket{\Psi(t)}=\displaystyle\frac{\alpha e^{\frac{iht}{2}}e^{-\frac{\gamma
t}{2}}\ket{\uparrow}+\beta e^{-\frac{iht}{2}}e^{\frac{\gamma
t}{2}}\ket{\downarrow}}{\sqrt{|\alpha|^{2}e^{-\gamma t}+|\beta|^{2}e^{\gamma
t}}},$ (24)

from which it follows that, as along as $\alpha$ and $\beta$ are different
from zero, the asymptotic state is $\ket{\downarrow}$ for $\gamma>0$ and
$\ket{\uparrow}$ for $\gamma<0$. This is for instance reflected by the time-
evolution of the projector on the $\ket{\downarrow}$ state:

$\langle\Psi(t)|\hat{S}^{-}\hat{S}^{+}|\Psi(t)\rangle=\frac{|\beta|^{2}}{|\beta|^{2}+|\alpha|^{2}e^{-2\gamma
t}}.$ (25)

The goal of the next sections is to use the spin-wave theory developed so far
in order to discuss the dynamics of the spin excitations in the system when
$J$ is different from zero.

Before dealing with the calculations, we can already anticipate that our
approach will not be able to describe the non-unitary dynamics at all times
for $\gamma>0$. Since the linear spin-wave theory is a good approximation when
the number of bosons $\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle\ll 1$,
the spin state must be close to the paramagnetic state
$\ket{\uparrow\uparrow...\uparrow\uparrow}$. Yet, as we have just seen, for
$\gamma>0$ the long-time dynamics drives the system towards the state
$\ket{\downarrow...\downarrow}$ where the spin-wave theory is not supposed to
work. The study that we are going to perform is thus expected to give a good
quantitative description only at short times. We will discuss this point
extensively. On the other hand, we expect this problem not to arise for
$\gamma<0$ because in the long-time limit the system remains close to the
state $\ket{\uparrow\uparrow...\uparrow\uparrow}$ which is well described in
terms of a small number of bosonic excitations, thus maintaining the system in
a regime where spin-wave theory works.

###  III.1 One-dimensional lattice: equations of motion and quench dynamics

In order to characterize the quench dynamics of our model, we derive the
equations of motion (EoM) associated to the quadratic bosonic correlators:

$G_{k}(t)=\langle\hat{a}^{{\dagger}}_{k}\hat{a}_{k}\rangle_{t},\qquad
F_{k}(t)=\langle\hat{a}_{k}\hat{a}_{-k}\rangle_{t},$ (26)

where $\langle...\rangle_{t}=\langle\Psi(t)|...|\Psi(t)\rangle$ with
$\ket{\Psi(t)}$ given at Eq. (23). Note that
$F_{k}(t)^{*}=\langle\hat{a}^{{\dagger}}_{-k}\hat{a}^{{\dagger}}_{k}\rangle_{t}=\langle\hat{a}^{{\dagger}}_{k}\hat{a}^{{\dagger}}_{-k}\rangle_{t}$
and
$\langle\hat{a}^{{\dagger}}_{k}\hat{a}_{q}\rangle_{t}=\delta_{k,q}G_{k}(t)$
due to the momentum conservation. We obtain the two non-linear coupled
differential equations (see Appendix C for details on the derivation):

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}F_{k}(t)=4\Im(A_{k})F_{k}(t)G_{k}(t)-2iA_{k}F_{k}(t)-iB_{k}\left(1+2G_{k}(t)\right);$
(27a)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}G_{k}(t)=-2B_{k}\Im(F_{k}(t))+2\Im(A_{k})\left(|F_{k}(t)|^{2}+G_{k}(t)+G_{k}(t)^{2}\right).$
(27b)

These equations of motion are valid for any value of $\gamma$, be it positive
or negative. In the rest of this section we will focus on the case $\gamma>0$,
where the breakdown of spin-wave theory is expected, and postpone the
discussion of the case $\gamma<0$ to Sec. III.5.

Figure 3: Space-time plot of the one-body correlation function $G_{R}(t)$ for
a sudden global quench on $\gamma$ of the TFIM in the paramagnetic phase. We
start from the ground state at $\gamma=0$ (Hermitian case) and evolve it in
time with $\gamma>0$ (non-Hermitian case) for the transverse Ising chain in
the $z$ polarized phase. (Top) The results are obtained by integrating the EoM
in Eq. (27). (Bottom) Theoretical guess of the correlation function
characterised by Eq. (29). The parameters of the simulation are: $N_{s}=200$,
$h=5$, $\gamma=0.2$, $J=1$. The slope of the superimposed lines are:
$V_{\mathrm{m}}=-5.2$ and $V_{\mathrm{CE}}=1.3$ (top panel) and
$V_{\mathrm{m}}=-5.8$ and $V_{\mathrm{CE}}=1.3$ (bottom panel).

We now focus on the spatial propagation of the quantum correlations induced by
the dissipative quench by analyzing the bosonic one-body density matrix:

$\displaystyle G_{R}(t)=$
$\displaystyle\langle\hat{a}^{{\dagger}}_{R}\hat{a}_{0}\rangle_{t}=\frac{1}{N_{s}}\sum_{k}\cos(kR)G_{k}(t),$
(28)

that should approximate the $\langle\hat{S}^{-}_{R}\hat{S}^{+}_{0}\rangle_{t}$
spin correlator. Fig. 3 presents the results of our numerical simulations,
where we clearly notice the emergence of a light-cone structure in the space-
time plot of $G_{R}(t)$, also called Lieb-Robinson light-cone [58]. A clear
correlation edge (CE) separates the causal region where the correlations have
non-zero values from the non-causal region where the quantum correlations are
zero. They propagate ballistically with a velocity $V_{\mathrm{CE}}$ which can
be deduced from the slope of the CE, highlighted by a green solid line in Fig.
3. In Appendix E we describe several techniques to numerically estimate
$V_{\rm CE}$.

Moreover, we observe that in the vicinity of the CE the $G_{R}(t)$ displays a
series of local minima and maxima which also propagate ballistically,
highlighted in Fig. 3(top) by several blue dashed lines. This behavior of the
space-time correlations for a non-Hermitian short-range interacting quantum
lattice model is reminiscent of that of Hermitian short-range and long-range
interacting quantum lattice models in the so-called local regime [44; 45; 46].
This regime refers to large values of $\alpha$ which corresponds to the power-
law exponent of the long-range algebraic interactions.

In order to give a clearer interpretation of the space-time pattern of
$G_{R}(t)$, we present a simplified theoretical model based on the following
analytical expression, see Appendix D for more details,

$\displaystyle
G_{R}(t)=\frac{1}{N_{s}}\sum_{k}F_{k}\cos(kR)\cos(2\Re(A_{k})t)\exp(2\Im(A_{k})t),$
(29)

where $A_{k}=h+(J/2)\cos(k)+i\gamma$ and the momentum-dependent function
$F_{k}$ corresponding to the sum of three Gaussians is given by

$\displaystyle F_{k}=$
$\displaystyle\exp\left(-\frac{(k-\pi)^{2}}{2}\right)+\exp\left(-\frac{(k+\pi)^{2}}{2}\right)$
$\displaystyle+\exp\left(-\frac{k^{2}}{2}\right),$ (30)

where $k\in[-\pi,\pi]$. Note that the numerical space-time pattern of
$G_{R}(t)$ deduced from the EoM, see top panel on Fig. 3, is qualitatively
very well reproduced by our theoretical guess, see bottom panel on Fig. 3.
Indeed, the twofold linear structure as well as the two characteristic
velocities $V_{\mathrm{CE}}$ and $V_{\mathrm{m}}$ are well captured.

Concerning the long-time quench dynamics of our model, we find a breakdown of
the linear spin-wave theory. By looking at the expressions of the EoM we could
anticipate an exponential growth of $G_{k}(t)$ as a function of the time $t$,
since its time derivative is proportional to $\Im(A_{k})=\gamma$, see Eq.
(27). This brings quickly the system to a regime where the linear
approximation breaks-down. For instance, for the set of parameters considered
in Fig. 3, we find a value of $\mathrm{max}(G_{k}(t))=\mathrm{max}(\langle
a^{{\dagger}}_{k}a_{k}\rangle_{t})$ significantly lower than $2s=1$. Hence,
for $t\leq t_{\mathrm{f}}$ where $t_{\mathrm{f}}$ refers to the observation
time, the linear spin-wave approximation is well justified. However, it is not
the case anymore when $t>t_{\mathrm{f}}$ because of
$\mathrm{max}(G_{k}(t))\gtrsim 2s$. In addition to this exponential growth we
find a genuine divergence of $G_{k}(t)$, on a time scale which is controlled
by $1/\gamma$. In order to rationalize this surprising result in the next
section we will simplify the problem and study the dynamics of a single
bosonic mode, corresponding to selecting only $k=0$ momentum in the spin wave
theory.

###  III.2 Dynamical Instability in a single-mode problem

To understand the origin of the dynamical instability discussed at the end of
previous section we now investigate the single-mode $k=0$ problem, which is
the simplest situation where we can test the predictions of the spin-wave EoM
developing a comparison with the numerics based on exact diagonalization (ED).
This study will provide a simple intuition concerning the increase of
$G(t)=\langle\hat{a}^{{\dagger}}\hat{a}\rangle_{t}$ and thus the breakdown of
the linear spin-wave approximation at large times, which makes the spin-wave
theory useless. In fact, we will observe that the equations predict a
mathematical singularity at finite times, a problem that will plague all our
studies of dynamics using spin-wave theory when $\gamma>0$.

Let us begin with the benchmark of the EoM of spin-wave theory agains ED; the
single-mode Hamiltonian of the non-Hermitian transverse Ising chain is given
by

$\hat{H}=\frac{A}{2}\left(\hat{a}^{{\dagger}}\hat{a}+\hat{a}\hat{a}^{{\dagger}}\right)+\frac{B}{2}\left(\hat{a}^{{\dagger}}\hat{a}^{{\dagger}}+\hat{a}\hat{a}\right),$
(31)

with $A=h+(J/2)+i\gamma$, $B=J/2$. $\hat{a}$ is a short-hand notation for
$\hat{a}_{k=0}$. We specify the EOMs in (27) to this case and integrate them
numerically using a fourth-order Runge-Kutta method. Similarly, we compute an
ED study of the dynamics of Hamiltonian (31) truncating the Hilbert space at
the maximal number of bosons $N_{\rm max}$. In both cases, the initial state
corresponds to the ground state of the Hamiltonian (31) for $\gamma=0$. The
results obtained with both methods are compared in Fig. 4, where we observe a
perfect agreement as long as we consider low bosonic occupation numbers,
$G(t)\ll 1$, here of the order of $10^{-2}$.

Figure 4: Single-mode two-body correlators as a function of time $t$ where a
sudden global quench has been considered on $\gamma$ from $\gamma=0$ to
$\gamma\neq 0$ for the non-Hermitian transverse Ising chain confined in the
$z$ polarized phase. (a) $G(t)=\langle\hat{a}^{{\dagger}}\hat{a}\rangle_{t}$,
(b) $\Re(F(t))=\Re(\langle\hat{a}\hat{a}\rangle_{t})$, (c)
$\Im(F(t))=\Im(\langle\hat{a}\hat{a}\rangle_{t})$. For each plot, we compare
the results found using exact diagonalization (see dashed red lines) and using
the equation of motion technique (see solid blue lines). Numerical parameters
: $h=5$, $\gamma=0.2$, $J=1$, $N_{\mathrm{max}}=10$.

We now inspect more closely the properties of the solution obtained
integrating the bosonic EoM. We study an initial state
$\ket{\Psi_{0}}=\ket{\Psi(0)}$ that is a single-mode squeezed vacuum (SMSV)
state i.e. $\ket{\Psi_{0}}=\ket{\mathrm{SMSV}}=\hat{S}(\xi)\ket{0}$, where
$\hat{S}(\xi)$ corresponds to the squeezing operator and $\xi=re^{i\phi}$ is
the squeezing parameter. It corresponds to the generic ground state of the
Hamiltonian in Eq. (31) for $\gamma=0$ and can be expanded as

$\ket{\Psi_{0}}=\frac{1}{\sqrt{\cosh(r)}}\sum_{n=0}^{+\infty}\left(-e^{i\phi}\tanh(r)\right)^{n}\frac{\sqrt{(2n)!}}{2^{n}n!}\ket{2n},$
(32)

on the bosonic Fock basis. Since we are interested in the effect of the local
dephasing noise on the quench dynamics, we focus on the simple model

$\hat{H}_{s}=i\frac{\gamma}{2}\hat{a}^{{\dagger}}\hat{a},$ (33)

where the prefactor $1/2$ is for convenience. It corresponds to the anti-
Hermitian part of the Hamiltonian (31).

With an explicit calculation, one easily obtains that the non-normalised state
$\ket{\tilde{\Psi}(t)}=e^{-i\hat{H}_{s}t}\ket{\Psi_{0}}$ reads

$\ket{\tilde{\Psi}(t)}=\frac{1}{\sqrt{\cosh(r)}}\sum_{n=0}^{+\infty}\left(-e^{i\phi}\tanh(r)e^{\gamma
t}\right)^{n}\frac{\sqrt{(2n)!}}{2^{n}n!}\ket{2n}.$ (34)

After some analytical calculations, we obtain an expression for its norm:

$||\ket{\tilde{\Psi}(t)}||^{2}=\frac{1}{\cosh(r)}\sum_{n=0}^{+\infty}\left(\frac{\tanh(r)e^{\gamma
t}}{2}\right)^{2n}\frac{(2n)!}{(n!)^{2}},$ (35)

which can be made more readable using the Stirling approximation $n!\sim
n^{n}$ for large $n$:

$||\ket{\tilde{\Psi}(t)}||^{2}=\frac{1}{\cosh(r)}\sum_{n=0}^{+\infty}\left(\tanh(r)e^{\gamma
t}\right)^{2n}.$ (36)

The latter converges only if $\tanh(r)e^{\gamma t}<1$ where $0<\tanh(r)<1$
since $r>0$. The previous inequality implies the existence of a typical time

$t_{\mathrm{f}}\approx\frac{1}{\gamma}\ln\left(\frac{1}{\tanh(r)}\right),$
(37)

at which the bosonic theory encounters a singularity that is responsible for
the divergence of $G(t)$. Not only thus the quadratic bosonic theory is
problematic because it cannot describe the physical fact that $0\leq\langle
S^{-}S^{+}\rangle_{t}\leq 1$, but it also predicts that $G(t)$ diverges at
finite time. We anticipate here that the divergence at finite time has been
observed in all the numerical simulations of the time evolution of the linear
spin-wave theory performed in this study with $\gamma>0$. In order to have a
large value of $t_{\rm f}$ we have always used a relatively small value of
$\gamma$, as detailed in the captions of the figures. We emphasize that this
divergence would be most likely cured by including non-linear terms in the
bosonic field which can provide a saturation mechanism to the amplification of
bosonic excitations. Still, we expect the bosonic theory to be particularly
inefficient in modelling the long time dynamics of this system.

###  III.3 Square lattice

Figure 5: One-body correlation function $G_{\mathbf{R}}(t)=G_{x,y}(t)$ for a
sudden global quench on $\gamma$ from $\gamma=0$ (Hermitian case) to
$\gamma\neq 0$ (non-Hermitian case) for the transverse Ising model on a square
lattice and confined in the $z$ polarized phase. Numerical results found using
the equation of motion technique to deduce $G_{\mathbf{k}}(t)$, see Eq. (38).
The parameters are : $N_{s}=100$, $h=5$, $\gamma=0.2$, $J=1$,
$t_{\mathrm{i}}=0$, $t_{\mathrm{f}}=10$, $\mathrm{steps}=300$. (first line
from left to right) $t=t_{\mathrm{i}}=0$, $t=t_{\mathrm{f}}/5=2$,
$t=2t_{\mathrm{f}}/5=4$ (second line from left to right)
$t=3t_{\mathrm{f}}/5=6$, $t=4t_{\mathrm{f}}/5=8$, $t=t_{\mathrm{f}}=10$.

We now turn on the two-dimensional case, i.e. the square lattice. Using the
notation employed in Sec. II.3, we introduce the functions
$G_{\mathbb{k}}(t)=\langle\hat{a}^{{\dagger}}_{\mathbf{k}}\hat{a}_{\mathbf{k}}\rangle_{t}$
and
$F_{\mathbb{k}}(t)=\langle\hat{a}_{\mathbb{k}}\hat{a}_{-\mathbb{k}}\rangle_{t}$
which obey the same EoM derived for the one-dimensional case in Eq. (27). The
one-body density matrix reads

$G_{\mathbf{R}}(t)=\frac{1}{N_{s}^{2}}\sum_{\mathbf{k}}\cos(\mathbf{k}\cdot\mathbf{R})G_{\mathbf{k}}(t),$
(38)

where $\mathbf{R}=(x\leavevmode\nobreak\ \leavevmode\nobreak\ y)^{T}$ and is
plotted at different times in Fig. 5 for a sudden global quench from
$\gamma=0$ to $\gamma\neq 0$.

From Fig. 5 dealing with the short-time quench dynamics of our 2D model, we
clearly notice that correlations propagate in the square lattice from the
central lattice site with time. This pattern is characterized an approximately
circular CE, whose radius increases linearly with time. The velocity can be
extracted by fitting the behaviour of the system both along the $x$ and $y$
direction. We track $x^{*}$ and $y^{*}$ corresponding to the maximal value
along the $x$ and $y$ direction of a non-zero correlation for several values
of the time $t$. The data are extracted from Fig. 5 and we find $V_{\rm
CE}\simeq 1.1$, see Appendix G for another numerical technique to determine
the correlation edge velocity $V_{\mathrm{CE}}$.

###  III.4 Fermionic approach for the one-dimensional lattice

In what follows, we consider a fermionic approach to deal with the quench
dynamics of the non-Hermitian transverse Ising chain confined in the $z$
polarized phase. This will permit to compare its predictions to those obtained
by using the bosonic approach discussed at Subsection III.1. Similarly as the
bosonic approach, the fermionic one consists in calculating analytically the
equations of motion associated to the two-body correlators, also called
second-order or quadratic correlators, expressed in terms of the fermionic
annihilation and creation operators denoted by $\hat{c}_{k}$ and
$\hat{c}^{{\dagger}}_{k}$. The two interesting correlators are
$\mathcal{G}_{k}(t)=\langle\hat{c}^{{\dagger}}_{k}\hat{c}_{k}\rangle_{t}$ and
$\mathcal{F}_{k}(t)=\langle\hat{c}_{k}\hat{c}_{-k}\rangle_{t}$.
$\langle...\rangle_{t}=\langle\Psi(t)|...|\Psi(t)\rangle$ and $\ket{\Psi(t)}$
is given at Eq. (23) where the Hamiltonian $\hat{H}$ is defined at Eq. (11).
We obtain the two following complex coupled differential equations [25; 26]

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{F}_{k}(t)=-8\Im(\mathcal{A}_{k})\mathcal{F}_{k}(t)\mathcal{G}_{k}(t)-4i\mathcal{A}_{k}\mathcal{F}_{k}(t)+2i\mathcal{B}_{k}(1-2\mathcal{G}_{k}(t)),$
(39)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{G}_{k}(t)=4\Im(\mathcal{A}_{k})(|\mathcal{F}_{k}(t)|^{2}+\mathcal{G}_{k}(t)-\mathcal{G}_{k}(t)^{2})+4i\mathcal{B}_{k}\Re(\mathcal{F}_{k}(t)),$
(40)

which can be solved numerically using a fourth-order Runge-Kutta method, see
Appendix F for more details about the calculation of the equations of motion.
In the following, two different observables, namely the spin-spin correlation
function denoted by $C^{zz}_{R}(t)$ and defined as

$C^{zz}_{R}(t)=\langle\hat{S}^{z}_{R}\hat{S}^{z}_{0}\rangle_{t}-\langle\hat{S}^{z}_{R}\rangle_{t}\langle\hat{S}^{z}_{0}\rangle_{t},$
(41)

and the local magnetization $\langle\hat{S}^{z}_{R}\rangle_{t}$, are
discussed.

#### Local magnetization -

We now turn to the second observable which is the local magnetization
$\langle\hat{S}^{z}_{R}\rangle_{t}$. Using a similar procedure as the one used
previously for the spin-spin correlations along the $z$ axis, we find for the
bosonic approach the following analytical expression

$\langle\hat{S}^{z}_{R}\rangle_{t}=\frac{1}{2}-\frac{1}{N_{s}}\sum_{k}G_{k}(t),$
(42)

and for the fermionic approach, we obtain a very similar form given by

$\langle\hat{S}^{z}_{R}\rangle_{t}=\frac{1}{2}-\frac{1}{N_{s}}\sum_{k}\mathcal{G}_{k}(t).$
(43)

Note that both theoretical expressions of the local magnetization are
independent of the distance $R$ due to the translational invariance of the
Hamiltonian associated to the transverse-field Ising chain. The theoretical
expression associated to the time-dependent local magnetization using the
bosonic and the fermionic approach is plotted in Fig. 6 for the transverse-
field Ising chain confined in the $z$ polarized phase where a sudden global
quench on $\gamma$ the dissipation strength is considered. We can clearly see
a quantitative agreement between both expressions at short, intermediate and
large times where the maximal relative error is very small. This means that
the bosonic and fermionic approaches provide very similar results as expected
and thus are equivalent even if the transverse-field Ising chain is confined
not very deep into the $z$ polarized phase. Note that at time $t=0$, the local
magnetization is close to $1/2$ due to the fact that the initial Hamiltonian
is confined relatively deep into the paramagnetic phase.

Figure 6: Local magnetization $\langle\hat{S}^{z}_{R}\rangle_{t}$ for a sudden
global quench on $\gamma$ from $\gamma=0$ to $\gamma\neq 0$ for the transverse
Ising chain in the $z$ polarized phase. For the profile associated to the
bosonic approach, see solid blue line and Eq. (42) and for the one associated
to the fermionic approach, see dashed red line and Eq. (43). The parameters
are : $N_{s}=200$, $h=5$, $\gamma=0.2$, $J=1$, $t_{\mathrm{i}}=0$,
$t_{\mathrm{f}}=13$, $\mathrm{steps}=300$.

#### Spin-spin correlation function -

Concerning the spin-spin correlation function $C^{zz}_{R}(t)$, the latter can
be expressed in terms of the two-body correlators $F_{k}(t)$ and $G_{k}(t)$
($\mathcal{F}_{k}(t)$ and $\mathcal{G}_{k}(t)$) when the bosonic (fermionic)
approach is considered. The procedure to follow in order to deduce the
analytical expression of $C^{zz}_{R}(t)$ is very similar for the two
theoretical approaches. Indeed, the first step consists in using the Holstein-
Primakoff (Jordan-Wigner) transformation. The second step is to use the
bosonic (fermionic) Wick theorem for the bosonic (fermionic) creation and
annihilation operators. Then, we rely on the Fourier transform in order to
manipulate the operators in momentum space and finally the momentum
conservation is used. For the bosonic approach, the theoretical expression of
$C^{zz}_{R}(t)$ is given by

$C^{zz}_{R}(t)=\frac{1}{N_{s}^{2}}\sum_{k,k^{\prime}}e^{i(k-k^{\prime})R}(F_{k}(t)^{*}F_{k^{\prime}}(t)+G_{k}(t)(1+G_{k^{\prime}}(t))),$
(44)

where $F_{k}(t)=\langle\hat{a}_{k}\hat{a}_{-k}\rangle_{t}$ and
$G_{k}(t)=\langle\hat{a}^{{\dagger}}_{k}\hat{a}_{k}\rangle_{t}$. Concerning
the fermionic approach, the spin-spin correlation function $C^{zz}_{R}(t)$ has
the following expression

$C^{zz}_{R}(t)=\frac{1}{N_{s}^{2}}\sum_{k,k^{\prime}}e^{i(k-k^{\prime})R}(\mathcal{F}_{k}(t)^{*}\mathcal{F}_{k^{\prime}}(t)+\mathcal{G}_{k}(t)(1-\mathcal{G}_{k^{\prime}}(t))),$
(45)

with $\mathcal{F}_{k}(t)=\langle\hat{c}_{k}\hat{c}_{-k}\rangle_{t}$ and
$\mathcal{G}_{k}(t)=\langle\hat{c}^{{\dagger}}_{k}\hat{c}_{k}\rangle_{t}$.

Figure 7: Space-time spin-spin correlation function $C^{zz}_{R}(t)$ defined at
Eq. (41) for a sudden global quench on $\gamma$ from $\gamma=0$ to $\gamma\neq
0$ for the transverse Ising chain in the $z$ polarized phase. Real part of
$C^{zz}_{R}(t)$ when using the bosonic approach, see Eq. (44) (top panel),
real part of $C^{zz}_{R}(t)$ when considering the fermionic approach, see Eq.
(45) (bottom panel). The parameters are : $N_{s}=200$, $h=5$, $\gamma=0.2$,
$J=1$, $t_{\mathrm{i}}=0$, $t_{\mathrm{f}}=13$, $\mathrm{steps}=300$,
$R_{\mathrm{i}}=0$, $R_{\mathrm{f}}=15$. Note that the logarithmic value of
the real part of $C^{zz}_{R}(t)$ is considered and denoted by
$\mathrm{log}_{10}(\Re(C^{zz}_{R}(t)))$.

The spin-spin correlation function $C^{zz}_{R}(t)$ is plotted on Fig. 7 for
both theoretical approaches. Before discussing the space-time correlation
patterns, the Jordan-Wigner approximation implying
$\mathcal{G}_{k}(t)=\langle\hat{c}^{{\dagger}}_{k}\hat{c}_{k}\rangle_{t}\ll
2s$ and meaning that the quantum spin chain is well confined into the $z$
polarized phase has been verified. Concerning the correlation patterns, one
can clearly notice a qualitative agreement. Indeed, the order of magnitude of
the spin-spin correlations is very similar between the two plots. This means
that both theoretical approaches provide qualitatively the same results and
thus are consistent for the transverse Ising chain relatively deep into the
$z$ polarized phase.

###  III.5 Quench dynamics for negative $\gamma$

In the following, we investigate the sudden global quench dynamics of the TFIM
in the one-dimensional case from $\gamma=0$ to $\gamma<0$. As discussed in
Sec. II this situation can arise for example for a different choice of jump
operators, associated to a projection on the spin-down rather than a spin-up
component. As we have already written, this situation tends to favor the
quantum state $\ket{\uparrow}$ and we expect that the divergences discussed in
the previous sections should disappear.

We begin our study by performing the time-evolution with the single-mode
Hamiltonian in Eq. (31), where our results for $\gamma<0$ are presented in
Fig. 8. We observe a perfect agreement between the results found using the EoM
and the ED. As expected, the dynamics is stabilized and no divergence occurs
anymore. Note that the values of $G(t)$ are very small, which means that the
system is close to the $\ket{\uparrow}$ state.

Figure 8: Single-mode two-body correlators as a function of time $t$ where a
sudden global quench has been considered on $\gamma$ from $\gamma=0$ to
$\gamma<0$ for the non-Hermitian transverse Ising chain confined in the $z$
polarized phase. (a) $G(t)=\langle\hat{a}^{{\dagger}}\hat{a}\rangle_{t}$, (b)
$\Re(F(t))=\Re(\langle\hat{a}\hat{a}\rangle_{t})$, (c)
$\Im(F(t))=\Im(\langle\hat{a}\hat{a}\rangle_{t})$. For each plot, we compare
the results found using exact diagonalization (see dashed red lines) and using
the equation of motion technique (see solid blue lines). Numerical parameters
: $h=5$, $\gamma=-0.2$, $J=1$, $N_{\mathrm{max}}=10$.

We now turn to the multi-mode problem where the time-evolution of the local
magnetization for the 1D spin chain is investigated. In Fig. 9 we present our
results for the local magnetization computed both with the bosonic and the
fermionic approaches; we can observe that they are in qualitative agreement.
In both cases, the local magnetization has values close but slightly below
$1/2$ as time increases because we start and remain close to an almost fully-
polarized product state
$\ket{\Psi_{0}}\simeq\ket{\uparrow\uparrow...\uparrow\uparrow}$. The dynamics
occurs at small times and then stabilizes to a final value larger than the
initial one.

Figure 9: Local magnetization $\langle\hat{S}^{z}_{R}\rangle_{t}$ for a sudden
global quench on $\gamma$ from $\gamma=0$ (Hermitian case) to $\gamma<0$ (non-
Hermitian case) for the transverse Ising chain in the $z$ polarized phase. For
the profile associated to the bosonic approach, see solid blue line and Eq.
(42) and for the one associated to the fermionic approach, see solid red line
and Eq. (43). The parameters of the simulation are: $N_{s}=200$, $h=5$, $J=1$,
$t_{\mathrm{i}}=0$, $t_{\mathrm{f}}=13$, $\mathrm{steps}=300$. (Top)
$\gamma=-0.2$, (bottom) $\gamma=-0.9$.

Using the fermionic approach, the stationary or final value of the local
magnetization is given by the EoMs at Eq. (39) and (40) by fixing each
derivative to zero, i.e. $\mathrm{d}\mathcal{F}_{k}(t)/\mathrm{d}t=0$ and
$\mathrm{d}\mathcal{G}_{k}(t)/\mathrm{d}t=0$. Then, the system to solve reads:

$\displaystyle-8\Im(\mathcal{A}_{k})\mathcal{F}_{k}\mathcal{G}_{k}-4i\mathcal{A}_{k}\mathcal{F}_{k}+2i\mathcal{B}_{k}(1-2\mathcal{G}_{k})$
$\displaystyle=0,$ (46) $\displaystyle
4\Im(\mathcal{A}_{k})(|\mathcal{F}_{k}|^{2}+\mathcal{G}_{k}-\mathcal{G}_{k}^{2})+4i\mathcal{B}_{k}\Re(\mathcal{F}_{k})$
$\displaystyle=0,$ (47)

and it yields the following equation for the correlator $\mathcal{G}_{k}$:

$\displaystyle
16\mathcal{B}_{k}^{2}\Im(\mathcal{A}_{k})+64\Im(\mathcal{A}_{k})(|\mathcal{A}_{k}|^{2}-\mathcal{B}_{k}^{2})\mathcal{G}_{k}+512\Im(\mathcal{A}_{k})^{3}\mathcal{G}_{k}^{3}$
$\displaystyle+(64\Im(\mathcal{A}_{k})(\mathcal{B}_{k}^{2}-|\mathcal{A}_{k}|^{2})-256\Im(\mathcal{A}_{k})^{3})\mathcal{G}_{k}^{2}-256\Im(\mathcal{A}_{k})^{3}\mathcal{G}_{k}^{4}=0.$
(48)

When solving the previous equation for each momentum $k$ and by taking the
real and lowest solution for $\mathcal{G}_{k}$, the stationary value of the
local magnetization $\langle\hat{S}^{z}_{R}\rangle_{\mathrm{f}}$ (the index
’$\mathrm{f}$’ stands for ’fermionic’) reads as:

$\langle\hat{S}^{z}_{R}\rangle_{\mathrm{f}}=\frac{1}{2}-\frac{1}{N_{s}}\sum_{k}\mathcal{G}_{k}.$
(49)

Using the latter, we recover the value of the stationary local magnetization
on Fig. 9 when the same sets of parameters are considered. To extract the
stationary local magnetization when using the fermionic approach, see the
value of the local magnetization at large times on solid red lines.

Similarly, using the bosonic approach, the stationary value of the local
magnetization is given by the EoMs at Eq. (27a) and (27b) by fixing each
derivative to zero, i.e. $\mathrm{d}F_{k}(t)/\mathrm{d}t=0$ and
$\mathrm{d}G_{k}(t)/\mathrm{d}t=0$. The system to solve reads:

$\displaystyle 4\Im(A_{k})F_{k}G_{k}-2iA_{k}F_{k}-iB_{k}(1+2G_{k})$
$\displaystyle=0,$ (50)
$\displaystyle-2B_{k}\Im(F_{k})+2\Im(A_{k})(|F_{k}|^{2}+G_{k}+G_{k}^{2})$
$\displaystyle=0,$ (51)

which yields for $G_{k}$ the following equation:

$\displaystyle-2B_{k}^{2}\Im(A_{k})+8\Im(A_{k})(|A_{k}|^{2}-B_{k}^{2})G_{k}+64\Im(A_{k})^{3}G_{k}^{3}$
$\displaystyle+8\Im(A_{k})(|A_{k}|^{2}-B_{k}^{2}+4\Im(A_{k})^{2})G_{k}^{2}+32\Im(A_{k})^{3}G_{k}^{4}=0.$
(52)

When solving the previous equation for each momentum $k$ and by taking the
real and lowest positive solution for $\mathcal{G}_{k}$, the stationary value
of the local magnetization $\langle\hat{S}^{z}_{R}\rangle_{\mathrm{b}}$ (the
index ’$\mathrm{b}$’ stands for ’bosonic’) reads as:

$\langle\hat{S}^{z}_{R}\rangle_{\mathrm{b}}=\frac{1}{2}-\frac{1}{N_{s}}\sum_{k}G_{k}.$
(53)

Using the latter, we recover the value of the stationary local magnetization
on Fig. 9 when the same sets of parameters are considered. To extract the
stationary local magnetization when using the bosonic approach, see the value
of the local magnetization at large times on solid blue lines.

We now focus on the spatial propagation of the quantum correlations induced by
the dissipative quench for $\gamma<0$ by analyzing the bosonic one-body
correlation function defined at Eq. (28). We remind that the latter
approximates the $\langle\hat{S}^{-}_{R}\hat{S}^{+}_{0}\rangle_{t}$ spin
correlator. Fig. 10 presents the results of our numerical simulations, where
we clearly notice the emergence of a light-cone structure in the space-time
plot of $G_{R}(t)$. A clear CE propagating ballistically with a positive
velocity $V_{\mathrm{CE}}$ is visible and highlighted by a green solid line in
Fig. 10. Besides, in the vicinity of the CE the correlation function
$G_{R}(t)$ displays a series of local minima and maxima which also spread
ballistically at the negative velocity $V_{\mathrm{m}}$, highlighted in Fig.
10 by several blue dashed lines. This behavior of the space-time correlations
is reminiscent of the one found previously at Fig. 3 for a dissipative quench
with $\gamma>0$.

Figure 10: Space-time plot of the one-body correlation function $G_{R}(t)$ for
a sudden global quench on $\gamma$ of the TFIM in the paramagnetic phase. We
start from the ground state at $\gamma=0$ (Hermitian case) and evolve it in
time with $\gamma<0$ (non-Hermitian case) for the transverse Ising chain in
the $z$ polarized phase. (Top) The results are obtained by integrating the EoM
in Eq. (27). (Bottom) Theoretical guess of the correlation function
characterised by Eq. (29). The parameters of the simulation are: $N_{s}=200$,
$h=5$, $\gamma=-0.2$, $J=1$. The slope of the superimposed lines are:
$V_{\mathrm{m}}=-4$ and $V_{\mathrm{CE}}=1$ (top panel) and
$V_{\mathrm{m}}=-5.8$ and $V_{\mathrm{CE}}=1$ (bottom panel).

##  IV Conclusion

In this work we have discussed the linear spin-wave theory for the non-
Hermitian TFIM in dimension $d=1,2,3$. Specifically we have first introduced
the Holstein-Primakoff transformation and then performed a linear
approximation to reduce the TFIM to a gas of non-Hermitian bosonic modes,
diagonalizable with a complex bosonic Bogolyubov transformation. We have shown
that this approach works well for the spectrum of the non-Hermitian
quasiparticles. For the one-dimensional case, the quasiparticle dispersion
relation found using the bosonic approach is in very good agreement with the
exact one, that can be computed using a fermionic formalism based on the
Jordan-Wigner transformation. This motivates the use of our theory for higher-
dimensional lattices where exact solutions are not known.

We have then applied the linear spin wave theory to describe the non-Hermitian
dynamics after a quench of the dissipative (non-Hermitian) coupling at time
$t=0$. We have shown that for certain initial states and values of the
dissipation the linear spin-wave approximation breaks down. While one could
have anticipated an amplification of quasiparticle occupation numbers due to
the non-Hermitian dynamics, leading to a violation of the harmonic
approximation, we have shown that the breakdown is more severe and takes the
form of a genuine finite-time divergence that we have explained using a simple
analytical model. We emphasize that higher order terms involving scattering
between bosonic modes, which are disregarded within the linear spin wave
approximation, would likely cure this divergence even though might be still
not sufficient to describe faithfully the non-Hermitian dynamics within a
purely bosonic model. We leave a detailed analysis of this point to the
future.

We used the developed theory at short times to study the dissipative Lieb-
Robinson light-cone of quantum correlations induced by the non-Hermitian
quench both in the one-dimensional and two-dimensional cases.
This work opens a number of questions about the possibility of applying spin-
wave theory to study the dynamics of open quantum systems. In this simple non-
Hermitian context we have shown that a naive application of spin-wave theory
is limited to short times when $\gamma>0$. Finding techniques for extending
the results to longer times is an intriguing open project.

###### Acknowledgements.

This research was supported by LabEx PALM (ANR-10-LABX-0039-PALM) in Orsay and
by Region Ile-de-France in the framework of DIM QuanTip.

##  Appendix A Bosonic and fermionic Bogolyubov transformations

Concerning the complex bosonic Bogolyubov transformation for the one-
dimensional case, the momentum-dependent function $\theta_{k}$ has the
following expression

$\theta_{k}=\mathrm{arctanh}\left(-B_{k}/A_{k}\right),$ (54)

where $A_{k}=h+(J/2)\cos(k)+i\gamma$ and $B_{k}=(J/2)\cos(k)$, see Fig. 11 for
an example. For the complex fermionic Bogolyubov transformation, $\theta_{k}$
is given by

$\theta_{k}=\mathrm{arctan}\left(i\mathcal{B}_{k}/\mathcal{A}_{k}\right),$
(55)

where $\mathcal{A}_{k}=(J/4)\cos(k)+(h+i\gamma)/2$ and
$\mathcal{B}_{k}=-i(J/4)\sin(k)$, see Fig. 11 for an example.

Figure 11: Momentum-dependent function $\theta_{k}$ where $\Re(\theta_{k})$
(solid cyan line) and $\Im(\theta_{k})$ (solid orange line) are plotted for
the complex bosonic Bogolyubov transformation (left) and the complex fermionic
Bogolyubov transformation (right). The numerical parameters are : $J=1$,
$h=20$ and $\gamma=10$.

##  Appendix B Fourier transform of the bosonic and fermionic operators

The convention of the Fourier transform is the following for the bosonic and
fermionic operators respectively for a $D$-dimensional hypercubic lattice

$\displaystyle\hat{a}^{{\dagger}}_{\mathbf{R}}=\frac{1}{N_{s}^{D/2}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\mathbf{R}}\hat{a}^{{\dagger}}_{\mathbf{k}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{a}_{\mathbf{R}}=\frac{1}{N_{s}^{D/2}}\sum_{\mathbf{k}}e^{i\mathbf{k}\mathbf{R}}\hat{a}_{\mathbf{k}},$
(56)
$\displaystyle\hat{c}^{{\dagger}}_{\mathbf{R}}=\frac{1}{N_{s}^{D/2}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\mathbf{R}}\hat{c}^{{\dagger}}_{\mathbf{k}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{c}_{\mathbf{R}}=\frac{1}{N_{s}^{D/2}}\sum_{\mathbf{k}}e^{i\mathbf{k}\mathbf{R}}\hat{c}_{\mathbf{k}}.$
(57)

##  Appendix C Equations of motion - bosonic approach

Figure 12: Numerical momentum-time correlator
$G_{k}(t)=\langle\hat{a}^{{\dagger}}_{k}\hat{a}_{k}\rangle_{t}$ for a sudden
global quench on $\gamma$ from $\gamma=0$ (Hermitian case) to $\gamma\neq 0$
(non-Hermitian case) for the transverse Ising chain in the $z$ polarized
phase. The numerical parameters are : $N_{s}=1000$, $h=5$, $\gamma=0.2$,
$J=1$.

Starting from the time-evolved quantum state $\ket{\Psi(t)}$ given by

$\ket{\Psi(t)}=\frac{e^{-i\hat{H}t}\ket{\Psi_{0}}}{||e^{-i\hat{H}t}\ket{\Psi_{0}}||},$
(58)

we calculate its derivative with respect to the time $t$ which can be written
as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\ket{\Psi(t)}=-i\hat{H}\ket{\Psi(t)}-\frac{i}{2}\langle\Psi(t)|\hat{H}^{{\dagger}}-\hat{H}|\Psi(t)\rangle\ket{\Psi(t)}.$
(59)

Using the latter and the equation below

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bra{\Psi(t)}=\frac{\mathrm{d}}{\mathrm{d}t}\ket{\Psi(t)}^{{\dagger}}=\left(\frac{\mathrm{d}}{\mathrm{d}t}\ket{\Psi(t)}\right)^{{\dagger}},$
(60)

we can deduce the equation of motion for the time-dependent expectation value
of the observable $\hat{A}$ denoted by
$G_{A}(t)=\langle\Psi(t)|\hat{A}|\Psi(t)\rangle=\langle\hat{A}\rangle_{t}$ and
given by

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}G_{A}(t)=i\langle\hat{H}^{{\dagger}}\hat{A}-\hat{A}\hat{H}\rangle_{t}+i\langle\hat{H}-\hat{H}^{{\dagger}}\rangle_{t}G_{A}(t).$
(61)

The Hamiltonian $\hat{H}$ used to perform the time evolution is given at Eq.
(4). After some analytical calculations where the bosonic Wick theorem
together with the momentum conservation are used, we find the coupled complex
differential equations at Eq. (27a) and (27b) for the correlator
$F_{k}(t)=\langle\hat{a}_{k}\hat{a}_{-k}\rangle_{t}$ and
$G_{k}(t)=\langle\hat{a}^{{\dagger}}_{k}\hat{a}_{k}\rangle_{t}$ respectively.
The initial conditions $F_{k}(0)=\langle\hat{a}_{-k}\hat{a}_{k}\rangle_{0}$
and $G_{k}(0)=\langle\hat{a}^{{\dagger}}_{k}\hat{a}_{k}\rangle_{0}$, with
$\langle...\rangle_{0}=\langle\Psi_{0}|...|\Psi_{0}\rangle$ and
$\ket{\Psi_{0}}=\ket{\Psi(0)}=\ket{\mathrm{GS}(\hat{H}_{i})}$ where the
initial (pre-quench) Hermitian Hamiltonian $\hat{H}_{i}$ has the following
expression:

$\displaystyle\hat{H}_{i}=\frac{1}{2}\sum_{k}A_{k,i}\left(\hat{a}^{{\dagger}}_{k,i}\hat{a}_{k,i}+\hat{a}_{-k,i}\hat{a}_{-k,i}^{{\dagger}}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+B_{k}\left(\hat{a}^{{\dagger}}_{k,i}\hat{a}^{{\dagger}}_{-k,i}+\hat{a}_{k,i}\hat{a}_{-k,i}\right),$
(62)

where $A_{k,i}=h+(J/2)\cos(k)$ and $B_{k}=(J/2)\cos(k)$, can be deduced by
following the procedure below. In order to calculate the correlators at time
$t=0$, we rely on the following bosonic Bogolyubov transformation

$\displaystyle\hat{a}_{k,i}=u_{k,i}\hat{\beta}_{k,i}+v_{-k,i}\hat{\beta}^{{\dagger}}_{-k,i},$
(63)
$\displaystyle\hat{a}_{k,i}^{{\dagger}}=u_{k,i}\hat{\beta}^{{\dagger}}_{k,i}+v_{-k,i}\hat{\beta}_{-k,i},$
(64)

where the real momentum-dependent functions $u_{k,i}$ and $v_{-k,i}$ are
parametrized by $\alpha_{k,i}$ and given by

$\displaystyle u_{k,i}=\cosh(\alpha_{k,i}),$ (65) $\displaystyle
v_{-k,i}=\sinh(\alpha_{k,i}),$ (66)

where $\alpha_{k,i}$ is defined as

$\alpha_{k,i}=\frac{1}{2}\mathrm{arctanh}\left(-\frac{B_{k}}{A_{k,i}}\right).$
(67)

By replacing the bosonic annihilation and creation operators $\hat{a}_{k}$ and
$\hat{a}^{{\dagger}}_{k}$ by $\hat{a}_{k,i}$ and $\hat{a}^{{\dagger}}_{k,i}$
respectively in the correlators and by injecting the expression of
$\hat{a}_{k,i}$ and $\hat{a}^{{\dagger}}_{k,i}$ in terms of the bosonic
Bogolyubov operators denoted by $\hat{\beta}_{k,i}$ and
$\hat{\beta}^{{\dagger}}_{k,i}$, we find:

$\displaystyle F_{k}(0)=u_{k,i}v_{k,i},$ (68) $\displaystyle
G_{k}(0)=v_{k,i}^{2}.$ (69)

##  Appendix D Theoretical guess for the one-body correlation function
$G_{R}(t)$

In order to unveil a theoretical expression for $G_{R}(t)$, it requires to
make an analytical guess for $G_{k}(t)$ according to Eq. (28). From Fig. 12
where the momentum-time behavior of $G_{k}(t)$ is displayed, we can make the
following theoretical guess:

$\displaystyle G_{k}(t)=F_{k}\cos(2\Re(E_{k})t)\exp(2\Im(E_{k})t),$ (70)

where the momentum-dependent $F_{k}$ is defined as

$\displaystyle
F_{k}=\exp\left(-\frac{(k-\pi)^{2}}{2\sigma^{2}}\right)+\exp\left(-\frac{(k+\pi)^{2}}{2\sigma^{2}}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\exp\left(-\frac{k^{2}}{2\sigma^{2}}\right).$ (71)

Note that the parameter $\sigma$ is deduced directly and qualitatively from
the numerical result of $G_{k}(t)$ on Fig. 12 by using the relation
$\mathrm{FWHM}\simeq 2.35\sigma$. We obtain $\sigma=1$. As a consequence, the
theoretical guess of $G_{R}(t)$ is given by

$\displaystyle
G_{R}(t)=\frac{1}{N_{s}}\sum_{k}F_{k}\cos(kR)\cos(2\Re(E_{k})t)\exp(2\Im(E_{k})t),$
(72)

where $E_{k}$ denotes the final (post-quench) excitation spectrum
corresponding to the quasiparticle dispersion relation of the non-Hermitian
transverse Ising chain, see Eq. (8). In order to simplify Eq. (72), we
consider the approximation $E_{k}\approx A_{k}$ valid for $h\gg J$ which leads
to the final expression:

$\displaystyle
G_{R}(t)=\frac{1}{N_{s}}\sum_{k}F_{k}\cos(kR)\cos(2\Re(A_{k})t)\exp(2\Im(A_{k})t),$
(73)

where the momentum-dependent function $A_{k}$ is given by
$A_{k}=h+(J/2)\cos(k)+i\gamma$. Note that the approximation $E_{k}\approx
A_{k}$ works very well in the regime where $h\gtrsim J$ according to Fig. 13
where the ratio $h/J=5$ is considered and for which the two space-time
patterns of $G_{R}(t)$ are identical.

Figure 13: Space-time one-body correlation function $G_{R}(t)$, see Eq. (28),
for a sudden global quench on $\gamma$ from $\gamma=0$ (Hermitian case) to
$\gamma\neq 0$ (non-Hermitian case) for the transverse Ising chain in the $z$
polarized phase. Plot of the theoretical guess at Eq. (72)(left) and the
theoretical guess with the approximation $E_{k}\approx A_{k}$ (right), see Eq.
(73). The parameters are : $N_{s}=200$, $h=5$, $\gamma=0.2$, $J=1$,
$t_{\mathrm{i}}=0$, $t_{\mathrm{zoom}}=2.5$, $t_{\mathrm{f}}=13$,
$\mathrm{steps}=300$, $R_{\mathrm{i}}=8$, $R_{\mathrm{f}}=20$, $\sigma=1$.

##  Appendix E Numerical technique to track the correlation edge

The correlation edge is found by tracking for each lattice $R$ the time
$t_{R}$ also called activation time where the correlation reaches a very small
percentage of the maximal value of the correlations, see black dots on Fig.
14. The different points $(R,t_{R})$ form a linear profile, see solid green
line on Fig. 14, where the corresponding velocity $V_{\mathrm{CE}}$ is given
by the slope.

Figure 14: Space-time one-body correlation function $G_{R}(t)$ for a sudden
global quench on $\gamma$ from $\gamma=0$ (Hermitian case) to $\gamma\neq 0$
(non-Hermitian case) for the transverse Ising chain in the $z$ polarized
phase. Numerical result found using the equation of motion technique to deduce
$G_{k}(t)$, see Eq. (28) (left) and theoretical guess, see Eq. (29) (right).
The parameters are : $N_{s}=200$, $h=5$, $\gamma=0.2$, $J=1$,
$t_{\mathrm{i}}=0$, $t_{\mathrm{zoom}}=2.5$, $t_{\mathrm{f}}=13$,
$\mathrm{steps}=300$, $R_{\mathrm{i}}=8$, $R_{\mathrm{f}}=20$, $\sigma=1$.
$V_{\mathrm{CE}}=1.3$ is found for each plot.

##  Appendix F Equations of motion - fermionic approach

We start from the equation of motion for the time-dependent expectation value
of the observable $\hat{A}$ denoted by
$G_{A}(t)=\langle\Psi(t)|\hat{A}|\Psi(t)\rangle=\langle\hat{A}\rangle_{t}$ and
given by

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}G_{A}(t)=i\langle\hat{H}^{{\dagger}}\hat{A}-\hat{A}\hat{H}\rangle_{t}+i\langle\hat{H}-\hat{H}^{{\dagger}}\rangle_{t}G_{A}(t).$
(74)

The latter is derived in Appendix C. The Hamiltonian $\hat{H}$ used to perform
the time evolution is given at Eq. (11). After some analytical calculations
where the fermionic Wick theorem together with the momentum conservation are
used, we find the two coupled complex differential equations at Eq. (39) and
(40) for the correlator
$\mathcal{F}_{k}(t)=\langle\hat{c}_{k}\hat{c}_{-k}\rangle_{t}$ and
$\mathcal{G}_{k}(t)=\langle\hat{c}^{{\dagger}}_{k}\hat{c}_{k}\rangle_{t}$
respectively. The initial conditions
$\mathcal{F}_{k}(0)=\langle\hat{c}_{k}\hat{c}_{-k}\rangle_{0}$ and
$\mathcal{G}_{k}(0)=\langle\hat{c}^{{\dagger}}_{k}\hat{c}_{k}\rangle_{0}$,
with $\langle...\rangle_{0}=\langle\Psi_{0}|...|\Psi_{0}\rangle$ and
$\ket{\Psi_{0}}=\ket{\Psi(0)}=\ket{\mathrm{GS}(\hat{H}_{i})}$ where the
initial (pre-quench) Hermitian Hamiltonian $\hat{H}_{i}$ has the following
expression

$\displaystyle\hat{H}_{i}=\sum_{k}\mathcal{A}_{k,i}\left(\hat{c}^{{\dagger}}_{k,i}\hat{c}_{k,i}-\hat{c}_{-k,i}\hat{c}^{{\dagger}}_{-k,i}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\mathcal{B}_{k}\left(\hat{c}^{{\dagger}}_{k,i}\hat{c}^{{\dagger}}_{-k,i}+\hat{c}_{k,i}\hat{c}_{-k,i}\right),$
(75)

where $\mathcal{A}_{k,i}=(J/4)\cos(k)+h/2$ and
$\mathcal{B}_{k}=-i(J/4)\sin(k)$, can be deduced by following the procedure
below. In order to calculate the correlators at time $t=0$, we rely on the
following fermionic Bogolyubov transformation

$\displaystyle\hat{c}^{{\dagger}}_{-k,i}=\cos\left(\frac{\theta_{k,i}}{2}\right)\hat{\eta}^{{\dagger}}_{-k,i}+i\sin\left(\frac{\theta_{k,i}}{2}\right)\hat{\eta}_{k,i},$
(76)
$\displaystyle\hat{c}_{k,i}=\cos\left(\frac{\theta_{k,i}}{2}\right)\hat{\eta}_{k,i}+i\sin\left(\frac{\theta_{k,i}}{2}\right)\hat{\eta}^{{\dagger}}_{-k,i},$
(77)

where the real momentum-dependent function $\theta_{k,i}$ is defined as

$\theta_{k,i}=\mathrm{arctan}\left(i\frac{\mathcal{B}_{k}}{\mathcal{A}_{k,i}}\right).$
(78)

By replacing the fermionic annihilation and creation operators $\hat{c}_{k}$
and $\hat{c}^{{\dagger}}_{k}$ by $\hat{c}_{k,i}$ and
$\hat{c}^{{\dagger}}_{k,i}$ respectively in the correlators and by injecting
the expression of $\hat{c}_{k,i}$ and $\hat{c}^{{\dagger}}_{k,i}$ in terms of
the fermionic Bogolyubov operators denoted by $\hat{\eta}_{k,i}$ and
$\hat{\eta}^{{\dagger}}_{k,i}$, we find

$\displaystyle\mathcal{F}_{k}(0)=-i\sin\left(\frac{\theta_{k,i}}{2}\right)\cos\left(\frac{\theta_{k,i}}{2}\right),$
(79)
$\displaystyle\mathcal{G}_{k}(0)=\sin^{2}\left(\frac{\theta_{k,i}}{2}\right).$
(80)

##  Appendix G Correlation edge velocity - square lattice

We propose here another numerical technique in order to determine the
correlation edge velocity $V_{\mathrm{CE}}$. The correlations are represented
as a function of $|R|$ corresponding to the distance from the central site of
the square lattice. The latter are reported on Fig. 15 for different times.
For each value of the time, the distance $|R|^{*}$ is tracked such that above
this value the correlations are zero corresponding to the black dots. These
dots have a linear profile fitted by the solid green line where the slope
corresponds to the solid green line. We find $V_{\mathrm{CE}}\simeq 1.2$.

Figure 15: One-body correlation function $G_{|R|}(t)$ for a sudden global
quench on $\gamma$ from $\gamma=0$ (Hermitian case) to $\gamma\neq 0$ (non-
Hermitian case) for the transverse Ising model on a square lattice and
confined in the $z$ polarized phase. $|R|$ represents the distance from the
central site of the square lattice. The parameters are : $N_{s}=100$, $h=5$,
$\gamma=0.2$, $J=1$, $t_{\mathrm{i}}=0$, $t_{\mathrm{f}}=10$,
$\mathrm{steps}=300$.

## References

  * [1] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F. Roos, _Quasiparticle engineering and entanglement propagation in a quantum many-body system_ , Nature 511, 202 (2014).
  * [2] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe, _Non-local propagation of correlations in quantum systems with long-range interactions_ , Nature 511, 198 (2014).
  * [3] M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauss, T. Fukuhara, C. Gross, I. Bloch, C. Kollath, and S. Kuhr, _Light-cone-like spreading of correlations in a quantum many-body system_ , Nature (London) 481, 484 (2012).
  * [4] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, _Colloquium: Nonequilibrium dynamics of closed interacting quantum systems_ , Rev. Mod. Phys. 83, 863 (2011).
  * [5] C. Gogolin and J. Eisert, _Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems_ , Reports on Progress in Physics 79(5), 056001 (2016).
  * [6] P. Calabrese and J. Cardy, _Evolution of entanglement entropy in one-dimensional systems_ , Journal of Statistical Mechanics: Theory and Experiment 2005(04), P04010 (2005).
  * [7] V. Eisler, _Crossover between ballistic and diffusive transport: the quantum exclusion process_ , Journal of Statistical Mechanics: Theory and Experiment 2011(06), P06007 (jun 2011).
  * [8] J.-S. Bernier, R. Tan, L. Bonnes, C. Guo, D. Poletti, and C. Kollath, _Light-cone and diffusive propagation of correlations in a many-body dissipative system_ , Physical Review Letters 120(2) (2018).
  * [9] X. Turkeshi and M. Schiró, _Diffusion and thermalization in a boundary-driven dephasing model_ , Phys. Rev. B 104, 144301 (Oct 2021).
  * [10] V. Alba and F. Carollo, _Noninteracting fermionic systems with localized losses: Exact results in the hydrodynamic limit_ , Physical Review B 105(5) (2022).
  * [11] L. Rosso, D. Rossini, A. Biella, and L. Mazza, _One-dimensional spin-1/2 fermionic gases with two-body losses: Weak dissipation and spin conservation_ , Phys. Rev. A 104, 053305 (2021).
  * [12] L. Rosso, A. Biella, J. De Nardis, and L. Mazza, _Dynamical theory for one-dimensional fermions with strong two-body losses: Universal non-hermitian zeno physics and spin-charge separation_ , Phys. Rev. A 107, 013303 (2023).
  * [13] G. Mazza and M. Schirò, _Dissipative dynamics of a fermionic superfluid with two-body losses_ , Phys. Rev. A 107, L051301 (May 2023).
  * [14] V. Alba and F. Carollo, _Spreading of correlations in markovian open quantum systems_ , Physical Review B 103(2) (2021).
  * [15] H.-P. Breuer and F. Petruccione, _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * [16] Y. Ashida, Z. Gong, and M. Ueda, _Non-hermitian physics_ , Advances in Physics 69(3), 249 (2020).
  * [17] A. J. Daley, _Quantum trajectories and open many-body quantum systems_ , Advances in Physics 63(2), 77 (2014).
  * [18] W. D. Heiss, _The physics of exceptional points_ , J. Phys. Math. Theor. 45, 444016 (2012).
  * [19] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, _Topological phases of non-hermitian systems_ , Phys. Rev. X. 8, 031079 (2018).
  * [20] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, _Exceptional topology of non-hermitian systems_ , Rev. Mod. Phys. 93, 015005 (2021).
  * [21] R. Couvreur, J. L. Jacobsen, and H. Saleur, _Entanglement in nonunitary quantum critical spin chains_ , Phys. Rev. Lett. 119, 040601 (2017).
  * [22] L. Herviou, N. Regnault, and J. H. Bardarson, _Entanglement spectrum and symmetries in non-hermitian fermionic non-interacting models_ , Scipost. Phys. 7 (2019).
  * [23] P.-Y. Chang, J.-S. You, X. Wen, and S. Ryu, _Entanglement spectrum and entropy in topological non-hermitian systems and nonunitary conformal field theory_ , Phys. Rev. Res. 2, 033069 (2020).
  * [24] S. Gopalakrishnan and M. J. Gullans, _Entanglement and purification transitions in non-hermitian quantum mechanics_ , Phys. Rev. Lett. 126, 170503 (2021).
  * [25] X. Turkeshi and M. Schiró, _Entanglement and correlation spreading in non-hermitian spin chains_ , Physical Review B 107(2) (2023).
  * [26] Y. L. Gal, X. Turkeshi, and M. Schirò, _Volume-to-area law entanglement transition in a non-hermitian free fermionic chain_ , SciPost Physics 14(5) (2023).
  * [27] K. Kawabata, T. Numasawa, and S. Ryu, _Entanglement phase transition induced by the non-hermitian skin effect_ , Phys. Rev. X 13, 021007 (2023).
  * [28] E. Granet, C. Zhang, and H. Dreyer, _Volume-law to area-law entanglement transition in a nonunitary periodic gaussian circuit_ , Phys. Rev. Lett. 130, 230401 (Jun 2023).
  * [29] L. Su, A. Clerk, and I. Martin, _Dynamics and phases of nonunitary floquet transverse-field ising model_ (2023), eprint 2306.07428.
  * [30] C. Zhang and E. Granet, _Antiunitary symmetry breaking and a hierarchy of purification transitions in floquet non-unitary circuits_ (2023), eprint 2307.07003.
  * [31] B. Dóra and C. P. Moca, _Quantum quench in $\mathcal{P}\mathcal{T}$-symmetric luttinger liquid_, Phys. Rev. Lett. 124, 136802 (2020).
  * [32] T. Banerjee and K. Sengupta, _Emergent conservation in the floquet dynamics of integrable non-hermitian models_ , Phys. Rev. B 107, 155117 (Apr 2023).
  * [33] N. Hatano and D. R. Nelson, _Localization transitions in non-hermitian quantum mechanics_ , Phys. Rev. Lett. 77, 570 (1996).
  * [34] S.-B. Zhang, M. M. Denner, T. c. v. Bzdušek, M. A. Sentef, and T. Neupert, _Symmetry breaking and spectral structure of the interacting hatano-nelson model_ , Phys. Rev. B 106, L121102 (2022).
  * [35] G. A. Starkov, M. V. Fistul, and I. M. Eremin, _Quantum phase transitions in non-hermitian pt-symmetric transverse-field ising spin chains_ , Annals of Physics p. 169268 (2023).
  * [36] S. Ghosh, K. Sengupta, and I. Paul, _Hilbert space fragmentation imposed real spectrum of a non-hermitian system_ (2023), eprint 2307.05679.
  * [37] K. Yamamoto, M. Nakagawa, M. Tezuka, M. Ueda, and N. Kawakami, _Universal properties of dissipative tomonaga-luttinger liquids: Case study of a non-hermitian xxz spin chain_ , Phys. Rev. B 105, 205125 (2022).
  * [38] P. Kattel, P. R. Pasnoori, and N. Andrei, _Exact solution of a non-hermitian -symmetric spin chain_ , Journal of Physics A: Mathematical and Theoretical 56(32), 325001 (2023).
  * [39] S. Dürr, J. J. García-Ripoll, N. Syassen, D. M. Bauer, M. Lettner, J. I. Cirac, and G. Rempe, _Lieb-liniger model of a dissipation-induced tonks-girardeau gas_ , Phys. Rev. A 79, 023614 (2009).
  * [40] T. Holstein and H. Primakoff, _Field dependence of the intrinsic domain magnetization of a ferromagnet_ , Phys. Rev. 58, 1098 (1940).
  * [41] F. J. Dyson, _General theory of spin-wave interactions_ , Phys. Rev. 102, 1217 (1956).
  * [42] M. Sandri, M. Schiró, and M. Fabrizio, _Linear ramps of interaction in the fermionic hubbard model_ , Phys. Rev. B 86, 075122 (Aug 2012).
  * [43] A. Lerose, B. Žunkovič, J. Marino, A. Gambassi, and A. Silva, _Impact of nonequilibrium fluctuations on prethermal dynamical phase transitions in long-range interacting spin chains_ , Phys. Rev. B 99, 045128 (Jan 2019).
  * [44] L. Cevolani, J. Despres, G. Carleo, L. Tagliacozzo, and L. Sanchez-Palencia, _Universal scaling laws for correlation spreading in quantum systems with short- and long-range interactions_ , Phys. Rev. B 98, 024302 (2018).
  * [45] J. Despres, L. Villa, and L. Sanchez-Palencia, _Twofold correlation spreading in a strongly correlated lattice bose gas_ , Scientific Reports 9(1) (2019).
  * [46] J. T. Schneider, J. Despres, S. J. Thomson, L. Tagliacozzo, and L. Sanchez-Palencia, _Spreading of correlations and entanglement in the long-range transverse ising chain_ , Physical Review Research 3(1) (2021).
  * [47] P. Barmettler, D. Poletti, M. Cheneau, and C. Kollath, _Propagation front of correlations in an interacting bose gas_ , Physical Review A 85(5) (2012).
  * [48] P. Hauke and L. Tagliacozzo, _Spread of correlations in long-range interacting quantum systems_ , Physical Review Letters 111(20) (2013).
  * [49] S. S. Natu and E. J. Mueller, _Dynamics of correlations in a dilute bose gas following an interaction quench_ , Physical Review A 87(5) (2013).
  * [50] L. Cevolani, G. Carleo, and L. Sanchez-Palencia, _Protected quasilocality in quantum systems with long-range interactions_ , Physical Review A 92(4) (2015).
  * [51] L. Cevolani, G. Carleo, and L. Sanchez-Palencia, _Spreading of correlations in exactly solvable quantum models with long-range interactions in arbitrary dimensions_ , New Journal of Physics 18(9), 093002 (2016).
  * [52] A. S. Buyskikh, M. Fagotti, J. Schachenmayer, F. Essler, and A. J. Daley, _Entanglement growth and correlation spreading with variable-range interactions in spin and fermionic tunneling models_ , Physical Review A 93(5) (2016).
  * [53] I. Frérot, P. Naldesi, and T. Roscilde, _Multispeed prethermalization in quantum spin models with power-law decaying interactions_ , Physical Review Letters 120(5) (2018).
  * [54] T. E. Lee and C.-K. Chan, _Heralded magnetism in non-hermitian atomic systems_ , Phys. Rev. X 4, 041001 (2014).
  * [55] T. E. Lee, F. Reiter, and N. Moiseyev, _Entanglement and spin squeezing in non-hermitian phase transitions_ , Physical Review Letters 113(25) (2014).
  * [56] J. M. Hickey, S. Genway, I. Lesanovsky, and J. P. Garrahan, _Time-integrated observables as order parameters for full counting statistics transitions in closed quantum systems_ , Phys. Rev. B 87, 184303 (2013).
  * [57] X. Z. Zhang and Z. Song, _Non-hermitian anisotropic $xy$ model with intrinsic rotation-time-reversal symmetry_, Phys. Rev. A 87, 012114 (2013).
  * [58] E. H. Lieb and D. W. Robinson, _The finite group velocity of quantum spin systems_ , Comm. Math. Phys. 28, 251 (1972).