Entanglement Spectrum Crossings Reveal non-Hermitian Dynamical Topology

The development of non-Hermitian topological band theory has led to observations of novel topological phenomena in effectively classical, driven and dissipative systems. However, for open quantum many-body systems, the absence of a ground state presents a challenge to define robust signatures of non-Hermitian topology. We show that such a signature is provided by crossings in the time evolution of the entanglement spectrum. These crossings occur in quenches from the trivial to the topological phase of a driven-dissipative Kitaev chain that is described by a Markovian quantum master equation in Lindblad form. At the topological transition, which can be crossed either by changing parameters of the Hamiltonian of the system or by increasing the strength of dissipation, the time scale at which the first entanglement spectrum crossing occurs diverges with a dynamical critical exponent of $\epsilon = 1/2$. We corroborate these numerical findings with an exact analytical solution of the quench dynamics for a spectrally flat postquench Liouvillian. This exact solution suggests an interpretation of the topological quench dynamics as a fermion parity pump. Our work thus reveals signatures of non-Hermitian topology which are unique to quantum many-body systems and cannot be emulated in classical simulators of non-Hermitian wave physics.


I. INTRODUCTION
Non-Hermitian topological band theory [1][2][3][4][5][6][7][8][9][10][11][12][13] allows to derive a topological classification of the complex spectra of quadratic Liouvillians [14][15][16][17][18][19][20], which describe the dynamics of noninteracting, driven and open quantum many-body systems. For complex spectra, the very notions of a ground state, as well as of filled and empty bands, are ill-defined. In the topological band theory for isolated and noninteracting fermionic systems, however, these concepts are central to predict robust topological features, such as the quantized response of a system in its ground state [21][22][23]. Moreover, the existence of a ground state allows to generalize topological phases to interacting many-body systems [23]. These observations raise the question that motivates our work: Which robust signatures of non-Hermitian topology survive in driven and open quantum many-body systems? Here we propose that the time evolution of the entanglement spectrum can provide such a signature.
The entanglement spectrum is a powerful tool to detect the presence of topological order in the pure ground state |ψ 0 of the Hamiltonian H 0 of an isolated Hermitian system [24][25][26][27][28][29][30][31][32]. For a subsystem A and its complement A c , the many-body entanglement spectrum is defined as the spectrum of the reduced density matrix ρ A = tr A c (ρ 0 ), where ρ 0 = |ψ 0 ψ 0 | is the density matrix that is associated with the pure state |ψ 0 . The bulk-edge correspondence of the entanglement spectrum [29,30] establishes a one-to-one correspondence between the low-energy edge states of the physical Hamiltonian at an open boundary of subsystem A, and the entanglement eigenstates of ρ A that have the highest entanglement eigenvalues, i.e., that are most entangled with the rest of the system. This correspondence holds both for noninteracting and interacting systems, and has even been applied to systems which are driven out of equilibrium: After a sudden quench from H 0 to H, the time-evolved state |ψ(t) = e −iHt |ψ 0 is the ground state of the time-dependent parent Hamiltonian H p (t) = e −iHt H 0 e iHt , and the entanglement spectrum of |ψ(t) reflects the topology of H p (t) [33][34][35][36].
Despite the strikingly universal applicability of the entanglement spectrum as a diagnostic for topological order, it is unclear whether the non-Hermitian topology of the complex spectrum of a LiouvillianL leaves signatures in the entanglement spectrum of an associated quantum state. For example, the time-evolved density matrix ρ(t) = eL t ρ 0 is, in general, a mixed state. Therefore, it cannot be interpreted as the ground state of any putative non-Hermitian parent Hamiltonian H(t). Intriguingly, establishing the connection between the entanglement spectrum of ρ(t) and the topology ofL could open up the possibility to characterize the non-Hermitian topology of both noninteracting and interacting open quantum many-body systems.
In this paper, we show that non-Hermitian dynamical topology of driven and open quantum many-body systems is revealed by crossings in the time evolution of the entanglement spectrum after a quench. We consider a Kitaev chain, which is subject to Markovian gain and dissipation [37][38][39][40][41] as illustrated in Fig. 1. The topology of the corresponding Liouvillian is characterized by a non-Hermitian winding number W [1,8]. A quench from a trivial to a topological phase leads to stable crossings in the entanglement spectrum of the time-evolved density matrix ρ(t) as shown in Fig. 1(a) Figure 1. Quench dynamics of a driven-dissipative Kitaev chain. At t = 0, the system is prepared in the topologically trivial ground state of the Kitaev chain, where Majorana fermions (spheres) pair locally (gray bonds) at each site of the chain (blue ellipses). Then, at t = 0 + , the system is driven out of its initial state by an abrupt change of the chemical potential. Simultaneously, the chain is connected to Markovian reservoirs, which induce dissipation at a rate γ. (a) For quenches to the topological phase of the postquench Liouvillian, the entanglement spectrum exhibits oscillatory decay with repeated zero crossings (only a single pair of entanglement eigenvalues ±ξ is shown as blue lines in the figure). At each zero crossing, the fermion parity of the entanglement ground state switches between even and odd. Red lines show the evolution of ±ξ for a quench closer to the topological phase boundary at γc. For δγ = |γ − γc| → 0, the time t1 at which the first crossing occurs diverges as t1 ∼ δγ − with a universal critical exponent = 1/2. (b) For quenches to a trivial phase, the entanglement eigenvalues ±ξ decay without crossing zero. illustrated in Fig. 1(b), the dynamics after a quench to a trivial phase does not feature such crossings. We show that topological entanglement spectrum crossings can be traced back to the reversal of the fermion parity in individual entanglement eigenstates. Thus, our work establishes that entanglement spectrum crossings and the concomitant pumping of fermion parity are robust signatures of non-Hermitian dynamical topology of the Liouvillian.
Surprisingly, we obtain these results for systems which heat up to a featureless infinite-temperature state at late times. By contrast, much previous work focused on inducing nontrivial topological features in the steady state ρ ss = lim t→∞ ρ(t), which can be achieved by properly designing purely dissipative dynamics [42][43][44][45][46]. Adopting this point of view, one is led to the conclusion that the systems we study are trivial. Instead, the perspective taken in this work, which is inspired by non-Hermitian band theory [1][2][3][4][5][6][7][8][9][10][11][12][13], focuses on the topological properties of the Liouvillian. As explained above, these properties are revealed in the dynamics, i.e., in the approach to the steady state. The time evolution of the entanglement spectrum reveals a sharp dynamical topological phase transition with associated dynamical criticality, reflected in a divergence of the time scale on which entanglement spectrum crossings occur as illustrated in Fig. 1(a). These features highlight the stark contrast between the two complementary approaches to the topology of driven and open quantum many-body systems: the one based on ρ ss [42][43][44][45][46], and the one that we propose, which is based on the entanglement spectrum.
We finally emphasize that the characterization of topology in terms of entanglement spectrum crossings is unique to quantum many-body systems and cannot be emulated by photonic [47][48][49][50][51][52][53][54][55][56][57] or mechanical [58] simulators of non-Hermitian wave physics. Moreover, by using the entanglement spectrum as a tool to diagnose non-Hermitian topology, our work paves the way to characterize non-Hermitian dynamical topology also of interacting open quantum many-body systems.
This paper is organized as follows: We present our main results in Sec. II. The model and formalism on which our study is based are introduced in Sec. III, and we summarize the non-Hermitian topological band theory of the driven-dissipative Kitaev chain in Sec. IV. The relationship between entanglement spectrum crossings and non-Hermitian topology of the Liouvillian depends on whether the jump operators, which describe the coupling of the system to Markovian reservoirs, are Hermitian or non-Hermitian. Section V is devoted to systems with Hermitian jump operators. In particular, we discuss the time evolution of entanglement spectra, the connection between entanglement spectrum crossings and parity pumping, and dynamical criticality. We consider systems non-Hermitian jump operators in Sec. VI. Our conclusions are presented in Sec. VII, together with an outlook on future perspectives. Technical details are deferred to Appendices A-G.

II. SUMMARY OF MAIN RESULTS
We consider the following model and quench protocol: At time t = 0, a Kitaev chain is prepared in its ground state |ψ 0 with a chemical potential µ 0 , which is then changed abruptly to µ 1 at t = 0 + . Simultaneously, the system is coupled with strength γ to local Markovian reservoirs. The resultant driven-dissipative dynamics is described by a quantum master equation for the density matrix ρ of the system, dρ/dt =Lρ. Here, the non-Hermitian Liouville superoperatorL comprises both coherent dynamics induced by the Hamiltonian H of the system, and the coupling to reservoirs through quantum jump operators. We study the time evolution of the entanglement spectrum, i.e., the spectrum of the reduced density matrix ρ A (t) = tr A c (ρ(t)), where ρ(t) = eL t ρ 0 is the state of the system at time t, and ρ 0 = |ψ 0 ψ 0 | is the initial state.
We first discuss systems with Hermitian jump operators, for which the steady state is the trivial fully mixed infinite-temperature state ρ ss = ρ ∞ = 1/D, where D is the Hilbert-space dimension. For such systems, our main results are as follows: Figure 2. Dynamical topological phase diagram and dynamical criticality of entanglement spectrum crossings. The non-Hermitian winding number W (for phases with real line gaps) and the global Berry phase Q distinguish six dynamical phases: the real-line gapped phases with nontrivial (Rtop, with W = Q = 1) and trivial (Rtr, with W = Q = 0) topology, the imaginary-line gapped phases with (I edge , with Q = 1) and without (Itr, with Q = 0) edge states, and the gapless phases with (G edge , with Q = 1) and without (Gtr, with Q = 0) edge states. Entanglement spectrum crossings occur for quenches from the trivial phase of the isolated Kitaev chain to Rtop. The logarithmic color scale in Rtop encodes the time t1 of the first entanglement spectrum zero crossing for a quench with µ0 = −3J and in a system of size N = 100. Inset: t1 diverges at the boundary of Rtop with a critical exponent of = 1/2. The data shown corresponds to approaching the phase boundary along the black arrows labeled by δµ and δγ in the main panel: µ1 = −2J + δµ, γ = 0 for the blue dots, and µ1 = −J, γ = 0.5J − δγ for the orange dots. In the former case, longer evolution times can be reached due to the absence of decay. The system size is N = 1000. For comparison, the gray solid line shows a square-root singularity.
Entanglement spectrum crossings reveal non-Hermitian dynamical topology.
We find that entanglement spectrum crossings occur exclusively for quenches from the ground state of the Kitaev chain in the topologically trivial phase to the nontrivial phase of the postquench LiouvillianL, which is designated as R top in the dynamical topological phase diagram in Fig. 2. In R top , the two bands of complex eigenvalues of the Liouvillian are separated by a real line gap, and each band can be characterized by a non-Hermitian winding number that takes the value W = 1 [1,8]. Entanglement spectrum crossings do not occur for quenches to R tr with W = 0, the gapless phases G edge and G tr , and the phases I edge and I tr with an imaginary line gap. Therefore, the presence of entanglement spectrum crossings is directly connected to the nontrivial non-Hermitian topology of the Kitaev chain, and can be used to map out the sharp phase boundary of the topological phase R top both at zero and nonzero strengths of dissipation.
A second topological invariant, the global Berry phase Q that pertains to the full Liouvillian and not to individual bands [59,60], can be defined even in the gapless phases G edge and G tr . The global Berry phase takes the value Q = 1 and predicts the existence of edge modes of the Liouvillian in the phases G edge , R top , and I edge , while Q = 0 in the phases G tr , R tr , and I tr without edge modes. Thus, the existence of edge modes can be understood as a global property of the Liouvillian [60], while, as we show in our work, the presence of entanglement spectrum crossings requires nontrivial non-Hermitian topology of individual bands.
Dynamical criticality of the entanglement spectrum at the topological phase boundary is governed by a universal critical exponent = 1/2. At the phase boundary of R top , the dynamics of the entanglement spectrum exhibits critical behavior. This is shown in Fig. 2, where the color scale within R top encodes the time t 1 of the first entanglement spectrum crossing, which diverges at the phase boundary as t 1 ∼ δγ − and t 1 ∼ δµ − . Through an exact analytical solution of the quench dynamics for µ 1 = 0, where the spectrum of the Liouvillian is flat, we show that the entanglement spectrum dynamical critical exponent takes the exact value = 1/2 when the phase boundary of R top is approached along the direction that is indicated with a black arrow labeled by δγ x . Numerical evidence indicates that this value is universal along the phase boundary: In the inset of Fig. 2, we show data which corresponds to approaching the phase boundaries along the directions δµ and δγ, and is compatible with = 1/2. This value is also compatible with numerical results of Ref. [61], which studied entanglement dynamics in the transverse-field Ising model, to which the isolated Hermitian Kitaev chain maps under a Jordan-Wigner transformation. The universality of the value = 1/2 is corroborated further by results for a generalized drivendissipative Kitaev chain with quasilocal jump operators, which we present in Appendix A. These findings highlight the phase boundary of R top as a sharp dynamical topological transition at finite dissipation. The system evolves towards the same trivial infinite-temperature state for all values of µ 1 and γ, in contrast to the clear dynamical distinction between the different phases that we find and show in Fig. 2.
Entanglement spectrum crossings in the drivendissipative Kitaev chain can be interpreted as a fermion parity pump. Our exact solution of the quench dynamics for µ 1 = 0 shows that with each crossing of entanglement eigenvalues, the fermion parity of many-body entanglement eigenstates is reversed and, therefore, the topological quench dynamics can be interpreted as a fermion parity pump. Thus, the exact solution generalizes earlier results [35] for the isolated Kitaev chain to nonzero Markovian dissipation and non-Hermitian topology. The reversal of the fermion parity is a robust signature which might be easier to detect than entanglement spectrum crossings, in particular, in interacting systems.
The non-Hermitian topological phase transition is associated with a distinct type of dynamical criticality: many-body critical damping. The dynamical criticality of the entanglement spectrum at the phase boundary between R top and G edge raises the question: Do "simple" two-time correlation functions also become critical at this dynamical topological transition?
To address this question, we focus on the retarded response function in the steady state. In the non-Hermitian models we consider here, the retarded response function exhibits two distinct types of dynamical criticality. At the phase boundary between R top and G edge , the response function shows what we dub many-body critical damping. We choose this term in analogy to the paradigmatic classical damped harmonic oscillator, in which critical damping marks the onset of the overdamped regime. Similarly, as we increase the rate of dissipation γ to approach the phase boundary, the period of oscillations of the response function diverges. The divergence is governed by the same value of the exponent of 1/2 as for the time scale t 1 of the first crossing in the entanglement spectrum. However, the oscillatory behavior and associated zero crossings of the response function occur in both the topological and the trivial phases with real line gaps, R top and R tr , respectively. Therefore, these zero crossings are not related to topology.
The second and more conventional type of dynamical criticality we find is critical relaxation, which occurs at the boundaries between the gapless phases G edge and G tr , and the imaginary-line gapped phases I edge and I tr . Upon approaching these phase boundaries, the relaxation time scales of both the response function and the entanglement spectrum diverge, leading to power-law decay exactly on the critical line. In stark contrast, the response function and the entanglement spectrum decay exponentially on the line of many-body critical damping, i.e., the phase boundary between R top and G edge . As for a damped harmonic oscillator, the decay rate is largest at the onset of the overdamped regime.
The phase boundaries between R top and G edge as well as G edge and G tr coincide for γ → 0, when the model becomes Hermitian. Therefore, the occurrence of two distinct types of dynamical critical behavior-many-body critical damping and critical relaxation-is unique to non-Hermitian dynamics.
The above results and, in particular, the direct connection between entanglement spectrum crossings and nontrivial values of the non-Hermitian winding number W , pertain to systems with Hermitian jump operators. Our key finding for systems with non-Hermitian jump operators reads as follows: To reveal the non-Hermitian topology of general quadratic Liouvillians through entanglement spectrum crossings, the jump operators have to be Hermitianized. While systems with Hermitian jump operators evolve towards a trivial infinite-temperature steady state, for non-Hermitian jump operators, the steady state is in general nontrivial. In particular, a well-defined value of the fermion parity of the entanglement ground state of the steady state can lead to the occurrence of entanglement spectrum crossings which are not related to the topology of the Liouvillian. However, the direct connection between non-Hermitian topology ofL and entanglement spectrum crossings can be restored by means of Hermitianization of the jump operators. Hermitianization is a continuous deformation of the jump operators to make them Hermitian, while the gap and the symmetries of the Liouvillian are preserved [14].

III. MODEL
The Kitaev chain of length N is described by the Hamiltonian [62] (1) where c † i and c i are, respectively, creation and annihilation operators for spinless Dirac fermions on lattice site i. These operators obey the cannonical anticommutation We set c N +1 = 0 and c N +1 = c 1 for a system with open and periodic boundary conditions, respectively. The parameters in the Hamiltonian (1) are the hopping matrix element J, the p-wave pairing amplitude ∆, and the chemical potential µ. We measure energy and time in units of J and /J, respectively, and we set = 1.
In what follows, we choose ∆ = J ∈ R >0 . For this choice, the Kitaev chain belongs to the Altland-Zirnbauer class BDI [63], which, in one spatial dimension, is characterized by a winding number W ∈ Z [23]. The groundstate phase diagram of the Kitaev chain comprises two topologically distinct phases. In the topologically nontrivial phase, which is realized for |µ| < 2J, an infinite bulk system is characterized by a nonvanishing winding number W = 1. By virtue of the bulk-edge correspondence, a finite system with open boundary conditions features two Majorana zero-energy modes that are localized at the ends of the chain. These edge modes are absent in the trivial phase at |µ| > 2J for which W = 0.
A. Driven-dissipative Kitaev chain We are interested in the topological properties of a driven-dissipative Kitaev chain which is realized by connecting the isolated Kitaev chain, Eq. (1), to Markovian baths. For driven-dissipative systems, the notion of a ground state is ill-defined, and, therefore, we focus on dynamical signatures of topology. Thus, we proceed to specify the dynamics of the driven-dissipative Kitaev chain, and the types of Markovian baths we consider in this paper.
The time evolution of the density matrix ρ of the driven-dissipative Kitaev chain is described by a master equation in Lindblad form, In this equation, the Liouvillian superoperatorL comprises both unitary dynamics generated by the Hamiltonian (1), and Markovian drive and dissipation, which are incorporated in the dissipatorD. The latter is given bŷ where the operators L i , termed Lindblad or quantum jump operators, describe the coupling between the system and its environment.
In this paper, we study the dynamics of noninteracting open quantum many-body systems that are described by quadratic Liouvillians. Such Liouvillians are defined in terms of a quadratic Hamiltonian, together with jump operators L i that are linear combinations of the fermionic field operators c i and c † i . For simplicity, we focus in most of our work on the following purely local jump operators [41]: For γ g = 0, these jump operators describe loss of particles at a rate γ l > 0, and for γ l = 0, particle gain with rate γ g > 0. The jump operators are Hermitian for balanced gain and loss rates, γ l = γ g = γ. As anticipated in Sec. II, nontrivial topology of the LiouvillianL is reflected in quench dynamics of the driven-dissipative Kitaev chain only if the jump operators L i are Hermitian. The fact that the jump operators (4) act locally on a single lattice site simplifies the problem and enables, as described in Sec. V, an exact analytical solution of the master equation (2) in a certain limit. To demonstrate the validity of our findings beyond the case of purely local jump operators, we also consider quasilocal and Hermitian jump operators which are given by where κ is the rate of Markovian dissipation. In the main text, we focus on the driven-dissipative Kitaev chain with the purely local jump operators defined in Eq. (4). Our results for the quasilocal jump operators (5) are summarized in Appendix A.

B. Third quantization
Quadratic Liouvillians are equivalent to a class of non-Hermitian Hamiltonians, and, therefore, the topology of quadratic Liouvillians can be studied by using tools and concepts of non-Hermitian topological band theory [14]. The equivalence between quadratic Liouvillians and non-Hermitian Hamiltonians is established within the formalism of third quantization [64,65], which we summarize in the following.
The formalism of third quantization builds on the representation of a quadratic Hamiltonian H and linear jump operators L i in terms of 2N Majorana fermions, which are defined as These definitions lead to the following general expressions for the Hamiltonian and the jump operators: where A is an antisymmetric real 2N × 2N matrix, and l is a complex matrix with dimensions N × 2N . For Hermitian jump operators, the matrix l is real. More generally, any operator can be expanded in a basis of products of Majorana operators. Each element w s1 1 · · · w s 2N 2N of this basis is fully determined by the 2N occupation numbers s i = 0, 1, which appear as exponents. Therefore, we use the shorthand notation |w s ⟫ = |w s1 1 · · · w s 2N 2N ⟫, where s = (s 1 , . . . , s 2N ), and the double angular brackets emphasize the interpretation of |w n ⟫ as a superket, i.e., a vector in the Fock space of operators on which superoperators such asL act.
We denote fermionic annihilation and creation superoperators over this Fock space byĉ i andĉ † i , respectively, where the hat symbol distinguishes superoperators from ordinary operators. The fermionic superoperators are defined by their action on basis vectors, c i |w n ⟫ = δ ni,1 |w i w n ⟫ andĉ † i |w n ⟫ = δ ni,0 |w i w n ⟫. These definitions make use of the following property of Majorana operators: The product w i w si i , which is contained in w i w s , is given by w i w si i = w i for s i = 0, and by w i w si i = w 2 i = w 0 i = 1 for s i = 1. Further, these definitions imply the canonical anticommutation relations In terms of the fermionic superoperators, the Liouvillian in Eq. (2) can be written aŝ where denotes the transpose of a matrix, and1 is the identity superoperator. The real matrices X and Y are defined in terms of the Hamiltonian matrix A in Eq. (7) and the bath matrix M = l l * as where M R = Re(M ) and M I = Im(M ). By definition, the bath matrix M is Hermitian. Consequently, its real and imaginary parts, M R and M I , are symmetric and antisymmetric, respectively. In view of the representation (8) of the Liouvillion as a quadratic form in the fermionic superoperatorsĉ i and c † i , the master equation (2), rewritten as can be interpreted as a Schrödinger equation for superkets |ρ⟫, with an effective quadratic Hamiltonian iL. For an isolated system with antisymmetric X = −A and Y = 0, this effective Hamiltonian is Hermitian. In contrast, the dynamics of an open system is described by a non-Hermitian effective Hamiltonian iL. Topological properties of the non-Hermitian superoperator iL are encoded in its spectrum and eigenmodes. The Liouvillian can be diagonalized as whereb i andb i are fermionic superoperators which obey the anticommutation relations [14]. The fact that the spectrum of L is determined by the matrix X alone follows from the block-triangular form of Eq. (8). In contrast, the matrix Y affects the shape of the eigenmodes ofL as well as the steady state of the master equation (2). The steady state is the eigenstate of the Liouvillian with eigenvalue zero, L|ρ ss ⟫ = 0. According to Eq. (11), it is the right vaccum state of the modesb i . Now that we have set out the formal framework in which the equivalence between quadratic Liouvillians and non-Hermitian Hamiltonians is established, we proceed to analyze the non-Hermitian topological band structure of the matrix Z, which determines the complex spectrum of the Liouvillian.

IV. NON-HERMITIAN BAND THEORY FOR THE DRIVEN-DISSIPATIVE KITAEV CHAIN
Since the matrix Z is derived from a Liouvillian that generates the dynamics of the density matrix ρ according to Eq. (2), the spectrum of Z obeys certain conditions. First, to ensure Hermiticity of the time-evolved density matrix ρ(t) = eL t ρ 0 , for each eigenvalue λ i , also its anticomplex-conjugate −λ * i must be an eigenvalue of Z [14]. Second, the fact that the evolution of the density matrix approaches a steady state implies that Im(λ i ) ≤ 0, where the negative imaginary parts − Im(λ i ) determine the rate of relaxation to the steady state.
For a translationally invariant system, the spectrum of the Liouvillian can be determined analytically. In particular, for the driven-dissipative Kitaev chain, which is described by the Hamiltonian (1) and the local jump operators (4), the matrix Z can be decomposed into translationally invariant 2 × 2 blocks, Therefore, the matrix Z is block-diagonal in momentum space. Its spectrum consists of two bands, which are given by the eigenvalues of the Fourier transform [66] where k ∈ (−π, π] is the quasimomentum. Further, σ = (σ x , σ y , σ z ) is a vector of Pauli matrices, and The two bands of complex eigenvalues of z k are given by Before we proceed to analyze the band structure, we note that for J = ∆, the shifted matrix is unitarily equivalent to a non-Hermitian Bloch Hamiltonian studied in Refs. [1,59,60]. The addition of the constant matrix i (γ l + γ g ) 1 corresponds to a shift of the origin in the complex plane of the eigenvalues λ ±,k . While such a shift does not affect the topology of z k , it has severe consequences for the dynamics. In particular, as we discuss in Sec. V, dynamical critical behavior is induced when the spectrum of z k includes the eigenvalue λ = 0 in the complex λ-plane.
Since the eigenvalues of the Liouvillian are complex, gapped and gapless phases can be defined by addressing two distinct questions. First, we can ask whether the two bands λ ±,k are separated in the complex plane. This question leads to the definition of real and imaginary line gaps [8] as specified below. Further, for a complex spectrum, the point at which two bands touch in a gapless phase does not have to be λ = 0. Therefore, we can introduce another notion of a spectral gap by asking: Is λ = 0 contained in the spectrum of z k ? If this is not the case, the system has a point gap at λ = 0 [8]. As we detail below and in Sec. V, closings of the two types of gaps, i.e., real or imaginary line gaps that separate the bands in the complex plane, or the point gap that separates both bands from λ = 0, are associated with different types of topological transitions and leave distinct signatures, in particular, in the dynamics of the system. We first discuss real and imaginary line gaps. If the expression under the square root in Eq. (15) is positive for all values of k, the bands are separated along the real axis by a real line gap, which is indicated by the blue shading in Fig. 3(a). In the phase diagram in Fig. 2, a real line gap occurs in the phases R top and R tr , in which the rates of dissipation are bounded by If the rates γ l and γ g are increased beyond this bound, the system enters one of the two gapless phases, which are denoted by G edge and G tr . In these phases, the spectrum features exceptional points at momenta ±k * , at which the expression under the square root in Eq. (15) vanishes and thus the eigenvalues and the corresponding eigenvectors of z k coincide. This is exemplified in Fig. 3(b). Finally, at strong dissipation, when the expression under the square root in Eq. (15) is negative for all values of k, the bands are separated by an imaginary line gap as indicated by the orange shading in Fig. 3(c). The corresponding phases I edge and I tr are delineated by A non-Hermitian winding number W for a single band can be defined only in the gapped phases, where the bands are separated, and the definition depends on the symmetries of the model. The driven-dissipative Kitaev chain belongs to class BDI † in the nomenclature of Ref. [8], and obeys the following forms of time-reversal, particle-hole, and chiral symmetries: For γ l = γ g = 0, when z k = z † k , the class BDI † reduces to the class BDI of the isolated Hermitian Kitaev chain (1) with ∆ ∈ R. The definition of the winding number W for a single band of the isolated system can be generalized to the entire phases R top and R tr with a real line gap [8]. In Ref. [1], the calculation of W for the shifted matrixz k in Eq. (16) is carried out, and leads to the result W = 1 in the topological phase R top with |µ| < 2J, and W = 0 in the trivial phase R tr with |µ| > 2J. As we show in Sec. V, non-Hermitian topology as characterized by the winding number W = 1 for an individual band is reflected in the presence of entanglement spectrum crossings in quench dynamics. In contrast to the phases R top and R tr with a real line gap, the topology of individual bands of the Liouvillian is always trivial in the phases I edge and I tr with an imaginary line gap [8]. Consistently, we observe no entanglement spectrum crossings for quenches into I edge and I tr .
A topological classification of the gapless phases can be given in terms of a global invariant, which is not a property of individual bands, but rather of the entire Liouvillian. In particular, the global Berry phase Q defined and calculated for the shifted matrixz k in Eq. (16) in Ref. [59], takes the values Q = 1 for |µ| < 2J and Q = 0 for |µ| > 2J. We stress that the global Berry phase Q is defined in the entire phase diagram in Fig. 2, including the gapless phases G edge and G tr . In contrast, the classification in terms of the non-Hermitian winding number W applies only to the phases with real line gaps.
In the driven-dissipative Kitaev chain, a nonzero value of Q implies the existence of edge modes [60] even within the gapless phase G edge and the phase I edge with an imaginary line gap [41], and not only in R top . Edge modes are shown in the spectra in Fig. 3 as black diamonds. A topological transition occurs at the boundaries between G edge and G tr as well as I edge and I tr , in which edge modes disappear because the value of Q changes from Q = 1 to Q = 0. In crossing this transition, and when γ l = γ g , the second type of band gap defined above, i.e., the point gap at λ = 0, closes. In Fig. 3(b), the point gap is indicated with a gray disk. Closings of this gap induce critical power-law relaxation in the dynamics of both the entanglement spectrum and response functions, as we discuss in Sec. V below.

V. SYSTEMS WITH HERMITIAN JUMP OPERATORS
The definition of a non-Hermitian winding number W for complex bands λ ±,k , which we discussed in Sec. IV, is purely formal, and raises the question about its physical implications. In the following, we show that the topological transition in the spectrum of the Liouvillian, in which W changes from W = 1 in R top to W = 0 or undefined in R tr and G edge , respectively, is associated with a dynamical phase transition in the time evolution of the entanglement spectrum.
We begin by introducing the tools which are required to track the time evolution of the entanglement spectrum, and which have been used to establish the entanglement spectrum bulk-edge correspondence in quench dynamics of isolated Hermitian systems [33][34][35]. Then, as one of the key results of our work, we generalize the entanglement spectrum bulk-edge correspondence to the drivendissipative Kitaev chain.
In this section, we consider the effective Schrödinger equation (10) with a non-Hermitian matrix Z = −iX and Y = 0, which corresponds to Hermitian jump operators, i.e., γ l = γ g = γ in Eq. (4). We discuss the general case with γ l = γ g and thus Y = 0 in Sec. VI.
A distinguishing feature of systems with Hermitian jump operators is that they generically heat up to a fully mixed infinite-temperature state, ρ ss = ρ ∞ = 1/D, where D is the Hilbert-space dimension. This can be seen by noting that for L i = L † i , the dissipator in Eq. (3) takes the form of a double commutator, Therefore, the trivial state ρ ∞ = 1/D, which satisfies is indeed a steady state. Typically, this steady state is also unique. Degeneracy of the steady state requires the existence of at least two zero modes of the matrix Z [16,41]. For the driven-dissipative Kitaev chain with Hermitian jump operators, the spectra in Fig. 3 display only a single zero mode which is localized at the right edge of the chain. Thus, ρ ss = ρ ∞ is the unique steady state.

A. Topology from quench dynamics
To explore the dynamical topology of the drivendissipative Kitaev chain, we study the time evolution of a pure state |ψ 0 , the ground state of the isolated Kitaev chain, after sudden changes of parameters. The quench protocol is illustrated in Fig. 1: At time t = 0 + , the chemical potential is varied abruptly from its prequench value µ 0 to a postquench value µ 1 , and at the same time the system is connected to Markovian baths.
Due to the coupling to Markovian baths, described by the termDρ in Eq. (2), the pure initial state ρ 0 = |ψ 0 ψ 0 | evolves into a mixed state ρ(t) = eL t ρ 0 . Since H andL are quadratic, both the initial ground state and the time-evolved density matrix ρ(t) are Gaussian states, i.e., they can be written as the exponential of a quadratic form in the fermionic operators c i and c † i . Gaussian states are fully determined by the real and antisymmetric covariance matrix [67,68], The covariance matrix obeys the equation of motion Here, Γ(0) is the covariance matrix of the initial state.
As stated above, we assume that the system is prepared in the ground state of the isolated Kitaev chain, and the corresponding covariance matrix is elaborated in Appendix B. For Y = 0, the second term in Eq. (22) vanishes, and Γ(t) decays to zero in the steady state ρ ss = ρ ∞ = 1/D.
In the following, we assume that the number of lattice sites N is even, and we consider an equal bipartition of the system into its left and right halves. There are in total 2N sites of Majorana fermions. We denote the set of site indices of Majorana fermions which belong to the left half of the system by A = {1, . . . , N }. Then, the reduced covariance matrix of subsystem A is defined as Its eigenvalues, which come in pairs ±ξ i of bounded values 0 ≤ ξ i ≤ 1 with i = 1, . . . , N/2, form the single-particle entanglement spectrum [70][71][72].
For Gaussian states, the equality ξ i = tanh(ε i /2) relates the entanglement eigenvalues ξ i to the singleparticle eigenvalues ε i of the entanglement Hamiltonian H A , which is determined by the reduced density matrix [70][71][72]. The many-body entanglement spectrum, which is defined as the spectrum of ρ A , is thus given by [26] where e = e 1 , . . . , e N/2 , and where the numbers e i are to be understood as occupations of single-particle entan- glement levels that take the values e i = 0, 1. Therefore, for Gaussian states, the single-particle and manybody entanglement spectra contain the same information. In particular, as detailed below, Eq. (23) implies that a zero crossing in the single-particle entanglement spectrum leads to the simultaneous crossing of pairs of many-body entanglement eigenvalues.

time evolution of entanglement spectra
Before we continue, it is instructive to recall previous results for isolated systems. Notably, Gong and Ueda [33] connected the occurrence of entanglement spectrum crossings in the quench dynamics and the topology of the prequench and postquench Hamiltonians, H 0 and H, respectively, of a one-dimensional isolated system. Their key observation is that the time-evolved state |ψ(t) = e −iHt |ψ 0 , where |ψ 0 is the ground state of H 0 , is the ground state of a time-dependent parent Hamiltonian defined as H p (t) = e −iHt H 0 e iHt . If H is assumed to be band-flattened, the time dependence of H p (t) is periodic, and, therefore, in the first-quantized Hamiltonian H p (k, t), the time t can be interpreted as a second quasimomentum variable. By virtue of the entanglement spectrum bulk-edge correspondence, the evolution of the en-tanglement spectrum over one period gives the exact edge spectrum of the (1+1)-dimensional band-flattened Bloch Hamiltonian H p (k, t) [26]. In particular, nontrivial topology of H p (k, t) is reflected in the presence of entanglement edge modes which cross the band gap as a function of t. For the Altland-Zirnbauer classes in which nontrivial topology of H p (k, t) is realized, the two-dimensional topological index of the parent Hamiltonian H p (k, t) is the difference between the one-dimensional topological indices of H 0 and H [33]. Therefore, entanglement spectrum crossings occur, in particular, for quenches from a trivial phase of H 0 to a topological phase of H, and serve as a proxy for nontrivial topology of H. As we discuss in the following, the connection between entanglement spectrum crossings and nontrivial topology of the generator of dynamics remains intact if, instead of unitary time evolution generated by a Hamiltonian H, we consider the driven-dissipative dynamics of an open quantum many-body system as described by the master equation (2) with a non-Hermitian superoperatorL.
Throughout this section, we take the system to be initialized in the topologically trivial ground state of the isolated Kitaev chain with Hamiltonian (1), where J = ∆ sets the energy scale. Our results do not depend qualitatively on the precise prequench value of the chemical potential. For concreteness, we choose µ 0 = −3J. Figure 4(a) shows a quench to the non-Hermitian topological phase R top with W = 1 at µ 1 = −J and γ = 0.2J. The parameters of the postquench Liouvillian are indicated with a blue inverted triangle in Fig. 2, and the spectrum of the Liouvillian for these parameters is shown in Fig. 3(a). As can be seen in the latter figure, the bulk modes, which are shown as blue and orange diamond symbols, decay with the same rate − Im(λ bulk ) = 2γ = 0.4J. Then, for the covariance matrix, and, consequently, its eigenvalues which form the single-particle entanglement spectrum, Eq. (22) implies a decay rate of 4γ = 0.8J. To best illustrate the evolution of the entanglement spectrum in Fig. 4, we account for this overall decay by rescaling the entanglement eigenvalues as ξ = ξe 4γt . With this rescalling, in Fig. 4(a), most of the single-particle entanglement eigenvalues, which are shown as blue lines, stay close to ±1. However, a pair of eigenvalues ±ξ i undergoes repeated zero crossings. To trace these crossings, we employ the Pfaffian of the reduced covariance matrix, The Pfaffian, which can be defined for any antisymmetric matrix, has the following properties: Its absolute value is given by the square root of the determinant, and, most importantly, it changes sign whenever a pair of entanglement eigenvalues ±ξ i crosses zero. In all panels of Fig. 4, the rescaled Pfaffian Pf = Pf e 4γN t is shown as a black line, and its sign is indicated with blue shaded areas. At long times, the covariance matrix (22), and, consequently, also the Pfaffian (24) decay to zero. The entanglement spectrum crossings disappear when the postquench parameters are within a trivial phase in which W = 0 or W is undefined. Figure 4(b) shows the entanglement spectrum dynamics following a quench to the real-line gapped trivial phase R tr with W = 0 at µ 1 = −2.5J and γ = 0.1J, as indicated by the red circle in Fig. 2. As above, the decay rates of all bulk modes are equal to 2γ. Therefore, rescaling the entanglement eigenvalues with e 4γt removes the overall decay of the entanglement spectrum and leaves all entanglement eigenvalues close to ±1. The evolution of the entanglement spectrum shows richer phenomenology for quenches to the gapless phases G edge and G tr , in which the winding number W is undefined. Figure 4(c) shows a quench to the phase G edge at µ 1 = −J and γ = 0.6J, as indicated by the green triangle in Fig. 2. The corresponding spectrum of the postquench Liouvillian is shown in Fig. 3(b). Decay rates of bulk modes are within a finite range of values, and, therefore, the rescaling of the entanglement spectrum with e 4γt leads to the occurrence of exponentially growing values. Thus, the rescaled Pfaffian in Fig. 4(c) soon exceeds one. Crucially, and in agreement with the fact that in a gapless phase the postquench Liouvillian is not characterized by a nontrivial non-Hermitian winding number, the Pfaffian does not change sign. Figure 4(d) presents the time evolution of the entanglement spectrum after a quench to the imaginary-line gapped phase I edge at µ 1 = 0 and γ = 1.1J, as indicated by the orange square in Fig. 2. For these postquench parameters, the two bands (15) are purely imaginary as shown in Fig. 3(c). This is reflected vividly in the purely decaying behaviour of the entenglement spectrum.
We proceed with a quantitative analysis of the time evolution of the entanglement spectrum for quenches to the topological phase R top . The color scale within the phase R top in Fig. 2 encodes the time t 1 at which the first zero crossing in the single-particle entanglement spectrum occurs. In Fig. 4(a), the time t 1 is indicated with a vertical black line. The values of t 1 diverge at the boundaries of the phase R top . This is illustrated in the inset of Fig. 2, where t 1 is shown for variations of the parameters of the postquench Liouvillian along the black arrows in the main panel of the figure, i.e., for δµ = µ 1 + 2J → 0 and γ = 0, and for µ 1 = −J and δγ = 0.5J − γ → 0, respectively. The numerical data is consistent with a square-root singularity, which is shown as a gray line in the inset of Fig. 2, and which corresponds to the behavior t 1 ∼ δµ − and t 1 ∼ δγ − , where = 1/2 can be interpreted as a dynamical critical exponent for entanglement spectrum crossings. In Sec. V A 2, we show that the value = 1/2 is exact for the specific case of a postquench Liouvillian with flat spectrum, in which an exact solution of the quench dynamics is possible. In the numerical data shown in the inset of Fig. 2, the singularity of t 1 is cut due to the finite system size. The values of t 1 saturate earlier if the phase boundary is approached by varying the strength of dissipation.
The singularity of t 1 at the boundaries of the topologi- cal phase R top raises the question, whether the time evolution of the entanglement spectrum exhibits also other signatures of dynamical criticality. To complement our analysis of oscillatory behavior of entanglement eigenvalues, which is characterized by the time scale t 1 , we illustrate in Fig. 5 the decay of the largest eigenvalue ξ max . It exhibits critical power-law relaxation at the phase boundary between the gapless phases G edge and G tr , where the global Berry phase jumps from Q = 1 to Q = 0 and, concomitantly, the edge modes of the Liouvillian disappear. This form of decay is also realized generically for quenches in isolated Hermitian systems [73]. However, the decay of ξ max follows a simple exponential form across the boundary between R top and G edge . In particular, in Fig. 5(a) we show the decay of ξ max for three different values of the rate of dissipation γ below, at, and above the critical value γ c = 0.5J for µ 1 = −J. For each value of γ, the approximately linear behavior of ξ max on the semilogarithmic scale indicates dominantly exponential relaxation. Indeed, we find good agreement of the numerical data with ξ max ∼ e −Γmaxt and ξ max ∼ e −Γmaxt / √ t for γ < γ c and γ > γ c , respectively. The numerically determined decay rate Γ max is shown as a function of γ in Fig. 5(b). It agrees well with the analytical expression for γ < γ c , which is inspired by the exact result for the decay rate of the retarded response function as elaborated in Sec. V B below. Interestingly, Γ max takes a maximum value of Γ max = 2J exactly at the critical point γ = γ c = 0.5J. This observation runs counter to the common expectation that relaxation to the steady state becomes exceedingly slow in the vicinity of a phase transition. In Sec. V B, we explain this unusual behavior through a mechanism of many-body critical damping.
The dynamics of the entanglement spectrum exhibits rather different behavior across the boundary between G edge and G tr . As discussed above, there is no qualitative change in the oscillatory dynamics: For quenches to both G edge and G tr , there are no zero crossings of entanglement eigenvalues. However, as illustrated in Fig. 5(c), the decay of the largest entanglement eigenvalue ξ max shows critical slowing down as the phase boundary is approached, and power-law relaxation ξ max ∼ (Jt) −1/2 exactly on the phase boundary. The origin of this behavior is the closing of the point gap at λ = 0. We provide a quantitative analysis of dynamical criticality in terms of two-time correlation functions in Sec. V B.

Driven-dissipative fermion parity pump
We proceed to complement our numerical results for the time evolution of entanglement spectra with an exact analytical solution, which is enabled by a judicious choice of the initial state and the parameters of the postquench Liouvillian. First, this provides us with an interpretation of topological quench dynamics as a fermion parity pump; and second, for specific parameter values, we are able to show explicitly that the entanglement spectrum dynamical critical exponent takes the exact value = 1/2.
To facilitate analytical progress, we assume that the system is initialized in the vacuum of Dirac fermions c i , which is the ground state of the Hamiltonian of the Kitaev chain for J = ∆ > 0 and µ 0 → −∞. Further, we choose the parameters of the postquench Liouvillian as J = ∆ > 0, µ 1 = 0, and γ l = γ g = γ. For this choice of postquench parameters, the bulk spectrum (15) for an infinite chain becomes flat with The spectrum exhibits a real line gap when γ < γ c = J, and the expression under the square root is positive. In contrast, for γ > J, the two flat bands are separated by an imaginary line gap. For a Liouvillian with flat bulk bands, also the spectrum for an open finite chain of length N can be obtained analytically, and is given by {λ ± , λ L , λ R }. In a finite system, the values λ ± given in Eq. (26) have degeneracy N − 1 each, and λ L = −i4γ and λ R = 0 are eigenvalues corresponding to edge modes which are localized on the left and the right end of the chain, respectively. For γ < J, the real part of the eigenvalues λ ± sets a time scale, for oscillations of the bulk modes. As we show in the following, this is the time scale on which the fermion parity of many-body entanglement eigenstates is reversed and, concomitantly, both single-particle and many-body entanglement eigenvalues cross. More precisely, as illustrated in Fig. 6 where the times mT , with m ∈ N 0 a nonnegative integer, are indicated with green vertical lines, entanglement spectrum crossings occur at times t m which obey mT < t m < (m + 1) T .
The key result, based on which we establish the connection between entanglement spectrum crossings and the reversal of the fermion parity at the edges of the Kitaev chain, is an exact expression for the reduced density matrix ρ A (mT ) at times t = mT with m ∈ N 0 , which we derive in Appendix C. As we show there, for the choice of initial state and postquench Liouvillian specified above, the reduced density matrix ρ A (mT ) is diagonal in a basis of Fock states with fixed number of particles on each lattice site, which we denote by |n A where n = n 1 , . . . , n N/2 is a vector of lattice-site occupation numbers with n i = 0, 1. In terms of projectors P n = |n A n| on these basis states, the reduced density matrix can be written as where the sum is over all combinations e = e 1 , . . . , e N/2 of single-particle entanglement occupation numbers e i = 0, 1, and the coefficients Ξ e,m are the eigenvalues of ρ A (mT ) and determine the many-body entanglement spectrum. We specify the relation between the latticesite occupation numbers n e,m and the entanglement occupation numbers e below.
For the exact solution (28) derived in Appendix C, the many-body entanglement eigenvalues Ξ e,m take the general form given in Eq. (23), with only two distinct single-particle entanglement eigenvalues, The first eigenvalue ξ 1,m is nondegenerate, whereas ξ 2,m has degeneracy N/2 − 1. Both single-particle entanglement eigenvalues decay exponentially with rates 6γ and Figure 6. (a) Time evolution of the single-particle entanglement spectrum for a quench from µ0 → −∞, to µ1 = 0, and γ l = γg = γ = 0.1J for N = 10. The solid black line is the Pfaffian, and the blue shading indicates the sign of the Pfaffian. Green vertical lines mark integer multiples of the time scale T in Eq. (27). Black and red crosses indicate, respectively, ±ξ1,m and ±ξ2,m, for the analytical results given in Eq. (29). (b) Fermion parity pump illustrated with the evolution of the many-body entanglement spectrum for a quench with the same parameters as in (a). Black and red lines correspond to many-body entanglement eigenstates with even and odd parity, respectively. The parity of the entanglement ground state is even (odd) in the regions with (without) gray shading. Crosses indicate analytical results at multiples of T marked by green vertical lines.
4γ, respectively. The exact analytical results for ±ξ 1,m and ±ξ 2,m are shown, respectively, as black and red crosses in Fig. 6(a). We find perfect agreement with the numerically calculated entanglement spectrum {±ξ i (t)}, which is indicated with blue lines, at times t = mT .
The reversal of the fermion parity at the ends of the chain is encoded in the relation between entanglement occupation numbers e and the site occupation numbers n e,m in Eq. (28). As we show in Appendix C, n e,m = e for m even, e for m odd, where vectors with an overbar are defined byē 1 = (e 1 + 1) mod 2, whereasē i = e i for i = 2, . . . , N/2. That is, e andē differ only in the first component.
To demonstrate the reversal of the fermion parity after each time step of duration T , we follow the evolution of the entanglement ground state, i.e., the state which corresponds to the smallest eigenvalue of the entanglement Hamiltonian H A . According to the relation between the reduced density matrix ρ A and the entanglement Hamiltonian, ρ A = e −H A /Z A where Z A = tr A (e −H A ), the entanglement ground state is the eigenstate of ρ A with the largest many-body entanglement eigenvalue Ξ e in Eq. (23). The latter equation implies that the largest many-body entanglement eigenvalue is obtained if all entanglement occupation numbers are equal to zero, e = 0. The corresponding lattice-site occupation numbers (30) are n 0,m = 0 for even values of m, and n 0,m =0 = (1, 0, . . . , 0) for odd values of m. The entanglement ground state is thus given by |0 A when m is even, and by |0 A when m is odd. These states have opposite fermion parity: as follows directly from the definition of the fermion parity operator for subsystem A, i.e., the left half of the chain, We proceed to connect the reversal of the fermion parity of the entanglement ground state to crossings in the many-body entanglement spectrum. Indeed, as we show in Appendix D, the entanglement eigenstates have definite parity at all times t. Therefore, a change of the fermion parity of the entanglement ground state between t = mT and t = (m + 1) T implies that there is a crossing in the entanglement spectrum in which the entanglement ground state and the first excited state, which have opposite fermion parity, switch places. Numerically, as shown in Fig. 6(b), we find that there is exactly one crossing between the entanglement ground state and the first excited state, which occurs at a time t m in the interval mT < t m < (m + 1) T .
A crossing between the entanglement ground state and the first excited state in turn implies a zero crossing in the single-particle entanglement spectrum. To see this, we consider the difference between the entanglement eigenvalues of the ground state and first excited state. As we have already noted, the ground state corresponds to e = 0. The first excited state is obtained by occupying the single-particle entanglement level with the smallest eigenvalues ξ 1 (t), where we order the single-particle entanglement eigenvalues ξ i (t) such that 0 ≤ ξ 1 (t) ≤ · · · ≤ ξ N/2 (t). Thus, for the first excited state, the entanglement occupation numbers are e =0. According to Eq. (23), which expresses the many-body entanglement spectrum in terms of the single-particle entanglement spectrum, the difference between the entanglement eigenvalues Ξ 0 (t) and Ξ0(t) is given by This difference is equal to zero when ξ 1 (t) = 0. That is, when the fermion parity for the many-body entanglement ground state is reversed, there is a zero crossing in the single-particle entanglement spectrum. Interestingly, a zero crossing of ξ 1 (t) leads to simultaneous crossings of all pairs of many-body entanglement eigenstates with entanglement occupation numbers e and e. The difference between the corresponding entanglement eigenvalues is given by and vanishes again when ξ 1 (t) = 0. Figures 6(a) and (b) illustrate this connection between zero crossings in the single-particle entanglement spectrum and simultaneous crossings of pairs of many-body entanglement eigenvalues. As can be seen in Fig. 6(b), there are further occasional crossings in the many-body entanglement spectrum that are not related to zero crossings in the singleparticle entanglement spectrum, and which are not of topological origin.
Finally, we demonstrate numerically that the simultaneous crossings in the many-body entanglement spectrum occur between pairs of states which have opposite parity, and that, therefore, a zero crossing in the singleparticle entanglement spectrum leads to the simultaneous reversal of the fermion parity in all entanglement eigenstates. To this end, we first note that according to Eq. (30), pairs of entanglement eigenstates with entanglement spectrum occupation numbers e andē have opposite parity at times t = mT . As discussed in Appendix D, parity is conserved during the dynamics. We can thus track the parity of all entanglement eigenstates at all times, even without calculating the full many-body density matrix ρ A . In Fig. 6(b), black and red lines correspond to even and odd parity, respectively. The parity of the entanglement ground state is even in regions with gray shading. At the boundaries of these regions, which agree with the boundaries of the blue-shaded regions in Fig. 6(a) that mark changes of the sign of the Pfaffian due to zero crossings in the single-particle entanglement spectrum, simultaneous crossings in the many-body entanglement spectrum occur between pairs of states with opposite parity.
The above considerations establish the interpretation of entanglement spectrum crossings as a fermion parity pump for a particular choice of parameters. In general, this interpretation remains valid for continuous variations of pre-and postquench parameters within the trivial phase of the isolated Kitaev chain and the non-Hermitian topological phase of the driven and open Kitaev chain, respectively. To convince oneself that this is the case, it is sufficient to note that, as detailed in Appendix D, fermion parity is always a good quantum number of entanglement eigenstates for the models and dynamics we consider here, i.e., quench dynamics starting from a state with definite parity, and where time evolution is generated by a quadratic Liouvillian.
The period T of the fermion parity pump diverges when the rate of dissipation approaches the critical value of γ c = J. Specifically, by inserting δγ = γ c − γ in Eq. (27), we obtain Consequently, also the times t m at which entanglement spectrum crossings occur exhibit a square-root singularity, t m ∼ δγ − with = 1/2. In the inset of Fig. 2, we present numerical evidence for the universality of the value of the entanglement spectrum critical exponent .
In particular, our numerical results indicate that = 1/2 when we (i) consider dispersive bands of the postquench Liouvillian, (ii) approach the phase boundary of R top from different directions, by tuning either γ or µ 1 , (iii) consider a more general topologically trivial initial state with J, ∆ = 0 and 2J < |µ 0 | < ∞, i.e., a trivial state that is different from the vacuum of Dirac fermions, and (iv) study different models such as the one with quasilocal jump operators defined in Eq. (5) and discussed in detail in Appendix A.

B. Dynamical critical behavior
In Sec. V A 1, we have encountered two distinct types of dynamical critical behavior of the entanglement spectrum: the time t 1 of the first entanglement spectrum crossing diverges with an exponent = 1/2 at the boundary of R top , while the relaxation of the largest entanglement eigenvalue ξ max follows a simple exponential form across the phase boundary; In contrast, there are no entanglement spectrum crossings on both sides of the phase boundary between G edge and G tr . However, ξ max exhibits critical power-law relaxation exactly at this phase boundary. In the following, we show that a unified perspective on these different forms of dynamical criticality can be obtained by studying two-time correlation functions.
We consider the Keldysh Green's function, together with the retarded and advanced response functions [74], in terms of which all two-time correlation functions can be specified: where θ(t) is the Heaviside step function. We omit the usual factors ∓i in the definitions of the retarded and advanced response functions to obtain real-valued quantities. The expectation values are taken in the steady state ρ ss and, therefore, it is sufficient to include a single time argument since two-time averages depend only on the time difference. For dynamics described by a master equation in Lindblad form, two-time averages can be calculated with the aid of the quantum regression theorem [75] as detailed in Appendix E. Further, we assume here an infinite bulk system such that the correlation functions are translationally invariant. The two-time correlation functions defined above are distinguished from the covariance matrix Eq. (21) in two aspects: first, the expectation value is taken in the steady state ρ ss in the two-time correlation functions, and not in the time-dependent state ρ(t) = e Lt ρ 0 as in the case of the covariance matrix; second, here the Majorana operators w i and w j are evaluated at different times. In the limit t → ∞, the covariance matrix approaches the equal-time Keldysh Green's function, Γ i,j (t) → −χ K i,j (0). We consider here a system with a trivial steady state ρ ss = 1/D, as is typically the case for Hermitian jump operators. Then, the trace of the commutator in Eq. (36) evaluates to zero, and the Keldysh Green's function vanishes identically. The retarded and the advanced response functions carry the same qualitative features in terms of divergent timescales and scaling behavior. For concreteness, in the following, we focus on the equalposition retarded response function χ R i,i (t). We present results for odd values of i. Even values of i lead to the same qualitative behavior.
As in our discussion of the interpretation of topological quench dynamics in terms of a fermion parity pump in Sec. V A 2, it is worthwhile to consider first the case of a flat-band Liouvillian, which enables exact analytics, and brings out the relation between the dynamics of the retarded response function and of a classical damped harmonic oscillator most clearly.
The bands λ ±,k in Eq. (15) are flat when the chemical potential is equal to zero. Then, the evolution equation for the retarded response function, as detailed in Appendix E, can be written in the form of the equation of motion of a damped harmonic oscillator which is kicked at t = 0: where δ(t) is the Dirac delta function. The solution which obeys the boundary condition χ R i,i (t) = 0 for t < 0 reads where ω = 2 J 2 − γ 2 . In the phase R top where γ < J, the response function exhibits underdamped oscillations. The boundary between the phases R top and I edge is at the point of critical damping where γ = J. As γ is increased beyond this value, ω becomes imaginary, and the response function ceases to oscillate.
While the dynamics of the response function for µ = 0 is completely analogous to that of a single damped harmonic oscillator, the many-body character of the response function is restored when the chemical potential takes nonzero values and, consequently, the bands λ ±,k in Eq. (15) are dispersive. For this case, the response function cannot be calculated exactly. However, its asymptotic late-time behavior, which is of main interest with regard to dynamical criticality, can still be obtained analytically as described in Appendix E.
The behavior of the retarded response function in the phase R top is illustrated in Fig. 7(a). In both phases with a real line gap, R top and R tr , the late-time asymptotic form of the response function reads That is, the response function is a superposition of two oscillatory contributions with frequencies ω ± , and damps out exponentially with a decay rate Γ R = 2γ. It is worth stressing that the response function exhibits oscillations and zero crossings in both the topological and the trivial phase. The amplitudes, phases, and frequencies of the oscillatory components with σ = ± are given by When γ is increased towards the critical value γ c = J + µ/2 which marks the boundary between the phases R top and G edge for µ < 0, the amplitude A + diverges and, therefore, the asymptotic behavior of the response function is dominated by the contribution which oscillates at the frequency ω + . This frequency vanishes at the phase boundary. In particular, for γ = γ c − δγ, we obtain ω + ∼ 2 (2J + µ) δγ. Thus, upon approaching δγ = 0, the oscillation period 2π/ω + ∼ δγ − diverges with the same exponent = 1/2 as t 1 . For µ > 0, the roles of σ = + and σ = − are interchanged: A − diverges and ω − vanishes at the phase boundary, while A + and ω + remain finite. For both a flat and a dispersive spectrum of the Liouvillian, which lead to Eqs. (40) and (41) for the response function, respectively, the period of oscillations of the response function diverges when the real line gap closes. However, there are notable differences between Eqs. (40) and (41). First, the latter equation does not reduce to Eq. (40) in the limit µ → 0, which shows that this limit does not commute with the late-time limit t → ∞ in which Eq. (41) is valid. The continuum of modes, which contribute to the response function in a manybody system, results in additional structure of Eq. (41) in comparison to Eq. (40). In particular, the exponential decay is augmented by an algebraic factor 1/ √ t in Eq. (41). Further, while the entire spectrum becomes purely damped with Re(λ) = 0 for critical damping with γ = J at µ = 0, for a generic point on the boundary between a real-line gapped phase and a gapless phase, the behavior of the many-body response function is determined by a continuum of modes λ. Even within the gapless phases there are modes for which Re(λ) = 0, as shown in Fig. 3(b). These modes give subleading oscillatory contributions with frequency Re(λ) to the response function, while the leading asymptotic late-time behavior is determined by modes with the smallest values of the damping rate − Im(λ). For these modes, the oscillation frequency vanishes, Re(λ) = 0. The damped oscillatory contribution is clearly visible for the line in Fig. 7(a) which corresponds to G edge . Due to these differences between Eqs. (40) and (41), we refer to the divergence of the period of oscillations of the response function at finite values of the chemical potential as many-body critical damping. Exactly on the boundary between R top and G edge , the retarded response function behaves as where γ c = J − |µ| /2 and Γ R = 2γ c . The late-time form of the response function in the gapless phase G edge is shown in Fig. 7(a). In fact, the leading asymptotic behavior of the response function is the same in two gapless phases G edge and G tr as well as the imaginary-line gapped phases I edge and I tr . This is because the modes with the smallest decay rates − Im(λ), which determine the late-time asymptotics, have the same behavior Re(λ) = 0 and − Im(λ) ∼ a + bk 2 for k → 0 with positive real coefficients a and b, as illustrated in Figs. 3(b) and 3(c) for G edge and I edge , respectively. As a result, we find Here, the decay rate is given by As can be seen in Fig. 7(a), the behavior of the response function in the gapless phase G edge and the imaginary-line gapped phase I edge differs by the presence of subleading oscillatory contributions in the former case. The above results for the late-time asymptotics of the retarded response function allow us to justify our finding in Sec. V A 1, that the numerically determined decay rate Γ max of the largest entanglement eigenvalue ξ max agrees with the analytical expression given in Eq. (25). Indeed, by combining Eqs. (41), (43), and (44), we see that decay rate Γ R of the retarded response function is related to Γ max by The factor of two can be explained by comparing Eqs. (22) and (E5) for the covariance matrix and the retarded response function, respectively. In the former equation, the matrix X appears twice in the exponent, and, therefore, the covariance matrix decays twice as fast.
The maximum of Γ R and Γ max at γ = γ c is again in agreement with the physics of a damped harmonic oscillator, which decays fastest for critical damping. Upon approaching the critical lines at µ = ±2J which separate the gapless phases G edge and G tr , the decay rate (45) vanishes as Γ R ∼ δµ 2 /(4γ), where we set µ = µ c − δµ for µ c = 2J and δµ > 0. That is, the response function exhibits critical slowing down with a divergent relaxation time scale. Exactly on the critical lines at µ = ±2J, the response function decays as a power law, As is illustrated in Fig. 7(b), subleading oscillations of the response function damp out for t → ∞.
In stark contrast to the critical relaxation of the response function at the boundary between the gapless phases G edge and G tr , according to Eq. (43), the response function exhibits exponential decay on the line of manybody critical damping, i.e., the phase boundary between R top and G edge . This difference of exponential vs. powerlaw decay can be traced backed to the different types of gap closings at the boundaries between R top and G edge , and G edge and G tr , respectively. In the latter case, the point gap at λ = 0 closes, and critical relaxation results from the existence of a continuum of modes which includes λ = 0: Technically, scaling behavior is due to the integral in Eq. (E6) over a continuum of modes, which develops a singularity when the continuum includes zero. In contrast, although the real line gap closes at the phase boundary between R top and G edge , the point gap at λ = 0 remains intact.
Finally, we note that the retarded response function can be calculated exactly at µ = ±2J and in the limit γ = 0 + when we assume that the system reaches the trivial infinite-temperature steady state even for vanishingly small dissipation. We find where J n (z) is a Bessel function of the first kind with late-time asymptotic behavior Thus, for γ = 0 + , oscillations of the response function persist indefinitely, as can be seen in Fig. 7(b).

VI. SYSTEMS WITH NON-HERMITIAN JUMP OPERATORS
So far, we have mainly focused on a special case of the driven-dissipative Kitaev chain: For γ l = γ g , the jump operators (4) are Hermitian, and the steady state of the time evolution described by the master equation (2) is a trivial infinite-temperature state. We showed that for Hermitian jump operators, non-Hermitian topology of the LiouvillianL is revealed through entanglement spectrum crossings in quench dynamics. Now we proceed to discuss the general case of non-Hermitian jump operators with γ l = γ g , where the interplay of unitary and dissipative dynamics results in a nontrivial steady state. As we explain in the following, in this general case, the presence or absence of entanglement spectrum crossings is determined not only by the topology ofL, but also by properties of the initial and steady states.
In particular, as discussed in Sec. V A 2, each zero crossing in the single-particle entanglement spectrum is accompanied by a reversal of the fermion parity of the entanglement ground state. Therefore, if the initial and steady states have opposite parities, an entanglement spectrum crossing must occur even if the topology ofL is trivial. However, as we show below, the connection between non-Hermitian topology ofL and entanglement spectrum crossings can be restored by means of Hermitianization of the jump operators [14].

A. Fermion parity and Pfaffian in the initial and steady states
For Hermitian jump operators, which lead to a trivial steady state ρ ss = ρ ∞ = 1/D, also the reduced density matrix is trivial, ρ A,ss = tr A c (ρ ss ) = 1 A /D A with D A = 2 N/2 . The corresponding many-body entanglement spectrum is flat, and states with even and odd fermion parity are degenerate. Consequently, the parity of the entanglement ground state is undetermined. Further, the covariance matrix of the steady state vanishes identically. Thus, also the sign of the steady-state Pfaffian Pf ss is undetermined. In contrast, for non-Hermitian jump operators with γ l = γ g , the entanglement spectrum of the steady state is nondegenerate, and the corresponding entanglement ground state has definite fermion parity p A,ss . The sign of the steady-state Pfaffian is given by sgn(Pf ss ) = p A,ss = (sgn(δ)) N/2 (50) with δ = γ l −γ g . As we show in Appendix F, this relation holds for J = ∆ > 0, and arbitrary postquench values of the chemical potential µ 1 . Notably, if δ < 0, Eq. (50) exhibits an even-odd half-system-size effect for sgn(Pf ss ). Here, as in Sec. V, we restrict our discussion to systems with an even number of lattice sites, N ∈ 2N, such that N/2 ∈ N is an integer.
Equation (50) should be compared with the corresponding relation for the initial state of the quench protocol, i.e., the ground state of the Hamiltonian (1) of the Kitaev chain. In the trivial phase with |µ 0 | > 2J, we find For values of µ 0 , δ, and N/2 that lead to disagreements between the right-hand sides of Eqs. (50) and (51), there must be a reversal of the sign of the Pfaffian in the course of the time evolution, and, consequently, also a crossing in the entanglement spectrum. This conclusion does not rely on the topology of the Liouvillian, and can even be in conflict with expectations based on the latter. In particular, for δ < 0 and µ 0 < 0, Eqs. (50) and (51) imply that an entanglement spectrum crossing must occur if N/2 is odd. By contrast, we should expect no crossing if N/2 is even. This even-odd half-system-size effect cannot be captured by the non-Hermitian winding number W , which is defined for the bands of an infinite bulk system. We proceed to demonstrate the occurrence of the evenodd effect explicitly in the time evolution of entanglement spectra, and we then show how the discrepancy between the entanglement spectrum crossings and the topology of the Liouvillian can be resolved through Hermitianization of the jump operators.

B. Quench dynamics with non-Hermitian and Hermitianized jump operators
The impact of the even-odd effect of the Pfaffians of the initial and steady states on the time evolution of the entanglement spectrum is illustrated in Fig. 8. In particular, Fig. 8(a) shows the evolution of the single-particle entanglement spectrum for a quench to the topological phase R top with γ l = 0.3J, γ g = 0.7J, and µ 1 = 0. For these parameter values, the bulk spectrum (15) of the Liouvillian becomes flat, and we might expect periodically recurring entanglement spectrum crossings. However, as shown in the figure, there is no entanglement spectrum zero crossing for N/2 = 50, and only a single crossing for N/2 = 51. These findings are in agreement with our expectations based on Eqs. (50) and (51).
In the course of the time evolution, the Pfaffian drops rapidly to very small values. Indeed, as discussed in Appendix F 2, the steady-state Pfaffian Pf ss is exponentially small in N . However, there we also show that the parity p A,ss , and, therefore, the sign of Pf ss are well-defined for arbitrarily large systems.
Quenches to the gapless trivial phase G edge with the same values of γ l and γ g as above and µ 1 = −1.5J are shown in Fig. 8(b). Again, there is no entanglement spectrum zero crossing for N/2 = 50, and only a single zero crossing for N/2 = 51.
The presence or absence of entanglement spectrum crossings can be reconciled with the topological properties of the LiouvillianL through Hermitianization of the jump operators, i.e., by continuously deforming the jump operators to make them Hermitian as described in Ref. [14]. Hermitianization of the jump operators corresponds to a transformation of the matrices X and Y that determine the Liouvillian through Eq. (8), in the course of which X = const., while Y → 0. That is to say, the spectral gap and the symmetries of the Liouvillian are preserved, whereas the steady state becomes trivial with Pf ss = 0. Consequently, for a master equation with Hermitianized jump operators, the even-odd half-system size effect is resolved.
In particular, for the quenches shown in Figs. 8(a) and 8(b), the time evolution of the entanglement spectrum with Hermitianized jump operators is shown in, respectively, the left and right panel of Fig. 8(c). Both for N/2 = 50 as shown in the figure and for N/2 = 51, en- The connection between entanglement spectrum crossings and topology is restored for the deformed Liouvillian with Hermitianized jump operators and a trivial steady state. The left and right panels correspond to quenches with the parameters of (a) and (b), respectively. We show the rescaled entanglement spectrum ξ = ξe 2(γ l +γg)t and Pfaffian Pf = Pf e 2(γ l +γg)N t for N/2 = 50. The corresponding plots for N/2 = 51 look essentially the same. tanglement spectrum crossings occur only when the non-Hermitian winding number of the postquench Liouvillian takes the value W = 1.
Hermitianization of the jump operators allows us to explore the entire topological phase diagram for γ l = γ g in terms of entanglement spectrum zero crossings. This is illustrated in Fig. 9. For the value µ 1 = 0 shown in the figure, the real-line gapped phase R top borders the phase I edge with an imaginary line gap. The time scale t 1 , at which the first entanglement spectrum crossing occurs, diverges at the boundary of the topological phase R top . In particular, as shown in the inset of Fig. 9, t 1 exhibits power-law behavior, t 1 ∼ δγ − . For the value of the exponent, we find = 1/2 in agreement with our results for Hermitian jump operators.

VII. CONCLUSIONS AND OUTLOOK
Our work establishes the entanglement spectrum as a tool to study non-Hermitian dynamical topology of driven and open quantum many-body systems. We circumvent the problem that complex spectra do not allow to define ground states [76] by considering the time evolution of entanglement spectra after a quench.
We show that for systems with Hermitian or Hermitianized jump operators, the presence of entanglement spectrum crossings reflects nontrivial topology of individual complex bands of the Liouvillian. The topology of individual bands is characterized by the non-Hermitian winding number W . In contrast, as discussed in Ref. [60], the presence of edge modes is tied to the global Berry phase Q, which pertains to the full spectrum, and is well-defined also in gapless phases. The existence of two distinct topological transitions, at which, respectively, the value of the non-Hermitian winding number W and the global Berry phase Q changes, is unique to non-Hermitian systems. Indeed, as shown in the dynamical topological phase diagram in Fig. 2, the correspond-ing phase boundaries between R top and G edge as well as G edge and G tr , coincide in the Hermitian limit where the rate of dissipation γ vanishes.
The characterization of non-Hermitian dynamical topology through the entanglement spectrum opens up interesting directions for future research. While we focus here on a driven and dissipative version of the Kitaev chain, which belongs to the non-Hermitian symmetry class BDI † [8], our approach can be applied directly to Liouvillian generalizations [14] of the Altland-Zirnbauer symmetry classes that accommodate nontrivial (1 + 1)dimensional dynamical topology [33].
Apart from different one-dimensional and noninteracting models, as well as other symmetry classes, a promising prospect is to generalize the entanglement spectrum bulk-edge correspondence to quench dynamics of drivendissipative systems in higher spatial dimensions, and of interacting open quantum many-body systems. Analogous generalizations were obtained for isolated Hermitian systems in Refs. [33,34]. For noninteracting Hermitian and non-Hermitian systems, the existence of a band structure is sufficient to define topological invariants, even without referring to a particular state. By contrast, for interacting Hermitian systems, a well-defined ground state commonly forms the basis of any discussion of topology [23]. This raises the following questions: How can the topology of the complex spectra of quartic Liouvillians, which describe interacting open quantum many-body systems, be classified? Do sharp dynamical topological transitions exist also for such systems? Our approach based on the entanglement spectrum provides a way to address, in particular, the latter question. For example, we expect that also for a generalization of an interacting Kitaev chain [26], which includes Markovian drive and dissipation, the time evolution of the manybody entanglement spectrum carries signatures of non-Hermitian dynamical topology of the quartic Liouvillian.
A substantial part of our work is concerned with a detailed investigation of dynamical criticality at the two distinct topological phase transitions which are related to the non-Hermitian winding number W and the global Berry phase Q. We find strikingly different behavior: While the transition that is described by Q exhibits conventional critical relaxation, we observe the phenomenon of many-body critical damping at the boundaries of the topological phase R top with W = 1.
Our analysis of dynamical criticality, and, in particular, of many-body critical damping, opens up new avenues for future research, which bear relevance beyond the field of non-Hermitian topology. An obvious question is whether it is possible to identify different universality classes of many-body critical damping. To address this question, it will be interesting to study, e.g., a driven and open Kitaev chain with long-range hopping and pairing [77][78][79], or to explore whether the dynamical critical exponent for entanglement spectrum crossings is renormalized due to interactions.
The focus of the present work is to explore theoreti-cally signatures of non-Hermitian topology in the time evolution of the entanglement spectrum. A natural next step is to devise experimental implementations of the physics discussed in our work. The required building blocks are available in synthetic quantum systems: First, realizations of Kitaev chain have been proposed in systems of cold atoms [80][81][82][83]. Alternatively, the transverse field Ising model, to which the Kitaev chain is mapped through the Jordan-Wigner transformation [84], can be implemented in quantum simulators of many-body spin systems such as trapped atomic ions [85], Rydberg atom arrays [86], or superconducting qubits [87,88]. Further, Markovian dissipation, which is described by a quantum master equation with Hermitian jump operators, can be induced by subjecting the system to classical noise. This approach has been used, e.g., to implement local dephasing of individual spins in a recent experiment with trapped ions [89]. Finally, the measurement of entanglement spectra in various experimental platforms has been developed in Refs. [90][91][92]. These developments open up the prospect to experimentally explore non-Hermitian topology of driven and open quantum many-body systems, and thus go fundamentally beyond the paradigm of classical non-Hermitian wave physics. A comprehensive understanding of the robust signatures of non-Hermitian topology in quantum many-body systems is central for this endeavour. Our work, which identifies entanglement spectrum crossings as such a signature, represents an important step in this direction. quasilocal jump operators, there is a direct phase boundary between the real-line gapped topological R top and trivial R tr phases, even at finite values of the dissipation rate κ. Second, the model with quasilocal jump operators does not exhibit phases with an imaginary line gap. We proceed with a concise summary of the key properties of the driven-dissipative Kitaev chain with quasilocal jump operators.
(A2) We obtain two complex bands which are given by For 4κ < |2J − µ|, the expression in brackets is positive for all values of the quasimomentum k, and the spectrum exhibits a real line gap. Spectra in the two phases R top and R tr with a real line gap are shown in Figs. 11(a) and 11(b), respectively. Figure 11. Liouvillian spectra of a driven-dissipative Kitaev chain with quasilocal jump operators and for J = ∆. The notation and symbols are the same as those in Fig. 3. The spectrum exhibits a real line gap in (a) the topological phase Rtop, shown here for κ = 0.4J and µ = −J, and (b) the trivial phase Rtr, shown for µ = −3J and κ = 0.4J. (c) Gapless spectrum in G edge at µ = −J and κ = 0.8J.
As illustrated in Fig. 11(c), the spectrum is gapless for 4κ > |2J − µ|. In the phase diagram in Fig. 10, the two gapless phases are designated by G edge and G tr . Unlike the case of local jump operators, the spectrum remains gapless even for high values of κ, and the present model does not have phases with an imaginary line gap.
The symmetries stated in Eq. (19) remain intact for quasilocal jump operators, and the system still belongs to the class BDI † . We omit details of the straightforward calculation of the non-Hermitian winding number W for brevity. This calculation yields W = 1 in R top and W = 0 in R tr . Further, the global Berry phase is given by Q = 1 for |µ| < 2J and Q = 0 for |µ| > 2J. Edge modes, which occur in R top and G edge where Q = 1, are indicated by black diamond symbols in Fig. 11. At the critical lines µ = ±2J, where the edge modes disappear, the point gap of the bands λ ±,k at λ = 0 closes. A notable difference between the spectra shown in Figs. 3 and 11 is the wide dispersion of imaginary parts Im(λ ±,k ) which occurs in the case of quasilocal jump operators. The coexistence of modes with widely disparate decay rates hampers the numerical analysis of the present model. In particular, the time evolution of entanglement spectra, which we discuss in the next section, is hard to resolve for values of κ close to κ = J within the topological phase R tr in Fig. 10. To mitigate this problem, we consider dynamics generated by a band-flattened Liouvillian. Band-flattening can be achieved through a continuous deformation of the original Liouvillian, and, therefore, does not affect its topology [8].
Our starting point is the third-quantized form of the LiouvillianL given in Eq. (8). As explained in Sec. III B, the spectrum of iL is determined by the matrix Z = −iX . To resolve the problem of disparate decay rates, it is sufficient to consider a half-flatttened Liouvillian, which we obtain by replacing Z = V ΛV −1 , with Λ = diag(λ 1 , . . . , λ 2N ) a diagonal matrix, by Z hf = V Re(Λ)V −1 , i.e., by removing the decay rates − Im(Λ) by hand. We note that the spectrum becomes completely flat if Z is replaced by Z fl = JV sgn(Re(Λ))V −1 , where the prefactor J is introduced as an arbitrary overall energy scale. However, Z fl is well-defined only if there are no eigenvalues of Z with vanishing real part. The timeevolved entanglement spectra which we discuss below, as well as the data for t 1 shown in Fig. 10, are obtained with a half-flattened Liouvillian.

time evolution of entanglement spectra
Entanglement spectrum crossings reveal non-Hermitian topology of the Liouvillian also for a driven-dissipative Kitaev chain with quasilocal Hermitian jump operators. This is illustrated in Figs. 12(a) and 12(b) for quenches to R top and R tr , respectively. For time evolution with a half-flattened Liouvillian, the entangement spectrum does not show significant decay, and most entanglement eigenvalues remain close to ±1.
A quantitative analysis of the time t 1 of the first entanglement spectrum zero crossing within the entire phase R top is shown in Fig. 10. The logarithmic color scale indicates that t 1 diverges at the boundaries of R top . This is illustrated further in the inset, which shows the increase of t 1 upon approaching the phase boundary along the black arrows in the main panel, i.e., for δµ = µ 1 +2J → 0 with κ = 0.4J, and for µ 1 = J with δκ = J/4 − κ → 0. The numerical data is consistent with a square-root singularity, and, therefore, corroborates the universality of the value = 1/2 of the dynamical critical exponent for entanglement spectrum crossings. To find the covariance matrix of the ground state of the Kitaev chain, we start from the Majorana representation of the Hamiltonian given in Eq. (7). The matrix A in this equation is real and antisymmetric. Therefore, there exists a real special orthogonal matrix O which brings A to the following canonical block-diagonal form: where i ≥ 0. The matrix O can be found by employing the algorithm described in Ref. [93]. Then, the covariance matrix of the ground state is given by [72] Appendix C: Exactly solvable quench dynamics Here we derive the exact solution in Eq. (28) for the time evolution of the reduced density matrix. We obtain this exact solution for a specific choice of pre-and postequench parameters as detailed in the following.

Initial state
We assume that the system is initialized in the ground state |ψ 0 of the Hamiltonian of the Kitaev chain for J = ∆ > 0 and µ 0 → −∞, which is the vacuum of Dirac fermions c i , i.e., |ψ 0 = |Ω with c i |Ω = 0 for all i = 1, . . . , N . The corresponding density matrix is Here and below, we use the following notation for projectors on states for which the lattice site i is empty and occupied, respectively: where the site occupation number operator is n i = c † i c i .

Postquench Liouvillian
For the choice of postquench parameters J = ∆, µ 1 = 0, and γ l = γ g = γ, the third-quantized Liouvillian (8) takes the form where the site occupation number superoperator is defined asn i =ĉ † iĉ i . Left and right edge modes are created by the superoperatorsĉ † L =ĉ † 1 andĉ † R =ĉ † 2N , respectively. The corresponding eigenvalues of iL are λ L = −i4γ and λ R = 0. Since Re(λ L ) = Re(λ R ) = 0, the dynamics of these modes is purely dissipative. In contrast, the dynamics of bulk modes has both oscillatory and dissipative components. The bulk spectrum is flat and consist of two values λ ± , which are determined by the eigenvalues of the matrix in the second term in Eq. (C3), and are given in Eq. (26).
To separate the dissipative and oscillatory contributions to the dynamics of the bulk, we decompose the Liouvillian asL =L 0 +L 1 , wherê and .
(C6) With the aid of anticommutation relation of superoperators stated in Sec. III B, it is straightforward to check thatL 0 andL 1 commute. Therefore, the time evolution superoperator factorizes as eL t = eL 1t eL 0t .
We proceed to analyze first the dissipative, and then the oscillatory component of the dynamics.

Dissipative dynamics
The LiouvillianL 0 in Eq. (C4) is a sum of commuting termsL 0,i , which describe the evolution of fermions on lattice site i. Further, since the initial density matrix (C1) factorizes in real space, we obtain and, therefore, we can consider each site i independently. We begin by collecting a few useful relations: First, by expanding the exponential in a power series and using thatn 2 Next, from the Majorana representation of the site occupation number operator, n i = 1 2 (1 − iw 2i−1 w 2i ), it follows that With these relations, it is straightforward to show that where the projectors P 0 i and P 1 i are defined in Eq. (C2). Thus, the evolution (C8) of the density matrix under the local LiouvilliansL 0,i is given by where γ 1 = 3γ, γ N = γ, and γ i = 2γ for i = 2, . . . , N − 1.
(C18) In particular, for t = mT , where m ∈ N 0 is a nonnegative integer, we obtain eL (1) and, therefore, which can also be written in the compact form The oscillatory evolution of the entire chain is thus given by the product whereP is the parity superoperator which is defined aŝ Finally, we note that the evolution superoperator can be simplified for even and odd values of m, respectively: eL 1mT = 1 for m even, 1 − 2n 1 1 − 2n 2N P for m odd. (C24)

Time evolution of the density matrix
We now combine Eq. (C13) with Eq. (C24) to find the time-evolved density matrix ρ(mT ). According to Eq. (C24), at times mT , not only the dissipative evolution underL 0 , but also the oscillatory dynamics generated byL 1 factorizes in real space: To determine ρ i (mT ) explicitly, we have to apply super-operators1 − 2n i , which occur in Eqs. (C23) and (C24), to the projectors P n i in Eq. (C13). This can be done with the aid of the relations which follow from Eqs. (C2) and (C10), and which immediately lead to for n = 0, 1. We thus find where m i = m for i = 1, N , and m i = 0 for i = 2, . . . , N − 1. Each of the reduced density matrices for a single lattice site i has trace one, tr i (ρ i (mT )) = 1. Therefore, the reduced density matrix for the left half of the chain is simply given by the product where the sum in the second equality is over all combinations n = n 1 , . . . , n N/2 of site occupation numbers, and ξ m,1 and ξ m,2 are given in Eq. (29). Finally, we obtain Eq. (28) by noting that P n = |n A n| = N/2 i=1 P ni i , and by replacing the sum over site occupation numbers n in Eq. (C29) by a sum over entanglement occupation numbers e = e 1 , . . . , e N/2 with e 1 = (n 1 + m) mod 2 and e i = n i for i = 2, . . . , N/2. Inverting the relation between e and n leads to Eq. (30).

Appendix D: Conservation of parity of entanglement eigenstates
In the following, we show that the fermion parity of any individual entanglement eigenstate is conserved. To this end, we note first that the fermion parity operator P = e iπ 2N i=1 ni is related to the parity superoperator (C23) viaP ρ = P ρP . This can be checked by expanding the density matrix ρ in the basis of products of Majorana operators specified in Sec. III B. The relation betweenP and P implies that the density matrix ρ = |ψ ψ| for any pure state |ψ which has definite positive or negative parity has positive superparity. In particular, this applies to the ground state of the Hamiltonian of the Kitaev chain (1), which is the initial state of the quench protocol we consider in the present work.
The superparity of the initial state is conserved under time evolution with a quadratic Liouvillian [64]. This follows from the equality [P ,L] = 0, which further implies that the density matrix can be diagonalized in a basis of states with definite fermion parity at all times, ρ = n c n |ψ n ψ n |, P |ψ n = p n |ψ n , with p n = ±1. Consequently, also the reduced density matrix ρ A = tr A c (ρ) for a subsystem A can be diagonalized in a basis of eigenstates of the corresponding fermion parity operator P A , which is defined in Eq. (32), This can be seen as follows: A single state |ψ n , which occurs in Eq. (D1), can be written as where the states |e A,A c α and |o A,A c α form a basis of states which contain, respectively, an even and odd number of fermions in A and its complement A c , and for concreteness we assume that |ψ n has even parity. Upon taking the trace over A c , we find That is, the reduced density matrix ρ n A is composed of two disconnected blocks corresponding to states with even and odd particle number, respectively. It follows that P A ρ n A P A = ρ n A . An analogous argument applies to states |ψ n with odd parity and, therefore, to all states in the sum in Eq. (D1). This proves Eq. (D2), according to which entanglement eigenstates have definite parity at all times.
averages can be obtained with the aid of the quantum regression theorem [75], which yields for t > 0 w i (t)w j = tr(w i e Lt (w j ρ ss )), w j w i (t) = tr(w i e Lt (ρ ss w j )), and, therefore, χ R i,j (t) = 1 2 θ(t) tr(w i e Lt {w j , ρ ss }) = θ(t) tr(w i e Lt w j ).
(E2) In the second equality, we specified the trivial steady state ρ ss = 1/D. A closed equation of motion can be derived for the response matrix, which collects all possible combinations of odd and even position indices, . (E3) The equation of motion, which results from taking the derivative of Eq. (E2), takes a particularly simple form in momentum space: where R k (t) = j∈Z e −ikj R i+j,i (t), x k = x k − 4 tr(M ), and x k = iz k . For the local dissipative model, z k is given in Eq. (13), and tr(M ) = γ. The solution to Eq. (E4) reads which yields the equal-position retarded response function as To obtain the late-time asymptotics of the retarded response function quoted in Sec. V B and illustrated in Fig. 7, we evaluate the integral over quasimomenta using standard techniques of asymptotic analysis [94]. The numerical results shown in Fig. 7 result from a direct numerical integration of Eq. (E6).
Appendix F: Even-odd half-system-size effect In this appendix, we derive Eqs. (50) and (51), which describe the even-odd half-system-size effect of the fermion parity of the entanglement ground state, and the sign of the Pfaffian.

Even-odd effect of the fermion parity
Before we confirm Eqs. (50) and (51) numerically for arbitrary values of µ 0 and δ, it is instructive to consider certain limiting cases, in which the fermion parities p A,0 initial steady vacuum filled δ → ∞ δ → −∞ pA,ss = 1 pA,ss = (−1) N/2 vacuum crossings are µ0 → −∞ no crossing present/absent if pA,0 = 1 N/2 is odd/even filled crossings are µ0 → ∞ present/absent if no crossing pA,0 = (−1) N/2 N/2 is odd/even Table I. Even-odd half-system-size effect for the parity in the initial and steady state, and, consequently, for the presence or absence of entanglement spectrum crossings. The initial and steady states are the vacuum state or the completely filled state, |Ω or N i=1 c † i |Ω , respectively, for the given limiting values of µ0 and δ. and p A,ss of the entanglement ground states of the initial and steady states, respectively, can be determined in a straightforward manner.
For J = ∆ > 0 and µ 0 → −∞, the ground state of the Hamiltonian (1) of the Kitaev chain is the vacuum state |Ω , which obeys c i |Ω = 0 for all i = 1, . . . , N ; for µ 0 → ∞, the ground state is the completely filled state N i=1 c † i |Ω . In the lattice-site occupation representation introduced in Sec. V A 2, the corresponding entanglement ground states are |0 A and |1 A , where 1 = (1, . . . , 1). For the vacuum state, the fermion parity p A,0 of the entanglement ground state is even, p A,0 = 1. In contrast, for the completely filled state, the parity depends on the size of the system as p A,0 = (−1) N/2 , i.e., the parity is even or odd if half of the system size is even or odd, respectively. As in the main text, we consider systems with an even number of lattice sites, N ∈ 2N, such that N/2 ∈ N is an integer.
With regard to the steady state, we consider the following limits of the value of δ = γ l − γ g : For δ → ∞, the steady state is the vacuum state with p A,ss = 1; for δ → −∞, the steady state is the completely filled state with p A,ss = (−1) N/2 .
If the difference of parities is finite, p A,0 − p A,ss = 0, in the course of the quench dynamics, there has to be at least one reversal of the parity of the entanglement ground state, and, concomitantly, a zero crossing in the single-particle entanglement spectrum. Our findings for p A,0 and p A,ss , and the resulting expectations for the presence or absence of entanglement spectrum crossings, are summarized in Table I.

Pfaffian in the initial and steady state
We proceed to confirm numerically and for generic values of pre-and postquench parameters, that the sign of the Pfaffian Pf (24) of the reduced covariance matrix of the initial and steady states exhibits the same half- system-size dependence as the fermion parity p A of the corresponding entanglement ground state. Indeed, in all cases we consider, we find sgn(Pf) = p A . The parity can be calculated as described in Appendix G. The even-odd half-system-size effect for the Pfaffian is illustrated in Fig. 13. In particular, Fig. 13(a) shows the Pfaffian of the reduced covariance matrix of the ground state of the Kitaev chain with open boundary conditions as a function of the chemical potential µ 0 . For increasing system size, the Pfaffian Pf 0 develops a sharp drop to zero at the critical values µ 0 = ±2J. In the topological phase of the Kitaev chain, for |µ 0 | < 2J, the Pfaffian Pf 0 is equal to zero, which reflects the presence of a zeroenergy edge mode at the entanglement cut. Indeed, since the square of the Pfaffian equals the determinant of the Figure 14. The gap ∆ξ of the single-particle entanglement spectrum approaches its thermodynamic value ∆ξ∞ exponentially. The orange line is a fit to the numerical data shown as blue dots. Parameters are γ l = 0.3J, γg = 0.7J, and µ1 = 0.
reduced covariance matrix, Pf 0 = 0 if one of the singleparticle entanglement eigenvalues ξ i is equal to zero. In the trivial phase with |µ 0 | > 2J, the sign of Pf 0 is given by Eq. (51).
In Appendix F 1, we argued that the fermion parity of the entanglement ground state in the steady state is always positive for δ → ∞, whereas it exhibits an even-odd half-system-size effect in the limit δ → −∞. As shown in Fig. 13(b), the Pfaffian in the steady state exhibits an analogous even-odd half-system-size effect also for finite values of the imbalance of dissipation rates δ = γ l − γ g . This confirms Eq. (50) for the values γ l + γ g = J and µ 1 = 0, chosen in the figure. We checked that Eq. (50) holds also for nonzero values of µ 1 .
The dependence of Pf ss on the half-system size N/2 is illustrated further in the inset of Fig. 13(b), which shows Pf ss for δ = −0.4J and a series of half-system sizes N/2 = 10, 11, 20, 21, . . . , 50, 51. While the linear behavior of Pf ss on the logarithmic scale clearly indicates exponential decay of the magnitude of the steady-state Pfaffian, the fermion parity of the entanglement ground state, and, therefore, the sign of the Pfaffian, is well-defined for arbitrarily large values of N . This follows from the fact that for δ = 0, the single-particle entanglement spectrum retains a finite gap ∆ξ > 0 for N → ∞, and that a finite gap implies a unique entanglement ground state with well-defined parity. In the limits δ → ±∞, the entanglement gap is ∆ξ = 1 for the vacuum and completely filled state, respectively, as follows from the explicit expressions Γ vac For generic parameter values, we present numerical evidence for a finite entanglement gap in Fig. 14.
Appendix G: Fermion parity of the entanglement ground state A numerical evaluation of Eq. (22) yields the reduced covariance matrix Γ A at arbitrary times t. Here we describe how the parity of the corresponding entanglement ground state can be determined. To this end, we use the following representation of the parity operator (32) Our goal is to calculate the expectation value of P A in the entanglement ground state. Since the entanglement Hamiltonian is quadratic, all of its eigenstates are Slater determinants. Therefore, expectations values of products of fermionic operators as in Eq. (G1) can be found by employing Wick's theorem, which yields [72] P A = pf(Γ 0 A ), where Γ 0 A is the covariance matrix of the entanglement ground state, i.e., the ground state of the entanglement Hamiltonian. The latter can be written as [72] where the real and antisymmetric matrix G can be expressed in the form Here, O is the real special orthogonal matrix which brings the reduced covariance matrix Γ A to its canonical blockdiagonal form, and ε i = 2 atanh(ξ i ). The covariance matrix Γ 0 A of the ground state of the entanglement Hamiltonian (G2) can be found as described in Appendix B for the original Hamiltonian (1), and yields the parity of the entanglement ground state as described above.