# Symmetries, conservation laws and entanglement in non-Hermitian fermionic lattices

Rafael Diogo Soares, Youenn Le Gal, Chun Leung, Dganit Meidan,

Alessandro Romito, Marco Schirò

► To cite this version:

Rafael Diogo Soares, Youenn Le Gal, Chun Leung, Dganit Meidan, Alessandro Romito, et al. Symmetries, conservation laws and entanglement in non-Hermitian fermionic lattices. SciPost Physics, 2025, 19 (4), pp.094. 10.21468/SciPostPhys.19.4.094. [https://arxiv.org/abs/2009.04468] (https://arxiv.org/abs/2009.04468)

HAL Id: hal-05408456

https://hal.science/hal-05408456v1

Submitted on 10 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Symmetries, conservation laws and entanglement in non-Hermitian fermionic lattices

Rafael D. Soares¹,², Youenn Le Gal¹, Chun Y. Leung³, Dganit Meidan⁴,⁵, Alessandro Romito³ and Marco Schirò¹

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

2. Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany

3. Department of Physics, Lancaster University, United Kingdom

4. Department of Physics, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel

5. SPEC, CEA, CNRS, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

## Abstract

Non-Hermitian quantum many-body systems feature steady-state entanglement transitions driven by the competition between unitary dynamics and dissipation. In this work, we reveal the fundamental role of conservation laws in shaping this competition. Focusing on translation-invariant non-interacting fermionic models with U(1) symmetry, we present a theoretical framework to understand the structure of the steady-state of these models and their entanglement content based on two ingredients: the nature of the spectrum of the non-Hermitian Hamiltonian and the constraints imposed on the steady-state single-particle occupation by the conserved quantities. These emerge from an interplay between Hamiltonian symmetries and initial state, due to the non-linearity of measurement back-action. For models with complex energy spectrum, we show that the steady state is obtained by filling single-particle right eigenstates with the largest imaginary part of the eigenvalue. As a result, one can have partially filled or fully filled bands in the steady-state, leading to an entanglement entropy undergoing a filling-driven transition between critical sub-volume scaling and area-law, similar to ground-state problems. Conversely, when the spectrum is fully real, we provide evidence that local observables can be captured using a diagonal ensemble, and the entanglement entropy exhibits a volume-law scaling independently on the initial state, akin to unitary dynamics. We illustrate these principles in the Hatano-Nelson model with periodic boundary conditions and the non-Hermitian Su-Schrieffer-Heeger model, uncovering a rich interplay between the single-particle spectrum and conservation laws in determining the steady-state structure and the entanglement transitions. These conclusions are supported by exact analytical calculations and numerical calculations relying on the Faber polynomial method.

<figure><img src="image_2.png" /></figure>

Copyright R. Diogo Soares et al.

This work is licensed under the Creative Commons

Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 2025-05-02

Accepted 2025-08-20

Published 2025-10-13

doi:10.21468/SciPostPhys.19.4.094

<figure><img src="image_3.png" /></figure>

Check for
updates

## Contents
<table><tbody><tr><td>2</td><td>Model, equations of motion and conservation laws</td><td>4</td></tr><tr><td>2.1</td><td>Impact of the initial state on conservation laws</td><td>5</td></tr><tr><td>3</td><td>Long-time limit, steady-state occupation and entanglement entropy</td><td>7</td></tr><tr><td>3.1</td><td>Steady-state structure: Complex spectrum</td><td>7</td></tr><tr><td>3.2</td><td>PT symmetric Hamiltonians</td><td>10</td></tr><tr><td>3.3</td><td>Entanglement entropy</td><td>11</td></tr><tr><td>4</td><td>Application: Hatano-Nelson model</td><td>13</td></tr><tr><td>5</td><td>Application: The non-Hermitian SSH Hamiltonian</td><td>15</td></tr><tr><td>5.1</td><td>The non-Hermitian SSH model</td><td>15</td></tr><tr><td>5.2</td><td>Steady-state single-particle distribution</td><td>17</td></tr><tr><td>5.3</td><td>Entanglement entropy</td><td>19</td></tr><tr><td>5.3.1</td><td>Class A</td><td>21</td></tr><tr><td>5.3.2</td><td>Class B</td><td>22</td></tr><tr><td>5.3.3</td><td>Class C</td><td>24</td></tr><tr><td>6</td><td>Conclusions</td><td>26</td></tr><tr><td>A</td><td>Numerical simulation of non-Hermitian lattices</td><td>28</td></tr><tr><td>B</td><td>Diagonal ensemble</td><td>29</td></tr><tr><td>B.1</td><td>Testing the diagonal ensemble</td><td>31</td></tr><tr><td>C</td><td>Analytical calculation of the entanglement entropy</td><td>31</td></tr><tr><td>C.1</td><td>The Hatano-Nelson model</td><td>32</td></tr><tr><td>C.2</td><td>The SSH chain</td><td>32</td></tr><tr><td>C.2.1</td><td>Class A</td><td>33</td></tr><tr><td>C.2.2</td><td>Class B</td><td>35</td></tr><tr><td>C.2.3</td><td>Class C</td><td>36</td></tr><tr><td>D</td><td>Scaling with subsystem size in the volume law phase</td><td>36</td></tr><tr><td></td><td>References</td><td>36</td></tr></tbody></table>

# 1 Introduction

Recent years have seen major progresses in the understanding of quantum dynamics in closed many-body systems [1, 2]. Here, the unitarity of time evolution puts strong constraints on the way correlations and entanglement spread through the system [3, 4], or local observables reach a stationary state and possibly approach thermal equilibrium at long-times [5].

A novel frontier of many-body quantum dynamics concerns the role of dissipation due to coupling to external environments, leading to a non-unitary time evolution [6]. An example of non-unitary dynamics which has attracted widespread interest is provided by non-Hermitian Hamiltonians [7], which arise naturally in the description of open quantum systems as no-click limit of a monitored evolution [8–10] and in the context of dynamics of systems with loss and gain [11–13].
Non-Hermitian quantum many-body systems have been the focus of several investigations, addressing their peculiar spectral and topological properties [14–16], correlation, and entanglement patterns [17–19]. An open question concerns the long-time limit of the system under non-Hermitian dynamics, namely, whether a generalized thermalization can be expected in these systems, and if so, what are the principles underpinning this process [20, 21]. Entanglement spreading has been also discussed in the context of non-Hermitian systems, leading to a wealth of entanglement transitions in different models, due to the skin effect [22, 23] or to spectral transitions [23–30]. However, a full understanding of the structure of entanglement of non-Hermitian systems is still lacking. Integrable or exactly solvable models have played a major role in driving the understanding of (generalized) thermalization in closed systems [31] and recently in open quantum systems described by the Lindblad master equation [32, 33]. Furthermore, they have been pivotal in the understanding of entanglement dynamics in unitary and monitored systems. It is therefore natural to investigate this class of systems in the non-Hermitian case.

In this work, we focus on non-interacting, translation-invariant, fermionic non-Hermitian models, augmented with a global U(1) symmetry. We first highlight a surprising consequence of the non-linearity of non-Hermitian dynamics, namely that the presence or absence of a conserved quantity associated to a symmetry depends nontrivially on the initial state. We then argue that the interplay between symmetries, conservation laws and spectral properties of the non-Hermitian Hamiltonian allow to completely specify the structure of the steady-state at long-times. Since for Gaussian systems the entanglement entropy is fully determined by the single-particle correlation matrix, this allows also to completely characterize the entanglement content of the steady-state and the associated entanglement transitions.

In particular, we show that the steady-state of a non-Hermitian free fermionic system with complex spectrum is obtained by filling single-particle orbitals with the highest imaginary part of the eigenvalue, compatibly with the constraints imposed by the symmetries of the non-Hermitian Hamiltonian and of the initial state. As a result, the scaling of the entanglement entropy in the steady-state obeys a logarithmic law or area law depending on whether the band of imaginary eigenvalues is partially or fully filled. A striking consequence of this statement is that non-Hermitian fermionic systems support entanglement transitions, which can be tuned by the filling of the system. This feature, which has no counterpart in unitary dynamics, resembles the behavior of entanglement in the ground-states of local Hamiltonians.

For a non-Hermitian Hamiltonian with purely real spectrum, such as in the presence of *PT*-symmetry [34] (or pseudo-hermiticity), we argue that the steady-state of local observables can be constructed from a diagonal ensemble similar to the one used for unitary dynamics. As a result, the entanglement entropy in the steady-state exhibits a volume-law scaling. Finally, if the spectrum is mixed with purely real and purely imaginary eigenvalues, such as in presence of *PT*-symmetry breaking, the steady-state depends nontrivially on the filling and symmetries of the initial state, and so, the entanglement entropy can display either a volume-law or logarithmic law scaling.

We test our theoretical scenario on two prototypical examples of non-Hermitian lattice models: the Hatano-Nelson [35, 36] and the Su-Heeger-Schrieffer model [27]. Our results and phase diagrams for the steady-state entanglement entropy are either obtained analytically or numerically using the Faber polynomial method [37], perfectly confirming our theoretical predictions.

The manuscript is organized as follows. In Sec. 2 we introduce our general model of non-Hermitian free fermionic systems and highlight the role of the initial state in the definition of conserved quantities for non-Hermitian systems. This will lead us to introduce three classes of initial states that will be used throughout the manuscript to discuss the results. In Sec. 3 we present our main general result concerning the structure of the steady-state occupation in our
models, both in presence of complex spectrum and of PT-symmetry, and its implication for the scaling of entanglement entropy. In Sec. 4 we present a first application of our framework to the Hatano-Nelson model and its entanglement entropy scaling in the steady-state. In Sec. 5 we discuss the non-Hermitian SSH model with its entanglement transitions and highlight the role of filling and initial state in controlling the structure of the phase diagram. Finally, Sec. 6 is devoted to the conclusions. Several appendices complete this work.

## 2 Model, equations of motion and conservation laws

In this work, we are interested in the dynamics generated by time-independent quadratic fermionic non-Hermitian Hamiltonians, that we can generally write as

$$ \mathcal{H}_{\text{nH}} = \sum_{ij} \sum_{\alpha\beta} c_{i\alpha}^{\dagger} h_{ij}^{\alpha\beta} c_{j\beta}, \qquad (1) $$

where i, j are sites of a d-dimensional lattice, α, β additional internal quantum numbers that can include orbital, sub-lattice or spin degrees of freedom and $h_{ij}^{\alpha\beta}$ is a non-Hermitian matrix that we can write as

$$ h_{ij}^{\alpha\beta} = h_{0,ij}^{\alpha\beta} + i\mathcal{V}_{ij}^{\alpha\beta}, \qquad (2) $$

with $h_{0,ij}^{\alpha\beta} = (h_{0,ji}^{\beta\alpha})^*$ and $\mathcal{V}_{ij}^{\alpha\beta} = (\mathcal{V}_{ji}^{\beta\alpha})^*$ respectively the real and imaginary part of the single-particle Hamiltonian. In the following, we will consider a one-dimensional lattice for simplicity, even though many of our conclusions apply more generally to non-interacting non-Hermitian fermions in higher dimensions. For our purposes, the non-Hermitianity of $\mathcal{H}_{\text{nH}}$ arises from dissipative processes such as gain and loss or more generally back-action terms associated to quantum measurements and postselection [10,38], see below for specific examples. For this reason, in the following, we will refer to the energy scale related to non-Hermiticity as dissipation or measurement back-action.

We assume the Hamiltonian to be translational-invariant, $h_{i,j}^{\alpha\beta} \equiv h_r^{\alpha\beta}$ with $r = |i - j|$, so that we can rewrite it in momentum space as

$$ \mathcal{H}_{\text{nH}} = \sum_k \sum_{\alpha\beta} c_{k\alpha}^{\dagger} h_k^{\alpha\beta} c_{k\beta}, \qquad (3) $$

with $h_k^{\alpha\beta} = \sum_r e^{ikr} h_r^{\alpha\beta}$ the single-particle Hamiltonian in momentum space. This Hamiltonian can be decomposed into its Hermitian part, $h_{0,k}$, which is responsible for the unitary evolution in the absence of dissipative processes, and the dissipative contribution, $i\mathcal{V}_k$,

$$ h_k^{\alpha\beta} = h_{0,k}^{\alpha\beta} + i\mathcal{V}_k^{\alpha\beta}, \qquad (4) $$

with $h_{0,k} = h_{0,k}^{\dagger}$. We assume that this Hamiltonian can be diagonalised with eigenvalues $\epsilon_{k,a} = \text{Re}\,\epsilon_{k,a} + i\text{Im}\,\epsilon_{k,a}$. This describes the complex energies of the non-Hermitian quasi-particles of the problem which have a finite decay rate $\gamma_{k,a} = \text{Im}\,\epsilon_{k,a}$.

The system evolution under $\mathcal{H}_{\text{nH}}$ reads

$$ |\Psi(t)\rangle = \frac{e^{-i\mathcal{H}_{\text{nH}}t} |\Psi(0)\rangle}{\|e^{-i\mathcal{H}_{\text{nH}}t} |\Psi(0)\rangle\|}, \qquad (5) $$

where the denominator properly normalizes the state, compensating the non-conservation of the state's norm due to the absence of unitarity in non-Hermitian systems. The equation above can be obtained as a no-click limit of the quantum jump dynamics of monitored systems [6].
Normalizing the state after the evolution is equivalent to solving a non-linear Schrödinger equation

$$d|\Psi(t)\rangle = -i\mathcal{H}_{\text{nH}}dt|\Psi(t)\rangle - i\frac{dt}{2}\langle\mathcal{H}_{\text{nH}} - \mathcal{H}_{\text{nH}}^{\dagger}\rangle_{t}|\Psi(t)\rangle, \quad (6)$$

where the last term, proportional to the imaginary part of the Hamiltonian, depends non-linearly on the state itself.

In this manuscript, we focus on the dynamics described by Eq. (5), (6) and discuss the steady-state properties of the system and its entanglement content. The non-Hermitian nature of the problem has direct consequences on the Heisenberg equations of motion of any operator, $\mathcal{O}$, evolving under $\mathcal{H}_{\text{nH}}$. Indeed, the equation of motion is given by

$$\partial_t \langle \mathcal{O} \rangle = i \langle \mathcal{H}_{\text{nH}}^{\dagger} \mathcal{O} - \mathcal{O} \mathcal{H}_{\text{nH}} \rangle - i \langle \mathcal{H}_{\text{nH}}^{\dagger} - \mathcal{H}_{\text{nH}} \rangle \langle \mathcal{O} \rangle. \quad (7)$$

The last term arises explicitly from the measurement-back action (i.e. from the non-Hermitian part of the Hamiltonian) and vanishes in Hermitian systems, thereby recovering the standard Heisenberg equation of motion. This term induces a non-linear, state dependent form to the equations of motion, even for a quadratic non-interacting systems such as in Eq. (1). To appreciate this point, it is useful to write down explicitly the equations of motion for the correlation matrix $G_{ij}^{\alpha\beta} = \langle c_{ia}^{\dagger} c_{j\beta} \rangle$, which for a fermionic gaussian problem whose dynamics conserves particle number encodes all relevant information. In momentum space $G_{kq}^{\alpha\beta} = \langle c_{ka}^{\dagger} c_{q\beta} \rangle$, one can show using Wick's theorem and the decomposition in Eq. (4) that the diagonal part satisfies the following equation of motion,

$$\partial_t G_{kk}^{\alpha\beta} = i \left( (h_k^*)^{\alpha\gamma} G_{kk}^{\gamma\beta} - G_{kk}^{\alpha\gamma} h_k^{\beta\gamma} \right) - 2 \sum_q G_{kq}^{\alpha\xi} (V_q^*)^{\xi\gamma} G_{qk}^{\gamma\beta}, \quad (8)$$

where repeated indices are summed. We note that the equation is non-linear due to the non-Hermitian part of the Hamiltonian. Furthermore, even for a translation-invariant Hamiltonian, the diagonal elements of the correlation matrix are coupled to the off-diagonal ones due to the non-Hermitian term in the Hamiltonian. This is not the case for a Hermitian system, as the last term is absent. This already suggests that symmetries and conserved quantities play a special role for non-Hermitian systems. Finally, for what concerns the entanglement content of the steady-state, we characterize it using the von Neumann entanglement entropy, defined as [39, 40]

$$S(t) = -\text{tr}_A[\rho^A(t) \ln \rho^A(t)], \quad (9)$$

where we have introduced a bipartition $A \cup B$ in the system, with the reduced density matrix $\rho^A(t) = \text{tr}_B|\Psi(t)\rangle\langle\Psi(t)|$. The von Neumann entropy correctly measures entanglement, as the system is initially in a pure state and purity is conserved throughout the time evolution, as seen by Eq. (5). As we focus on quadratic models, the entanglement entropy can be directly obtained from the correlation matrix [41].

## 2.1 Impact of the initial state on conservation laws

We now discuss the definition of a conserved quantity for non-Hermitian dynamics and the influence of the initial state, which play a crucial role in the dynamical evolution.

In a Hermitian system with Hamiltonian $\mathcal{H}$ a time-independent observable $\mathcal{Q}$ is conserved if $[\mathcal{H}, \mathcal{Q}] = 0$, as directly given by the Heisenberg equation of motion: $\partial_t \langle \mathcal{Q} \rangle = i[\mathcal{H}, \mathcal{Q}]$. In this case, the expectation value of $\mathcal{Q}$ is fixed by the initial state, $\langle \Psi_0 | \mathcal{Q} | \Psi_0 \rangle$. All symmetries of the Hamiltonian correspond to a conserved quantity, and vice versa, as encapsulated by the Noether theorem, irrespectively of whether the initial state respects the symmetry.
Non-Hermitian systems, on the other hand, exhibit peculiar behavior due to their non-linear state-dependent evolution equation. In particular, the conservation laws in non-Hermitian systems are not solely determined by the symmetries of the Hamiltonian. As in other settings of open quantum systems, such as the Lindblad master equation, the existence of a (weak) symmetry in the dynamics does not imply the presence of a conserved quantity [6]. For a non-Hermitian system, in addition to the symmetries, the initial state also plays a key role in determining the existence or not of a conserved quantity. Indeed, as one can see from Eq. (7) for non-Hermitian systems in order for the expectation value of an operator $Q$ to be conserved by the non-unitary evolution, two conditions are required, namely: (i) the operator has to commute with $\mathcal{H}_{nH}$, i.e., $[\mathcal{H}_{nH}, Q] = 0$ and (ii) the initial state must be an eigenstate of the operator $Q$, i.e. $Q|\Psi_0\rangle = q_0|\Psi_0\rangle$. These two conditions together guarantee that the state remains an eigenstate of the conserved quantity throughout the evolution. In fact, the equation of motion is identically zero because the following identity holds $\langle\mathcal{H}_{nH}^\dagger Q - Q\mathcal{H}_{nH}\rangle = \langle\mathcal{H}_{nH}^\dagger - \mathcal{H}_{nH}\rangle\langle Q\rangle$. In other words, we can write

$$[\mathcal{H}_{nH}, Q] = 0, \quad \text{and} \quad Q|\Psi_0\rangle = q_0|\Psi_0\rangle \Rightarrow \partial_t\langle Q\rangle = 0. \qquad (10)$$

To appreciate the role of the initial state in determining the conservation laws of our system, it is useful to note that our non-Hermitian Hamiltonian in Eq. (1) admits both a U(1) symmetry, $c_{i\alpha} \to e^{i\phi}c_{i\alpha}$ and a symmetry under spatial translation. As a result, we can identify two operators that commute with $\mathcal{H}_{nH}$: the number of fermions at momentum $k$,

$$\hat{n}_k = \sum_{\alpha} c_{k,\alpha}^{\dagger} c_{k,\alpha}, \qquad (11)$$

as well as, the total number of particles,

$$\hat{N} = \sum_k \hat{n}_k. \qquad (12)$$

We note that, conversely, the occupation of the eigenmodes of $\mathcal{H}_{nH}$ are not in general Hermitian operators and, therefore, not of direct utility for our purpose here. To assess whether to these operators we can associate genuine conserved quantities, we will have to consider three distinct classes of initial states $|\Psi_0\rangle$.

* class A: $|\Psi_0\rangle$ is an eigenstate of $\hat{n}_k$, for all $k$. In general, this condition gives a Slater determinant of the form

$$|\Psi_0\rangle = \prod_{k \in \text{occupied}} \eta_k^{\dagger} |\text{vac}\rangle, \qquad (13)$$

with $\eta_k = \sum_{\beta} U_{k,\beta} c_{k,\beta}$ with $U_k$ satisfying $\sum_{\beta} |U_{k,\beta}|^2 = 1$. For initial states in this class, the expectation value of $\hat{n}_k$ is conserved even under the non-Hermitian evolution. In other words, states in this class conserve the total particle number and are translationally invariant.

* class B: $|\Psi_0\rangle$ is an eigenstate of $\hat{N}$ but not an eigenstate of $\hat{n}_k$, hence the expectation value of $\hat{n}_k$ admits a nontrivial dynamics, while the total particle number is a genuine conserved quantity of the non-Hermitian dynamics. An example of initial state in this class, that we will use later on in the manuscript, read:

$$|\Psi_0\rangle = \prod_{l=0}^{N-1} c_{\sigma(l)}^{\dagger} |\text{vac}\rangle, \qquad (14)$$

where $\sigma$ is a permutation of $\{1, ..., L\}$ with $L$ the size of the chain.
* class C: $|\Psi_0\rangle$ is neither an eigenstate of $\hat{n}_k$ nor of the particle number $\hat{N}$, therefore, the associated expectation values are not conserved during the dynamics. Typically, such a state can be an eigenstate of Bogoliubov-de-Gennes Hamiltonians and can be written in the form

$$|\Psi_0\rangle = \prod_{l=0}^{L/2-1} [u_l + v_l c_{2l}^{\dagger} c_{2l+1}^{\dagger}] |\text{vac}\rangle, \quad (15)$$

where $|u_l|^2 + |v_l|^2 = 1$ to assure the correct normalization. One useful quantity is then the average number of particles, $\langle \hat{N} \rangle$, and the fluctuations of the total particle number $\langle \delta \hat{N}^2 \rangle = \langle \hat{N}^2 \rangle - \langle \hat{N} \rangle^2$, which for the initial state written above are given by,

$$\langle \hat{N} \rangle = 2 \sum_{l=0}^{L/2-1} |v_l|^2, \quad \langle \delta \hat{N}^2 \rangle = 4 \sum_{l=0}^{L/2-1} |u_l|^2 |v_l|^2, \quad (16)$$

respectively.

In the following, we will discuss the general features of the long-time dynamics of non-Hermitian free fermionic systems, starting from initial states in these three classes.

## 3 Long-time limit, steady-state occupation and entanglement entropy

We now present our main result concerning the long-time limit of a free non-Hermitian fermionic system and its entanglement structure, which we will then discuss in detail for specific examples in the remainder of the manuscript (see Sections 4 and 5). Specifically, here we aim at characterizing the steady-state of the system in terms of the occupation of single-particle states $\langle \hat{n}_k \rangle$, whose dynamics as discussed above depends strongly on the initial condition. This will allow us to understand the role that the filling of the system can have and the consequences we can observe on the scaling of entanglement entropy.

### 3.1 Steady-state structure: Complex spectrum

To make general statements about the steady-state structure, we assume the non-Hermitian Hamiltonian to be diagonalizable in a bi-orthogonal basis,

$$\mathcal{H}_{\text{nH}} = \sum_{n} (E_n + i\Gamma_n) |\Psi_n^R\rangle\langle\Psi_n^L|, \quad (17)$$

where $|\Psi_n^{R/L}\rangle$ is the many-body right/left $n^{th}$ eigenstate and $(E_n + i\Gamma_n)$ the associated eigenvalue that we take here to be complex, i.e. we will assume $\Gamma_n$ to be non-zero for any $n$. In general, $\Gamma_n$ can be either positive (describing growing modes) or negative (describing decaying or damped modes). However, due to the normalization in the non-Hermitian evolution, the spectrum is defined up to a global shift, so its overall sign is irrelevant: only differences among eigenvalues matter.

The evolution of the state (given in Eq. (5)) can be written using this bi-orthogonal basis by inserting the identity $1 = \sum_n |\Psi_n^R\rangle\langle\Psi_n^L|$ between the initial state and the time evolution operator,

$$|\Psi(t)\rangle = \frac{e^{\Gamma_{ss}t}}{\mathcal{N}(t)} \left( e^{-iE_{ss}t} \langle\Psi_{ss}^L|\Psi(0)\rangle |\Psi_{ss}^R\rangle + \sum_{n \neq ss} e^{(\Gamma_n - \Gamma_{ss})t} e^{-iE_n t} \langle\Psi_n^L|\Psi(0)\rangle |\Psi_n^R\rangle \right), \quad (18)$$
<figure><img src="image_4.png" /><figcaption>Figure 1: Sketch of the steady-state occupation of single-particle energy levels by class. For simplicity, we consider the case of a two-band single-particle Hamiltonian. The black balls represent the initial distribution. This should be interpreted as the expectation value of $\hat{n}_k$ for class B and class C.</figcaption></figure>

where $N(t)$ properly normalizes the state, and we identified with the index $n = ss$ the pair of right and left eigenstates $|\Psi_n^{R/L}\rangle$ that have the highest $\Gamma_n$, the imaginary part of the eigenvalue, among those with a non-zero overlap with the initial state, $\langle\Psi_{ss}^L|\Psi(0)\rangle \neq 0$. For simplicity, we have assumed that the state with the highest $\Gamma_n$ is non-degenerate. In the long time limit, we obtain

$$ \lim_{t \to +\infty} |\Psi(t)\rangle = e^{-iE_{ss}t} |\Psi_{ss}^R\rangle. \qquad (19) $$

The interpretation of this statement is simple: the state which is populated at long times is the one with highest imaginary part of the eigenvalue, corresponding to the fastest growth rate (if $\Gamma_{ss} > 0$) or slowest decay rate (if $\Gamma_{ss} < 0$). We can say more about the structure of the steady-state by taking into account the symmetries of the non-Hermitian Hamiltonian. Indeed, both the left and right eigenvectors are eigenstates of the Hamiltonian's symmetries. For the non-Hermitian quadratic Hamiltonian considered here, this implies that each eigenstate is labeled by the occupation number, $n_k$, for each momentum $k$, and the total particle number, $N$, namely we have that

$$ \begin{aligned} \forall k, \hat{n}_k |\Psi_m^{R/L}\rangle &= (n_k)_m |\Psi_m^{R/L}\rangle, \\ \hat{N} |\Psi_m^{R/L}\rangle &= N_m |\Psi_m^{R/L}\rangle, \end{aligned} \qquad (20) $$

where $(n_k)_m$ and $N_m$ are the eigenvalues corresponding to $m^{th}$ eigenstate. Finally, the detailed structure of the steady-state differs for the various classes of initial states considered in this work, as we now discuss.

* An initial state in class A is also an eigenstate of the occupation number $\hat{n}_k$, for each momentum $k$ with $\hat{n}_k|\Psi_0\rangle = (n_k)_{t=0}|\Psi_0\rangle$, and the total particle number, $\hat{N}$ with $\hat{N}|\Psi_0\rangle = N_{t=0}|\Psi_0\rangle$. Consequently, as imposed by Eq. (7), throughout the dynamics, the average occupation for each momentum, $\langle\hat{n}_k\rangle$ remains conserved and fixed at its initial value, $(n_k)_{t=0}$. As a result, the steady-state in this case corresponds to the right eigen-state of the Hamiltonian characterized by $(n_k)_{t=0}$ quantum numbers, with the highest imaginary eigenvalue, $\Gamma_{ss}$. This component can be written using the single-particle energies as follows

$$ \Gamma_{ss} = \sum_k \sum_{\alpha=0}^{(n_k)_{t=0}-1} \mathrm{Im} \, \varepsilon_{k,\alpha}, \qquad (21) $$
where $(n_k)_{t=0}$ is the eigenvalue of the $\hat{n}_k$ operator in the initial state, $\alpha$ corresponds to a band index of the single-particle spectra, and we have assumed the ordering convention $\mathrm{Im} \, \varepsilon_{k,0} > \mathrm{Im} \, \varepsilon_{k,1} > \cdots$. This is represented in the scheme of Fig. 1. The single-particle states with the largest imaginary part of the eigenvalue are filled, preserving the total particle number and initial momentum occupation. As seen in the example of Fig. 1, the number of red balls matches the number of black ones, and both occupy the same momentum values.

* For an initial state in class B the situation is distinct. In this case, the initial state only has a well-defined particle number. The expectation value $\langle \hat{n}_k \rangle$ is no longer conserved under the evolution, as the initial state as a non-zero overlap with eigenstates that have different values of $(n_k)_m$. In this case, the steady state is the right eigenstate that maximizes the imaginary part of the eigenvalue, while maintaining the same particle number as the initial state, $\langle \hat{N} \rangle_{t=0}$.

The steady state is then obtained by filling the single-particle states with the largest imaginary parts of the eigenvalue, subject to the total particle number constraint. Given this,

$$ \Gamma_{ss} = \sum_{m=0}^{\langle \hat{N} \rangle_{t=0}-1} \mathrm{Im} \, \varepsilon_m, \qquad (22) $$

where $\langle N \rangle_0$ is the eigenvalue of total particle number operator correspondent to the initial state, $m = (k, \alpha)$ and we assume the ordering convention $\mathrm{Im} \, \varepsilon_0 > \mathrm{Im} \, \varepsilon_1 > \cdots$. This last equation also tells us the value of the distribution $(n_k)_{ss}$ in the steady-state. As seen in the scheme of Fig. 1, the number of occupied final states (depicted in green) matches the number of initial occupied states (represented in black), but they do not occupy the same momentum values.

* In class C, the initial states are neither eigenstates of the occupation in momentum nor of the total number operator. Generically, the initial state overlaps with all particle number sectors. Thus, in general, the steady state of the system corresponds to the right eigenstate with the highest imaginary part of the eigenvalue in the full many-body spectrum. This corresponds to filling all single-particle states that have a positive imaginary eigenvalue, corresponding to amplifying modes

$$ \Gamma_{ss} = \sum_m \mathrm{Im} \, \varepsilon_m, \quad \mathrm{Im} \, \varepsilon_m \ge 0, \qquad (23) $$

where $m = (k, \alpha)$ and we assume that the imaginary part of the single-particle spectrum has a zero center of mass, meaning the imaginary part of the sum of all single-particle eigenvalues is zero. In the example of Fig. 1, the number of blue balls does not match that of black ones.

We can summarize these considerations and express the steady-state occupation of single-particle modes $n_k$ as

$$ (n_k)_{ss} = \begin{cases} (n_k)_{t=0}, & \text{class A,} \\ \Theta(\mathrm{Im} \, \varepsilon_{k,\alpha} - \mu_{\mathrm{eff}}), & \text{class B,} \\ \Theta(\mathrm{Im} \, \varepsilon_{k,\alpha}), & \text{class C,} \end{cases} \qquad (24) $$

where $(n_k)_{ss}$ is the eigenvalue of the operator $\hat{n}_k$ in the steady state, $\Theta(\cdot)$ denotes the Heaviside step function, and $\mu_{\mathrm{eff}}$ is an effective chemical potential ensuring that the total particle number remains equal to that of the initial state, i.e. $\sum_k (n_k)_{ss} = N_0$. As before, we assume that the imaginary part of the single-particle spectrum has a zero center of mass.
We note that the steady state given by Eq. (18) possesses the same symmetries as the considered non-Hermitian Hamiltonian, although the initial state did not. Interestingly, this implies that time evolution allows for the restoration of the symmetries, even though they were initially explicitly broken by the state. For example, an initial state in class C is a coherent superposition of states with distinct values of the total particle number, but the steady state is an eigenstate of the total particle number. We will return to this point when discussing specific examples.

We conclude by noticing the formal analogy between what we discussed so far for non-Hermitian evolution concerning filling of modes with highest imaginary eigenvalue, and the long-time limit of a state evolved in imaginary time with a Hermitian Hamiltonian. In this case, the initial state converges to the ground-state of the Hamiltonian within the sector that has the same quantum numbers as the initial state and has the lowest (real) energy. This can also be viewed as real-time evolution under the non-Hermitian Hamiltonian $\mathcal{H}_{\text{nH}} = -i\hat{H}$, where $\hat{H}$ is a Hermitian Hamiltonian. Thus, minimizing the real energy in imaginary-time evolution corresponds to maximizing the imaginary energy in real-time evolution.

## 3.2 PT symmetric Hamiltonians

In the previous section we have assumed a complex energy spectrum. However, in non-Hermitian physics an important class of problems features PT-symmetry [34]. In this case, the non-Hermitian Hamiltonian $\mathcal{H}_{\text{nH}}$ commutes with the PT operator, the conjugation of the parity (P) and time-reversal (T) operators, that is, $[\mathcal{H}_{\text{nH}}, PT] = 0$. This symmetry ensures that the spectrum of the Hamiltonian is real whenever the eigenstates are also eigenvectors of the PT operator [34]. Depending on the values of the system's parameters, this condition may or may not hold. When it does, the system is said to be in a PT-symmetric phase. Conversely, if the condition is not satisfied, the spectrum becomes complex and the system is in a PT-broken phase, as the state spontaneously breaks the PT symmetry. The breaking of PT symmetry can also occur in stages, with only part of the spectrum developing imaginary eigenvalues, before the full spectrum becomes complex, as we will discuss in an example later in the manuscript.

In the PT-symmetric phase, the steady-state cannot be described by Eq. (18) as all eigenstates are real, and so, in principle, none is a priori more enhanced/damped than the others. Although the time evolution is still non-unitary due to the non-orthogonality between distinct right/left eigenstates, $\langle \Psi_n^R | \Psi_m^R \rangle \neq 0$, we can understand the structure of the steady-state using an analogy with the emergence of the diagonal ensemble under unitary dynamics in closed isolated quantum systems. As we show in the appendix B, the expectation value of local observables in the steady-state can be approximated by a diagonal ensemble of the form:

$$ \rho_{\text{DE}} = \frac{1}{Z_{\text{DE}}} \sum_{n} |\langle \Psi(0) | \Psi_{n}^{L} \rangle|^2 |\Psi_{n}^{R} \rangle \langle \Psi_{n}^{R} |, \qquad (25) $$

where $Z_{\text{DE}} = \sum_{n} |\langle \Psi(0) | \Psi_{n}^{L} \rangle|^2$ is the normalization factor. We note that the diagonal ensemble explicitly depends on the initial state of the system. The density matrix of Eq. (25) is constructed only with the right eigenstates that have the same quantum numbers as the initial state, since otherwise the overlap $\langle \Psi(0) | \Psi_{n}^{L} \rangle$ is zero. We stress that, as its unitary case counterpart, Eq. (25) should only be understood as valid when evaluating the long-time expectation value of sufficiently local operators. The full wave function remains in a pure state under non-Hermitian evolution, yet local observables can reach a steady-state via many-body dephasing (see Appendix B). Finally, we mention that a similar construction of the diagonal ensemble could be useful for other types of pseudo-Hermitian Hamiltonians (among which PT-symmetric Hamiltonians represent one possible example).
If instead the system breaks spontaneously PT symmetry and develops imaginary eigenvalues (i.e. in the PT-broken phase), the structure of the steady-state can be still understood using the argument in previous section, at least when all the spectrum becomes complex. In this scenario, the steady state corresponds to the least damped mode subject to the activated conservation laws. We emphasize that for certain values of the Hamiltonian parameters, exceptional points may arise, preventing us from diagonalizing the full Hamiltonian. In this case, the construction of Eq. (18) cannot be made. Typically, these exceptional points appear in the spectrum along with spontaneous symmetry breaking of PT symmetry [24]. We discuss this pathological case later in the manuscript.

### 3.3 Entanglement entropy

We now use the knowledge of the steady-state occupation discussed above to understand the structure of the entanglement entropy in the steady state, $S(+\infty) = \lim_{t \to +\infty} S(t)$, where $S(t)$ is defined in Eq. (9). As in previous sections, we focus here on one-dimensional fermions, although our conclusions can be easily extended to higher dimensions. Moreover, since we are considering Gaussian fermionic states, knowledge of the single-particle correlation matrix is sufficient to obtain the entanglement entropy, at least numerically. Remarkably, we argue that one can obtain exact analytical predictions for the scaling of the entanglement entropy by knowing the steady-state occupation. Quite generally we can write the steady-state entanglement entropy of a bi-partition of linear size $l$ in the form

$$S(+\infty) = a_1 l + a_2 \ln l + \mathcal{O}(1). \qquad (26)$$

The coefficient $a_1$ corresponds to the extensive volume-law contribution to the entanglement entropy, $a_2$ controls the sub-extensive term characteristic of critical behavior, and finally we have omitted the constant term corresponding to the area-law scaling. While Eq. (26) encompasses all cases of interest for our models, we show in Sec. 5 and Appendix C how the coefficients $a_1, a_2$ can be computed analytically in certain cases.

To see how Eq. (26) comes about, let us start by discussing the case when the spectrum is complex and our discussion in Sec. 3.1 applies (see Eq. (18)). If this is the case, then the entanglement entropy in the steady state is fully determined by the distribution $\langle \hat{n}_k \rangle$ and the structure of the imaginary spectrum of the single-particle Hamiltonian. In this regime, the steady-state correlation matrix is diagonal in the momentum index, that is, $G_{kq}^{\alpha\beta} = 0$ for $k \neq q$. For a given distribution of $\langle \hat{n}_k \rangle$ and an imaginary single-particle spectrum, the steady state resembles a Slater-determinant characteristic of fermions with a partially or completely filled band of imaginary eigenvalues. Consequently, if the eigenstates of the non-Hermitian Hamiltonian can be analytically determined, the steady-state entanglement entropy can also be explicitly obtained for all classes of initial states using the Szegö theorem [3] and the Fisher-Hartwig conjecture [42–44]. When the steady state exhibits a fully occupied band of imaginary eigenvalues, and the final occupied band is separated by an energy gap from the next band in the imaginary part of the spectrum, the entanglement entropy obeys an area law. In contrast, when the spectrum is gapless, it exhibits a critical sub-volume logarithmic scaling. Generically, if it exhibits metallic-like behavior, i.e. a partially filled band, the system’s entanglement entropy will grow logarithmically with subsystem size.

This behavior mirrors the scaling of the entanglement entropy observed in the ground-states of models with short-range interactions and hopping [45, 46]. However, for the quadratic non-Hermitian Hamiltonians considered here, the structure of the imaginary part of the single-particle spectrum (i.e. whether it is gapped or gapless) alone does not determine the entanglement entropy scaling. Instead, the momentum-space occupation, which depends on the initial state and on the initial filling, also plays a crucial role.
<table><thead><tr><th>Spectrum</th><th>Steady State</th><th>Purely Real Energy Eigenstates</th><th>Complex Energy Band</th><th>Band Occupation</th><th>Entanglement scaling</th></tr></thead><tbody><tr><td rowspan="3">Complex</td><td rowspan="3">|Ψ<sup>R</sup><sub>ss</sub>⟩</td><td>N.A.</td><td>Fully filled</td><td></td><td>S ∝ O(1)</td></tr><tr><td>N.A.</td><td>Partially filled</td><td></td><td>S ∝ ln ℓ</td></tr><tr><td>N.A.</td><td>Fully filled</td><td></td><td>S ∝ ln ℓ</td></tr><tr><td rowspan="3">Mixed</td><td rowspan="3"></td><td>Empty</td><td>Partially filled</td><td></td><td>S ∝ ln ℓ</td></tr><tr><td>Occupied</td><td>Partially filled</td><td></td><td>S ∝ ℓ</td></tr><tr><td>Fully Occupied</td><td>Partially filled</td><td></td><td>S ∝ ln ℓ</td></tr><tr><td>Purely Real</td><td>ρ<sub>DE</sub></td><td>Occupied</td><td>N.A.</td><td>●●●●●</td><td>S ∝ ℓ</td></tr></tbody></table>

○ Empty mode ● Filled Mode ● Occupied Purely Real Mode ● Fully Occupied Purely Real Mode

<caption>Figure 2: Scaling of steady-state entanglement entropy, depending on the spectrum of the non-hermitian Hamiltonian and the structure of the steady state. In the *band* occupation column, we illustrate the case of a two-band model for simplicity.</caption>

Let us now consider the case of a non-Hermitian Hamiltonian with *PT* symmetry. If this symmetry is unbroken, all eigenstates share the same imaginary eigenvalue (zero), leading to a high degree of degeneracy. In these cases, Eq. (18) is no longer applicable. Instead, the steady state is typically a linear combination of the Hamiltonian's eigenstates, and the correlation matrix is generically not diagonal in momentum space. As a result, the Szegö theorem and the Fisher-Hartwig conjecture cannot be directly used. Based on our discussion of the steady-state occupation and the emergence of a diagonal ensemble, we can expect the entanglement entropy to satisfy a volume-law scaling, in analogy with closed systems. This can be checked explicitly if the initial state belongs to class A, when the correlation matrix remains diagonal in momentum space. In this case, the computation is feasible, provided that the equations of motion are solved exactly and the steady-state correlation matrix is determined, as we will show in Sec. 5 for the SSH model.

Finally, if the spectrum contains both purely real and purely imaginary eigenvalues, such as when the *PT* symmetry is broken, the entanglement entropy depends on the filling of the imaginary and real eigenstates, with possible transition between volume-law and subvolume (logarithmic) scaling. We will discuss a specific example of this case in Sec. 5.

In the table of Fig. 2, we present a summary of the possible steady-state entanglement scalings as a function of the distribution of single-particle states in the steady state. As we have explained, the specific distribution depends on both the spectrum of the non-Hermitian Hamiltonian and the conservation laws imposed by the initial state.
## 4 Application: Hatano-Nelson model

We start by analyzing the paradigmatic Hatano-Nelson model [35, 36] with periodic boundary conditions, a one-dimensional chain of spinless fermions with an asymmetric hopping, whose real-space Hamiltonian is given by

$$ \mathcal{H}_{\text{nH}} = -\frac{1}{2} \sum_{l=0}^{L-1} ([J + \gamma] c_l^\dagger c_{l+1} + [J - \gamma] c_{l+1}^\dagger c_l) = \sum_k \varepsilon_k c_k^\dagger c_k, \quad (27) $$

where $J$ is the coherent hopping strength to the first-neighbor sites, $\gamma \in \mathbb{R}$ induces an imbalance in the charge hopping and $c_l^\dagger(c_l)$ is the creation (annihilation) operator which creates (destroys) a fermion on site $l$. The Hatano-Nelson Hamiltonian is diagonalized by Fourier transforming the fermionic operators, and as such the Hamiltonian is diagonal in the momentum occupation $\hat{n}_k = c_k^\dagger c_k$ with a complex single-particle dispersion relation

$$ \varepsilon_k = -J \cos(ka) + i\gamma \sin(ka). \quad (28) $$

We see that taking $\gamma > 0$, the imaginary part of the spectrum is positive for $0 < k < \pi$ and negative otherwise.

We start considering an initial state in class A, corresponding to

$$ |\Psi_0\rangle = \prod_{k \in \Omega} c_k^\dagger |\text{vac}\rangle, \quad (29) $$

where $\Omega$ a collection of the values of momenta $k$. For this specific case, the time evolution is trivial, as the initial state is an eigenstate of the Hamiltonian, and so the state does not evolve in time

$$ |\psi(t)\rangle = e^{-i\text{Re}(E_k)t} |\psi_0\rangle, \quad (30) $$

where $E_k = \sum_{k \in \Omega} \varepsilon_k$. As expected, the state maintains the expectation value of $\langle \hat{n}_k \rangle$ for all $k$. Moreover, the entanglement entropy in the steady state is exactly the same as in the initial state.

On the other hand, the state evolves nontrivially for an initial state in either class B or C. Using the considerations discussed in the previous section and Eq. (18), the steady state for an initial state in class B is a Slater determinant obtained by filling the single-particle eigenstates with the largest imaginary part. For the single-particle dispersion relation in Eq. (28) and for $\gamma > 0$ this corresponds to filling the single-particle eigenstates around $k = \pi/2$,

$$ |\Psi_{ss}^B\rangle = \prod_{k \in [\frac{\pi}{2}-k_f, \frac{\pi}{2}+k_f]} c_k^\dagger |\text{vac}\rangle, \quad (31) $$

where $k_f$ is determined by the total particle number of the initial state. With this in mind, the distribution is given by

$$ (n_k)_{ss} = \begin{cases} 1, & k \in [\frac{\pi}{2} - k_f, \frac{\pi}{2} + k_f], \\ 0, & \text{otherwise.} \end{cases} \quad (32) $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Time evolution of the expectation value of $\langle \hat{n}_k \rangle$ distribution. In panel a), the initial state corresponds to a charge density wave, $|\Psi_0\rangle = |101 \cdots 101\rangle$, and so it is in class B (corresponding to an initial filling equal to 1/2), while in panel b) the initial state belongs to class C. Other parameters: $L = 128$ and $\gamma = 0.4J$.</figcaption></figure>

We show this in the panel a) of Fig. 3, where we plot the dynamics of the occupation number $\langle \hat{n}_k \rangle$ for increasing values of time. We see that for an initial state at half filling, $\langle \hat{n}_k(0) \rangle = 1/2$, the non-Hermitian dynamics fills up the states around $k = \pi/2$ in the long-time limit, while depleting states with negative imaginary part of the eigenvalues, corresponding to decaying modes, in such a way that the total number of particles remains preserved. For an initial state in class C, based on our general framework, we expect the steady state to be of the form

$$|\Psi_{ss}^C\rangle = \prod_{k \in [0, \pi]} c_k^\dagger |\text{vac}\rangle, \qquad (33)$$

corresponding to a Slater determinant with the highest imaginary eigenvalue. Since for initial states in class C there is no particle number conservation, this steady state corresponds to filling all single-particle states with positive energy. For $\gamma > 0$, this corresponds to fill $k \in [0, \pi]$, thus

$$(n_k)_{ss} = \begin{cases} 1, & k \in [0, \pi], \\ 0, & \text{otherwise.} \end{cases} \qquad (34)$$

We observe this in the numerical results shown in the panel b) of Fig. 3. The expectation value of $\langle \hat{n}_k \rangle$ evolves until all states with a positive imaginary energy are occupied. In contrast with the state in class B the area under the curve, corresponding to the total particle number, is not conserved throughout evolution.

Given the structure of the steady-state occupation, we can immediately infer the steady-state entanglement entropy for the Hatano-Nelson model, which can also be computed analytically by using the Fisher-Hartwig conjecture (see the Appendix C). Indeed, for both classes B and C of initial states, the steady-state of the non-Hermitian dynamics takes the form of a Slater determinant describing fermions in a single imaginary-energy band, which is partially filled. This leads to an entanglement entropy scaling as the logarithm of the subsystem size, i.e.

$$S^{B/C}(+\infty) = \frac{1}{3} \ln \ell + \mathcal{O}(1), \qquad (35)$$

with a prefactor of the log term equal to 1/3 as expected for gapless Dirac fermions described by a (1+1)D conformal field theory with a central charge $c = 1$ [47]. This result is independent
<figure><img src="image_6.png" /><figcaption>Figure 4: Time evolution of the entanglement entropy for different system sizes, starting from initial states in class B (panel a)) and class C (panel b)). The inset in both panel shows the steady-state entanglement entropy as a function of the total system size. The scaling behavior is found to be $S(+\infty) \propto a_2 \ln(L) + \mathcal{O}(1)$, with $a_2 = 0.321 \pm 0.05$ in panel a and $a_2 = 0.3336 \pm 0.0005$ in panel b). This is consistent with the analytical prediction $a_2 = 1/3$ for the thermodynamic limit. Other parameters: $\gamma = 0.4J$.</figcaption></figure>

of the state's filling.<fnref n="1" /> The scaling of the entanglement entropy is confirmed by the numerical results shown in Fig. 4 for both class B and C. While the dynamics presents a rather different behavior, the long-time steady-state entanglement converges to the expected scaling with subsystem size, here for a cut $l = L/4$.

# 5 Application: The non-Hermitian SSH Hamiltonian

In this section, we examine the non-Hermitian Su-Schrieffer-Heeger (SSH) model [15, 17, 18, 27, 48–50], which exhibits richer spectral properties than the Hatano-Nelson model, including PT symmetric and broken phases. We first review the SSH spectral properties and then analyze the steady-state momentum occupation and its connection to the entanglement entropy. In particular, we find that the scaling of the steady-state entanglement entropy depends on the class of the initial state, whether it spontaneously breaks the PT symmetry, and the filling in classes A and B.

## 5.1 The non-Hermitian SSH model

The SSH model describes a spinless fermionic chain with 2 atoms per unit cell, giving rise to two sublattices A, B, and dimerized hopping $-J \pm h/2$, see Figure 5. We consider a non-Hermitian extension of the SSH model where fermions in the A sublattice experience gain, while those in the B sub-lattice are damped. In real-space the non-Hermitian Hamiltonian reads [27, 51]

$$ \mathcal{H}_{\text{nH}} = -\sum_{l=0}^{L-1} \left[ \left( J - \frac{h}{2} \right) c_{B,l}^{\dagger} c_{A,l+1} + \left( J + \frac{h}{2} \right) c_{A,l}^{\dagger} c_{B,l} + \text{h.c.} \right] + i\gamma \sum_{l=0}^{L-1} \left( c_{A,l}^{\dagger} c_{A,l} - c_{B,l}^{\dagger} c_{B,l} \right), \quad (36) $$

which can be written in the form of our original Eq. (1) and contains now two bands corresponding to the sublattice index A, B. In momentum space, it is given by a two band Hamilto-

<footnote n="1">¹The difference with the filling is reflected in $\mathcal{O}(1)$ term.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 5: Scheme of the non-Hermitian SSH model: Each unit cell is enclosed within a grey area. The intra-cell hopping is given by $-J-h/2$, and the inter-cell hopping is given by $-J+h/2$. The red sites (A sublattice) have a local gain term, $i\gamma$, while the blue sites (B sublattice) have local damping, $-i\gamma$.</figcaption></figure>

nian of the form,

$$ \mathcal{H} = \sum_k \Psi_k^\dagger \mathcal{H}_k \Psi_k, \quad \mathcal{H}_k = \begin{pmatrix} i\gamma & J_k \\ J_k^* & -i\gamma \end{pmatrix}, \qquad (37) $$

where $J_k = (-J + h/2)e^{-ika} - (J + h/2)$ and $\Psi^\dagger = (c_{k,A}^\dagger \ c_{k,B}^\dagger)$. The non-Hermitian SSH model is obtained as the no-click limit of a monitored SSH chain, where the local hole density in the A sublattice, $c_{l,A}c_{l,A}^\dagger$, and the local fermionic density in the B sublattice, $c_{l,B}^\dagger c_{l,B}$ are monitored [27].

The Hamiltonian of this non-Hermitian system possesses PT symmetry, as discussed in subsection 3.2. For each momentum $k$, the single-particle eigenstate of the Hamiltonian can either be an eigenstate of PT or spontaneously break this symmetry. As previously emphasized, if PT symmetry is not spontaneously broken, the eigenvalues $\epsilon_{k,\alpha}$ are real. However, if the symmetry is broken, the spectrum becomes complex. For this specific model, one can show that in this case $\epsilon_{k,+} = \epsilon_{k,-}^*$. This can be directly seen from the expression for the single-particle dispersion relation

$$ \epsilon_{k,\pm} = \pm \epsilon_k = \pm \sqrt{h^2 - \gamma^2 + (4J^2 - h^2) \cos^2\left(\frac{k}{2}\right)}. \qquad (38) $$

We define three possible phases based on the total number of real eigenvalues:

* The PT-symmetric phase, where the eigenvalues are real for all momenta $k$;

* The PT-fully-broken phase, where the spectrum is purely imaginary for all momenta $k$;

* The PT-mixed phase, where purely real and purely imaginary modes coexist at different momenta $k$.

In most of the literature [24], both the PT-Fully Broken and PT-Mixed phases are referred indistinguishably as the PT-broken phase. However, this distinction is important for the remainder of our discussion. The three phases are depicted in the panel a) of Fig. 6, where we plot the spectral density of the real eigenvalues defined, in the thermodynamic limit, by the integral $\phi = \int_{-\pi}^{+\pi} dk/2\pi \Theta(\epsilon_k^2)$, where $\Theta(x)$ is the Heaviside step function.

For $\gamma$ smaller than $h$ and $2J$, the entire spectrum is real and gapped, as seen in panel b) of Fig. 6. As $\gamma/J$ increases, the real energy gap decreases and vanishes at $\gamma = h$. At this point, two exceptional points emerge at $k = \pm\pi$, and the system transitions from the PT-symmetric phase to the PT-mixed phase. In this phase, the single-particle energy spectrum is gapless in both its real and imaginary components, as the two bands touch at

$$ k_* = \pm 2 \arccos \sqrt{\frac{\gamma^2 - h^2}{4J - h^2}}, \qquad (39) $$
<figure><img src="image_8.png" /><figcaption>Figure 6: (a) Spectral density of the real eigenstates as a function of the system parameters. The transition lines ($\gamma = 2J$ and $h = \gamma$) are marked in black. The single-particle dispersion relation is shown for the $PT$-symmetric phase in panel (b), the $PT$-mixed phase in panel (c), and the $PT$-fully broken phase in panel (d). The blue line corresponds to the real part of the spectrum, while the red line corresponds to the imaginary part. Other parameters: $h = J$.</figcaption></figure>

for $h < 2J$, as shown in panel c) of Fig. 6. For further convenience, we define $k_* = \pm \pi$ in the $PT$-symmetric phase and $k_* = 0$ in the $PT$-fully-broken phase. With a further increase in $\gamma/J$, another transition occurs at $\gamma_c = 2J$, where the system enters the $PT$-fully broken phase, characterized by a purely imaginary spectrum with a gap, as seen in panel d) of Fig. 6.

## 5.2 Steady-state single-particle distribution

We begin by discussing the steady-state distribution of the single-particle occupation numbers

$$ \langle \hat{n}_k \rangle = \sum_{\alpha=A,B} \langle \hat{n}_{k,\alpha} \rangle, \qquad (40) $$

for different classes of initial states and system parameters. For initial states in class A, $\langle \hat{n}_k \rangle$ remains strictly conserved throughout time evolution, regardless of the value of $\gamma/J$. In contrast, for classes B and C, the distribution is nontrivial. This behavior is illustrated in Fig. 7, where we show the long-time limit of the dynamics for $\langle \hat{n}_k \rangle$ for different values of $\gamma/J$ at fixed $h$, allowing us to examine its dependence across the phase diagram.

In the $PT$-symmetric phase and for both class B and C the steady-state occupation resembles strongly the one at initial times, i.e. a uniform distribution in momentum space, with small fluctuations, indicating that all $k$ modes are equally occupied. This can be understood by recalling that in this phase the many-body spectrum is entirely real, meaning that no eigen-state is, in principle, more damped or enhanced than the others. This phase can be regarded as
<figure><img src="image_9.png" /><figcaption>Figure 7: Steady-state expectation value of $\hat{n}_k$ in the three spectral regimes for different values of the fillings. The values $\gamma$ are 0.5J (panel a) and b)), 1.5J (panel b) and e)), and 2.5J (panel c) and f)). The upper row corresponds to an initial state in Class B, while the bottom row corresponds to an initial state in the Class C. Other parameters: $h = J$.</figcaption></figure>

a minor perturbation of the Hermitian evolution, preserving the system's dynamical stability. Consequently, $\langle \hat{n}_k \rangle$ is almost conserved, becoming strictly conserved in the limit $\gamma \to 0$.

In the *PT*-mixed phase and for class **B**, on the other hand, the physics is highly sensitive to the filling of the initial state. The steady state in this regime is constructed by populating the single-particle modes according to the hierarchy of imaginary energies. When the filling is less than the total number of purely imaginary modes, the steady state becomes an eigenstate of $\hat{n}_k$, as illustrated in Fig. 7 for $\nu = 1/8$. In contrast, when the filling exceeds the number of purely imaginary modes, the steady-state distribution becomes more complex, as presented for the case of $\nu = 1/4$. Here, the quasimomentum values $k$ corresponding to purely imaginary energy modes are fully occupied while those associated with real energy modes exhibit a nontrivial $\langle \hat{n}_k \rangle$ distribution, which is harder to predict. We have verified that this distribution, associated with the real modes, depends on the specific choice of the initial state within class **B**.

For class **C** in the *PT*-mixed phase, the steady-state expectation value of $\hat{n}_k$ is unity for values of $k$ where the spectrum is purely imaginary. In contrast, for $k$ values where the single-particle spectrum remains real, $\langle \hat{n}_k \rangle$ exhibits a nontrivial distribution that depends on the initial state. Figure 7 confirms this behavior, showing that, regardless of the initial expectation value of the filling, the steady state satisfies $\langle \hat{n}_k \rangle = 1$ for all $k$ where the single-particle spectrum is purely imaginary.

Finally, in the *PT*-fully broken phase, we can apply Eq. (18), leading to a steady state that corresponds to the right eigenstate with the highest imaginary energy, subject to the conservation laws imposed by the initial state. For class **B**, where the total number of particles is conserved, the $k$ modes with the largest imaginary energies are filled up to a value $k_F$ which is determined by the total particle number. This is evident in Fig. 7, where the modes around $k = \pm \pi$, corresponding to the highest imaginary eigenenergies, are occupied first as the total particle number increases.
<figure><img src="image_10.png" /></figure>

Figure 8: Steady-state entanglement phase diagram as a function of $\gamma$ and $\nu$ for the different classes of initial states considered. The blue region indicates a volume-law scaling of the steady-state entanglement entropy, the red region corresponds to logarithmic scaling, and the white region represents area-law scaling. For class C, $\nu_0$ denotes the expectation value of the filling in the initial state. The orange lines mark the phase boundaries. In class A we observe different boundaries depending on the chosen central $k_0$ from which the state is filled, the boundary is given by Eq. (45). In class B the boundary is given by Eq. (46) and is the same as the one for $k_0 = \pm\pi$ in class A. The diagram highlights the dominant scaling behavior in the thermodynamic limit. Other parameters: $h=J$.

On the other hand, for class C, where the particle number is not fixed, the least damped right eigenstate consists of all $k$ states being occupied by exactly one particle. This results in a steady state where the entire upper imaginary band is populated, as shown in panel d) of Fig. 7. Thus, independently of the initial conditions, the steady state in this phase is always half-filled. We note that in both cases of class B and C, the steady state has restored the symmetry that was broken by the initial state, respectively translation invariance (for class B) and conservation of particle number (for class C). Our results are therefore consistent and in-line with recent ones exploring symmetry restoration in a related non-Hermitian SSH model [52].

## 5.3 Entanglement entropy

We now discuss the properties of the steady-state entanglement entropy as a function of system parameters and the class of the initial state. Our main results are summarized in the entanglement phase diagram shown in Fig. 8, as a function of dissipation $\gamma$ and filling $\nu$, for the three different classes of initial states. We briefly highlight its main features, which will be discussed in more detail in the following.

For initial states in class A the steady-state entanglement entropy displays either a volume law scaling at weak dissipation or a logarithmic scaling at large dissipation, except for half-filling where the logarithmic phase turns into an area law phase. Due to the conservation of the distribution of $n_k$, the phase diagram in this case strongly depends on the initial state, both in terms of filling $\nu$ and on the offset momentum value $k_0$ used to construct the initial distribution $(n_k)_{t=0}$. As we see in the left panel of Fig. 8, depending on the value of $k_0$ used to construct the initial state, we have different phase boundaries between the volume law and the logarithmic-law phase, including possibly a phase boundary that does not depend on the filling with the exception of the half-filled case where the transition from the volume-law phase is toward an area-law phase.
<figure><img src="image_11.png" /><figcaption>Figure 9: a) Pre-factor of the volume term in the entanglement entropy for an initial state with $k_0 = 0$. This term vanishes at $\gamma_c = 2J$. b) Pre-factor of the logarithmic term for the initial state with $k_0 = 0$. The orange dashed lines indicate the change of $a_2$ into $1/3$. At half-filling the component $a_2$ is always zero. c) Pre-factor of the volume term in the entanglement entropy for a fixed filling $v = 1/8$ and initial states with different $k_0$. d) Pre-factor of the logarithmic term in the entanglement entropy for a fixed filling $v = 1/8$ and initial states with different $k_0$. For $\gamma$ large enough we always find $a_2 = 1/3$. Other parameters: $h = J$.</figcaption></figure>

For initial states in class B, the phase diagram also displays three phases characterized, respectively, by a volume law scaling (at weak dissipation), a logarithmic scaling of the entanglement entropy at strong dissipation and a line with area-law scaling for an initial filling equal to $v = 1/2$. In contrast to the previous case, the steady-state scaling of the entanglement entropy does not depend on the specific form of the initial states within this class. The boundary of the entanglement transition between the volume law phase and logarithmic law phase exhibits a nontrivial dependence on the filling, making it possible to drive the transition by tuning either the filling or the dissipation.

Finally, in class C, the phase diagram features only two phases: a volume-law phase under weak dissipation, corresponding to the *PT*-symmetric and *PT*-mixed phases, and an area-law phase under strong dissipation, corresponding to the *PT*-fully broken phase. These phases emerge independently of the expectation value of the filling in the initial state. We now provide analytical and numerical evidence in support of the general picture drawn so far on the entanglement structure of the non-Hermitian SSH.
### 5.3.1 Class A

For initial states in this class, when the occupation number $\langle \hat{n}_k \rangle$ is conserved and the correlation matrix $G_k^{\alpha\beta}(+\infty)$ diagonal in momentum, it is possible to solve the dynamics and determine the steady-state entanglement entropy analytically, throughout the phase diagram. We note that for this class the specific form of the initial state is crucial, as it determines and fixes the expectation values of $\hat{n}_k$ for all quasimomenta $k$. Here we consider initial states defined as Slater determinants of single-particle momentum eigenstates filled from $k_0$. Specifically, we take

$$|\Psi_0\rangle = \prod_{|k-k_0|<2\pi\nu} D_k^\dagger |\text{vac}\rangle, \quad (41)$$

with $|k-k_0|$ defined modulo $2\pi$, and $D_k^\dagger$ is a compact notation to consider the possibility that a given filled momentum $k$ can be single or doubly occupied

$$D_k^\dagger = \begin{cases} \frac{c_{kA}^\dagger + c_{kB}^\dagger e^{ik/2}}{\sqrt{2}} & \text{(single occupancy),} \\ \left( \frac{c_{kA}^\dagger - c_{kB}^\dagger e^{-ik/2}}{\sqrt{2}} \right) \left( \frac{c_{kA}^\dagger + c_{kB}^\dagger e^{ik/2}}{\sqrt{2}} \right) & \text{(double occupancy).} \end{cases} \quad (42)$$

For $k_0 = 0$, this corresponds to the ground-state of the hermitian SSH model ($\gamma = 0$) with $h = 0$. Despite this specific choice of initial state, all subsequent considerations and discussions remain valid and the results can be mapped to other initial distributions.

After the analytical manipulations, that we report in Appendix C for the interested reader, we obtain two distinct contributions for the steady-state entanglement entropy, besides the irrelevant area law coefficient,

$$S_{\ell}(+\infty) = a_1 \ell + a_2 \ln \ell + \mathcal{O}(1). \quad (43)$$

As shown in Appendix C, the prefactors $a_1$ and $a_2$, controlling the volume-law and logarithmic-law term respectively, can be computed using the Szegö lemma and the Fisher-Hartwig conjecture [42] through the eigenvalues $\mu_{k,\pm}$ of the steady-state correlation matrix $G_k^{\alpha\beta}(+\infty)$. We obtain the following expression,

$$a_1 = \int_V \frac{dk}{2\pi} [s(\mu_{k,+}) + s(\mu_{k,-})], \quad (44)$$

where the integration is performed over the domain $V = [-k_*, k_*] \cap [k_0 - 2\pi\nu, k_0 + 2\pi\nu]$ with $k_*$ defined in Eq. (39), and the binary entropy function is defined as

$$s(x) = -(1-x)\ln(1-x) - x \ln(x).$$

As obtained from the exact calculation, only the eigenvalues of the correlation matrix, $G_k^{\alpha\beta}(+\infty)$, corresponding to a single-occupied quasimomentum $k$ with real energy, contribute to the volume coefficient $a_1$.

Panel a) of Fig. 9 shows the numerical value of $a_1$ for an initial state with $k_0 = 0$. In this case, for all fillings, the coefficient $a_1$ vanishes only when the system transitions to the PT-Fully Broken phase, i.e. for $\gamma = 2J$, which therefore signals the transition out of the volume-law phase for $k_0 = 0$ (see the phase diagram in Fig. 8). The filling does not affect the transition for this particular initial state ($k_0 = 0$), but in general it does it for other initial states as we show in Fig. 8(a).

In the panel c) of Fig. 9, we show, for a fixed value of the filling $\nu = 1/8$, the dependence of the coefficient $a_1$ with $\gamma$ for initial states with distinct $k_0$. For this filling and $k_0 \neq 0$ we
see that $a_1$ vanishes at a dissipation value $\gamma < 2J$, i.e. the transition out of the volume law phase does not coincide any longer with the full breaking of the PT-symmetry, as was the case for the half-filled SSH model discussed previously [27]. We can derive an expression for the critical value of $\gamma$ at which $a_1$ vanishes, corresponding to the boundary between the phase with volume-law scaling of the entanglement and the phase exhibiting logarithmic or area-law scaling. The critical value of dissipation reads

$$ \gamma_c(k_0, \nu) = \left( h^2 + (4J^2 - h^2) \cos \left[ \frac{\max(|k_0| + \pi|\nu - 1| - \pi, 0)}{2} \right] \right)^{1/2} . \quad (45) $$

These boundaries are plotted in the panel a) of Fig. 8 for different values of $k_0$, ranging from 0 to $\pi$ (the expression is symmetric in $k_0$). The transition is only independent of the filling when $k_0 = 0$, while for other values, the phase boundary shows a strong dependence on $\nu$, especially away from half-filling.

Returning to the general result in Eq. (43) we see that once $a_1$ vanishes the entanglement entropy scales as the logarithm of the subsystem size, whenever the coefficient $a_2 \neq 0$. Interestingly, this sub-leading term only appears when the system is doped away from half-filling. In fact $a_2$ originates from discontinuities in the two-point correlation matrix $G_k^{\alpha\beta}$ near $k_0 \pm 2\pi\nu$. The integral yielding $a_2$ is detailed in Appendix C.

In panel b) of Fig. 9, we show the numerical values of $a_2$ for different fillings and $\gamma/J$ for the initial state with $k_0 = 0$. We note that $a_2$ is different from zero throughout the phase diagram, except along the line $\nu = 1/2$, however, its effect is visible only in regions where the leading term $a_1 = 0$. We can understand the origin of the logarithmic phase and its transition into an area law at half-filling from the perspective of our general framework. Indeed, whenever the spectrum is purely imaginary (as for $\gamma \ge 2J$) the steady state is represented by a Slater determinant obtained by filling the single-particle states with the highest imaginary energies. When the state is away from half filling, we thus get the purely imaginary spectrum to be partially filled, which explains the logarithmic scaling of the entanglement entropy. We also note that the logarithmic contribution to the entanglement entropy features a prefactor of $1/3$, as expected from a free-fermionic Conformal Field Theory. In contrast, when the system is exactly at half-filling, one of the imaginary bands is completely filled while the other remains empty, and the spectrum in that case is effectively gapped. This scenario is analogous to the ground-state of an insulating band model, resulting in an area-law behavior for the entanglement entropy.

Finally, in panel d) of Fig. 9, we show cuts of the coefficient $a_2$ at fixed filling $\nu = 1/8$ and for different initial states parametrized by $k_0$. We can then observe that the volume law to area law transition occurs at a different $\gamma_*(\nu = 1/4)$ depending on $k_0$. This can be understood from Eq. (44) since $\nu$ is based on an intersection that changes depending on $k_0$ and can more easily be null for small and large filling.

### 5.3.2 Class B

For the initial states in class B, we investigate the steady-state scaling of the entanglement entropy by propagating the initial state using the Faber polynomial method [37] which is detailed in appendix A. In the following discussion, we fix $h = J$, and tune $\gamma/J$ across the phase diagram in Fig. 10 and the filling of the initial state.

In the PT-symmetric phase, that is, for $\gamma < \gamma_c(h = 1) = J$, the entanglement entropy scales linearly with the size of the system, $S_\ell \sim \ell$, as shown in panel a) of Fig. 10 for a cut $\ell = L/4$. This result is independent of the system's filling, which controls only the prefactor of the linear term. We can understand this result by noticing that PT-symmetric phase is
<figure><img src="image_12.png" /><figcaption>Figure 10: Steady-state entanglement entropy for an initial state in the class B, across the three different spectral regimes for various fillings $v$. The values of $\gamma$ are 0.5$J$ (a)), 1.5$J$ (b)) and 2.5$J$ (c)). The parameter $L$ corresponds to the total number of unit cells. Other parameters: $h = J$.</figcaption></figure>

similar to and adiabatically connected<fnref n="2" /> to the unitary evolution, where a volume-law scaling of the entanglement entropy is generically observed after a global quantum quench in free fermionic systems [3]. Furthermore, as discussed in Sec. 3, in this phase the system's steady state can be described by a diagonal ensemble, which usually encodes effectively thermal or (generalized) equilibrium states. This further corroborates the volume-law scaling of the entanglement entropy.

Upon increasing the dissipation and entering the *PT*-mixed phase, volume-law scaling of the entanglement entropy becomes strongly dependent on the initial filling. In particular, we see in panel *b*) of Fig. 10 that by varying the filling of the initial state, one can tune an entanglement transition from volume scaling to logarithmic scaling. This result can be understood within our general framework: volume-law scaling is observed whenever the state filling exceeds the number of quasimomenta $k$ with purely imaginary energies (see panel *b*) of Fig. 10) and the modes with real energies are singly occupied. Conversely, when the filling is less than the number of single-particle modes with purely imaginary energies, the entanglement entropy exhibits logarithmic scaling. A similar behavior occurs if the filling is such that all modes associated with real energies are doubly occupied, as illustrated in the phase diagram of Fig. 8 and seen explicitly for the case $v = 7/8$ in Fig. 10.

The results above and in particular the filling-driven entanglement transition between logarithmic and area-law phase are fully consistent, at least at the qualitative level, with our general framework. Technically, however, the existence of a pair of exceptional points prevents us from fully diagonalizing the Hamiltonian throughout the Brillouin zone. Indeed, as discussed in Sec. 5.2 in this *PT*-mixed phase the steady state is a Slater determinant constructed by occupying the single-particle states with the highest imaginary energies. Consequently, it is possible to apply the Szegö lemma and the Fisher-Hartwig conjecture to corroborate the numerical results shown in Fig. 7.

In the *PT*-fully broken phase, the spectrum is entirely imaginary, and so we can apply our general framework to make analytical progress. In this scenario, the steady state is described by a Slater determinant formed from the slowest decaying single-particle states. Using the Fisher-Hartwig conjecture, the steady-state entanglement entropy exhibits logarithmic scaling with a $1/3$ prefactor in the gapless phase, while in the specific case of half-filling, it adheres to an area law (see Appendix C for details). This behavior is supported by finite-sized system numerical simulations, as illustrated in Fig. 10. In the thermodynamic limit, the steady state

<footnote n="2">²When tuning $\gamma \to 0$ no gap is open in the single-particle spectrum.</footnote>
can once again be understood as a Slater determinant, interpreted as the low-energy state of a free-fermionic conformal field theory. The prefactor of the logarithmic scaling corresponds to a central charge $c = 1$, whenever $\nu \neq 1/2$.

We can obtain an analytical expression for the phase boundary in Fig. 8 separating the volume-law phase from the logarithmic-law phase, as a function of the initial filling. This boundary can be determined by examining whether, for a given filling, any quasimomenta $k$ associated with real energy are singly populated. If this is the case, one obtains a volume law scaling for the entanglement entropy. Conversely, if only purely imaginary modes are occupied and the real modes are empty or doubly occupied, a logarithmic scaling occurs, except for $\nu = 1/2$, which follows an area law scaling. The critical value of $\gamma$ is then given by

$$ \gamma_c = \left[ (4J^2 - h^2) \cos^2\left(\frac{\pi}{2}(1-2\nu)\right) + h^2 \right]^{1/2}. \qquad (46) $$

This is the same boundary found for states considered in class A with $k_0 = \pi$. This is the case in class A since filling from $k_0 = \pi$ corresponds to filling according to the highest imaginary single-particle energies.

### 5.3.3 Class C

In this last section, we discuss the steady-state entanglement structure for the initial state in class C. In Fig. 11, we plot the scaling of the entanglement entropy for a partition $l = L/4$ for three values of $\gamma$, respectively in the PT-symmetric, mixed and fully-broken phase and for different values of the expectation value for the initial filling $\nu_0$. In the PT-symmetric and mixed phase, we consistently observe a volume-law scaling of the steady-state entanglement entropy, see Fig. 11(a, b), with the initial filling controlling the slope of the entanglement entropy, i.e. the coefficient of the volume-law term. In contrast, in the fully broken PT phase, the steady-state entanglement is independent of system size, i.e. the system enters the area law phase. In other words, for initial states in class C the phase diagram displays a volume-to-area law entanglement transition, which we locate numerically at $\gamma_c = 2J$, independently on the initial filling $\nu_0$. This value coincides with the spectral transition into the fully broken PT-phase.

<figure><img src="image_13.png" /><figcaption>Figure 11: Steady-state entanglement entropy for different system parameters for the states in class C. The values of $\gamma$ are 0.5J (a)), 1.5J (b)) and 2.5J (c)). Other parameters: $h = J$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: Fluctuations of the expectation value of the total particle number for states in class C and a system with linear size $L = 512$. a) Dynamics of the intensive fluctuations for different values of $\gamma$. b) Steady-state value of the intensive particle number fluctuations as a function of $\gamma/J$ and different linear system sizes. We see that charge-sharpening is not effective for $\gamma < J$ and it is fully complete for $\gamma > 2J$, leaving in between a phase where charge is partially sharpened. In the top of panel b), we identify the corresponding spectral phase for this specific ratios of $\gamma/J$ with $h=J$.</figcaption></figure>

The above results are again fully consistent with our general framework for the steady-state occupancy. Indeed, as discussed in Sec. 5.2, for initial states in class C the steady state always fills up states with zero imaginary eigenvalue (i.e. purely real eigenvalues) whenever available and therefore the leading contribution to the entanglement entropy displays volume law scaling in the $PT$ and $PT$-mixed phases. In the fully broken $PT$ phase, when the spectrum is complex, the system relaxes into the right eigenstate of the many-body Hamiltonian. As discussed earlier, this corresponds to the Slater determinant obtained by filling all right single-particle states with the highest imaginary energy. Consequently, the correlation matrix is diagonal in the $k$ basis, leading to an area-law scaling of the entanglement entropy. The steady state resembles that of an insulator as it consists of a fully occupied upper imaginary band.

Finally, it is interesting to note that initial states in class C which are linear superposition of states with different particle numbers allow to explore charge-sharpening dynamics in a non-Hermitian setting. For monitored random circuits with U(1) symmetry, it is known that a charge-sharpening transition occurs within the volume-law phase [53]. This transition distinguishes a phase where initial charge superpositions are not effectively collapsed by the monitoring protocol, remaining fuzzy along individual trajectories for times $t \sim L$, from a regime where this process occurs rapidly: the system is in a well-defined charge eigenstate well before becoming fully pure. It is therefore tempting to investigate whether a similar behavior occurs within the volume-law phase of our non-Hermitian SSH system, where the charge associated with the U(1) symmetry here is the total particle number. In Fig. 12 a), we plot the dynamics of the fluctuations of the total particle number $\hat{N} = \sum_{i\alpha} \hat{n}_{i\alpha}$, for different values of $\gamma$. We see that for weak dissipation the fluctuations are almost constant in time and equal to the ones encoded in the initial state. Upon further increasing the dissipation rate, we observe a nontrivial dynamics that reaches a finite steady-state value, i.e. the non-Hermitian evolution has sharpened the charge but not completely. Finally, for large values of $\gamma$ the dynamics collapse the initial superposition and fluctuations at long times vanish. In Fig. 12 b), we plot the steady-state value of the intensive fluctuations, which clearly distinguishes three regimes: for $\gamma < J$ the fluctuations are almost unaffected by the dissipation; for $J < \gamma < 2J$ the residual fluctuations are small and decrease with $\gamma$, and for $\gamma > 2J$ the system at long-times
is in eigenstate of the total particle number, that is, the symmetry is restored and the particle-
number fluctuations vanish. Quite interestingly, the three regimes above coincide with the
spectral transition where PT-symmetry is first partially broken ($\gamma = J$) and then fully broken
($\gamma = 2J$). This can be understood within our framework: once PT-symmetry breaks, the sys-
tem develops purely imaginary eigenmodes. These will be populated at long times, reducing
the fluctuations of the occupation numbers $n_k$, and, consequently, of the total particle number.
This process continues upon increasing $\gamma$ and completes once PT is fully broken, at which
point the charge fluctuations completely vanish.

Although our simple model of non-interacting fermions is not able to reproduce the full
features of the charge sharpening criticality observed in random circuits, in particular the
different behavior of the time-scales controlling sharpening and purification, it is already quite
remarkable that such phenomenology is observed in the no-click limit of purely non-Hermitian
evolution. This is even more so considering that our model is non-interacting: for comparison,
non-interacting monitored fermions with U(1) symmetry do not sustain a volume-law phase,
and so, they also lack a charge sharpening phase. Future work is required to investigate this
non-Hermitian charge sharpening transition in a fully interacting and chaotic case.

6 Conclusions

In this work, we have introduced a general theoretical framework to understand the steady-
state structure of non-interacting fermionic non-Hermitian lattice models and their entangle-
ment content. We have highlighted in particular the role of symmetries and conserved quan-
tities in constraining the dynamics. A unique feature of non-Hermiticity which directly arises
from the nonlinearity of measurement back-action is that conserved quantities arise from an
interplay between Hamiltonian symmetries and initial states breaking or preserving such sym-
metry. Focusing on translation-invariant free fermionic models with U(1) symmetry, we have
identified three classes of initial states (class A, B, and C) which break at different levels the
symmetries of the problem and therefore have a strong impact on the steady-state properties.

We have shown that for generic non-interacting non-Hermitian Hamiltonian with complex
spectrum, the steady state can be obtained by filling single-particle states with the highest
imaginary eigenvalue, compatibly with the constraints introduced by the initial states. Initial
states in class A lead to a dynamic that conserves the average occupation ⟨$\hat{n}$<sub>k</sub>⟩, leaving strong
imprints on the structure of the steady state. Initial states in class B are characterized by a con-
served total number of particles, and as a result, the initial filling plays a key role in shaping the
steady state. In particular, we have shown that for initial states in this class, the steady state
can feature either partially or completely filled bands of imaginary eigenvalues. Finally, the
initial states in class C are characterized by breaking the particle number conservation, lead-
ing to a filling of the imaginary eigenvalue bands without additional constraints. Interestingly,
while the initial states break the symmetries of the non-Hermitian Hamiltonian, the dynamics
restores them in the steady state. Since for fermionic gaussian states the knowledge of the
correlation matrix determines the structure of the entanglement entropy, our results allow to
conclude that for non-Hermitian fermions with complex energy spectrum, the entanglement
scaling in the steady-state is either area law, if the imaginary eigenvalue bands are gapped and
half-filled, or scaling with the logarithm of system size for partially filled bands. This imme-
diately implies that the structure of the steady-state entanglement can be strongly dependent
on filling, a feature that is remarkably different from the unitary case and resembles the case
of ground-state problems.
We have also discussed the situation in which, due to *PT*-symmetry, the energy spectrum is either fully real (*PT*-symmetric phase) or partially real and partially complex (so-called *PT*-mixed phase). In the first case, we have argued that the steady-state density matrix can be described in terms of a diagonal ensemble, and due to the purely real-spectrum, the entanglement entropy shows volume law scaling. In the *PT*-mixed phase, the steady-state occupation depends nontrivially on filling: if the purely real energy is partially occupied, the entanglement entropy still shows a volume law scaling.

We have applied our framework to two relevant models of non-Hermitian fermionic lattices, the Hatano-Nelson and the Su-Schrieffer-Heeger models. In the former case, we have shown how the dynamics of single-particle states always leads to partially filled bands, giving rise to a steady-state entanglement that is logarithmic in the system size as in critical systems. For the SSH model, the landscape is richer because of its spectral properties, involving *PT* symmetry and its breaking. This has direct consequences on the structure of the steady-state entanglement. In particular, we have shown that in the *PT*-symmetric phase the dynamics for all initial states resembles the unitary evolution. In particular, local observables thermalize to a diagonal ensemble, and the entanglement entropy displays a volume-law scaling, independently of the initial state.

On the other hand, when the dissipation is large enough such that *PT*-symmetry is fully broken and the spectrum is purely imaginary, the initial state strongly affects the dynamics. We have shown in Sec. 3 that the steady state is obtained by filling single-particle orbitals with the highest imaginary eigenvalue, compatibly with the initial filling for states in class **B**, or independently of the filling for class **C**. The entanglement entropy is therefore either area law or logarithmic. In contrast, in the *PT*-mixed phase, the situation is richer for classes **A** and **B**, with a steady-state occupation that depends nontrivially on filling. As such, the entanglement entropy shows a transition driven by filling between logarithmic scaling and volume law scaling. However, for initial states in class **C** in the *PT*-mixed phase, the entanglement entropy still shows the volume law. Nevertheless, we have unveiled a nontrivial dynamics of the fluctuations of the total particle number which displays a phenomenology which resembles the charge-sharpening of random circuits. To summarize, our results for the non-Hermitian SSH show that by combining conservation laws, *PT*-symmetry, and the initial state's filling one can generate rich entanglement patterns.

Our work opens up several directions for further investigation. Relaxing the assumption of translation invariance and periodic boundary conditions could allow to discuss entanglement patterns in disordered free fermionic systems or in systems with open boundary conditions, where the skin-effect is expected to play an important role. It would be interesting to explore whether our framework could be generalized to these cases as well. Similarly, bosonic non-Hermitian systems have been shown to display rich physics and could represent another arena to explore the connection between conserved quantities, steady-state occupations, and entanglement. The major challenge is represented by interacting non-Hermitian many-body systems, which would be exciting to address using extensions of our theoretical framework. Finally, it would be interesting to explore whether these features could persist in the context of monitored systems. It has been numerically shown that in certain models, the monitored dynamics of the entanglement entropy closely agrees with the no-click limit [51, 54]. This also depends on how the symmetries of the quantum jump operators interplay with those of the Hamiltonian [54], so it would be interesting to explore whether, by appropriately choosing the quantum jump operators, there could exist a dependence of the measurement-induced entanglement transition on the filling of the state. Similarly, it would be interesting to explore whether some of our results can be generalized to the dynamics of entanglement of the averaged mixed state, using different entanglement estimators such as the negativity. In this respect it is known that for non-interacting driven-dissipative systems the dynamics of the cor-
relation matrix is also controlled by a non-Hermitian Hamiltonian, yet generally different from the no-click one [13]. One could then use similar principles to the ones used in this work to make progress in this problem, a direction that we leave for future work.

## Acknowledgments

**Funding information** R.D.S acknowledges funding from Erasmus+ (Erasmus Mundus program of the European Union). We acknowledge Collège de France IPH cluster where the numerical calculations were carried out.

## A Numerical simulation of non-Hermitian lattices

In this appendix, we describe the Faber Polynomial method [37] that we use throughout this work to numerically perform the time-evolution of the states in classes **B** and **C**. In practice, the time-evolution operator in Eq. (5) is expanded in a Faber Polynomial series,

$$e^{-i\mathcal{H}_{\text{nH}}t} |\Psi(t_0)\rangle = \sum_{n=0}^{+\infty} c_n(t) F_n\left(\frac{\mathcal{H}_{\text{nH}}}{\lambda}\right) |\Psi(t_0)\rangle, \qquad (A.1)$$
$$c_n(t) = e^{-i\lambda\gamma_0 t} \left(\frac{-i}{\sqrt{\gamma_1}}\right)^n J_n(2\sqrt{\gamma_1}\lambda t),$$

where $F_n(\cdot)$ is the n^th Faber Polynomial, $J_n(\cdot)$ is the n^th Bessel function of the first kind, $\lambda$ is a rescaling parameter and $\gamma_0$ and $\gamma_1$ are parameters related with the bounds of the spectrum of the non-Hermitian Hamiltonian, consult [37] for the specific details. In practice, the series is truncated up to a finite order, and the action of $F_n(\mathcal{H}_{\text{nH}}/\lambda) |\Psi(t_0)\rangle$ is done efficiently by using the recurrence relation satisfied by the Faber Polynomials, namely,

$$\begin{aligned}
|\Psi_0\rangle &= |\Psi(t_0)\rangle, \\
|\Psi_1\rangle &= (\tilde{\mathcal{H}}_{\text{nH}} - \gamma_0) |\Psi_0\rangle, \\
|\Psi_2\rangle &= (\tilde{\mathcal{H}}_{\text{nH}} - \gamma_0) |\Psi_1\rangle - 2\gamma_1 |\Psi_0\rangle, \\
|\Psi_{n+1}\rangle &= (\tilde{\mathcal{H}}_{\text{nH}} - \gamma_0) |\Psi_n\rangle - \gamma_1 |\Psi_{n-1}\rangle, \quad n > 2,
\end{aligned} \qquad (A.2)$$

where $|\Psi_n\rangle = F_n(\mathcal{H}_{\text{nH}}/\lambda) |\Psi(t_0)\rangle$. This approach is highly memory efficient, as it requires storing only the last two vectors in memory, $|\Psi_{n-1}\rangle$ and $|\Psi_n\rangle$, to compute the next one, $|\Psi_{n+1}\rangle$. After each time step, the state is normalized. In practice, the expansion in Eq. (A.1) is done at the level of the single-particle Hamiltonian.

For the states in class **B**, as these have a well-defined particle number M, the many-body state can always be represented in the form

$$|\Psi(t)\rangle = \prod_{n=0}^{M-1} \left[ \sum_{l=0}^{L-1} U_{ln}(t) c_l^\dagger \right] |\text{vac}\rangle, \qquad (A.3)$$

with M the total number of particles and L the dimension of the single-particle Hilbert space. The time evolution is then reduced to the following equation

$$U_n = e^{-ith} U_n, \qquad (A.4)$$

where $U_n$ is the n^th column vector and h is the single-particle Hamiltonian matrix, a L x L matrix. In this case, we proceed by expanding $\exp(-ith)$ in a Faber series as described above.
The Faber polynomial method can also be used to deal with the states in class C [54]. In this case, the fermionic state is parameterized according to

$$|\Psi(t)\rangle = \mathcal{N}(t) \exp\left(-\frac{1}{2} \sum_{m,n} \left[(U_t^\dagger)^{-1} V_t^\dagger\right]_{m,n} c_m^\dagger c_n^\dagger\right) |\text{vac}\rangle, \quad (A.5)$$

where $\mathcal{N}$ enforces the correct normalization and $U$ and $V$ are L × L matrices, where L is the dimension of the single-particle Hilbert space. These matrices evolve according to

$$\begin{pmatrix} U_t \\ V_t \end{pmatrix} = e^{-2i\mathcal{H}_{\text{nH}}} \begin{pmatrix} U(0) \\ V(0) \end{pmatrix}, \quad (A.6)$$

where $\mathcal{H}_{\text{nH}}$ is the single-particle non-Hermitian Hamiltonian in the Nambu representation. So, for the states in class C, we develop $\exp(-it\mathcal{H}_{\text{nH}})$ in a Faber series as described above.

## B Diagonal ensemble

In this Appendix, we derive the diagonal ensemble of Eq. (25). We demonstrate that, analogously to Hermitian systems, one can construct a diagonal ensemble that accurately describes the long-time behavior of the expectation value of a local observable $\mathcal{O}$ in non-Hermitian models. Importantly, this construction is valid only for non-Hermitian systems with a real spectrum. So, it should only be applicable to systems driven by $\mathcal{PT}$-symmetric Hamiltonians in the $\mathcal{PT}$-symmetric phase, or to generic pseudo-Hermitian Hamiltonians.

To construct the diagonal ensemble, we first perform an eigendecomposition of the Hamiltonian in a bi-orthogonal basis

$$\mathcal{H}_{\text{nH}} = \sum_n E_n |\Psi_n^R\rangle\langle\Psi_n^L|, \quad (B.1)$$

where $|\Psi_n^{R/L}\rangle$ is the many-body right/left $n^{th}$ eigenstate and $E_n \in \mathbb{R}$ the associated eigenenergy. We processed by expressing the long-time expectation value of $\mathcal{O}$ as

$$\overline{\langle\mathcal{O}\rangle} = \lim_{\tau \to +\infty} \frac{1}{\tau} \int_0^{\tau} dt \, \langle\Psi(t)|\mathcal{O}|\Psi(t)\rangle. \quad (B.2)$$

By making use of the eigendecomposition of the Hamiltonian, this is equivalent to

$$\overline{\langle\mathcal{O}\rangle} = \sum_{n,m} \lim_{\tau \to +\infty} \frac{1}{\tau} \int_0^{\tau} dt \frac{\langle\Psi(0)|\Psi_n^L\rangle \langle\Psi_m^L|\Psi(0)\rangle \langle\Psi_n^R|\mathcal{O}|\Psi_m^R\rangle}{\mathcal{N}^2(t)} e^{i(E_n - E_m)t}, \quad (B.3)$$

with the normalization factor in the denominator equal to

$$\mathcal{N}^2(t) = \sum_{n,m} \langle\Psi(0)|\Psi_n^L\rangle \langle\Psi_m^L|\Psi(0)\rangle \langle\Psi_n^R|\Psi_m^R\rangle e^{i(E_n - E_m)t}. \quad (B.4)$$

The overlap $\langle\Psi_n^R|\Psi_m^R\rangle$ is generically non-vanishing for $n \neq m$. In the case of the Hamiltonian of the type of Eq. (1), the overlap between many-body eigenstates with different quantum numbers occurs only for a very small subset of states. This is particularly evident in the one-body sector, where each right eigenstate $|\Psi_n^R\rangle$ has a non-zero overlap with precisely $n-1$ other right eigenstates, with $n$ being the total number of bands.
<figure><img src="image_15.png" /><figcaption>Figure 13: Comparison of the expectation value predicted by the diagonal ensemble with the exact numerical evolution for a system of size L = 512 unit cells. The effects of finite system size are evident in the small oscillations around the steady-state value. Other parameters: $\gamma = 0.5J$ and $h = J$.</figcaption></figure>

Assuming that we can perform the time-average of the norm and the state independently, we obtain the following expressing for the long-time limit of $\mathcal{O}$ as

$$ \overline{\langle \mathcal{O} \rangle} = \frac{1}{\mathcal{N}^2} \left[ \sum_{\alpha} |\langle \Psi(0) | \Psi_{\alpha}^{L} \rangle|^2 \langle \Psi_{\alpha}^{R} | \mathcal{O} | \Psi_{\alpha}^{R} \rangle + \lim_{\tau \to +\infty} \sum_{\alpha \neq \beta} \langle \Psi(0) | \Psi_{\alpha}^{L} \rangle \langle \Psi_{\beta}^{L} | \Psi(0) \rangle \langle \Psi_{\alpha}^{R} | \mathcal{O} | \Psi_{\beta}^{R} \rangle \frac{e^{i(E_{\alpha} - E_{\beta})\tau} - 1}{i\tau(E_{\alpha} - E_{\beta})} \right]. \quad (B.5) $$

Assuming that phase coherence between the different exponential contributions is lost, the second term in the equation should decay to zero in the long-time limit. Consequently, the expectation value reduces to a weighted average of $\mathcal{O}$ evaluated in the right eigenstates of the Hamiltonian,

$$ \overline{\langle \mathcal{O} \rangle} = \frac{1}{\mathcal{N}^2} \sum_{\alpha} |\langle \Psi(0) | \Psi_{\alpha}^{L} \rangle|^2 \langle \Psi_{\alpha}^{R} | \mathcal{O} | \Psi_{\alpha}^{R} \rangle, \quad (B.6) $$

where $\mathcal{N}^2 = \sum_{\alpha} |\langle \Psi(0) | \Psi_{\alpha}^{L} \rangle|^2$. So, the density matrix predicted from this construction is given by

$$ \rho_{DE} = \frac{1}{\mathcal{N}^2} \sum_{\alpha} |\langle \Psi(0) | \Psi_{\alpha}^{L} \rangle|^2 |\Psi_{\alpha}^{R} \rangle \langle \Psi_{\alpha}^{R} |. \quad (B.7) $$

This result is similar to the one derived for an Hermitian model [5]. The main difference lies in the existence of the right and left eigenstates.
## B.1 Testing the diagonal ensemble

In this subsection, we make a comparison between the results predicted by the diagonal ensemble $\rho_{DE}$, with those given by the exact numerical simulation of the lattice mode. We take the SSH Hamiltonian of Eq. (37) in the $PT$-symmetric regime. In this comparison, we only consider observables that are local in real space degrees of freedom such as the onsite particle density or the two-point correlation function between neighboring sites. We expect that only the expectation value of local observables to relax to the prediction of the diagonal ensemble, as the system as a whole remains pure at all times. Only the reduced density matrix describing a given small region $A$, $\rho_A = \text{Tr}_{\bar{A}}(|\Psi(t)\rangle\langle\Psi(t)|)$, is expected to reach a steady state that can be described by the diagonal ensemble

$$ \lim_{t \to +\infty} \rho_A(+\infty) = \rho_{DE}. \qquad (B.8) $$

In the following, we take the initial state to be

$$ |\Psi_0\rangle = \sum_{l=0}^{N-1} \psi_l c_l^\dagger |\text{vac}\rangle, \qquad (B.9) $$

with $\sum_l |\psi_l|^2 = 1$. In this way, the initial state has a significant overlap with a large number of left eigenstates of the Hamiltonian.

In Fig. 13, we compare the expectation values of the two-point evolving fermionic operators within each unit cell and between adjacent cells considering the initial state $|\Psi_0\rangle = \prod_{n=0}^{L} c_{n,A}^\dagger |0\rangle$. Our results indicate that the diagonal ensemble accurately captures the steady-state value of these local operators.

## C Analytical calculation of the entanglement entropy

In this appendix, we explain how to calculate the entanglement entropy using the Szegö lemma [3] and the Fisher-Hartwig conjecture [42–44].

For fermionic Gaussian states, the entanglement entropy can be computed using the two-point correlation matrix truncated to the sites within the region of interest. The entanglement entropy $S_\ell$ of a region $\ell$ is given by the expression

$$ S_{\ell} = \frac{1}{2\pi i} \lim_{\epsilon \to 0^{+}} \oint_{\Gamma} d\lambda s(\epsilon, \lambda) \frac{d}{d\lambda} \ln \det \mathcal{D}_{\lambda}^{\ell}, \qquad (C.1) $$

with

$$ \mathcal{D}_{\lambda}^{\ell} = \lambda \mathcal{I}_{\ell} - \mathcal{G}_{\ell}, \quad \text{and} \quad s(\epsilon, x) = -(x + \epsilon)\ln(x + \epsilon) - (1 + \epsilon - x)\ln(1 + \epsilon - x), $$

where $\lambda \in \mathbb{C}$, $\mathcal{I}_\ell$ is the identity matrix corresponding to the region $\ell$, and $\mathcal{G}_{i,j}^{\alpha\beta}|_{i,j \in \ell} = \langle c_{i\alpha}^{\dagger} c_{j\beta} \rangle$ is the real-space correlation matrix truncated to $\ell$. The integral is performed along a contour $\Gamma$, which encircles the interval $[0, 1]$ and avoids the branch cuts $(-\infty, -\epsilon) \cup (1+\epsilon, +\infty)$ [42].

In general, calculating the determinant of $D_\lambda^\ell$ analytically is challenging. However, analytical progress can be done if the correlation matrix is a Toeplitz matrix, where $\mathcal{G}_{n,m}^{\alpha\beta} = \mathcal{G}_{n-m}^{\alpha\beta}$, and thus it is diagonal in the momentum sector:

$$ \mathcal{G}_{n,m}^{\alpha\beta} = \int_{-\pi}^{+\pi} \frac{dk}{2\pi} e^{ik(n-m)} \mathcal{G}_k^{\alpha\beta}. \qquad (C.2) $$
In this case, for $\ell \gg 1$, the determinant of $\mathcal{D}_{\lambda}^{\ell}$ can be calculated using the Szegö lemma and the Fisher-Hartwig conjecture:

$$ \ln \det \mathcal{D}_{\lambda}^{\ell} = \ell \int_{-\pi}^{\pi} \ln \det \mathcal{D}_k \frac{dk}{2\pi} + \frac{\ln \ell}{4\pi^2} \sum_{r=1}^{R} \mathrm{Tr} \left[ \ln \left( \mathcal{D}_{k,r}^{-} \left( \mathcal{D}_{k,r}^{+} \right)^{-1} \right)^2 \right] + \mathcal{O}(1), \quad (C.3) $$

where $\mathcal{D}_{\lambda,k} = \lambda \mathcal{I}_k - G_k$. This identity holds under the assumptions that: (i) $\det \mathcal{D}_k \neq 0$ for all $k$, and (ii) $\det \mathcal{D}_k$ is a piecewise continuous function with jumps at $k = k_1, \dots, k_R$. The boundary term can also be obtained (see, for instance [55]), but it is not relevant for our purposes.

Combining Eq. (C.1) with Eq. (C.3) the entanglement entropy is given as

$$ S_{\ell} = a_1 \ell + a_2 \ln \ell + \mathcal{O}(1), \quad (C.4) $$

with

$$ a_1 = \lim_{\epsilon \to 0^+} \oint_{\Gamma} \frac{s(\epsilon, \lambda)}{4\pi^2 i} \int_{-\pi}^{+\pi} \frac{d}{d\lambda} \ln \det \mathcal{D}_{\lambda,k} dk d\lambda, \quad (C.5) $$

$$ a_2 = \lim_{\epsilon \to 0^+} \oint_{\Gamma} \frac{s(\epsilon, \lambda)}{8\pi^3 i} \sum_{r=1}^{R} \frac{d}{d\lambda} \mathrm{Tr} \left[ \ln \left( \mathcal{D}_{k,r}^{-} \left( \mathcal{D}_{k,r}^{+} \right)^{-1} \right)^2 \right] d\lambda. \quad (C.6) $$

## C.1 The Hatano-Nelson model

In this subsection, we calculate the steady-state entanglement entropy for an initial state in class B and C evolving under the Hatano-Nelson model. In these cases, the steady-state correlation matrix is diagonal in the momentum space,

$$ \langle c_k^\dagger c_q \rangle_{ss} = \begin{cases} \delta_{kq}, & k \in [\pi/2 - k_f, \pi/2 + k_f], \\ 0, & \text{otherwise,} \end{cases} \quad (C.7) $$

where $k_f$ is defined by the total particle number of the initial state in class B, while $k_f = \pi/2$ for the states in class C.

The determinant of the matrix $\mathcal{D}_{\lambda,k}$ is discontinuous in $\pi/2 - k_f$ and $\pi/2 + k_f$. Using the Fisher-Hartwig conjecture, we see that there is no volume term, since

$$ a_1 = \lim_{\epsilon \to 0^+} \oint_{\Gamma} \frac{s(\epsilon, \lambda)}{4\pi^2 i} \left[ \int_A dk \frac{1}{\lambda - 1} + \int_{\bar{A}} dk \frac{1}{\lambda} \right] = 0, \quad (C.8) $$

where $A = [\pi/2 - k_f, \pi/2 + k_f]$, and $\bar{A}$ the complement of $A$ in $[-\pi, \pi]$. In contrast, the coefficient of the logarithmic contribution is non-null and is given by

$$ a_2 = -\frac{1}{4\pi^3 i} \lim_{\epsilon \to 0^+} \oint_{\Gamma} \frac{d}{d\lambda} s(\epsilon, \lambda) \ln \left( \frac{\lambda}{\lambda-1} \right)^2 d\lambda = \frac{1}{\pi^2} \int_0^1 \ln \left( \frac{1-\lambda}{\lambda} \right)^2 d\lambda = \frac{1}{3}. \quad (C.9) $$

## C.2 The SSH chain

In this subsection, we focus on the non-Hermitian SSH model. We describe the main steps involved in the analytical calculation of the steady-state entanglement entropy for the system parameters for which this is possible in the three different classes.
### C.2.1 Class A

As demonstrated, the expectation value of the operator $\hat{n}_k$ is conserved if and only if the initial state is an eigenstate of the operator. With this in mind, it is possible to explicitly show that the equations of motion satisfied by the two-point correlation matrix factorize for each momentum label. This makes the problem analytically tractable, as the dynamic is confined to the orbital degrees of freedom within each k sector, which in this case are only two. To make it concrete, we assume that the initial state is defined as Eq. (41) of the main text with a filling $v$. We start by defining the following basis that simplifies the calculations,

$$ \begin{pmatrix} d_{k,+} \\ d_{k,-} \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -e^{ik/2} \\ e^{-ik/2} & 1 \end{pmatrix} \begin{pmatrix} c_{A,k} \\ c_{B,k} \end{pmatrix}, \quad g_k = \begin{pmatrix} \langle d_{k,+}^{\dagger} d_{k,+} \rangle \\ \langle d_{k,-}^{\dagger} d_{k,+} \rangle \\ \langle d_{k,-}^{\dagger} d_{k,-} \rangle \end{pmatrix}. \quad (C.10) $$

Depending on the eigenvalue $\hat{n}_k = \sum_\alpha d_{k,\alpha}^\dagger d_{k,\alpha}$, which shows if the k-mode starts empty, singly occupied, or doubly occupied, we have three initial conditions: $g_k = \mathbb{0}_{2 \times 2}$, $g_k = (\mathcal{I}_{2 \times 2} - \sigma_z)/2$, or $g_k = \mathcal{I}_{2 \times 2}$, respectively. Where $\mathbb{0}_{2 \times 2}$ is the two by two zero matrix and $\sigma_z$ the z Pauli matrix. After some lengthy algebra using the expression Eq. (7), the dynamics in this basis can be integrated and then rotated back to get the two-point function for a specific k in the original basis [27],

$$ G_k = \begin{pmatrix} \langle c_{k,A}^{\dagger} c_{k,A} \rangle & \langle c_{k,A}^{\dagger} c_{k,B} \rangle \\ \langle c_{k,B}^{\dagger} c_{k,A} \rangle & \langle c_{k,B}^{\dagger} c_{k,B} \rangle \end{pmatrix}. \quad (C.11) $$

The average value in the steady state is obtained by using a stationary phase approximation, which gives the following results:

* For $\langle \hat{n}_k \rangle = 0$, $G_k = \mathbb{0}_{2 \times 2}$.

* For $\langle \hat{n}_k \rangle = 1$, the correlation matrix is nontrivial and given if $\varepsilon_k^2 \ge 0$ by

$$ G_k = \frac{1}{2} \begin{pmatrix} 1 & e^{-ik/2} \\ e^{ik/2} & 1 \end{pmatrix} + \begin{pmatrix} 0 & e^{-ik/2}(a_k + 2Jib_k)\chi_k \\ e^{ik/2}(a_k - 2Jib_k)\chi_k & 0 \end{pmatrix}, \quad (C.12) $$

where

$$ \begin{aligned} \chi_k &= \frac{A_k}{(1+a_k)A_k} \left( 1 - \frac{1}{\sqrt{2(1+a_k)A_k + 1}} \right), & A_k &= \frac{\gamma^2 - h^2 \sin(k/2)^2}{2\varepsilon_k^2}, \\ a_k &= \frac{\gamma + h \sin(k/2)}{\gamma - h \sin(k/2)}, & b_k &= \frac{\cos(k/2)}{\gamma - h \sin(k/2)}. \end{aligned} $$

Otherwise if $\varepsilon_k^2 < 0$, we get

$$ G_k = \frac{1}{2} \begin{pmatrix} 1 & e^{-ik/2} \\ e^{ik/2} & 1 \end{pmatrix} + \frac{1}{2\gamma} \begin{pmatrix} i|\varepsilon_k| & e^{-ik/2}c_k \\ e^{-ik/2}c_k^* & -i|\varepsilon_k| \end{pmatrix}, \quad (C.13) $$

with $c_k = \gamma + h \sin(k/2) + 2iJ \cos(k/2)$.

* For $\langle \hat{n}_k \rangle = 2$, $G_k = \mathcal{I}_{2 \times 2}$.

Then we can use the formula Eq. (C.4) which depends on the eigenvalues $\mu_{k,\pm}$ of $G_k$ to calculate the coefficient $a_1$ and $a_2$. For $\langle \hat{n}_k \rangle = 0$ or 2, eigenvalues are trivial from $G_k$, in the case of $\langle \hat{n}_k \rangle = 1$ we obtain

$$ \left\{ \begin{array}{ll} \mu_{k,\pm} = \frac{1 \pm \xi_k}{2} = \frac{1 \pm \sqrt{1 - 4a_k \chi_k + 4(a_k^2 + 4b_k^2)\chi_k^2}}{2}, & \text{if } \varepsilon_k^2 \ge 0, \\ \mu_{k,+} = 1, \quad \mu_{k,-} = 0, & \text{if } \varepsilon_k^2 < 0. \end{array} \right. \quad (C.14) $$
<figure><img src="image_16.png" /><figcaption>Figure 14: Examples of functions $\mu_{k,\pm}$ for $\nu < 1/2$ and $\nu > 1/2$ - a) Case $\nu < 1/2$ where $k$ modes are either empty or singly occupied. Discontinuities are in $k_{f,1} = k_0 - 2\pi\nu$ and $k_{f,2} = k_0 + 2\pi\nu$ and functions are plotted for $\nu = 0.2$. b) Case $\nu > 1/2$ $k$ modes are either singly or doubly occupied. Discontinuities are in $k_{f,1} = k_0 - 2\pi\nu$ and $k_{f,2} = k_0 + 2\pi\nu$, functions are plotted for $\nu = 0.75$. Other parameters $\gamma = 1.7J$, and $k_0 = 1/2$.</figcaption></figure>

We also define

$$k_{f,1} = k_0 - 2\pi\nu, \quad \text{and} \quad k_{f,2} = k_0 + 2\pi\nu, \qquad (C.15)$$

the borders of the fermi sea of the initial state we consider which are depicted in Fig. 14.

Thus $a_1$ is simply computed as shown in the main text, with modes empty or doubly occupied that do not contribute to the linear prefactor. We rewrite the result here for ease,

$$a_1 = \int_V \frac{dk}{2\pi} [s(\mu_{k,+}) + s(\mu_{k,-})], \quad V = [-k_*, k_*] \cap [k_{f,1}, k_{f,2}]. \qquad (C.16)$$

where $s(x) = s(0, x)$ and which directly comes from Eq. (C.6).

To compute $a_2$, we need to understand the discontinuities of $G_k$ depending on initial filling $\nu$ and $\gamma$.

First, we describe the situation for $\nu < 1/2$, in this case we have only empty or singly occupied $k$ modes, as shown in Fig. 14 a). We obtain the following two cases:

* If $k_{f,1/2} \in [-\pi, -k_*] \cup [k_*, \pi]$ (i.e. in the imaginary part of the spectrum) the discontinuity comes from $\mu_{k,\pm} = 0$ (empty modes) jumping to $\mu_{k,\pm} = 0, 1$ (singly filled modes in the imaginary part of the spectrum). In that case, only one of the eigenvalues is discontinuous and contributes to $a_2$ in Eq. (C.6).

* If $k_{f,1/2} \in [-k_*, k_*]$ (i.e. in the real part of the spectrum) the discontinuity comes from $\mu_{k,\pm} = 0$ (empty modes) jumping to $\mu_{k,\pm} = (1 \pm \xi_k)/2$ (singly filled modes in the real part of the spectrum). Both eigenvalues are discontinuous.

We point out that these are the only discontinuities because for singly filled modes,

$$\mu_{k,\pm}[\varepsilon_k^2 = 0^+] = \mu_{k,\pm}[\varepsilon_k^2 = 0^-]. \qquad (C.17)$$

We now describe the situation for $\nu > 1/2$ as shown in Fig. 14 b). The reasoning is similar, we have two cases:

* If $k_{f,1/2} \in [-\pi, -k_*] \cup [k_*, \pi]$ (i.e. in the imaginary part of the spectrum) the discontinuity comes from $\mu_{k,\pm} = 1$ (doubly occupied modes) jumping to $\mu_{k,\pm} = 0, 1$ (singly filled modes in the imaginary part of the spectrum). In that case, only one of the eigenvalues is discontinuous and contributes to $a_2$ in Eq. (C.6).
• If $k_{f,1/2} \in [-k_*, k_*]$ (i.e. in the real part of the spectrum) the discontinuity comes from $\mu_{k,\pm} = 1$ (doubly occupied modes) jumping to $\mu_{k,\pm} = (1 \pm \xi_k)/2$ (singly filled modes in the real part of the spectrum). Both eigenvalues are discontinuous.

We then obtain $a_2$ through Eq. (C.6) using the discussed discontinuities in matrices $\overline{D}_{k,r}^+$ ($\overline{D}_{k,r}^-$) for the left (right) limit around the discontinuity.<fnref n="3" />

### C.2.2 Class B

For an initial state in class B we can compute the entanglement entropy rigorously in the PT fully broken phase, where Eq. (18) holds. As such, the steady state is a Slater determinant, constructed by filling the single-particle states with the largest imaginary eigenvalues,

$$|\Psi_{ss}\rangle = \frac{1}{\mathcal{N}} \prod_{|k-\pi|<2\pi\nu} \tilde{B}_{k,+}^{\dagger} |\text{vac}\rangle, \qquad (C.18)$$

where $\mathcal{N}$ is a normalization factor and,

$$\tilde{F}_k^\dagger = \begin{cases} \tilde{f}_{k+}^\dagger & (\text{single occupancy}), \\ \tilde{f}_{k-}^\dagger \tilde{f}_{k+}^\dagger & (\text{double occupancy}). \end{cases} \qquad (C.19)$$

These fermionic modes appear when diagonalizing the Hamiltonian (Eq. (37)) in a biortogonal basis,

$$\mathcal{H}_{\text{nH}} = \sum_k (\varepsilon_{k,+} \tilde{f}_{k+}^\dagger f_{k+} - \varepsilon_{k,-} \tilde{f}_{k-}^\dagger f_{k-}), \qquad (C.20)$$

where $f_{k\alpha}$ annihilates a fermion in the left eigenstate, while $\tilde{f}_{k\alpha}^\dagger$ creates a fermion in the right eigenstate of the single-particle Hamiltonian. These operators obey unconventional anticommutation relations:

$$\{f_{q\alpha}, \tilde{f}_{k\beta}^{\dagger}\} = \delta_{k,q} \delta_{\alpha,\beta}, \quad \{f_{q\alpha}, f_{k\beta}^{\dagger}\} = \delta_{k,q} V_{\alpha,\beta}, \quad \{\tilde{f}_{q\alpha}, \tilde{f}_{k\beta}^{\dagger}\} = \delta_{k,q} P_{\alpha,\beta}, \qquad (C.21)$$

with $V_{\alpha,\beta}$ and $P_{\alpha,\beta}$ determined explicitly by the microscopic parameters of the Hamiltonian.

Depending on whether the correlation matrix is singly or doubly occupied, it has the eigenvalues:

$$\mu_{k,\pm} = \begin{cases} 0 \text{ and } 0 & (\text{empty}), \\ 0 \text{ and } 1 & (\text{single occupancy}), \\ 1 \text{ and } 1 & (\text{double occupancy}). \end{cases} \qquad (C.22)$$

As explained in the previous subsection, this implies that $a_1 = 0$ and $a_2 = 1/3$.

We note from the numerical results that in the PT-Mixed phase, the steady-state is also a Slater determinant constructed from the single-particle modes with the highest imaginary energy, whenever the total number of particles is less than the number of purely imaginary eigenmodes of the Hamiltonian. In this case, we can also apply Segzö lemma and the Fisher-Hartwig conjecture, obtaining $a_1 = 0$ and $a_2 = 1/3$.

We stress that we cannot analytically prove that the steady-state is a Slater determinant constructed from the single-particle modes with the highest imaginary energy. In this case, we cannot apply Eq. (18) as the Hamiltonian contains two values of the quasimomentum exceptional points where it cannot be diagonalized explicitly:

$$\mathcal{H}_{\text{nH}} = \sum_{k \in [-\pi, \pi] \setminus \pm k_{*}, \alpha} \varepsilon_{k,\alpha} \tilde{f}_{k\alpha}^{\dagger} f_{k\alpha} + \sum_{k \in \pm k_{*}, \alpha, \beta} h_{k}^{\alpha\beta} c_{k\alpha}^{\dagger} c_{k\beta}. \qquad (C.23)$$

<footnote n="3">³We recall that the $\pm$ symbol in $\mu_k$ is due a unit cell of 2 sites, A and B, whereas $\pm$ in $\overline{D}_{k,r}$ indicates the left and right limits around a discontinuity.</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 15: Steady-state entanglement entropy as a function of the subsystem fraction $ℓ/L$ in the volume-law phase for class B (a) and class C (b). Other parameters: $v = 1/4$, $v_0 = 1/4$, $\gamma = 0.5J$, and $h = J$.</figcaption></figure>

### C.2.3 Class C

In the PT-Fully Broken phase, the steady-state correspondent to an initial state in class C is given by filling all the single-particle states with positive imaginary eigenvalues,

$$|\Psi_{ss}\rangle = \prod_{k \in [-\pi, \pi)} \tilde{f}_{k,+}^{\dagger} |\text{vac}\rangle. \qquad (C.24)$$

Following the same arguments as in the previous subsection, the entanglement entropy follows an area law since the eigenvalues of the correlation matrix for all $k$ are either zero or one, and the determinant of the correlation matrix is continuous. Given this $a_1 = a_2 = 0$.

## D Scaling with subsystem size in the volume law phase

Here, we present additional results characterizing entanglement in the volume-law phase. Specifically, we analyze the dependence of the steady-state entanglement entropy on the subsystem size $ℓ$ for finite chains of total size $L$, as shown in Fig. 15. The results obtained for both class B and class C indicate that the entanglement entropy follows the scaling law

$$S(\ell/L) = L \cdot f(\ell/L) + b_0, \qquad (D.1)$$

where $f(\cdot)$ is an $O(1)$ function that depends only on the ratio $ℓ/L$ (and on the system parameters), and does not scale with the total system size. The term $b_0$ represents subleading contributions; these are negligible in the limit $L \gg 1$, and only affect smaller chains. Thus, the leading linear scaling of the entanglement entropy with total system size is independent of the specific value of $ℓ$ chosen.

## References

[1] A. Polkovnikov, K. Sengupta, A. Silva and M. Vengalattore, *Colloquium: Nonequilibrium dynamics of closed interacting quantum systems*, Rev. Mod. Phys. **83**, 863 (2011), doi:10.1103/RevModPhys.83.863 [preprint doi:10.48550/arXiv.1007.5331].
[2] C. Gogolin and J. Eisert, *Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems*, Rep. Prog. Phys. **79**, 056001 (2016), doi:10.1088/0034-4885/79/5/056001.

[3] P Calabrese and J. Cardy, *Evolution of entanglement entropy in one-dimensional systems*, J. Stat. Mech.: Theory Exp. P04010 (2005), doi:10.1088/1742-5468/2005/04/P04010.

[4] H. Kim and D. A. Huse, *Ballistic spreading of entanglement in a diffusive nonintegrable system*, Phys. Rev. Lett. **111**, 127205 (2013), doi:10.1103/PhysRevLett.111.127205 [preprint doi:10.48550/arXiv.1306.4306].

[5] M. Rigol, V. Dunjko and M. Olshanii, *Thermalization and its mechanism for generic isolated quantum systems*, Nature **452**, 854 (2008), doi:10.1038/nature06838 [preprint doi:10.48550/arXiv.0708.1324].

[6] R. Fazio, J. Keeling, L. Mazza and M. Schirò, *Many-body open quantum systems*, SciPost Phys. Lect. Notes 99 (2025), doi:10.21468/SciPostPhysLectNotes.99 [preprint doi:10.48550/arXiv.2409.10300].

[7] Y. Ashida, Z. Gong and M. Ueda, *Non-Hermitian physics*, Adv. Phys. **69**, 249 (2020), doi:10.1080/00018732.2021.1876991 [preprint doi:10.48550/arXiv.2006.01837].

[8] J. Dalibard, Y. Castin and K. Mølmer, *Wave-function approach to dissipative processes in quantum optics*, Phys. Rev. Lett. **68**, 580 (1992), doi:10.1103/PhysRevLett.68.580.

[9] A. J. Daley, *Quantum trajectories and open many-body quantum systems*, Adv. Phys. **63**, 77 (2014), doi:10.1080/00018732.2014.933502 [preprint doi:10.48550/arXiv.1405.6694].

[10] H. M. Wiseman and G. J. Milburn, *Quantum measurement and control*, Cambridge University Press, Cambridge, UK, ISBN 9780521804424 (2009), doi:10.1017/CBO9780511813948.

[11] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter and D. N. Christodoulides, *Non-Hermitian physics and PT symmetry*, Nat. Phys. **14**, 11 (2018), doi:10.1038/nphys4323.

[12] M.-A. Miri and A. Alù, *Exceptional points in optics and photonics*, Science **363**, eaar7709 (2019), doi:10.1126/science.aar7709.

[13] A. McDonald, R. Hanai and A. A. Clerk, *Nonequilibrium stationary states of quantum non-Hermitian lattice models*, Phys. Rev. B **105**, 064302 (2022), doi:10.1103/PhysRevB.105.064302 [preprint doi:10.48550/arXiv.2103.01941].

[14] 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), doi:10.1103/PhysRevX.8.031079 [preprint doi:10.48550/arXiv.1802.07964].

[15] L. Herviou, N. Regnault and J. H. Bardarson, *Entanglement spectrum and symmetries in non-Hermitian fermionic non-interacting models*, SciPost Phys. **7**, 069 (2019), doi:10.21468/SciPostPhys.7.5.069 [preprint doi:10.48550/arXiv.1908.09852].

[16] E. J. Bergholtz, J. C. Budich and F. K. Kunst, *Exceptional topology of non-Hermitian systems*, Rev. Mod. Phys. **93**, 015005 (2019), doi:10.1103/RevModPhys.93.015005 [preprint doi:10.48550/arXiv.1912.10048].
[17] Y. Ashida and M. Ueda, Full-counting many-particle dynamics: Nonlocal and chiral propagation of correlations, Phys. Rev. Lett. **120**, 185301 (2018), doi:10.1103/PhysRevLett.120.185301 [preprint doi:10.48550/arXiv:1709.03704].

[18] 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), doi:10.1103/PhysRevResearch.2.033069 [preprint doi:10.48550/arXiv:1909.01346].

[19] G. Cipolloni and J. Kudler-Flam, Entanglement entropy of non-Hermitian eigenstates and the Ginibre ensemble, Phys. Rev. Lett. **130**, 010401 (2023), doi:10.1103/PhysRevLett.130.010401 [preprint doi:10.48550/arXiv:2206.12438].

[20] S. S. Roy, S. Bandyopadhyay, R. Costa de Almeida and P. Hauke, Unveiling eigen-state thermalization for non-Hermitian systems, Phys. Rev. Lett. **134**, 180405 (2025), doi:10.1103/PhysRevLett.134.180405 [preprint doi:10.48550/arXiv:2309.00049].

[21] G. Cipolloni and J. Kudler-Flam, Non-Hermitian Hamiltonians violate the eigenstate thermalization hypothesis, Phys. Rev. B **109**, L020201 (2024), doi:10.1103/PhysRevB.109.L020201 [preprint doi:10.48550/arXiv:2303.03448].

[22] K. Kawabata, T. Numasawa and S. Ryu, Entanglement phase transition induced by the non-Hermitian skin effect, Phys. Rev. X **13**, 021007 (2023), doi:10.1103/PhysRevX.13.021007 [preprint doi:10.48550/arXiv:2206.05384].

[23] T. Orito and K.-I. Imura, Entanglement dynamics in the many-body Hatano-Nelson model, Phys. Rev. B **108**, 214308 (2023), doi:10.1103/PhysRevB.108.214308 [preprint doi:10.48550/arXiv:2308.03078].

[24] S. Gopalakrishnan and M. J. Gullans, Entanglement and purification transitions in non-Hermitian quantum mechanics, Phys. Rev. Lett. **126**, 170503 (2021), doi:10.1103/PhysRevLett.126.170503 [preprint doi:10.48550/arXiv:2012.01435].

[25] S.-K. Jian, Z.-C. Yang, Z. Bi and X. Chen, Yang-Lee edge singularity triggered entanglement transition, Phys. Rev. B **104**, L161107 (2021), doi:10.1103/PhysRevB.104.L161107 [preprint doi:10.48550/arXiv:2101.04115].

[26] X. Turkeshi and M. Schiró, Entanglement and correlation spreading in non-Hermitian spin chains, Phys. Rev. B **107**, L020403 (2023), doi:10.1103/PhysRevB.107.L020403 [preprint doi:10.48550/arXiv:2201.09895].

[27] Y. Le Gal, X. Turkeshi and M. Schirò, Volume-to-area law entanglement transition in a non-Hermitian free fermionic chain, SciPost Phys. **14**, 138 (2023), doi:10.21468/SciPostPhys.14.5.138 [preprint doi:10.48550/arXiv:2210.11937].

[28] C. Zerba and A. Silva, Measurement phase transitions in the no-click limit as quantum phase transitions of a non-hermitean vacuum, SciPost Phys. Core **6**, 051 (2023), doi:10.21468/SciPostPhysCore.6.3.051.

[29] 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 (2023), doi:10.1103/PhysRevLett.130.230401 [preprint doi:10.48550/arXiv:2212.10584].

[30] L. Su, A. Clerk and I. Martin, Dynamics and phases of nonunitary Floquet transverse-field Ising model, Phys. Rev. Res. **6**, 013131 (2024), doi:10.1103/PhysRevResearch.6.013131 [preprint doi:10.48550/arXiv:2306.07428].
[31] L. Vidmar and M. Rigol, Generalized Gibbs ensemble in integrable lattice models, J. Stat. Mech.: Theory Exp. 064007 (2016), doi:10.1088/1742-5468/2016/06/064007.

[32] E. Starchl and L. M. Sieberer, Relaxation to a parity-time symmetric generalized Gibbs ensemble after a quantum quench in a driven-dissipative Kitaev chain, Phys. Rev. Lett. 129, 220602 (2022), doi:10.1103/PhysRevLett.129.220602 [preprint doi:10.48550/arXiv.2203.14589].

[33] E. Starchl and L. M. Sieberer, Quantum quenches in driven-dissipative quadratic fermionic systems with parity-time symmetry, Phys. Rev. Res. 6, 013016 (2024), doi:10.1103/PhysRevResearch.6.013016 [preprint doi:10.48550/arXiv.2304.01836].

[34] C. M. Bender and D. W. Hook, *PT*-symmetric quantum mechanics, Rev. Mod. Phys. **96**, 045002 (2024), doi:10.1103/RevModPhys.96.045002.

[35] N. Hatano and D. R. Nelson, Localization transitions in non-Hermitian quantum mechanics, Phys. Rev. Lett. **77**, 570 (1996), doi:10.1103/physrevlett.77.570 [preprint doi:10.48550/arXiv.cond-mat/9603165].

[36] N. Hatano and D. R. Nelson, Vortex pinning and non-Hermitian quantum mechanics, Phys. Rev. B **56**, 8651 (1997), doi:10.1103/physrevb.56.8651 [preprint doi:10.48550/arXiv.cond-mat/9705290].

[37] R. Diogo Soares and M. Schirò, Non-unitary quantum many-body dynamics using the Faber polynomial method, SciPost Phys. **17**, 128 (2024), doi:10.21468/SciPostPhys.17.5.128 [preprint doi:10.48550/arXiv.2406.10135].

[38] K. Jacobs, Quantum measurement theory and its applications, Cambridge University Press, Cambridge, UK, ISBN 9781107025486 (2014), doi:10.1017/CBO9781139179027.

[39] P Calabrese and J. Cardy, Entanglement entropy and quantum field theory, J. Stat. Mech.: Theory Exp. P06002 (2004), doi:10.1088/1742-5468/2004/06/P06002.

[40] L. Amico, R. Fazio, A. Osterloh and V. Vedral, Entanglement in many-body systems, Rev. Mod. Phys. **80**, 517 (2008), doi:10.1103/RevModPhys.80.517.

[41] I. Peschel and V. Eisler, Reduced density matrices and entanglement entropy in free lattice models, J. Phys. A: Math. Theor. **42**, 504003 (2009), doi:10.1088/1751-8113/42/50/504003.

[42] B.-Q. Jin and V. E. Korepin, Quantum spin chain, Toeplitz determinants and the Fisher-Hartwig conjecture, J. Stat. Phys. **116**, 79 (2004), doi:10.1023/b:joss.0000037230.37166.42.

[43] D. A. Ivanov, A. G. Abanov and V. V. Cheianov, Counting free fermions on a line: A Fisher-Hartwig asymptotic expansion for the Toeplitz determinant in the double-scaling limit, J. Phys. A: Math. Theor. **46**, 085003 (2013), doi:10.1088/1751-8113/46/8/085003.

[44] S. Groha, F. Essler and P. Calabrese, Full counting statistics in the transverse field Ising chain, SciPost Phys. **4**, 043 (2018), doi:10.21468/SciPostPhys.4.6.043 [preprint doi:10.48550/arXiv.1803.09755].

[45] M. B. Hastings, An area law for one-dimensional quantum systems, J. Stat. Mech.: Theory Exp. P08024 (2007), doi:10.1088/1742-5468/2007/08/P08024.
[46] M. M. Wolf, Violation of the entropic area law for fermions, Phys. Rev. Lett. **96**, 010404 (2006), doi:**10.1103/PhysRevLett.96.010404** [preprint doi:**10.48550/arXiv.quant-ph/0503219**].

[47] P Di Francesco, P. Mathieu and D. Sénéchal, Conformal field theory, Springer, New York, USA, ISBN 9781461274759 (1997), doi:10.1007/978-1-4612-2256-9.

[48] Á. Bácsi and B. Dóra, Dynamics of entanglement after exceptional quantum quench, Phys. Rev. B **103**, 085137 (2021), doi:**10.1103/PhysRevB.103.085137** [preprint doi:**10.48550/arXiv.2011.11979**].

[49] N. Chakrabarti, N. Nirbhan and A. Bhattacharyya, Dynamics of monitored SSH model in Krylov space: From complexity to quantum Fisher information, J. High Energy Phys. **07**, 203 (2025), doi:**10.1007/JHEP07(2025)203** [preprint doi:**10.48550/arXiv.2502.03434**].

[50] E. Medina-Guerra, I. V. Gorny i and Y. Gefen, Phase transitions in a non-Hermitian Su-Schrieffer-Heeger model via Krylov spread complexity, Phys. Rev. B **112**, 035427 (2025), doi:**10.1103/6lvg-7qdn** [preprint doi:**10.48550/arXiv.2503.18936**].

[51] Y. Le Gal, X. Turkeshi and M. Schirò, Entanglement dynamics in monitored systems and the role of quantum jumps, PRX Quantum **5**, 030329 (2024), doi:**10.1103/PRXQuantum.5.030329** [preprint doi:**10.48550/arXiv.2312.13419**].

[52] G. Di Giulio, X. Turkeshi and S. Murciano, Measurement-induced symmetry restoration and quantum Mpemba effect, Entropy **27**, 407 (2025), doi:**10.3390/e27040407**.

[53] U. Agrawal, A. Zabalo, K. Chen, J. H. Wilson, A. C. Potter, J. H. Pixley, S. Gopalakrishnan and R. Vasseur, Entanglement and Charge-sharpening transitions in $U(1)$ symmetric monitored quantum circuits, Phys. Rev. X **12**, 041002 (2022), doi:**10.1103/PhysRevX.12.041002** [preprint doi:**10.48550/arXiv.2107.10279**].

[54] R. D. Soares, Y. Le Gal and M. Schirò, Entanglement transition due to particle losses in a monitored fermionic chain, Phys. Rev. B **111**, 064313 (2025), doi:**10.1103/PhysRevB.111.064313** [preprint doi:**10.48550/arXiv.2408.03700**].

[55] P. Calabrese and F. H. L. Essler, Universal corrections to scaling for block entanglement in spin-1/2 XX chains, J. Stat. Mech.: Theory Exp. P08029 (2010), doi:**10.1088/1742-5468/2010/08/p08029**.