Quantum quenches in driven-dissipative quadratic fermionic systems with parity-time symmetry

We study the quench dynamics of noninteracting fermionic quantum many-body systems that are subjected to Markovian drive and dissipation and are described by a quadratic Liouvillian which has parity-time (PT) symmetry. In recent work, we have shown that such systems relax locally to a maximum entropy ensemble that we have dubbed the PT-symmetric generalized Gibbs ensemble (PTGGE), in analogy to the generalized Gibbs ensemble that describes the steady state of isolated integrable quantum many-body systems after a quench. Here, using driven-dissipative versions of the Su-Schrieffer-Heeger (SSH) model and the Kitaev chain as paradigmatic model systems, we corroborate and substantially expand upon our previous results. In particular, we confirm the validity of a dissipative quasiparticle picture at finite dissipation by demonstrating light cone spreading of correlations and the linear growth and saturation to the PTGGE prediction of the quasiparticle-pair contribution to the subsystem entropy in the PT-symmetric phase. Further, we introduce the concept of directional pumping phases, which is related to the non-Hermitian topology of the Liouvillian and based upon qualitatively different dynamics of the dual string order parameter and the subsystem fermion parity in the SSH model and the Kitaev chain, respectively: Depending on the postquench parameters, there can be pumping of string order and fermion parity through both ends of a subsystem corresponding to a finite segment of the one-dimensional lattice, through only one end, or there can be no pumping at all. We show that transitions between dynamical pumping phases give rise to a new and independent type of dynamical critical behavior of the rates of directional pumping, which are determined by the soft modes of the PTGGE.


I. INTRODUCTION
Quantum quenches are a paradigmatic scenario for inducing and studying far-from-equilibrium dynamics in isolated quantum many-body systems.In a quantum quench, an initial state, often chosen to be the ground state of a prequench Hamiltonian, is evolved in time with a postquench Hamiltonian [1,2].Of particular interest are quenches across phase boundaries, which can lead to a variety of intriguing phenomena, especially in systems with nontrivial topology [3][4][5][6][7][8][9].But the two basic questions that underlie the interest in nonequilibrium dynamics induced by quantum quenches read as follows: Does a given system equilibrate after a quench, i.e., do expectation values of local observables reach a steady state?And what is the nature of this steady state?It is well established that integrable quantum many-body systems, which are characterized by an extensive set of integrals of motion, relax locally to a statistical ensemble determined by the principle of maximum entropy [10]-the generalized Gibbs ensemble (GGE) [11][12][13][14][15][16][17].The GGE can be regarded as being universal in the sense that its general structure is model-independent, whereas specific model-dependent details enter the GGE only through the form of the integrals of motion and as Lagrangian multipliers, the values of which are determined by the state in which the system is prepared before the quench.In this way, the GGE conserves an extensive amount of information about the initial state.Relaxation to the GGE is accompanied by a number of universal characteristic features, such as light cone propagation of correlations [16,18], linear growth and volume-law saturation of the subsystem entropy [19][20][21][22] and the equilibration of local order parameters [13][14][15][16][17].
In stark contrast, open many-body quantum systems, which are subjected to Markovian drive and dissipation, generally approach a highly model-dependent steady state that is determined by the interplay between internal Hamiltonian dynamics and the coupling to external reservoirs [23][24][25][26][27][28][29][30][31][32].Through these couplings, the system can exchange energy and particles with the reservoirs, breaking conservation laws the system would have in isolation.In particular, for integrable systems, this generically means that local integrals of motions of an isolated system are not conserved anymore if the system is open.As a result, after a quench, all memory of the initial state fades away with time, and the steady state takes the form of a GGE only in the limit of weak coupling γ between the system and its environment [33][34][35][36].Yet, as we have shown in Ref. [37], the universal principles governing generalized thermalization after a quantum quench in an isolated integrable system do apply-in suitably generalized formto specific driven-dissipative systems even for finite systembath coupling γ; and further, key signatures, such as a drivendissipative quasiparticle picture [38][39][40] and equilibration of local observables, still accompany relaxation to a maximum entropy ensemble.This unexpected robustness of generalized thermalization is ensured by parity-time (PT) symmetry of the Liouvillian that generates the dynamics of these drivendissipative systems.To highlight the fundamental roles that are played by PT symmetry and the principle of maximum entropy, we have dubbed the ensemble such systems locally relax to the PT-symmetric GGE (PTGGE) [37].
Originally, PT symmetry has been studied as a framework for a non-Hermitian generalization of quantum mechanics.In conventional quantum mechanics, physical observables are represented by Hermitian operators, mainly because these operators have real spectra.However, as shown in seminal work by Bender et al. [41], symmetry under the combination of spatial inversion or parity and time reversal can also lead to real spectra in non-Hermitian operators.In particular, an eigenvector of a PT-symmetric non-Hermitian operator is associated with a real eigenvalue if also the eigenvector itself is PT-arXiv:2304.01836v2[cond-mat.stat-mech]19 Feb 2024 symmetric, i.e., invariant under the PT transformation.Eigenvectors that are not invariant under the PT transformation are said to spontaneously break PT symmetry and are associated with complex eigenvalues.Typically, if a PT-symmetric operator depends on a parameter γ that measures the degree of non-Hermiticity such that the operator is Hermitian for γ = 0, all eigenvectors are PT-symmetric and the spectrum is entirely real if γ is sufficiently small; this situation is referred to as the PT-symmetric phase.At intermediate values of γ, in the PT-mixed phase, PT-symmetric and PT-breaking eigenvectors coexist.And for large values of γ, typically all eigenvectors spontaneously break PT symmetry.Non-Hermitian generalizations of quantum mechanics, as envisioned originally in Refs.[41], are restricted to the PT-symmetric phase with a real spectrum.More recently, PT-symmetric non-Hermitian versions of paradigmatic models of condensed matter theory such as the Kitaev chain [42] have received great interest due to their unconventional topological properties, in particular, in connection with the spontaneous breaking of PT symmetry that leads to the occurrence of complex eigenvalues [43][44][45][46][47][48][49][50][51].
Here, however, we are interested in PT symmetry of a special type of non-Hermitian operator: the Liouvillian that generates the dynamics of an open quantum many-body system [52][53][54][55][56][57][58][59][60][61][62][63][64].The coupling to external reservoirs generically induces exponential decay, which is reflected in the spectrum of the Liouvillian being complex.On the one hand, this implies that not only the PT-symmetric but also the PT-mixed and PT-broken phases can describe quantum dynamics of driven-dissipative systems; on the other hand, this puts into question the very existence of a PT-symmetric phase, since all eigenvalues of a generic Liouvillian become real only in the limit of vanishing dissipation-we choose here the convention that the real and imaginary parts of an eigenvalue of a Liouvillian determine, respectively, the oscillation frequency and decay rate of the corresponding eigenmode.However, after a shift by a constant decay rate, the single-particle spectrum of a quadratic Liouvillian, which describes a noninteracting open quantum many-body system, can become entirely real-a property that is known as passive PT symmetry [65,66].
Passive PT symmetry of quadratic Liouvillians is the key feature that underlies local relaxation to the PTGGE of the translationally invariant driven-dissipative fermionic lattice systems that we study below.In the PT-symmetric phase, the single-particle eigenvalues of such Liouvillians form two bands that are given by λ ±,k = −iκ ± ω k , where k is the quasimomentum.The momentum-independent decay rate κ results in temporally uniform overall exponential relaxation.Independently from that, dephasing of modes with different frequencies ω k ω k ′ leads to local relaxation to the PTGGEthis is completely analogous to isolated two-band models with single-particle dispersion ±ε k , where generalized thermalization to the GGE is induced by dephasing of purely oscillatory modes with ε k ε k ′ [16,67].However, note that the PTGGE is intrinsically time-dependent due to the exponential decay at the rate κ.Relaxation of an isolated system to the GGE and of a driven-dissipative system to the PTGGE, and the underlying single-particle spectra, are illustrated in Figs.1(a) and (b), respectively.After a quench in an isolated Figure 1.(a) Schematic representation of (left) the single-particle spectrum λ k of an isolated system (blue lines) with κ = 0 and (right) relaxation of an observable ⟨O ℓ ⟩ acting on ℓ sites to the GGE (red, dashed line) on a time scale t F (purple, dashed line).(b) For a drivendissipative system with PT symmetry, (left) the spectrum acquires a global shift by −iκ, and (right) relaxation to the PTGGE (red, dashed line) is revealed by rescaling ⟨O ℓ ⟩ with e 2ℓκt .(c) Quenches from the trivial to the topological phase are accompanied by oscillatory dynamics of a topological disorder parameter are shown for (left) an isolated system, with zeros of ⟨O ℓ ⟩ at multiples of a single soft-mode period t s (green vertical lines), and (right) a driven-dissipative system, where directional pumping is characterized by two soft-mode periods t s,− (dark green) and t s,+ (light green).system, a suitably chosen observable ⟨O ℓ ⟩ that acts on ℓ contiguous lattice sites equilibrates to the value predicted by the GGE after the Fermi time t F ∼ ℓ [14].In contrast, after a quench to the PT-symmetric phase of a driven-dissipative system, the global shift of the spectrum by −iκ results in exponential decay, ⟨O ℓ ⟩ ∼ e −2ℓκt .Relaxation to the PTGGE through dephasing can thus be revealed by considering the rescaled expectation value e 2ℓκt ⟨O ℓ ⟩.For stronger coupling to the environment, PT symmetry is broken spontaneously and the imaginary part of λ ±,k becomes dispersive.The long-time dynamics is then determined not by dephasing of all single-particle eigenmodes but rather by the single slowest-decaying mode that corresponds to the smallest value of − Im(λ ±,k ).Therefore, the spontaneous breaking of passive PT symmetry defines a sharp dynamical transition that delimits relaxation to the PTGGE and thus the validity of the principle of maximum entropy in driven-dissipative systems.An important caveat is that this argument for local relaxation to the PTGGE in the PTsymmetric phase is based solely on the form of the spectrum of the Liouvillian.Indeed, relaxation to the PTGGE applies as long as the overall decay with constant rate κ causes the system to heat up to infinite temperature.Instead, for a nontrivial steady state, local relaxation to the PTGGE occurs only transiently up to a sharply-defined crossover time t × [37].
In this work, we perform a detailed study of the quench dynamics and relaxation to the PTGGE of driven-dissipative fermionic many-body systems with PT symmetry, corroborat-ing and substantially expanding upon our earlier results that we have presented in Ref. [37].We consider two models, driven-dissipative versions of the Su-Schrieffer-Heeger (SSH) model [68] and the Kitaev chain [42], which are representative for broad classes of quadratic fermionic systems.Studying the time evolution of the dual string order parameter and the subsystem fermion parity, which serve as topological disorder parameters for the SSH model and the Kitaev chain, respectively, leads us to introduce directional pumping phases that are characterized by qualitatively different dynamics of the topological disorder parameters as illustrated in Fig. 1(c).
Using the example of a Kitaev chain with long-range hopping and pairing, we show that directional pumping phases and the dynamical critical behavior at transitions between these phases are in general independent from the phases and the associated criticality that are defined in terms of PT symmetry and gap closings of the postquench Liouvillian.
To make the differences between our previous [37] and current work more explicit, let us briefly summarize the additional results presented in this work that have not been presented before.We extend our previous results in two major ways.First, we present a detailed study of the quench dynamics of the driven-dissipative SSH model, which has been mentioned only briefly in Ref. [37].In particular, we derive the corresponding expression for the PTGGE and study the evolution of the subsystem entropy for this model.We complement our previous results for the Kitaev chain by studying the spreading of correlations.Second, we provide an extensive analysis of the dynamics and topological zero crossings of topological disorder parameters in both models.The most important new concepts and results in this context are directional pumping phases with dynamical critical behavior and a notion of universality at the phase boundaries.We further demonstrate modifications of this critical behavior induced by changes in the long-wavelength description due to long-range couplings.Finally, we compare the dynamical critical behavior of topological disorder parameters to the critical behavior of the connected density autocorrelation function.This paper is organized as follows: In Sec.II, we summarize our key results.The models we study, which are drivendissipative versions of the SSH and Kitaev chains, are introduced in Sec.III.We then establish the PTGGE as the maximum entropy ensemble for these models in Sec.IV.The spreading of correlations after a quench is studied in Sec.V, which is followed in Sec.VI by a discussion of the time evolution of the subsystem entropy.We study the dynamics of the dual string order parameter and the subsystem parity in Sec.VII, where we also introduce directional pumping phases and characterize the associated critical behavior.In Sec.VIII, we investigate directional pumping phases and phase transitions in a Kitaev chain with long-range hopping and pairing.Open research questions are presented in Sec.IX, and technical details are described in the Appendices A-F.

II. KEY RESULTS
We consider the following quench protocol: The system is initialized in the ground state |ψ 0 ⟩ of the Hamiltonian H 0 that describes an isolated SSH model or Kitaev chain, where we focus on quenches starting from the topologically trivial phases of these models.At t = 0, a parameter of the Hamiltonian is changed abruptly; simultaneously, the system is coupled to Markovian reservoirs.The ensuing dynamics is generated by a quantum Liouvillian L in Lindblad form, such that the state of the system at time t is given by ρ(t) = e −iLt ρ 0 , where ρ 0 = |ψ 0 ⟩⟨ψ 0 | is the initial state.Our main results can be summarized as follows: After quenches to the PT-symmetric phase, drivendissipative free fermionic models relax to a maximum entropy ensemble, the parity-time symmetric generalized Gibbs ensemble (PTGGE).We have introduced the PTGGE as the maximum entropy ensemble that describes relaxation of a driven-dissipative Kitaev chain and have briefly discussed quench dynamics of an SSH model with incoherent loss and gain in Ref. [37].Here, we present a detailed derivation of the PTGGE for the SSH model.As noted in the Introduction, dephasing of modes with different oscillation frequencies, but with a decay rate that is guaranteed to be equal for all modes by PT symmetry, is the fundamental process that underlies relaxation to the PTGGE in free fermionic systems.Due to the common decay rate of all modes, the PTGGE is inherently time-dependent.We choose the SSH model and the Kitaev chain as representatives of the two fundamental classes of quadratic fermionic systems: The isolated SSH model has a U(1) symmetry associated with the conservation of the number of particles, whereas due the presence of pairing terms in the isolated Kitaev chain only the fermion parity is conserved and the symmetry group is reduced to Z 2 .The driven-dissipative generalizations of these models combine their quadratic Hamiltonians with linear Lindblad operators, which enable the systems to exchange particles with external reservoirs and, therefore, always break particle number conservation.However, as detailed in Sec.III, we can still choose a particular form of dissipation so as to preserve a weak U(1) symmetry in the driven-dissipative SSH model [69].Then, as in the isolated SSH model, no anomalous correlations are generated in the course of the dynamics.Therefore, the drivendissipative versions of the SSH model and the Kitaev chain we consider in this work can be regarded as natural open-system generalizations of topological insulators and superconductors.
We illustrate relaxation to the PTGGE by considering the dynamics of topological disorder parameters, which take finite expectation values in the ground state in the trivial phase and vanish in the topological phase.Topological disorder parameters for the SSH model and the Kitaev chain are the dual string order parameter [70] and the subsystem fermion parity [37], respectively.The corresponding operators act nontrivially on ℓ contiguous lattice sites.As illustrated schematically in Fig. 1, relaxation of the topological disorder parameters to the PTGGE happens on a characteristic time scale given by the Fermi time t F ∼ ℓ [14].We further provide analytical conjectures for the evolution of the topological disorder pa-rameters, which we find to be in excellent agreement with our exact numerical results.These conjectures generalize analytical results for the quench dynamics of the isolated transverse field Ising model in the space-time scaling limit ℓ, t → ∞ with ℓ/t fixed [13][14][15], and show clearly how the dynamics are affected by drive and dissipation: Apart from the abovementioned uniform exponential decay, both the dispersion relation and the statistics of the eigenmodes of the adjoint Liouvillian, which generates the time evolution of operator expectation values, are modified, leading to pronounced quantitative differences even after rescaling of expectation values to compensate the exponential decay.
For the models we consider, a description of the late-time dynamics in terms of the PTGGE applies up to arbitrarily long times for balanced loss and gain, which leads to a steady state at infinite temperature.For a small imbalance, the PTGGE still provides an accurate description on intermediate time scales, up to a crossover time scale t × that scales logarithmically with the difference between loss and gain rates.However, some observables such as the dual string order parameter are not affected at all by an imbalance between loss and gain.
Correlations show ballistic light cone spreading in the PTsymmetric phase.In contrast, after quenches to the PT-mixed and PT-broken phases, correlations spread diffusively.The single-particle spectra for isolated and PT-symmetric drivendissipative systems shown in Figs.1(a) and (b), respectively, suggest that PT-symmetric quadratic Liouvillians admit a notion of quasiparticles that propagate coherently with a velocity v k = dω k /dk, as is also the case for quasiparticle excitations of an isolated system, but have a finite lifetime ∼ 1/κ.Based on this dissipative quasiparticle picture, we can expect many characteristic features of the dynamics of isolated systems to carry over to PT-symmetric driven-dissipative systems in the PT-symmetric phase.In particular, we find that the spreading of correlations after quenches to the PT-symmetric phase is described by a clear light cone structure.Interestingly, the speed at which correlations propagate is increased as compared to isolated systems; however, the finite lifetime of quasiparticles, which leads to an overall exponential decay of correlations, indicates that, in fact, it is not the case that in open systems more information is transported in a shorter time.
Light cone spreading of correlations is restricted to the PTsymmetric phase.But also after quenches to the PT-mixed and PT-broken phases there is a pronounced peak of correlations that propagates through the system.However, the peak position evolves diffusively rather than ballistically as in the PT-symmetric phase.
The growth and saturation of the subsystem entropy obeys the quasiparticle picture, adapted to driven-dissipative systems.For isolated systems, the quasiparticle picture leads to quantitative predictions for the full time evolution of the entropy of finite subsystems [19][20][21][22].One can regard the initial ground state as a source of pairs of entangled quasiparticles with opposite momenta, which propagate through the system with different velocities of at most v max .If one of the two quasiparticles that form a pair is located within the subsystem while the other one is outside of the subsystem, then this pair contributes to the entanglement between the subsystem and its complement and, therefore, to the subsystem entropy.Based on this picture, knowledge of the quasiparticle velocity and the stationary value of the subsystem entropy that is reached for t → ∞ is sufficient to determine the full time evolution of the subsystem entropy in the space-time scaling limit ℓ, t → ∞, where ℓ is the size of the subsystem.The quasiparticle picture has been extended to open systems, where the requirement of ballistically propagating quasiparticles restricts its applicability to the limit of weak dissipation γ → 0 [38,39], and an additional contribution to the subsystem entropy due to the mixedness of the state has to be accounted for.If, however, the system under consideration has PT symmetry, then, as explained above, ballistically propagating quasiparticles exist within the entire PT-symmetric phase.Based on this observation, in Ref. [37], we have proposed an analytical conjecture for the quasiparticle-pair contribution to the subsystem entropy of a PT-symmetric Kitaev chain in the space-time scaling limit and for finite dissipation strength γ.Here, we provide further evidence for broad validity of our conjecture by applying it to the SSH model with incoherent loss and gain, where we again find excellent agreement with numerical data.
Quantum quenches in driven-dissipative systems can give rise to the unique phenomenon of directional pumping of topological disorder parameters.The time scales of directional pumping are determined by the soft modes of the PTGGE.As detailed in Sec.VII, in the isolated SSH model and the Kitaev chain, crossing the boundary between the trivial and the topological phase in the course of a quench is reflected in pumping of topological disorder parameters.That is, the dual string order parameter and the fermion parity of a subsystem of size ℓ exhibit oscillatory decay, crossing zero at multiplies of a time scale t s as illustrated in the left panel in Fig. 1(c).The period of zero crossings t s is determined by the momenta at which the GGE Hamiltonian vanishes-soft modes of the GGE [14,16].In contrast, for quenches within the trivial phase, the disorder parameters show nonoscillatory decay.Since the Hamiltonians of the SSH model and the Kitaev chain commute with the respective topological disorder parameters, processes that change the dual string order parameter and the fermion parity occur not in the bulk of the subsystem but at the interfaces between the subsystem and its complement.Interestingly, for quenches in driven-dissipative systems, the rates at which topological disorder parameters are pumped through the left and right ends of a subsystem are different.Therefore, as illustrated in the right panel in Fig. 1(c), there are two distinct time scales t s,± for zero crossings.These time scales are determined by soft modes of the PTGGE.As shown in Ref. [37] for the driven-dissipative Kitaev chain, necessary conditions for the phenomenon of directional pumping to occur are opensystem dynamics leading to mixed states, and the breaking of inversion symmetry by the coupling to reservoirs.
A dynamical phase diagram can be defined in terms of directional pumping.The resulting phase boundaries and the dynamical critical behavior at these phase boundaries are, in general, different and independent from the corresponding properties defined in terms of gap closings and PT symmetry of the postquench Liouvillian.PT symmetry of the Liouvillian leads to the distinction between PT-symmetric, PT-breaking, and PT-mixed phases.For the models we consider here, the transition from the PT-symmetric to the PT-mixed phase occurs when the gap between the bands λ ±,k = −iκ ± ω k closes.Note that this corresponds to a purely dynamical phase transition that marks a qualitative change of only the coherent dynamics: As we demonstrate using the example of the density autocorrelation function, upon approaching the phase boundary from the PT-symmetric phase, a characteristic period of oscillations of local correlation functions diverges; in the PTmixed phase, the decay of the density autocorrelation function is overdamped.Further, this transition is purely dynamical in the sense that it does not affect the steady state.Indeed, for this part of our analysis, we focus on balanced loss and gain, such that the steady state is always at infinite temperature.
Here, we show that a different and independent characterization of dynamical phases can be given in terms of directional pumping: As explained above, the pumping of topological disorder parameters for quenches from the trivial to the topological phase in isolated systems becomes directional in open systems.In particular, the topological phases of the isolated SSH model and Kitaev chain are continuously connected to phases with pumping at different rates through both ends of a subsystem.Surprisingly, we find transitions to phases with pumping through only the left or the right end of a subsystem.As the transition to such a phase is approached, one of the time scales t s,± of zero crossings diverges.
The Kitaev chain with long-range hopping and pairing represents an especially interesting model to study directional pumping phases and the associated critical behavior.In particular, we find that long-range couplings can modify the critical exponents that govern the divergence of the time scales t s,± , whereas the exponents that describe the divergence of the period of oscillations of the density autocorrelation function at the gap closing at the transition from the PT-symmetric to the PT-mixed phase remain unchanged.Moreover, in the presence of long-range couplings, directional pumping phase boundaries do not always coincide with phase boundaries that are determined by gap closings, which implies that a divergence of one of the time scales t s,± does not require a gap closing in the spectrum of the Liouvillian.These findings establish directional pumping phases, phase transitions, and the associated critical behavior as new and independent concepts that are unique to quantum quenches in driven-open systems.

III. MODELS
For the quench protocol outlined above, the postquench time evolution of the system density matrix ρ is described by a Markovian quantum master equation in Lindblad form [71], where the Liouvillian L incorporates both unitary dynamics generated by the Hamiltonian superoperator H and Markovian drive and dissipation described by the dissipator D. These superoperators are defined through their action on the density matrix ρ, where L l are the quantum jump operators, describing the coupling between system and environment.In this work, we focus on noninteracting fermionic lattice models that are described by a Hamiltonian H and jump operators L l that are quadratic and linear in fermionic operators, respectively, leading to a quadratic Liouvillian L. For the Hamiltonian H, we study two models: the SSH model and the Kitaev chain, which are paradigmatic examples for a one-dimensional topological insulator and superconductor, respectively.We specify the Hamiltonians and jump operators in the following, and give a detailed description of the symmetries, spectrum, and time evolution for SSH model.For a detailed account of the drivendissipative Kitaev chain, we refer to Ref. [37].

A. Driven-dissipative SSH model
The SSH model with L unit cells is described by the following many-body Hamiltonian [68]: where c s,l and c † s,l are, respectively, annihilation and creation operators for fermions on sublattice s ∈ {A, B} at lattice site l.Further, J 1 , J 2 ∈ R >0 are the hopping amplitudes within and between unit cells.Unless stated otherwise, we assume periodic boundary conditions (PBC) with c s,L+1 = c s,1 ; open boundary conditions (OBC) are implemented by setting c s,L+1 = 0. We often find it convenient to parameterize J 1 and J 2 in terms of total and relative hopping amplitudes, As detailed below, the SSH model belongs to the Altland-Zirnbauer class BDI [72].Therefore, the topology of the SSH model is characterized by an integer-valued invariant, the winding number W, and the ground state of the SSH model is topologically trivial and nontrivial for ∆J > 0 and ∆J < 0, respectively [73][74][75].In this work, we study quench dynamics whereby the initial state |ψ(t)⟩ at t = 0 is chosen to be the ground state |ψ 0 ⟩ for prequench parameters J 1,0 > 0 and J 2,0 = 0, corresponding to the topologically trivial phase.The ground state |ψ 0 ⟩ is Gaussian, i.e., the density matrix ρ 0 = |ψ 0 ⟩⟨ψ 0 | can be written as the exponential of a quadratic form in fermionic operators; further, since the Hamiltonian Eq. ( 4) is quadratic, also the time-evolved state ρ(t) = e −iHt ρ 0 e iHt is Gaussian [76].Therefore, after a quench in the isolated SSH model, the state of the system is Gaussian at all times, and thus fully determined by the covariance matrix, for s, s ′ ∈ {A, B} and l, l ′ ∈ {1, . . ., L}, and where In particular, anomalous correlations vanish, Turning now to a drivendissipative generalization of the SSH model, we note first that the dynamics generated by a quadratic Liouvillian L preserves the Gaussianity of the time-evolved mixed state ρ(t) = e −iLt ρ 0 [76,77].The vanishing of anomalous correlations, which for the isolated SSH model is a consequence of particle number conservation, can be ensured in an open SSH model by imposing the weak U(1) symmetry that is defined by the relation L(S ρS † ) = S L(ρ)S † for S = e iθN where the particle number operator is N = s∈{A,B} L l=1 c † s,l c s,l and θ ∈ [0, 2π) [69].This weak symmetry condition guarantees that no correlations are established between Hilbert space sectors with different numbers of particles, and, therefore, anomalous correlations vanish.For the Hamiltonian part of the time evolution in Eq. ( 1), the weak symmetry is again a consequence of particle number conservation of the SSH Hamiltonian Eq. ( 4), [H, N] = 0. Note that the weak U(1) symmetry does not affect the integrability of the model, which is guaranteed by the Liouvillian being quadratic.The weak U(1) symmetry only ensures the vanishing of anomalous correlations, which leads to a simplification of the dynamics.A specific choice of dissipation that respects the weak U(1) is given by incoherent loss and gain as described by the dissipator D = D l + D g with 2L g,s,l ρL † g,s,l − {L † g,s,l L g,s,l , ρ} .
The jump operators corresponding to local loss and gain are given by where we consider loss and gain rates γ l,s and γ g,s , respectively, that depend on the sublattice index.We parameterize these rates in terms of their mean and difference, , .
Below it will prove convenient to express the loss and gain rates in terms of the four independent parameters γ, δ, ∆γ, and ∆, which are defined by To determine the dynamics of the covariance matrix Eq. ( 6), let us first represent the jump operators in the general form where we collect sublattice and lattice indices in a single index such that c A,l = c 2l−1 and c B,l = c 2l .With the L × 2L matrices B l and B g , we can then introduce the bath matrices Finally, the Liouvillian dynamics of a Gaussian state is described by the equation of motion for the covariance matrix which reads [78] where the generator of the dynamics Z can be interpreted as a non-Hermitian Hamiltonian and is defined by The derivation of Eq. ( 13) is presented in Appendix A. Before we proceed with the analysis of the drivendissipative SSH model, let us briefly comment on the relation between the formalism of Lindblad master equations on which our work is based, and the framework of non-Hermitian Hamiltonians [79].Formally, the matrix Z defined in Eq. ( 14) is equivalent to non-Hermitian SSH model of Refs.[48,80,81] up to a shift by −i2γ1.Such a shift does not affect the topological properties that are considered in these references.However, the shift guarantees that the eigenvalues of Z are negative semidefinite, which is a necessary requirement for the dynamics of the covariance matrix generated by Z to be stable and thus physically meaningful.More generally, and in contrast to non-Hermitian Hamiltonians, the master equation (1), from which the evolution equation (13) follows, provides a full description of the time evolution of an open quantum system.In particular, the time evolution generated by the Liouvillian L is completely positive and trace preserving.This guarantees that the time-evolved density matrix ρ(t) represents, at all times, a physical state, from which the expectation value of an observable O can be evaluated as ⟨O(t)⟩ = tr(Oρ(t)).In Appendix B we discuss how the left and right eigenvectors of Z enter ⟨O(t)⟩.For non-Hermitian Hamiltonians, these eigenvectors are used to define biorthogonal expectation values, and are of key importance for the characterization of topology [82].
In the following, we discuss some fundamental properties of the driven-dissipative SSH model.In particular, we examine symmetries of the isolated and driven-dissipative system, the corresponding spectrum and mode structure, non-Hermitian topology, and finally the dynamics of the covariance matrix.

Symmetries of the isolated and driven-dissipative SSH model
Symmetries of driven-dissipative systems are of fundamental importance not only for their topological classification but also-as highlighted, in particular, in our work [37]-for their quench dynamics.In our discussion of symmetries of the driven-dissipative SSH model, we will exploit translational invariance and work in momentum space.To that end, we first rewrite the Hamiltonian in Eq. ( 4) in terms of spinors with the 2 × 2 matrix h l = h l • σ where σ = (σ x , σ y , σ z ) ⊺ is a vector of Pauli matrices.The momentum-space representation of the SSH model is given by the Bloch Hamiltonian h k that is defined as where Similarly, also the bath matrices in Eq. ( 12) can be expressed in terms of translationally invariant 2 × 2 blocks: Combining the block representations of the Hamiltonian and the bath matrices, we can write the matrix Z given in Eq. ( 14) in the form Then, as in Eq. ( 16), we obtain the non-Hermitian Bloch Hamiltonian z k that generates the time evolution of the covariance matrix: with where êz = (0, 0, 1) ⊺ is a unit vector along the z axis.a. Isolated system.We first discuss symmetries of the isolated SSH model that is described by the Bloch Hamiltonian in Eq. ( 16).As stated above, the SSH model belongs to the Altland-Zirnbauer class BDI and has particle-hole symmetry (PHS), time reversal symmetry (TRS), and chiral symmetry (CS).Due to the Hermiticity of the Bloch Hamiltonian, h k = h † k , each of these symmetries can be expressed in two equivalent ways: PHS: which can be confirmed by using Another symmetry that will play a key role in the following is inversion symmetry (IS).In real space and when expressed as a transformation of the fermionic operators c s,l , inversion amounts to an exchange of the sublattices and a reflection across the middle of the chain, c A/B,l → c B/A,L+1−l .The combination of IS with TRS yields PT symmetry (PTS).In momentum space, IS and PTS are described by IS: b. Driven-dissipative system.For the driven-dissipative SSH model, we can regard the non-Hermitian matrix z k z † k in Eq. (20) as the generalization of the Bloch Hamiltonian to open systems.Of the two equivalent versions of the symmetries of the Bloch Hamiltonian h k stated in Eq. ( 22), only one of each applies to z k : Therefore, in the nomenclature of Ref. [83], the drivendissipative SSH model belongs to the class BDI † .Further, inversion symmetry in the form given in Eq. ( 23) is broken by the dissipative contributions to z k .However, an inversion symmetry IS † still applies to the traceless part with z 1 given in Eq. ( 21); and by combining IS † with TRS † we obtain PTS of z ′ k : This form of a PT symmetry which applies after a shift that renders the generator of the dynamics traceless is called passive PT-symmety [62,65,66].
Let us now discuss the implication of PTS for the mode structure of the Liouvillian [41].The eigenvalues and eigenvectors of the shifted matrix z ′ k are denoted by λ ′ ±,k = ± √ z k • z k and |ψ ±,k ⟩, respectively.Since z ′ k and z k are related by a shift z 1 1, they have the same eigenvectors, and their eigenvalues are related by λ ±,k = λ ′ ±,k + z 1 with z 1 = −i2γ.The PTS condition in Eq. ( 25) then leads to This implies that for any eigenvalue λ ′ ±,k of z ′ k correpsonding to an eigenvector |ψ ±,k ⟩, there is also an eigenvalue λ ′ * ±,k with eigenvector σ x |ψ ±,k ⟩ * .Since for a given value of k there are only two eigenvalues which are related by λ ′ +,k = −λ ′ −,k , there are two possibilities: (i) PT-symmetric eigenmodes with λ ′ ±,k = λ ′ * ±,k ∈ R and σ x |ψ ±,k ⟩ * = |ψ ±,k ⟩.For the eigenvalues λ ±,k of the full matrix z k this implies that Im(λ ±,k ) = −2γ and Re(λ +,k ) = − Re(λ −,k ).(ii) Alternatively, there exist PTbreaking eigenmodes with λ ′ ±,k = −λ ′ * ±,k ∈ iR and eigenvectors obeying σ x |ψ ±,k ⟩ * = |ψ ∓,k ⟩ such that Re(λ ±,k ) = 0 and Im(λ +,k ) + 2γ = − Im(λ −,k ) + 2γ .That is, the eigenvalues that are associated with PT-breaking eigenmodes are related by reflection across the line −i2γ.Depending on the symmetry of the eigenvectors |ψ ±,k ⟩ under the PT transformation, we can distinguish three different phases: (i) In the PT-symmetric phase, all eigenmodes are PTsymmetric; then, eigenvalues have a constant imaginary value, while the real part is dispersive.This is shown in Fig. 2(b).(ii) In the PT-broken phase, all eigenmodes are PT-breaking, and the corresponding eigenvalues are purely imaginary and dispersive as illustrated in Fig. 2(d).(iii) Finally, in the PTmixed phase, PT-symmetric and PT-breaking modes exist simultaneously.Therefore, as shown in Fig. 2(c), there are two sets of eigenvalues: one that is dispersive only in the real part, and one that is dispersive only in the imaginary part.The resulting dynamical phase diagram of the driven-dissipative SSH model is shown in Fig. 2(a).

Spectrum of the Liouvillian
We next consider the spectra of the isolated and drivendissipative SSH models with PBC.As shown in Appendix C, the Bloch Hamiltonian h k Eq. ( 16) can be diagonalized as where the unitary matrix U k is given by Eq. (C1) after setting γ = ∆γ = 0, and with the single particle dispersion relation ε k of the isolated SSH model defined by the magnitude of the Bloch vector, For future reference, we note that in terms of the fermionic operators d s,k defined through the Hamiltonian takes the following diagonal form: where the Brillouin zone is BZ = {−π + ∆k, −π + 2∆k, . . ., π} with ∆k = 2π/L.The ground state |ψ 0 ⟩ is obtained by filling the band with negative energy, where |Ω⟩ is the vacuum of particles.
Applying the formalism of third quantization [84], one can show that the spectrum of the Liouvillian is determined by the eigenvalues of the matrix z k in essentially the same way as the spectrum of the second-quantized Hamiltonian Eq. ( 4) is obtained by occupying single-particle states with energies ±ε k .The eigenvalues of z k , which thus can be regarded as forming the single-particle spectrum of the Liouvillian, form two bands and are given by with the dispersion relation where γ and ∆γ are defined in Eq. ( 10).The dispersion relation ω k determines the dynamical phase diagram of the driven-dissipative SSH model in the ∆J-∆γ plane depicted in Fig. 2(a).In particular, different phases can be defined in terms of the gap structure of the bands λ ±,k [83]-as we explain next, this is equivalent to the distinction of phases in terms of PT symmetry: For small values of ∆γ, the bands are separated by a real line gap (blue, red regions in the figure); that is, ω k > 0 for all values of k ∈ BZ, with a typical complex band structure shown in Fig. 2(b).This is the PTsymmetric phase.Upon increasing the value of ∆γ, the real line gap closes, and the spectrum of z k is gapless in a finite region of the ∆J-∆γ plane (green, yellow).Then, for a range of momenta, ω k is real; for all other values of k, ω k = iκ k is purely imaginary, with A typical band structure in this gapless or PT-mixed phase is illustrated in Fig. 2(c).The phase boundary between PTsymmetric and PT-mixed phases is determined by a critical value ∆γ c .For |∆J| < J, the critical value is ∆γ c = |∆J|, and for |∆J| > J, it is given by ∆γ c = J.At |∆J| = J, the dispersion is flat, ω k = 2J, which enables a direct transition between the phase with a real line gap and the phase with an imaginary line gap or PT-broken phase (orange, purple), where κ k > 0 for all values of k ∈ BZ.The band structure in this phase is exemplified in Fig. 2(d).

Non-Hermitian topology
We have already mentioned that the SSH model is in the Altland-Zirnbauer class BDI.The class BDI in one spatial dimension is topologically nontrivial and characterized by an integer-valued invariant, called the winding number W. For models with a real line gap, the Altland-Zirnbauer classification can be extended to non-Hermitian systems, where in the nomenclature of Ref. [83], the driven-dissipative SSH model belongs to the class BDI † .Accordingly, within the blue and red regions in Fig. 2(a), which correspond to the phase with a real line gap or, equivalently, the PT-symmetric phase, the topology of the driven-dissipative Kitaev chain is characterized by a non-Hermitian generalization of the winding number.To calculate the non-Hermitian winding number, we first find the left and right eigenstates of z k , which we denote by |ψ L ±,k ⟩ and |ψ R ±,k ⟩, respectively.By employing the representation of z k given in Eq. (C4), we find the eigenvalue equations where |±⟩ are the eigenvectors of σ z with eigenvalues ±1.The left and right eigenvectors of z k are thus given by respectively, where U − † k is the short-hand notation for the inverse of the Hermitian adjoint of U k .We can now define a projector on the − band: and determine the non-Hermitian winding number by where Γ = σ z is the matrix that appears in the chiral symmetry condition Eq. ( 24) [83], and we have replaced summation of k ∈ BZ by an integral over k ∈ [−π, π] in the thermodynamic limit L → ∞.With {Q k , σ z } = 0, following from chiral symmetry, and and the calculation reduces to In the PT-symmetric phase, we find q k = − J 1 + J 2 e −ik ω k , and, therefore, The second term on the right-hand side of this equation is an odd function of k and does not contribute to the symmetric integral in Eq. ( 40); the first term yields which is solved by substituting z = e ik and integrating over the unit circle |z| = 1.Simple poles of the integrand are located at z 1 = 0 and z 2 = −J 2 /J 1 .By using the residue theorem, we obtain the result This is essentially the same result as for the isolated SSH chain, but extended to the whole PT-symmetric phase.The topological PT-symmetric phase is indicated in Fig. 2(a) by the blue area with the label "TOP," and the trivial phase by the red area with the label "TRIV."Interestingly, while the definition of the non-Hermitian winding number is restricted to the PT-symmetric phase with a real line gap, edge modes with λ edge = 0, −i4∆γ occur in a chain with OBC when ∆J < 0 and for arbitrarily large values of ∆γ, including in the PTmixed and PT-broken phases.In Figs.2(b), (c), and (d), edge modes are indicated with black dots.The existence of these edge modes can be understood to be a consequence of chiral symmetry [85]: Due to chiral symmetry, each sublattice supports one edge mode.But on a given sublattice, the non-Hermitian contribution to Z in Eq. ( 14) is constant, and its presence does not affect the eigenvectors with support on that sublattice.Therefore, the edge modes of Z are identical to the edge modes of the Hamiltonian H, and their existence is determined by the topology of H.This concludes our discussion of the static properties of the SSH model, and we turn now to quench dynamics, which are the main interest of our work.

Dynamics of the driven-dissipative SSH model
As explained at the beginning of Sec.III A, in this work, we consider quench dynamics with the system initially prepared in the ground state given in Eq. ( 31), for prequench parameters J 1,0 > 0 and J 2,0 = 0 corresponding to the trivial phase.Gaussianity of this state is preserved in the time evolution generated by a quadratic Liouvillian, and, therefore, the state of the system is fully determined by the covariance matrix G defined in Eq. ( 6), whose dynamics are described by Eq. ( 13).For PBC, the covariance matrix is a 2L × 2L block Toeplitz matrix, which is built from 2 × 2 blocks given by Further, for PBC and due to translational invariance of the Hamiltonian, the Markovian baths and the chosen initial state, all matrices in Eq. ( 13) are block-circulant Toeplitz matrices.Therefore, the representation of the blocks of the covariance matrix in momentum space, which is obtained through a discrete Fourier transform, obeys the following equation of motion: with where m l,k and m g,k are the momentum-space representations of the blocks of the bath matrices defined in Eq. ( 18).The time-evolved covariance matrix in momentum space can be split into two contributions, where Here, g 1,k (t) encodes information of the initial condition g k (0); the second contribution g 2,k (t) has no counterpart in isolated systems and describes the approach to the steady state ρ SS for t → ∞, where L(ρ SS ) = 0.
To evaluate the initial value g k (0), let us first define spinors , where the operators c s,k and d s,k are defined in Eq. (29).Further, by rearranging Eq. ( 27) and defining an initial unit vector n0,k = h 0,k /ε 0,k , where h 0,k and ε 0,k are, respectively, the Bloch Hamiltonian and the single-particle dispersion relation for prequench parameters J 1,0 and J 2,0 , we find n0,k for the ground state |ψ 0 ⟩ in Eq. ( 31), we obtain In particular, for our choice of prequench parameters given by J 1,0 > 0 and J 2,0 = 0, we find With this result for the initial value, and after some algebraic simplifications, a general form for the time-dependent covariance matrix can be found [37].We first split Eqs. ( 49) and ( 50) into two contributions that are proportional to the identity and traceless, respectively, In the PT-symmetric phase, the components of g 1,k (t) read and where we have introduced vectors h ∥,k and h ⊥,k which are parallel and perpendicular to nk = h k /ε k , respectively, and the out-of-plane vector h o,k that is orthogonal to both n0,k and nk : , and semi-minor axis ε k /ω k along h o,k ; for an isolated system with γ = ∆γ = 0, this reduces to the wellknown precession of g 1,k (t) around h k .
The components of g 1,k (t) given in Eqs. ( 54) and ( 55) decay exponentially at a rate 4γ.Therefore, as indicated above, the steady state, which is determined by the limit t → ∞ of g k (t), is described by g 2,k (t) Eq. ( 50).In turn, g 2,k (t) is proportional to y k , which vanishes for balanced loss and gain, δ = ∆ = 0.The vanishing of all correlations indicates that balanced loss and gain lead to a steady state at infinite temperature, ρ SS = ρ ∞ = 1/2 2L .In contrast, for δ, ∆ 0, also g 2,k (t) 0 and the steady state is nontrivial.The integral in the expression for g 2,k (t) in Eq. ( 50) can be solved by elementary means, but we omit the lengthy result.
Finally, for future reference, we briefly discuss the consequences of PHS for g k (t).In particular, PHS of the initial state and PHS † of z k imply that where we omit the time argument to shorten the notation.In contrast, y k in Eq. ( 47) breaks PHS, but has TRS and commutes with σ z .Combined with PHS † of z k this leads to By inversion of Eq. ( 45), one then immediately finds

B. Driven-dissipative Kitaev chain
In Ref. [37], we have presented a detailed study of the quench dynamics of a driven-dissipative Kitaev chain with short-range hopping and pairing.We summarize key properties of this model in the following.These properties will form the basis for our discussion of new results that concern the spreading of correlations and the effects of long-range hopping and pairing in Secs.V and VIII, respectively.
The Hamiltonian of a Kitaev chain [42] of length L, with hopping matrix element J, pairing amplitude ∆, and chemical potential µ, reads where the fermionic annihilation and creation operators at lattice site l are c l and c † l , respectively, with canonical anticommutation relations {c l , c † l ′ } = δ l,l ′ and {c l , c l ′ } = {c † l , c † l ′ } = 0. Depending on the observable under study, we consider both periodic boundary conditions with c L+1 = c 1 and open boundary conditions with c L+1 = 0, which will be indicated accordingly.We assume that J and ∆ are positive and real, such that the Kitaev chain has TRS and belongs to the Altland-Zirnbauer class BDI.In the following, we keep J and ∆ as distinct parameters in most expressions.However, our main results are obtained for J = ∆.For this choice, the ground state of the Kitaev chain is topologically nontrivial for |µ| < 2J.In Sec.VIII, we will also study a generalization of the Kitaev chain that incorporates long-range hopping and pairing.
We now subject the Kitaev chain to Markovian drive and dissipation in the form of local particle loss and gain as described by the jump operators [85][86][87] with loss and gain rates γ l and γ g , respectively.A convenient parameterization of the coupling to Markovian reservoirs is obtained by introducing the mean and relative rates, The mean rate γ determines the overall strength of dissipation; and the relative rate δ can be interpreted as an effective inverse temperature in the sense that for δ = 0, the steady state ρ SS is the jump operators Eq. ( 61) describe particle loss, leading to a pure steady state, ρ SS = |Ω⟩⟨Ω| [37].
In this work, we consider quenches originating from the trivial phase of the isolated chain with µ 0 → −∞.The initial state is then given by the vacuum of fermions |Ω⟩.As in the case of the SSH model discussed above, Gaussianity of the state is preserved under the evolution generated by the quadratic Hamiltonian Eq. ( 60) and the linear jump operators Eq. ( 61), and the state is fully determined by two-point correlations.However, for the Kitaev chain, there are also nonvanishing anomalous correlations, ⟨c l c l ′ ⟩, ⟨c † l c † l ′ ⟩ 0, and, therefore, it is convenient to describe correlations by 2L real Majorana fermions instead of L complex Dirac fermions.The transformation to Majorana operators w l reads These operators obey the anticommutation relation {w l , w l ′ } = 2δ l,l ′ .The covariance matrix can then be defined using Majorana operators as where Quench dynamics can thus be described by the time evolution of Γ(t).For details we refer to Ref. [37].Finally, analogous to the phases of the SSH model discussed in Sec.III A 1, the driven-dissipative Kitaev chain features PT-symmetric, PT-breaking, and PT-mixed phases.
The phase boundaries can be defined in terms of the spectrum of the Liouvillian, determined by λ ±,k = −i2γ ± ω k with the dispersion relation where Of particular interest for this work is the PT-symmetric phase where ω k ∈ R >0 for all k ∈ BZ, realized for while the PT-breaking phase, with For values of 2 √ γ l γ g between the boundaries given in Eqs. ( 67) and ( 68), the system is in the PT-mixed phase.

IV. PT-SYMMETRIC GENERALIZED GIBBS ENSEMBLE
After the technical preliminaries of the previous section, let us now focus on the main purpose of this work, which is to study the quench dynamics and relaxation of noninteracting driven-dissipative fermionic models.We start by deriving the maximum entropy ensemble [10] for PT-symmetric free fermionic systems in the PT-symmetric phase, the PTsymmetric generalized Gibbs ensemble (PTGGE), using two different approaches: First, we present a derivation based on the dephasing of the covariance matrix, which has the benefit of being more practicable for analytical and numerical purposes; then, we show how the PTGGE can be derived from the quadratic eigenmodes of the adjoint Liouvillian.This, in comparison, is a more physically insightful approach, offering a better understanding in terms of the evolution of Liouvillian eigenmodes with modified dynamics and statistics, and connects neatly to the structure found in isolated systems.
We have discussed the derivation of the PTGGE for the driven-dissipative Kitaev chain in depth in Ref. [37].Therefore, in Sec.IV A, we present a detailed derivation of the PTGGE for the SSH model with incoherent loss and gain, and we provide only a brief summary of the corresponding results for the Kitaev chain in Sec.IV B.

A. Driven-dissipative SSH model
As stated above, we commence the derivation of the PTGGE for the SSH model in terms of the covariance matrix.This derivation can be divided into three steps: The first step, which we have taken in Sec.III A 4, is to find an analytical result for the evolution of the covariance matrix.Then, we obtain the long-time asymptotic behavior, determined by dephasing of momentum modes.Finally, we relate the dephased covariance matrix to the density matrix of the system.The last two steps are presented in Sec.IV A 1 below.After this we proceed with the derivation by finding quadratic eigenmodes of the Liouvillian in Sec.IV A 2.

Derivation of the PTGGE from dephasing of the covariance matrix
In isolated integrable models that can be mapped to noninteracting fermions, generalized thermalization to a maximum entropy ensemble following a quantum quench happens through dephasing of momentum modes that oscillate at different frequencies 16,67].Consequently, local observables take on stationary expectation values predicted by the generalized Gibbs ensemble.In drivendissipative integrable systems, the same mechanism is responsible for relaxation to a maximum entropy ensemble inside the PT-symmetric phase [37], yet there are fundamental differences which we will illustrate in the following.To that end, let us study a general block g l (t) of elements of the covariance matrix, determined by the momentum-space representation g k (t) given in Eq. ( 48) through inversion of Eq. ( 45), For the moment, let us consider balanced loss and gain δ = ∆ = 0 such that g 2,k (t) = 0; explicit expressions for g 1,k (t) are given in Eqs. ( 54) and (55).Note that g 1,k (t) depends on time through an overall factor e −4γt , and through oscillating factors sin(2ω k t) and cos(2ω k t).According to the Riemann-Lebesgue lemma, to obtain the behavior of g l (t) for t → ∞, we have to drop these oscillating terms [16].Stated in terms of g 1,k (t), in the limit t → ∞, we may write where the subscript "d" denotes the dephased value obtained by omitting oscillatory contributions, and where g ′ d,1,k is timeindependent and explicitly given by where and with the angle ∆ϕ k defined by The asymptotic behavior of g l (t) is thus given by Since for a Gaussian state any expectation value can be expressed in terms of the covariance matrix by using Wick's theorem, this result shows that after appropriate rescaling to compensate overall exponential decay, local observables relax to stationary values that are determined by the dephased covariance matrix.In particular, the expectation value of a product O ℓ of ℓ fermionic operators becomes stationary after rescaling with a factor of e 2ℓγt .We note that this applies only when ⟨O ℓ ⟩ SS = tr(O ℓ ρ SS ) = tr(O ℓ )/2 2L = 0. Otherwise, one should consider the tracless part O ℓ − tr(O ℓ )/2 2L , the expectation value of which measures actual correlations and vanishes in the steady state at infinite temperature.Further, similarly to isolated systems, in the equilibrated values of rescaled local observables, memory of the initial state is preserved through the angle defined in Eq. ( 74).
As a consequence of Gaussianity, not only arbitrary expectation values but also the full state of the system is determined by the covariance matrix.In particular, the density matrix that describes the PTGGE corresponding to the dephased covariance matrix, is given by [88] ρ where Z(t) is a normalization such that tr(ρ PTGGE (t)) = 1.We note that while our derivation of ρ PTGGE (t) is based on explicit results for the covariance matrix for the drivendissipative SSH model, our considerations clearly generalize to other open fermionic systems that are described by number-conserving quadratic fermionic Hamiltonians and are subjected to incoherent loss and gain.Even though the same mechanism of dephasing underlies both relaxation of isolated systems to the GGE and of PTsymmetric systems to the PTGGE, there are several important differences: (i) The PTGGE is intrinsically time-dependent due to the exponential decay of correlations at a rate determined by the imaginary part of the eigenvalues λ ±,k of z k .Crucially, in the PT-symmetric phase, the decay rate is identical for all momentum modes.Relaxation of local observables to the PTGGE is then visualized best by factoring out the overall exponential decay.(ii) The oscillation frequencies in driven-dissipative models are given by the Liouvillian dispersion ω k and not by the bare Hamiltonian dispersion relation ε k .This affects the characteristic time scale of relaxation to the PTGGE.(iii) Until now, we have considered balanced loss and gain rates with δ = ∆ = 0.For finite values of δ and ∆, also g 2,k (t) 0 in Eq. ( 69) is nonzero, and this contribution to the covariance matrix has no counterpart in isolated systems.While in driven-dissipative systems information about the initial state is incorporated in g 1,k (t), the contribution g 2,k (t) describes the approach to the steady state, which is nontrivial for δ, ∆ 0. Hence, for finite δ and ∆, the state of the system will deviate from the PTGGE for t → ∞.However, as we explain in the following, for sufficiently small values of δ and ∆, one can clearly observe transient relaxation to the PTGGE.
Similarly to g 1,k (t), the contribution g 2,k (t) dephases in the long-time limit and contains terms that decay exponentially.However, as mentioned before, the integral in Eq. ( 50) also includes a time-independent steady-state contribution, where and with the steady state contribution given by Hence, for δ, ∆ 0, Eq. ( 75) is replaced by Clearly, at some point the steady state contribution will dominate over the exponentially decaying terms, and the covariance matrix will assume its steady state form g l (t) ∼ g SS,l .Yet, for sufficiently small values of δ and ∆, the PTGGE still gives an accurate description for relaxation of local observables on intermediate time scales, up to the crossover time t × at which g SS,l becomes dominant.The explicit value of t × depends on the observable under consideration.For the example of correlations within a unit cell as measured by ⟨c B,l c † A,l (t)⟩ = G B,A l,l (t)/2 = g 2,1 0 (t)/2, we can estimate the crossover time as follows: Assuming that dephasing happens on a time scale that is shorter than t × , the crossover time t × is determined by the condition that the absolute value of the exponentially decaying term in Eq. ( 81) is equal to the steady state contribution, which leads to Since g SS,l is proportional to y k given in Eq. ( 47), this estimate implies that t × diverges logarithmically for δ, ∆ → 0. The time evolution of ⟨c B,l c † A,l (t)⟩ and the logarithmic divergence of t × are illustrated in Fig. 3.
Finally, we want to contrast relaxation to the PTGGE with the asymptotic behavior of the covariance matrix in the PTmixed and PT-breaking phases, again based on the dynamics of ⟨c B,l c † A,l (t)⟩, which, for δ = ∆ = 0, is given by In the PT-breaking phase, the dispersion is purely imaginary, ω k = iκ k with κ k ∈ R >0 given in Eq. (34).After some algebraic simplifications of g 1,k (t), the above integral reads [37] ⟨c Clearly, the dominant contribution is due to momenta in the vicinity of the maximum of κ k .In the PT-mixed phase, the integration in Eq. ( 83) can be split into a contribution only consisting of PT-symmetric modes and another contribution due to PT-breaking modes, and again the dominant contribution at late times is due to PT-breaking modes in the vicinity of the maximum of κ k .For example, for |∆J| > J, the dispersion κ k has a maximum at k max = 0. Using standard asymptotic expansion techniques [89], we obtain That is, in the PT-mixed and PT-breaking phases, there is a single momentum mode k max that maximizes κ k and dominates the dynamics, and the continuum of modes in the vicinity of k max leads to additional algebraic decay.In contrast, all momenta k ∈ BZ contribute to the result for the PT-symmetric phase given in Eq. ( 75).

Derivation of the PTGGE from the principle of maximum entropy
As we show next, the PTGGE can also be derived from the principle of maximum entropy [10] by properly taking into account the modified statistics of Liouvillian eigenmodes and the expectation values of their commutators.This approach is explained best by considering first the isolated SSH model.The dynamics of the isolated SSH model are generated by the Hamiltonian superoperator H, whose action is defined in Eq. ( 2), and with eigenmodes d s,k given in Eq. ( 29).Since we are concerned with Gaussian states that are by definition fully determined by two-point functions, let us consider also quadratic forms of operators.Any quadratic form of eigenmodes d s,k can be expressed in terms of commutators and anticommutators through the decomposition For fermions, statistics are encoded in anticommutators, 87) By contrast, commutators describe the dynamics.In particular, commutators of the modes d s,k are also eigenmodes of the Hamiltonian superoperator H: That is, mode-diagonal commutators are conserved, while mixed-index commutators oscillate with frequency 2ε k and dephase.In quadratic fermionic models or in integrable models that can be mapped to noninteracting fermions, the GGE is usually stated as the maximum entropy ensemble that is compatible with conserved mode occupation numbers n s,k = d † s,k d s,k [16,17].According to our discussion, we can equivalently define the GGE as the maximum entropy ensemble that is consistent with canonical anticommutations relations Eq. ( 87) and the conservation of mode-diagonal commutators Eq. (88).Crucially, this latter definition generalizes to the PTGGE, however, with some important differences: First, the adjoint Liouvillian L † that generates the dynamics of operators in the driven-dissipative setting, and is specified below for the SSH model, is a non-Hermitian operator.Therefore, its eigenmodes obey modified noncanonical anticommutation relations.Second, mode-offdiagonal commutators of eigenmodes of L † oscillate at modified frequencies 2ω k , directly affecting the dynamics.And third, mode-diagonal commutators are not conserved; instead, these commutators decay exponentially.But crucially, they do not oscillate and are, therefore, not affected by dephasing.Based on the above definition of the PTGGE, in the following, we present a detailed derivation of the PTGGE for the driven-dissipative SSH model.The computation consists of three steps: (i) Specifying the generator of operator dynamics L † and the corresponding eigenvalue equation.(ii) Solving the eigenvalue equation to obtain (ii.a) nonoscillatory, mode-diagonal and (ii.b) oscillatory, mode-offdiagonal commutators of eigenmodes of L † .(iii) Constructing the PTGGE as the maximum entropy ensemble that is compatible with the statistics of the eigenmodes of the adjoint Liouvillian L † and the expectation values of nonoscillatory commutators.
(i) Adjoint Liouvillian.As detailed in Appendix A, for a density matrix ρ evolving according to a Liouvillian superoperator L, the expectation value of an operator O follows the equation of motion given by with the adjoint Liouvillian Hermitian conjugation is defined here with respect to the Hilbert-Schmidt scalar product, leading to H = H † , and According to Eq. ( 89), the adjoint Liouvillian L † generates the dynamics of operator expectation values.Suppose now that O is an eigenmode of L † in the following sense: where O SS is a number.Then, the equation of motion of the expectation value ⟨O⟩ reduces to Dynamical stability requires the imaginary part of the eigenvalue λ to be negative, such that the expectation value ⟨O⟩ approaches O SS for t → ∞.
(ii) Eigenmodes of the adjoint Liouvillian.We want to solve the eigenvalue equation ( 92), where we consider quadratic eigenmodes of the adjoint Liouvillian.In particular, we seek eigenmodes in the form of mode-diagonal and mode-offdiagonal commutators, We note that only the commutators η k and χ k and not the modes d s,k themselves satisfy the eigenvalue equation (92).Nevertheless, for simplicity, we refer to both the modes d s,k and the commutators η k and χ k as eigenmodes of L † .In analogy to Eq. ( 88), we anticipate that the expectation values of the diagonal commutators η k are nonoscillatory whereas the expectation values of the offdiagonal commutators χ k are oscillatory and, therefore, subject to dephasing.
(ii.a) Nonoscillatory eigenmodes.To find the nonoscillatory eigenmodes of the adjoint Liouvillian, we use a general bilinear ansatz given by where the goal is to find Q such that η satisfies the eigenvalue equation ( 92).This approach does not rely on translational invariance and, therefore, we omit the momentum index for the time being.Plugging the ansatz into the eigenvalue equation (92), we obtain two contributions on the left-hand side that are due to the Hamiltonian superoperator H and the adjoint dissipator D † , respectively.The Hamiltonian part reads and the action of the dissipator is given by We write the eigenvalue λ in Eq. ( 92) as λ = −iκ with an undetermined real parameter κ ∈ R, and we identify the steadystate value η SS with the last term in Eq. ( 97), Then, the eigenvalue equation takes the form Multiplying this equation with c m c † m ′ and using we obtain two equations for Q and its transpose Q ⊺ , respectively, where Z is defined in Eq. ( 14).Since Z = Z ⊺ , the two equations above are actually equivalent.To make further progress, we use translational invariance of the driven-dissipative SSH model, which implies that the eigenmodes η k are labeled by a momentum k, and that Q k are block Toeplitz matrices with 2 × 2 blocks such that the commutator of Liouvillian eigenmodes in Eq. ( 94) now takes the form with Defining the discrete Fourier transformations of the blocks q k,l and the spinors C l as we can recast Eq. ( 102) as and write the commutator in momentum space, In the PT-symmetric phase, z k can be diagonalized as stated in Eq. (C4), and the spectrum of z k is given by σ(z k ) = {−i2γ ± ω k } where ω k ∈ R >0 .Using this representation of z k in Eq. ( 106), we obtain with the shorthand notation Then, identifying κ = 4γ, we can rewrite this equation as a commutator, Therefore, for η k to be an eigenmode of the adjoint Liouvillian with eigenvalue λ = −i4γ, q k,k ′ has to satisfy the above commutation relation.The general solution for q k,k ′ reads where α 1,k,k ′ and α z,k,k ′ are undetermined parameters.To obtain the solution anticipated in Eq. ( 94), we choose Then, Eq. ( 107) takes the form where Defining now the Liouvillian eigenmodes as with we obtain the final form of η k : The angle ϕ k is defined in Eq. (C2) and ψ k is determined by the relation In the transformation to the Liouvillian eigenmodes in Eq. ( 114), we have chosen the normalization such that {d A,k , d † A,k } = 1.This leads to with θ k given in Eq. (C3) and related to ψ k by Note that the alternative choice Further, the transformation V k given in Eq. ( 114) determines the statistics of the Liouvillian eigenmodes.Since V k is not unitary, the operators d s,k do not obey the usual canonical anticommutation relations.The statistics of these modes are encoded in the anticommutators collected in where Next, having specified the solution for the matrix Q that leads to eigenmodes η k in the form given in Eq. ( 115), we can determine the steady state contribution in Eq. ( 98).We find Finally, to fully determine the time evolution of the expectation value ⟨η k ⟩, we have to calculate its initial value.Relating the expectation value of the commutator to the covariance matrix by we find the initial value, determined by the covariance matrix g k (0) in Eq. ( 51), to be given by where the angle ∆ϕ k is defined in Eq. (74).Therefore, by solving Eq. ( 93) for the commutator in Eq. (115), we obtain the time evolution of the expectation value of diagonal commutators: with the initial and steady-state values given in Eq. ( 123) and Eq. ( 121), respectively.
(ii.b) Oscillatory eigenmodes.We proceed by showing that the offdiagonal commutators ] with mode operators d s,k given in Eq. ( 113) are oscillatory eigenmodes of the Liouvillian.To that end we write χ k in the form where we define and σ ± = (σ x ± iσ y )/2.The fact that the offdiagonal commutators χ k are eigenmodes of L † with eigenvalue λ = −2 (ω k + i2γ) can be confirmed by noting that p k,k ′ satisfies the eigenvalue equation that is similar to Eq. ( 106): The time evolution of expectation values of χ k thus reads where the initial and steady-state values are given by (iii) Construction of the PTGGE.We consider now the case of balanced loss and gain with δ = ∆ = 0, such that the steady state contributions to Eqs. ( 124) and ( 128) vanish, Further, due to the oscillatory exponent in ( 128), the contribution of offdiagonal commutators to local observables vanishes at late times due to dephasing between offdiagonal commutators with different momenta Therefore, the state of the system at late times is determined by the expectation values of diagonal commutators given in Eq. ( 124) and the statistics of the modes d s,k that follow from the nonunitary transformation V k Eq. (114) and are encoded in the matrix f k in Eq. ( 120).The PTGGE can be defined as the maximum entropy ensemble that is compatible with these requirements.We can also collect the same information in the dephased covariance given in Eq. ( 76), which can be related to the eigenmodes by where Then for the given two-point function g PTGGE,k (t), the entropy is maximized for the Gaussian state that is uniquely determined through Eq. ( 77).In terms of the eigenmodes of the adjoint Liouvillian, the PTGGE can be written as which shows most transparently both the dependence of the PTGGE on the initial conditions through ζ k (t) and the effect of the noncanonical anticommutation relations of the Liouvillian eigenmodes described by f k .

B. Driven-dissipative Kitaev chain
Having presented the derivation of the PTGGE for the SSH model, we now turn to the driven-dissipative Kitaev chain.Since the latter has been discussed in detail in Ref. [37], we restrict ourselves here to highlighting the aspects in which the derivation for the Kitaev chain differs from the one for the SSH model.In particular, for the Kitaev chain, we define bilinear forms of eigenmodes as which can be regarded as normal and anomalous commutators, respectively, and where the Liouvillian eigenmodes d k are related to the original fermionic operators c k = with angles θ k and ϕ k that are defined through the relations The Liouvillian eigenmodes d k obey noncanonical anticommutation relations, where In analogy to the modediagonal commutators for the SSH model, the normal commutators η k are nonoscillatory, where k ]⟩ SS contains the steady state contribution.The nonoscillatory normal commutators decay with an overall decay rate γ defined in Eq. ( 62), and are not affected by dephasing.In contrast, and analogously to the mode-offdiagonal commutators of the SSH model, the anomalous commutators χ k are oscillatory, The anomalous commutators oscillate with frequency ω k and are thus affected by dephasing.As above, the PTGGE describes the late-time relaxation dynamics for vanishing steady-state contributions accomplished by balanced loss and gain rates δ = γ l − γ g = 0, and after dephasing of the contributions due to anomalous commutators.The explicit form of the PTGGE is given by where The key difference between Eqs. ( 141) and ( 134) for the Kitaev chain and the SSH model, respectively, is that in general anomalous correlations such as ⟨c l c l ′ ⟩ do not vanish for the Kitaev chain; in contrast, the weak U(1) symmetry of the driven-dissipative SSH model ensures that anomalous correlations like ⟨c s,l c s ′ ,l ′ ⟩ vanish at all times.

V. SPREADING OF CORRELATIONS
Now that we have derived the ensemble that describes the late-time dynamics for quenches to the PT-symmetric phase, we turn to a detailed study of relaxation to the PTGGE.In the following sections, both for the SSH model and the Kitaev chain, we focus on balanced loss and gain with δ = ∆ = 0, such that the steady state is at infinite temperature, and a description of the dynamics in terms of a PTGGE applies on arbitrarily long time scales.
A well-established property of isolated integrable systems that exhibit generalized thermalization after a quantum quench is the ballistic propagation of correlations through the system [16,18].This is also referred to as light cone spreading of correlations, since correlations outside of the light cone that is defined by the group velocity v max are suppressed exponentially [90].As we demonstrate in the following, light cone propagation of correlations is not unique to isolated systems [91,92].We illustrate this behavior for the driven-dissipative Kitaev chain, where we consider the evolution of the normal commutators C l−l ′ (t) = ⟨[c l , c † l ′ ](t)⟩ for δ = ∆ = 0.In the PT-symmetric phase, we find ballistic propagation of quasiparticles which, however, have a finite lifetime, and a maximum velocity adjusted to the modified dispersion relation ω k ; in contrast, in the PT-mixed and PT-broken phases, we observe diffusive spreading of correlations.
We conclude that the spreading of correlations with finite velocity is maintained also in the driven-dissipative model throughout the whole PT-symmetric phase.For a finite rate of dissipation γ, the light cone velocity with dispersion ω k is increased compared to the velocity corresponding to the isolated system with dispersion ε k .This is further verified in Fig. 4(c), where the position of the light cone boundary is traced numerically and compared to the analytical prediction.However, the increased speed at which correlations propagate is offset by the exponential decay of correlations at rate 4γ.Further, the phase velocity defined as [93] v PH = |ω k /k| k=k max , where k max is the momentum for which v k = v max , decreases with increasing reservoir coupling.While v max defines the boundary beyond which correlators become exponentially suppressed, the phase velocity describes the propagation of local peaks within the light cone.At the boundaries of the PTsymmetric phase, the phase velocity and maximum particle velocity line up with v max = v PH = 2J |µ| for µ < 0, while for µ > 0 we still have v max = 2J |µ| but the phase velocity v PH = 0.This behavior is depicted in Fig. 5, where v max and v PH are plotted as a function of µ.For the value µ = −0.5Jchosen in Figs.4(a) and (b), we obtain v PH > v max .However, as can be seen in Fig. 5, also the opposite order v PH < v max is realized for different values of µ.

B. PT-mixed and PT-broken phases
In the PT-mixed and PT-broken phases, where γ > γ c = |J − |µ| /2|, the spreading of correlations is dominated by the single slowest-decaying mode with decay rate Accordingly, we define rescaled normal commutators as C ′ l (t) = e 4γ s t C l (t).In Fig. 6(a) and (b), the absolute value of rescaled normal commutators is shown for the PT-mixed and PT-broken phase, respectively.It is worthwhile to first discuss the qualitative differences of correlations in the phases with PT-symmetric and PT-breaking modes.In the PT-symmetric phase, at any time, the normal commutators show multiple oscillations inside the light cone, with a peak at the boundary of the light cone and ensuing suppression of correlations outside the light cone.This creates the sharp boundaries in Fig. 4. In contrast, in the presence of PT-breaking modes, normal commutators show a single peak without oscillations and long decaying tails, giving the unstructured appearance inside the boundary in Fig. 6.Crucially, a correlation boundary can still be defined through this single peak.The correlation boundary, however, spreads diffusively according to ℓ = 2(Dt) 1/2 , where the diffusion constant is given by with κ k defined in Eq. ( 34) and k max the momentum maximizing κ k .The diffusive evolution of the peak position is illustrated further in Fig. 6(c).
Let us briefly comment on the spreading of correlations in the driven-dissipative SSH model.Based on the analytical form of the covariance matrix derived in Sec.III A 4, which is structurally similar to the covariance matrix of the Kitaev chain, we can expect that the dynamics of correlations is qualitatively the same in both models.This expectation is confirmed by numerical results, which we do not show here.

VI. TIME EVOLUTION OF THE SUBSYSTEM ENTROPY
The linear growth and volume-law saturation of the von Neumann entropy of a finite subsystem is a key signature of thermalization in isolated systems [19][20][21][22].As we show in the following, the exact same phenomenology can be observed after quenches to the PT-symmetric phase of the drivendissipative SSH model-if we consider the contribution to the entropy that measures the spreading of correlations due to the propagation of pairs of entangled quasiparticles, and after appropriate rescaling to compensate for exponential decay.The corresponding analysis for the driven-dissipative Kitaev chain is provided in Ref. [37].
We consider a subsystem of the SSH chain that contains ℓ contiguous unit cells.The corresponding reduced density matrix is obtained by tracing out the remaining L−ℓ unit cells, ρ ℓ = tr L−ℓ (ρ), and the von Neumann subsystem entropy is defined by S vN,ℓ = − tr(ρ ℓ ln(ρ ℓ )).For a Gaussian state, the subsystem entropy can be calculated as [88] where and where ±ξ l with 0 ≤ ξ l ≤ 1 and l ∈ {1, . . ., ℓ} are the eigenvalues of the reduced covariance matrix G ℓ defined by with G s,s ′ l,l ′ given in Eq. ( 6).The set of eigenvalues {±ξ l } of G ℓ forms the single-particle entanglement spectrum [88].In pure states of isolated systems, the subsystem entropy measures the entanglement between the subsystem and its compliment.After a quench, the subsystem entropy grows linearly in t before it saturates to a volume-law value, S vN,ℓ ∝ ℓ.This behavior can be explained through a quasiparticle picture, which does not only provide a qualitative interpretation but also quantitative predictions for the full time evolution after a quench in the space-time scaling limit [19][20][21][22].In this picture, the initial state acts as a source of pairs of entangled quasiparticles with opposite momenta.After creation, these quasiparticles move ballistically with a velocity of at most v max through the system.All pairs whose members are separated by the boundary between the subsystem and its compliment contribute to the subsystem entropy.Consequently, the subsystem entropy starts to saturate when all maximum-velocity pairs generated in the subsystem have left the subsystem.
In open systems, the time-evolved state is no longer pure and the subsystem entropy involves two contributions [38][39][40]94]: S QP vN,ℓ , which describes correlations due to quasiparticle pairs as in an isolated system; and S stat vN = S vN,L , the statistical entropy due to the mixedness of the state.The quasiparticlepair contribution is thus given by the difference In the following, we consider quenches to the PT-symmetric phase of the driven-dissipative SSH model.Building upon findings of Refs.[38,39], in Ref. [37], we have proposed an analytical conjecture for the evolution of the quasiparticle-pair contribution S QP vN,ℓ in the space-time scaling limit: where for the SSH model ζ k (t) and f k are defined in Eqs.(133) and (120), respectively, and g k (t) is given in Eq. ( 48) with the subscript "d," which stands for dephasing, indicating that only nonoscillatory components contribute.For long times, γt ≫ 1, we can expand the quasiparticle-pair entropy in the exponentially decaying factors ζ k (t), g k (t) ∼ e −4γt .To lowest nontrivial order, we find S (ξ) ∼ ln(2) − ξ 2 /2.Then, the constant contributions in Eq. ( 151) cancel, and we find S QP vN,ℓ (t) ∼ e −8γt .Therefore, relaxation to the PTGGE can be revealed by considering the rescaled quasiparticle-pair entropy e 8γt S QP vN,ℓ (t), which is shown in Fig. 7.The rescaled quasiparticle-pair entropy grows linearly up to the Fermi time [14], where saturation to the value predicted by the PTGGE sets in.As shown in the inset, the difference between the numerical results and the analytical conjecture Eq. ( 150) vanishes as 1/ℓ.The numerical method we have used to obtain this data is described in Ref. [37].

VII. TIME EVOLUTION OF STRING ORDER AND SUBSYSTEM FERMION PARITY
Generalized thermalization after a quantum quench in isolated systems is defined as relaxation of the averages of local observables to stationary values that are determined by the  151) for the space-time scaling limit.The approach of the numerical results to the conjecture with increasing subsystem size is shown in the inset: For the topological quench at t = 2t F , the difference between the numerical data and Eq. ( 151) (blue dots) vanishes as 1/ℓ (orange line).
GGE.Similarly, based on the arguments given in Sec.IV, after quenches to the PT-symmetric phases of the driven-dissipative Kitaev and SSH chains, we expect to observe local relaxation to the PTGGE.To confirm this expectation, we proceed to study the time evolution of the dual string order parameter for the SSH model and the subsystem parity for the Kitaev chain, which are defined below in Secs.VII A and VII B, respectively.These observables are of special interest due to their connection to topology: On the one hand, in the ground state of an isolated system, the dual string order parameter and the subsystem parity function as topological disorder parameters.That is, they are finite in the trivial phase and vanish in the topological phase.On the other hand, after quenches originating from the trivial phase, the decay of dual string order and subsystem parity carries a robust signature of the topology of the postquench Hamiltonian: For quenches across the topological phase boundary, the topological disorder parameters exhibit periodically recurring zero crossings; in contrast, they do not cross zero for quenches within the trivial phase.Physically, these zero crossings correspond to pumping of spin order and fermion parity between a finite subsystem and its complement.As we discuss in the following, in the presence of drive and dissipation, this phenomenology acquires qualitative modifications that are unique to open systems.

A. Dual string order parameter
The concept of string order has originally been introduced and is mostly discussed in the context of spin chains [95,96].For the SSH model Eq. ( 4), an equivalent formulation in terms of a dimerized spin-1 /2 chain can be obtained through the Jordan-Wigner transformation [97], which yields where S µ l with µ ∈ {x, y, z} are spin-1 /2 operators that act on the spin at site l with l ∈ {1, . . ., 2L}.The Jordan-Wigner transformation, which maps the set of 2L fermionic operators c l and c † l with c A,l = c 2l−1 and c B,l = c 2l to spin operators S ± l = S x l ± iS y l , is given by where n l = c † l c l .For ∆J < 0, the ground state of the dimerized spin chain is topologically ordered.Indeed, the topologically ordered phase of the dimerized spin-1 /2 chain is continuously connected to the Haldane phase of the antiferromagnetic spin-1 Heisenberg chain [98].The topological and trivial phases of the dimerized spin chain can be distinguished through the string order parameter, which takes a finite expectation value in the topological phase, and vanishes in the trivial phase.Equivalently, the dual string order parameter is nonzero in the trivial phase, and vanishes in the topological phase [70].The designation as dual refers here to the self duality of the dimerized spin chain: Under a translation by one lattice site, S µ l → S µ l+1 , the Hamiltonian Eq. ( 153) maps to itself with ∆J → −∆J.Consequently, also the trivial and topological phases are exchanged, and, therefore, the dual string order parameter, which is obtained by performing the translation by one lattice site on the usual string order parameter, serves as a topological disorder parameter.The dual string order parameter for the dimerized spin-1 /2 chain is given by [70,99] with µ ∈ {x, y}.We note that the invariance of the Hamiltonian Eq. ( 153) under rotations around the z-axis implies that ⟨O x ℓ ⟩ = ⟨O y ℓ ⟩.An alternative interpretation of string order is obtained through a nonlocal unitary transformation that maps the dimerized spin chain Eq. ( 153) to two coupled Ising models [100], which have the doubled Ising symmetry Z 2 × Z 2 .This symmetry is broken in the topologically ordered phase, and correlations of the Ising order parameters map to the string order parameters O x ℓ and O y ℓ .To compute the dual string order parameter for the SSH model, it is convenient to introduce 4L Majorana operators w l according to Eq. (63).For concreteness, we consider the dual string order parameter with µ = x, which can be written in terms of Majorana fermions as Since the quench dynamics we consider here only involve Gaussian states, we can employ Wick's theorem [101][102][103], according to which the expectation value of the dual string order parameter is given by the Pfaffian of a submatrix of the covariance matrix Γ for Majorana fermions defined in Eq. ( 64), The submatrix Γ ℓ consists of the elements of Γ that correspond to the Majorana operators appearing in Eq. ( 156), As shown in Appendix D, up to a factor of i, Γ ℓ is unitarily equivalent to the contribution G 1,ℓ to the reduced covariance matrix of complex fermions defined in Eq. ( 149), where The decomposition of the reduced covariance matrix as G ℓ = G 1,ℓ +G 2,ℓ employed here is analogous to that of g k in Eq. ( 48).

Relaxation of the dual string order parameter
To illustrate relaxation to the PTGGE, we consider the time evolution of the dual string order parameter after quenches to the trivial and topological PT-symmetric phases.As discussed in Sec.IV A 1, the expectation value of a product of ℓ fermionic operators becomes stationary after rescaling with a factor of e 2ℓγt .Accordingly, in Fig. 8, we show the evolution of the rescaled dual string order parameter e 4ℓγt ⟨O x ℓ ⟩ .For quenches to both the topological (blue) and the trivial (green) phase, we observe fast decay up to the Fermi time t F Eq. ( 152), followed by much slower relaxation to the value predicted by the PTGGE.As an aside, it is interesting to note that, since G 2,ℓ does not contribute to Γ ℓ in Eq. ( 159), relaxation of the dual string order parameter is described by the PTGGE even if δ, ∆ 0. The numerical data for the quench to the trivial phase is in excellent agreement with the following analytical form: where ζ ′ k and f k are defined in Eqs. ( 133) and ( 120), respectively, and the prefactor O 0 is determined by fitting the late-time asymptotic behavior to the PTGGE prediction.Equation ( 161) is an analytical conjectures we have proposed originally for the dynamics of the subsystem parity in the drivendissipative Kitaev chain [37], and which is based on analytical results obtained by Calabrese et al. [13][14][15] for the relaxation of order parameter correlations in the transverse field Ising model and in the space-time scaling limit ℓ, t → ∞ with ℓ/t fixed.As mentioned above, the SSH model can be mapped to two coupled Ising models by a nonlocal unitary transformation [100], which leads us to expect that our conjecture also applies to the dual string order parameter.Numerically, this expectation is confirmed in Fig. 8, where the solid lines indicate the analytical predictions.
The analytical conjecture Eq. ( 161) shows that as compared to the isolated SSH model, there are two important modifications: (i) The dynamics are determined the Liouvillian dispersion relation ω k rather than the Hamiltonian dispersion relation ε k , leading here to a shorter Fermi time t F ; (ii) noncanonical anticommuation relations of Liouvillian quasiparticles encoded in f k 1 affect the result quantitatively.
For quenches to the topologically nontrivial phase with ∆J < 0, the dual string order parameter exhibits additional oscillatory dynamics for t < t F , well described in the spacetime scaling limit by where ⟨O x ℓ (t)⟩ nonosc is the nonoscillatory evolution described by Eq. ( 161), α ± are numerically determined phase shifts, and k s,± are the soft modes of the PTGGE.These modes are determined by the condition that the matrix in the exponent in Eq. ( 77) has an eigenvalue that is equal to zero.According to Eq. ( 132), this condition can be stated as or, equivalently, cos(∆ϕ k ± ψ k ) = 0.As we show below, cos(∆ϕ k − ψ k ) = 0 has two solutions k = k s,± .The solutions to cos(∆ψ k + ψ k ) = 0 are then given by −k s,± : Indeed, IS of the Bloch Hamiltonian h k and IS † of the non-Hermitian Bloch Hamiltonian z k imply that ∆ϕ k = −∆ϕ k and ψ k = ψ −k , respectively [37], and, therefore, cos For our choice of initial state with vanishing intercell hopping, J 2,0 = 0, this condition leads to The two solutions are given by We have chosen the designation of each of the solutions of Eq. ( 164) as either k s,+ or k s,− such that k s,− and k s,+ are associated with the pumping of string order through the right and left boundary of a subsystem, respectively, as detailed below.Note that the designation is reversed when J 1 = J + ∆J changes sign.Crucially, for this choice, the frequencies ω k s,± are continuous functions of ∆J and ∆γ.For ∆γ 0, the PTGGE has two distinct soft modes k s,+ −k s,− .In contrast, for the GGE, which is obtained for ∆γ → 0, the soft modes are locked onto each other by inversion symmetry, k s,+ = −k s,− .Associated with the soft modes are two time scales [104], which determine the oscillation periods of the dual string order parameter, and, in particular, the periodicity of the recurring zero crossings of the dual string order parameter.Note that since ω k = ω −k , there is only a single time scale t s = t s,+ = t s,− for isolated systems.As shown in Fig. 8, the periodicity of zero crossings persists up to t F , and the subsequent relaxation of the dual string order parameter for t > t F approximately follows the prediction given in Eq. ( 161).
The oscillatory dynamics of the dual string order parameter for quenches to the topological phase are connected to the topological zero crossings in the entanglement spectrum.We explore this connection in more detail in the following section.

Topological zero crossings
For quantum quenches in the isolated SSH model, zero crossings of the dual string order parameter are a robust signature of the topology of the postquench Hamiltonian.To explain why that is the case, we observe first that according to Eq. (157) zero crossings of the dual string order parameter occur whenever an eigenvalue of Γ ℓ , or, equivalently as per Eq. ( 159), of the reduced covariance matrix G ℓ , crosses zero.The coincidence of zero crossings of the dual string order parameter and in the single-particle entanglement spectrum applies to all Gaussian states, and is illustrated for a quench to the topological PT-symmetric phase in Fig. 9. Furthermore, the dynamical entanglement spectrum bulk-boundary correspondence for one-dimensional lattice models [4] that belong to the class BDI [105] states that for quenches originating from the trivial phase, zero crossings in the entanglement spectrum occur if and only if the postquench Hamiltonian is topologically nontrivial.Consequently, also zero crossings of the dual string order parameter occur if and only if the Figure 10.Directional string order pumping for a quench to the topological PT-symmetric phase with ∆J = −0.5J,γ = ∆γ = 0.2J, δ = ∆ = 0, and ℓ = 40.For PBC, the numerical data (black line) crosses zero at multiples of both t s,+ (light green) and t s,− (dark green), with blue shading indicating the sign of the dual string order parameter.The analytical conjecture Eq. ( 162) agrees well with the numerics (red line).For OBC, zero crossings occur at multiples of t s,− or t s,+ depending on whether the subsystem is at the left (L ℓ , violet line) or the right end of the chain (R ℓ , blue line).Additional rescaling with e ±2γt of the data for L ℓ and R ℓ compensates for exponential decay and growth, respectively, due to edge modes.
postquench Hamiltonian of the isolated SSH model is topologically nontrivial.As we discuss in more detail below, for the Kitaev chain, the role of the dual string order parameter is played by the subsystem fermion parity [37].
Turning now to the driven-dissipative SSH model, the continuity of non-Hermitian real-line-gap topology suggests that the relation between zero crossings of the dual string order parameter and nontrivial topology of the generator of the postquench dynamics remains valid through the entire PTsymmetric phase.An important qualitative modification due to drive and dissipation has already been mentioned above: There are now two distinct time scales for zero crossings of the dual string order parameter, which are determined by soft modes of the PTGGE.In Ref. [37], for the driven-dissipative Kitaev chain, we have provided a physical interpretation of these two time scales in terms of directional parity pumping: The two time scales correspond to different rates of the exchange of parity between a subsystem and its complement through the left and right ends of the subsystem.This effect requires the breaking of inversion symmetry and mixedness of the time-evolved state, and is thus unique to driven-dissipative systems.As we illustrate in Fig. 10, the driven-dissipative SSH model exhibits an analogous effect of directional string order pumping.The figure shows the dynamics of the rescaled dual string order parameter for a quench to the PT-symmetric topological phase.After a short initial period, the numerical data for PBC (black line) are well described by the analytical prediction Eq. (162) (red line), and exhibit zero crossings at multiples of both t s,− and t s,+ .In contrast, the dual string order parameter for subsystems L ℓ = {1, . . ., ℓ} at the left end and R ℓ = {L − ℓ + 1, . . ., L} at the right end of a chain with OBC, crosses zero at multiples of only t s,− and t s,+ , respectively.These observations confirm that string order pumping does indeed occur at different rates through the left and right ends of a subsystem.We note that a previous study of topological entanglement spectrum crossings in the driven-dissipative Kitaev chain has focused on subsystems of type L ℓ and has, therefore, not observed directional pumping [85].
As explained above, our expectation for the occurrence of zero crossings of the dual string order parameter in the topological PT-symmetric phase is based on the continuity of topological properties throughout this phase.Indeed, for the driven-dissipative SSH model, two soft modes k s,± of the PTGGE given by Eq. ( 165) exist within the entire topological PT-symmetric phase.The corresponding frequencies ω k s,± are shown in Fig. 11.Upon increasing ∆γ, one of the soft mode frequencies, ω s,− , vanishes at the transition from the PT-symmetric to the PT-mixed phase.In contrast, ω k s,+ remains finite within the entire PT-mixed phase, and vanishes at the boundaries separating the PT-mixed from the trivial PTsymmetric and the PT-broken phases.This observation suggests that zero crossings of the dual string order parameter can also occur in the PT-mixed phase-an expectation that is confirmed in Fig. 12.The numerical analysis of zero crossings of the dual string order parameter for quenches to the PTmixed phase is more challenging than for the PT-symmetric phase: In the PT-symmetric phase, through a simple rescaling, the overall exponential decay of the covariance matrix can be removed analytically before evaluating the covariance matrix numerically [37].In contrast, due to the finite bandwidth of decay rates in the PT-mixed phase, exponential decay cannot be fully accounted for beforehand, leading at late times to extremely small numbers below the numerical precision.The data shown in Fig. 12 is rescaled by e (4ℓγ s +c)t , where 2γ s is the decay rate of the slowest-decaying PT-breaking mode given in Eq. ( 145), and we have chosen c = 51 to best present the numerical results.In agreement with our expectations based on the values of ω k s,± shown in Fig. 11, the dual string order parameter for the subsystem L ℓ does not exhibit zero crossings, while for R ℓ there are zero crossings at multiples of t s,+ .Note that the numerical data for R ℓ stops abruptly due to the occurrence of exponentially small numbers in the calculation of the Pfaffian.For PBC we, further observe a transition from a fast initial decay to a much slower relaxation behavior, in analogy to but much smoother than the transition at t F for the quenches to the PT-symmetric phase shown in Fig. 8.
Our findings lead us to speculate that the non-Hermitian topology of a suitably defined restriction of the matrix z k to its Figure 12.Directional string order pumping for quenches to the PTmixed phase with ∆J = −5J, ∆γ = γ = 1.1J, and ℓ = 30, for PBC (black line) and for OBC with subsystems located at the left (L ℓ , purple line) and the right edge of the chain (R ℓ , blue line).The blue shading indicates the sign of the dual string order parameter for the subsystem R ℓ .Zero crossings of the dual string order parameter for PBC and the subsystem R ℓ for OBC fit well to the soft mode time scale t s,+ .We include an additional exponential rescaling e ct with c = 51 to improve the presentation of the data.
PT-symmetric eigenmodes remains intact even upon crossing the transition to the PT-mixed phase, and that a finite gap is not essential for such a construction.In physical terms, our results suggest to define directional pumping phases based on the existence of soft modes.For the driven-dissipative SSH model, as can be seen in Fig. 11, the phase boundaries of directional pumping phases are also phase boundaries in terms of the gap structure and PT symmetry.However, as we discuss for the example of a driven-dissipative Kitaev chain with long-range hopping and pairing in Sec.VIII below, this does not have to be the case.

Dynamical criticality
Having introduced the concept of directional pumping phases, we may ask what kind of critical behavior occurs at the transitions between these phases, and whether this behavior is universal.Usually, dynamical critical phenomena in driven-dissipative systems are associated with symmetry-breaking second order phase transitions in the steady state [24][25][26][27][28][29][30][31][32].Then, critical behavior is induced by the closing of the dissipative gap, given by the decay rate of the slowest-decaying eigenmode of the Liouvillian, at the critical point.In particular, the correlation length diverges as ξ ∼ δ g −ν , where δ g = g − g c is the deviation of the parameter g from its critical value g c , and the divergence of the relaxation time scale τ ∼ δ g −νz , where z is the dynamical critical exponent, leads to the phenomenon of critical slowing down.Contrary to that, directional pumping phase transitions mark a qualitative change in the oscillatory dynamics of the dual string order parameter and are not associated with any changes in the steady state.Indeed, we consider here dynamics which, for δ = ∆ = 0, always lead to a steady state at infinite temperature with a vanishing correlation length.Further, the rate of relaxation to the steady state is finite for any value of γ > 0. Instead, directional pumping phase transitions are characterized by divergences of the soft-mode time scales t s,± , which determine the periodicity of zero crossings of the dual string order parameter.The soft modes k s,± and the associated time scales t s,± are genuinely dynamical quantities and depend on both pre-and postquench parameters.However, as we discuss in the following, the critical exponents that govern the powerlaw behavior of k s,± and t s,± at the phase boundaries do not depend on the specific choice of parameters.These exponents are determined by an effective long-wavelength description, which indicates that their values are indeed universal.At the boundaries of the directional pumping phases, the soft modes k s,± approach values k s,±,c .The frequencies ω k s,± = J 2 |sin(k s,± )| vanish at the phase boundaries, implying that k s,±,c = 0, ±π.In analogy to the correlation length exponent ν and the dynamical exponent z, we define critical exponents ν ′ and z ′ through the scaling behavior for δ g = g−g c,± → 0, leading to a divergence of the soft-mode time scales t s,± ∼ δ g . For the SSH model, we consider the parameters g = ∆J and g = ∆γ.Let us assume that the values ∆J c,± and ∆γ c,± correspond to a particular point on the phase boundary of ω k s,± , where k s,± takes the value k s,±,c .To obtain the behavior of k s,± and ω k s,± in the vicinity of this critical point, we set ∆J = ∆J c,± + δ ∆J and ∆γ = ∆γ c,± + δ ∆γ in Eq. ( 164), and expand in k around k s,±,c .We thus find for both g = ∆J and g = ∆γ.Therefore, the time scales t s,± exhibit a square-root divergence with critical exponents irrespective of the direction from which the phase boundary is approached in the ∆J-∆γ plane.

B. Subsystem fermion parity
As mentioned above, the subsystem fermion parity serves as a topological disorder parameter for the Kitaev chain.This can be seen by noting that a combination of the Jordan-Wigner [97] and Kramers-Wannier [106,107] transformations maps the Kitaev chain to the transverse field Ising model such that the trivial phase of the Kitaev chain corresponds to the ferromagnetically ordered phase of the Ising model.Then, orderparameter correlations of the Ising model over a distance ℓ are equivalent to the fermion parity P ℓ of a subsystem of size ℓ of the Kitaev chain [37], which is defined by For a Gaussian state, the expectation value ⟨P ℓ ⟩ = pf(Γ ℓ ) is given by the Pfaffian of the reduced covariance matrix Γ ℓ defined in Eq. ( 64) [101,102].

Topological zero crossings
In complete analogy to our the discussion for the SSH model in Sec.VII A 2, the subsystem parity can be seen to exhibit zero crossings for quenches of the isolated Kitaev chain from the trivial to the topological phase.As shown in Refs.[37], the connection between zero crossings and topology remains valid throughout the entire PT-symmetric phase of the driven-dissipative Kitaev chain.The periodicities t s,± of these zero crossings are determined by the soft modes k s,± , which, for µ 0 → −∞, are given by the solutions to [37] 2J cos(k) + µ = sgn(sin(k))2γ, and the corresponding frequencies ω k s,± determine the directional pumping phases shown in Fig. 13.As in the case of the SSH model, the frequency ω k s,− , which describes parity pumping through the right end of a subsystem, is nonzero only within the topological PT-symmetric phase; in contrast, the frequency ω k s,+ , corresponding to parity pumping through the left end of a subsystem, remains finite within the entire PT-mixed phase.The phase boundaries defined by directional parity pumping coincide with those that are related to PT symmetry.However, as we show in Sec.VIII A where we consider a Kitaev chain with long-range hopping and pairing, this does not have to be the case.

Dynamical criticality
The critical exponents ν ′ and z ′ , which describe the behavior of the soft modes k s,± and the frequencies ω k s,± in the vicinity of transitions between directional pumping phases according to Eq. ( 167), take on the same values for the Kitaev chain as for the SSH model.To obtain this result, we insert µ = µ c + δ µ and γ = γ c + δ γ in Eq. ( 171) and expand in k around k s,±,c = 0, ±π.We obtain for both g = µ and g = γ, which leads to the values of ν ′ and z ′ given in Eq. ( 169).As we show in the next section, the values of these exponents can be modified in the presence of long-range hopping and pairing.

VIII. DRIVEN-DISSIPATIVE KITAEV CHAIN WITH LONG-RANGE HOPPING AND PAIRING
We have introduced the concept of directional pumping phases to distinguish parameter regions with qualitatively different dynamics of string order and fermion parity.Transitions between directional pumping phases are characterized by divergent time scales t s,± for string order and parity pumping, and there is evidence for universality of the exponents that govern the critical behavior of t s,± : We have found the same exponents given in Eq. ( 169) for two models that differ by the presence of a weak U(1) symmetry but belong to the same Altland-Zirnbauer class; and these exponents can be obtained from an expansion in momenta around critical values k s,±,c , which indicates that models with the same lowmomentum or, equivalently, long-wavelength description will have the same exponents, while microscopic details that require the entire Brillouin zone for their description are irrelevant.However, in the examples we have considered so far, the boundaries of directional pumping phases coincide with gap closings of the Liouvillian single-particle spectrum in the complex plane, suggesting that also directional pumping phase transitions are a mere manifestation of gap closings, and that the exponents that describe the critical behavior of t s,± are determined by those that govern the divergences of oscillation periods of low-order correlation functions [85].To better understand the relation between directional pumping phases and the more elementary notion of dynamical phases defined in terms of gap closings or PT symmetry, we now consider quench dynamics in a Kitaev chain with long-range hopping and pairing [108][109][110][111][112][113][114][115][116][117][118][119][120][121][122][123][124].Long-range couplings are known to modify critical properties at gap closings in isolated systems.As we discuss in the following, for a Kitaev chain with Markovian drive and dissipation, the presence of long-range couplings can likewise lead to modifications of the critical exponents that are associated with the pumping time scales t s,± , but does not affect the exponents that govern the divergence of the oscillation period of the density autocorrelation function.Furthermore, even the boundaries of the directional pumping phases of the driven-dissipative long-range Kitaev chain do not always coincide with gap closings.These results indicate that directional pumping phases and the associated critical behavior are indeed new and independent concepts.

A. Long-range Kitaev chain
We consider a Kitaev chain with long-range hopping and pairing as described by the Hamiltonian [122] where both the hopping matrix element J r = J/(N α r α ) and the pairing amplitude ∆ r = ∆/(N α r α ) decay with distance as a power law with exponent α > 1, and the Kac normalization factors are defined as N α = ⌊L/2⌋ r=1 r −α [125].As with the short-range Kitaev chain, we assume that the coupling to Markovian reservoirs is described by the jump operators given in Eq. ( 61), where we focus on γ l = γ g in this section.Consequently, the single-particle dispersion relation of the Liouvillian takes the form given in Eq. ( 65), but with the Hamiltonian dispersion relation ε k given by where, in the thermodynamic limit L → ∞, Explicit expressions for the polylogarithm Li α (z) and the Riemann zeta function ζ(α) are provided in Appendix E.

B. Phase diagram of the long-range Kitaev chain
Long-range couplings do not affect the symmetry properties of the Kitaev chain; in particular, the driven-dissipative longrange Kitaev chain is PT-symmetric.The phase diagram of the long-range Kitaev chain, determined by the spontaneous breaking of PT symmetry, is shown in Figs.14(a) and (c) for α = 3.3 and α = 1.5, respectively.In comparison to the phase diagram of the short-range Kitaev chain [37,85], the key qualitative differences are the absence (i) of mirror symmetry with respect to axis µ = 0 and (ii) of a direct transition between the PT-symmetric and the PT-broken phase, and (iii) a segment of the boundary of the topological PT-symmetric phase being curved instead of described by straight lines.
To understand these modifications, let us first discuss how the corresponding properties come about in the short-range Kitaev chain: (i) The unitary transformation c l → (−1) l c l maps the Hamiltonian H Eq. ( 60) and the chemical potential µ to −H and −µ, respectively.Therefore, gap closings that determine phase boundaries occur symmetrically with respect to µ = 0.This symmetry of the isolated Kitaev chain extends to the driven-dissipative model.However, in the isolated longrange Kitaev chain, the mapping c l → (−1) l c l does not result in a simple sign change and, therefore, the phase diagram does not have reflection symmetry with respect to µ = 0.The We note that also the midpoint between µ c,< and µ c,> , given by µ m = − (J π + J 0 ), does not describe a mirror symmetry of the phasediagram, but will prove to be convenient in the characterization of the long-range Kitaev chain.
(ii) At µ = 0, the dispersion relation of the short-range Kitaev chain is flat, Therefore, upon increasing γ, the Liouvillian dispersion relation Eq. ( 65) becomes imaginary simultaneously for all k ∈ BZ, leading to a direct transition from the PT-symmetric to the PT-broken phase.In contrast, for 1 < α < ∞, the dispersion relation Eq. ( 175) of the long-range Kitaev chain is never flat.The resulting absence of a direct transition the PT-symmetric and PT-broken phase is clearly visible in Fig. 14. (iii) Spontaneous breaking of PT symmetry, which determines the boundary of the PT-symmetric phase, occurs when ω k Eq. ( 65) where γ = γ l = γ g with ε k given in Eq. (175) becomes imaginary for some k.This happens when min k∈BZ (ε k ) = 2γ.For the short-range Kitaev chain, the dispersion relation ε k takes its minimum value at k min = 0 when µ < 0, and at k min = π when µ > 0; the jump of k min from 0 to π occurs at µ = 0 where the dispersion is flat.In contrast, in the long-range Kitaev chain with 1 < α < ∞, there are two bifurcation points of the minimum of ε k at µ b,≶ : There is a single minimum at k = 0 for µ < µ b,< ; then, in the range µ b,< < µ < µ b,> , there are two degenerate minima at momenta ±k b with k b increasing monotonically from 0 to π for µ increasing from µ b,< to µ b,> ; finally, for µ b,> < µ, there is again a single minimum at k = π.This is illustrated for α = 3.3 and α = 1.5 in Figs.14(b) and (d), respectively.To determine the precise shape of the phase boundary, we note that ∆ k defined in Eq. ( 176) vanishes at k = 0, π.Therefore, for µ < µ b,< and µ b,> < µ, the boundary of the PTsymmetric phase is determined by ε 0 = |2J 0 + µ| = 2γ and ε π = |2J π + µ| = 2γ, respectively.These conditions describe straight lines and are symmetric with respect to µ m defined above.But for µ b,< < µ < µ b,> , PT symmetry breaking occurs for ε k b = 2γ, yielding a smaller critical value of γ as compared to an extension of the straight boundaries over the entire range of values of µ.
The momenta ±k b at which the dispersion relation takes on its minimum value for µ b,< < µ < µ b,> can only be found numerically.However, the bifurcation points µ b,≶ can be found analytically by using the series expansions of the polylogarithm given in Appendix E: To determine the right bifurcation point µ b,> , we employ the expansion of ε k around k = π given by ε k ∼ ε π + (k − π) 2 /(2m π ), and solve the equation 1/m π = 0 for µ.Using Eq. (E4), for 1 < α < ∞, we find where primes denote derivatives with respect to k.The left bifurcation point µ b,< is determined by ε k ∼ ε 0 + k 2 /(2m 0 ) and 1/m 0 = 0, which holds for 3 < α < ∞, where the term ∼ k α−1 in Eq. (E2) can be neglected.We obtain For 1 < α < 3, using the same idea as above but keeping only the relevant terms in the expansion in Eq. (E2), we are led to where c 1 , c 2 ∈ R are constants.Therefore, In Fig. 15, our analytical expressions for the bifurcation points are shown to agree with numerical solutions.

C. Subsystem fermion parity
The time evolution of the subsystem parity after quenches to the topological PT-symmetric phase for α = 2.5 and α = 1.5 is shown in Fig. 16.In comparison to the short-range Kitaev chain with Markovian drive and dissipation, there are more pronounced oscillations after t = t F .However, the stationary values of the rescaled subsystem parity are still well described by the PTGGE.

Topological zero crossings
The time scale of topological zero crossings of the subsystem parity is determined by soft modes of the PTGGE.In the presence of long-range couplings, the conditions Eq. (171), which determine the soft modes and associated frequency scales for the short-range Kitaev chain, are generalized to To solve these equations, we use the following properties of function J k and ∆ k defined in Eq. ( 176), which hold for 1 < α < ∞: (i) J k = J −k and ∆ k = −∆ −k are even and odd, respectively.(ii) J k is monotonic for k ∈ [0, π], and can thus be inverted on that interval.We denote the inverse by k where we have introduced a value µ s of the chemical potential that depends on γ and at which the designation of the solutions of Eq. (182) as k s,± is reversed.As in the case of the Kitaev chain with short-range couplings, we have chosen k s,− and k s,+ to describe the pumping of parity through the right and left ends of a subsystem, respectively.The value of µ s In contrast, for OBC, zero crossings occur at multiples of t s,− or t s,+ for a subsystem located at the left (L ℓ , violet line) and right end of the chain (R ℓ , blue line).Factors e ±2γt compensate for additional exponential decay and growth due to edge modes.The analytical conjecture Eq. ( 162) agrees well with the numerics also for the longrange model after stronger initial discrepancies (red line).can be found numerically be requiring the frequencies ω k s,± to be continuous functions of µ and γ.As illustrated in Fig. 17 for the long-range Kitaev chain with α = 2.5, the soft mode periods t s,± = π/ω k s,± agree with the zero crossings of the subsystem parity, both for PBC and OBC.
A unique property of the long-range model is the existence of a finite region in the µ-γ plane which lies outside of the PT-symmetric phase but in which the PTGGE has two soft modes as shown in Figs.18(a), (b) for α = 2.5 and (c), (d) for α = 1.5, where the boundary of the PT-symmetric phase is indicated by a black line.This finding demonstrates that directional pumping phases are in fact not bound to dynamical phases determined by the gap structure of the matrix z k or PT symmetry.Moreover, this finding rules out a relation between topological zero crossings and exceptional points, which are present in the spectrum of the matrix z k in the gapless PTmixed phase at the crossing of the bands λ ±,k .For the SSH model, the band crossing is illustrated in Fig. 2(c).

Dynamical criticality
The product ν ′ z ′ of critical exponents describes how the soft-mode time scales t s,± diverge at the boundaries of directional pumping phases.For the long-range Kitaev chain, this divergence is the same as in the short-range Kitaev chain for the phase boundaries that are shown as green lines in Fig. 18; but for the phase boundaries that are shown as red lines, long-range hopping and pairing lead to modified critical behavior.To derive the corresponding exponents, we note that according to Eq. (182), at directional pumping phase boundaries, ∆ k s,± = ω k s,± /2 vanishes, which is the case for k s,±,c = 0, ±π.We can then determine the scaling behavior of k s,± as in the short-range Kitaev chain by inserting µ = µ c + δ µ and γ = γ c + δ γ in Eq. ( 182) and expanding in k around k s,±,c = 0, ±π.Modified scaling behavior of k s,± as compared to the short-range model occurs if the expansions of J k and ∆ k are modified as compared to J cos(k) and ∆ sin(k), respectively.This is the case for the expansions around k = 0: Depending on the value of α, different exponents dominate in Eq. (E2), leading to and For the phase boundaries with modified scaling behavior, which are indicated by red lines in Fig. 18, we thus find Otherwise, when k s,±,c = ±π which is the case for the phase boundaries that are shown as green lines, the leading powers in the expansions of J k and ∆ k and, therefore, the critical exponents ν ′ and z ′ are the same as in the short-range Kitaev chain and the SSH model and given in Eq. ( 169).
An efficient way to probe the modified exponents numerically is to perform quenches for a range of values of µ and γ close to a phase boundary and measure the time t 1 at which the first zero crossing of the subsystem fermion parity occurs.In Fig. 19(a), we set µ = µ c,< + δ µ and γ = 0, and we consider the fermion parity of the left half of the system L L/2 = {1, . . ., L/2} in a chain with OBC and for different values of α.By decreasing the value of δ µ , we expect to observe a scaling behavior t 1 ∼ δ −ν ′ z ′ µ with exponents given in Eq. ( 186).Indeed, we find good agreement between the predicted exponents and the numerical data for α > 2. For α < 2, due to the limited system sizes in our simulations, we cannot fully reach the scaling regime, we still observe a clear tendency toward the analytically predicted scaling behavior.In Fig. 19(b), we present an analogous analysis for quenches to the PT-mixed phase, where we set µ = µ c,< and γ = δ γ .Here, pumping of fermion parity occurs only through the left end of a subsystem and, therefore, we consider a subsystem R L/2 = {L/2 + 1, . . ., L} corresponding to the right half of the chain.Even though numerics in the PT-mixed phase are restricted to smaller system sizes, we again find compelling agreement with the analytical exponents for α > 2, and a clear trend toward the analytical prediction for α < 2.

D. Connected density autocorrelation function
For each momentum mode k, the Liouvillian single-particle eigenvalues λ ±,k = −i2γ ± ω k are formally identical to the eigenfrequencies of a damped harmonic oscillator with undamped natural frequency ε k and damping rate 2γ.When the damping rate is increased, the mode k undergoes a transition from under-to overdamped oscillations-or, equivalently, from PT-symmetric to PT-breaking-when ε k = 2γ in Eq. (65).It is interesting to compare how different observables are affected by this overdamping transition.The oscillation frequencies of the topological disorder parameters of the SSH model and the Kitaev chain are determined by the soft modes k s,± , and directional pumping transitions occur when either of these modes becomes overdamped.According to our definition of these modes, k s,− is the first to become overdamped, as shown in the phase diagrams in Figs.11, 13, and 18.But even when k s,− is overdamped while k s,+ is not, for a system with PBC or a subsystem that is not on the left end of a system with OBC, the topological disorder parameters exhibit underdamped oscillatory decay as illustrated in Fig. 12.To observe this behavior, it seems to be crucial that the disorder parameters act nontrivially on sufficiently large subsystems of size ℓ.Indeed, underdamped oscillatory decay persist until leveling, caused by the finiteness of ℓ, sets in.As an important example of a local observable with support on only a small number of lattice sites, we consider the density autocorrelation function in the steady state of the driven-dissipative long-range Kitaev chain.As we show in the following, upon increasing the strength of dissipation γ, the decay of the density autocorrelation function becomes overdamped as soon as the first momentum mode becomes overdamped, i.e., at the transition from the PT-symmetric to the PT-broken phase.Concomitantly, one of the base oscillation periods of the density autocorrelation function diverges.However, the exponent that governs this divergence is fully determined by PT symmetry, and its value is not modified by long-range couplings.
The density autocorrelation function is defined as where n l = c † l c l and the normalization is chosen to obtain a simple expression in terms of Majorana fermions.We evaluate the density autocorrelation function in the steady state, which for γ l = γ g is at infinite temperature.In Appendix F, we derive and solve the equation of motion for the density autocorrelation function.The solution can be expressed in terms of sums over contributions from momentum modes, each containing a factor e −iλ ±,k t .This form immediately implies overdamped late-time decay in PT-mixed and PT-broken phases: In these phases, the dominant contribution stems from the mode with the smallest decay rate, given by γ s in Eq. (145).But this mode is nonoscillatory.In contrast, in the PT-symmetric phase, the density autocorrelation function shows oscillatory decay.As detailed in Appendix F, in the thermodynamic limit, the rescaled density autocorrelation function e 4γt A l (t) can be written as a sum of squares of integrals of the form where f k = 2J k +µ or f k = 1.We obtain the asymptotic behavior of the density autocorrelation function for t → ∞ by evaluating these integrals in a stationary phase approximation [89].
To apply this method, we need to identify the stationary points k 0 determined by ω ′ k | k=k 0 = ε ′ k | k=k 0 = 0, where primes denote derivatives with respect to k.According to our discussion of the dispersion relation in Sec.VIII B, the stationary points in the interval 0 ≤ k 0 ≤ π are given by k 0 ∈ {0, k b , π}.The asymptotic behavior of the integral in Eq. ( 188) is determined by momenta close to the stationary points, and we obtain where the terms I s (k 0 ) oscillate at the frequencies ω k 0 , For α < 3, the second derivative ω ′′ k is undefined at k 0 = 0, and the corresponding contribution to I has to be dropped.We omit the lengthy explicit asymptotic expression for A l (t).But the asymptotic result for I already shows that the density autocorrelation function exhibits oscillatory decay with the base frequencies ω 0 , ω k b , and ω π .
As discussed in Sec.VIII B, the minimum of ε k is at 0, k b , and π, for µ < µ b,< , µ b,< < µ < µ b,> , and µ b,> < µ, respectively.Therefore, upon increasing γ, gap closings with ω k 0 → 0, which mark the transition to the PT-mixed phase, occur at these momenta, similarly to ω k s,± going to zero at directional pumping phase transitions.However, while the exponents that describe the critical behavior of ω k s,± are modified in the presence of long-range couplings as described by Eq. ( 186) and illustrated in Fig. 19, the base frequencies of the density autocorrelation function always vanish as a square root of the deviation δ γ = γ − γ c of the dissipation rate γ from the critical value for the transition to the PT-mixed phase, The value of 1/2 of the exponent is, in fact, determined by PT symmetry.For the Kitaev chain, the PT symmetry condition takes the same form as given in Eq. ( 25) for the SSH model, but for the matrix x ′ k = x k • σ corresponding to the traceless part of the matrix x k defined in Appendix F, and given explicitly by x k = (2∆ k , −2J k − µ, −i2γ) [37].Without resorting to the explicit form of x k for the Kitaev chain, PTS implies that x x,k , x y,k ∈ R whereas x z,k ∈ iR, such that At the transition to the PTmixed phase, the difference under the square root vanishes for a particular momentum k 0 .For the example of the Kitaev chain, ix z,k = 2γ, and a gap closing occurs at γ = γ c , where x 2 x,k 0 + x 2 y,k 0 − 4γ 2 c = 0.If we now set γ = γ c + δ γ and approach the gap closing by taking the limit δ γ → 0, we find in agreement with the analysis of response functions in Ref. [85].We stress that this result follows directly from the reality conditions imposed upon x x,k , x y,k , and ix z,k by PT symmetry, and holds also for more general forms of x z,k .Finally, we note that for an isolated system with x z,k = 0 such that ω k = x 2 x,k + x 2 y,k , gap closings require that x x,k 0 = x y,k 0 = 0. Considering again the Kitaev chain with x x,k = 2∆ k and x y,k = −2J k − µ, gap closings occur at k 0 = 0, π where x x,k 0 = 2∆ k 0 = 0.Then, with µ = µ c,≶ + δ µ , we obtain ω k 0 ∼ δ µ .In contrast, the soft-mode frequencies vanish for γ = 0 as ω k s,± ∼ δ µ 1/2 , which provides further evidence for the independence of dynamical criticality at the boundaries of directional pumping phases.

IX. CONCLUSIONS AND OUTLOOK
Through the study of quantum quenches in drivendissipative many-body systems, our work determines PT symmetry as the principal driver of local relaxation to a PTsymmetric generalized Gibbs ensemble for two fundamental classes of quadratic fermionic models.In this way, we substantially extend the field of quantum quenches and relaxation of many-body systems by establishing local equilibration also for open systems with finite coupling to external reservoirs.
We have presented the theoretical framework of the PTGGE for driven-dissipative versions of the SSH model and the Kitaev chain, which can be regarded as natural open-system generalizations of paradigmatic examples of one-dimensional topological insulators and superconductors.These models differ by the presence of a weak U(1) symmetry [69], which leads to the vanishing of anomalous correlations in the open SSH model, even though the coupling to reservoirs breaks particle number conservation.Our analysis shows that PT symmetry of the quadratic Liouvillian is the fundamental property leading to coherent local relaxation in the presence of temporally uniform and spatially global exponential decay.
After rescaling observables to compensate this exponential decay, key features of the quench dynamics of isolated systems are revealed to persist in the PT-symmetric phase.This includes the light cone spreading of correlations and the linear-growth and volume-law saturation of the contribution to the subsystem entropy due to the propagation of pairs of entangled quasiparticles-however, with modified quasiparticle dynamics and statistics.Based on a dissipative quasiparticle picture [38][39][40], we have proposed an analytical conjecture for the time evolution of the quasiparticle-pair contribution to the subsystem entropy in the space-time scaling limit, which is in excellent agreement with our numerical results.
Furthermore, we have provided a detailed analysis of the dynamics of topological disorder parameters-the dual string order parameter and the subsystem fermion parity for the SSH model and the Kitaev chain, respectively.In the isolated versions of these models, the topological disorder parameters show oscillatory decay for quenches from the trivial to the topological phase.This pumping phenomenon becomes directional in driven-dissipative systems, i.e., the pumping of string order and fermion parity happens at different rates through the left and right ends of a subsystem.The pumping rates are determined by soft modes of the PTGGE, and based on this insight, we have formulated analytical conjectures for the time evolution of the topological disorder parameters, which we have found to match numerical simulations very well.Interestingly, there are parameter regimes in which there is pumping through only one end of a subsystem.This has led us to introduce the notion of directional pumping phases.As we have demonstrated using the example of a Kitaev chain with long-range hopping and pairing, directional pumping phases do, in general, not coincide with the phases determined by the breaking of PT symmetry.Moreover, we have identified a distinct form of dynamical criticality at the transitions between directional pumping phases, and we have shown that the critical exponents that govern the divergences of pumping rates are modified when the effective long-wavelength description is modified due to the presence of long-range couplings.
Our work opens up interesting prospects for future research.An important next step is to establish the generality of our results, in particular, of relaxation to the PTGGE.Our reasoning leading to the PTGGE applies to all fermionic many-body systems described by a quadratic Liouvillian with a fully PTsymmetric phase, and can be extended straightforwardly to more general quench protocols, such as for systems prepared in excited or mixed states; and we have presented first results for generalizations of the PTGGE to quadratic bosonic systems and to noninteracting fermionic systems with quadratic Hermitian jump operators in Ref. [37].However, it remains to be seen whether relaxation to a suitably defined PTGGE occurs more generally in PT-symmetric integrable drivendissipative systems [53,63,[126][127][128][129][130][131][132][133][134][135][136][137].
Furthermore, a possible topological origin of the zero crossings of the topological disorder parameters that occur after quenches to the gapless PT-mixed phase warrants further investigation.In Ref. [85], zero crossings of the subsystem parity were studied in a Kitaev chain with OBC and for a subsystems given by the left half of the chain.Then, crossings occur only within the gapped fully PT-symmetric phase, and have been attributed to the nontrivial non-Hermitian topology of the postquench Liouvillian.Our finding, that zero crossings occur for a subsystem at the right end of a chain with OBC and or in a chain with PBC also within the PT-mixed phase, raises the question whether the dynamical entanglement-spectrum bulk-boundary correspondence [4,105,138] can be extended to such phases and how it can be refined to distinguish between left and right entanglement cuts.
In addition to topological disorder parameter pumping or, equivalently, entanglement spectrum crossings, nonanalyticities of the Loschmidt echo, which have been dubbed dynamical quantum phase transitions [8,139,140], have also been proposed as a dynamical signature of topology [6].Indeed, for the isolated Kitaev chain, the soft-mode period t s coincides with the period of singularities of the Loschmidt echo.It will be interesting to see whether a suitable generalization of the Loschmidt echo to open systems [141][142][143] can capture the existence of two distinct soft-mode periods t s,+ t s,− in the driven-dissipative Kitaev chain.expectation value can be expressed in terms of the covariance matrix as ⟨O(t)⟩ = 1 2 l,l ′ O l,l ′ δ l,l ′ − G l ′ ,l (t) .(B1) Further, we focus on balanced loss and gain rates, such that M l − M g = 0 in the evolution equation (13).Then, the solution to the evolution equation that satisfies the initial condition G(0) = G 0 reads (B2) where we have employed the spectral decomposition Z = 2L l=1 λ l v R,l v † L,l .Here, v L,l and v R,l are left and right eigenvectors of Z, respectively, corresponding to the eigenvalues λ l .The left and right eigenvectors obey the biorthogonality condition v † L,l v R,l ′ = δ l,l ′ .We can thus write the evolution of the expectation value of the observable O as (B3) Therefore, the evolution of a quadratic operator can be expressed in terms of matrix elements of the initial covariance matrix G 0 and the observable O between left and right eigenvectors of Z, respectively.For Gaussian states, expectation values of observables that involve products of more than two fermionic fields can be reduced to expectation values of quadratic operators by employing Wick's theorem.

Appendix C: Diagonalization of z k
The non-Hermitian Bloch Hamiltonian of the drivendissipative SSH model, z k = z 1 1 + z k • σ defined in Eq. (20), can be diagonalized in the PT-symmetric phase by applying two consecutive rotations by a real and an imaginary angle, respectively, as described by the nonunitary matrix U k = U z,−ϕ k U y,−iθ k +π/2 = e −i(ϕ k /2)σ z e [(θ k −iπ/2)/2]σ y , (C1) where the angles ϕ k and θ k are defined by the relations ε k e iϕ k = J 1 + J 2 cos(k) + iJ 2 sin(k), (C2) This leads to the following representation of z k : with the dispersion relation of the driven-dissipative SSH model ω k given in Eq. (33).For the isolated SSH model with ∆γ = 0, the above relation reduces to Eq. (27).As explained in Sec.III A, due to the weak U(1) symmetry of the driven-dissipative SSH model, if there are no anomalous correlations in the initial state, no such correlations will be generated in the dynamics, and, therefore, the state of the system is at all times fully determined by the covariance matrix defined in Eq. ( 6).Consequently, describing the state in terms of a covariance matrix for Majorana fermions appears to be unnatural and inefficient.However, as stated in Eq. ( 156), the dual string order parameter has a simple representation in terms of Majorana fermions, which, by applying Wick's theorem, immediately leads to the expression Eq. (157) for the expectation value of the dual string order parameter as the Pfaffian of a submatrix Γ ℓ of the Majorana covariance matrix.In this appendix, we provide a proof for Eq. ( 159) which relates Γ ℓ to the reduced covariance matrix G 1,ℓ defined in terms of complex fermions-confirming that there is no information in Γ ℓ that is not already contained in G 1,ℓ .For simplicity, we consider a translationally invariant state of a system with PBC, but these assumptions are not crucial.
We begin by defining a complex covariance matrix that includes anomalous correlations ⟨c l c l ′ ⟩ and ⟨c † l c † l ′ ⟩, even though these vanish for the SSH model: with . The Majorana covariance matrix for the SSH model, introduced in Sec.VII A, can then be written as where the matrix converting between Majorana and complex fermions is given by We want to understand the structure of the submatrix Γ ℓ introduced in Eq. (158).To that end, we first consider the structure of F. For a translationally invariant state, F can be decomposed into 4 × 4 blocks f l which can be expressed in terms of the elements of the covariance matrix Eq. ( 6) as where the zero entries are due to the vanishing of anomalous correlations.Next, we consider the transformation to Majorana fermions in Eq. (D2).For a 4 × 4 block γ l of Γ, we obtain To calculate the dual string order parameter according to Eq. (157), we require not the full matrix Γ but rather the submatrix Γ ℓ specified in Eq. ( 158).This submatrix is composed in terms of which Eq. (F7) can be recast as  (F14) The asymptotic behavior of these integrals for t → ∞ can be obtained by means of a stationary phase approximation as described in Sec.VIII D.

Figure 2 .
Figure 2. (a) Dynamical phase diagram of the driven-dissipative SSH model with topological phase for ∆J < 0 and trivial phase for ∆J > 0. The model features PT-symmetric (blue, red), PT-broken (orange, purple), and PT-mixed phases (green, yellow).Examples of the mode structure for the three topological phases are shown for (b) the PT-symmetric (PTS), (c) the PT-mixed (PTM), and (d) the PT-broken (PTB) phases.Positive (negative) bands are represented by orange (blue) dots and edge modes for a system with OBC by black dots.Star, circle and diamond markers in (a) indicate parameters corresponding to the spectra shown in (b)-(d).

Figure 3 .
Figure 3. Relaxation of local correlations for a PT-symmetric quench with ∆J = −0.5J,∆γ = γ = 0.2J, and with small imbalance between loss and gain.(a) The dashed line corresponds to an imbalance of ∆ = δ/2 = 10 −9 J, while the solid line shows the evolution for balanced loss and gain, ∆ = δ = 0.For finite imbalance, the rescaled correlation function relaxes to the PTGGE prediction on intermediate time scales up to t × , indicated by a green vertical line.(b) The analytical prediction from Eq. (82) is in excellent agreement with the numerically determined values for a range of values of δ = 2∆.Numerically, we define t × as the time at which the deviation from the PTGGE value equals one.
A. PT-symmetric phase In Figs.4(a) and (b), we show the absolute value of rescaled normal commutators C ′ ℓ (t) = e 4γt C ℓ (t) for the isolated and the driven-dissipative Kitaev chain in the PT-symmetric phase, respectively.A clear light cone structure is visible in both figures.The peak of correlations, defining the boundary of the light cone and found at position ℓ = 2v max t, is well described by the ballistic propagation of quasiparticles with finite lifetime ∼ 1/γ and velocity

Figure 5 .
Figure 5. Maximum of the quasiparticle velocity v max and the phase velocity v PH as a function of the chemical potential µ for (a) the isolated system (γ = 0) and (b) the driven-dissipative system (γ = 0.3J).For µ < 0, at the phase boundary, v max = v PH = 2J |µ|, while for µ > 0, we find v max = 2J |µ| and v PH = 0.

Figure 7 .
Figure 7. Quasiparticle-pair contribution to the subsystem entropy in the driven-dissipative SSH model after quenches to the trivial (green, ∆J = 0.5J) and topological (blue, ∆J = −0.5J)PT-symmetric phases for γ = ∆γ = 0.1J, δ = ∆ = 0, and ℓ = 20.Dashed lines show numerical data, and solid lines correspond to the analytical conjecture Eq. (151) for the space-time scaling limit.The approach of the numerical results to the conjecture with increasing subsystem size is shown in the inset: For the topological quench at t = 2t F , the difference between the numerical data and Eq.(151) (blue dots) vanishes as 1/ℓ (orange line).

Figure 8 .
Figure 8. Dual string order parameter after quenches to the trivial (green, ∆J = 0.5J) and topological (blue, ∆J = −0.5J)PTsymmetric phases for γ = ∆γ = 0.3J, δ = ∆ = 0, and ℓ = 20.Dashed lines correspond to numerical data, and the solid lines show the analytical conjectures in Eqs.(161) and (162) with α + = α − = 0.15.We indicate the asymptotic PTGEE prediction by straight horizontal red lines and the the change of dynamics at the characteristic time scale t = t F by the purple vertical line.The system size L is chosen sufficiently large to avoid finite-size effects.

Figure 9 .
Figure 9. Simultaneous zero crossings of the dual string order parameter (black line) and in the single-particle entanglement spectrum (blue lines) for a quench to the topological phase, ∆J = −1.5Jand γ = ∆γ = 0.3J, with balanced loss and gain δ = ∆ = 0, and for a subsystem of size ℓ = 10.

Figure 11 .
Figure 11.Directional string order pumping phase diagram of the driven-dissipative SSH model determined by the soft-mode frequencies (a) ω ks,− and (b) ω ks,+ .The phases determined by PT symmetry, introduced in Fig. 2, are shown in the background.

Figure 13 .
Figure 13.Directional parity pumping phases determined by the softmode frequencies (a) ω ks,− and (b) ω ks,+ for the driven-dissipative Kitaev chain.The background colors indicate the phases defined by PT symmetry, with PT-symmetric phases (blue, red), PT-mixed phases (green, yellow) and a PT-broken phase (orange) [37, 85].For better visual differentiation the saturation of background colors in (b) has been increased.

Figure 16 .
Figure 16.Relaxation of the subsystem parity after quenches in the driven-dissipative Kitaev chain with long-range couplings.Both quenches are to the topological phase, where the blue line corresponds to α = 2.5, µ = −J, γ = 0.3J, and ℓ = 20, and the orange line to α = 1.5, µ = −0.5J,γ = 0.15J, and ℓ = 80.Dashed and solid lines show numerical data and the analytical conjecture in Eq. (162), respectively.For smaller values of α corresponding to longer-range couplings, finite-size effects are more strongly pronounced and larger system sizes are required.The purple vertical line indicates the characteristic relaxation time scale t F and the red horizontal lines show the PTGGE predictions.

Figure 17 .
Figure 17.Directional pumping of subsystem parity for a quench to the topological PT-symmetric phase with α = 2.5, µ = −0.5J,γ = 0.3J, δ = 0, and ℓ = 20.For PBC, the subsystem parity (black line) crosses zero at multiples of both t s,+ (light green) and t s,− (dark green), with blue shading indicating the sign of the subsystem parity.In contrast, for OBC, zero crossings occur at multiples of t s,− or t s,+ for a subsystem located at the left (L ℓ , violet line) and right end of the chain (R ℓ , blue line).Factors e ±2γt compensate for additional exponential decay and growth due to edge modes.The analytical conjecture Eq. (162) agrees well with the numerics also for the longrange model after stronger initial discrepancies (red line).

Figure 18 .
Figure 18.Directional parity pumping phase diagrams of the longrange Kitaev chain determined by the soft-mode frequencies ω ks,± for (a), (b) α = 2.5 and (c), (d) α = 1.5.Background colors indicate phases defined by PT symmetry as in Fig. 14.The boundary of the PT-symmetric topological phase is shown as a black line.Red and green lines indicate phase boundaries at which the critical exponents are and are not modified, respectively, in comparison to the shortrange Kitaev chain.

Figure 19 .
Figure 19.Scaling behavior of the first crossing time t 1 for the subsystem parity in the long-range Kitaev chain.The connected dots show numerical results for α = 4.9 (blue), α = 2.5 (green) and α = 1.9 (orange), and the expected scaling behavior according to Eq. (186) is indicated by solid lines, with scaling exponents given in the left lower corner of each panel.We consider quenches to (a) the PT-symmetric (PTS) phase with µ = µ c,< + δ µ and γ = 0, and (b) the PT-mixed (PTM) phase with µ = µ c,< and γ = δ γ , for a chain with OBC and subsystems located at (a) the left (L ℓ with ℓ = L/2 = 500) and (b) the right end of the chain (R ℓ with ℓ = L/2 = 300).