MnLargeSymbols’164 MnLargeSymbols’171

# Supplemental Material: Relaxation to a Parity-Time Symmetric Generalized
Gibbs
Ensemble after a Quantum Quench in a Driven-Dissipative Kitaev Chain

Elias Starchl  Lukas M. Sieberer<EMAIL_ADDRESS>Institute for
Theoretical Physics, University of Innsbruck, 6020 Innsbruck, Austria

###### Abstract

In the Supplemental Material, we provide technical details of our analytical
and numerical studies of quench dynamics in a driven-dissipative Kitaev chain,
and we present additional results for PT-symmetric relaxation in quench
dynamics of fermionic and bosonic SSH models with loss and gain and a
dimerized XY spin chain with dephasing. We first discuss relevant properties
of the isolated Kitaev chain, in particular, its symmetries and the mapping to
the transverse field Ising model employed in the main text. Then, we summarize
how these properties are affected by Markovian drive and dissipation, and how
quench dynamics in the driven-dissipative Kitaev chain can be described
analytically and numerically in terms of the time evolution of the covariance
matrix. Focusing on the PT-symmetric phase, we present two complementary
derivations of the PT-symmetric generalized Gibbs ensemble (PTGGE): The first
one is based on the dephased form of the covariance matrix, and the second
one, which is outlined in the main text, benefits from the clear picture of
quench dynamics that is provided by a description in terms of eigenmodes of
the Liouvillian. We then describe how the numerical results for the subsystem
parity and entropy shown in the main text are obtained, and we motivate our
conjectures for the time-dependence of these quantities in the space-time-
scaling limit, which are summarized in Eqs. (9), (10), and (11) in the main
text. Finally, we present our results for relaxation to the PTGGE in other
models.

###### Contents

  1. I Isolated Kitaev chain
    1. I.1 Hamiltonian
    2. I.2 Covariance matrix
    3. I.3 Diagonalization and ground state
      1. I.3.1 Numerical diagonalization
      2. I.3.2 Analytical diagonalization for periodic boundary conditions
    4. I.4 Symmetries of the isolated Kitaev chain
    5. I.5 Relation to the transverse field Ising model
  2. II Driven-dissipative Kitaev chain
    1. II.1 Master equation and Liouvillian
    2. II.2 Time evolution of the covariance matrix
      1. II.2.1 Numerical time evolution
      2. II.2.2 Analytical time evolution for periodic boundary conditions
      3. II.2.3 Spectrum of the covariance matrix
    3. II.3 Time evolution of the complex covariance matrix
    4. II.4 Symmetries of the Liouvillian
    5. II.5 Spectrum of the Liouvillian
  3. III Quench dynamics in the PT-symmetric phase and relaxation to the PTGGE
    1. III.1 Dephasing and relaxation to the PTGGE
      1. III.1.1 The covariance matrix in the PT-symmetric phase
      2. III.1.2 Late-time limit and dephasing
      3. III.1.3 From the covariance matrix to the density matrix
    2. III.2 Derivation of the PTGGE from the principle of maximum entropy
      1. III.2.1 Eigenmodes of the adjoint Liouvillian
      2. III.2.2 Maximum entropy ensemble
    3. III.3 Subsystem parity
      1. III.3.1 Numerical time evolution
      2. III.3.2 Long-time behavior
      3. III.3.3 Space-time scaling limit
      4. III.3.4 Condition to observe relaxation to the PTGGE
      5. III.3.5 Oscillatory dynamics from zero modes
      6. III.3.6 Necessary conditions for $t_{s,+}\neq t_{s,-}$
    4. III.4 Subsystem entropy
      1. III.4.1 Subsystem entropy from the covariance matrix
      2. III.4.2 Statistical entropy
      3. III.4.3 Numerical time evolution
      4. III.4.4 Late-time behavior
      5. III.4.5 Space-time scaling limit
  4. IV Universality of PT-symmetric relaxation
    1. IV.1 Fermionic SSH model with loss and gain
    2. IV.2 Bosonic SSH model with loss and gain
    3. IV.3 Dimerized XY model with dephasing
  5. A Time evolution induced by a non-Hermitian $2\times 2$ Hamiltonian

##  I Isolated Kitaev chain

Here, we summarize aspects of the isolated Kitaev chain that are relevant for
our study of quench dynamics in a driven-dissipative Kitaev chain. This
includes the diagonalization and ground state of the Kitaev chain for open and
periodic boundary conditions, symmetries, and the mapping between the Kitaev
chain and the transverse field Ising model.

###  I.1 Hamiltonian

The Hamiltonian of the isolated Kitaev chain, as given in the main text, reads

$H=\sum_{l=1}^{L}\left(-Jc^{\dagger}_{l}c_{l+1}+\Delta
c_{l}c_{l+1}+\mathrm{H.c.}\right)-\mu\sum_{l=1}^{L}\left(c^{\dagger}_{l}c_{l}-\frac{1}{2}\right).$
(1)

All of our results for quench dynamics presented in the main text are obtained
for $J=\Delta$, but we keep $J$ and $\Delta$ as distinct parameters in most
expressions below. However, we always assume that
$J,\Delta\in\mathbb{R}_{>0}$. Further, we consider both periodic boundary
conditions with $c_{L+1}=c_{1}$, and open boundary conditions with
$c_{L+1}=0$.

An alternative form of the Hamiltonian, which puts emphasis on translational
invariance and, therefore, proves especially convenient when we consider
periodic boundary conditions, is given by

$H=\frac{1}{2}\sum_{l,l^{\prime}=1}^{L}C_{l}^{\dagger}h_{l-l^{\prime}}C_{l^{\prime}},$
(2)

where we collect fermionic annihilation and creation operators in spinors,

$C_{l}=\left(c_{l},c_{l}^{\dagger}\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(3)

and where

$h_{l}=\mathbf{h}_{l}\cdot\boldsymbol{\sigma},$ (4)

with the vector of Pauli matrices,

$\boldsymbol{\sigma}=\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(5)

and with coefficients given by

$\mathbf{h}_{l}=\left(0,i\Delta\left(\delta_{l,1}-\delta_{l,-1}\right),-J\left(\delta_{l,1}+\delta_{l,-1}\right)-\mu\delta_{l,0}\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}.$
(6)

Here and in the following, for periodic boundary conditions, each index of a
Kronecker delta is to be understood modulo the length of the chain $L$, e.g.,
$\delta_{l,l^{\prime}+L}=\delta_{l,l^{\prime}}$.

In addition to the above formulation of the Kitaev chain in terms of $L$
complex or Dirac fermions $c_{l}$ with $l\in\\{1,\dotsc,L\\}$, it is often
convenient to work with $2L$ Majorana fermions $w_{l}$, where
$l\in\\{1,\dotsc,2L\\}$, which are defined by

$w_{2l-1}=c_{l}+c^{\dagger}_{l},\qquad
w_{2l}=i\left(c_{l}-c^{\dagger}_{l}\right),$ (7)

and obey the anticommutation relation
$\\{w_{l},w_{l^{\prime}}\\}=2\delta_{l,l^{\prime}}$. Further, it is often
convenient to express the transformation to Majorana fermions in terms of
spinors of operators as follows:

$C_{l}=\frac{1}{\sqrt{2}}rW_{l},$ (8)

where $C_{l}$ is defined in Eq. (3),

$W_{l}=\left(w_{2l-1},w_{2l}\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(9)

and the unitary matrix $r$ is given by

$r=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-i\\\ 1&i\end{pmatrix}.$ (10)

In terms of Majorana fermions, the Hamiltonian of the isolated Kitaev chain in
Eq. (1) can be written as

$H=\frac{i}{4}\sum_{l,l^{\prime}=1}^{L}W_{l}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}a_{l-l^{\prime}}W_{l^{\prime}},$
(11)

where

$a_{l}=-ir^{\dagger}h_{l}r=\mathbf{a}_{l}\cdot\boldsymbol{\sigma},$ (12)

with

$\mathbf{a}_{l}=\left(\Delta\left(\delta_{l,1}-\delta_{l,-1}\right),i\left[J\left(\delta_{l,1}+\delta_{l,-1}\right)+\mu\delta_{l,0}\right],0\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}.$
(13)

Here and in the following, the transformations of various quantities between
the bases of complex and Majorana fermions, such as the transformation from
$h_{l}$ to $a_{l}$, can be performed efficiently by noting that they
correspond merely to permutations of the components in expansions in terms of
Pauli matrices, which are here given by $\mathbf{h}_{l}$ and $\mathbf{a}_{l}$,
respectively. In particular, to calculate

$r^{\dagger}h_{l}r=\mathbf{h}_{l}\cdot\left(r^{\dagger}\boldsymbol{\sigma}r\right),$
(14)

we use

$r^{\dagger}\begin{pmatrix}\sigma_{x}\\\ \sigma_{y}\\\
\sigma_{z}\end{pmatrix}r=\begin{pmatrix}\sigma_{z}\\\ \sigma_{x}\\\
\sigma_{y}\end{pmatrix}=P^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\begin{pmatrix}\sigma_{x}\\\
\sigma_{y}\\\ \sigma_{z}\end{pmatrix},$ (15)

where the permutation matrix $P$ is given by

$P=\begin{pmatrix}0&1&0\\\ 0&0&1\\\ 1&0&0\end{pmatrix}.$ (16)

We thus obtain Eqs. (12) and (13):

$r^{\dagger}h_{l}r=\mathbf{h}_{l}\cdot\left(P^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\boldsymbol{\sigma}\right)=\left(P\mathbf{h}_{l}\right)\cdot\boldsymbol{\sigma}=i\mathbf{a}_{l}\cdot\boldsymbol{\sigma}=ia_{l}.$
(17)

For future reference, we define a real antisymmetric matrix $A$ as the block
Toeplitz matrix built from the $2\times 2$ blocks $a_{l}$,

$A=\begin{pmatrix}a_{0}&a_{-1}&\cdots&a_{1-L}\\\ a_{1}&a_{0}&&\vdots\\\
\vdots&&\ddots&\vdots\\\ a_{L-1}&\cdots&\cdots&a_{0}\end{pmatrix}.$ (18)

In terms of the matrix $A$, the Hamiltonian of the Kitaev chain can be written
as

$H=\frac{i}{4}\sum_{l,l^{\prime}=1}^{2L}w_{l}A_{l,l^{\prime}}w_{l^{\prime}}.$
(19)

###  I.2 Covariance matrix

In our studies of quench dynamics, we always deal with Gaussian states $\rho$:
The initial state $\rho_{0}=\lvert\psi_{0}\rangle\langle\psi_{0}\rvert$, where
$\lvert\psi_{0}\rangle$ is the ground state of the isolated Kitaev chain, is
Gaussian, and the quadratic Liouvillian we consider below preserves
Gaussianity. Gaussian states are fully determined by two-point functions,
which are collected in the real and antisymmetric covariance matrix,

$\Gamma_{l,l^{\prime}}=\frac{i}{2}\langle[w_{l},w_{l^{\prime}}]\rangle=i\left(\langle
w_{l}w_{l^{\prime}}\rangle-\delta_{l,l^{\prime}}\right),$ (20)

where $\langle\dotsc\rangle=\mathop{\mathrm{tr}}(\dotsc\rho)$. We often find
it convenient to work with $2\times 2$ blocks of $\Gamma$, which can be
specified in terms of the Majorana spinors defined in Eq. (9):

$\gamma_{l,l^{\prime}}=\begin{pmatrix}\Gamma_{2l-1,2l^{\prime}-1}&\Gamma_{2l-1,2l^{\prime}}\\\
\Gamma_{2l,2l^{\prime}-1}&\Gamma_{2l,2l^{\prime}}\end{pmatrix}=i\left(\langle
W_{l}W_{l^{\prime}}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\rangle-\delta_{l,l^{\prime}}\mathbbm{1}\right).$
(21)

Switching to the basis of complex fermions, we define a matrix $G$ through the
relation

$G=-iR\Gamma R^{\dagger},$ (22)

where the matrix $R$ is given in terms of the $2\times 2$ matrix $r$ in Eq.
(10) by

$R=\bigoplus_{l=1}^{L}r.$ (23)

We refer to $G$ as covariance matrix in the complex basis, or simply complex
covariance matrix. As done in Eq. (21) for the covariance matrix in the
Majorana basis, we define blocks of the complex covariance matrix,

$g_{l,l^{\prime}}=\begin{pmatrix}G_{2l-1,2l^{\prime}-1}&G_{2l-1,2l^{\prime}}\\\
G_{2l,2l^{\prime}-1}&G_{2l,2l^{\prime}}\end{pmatrix}=-ir\gamma_{l,l^{\prime}}r^{\dagger}.$
(24)

In terms of the complex spinors given in Eq. (3), the blocks
$g_{l,l^{\prime}}$ can be written as

$g_{l,l^{\prime}}=2\langle
C_{l}C_{l^{\prime}}^{\dagger}\rangle-\delta_{l,l^{\prime}}\mathbbm{1}=\begin{pmatrix}\langle[c_{l},c_{l^{\prime}}^{\dagger}]\rangle&\langle[c_{l},c_{l^{\prime}}]\rangle\\\
\langle[c_{l^{\prime}}^{\dagger},c_{l}^{\dagger}]\rangle&\langle[c_{l^{\prime}}^{\dagger},c_{l}]\rangle\end{pmatrix}.$
(25)

###  I.3 Diagonalization and ground state

We proceed to explain how the Hamiltonian of the Kitaev chain in Eq. (1) can
be diagonalized and how its ground state can be found. This can be done
numerically for finite systems with open or periodic boundary conditions.
Analytical progress is possible for a system with periodic boundary
conditions. The ground state serves as the initial state in the quench
protocol discussed in the main text and further below.

####  I.3.1 Numerical diagonalization

The excitation spectrum and ground state of the Kitaev chain can be found
numerically as follows: By employing the algorithm described in Ref. [1] (see
also Ref. [2] for a discussion in the context of fermionic Gaussian systems),
one obtains a real special orthogonal matrix $O$, which brings $A$ in Eq. (18)
to the following canonical block-diagonal form:

$OAO^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}=\bigoplus\limits_{l=1}^{L}\begin{pmatrix}0&-\varepsilon_{l}\\\
\varepsilon_{l}&0\end{pmatrix},$ (26)

where $\varepsilon_{l}\geq 0$ are the energies of elementary excitations.
Further, the covariance matrix of the ground state is given by [3]

$\Gamma_{0}=O^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\left(\bigoplus\limits_{l=1}^{L}i\sigma_{y}\right)O.$
(27)

Below and in the main text, we focus on quenches from the ground state of the
Kitaev chain for $\mu_{0}\to-\infty$, i.e., the vacuum state
$\lvert\Omega\rangle$, for which $c_{l}\lvert\Omega\rangle=0$ for all
$l\in\\{1,\dotsc L\\}$. Then, $O$ is simply the $2L\times 2L$ identity matrix.

####  I.3.2 Analytical diagonalization for periodic boundary conditions

For periodic boundary conditions, analytical progress can be made most
conveniently in the complex basis. In particular, we express the fermionic
operators $c_{l}$ in terms of momentum-space operators $c_{k}$ as

$c_{l}=\frac{1}{\sqrt{L}}\sum_{k\in\mathrm{BZ}}e^{ikl}c_{k},\quad
c_{k}=\frac{1}{\sqrt{L}}\sum_{l=1}^{L}e^{-ikl}c_{l},$ (28)

where the Brillouin zone $\mathrm{BZ}$ contains $L$ momenta with spacing
$\Delta k=2\pi/L$, that is, $\mathrm{BZ}=\\{\pi+\Delta k,-\pi+\Delta
k,\dotsc,\pi\\}$. The spinors in Eq. (3) transform to

$C_{k}=\begin{pmatrix}c_{k}\\\
c_{-k}^{\dagger}\end{pmatrix}=\frac{1}{\sqrt{L}}\sum_{l=1}^{L}e^{-ikl}C_{l},$
(29)

and the Hamiltonian of the Kitaev chain can be written as

$H=\sum_{k\geq 0}C_{k}^{\dagger}h_{k}C_{k},$ (30)

with the Bogoliubov-De Gennes (BdG) Hamiltonian

$h_{k}=\sum_{l=1}^{L}e^{-ikl}h_{l}=\mathbf{h}_{k}\cdot\boldsymbol{\sigma},$
(31)

where

$\mathbf{h}_{k}=\left(0,2\Delta\sin(k),-2J\cos(k)-\mu\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}.$
(32)

Below, we find it convenient to parameterize the vector $\mathbf{h}_{k}$ in
terms of its magnitude,

$\varepsilon_{k}=\left\lvert\mathbf{h}_{k}\right\rvert=\sqrt{\left(2J\cos(k)+\mu\right)^{2}+4\Delta^{2}\sin(k)^{2}},$
(33)

and a unit vector

$\hat{\mathbf{n}}_{k}=\frac{\mathbf{h}_{k}}{\varepsilon_{k}}=\left(0,\sin(\theta_{k}),\cos(\theta_{k})\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(34)

where the angle $\theta_{k}$ is defined implicitly by the relation

$\varepsilon_{k}e^{i\theta_{k}}=-2J\cos(k)-\mu+i2\Delta\sin(k).$ (35)

The diagonalization of the BdG Hamiltonian corresponds geometrically to a
rotation of $\mathbf{h}_{k}$ to the $z$ axis, which can be accomplished by a
rotation around the $x$ axis by the angle $\theta_{k}$: In particular, with

$U_{x,\theta_{k}}=e^{i\left(\theta_{k}/2\right)\sigma_{x}}=\begin{pmatrix}\cos(\theta_{k}/2)&i\sin(\theta_{k}/2)\\\
i\sin(\theta_{k}/2)&\cos(\theta_{k}/2)\end{pmatrix},$ (36)

we obtain

$U_{x,\theta_{k}}^{\dagger}h_{k}U_{x,\theta_{k}}=\varepsilon_{k}\sigma_{z}.$
(37)

The corresponding Bogoliubov transformation of fermionic operators reads

$D_{k}=\begin{pmatrix}d_{k}\\\
d_{-k}^{\dagger}\end{pmatrix}=U_{x,\theta_{k}}^{\dagger}C_{k}.$ (38)

In terms of the thusly defined fermionic modes $d_{k}$, the Hamiltonian
becomes

$H=\sum_{k\geq
0}\varepsilon_{k}\left(d_{k}^{\dagger}d_{k}-d_{-k}d_{-k}^{\dagger}\right),$
(39)

which shows that $\varepsilon_{k}$ determines the dispersion relation of
single-particle excitations, and that the ground state $\lvert\psi_{0}\rangle$
is given by the vacuum state of the operators $d_{k}$, i.e.,
$d_{k}\lvert\psi_{0}\rangle=0$ for all $k\in\mathrm{BZ}$. To obtain an
explicit expression for the ground state, we have to treat momenta
$k\in\\{0,\pi\\}$ and $k\in\mathrm{BZ}\setminus\\{0,\pi\\}$ separately. For
momenta $k\in\\{0,\pi\\}$, Eq. (35) implies

$\theta_{0}=\begin{cases}0&\text{for }\mu<-2J,\\\ \pi&\text{for
}\mu>-2J,\end{cases}$ (40)

and

$\theta_{\pi}=\begin{cases}0&\text{for }\mu<2J,\\\ \pi&\text{for
}\mu>2J.\end{cases}$ (41)

By inserting these results in the Bogoliubov transformation Eq. (38), we
obtain

$d_{0}=\begin{cases}c_{0}&\text{for }\mu<-2J,\\\ -ic_{0}^{\dagger}&\text{for
}\mu>-2J,\end{cases}$ (42)

and, with $c_{\pi}=c_{-\pi}$,

$d_{\pi}=\begin{cases}c_{\pi}&\text{for }\mu<2J,\\\
-ic_{\pi}^{\dagger}&\text{for }\mu>2J.\end{cases}$ (43)

Therefore, the conditions
$d_{0}\lvert\psi_{0}\rangle=d_{\pi}\lvert\psi_{0}\rangle=0$ are satisfied, if,
in the state $\lvert\psi_{0}\rangle$, the mode with momentum $k=0$ is occupied
when $\mu>-2J$, and the mode with momentum $k=\pi$ is occupied when $\mu>2J$.

For momenta $k\in\mathrm{BZ}\setminus\\{0,\pi\\}$, we employ the following
alternative form of the Bogoliubov transformation:

$d_{k}=\mathfrak{u}_{k}c_{k}\mathfrak{u}_{k}^{\dagger},$ (44)

where

$\begin{split}\mathfrak{u}_{k}&=e^{\left(\theta_{k}/2\right)\left(c_{-k}c_{k}-c_{k}^{\dagger}c_{-k}^{\dagger}\right)}\\\
&=e^{i\tan(\theta_{k}/2)c_{k}^{\dagger}c_{-k}^{\dagger}}e^{i\sin(\theta_{k}/2)\cos(\theta_{k}/2)c_{-k}c_{k}}e^{-\ln(\cos(\theta_{k}/2))\left(c_{k}^{\dagger}c_{k}+c_{-k}^{\dagger}c_{-k}-1\right)}.\end{split}$
(45)

From $c_{k}\lvert\Omega\rangle=0$ we obtain
$d_{k}\mathfrak{u}_{k}\lvert\Omega\rangle=0$, i.e.,
$\mathfrak{u}_{k}\lvert\Omega\rangle$ is the vacuum of the Bogoliubov modes
$d_{k}$. Further, since $\mathfrak{u}_{k}$ is unitary, the transformed vacuum
$\mathfrak{u}_{k}\lvert\Omega\rangle$ is normalized.

When we study quench dynamics below, we restrict ourselves to initial states
with $\mu<-2J$ such that, according to the above discussion, the modes with
$k\in\\{0,\pi\\}$ are not occupied. Then, the ground state is given by

$\begin{split}\lvert\psi_{0}\rangle&=\prod_{k\in\mathrm{BZ}\setminus\\{0,\pi\\}}\mathfrak{u}_{k}\lvert\Omega\rangle\\\
&=\prod_{k\in\mathrm{BZ}\setminus\\{0,\pi\\}}\left(\cos(\theta_{k}/2)+i\sin(\theta_{k}/2)c_{k}^{\dagger}c_{-k}^{\dagger}\right)\lvert\Omega\rangle.\end{split}$
(46)

Finally, we specify the covariance matrix in the ground state. To that end, we
first note that due to translational invariance in a system with periodic
boundary conditions, the $2\times 2$ block in Eq. (24) depends only on the
difference of its indices, $g_{l,l^{\prime}}=g_{l-l^{\prime}}$. Then, its
Fourier transform is given by

$g_{k}=\sum_{l=1}^{L}e^{-ikl}g_{l}=2\langle
C_{k}C_{k}^{\dagger}\rangle-\mathbbm{1}=\begin{pmatrix}\langle[c_{k},c_{k}^{\dagger}]\rangle&\langle[c_{k},c_{-k}]\rangle\\\
\langle[c_{-k}^{\dagger},c_{k}^{\dagger}]\rangle&\langle[c_{-k}^{\dagger},c_{-k}]\rangle\end{pmatrix}.$
(47)

To obtain an explicit expression for the complex covariance matrix $g_{k}$ in
the ground state, we use the following ground-state expectation values, which
follow from the condition $d_{k}\lvert\psi_{0}\rangle=0$ that defines the
ground state,

$\langle D_{k}D_{k}^{\dagger}\rangle_{0}=\begin{pmatrix}\langle
d_{k}d_{k}^{\dagger}\rangle_{0}&\langle d_{k}d_{-k}\rangle_{0}\\\ \langle
d_{-k}^{\dagger}d_{k}^{\dagger}\rangle_{0}&\langle
d_{-k}^{\dagger}d_{-k}\rangle_{0}\end{pmatrix}=\frac{\mathbbm{1}+\sigma_{z}}{2},$
(48)

where $\langle\dotsc\rangle_{0}=\langle\psi_{0}|\dotsc|\psi_{0}\rangle$.
Further, we use Eq. (37), according to which
$h_{k}/\varepsilon_{k}=\hat{\mathbf{n}}_{k}\cdot\boldsymbol{\sigma}=U_{x,\theta_{k}}\sigma_{z}U_{x,\theta_{k}}^{\dagger}$,
and Eq. (38). We thus find

$g_{0,k}=2U_{x,\theta_{k}}\langle
D_{k}D_{k}^{\dagger}\rangle_{0}U_{x,\theta_{k}}^{\dagger}-\mathbbm{1}=U_{x,\theta_{k}}\sigma_{z}U_{x,\theta_{k}}^{\dagger}=\hat{\mathbf{n}}_{k}\cdot\boldsymbol{\sigma}.$
(49)

When we rotate this expression to the Majorana basis, we find it convenient to
include a prefactor $-i$ relative to Eq. (24), i.e., we define the momentum-
space covariance matrix in the Majorana basis as

$\gamma_{k}=r^{\dagger}g_{k}r=-i\sum_{l=1}^{L}e^{-ikl}\gamma_{l},\qquad\gamma_{l}=\frac{i}{L}\sum_{k\in\mathrm{BZ}}e^{ikl}\gamma_{k}.$
(50)

With this convention, $\gamma_{k}=\gamma_{k}^{\dagger}$ is Hermitian. We note
that in terms of spinors of Majorana operators, which we define in analogy to
Eq. (29),

$W_{k}=\frac{1}{\sqrt{L}}\sum_{l=1}^{L}e^{-ikl}W_{l}=\sqrt{2}r^{\dagger}C_{k},$
(51)

the covariance matrix can be written as

$\gamma_{k}=\langle W_{k}W_{k}^{\dagger}\rangle-\mathbbm{1}.$ (52)

We obtain the covariance matrix in the ground state from Eq. (49) by using
Eqs. (15) and (16),

$\gamma_{0,k}=\hat{\mathbf{n}}_{k}\cdot\left(r^{\dagger}\boldsymbol{\sigma}r\right)=\hat{\mathbf{n}}_{k}\cdot\left(P^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\boldsymbol{\sigma}\right)=\left(P\hat{\mathbf{n}}_{k}\right)\cdot\boldsymbol{\sigma}=\hat{\mathbf{a}}_{k}\cdot\boldsymbol{\sigma},$
(53)

where $\hat{\mathbf{a}}_{k}=\mathbf{a}_{k}/\varepsilon_{k}$ and

$\mathbf{a}_{k}=P\mathbf{h}_{k}=\left(2\Delta\sin(k),-2J\cos(k)-\mu,0\right)^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}.$
(54)

For future reference, we define the Fourier transform of $a_{l}$ in Eq. (12)
including the imaginary unity as a prefactor:

$a_{k}=r^{\dagger}h_{k}r=i\sum_{l=1}^{L}e^{-ikl}a_{l}=\mathbf{a}_{k}\cdot\boldsymbol{\sigma}.$
(55)

Again, this convention leads to a Hermitian matrix $a_{k}=a_{k}^{\dagger}$.

###  I.4 Symmetries of the isolated Kitaev chain

We next discuss symmetries of the Kitaev chain, focusing on the BdG
Hamiltonian $h_{k}$ in Eq. (31). To prepare for a comparison to the driven-
dissipative Kitaev chain introduced below, we state two equivalent forms of
each symmetry, where the second form follows from Hermiticity of the BdG
Hamiltonian, $h_{k}=h_{k}^{\dagger}$. Particle-hole symmetry (PHS), time-
reversal symmetry (TRS), and chiral symmetry (CS), which results from the
combination of PHS and TRS, are expressed by

PHS: $\displaystyle h_{k}$
$\displaystyle=-\sigma_{x}h_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{x}=-\sigma_{x}h_{-k}^{*}\sigma_{x},$
(56) TRS: $\displaystyle h_{k}$
$\displaystyle=h_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}=h_{-k}^{*},$
(57) CS: $\displaystyle h_{k}$
$\displaystyle=-\sigma_{x}h_{k}\sigma_{x}=-\sigma_{x}h_{k}^{\dagger}\sigma_{x}.$
(58)

We note that TRS is ensured by our choice $J\in\mathbb{R}_{>0}$ (a nontrivial
phase of $\Delta$ can always be absorbed in a redefinition of the operators
$c_{l}$). Then, the time-reversal invariant Kitaev chain belongs to the
Altland-Zirnbauer class BDI [4], which, in one spatial dimension, is
characterized by an integer-valued winding number [5, 6, 7]. The ground state
of the Kitaev chain is topologically trivial with a vanishing winding number
for $\left\lvert\mu\right\rvert>2J$, and topologically nontrivial, as
indicated by a nonvanishing value of the winding number, for
$\left\lvert\mu\right\rvert<2J$.

In what follows, an important role is played by inversion symmetry (IS) and
its breaking due to Markovian drive and dissipation. For the isolated Kitaev
chain, inversion is represented by a unitary operator $I$ that is defined as
follows [8]:

open boundary conditions: $\displaystyle Ic_{l}I^{\dagger}$
$\displaystyle=ic_{L+1-l},$ (59) periodic boundary conditions: $\displaystyle
Ic_{l}I^{\dagger}$ $\displaystyle=ic_{-l}.$

The Hamiltonian in Eq. (1) is symmetric under inversion, $H=IHI^{\dagger}$. In
the transformation of the Hamiltonian, the factor $i$ on the right-hand side
of Eq. (59) cancels in the hopping and chemical potential terms, but is
required to leave the pairing term invariant. For periodic boundary
conditions, the spinors $C_{k}$ in Eq. (29) transform as

$IC_{k}I^{\dagger}=i\sigma_{z}C_{-k}.$ (60)

Applying this transformation to the momentum-space representation of the
Hamiltonian in Eq. (30) yields, for the BdG Hamiltonian Eq. (31), the
following conditions of IS and, combining IS with TRS, parity-time symmetry
(PTS):

IS: $\displaystyle h_{k}$
$\displaystyle=\sigma_{z}h_{-k}\sigma_{z}=\sigma_{z}h_{-k}^{\dagger}\sigma_{z},$
(61) PTS: $\displaystyle h_{k}$
$\displaystyle=\sigma_{z}h_{k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{z}=\sigma_{z}h_{k}^{*}\sigma_{z}.$
(62)

For future reference, we list symmetry properties that can be derived from the
symmetries of the BdG Hamiltonian. First, for $a_{k}$ defined in Eq. (55), we
find

PHS: $\displaystyle a_{k}$
$\displaystyle=-a_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}=-a_{-k}^{*},$
(63) TRS: $\displaystyle a_{k}$
$\displaystyle=\sigma_{z}a_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{z}=\sigma_{z}a_{-k}^{*}\sigma_{z},$
(64) CS: $\displaystyle a_{k}$
$\displaystyle=-\sigma_{z}a_{k}\sigma_{z}=-\sigma_{z}a_{k}^{\dagger}\sigma_{z},$
(65) IS: $\displaystyle a_{k}$
$\displaystyle=\sigma_{y}a_{-k}\sigma_{y}=\sigma_{y}a_{-k}^{\dagger}\sigma_{y},$
(66) PTS: $\displaystyle a_{k}$
$\displaystyle=\sigma_{x}a_{k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{x}=\sigma_{x}a_{k}^{*}\sigma_{x}.$
(67)

Further, PHS yields restrictions on the transformation $U_{x,\theta_{k}}$ in
Eq. (36) that diagonalizes the BdG Hamiltonian. By combining Eqs. (37) and
(56) we obtain

$\sigma_{x}U_{x,\theta_{k}}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{x}h_{-k}\sigma_{x}U_{x,\theta_{k}}^{*}\sigma_{x}=\varepsilon_{k}\sigma_{z}.$
(68)

That is, we can identify $\varepsilon_{k}=\varepsilon_{-k}$ and

$U_{x,\theta_{k}}=\sigma_{x}U_{x,\theta_{-k}}^{*}\sigma_{x},$ (69)

which implies $\theta_{k}=-\theta_{-k}$. Alternatively, $\theta_{k}$ being an
odd function of $k$ can also be seen as a consequence of IS Eq. (61), which
implies

$U_{x,\theta_{k}}=\sigma_{z}U_{x,\theta_{-k}}\sigma_{z}.$ (70)

Moreover, in analogy to the PHS of the BdG Hamiltonian $h_{k}$ and its
counterpart $a_{k}$ in the Majorana basis, one can check that also the
covariance matrices $g_{k}$ and $\gamma_{k}$, defined in Eqs. (47) and (50),
respectively, obey PHS,

$g_{k}=-\sigma_{x}g_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{x},\qquad\gamma_{k}=-\gamma_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}.$
(71)

In particular, when we expand the covariance matrix in the basis of Pauli
matrices as

$\gamma_{k}=\gamma_{\mathbbm{1},k}\mathbbm{1}+\boldsymbol{\gamma}_{k}\cdot\boldsymbol{\sigma},$
(72)

PHS implies

$\gamma_{\mathbbm{1},k}=-\gamma_{\mathbbm{1},-k},$ (73)

and

$\gamma_{x,k}=-\gamma_{x,-k},\quad\gamma_{y,k}=\gamma_{y,-k},\quad\gamma_{z,k}=-\gamma_{z,-k}.$
(74)

Finally, Eq. (60) yields, for an inversion-symmetric state $\rho=I\rho
I^{\dagger}$ in in Eqs. (47) and (52), the following conditions for the
correlation and covariance matrices:

$g_{k}=\sigma_{z}g_{-k}\sigma_{z},\qquad\gamma_{k}=\sigma_{y}\gamma_{-k}\sigma_{y}.$
(75)

###  I.5 Relation to the transverse field Ising model

Through the Jordan-Wigner transformation, the isolated Kitaev chain can be
mapped to the transverse field Ising model, which serves as a paradigmatic
model system in studies of quench dynamics and generalized thermalization [9,
10, 11]. To establish the connection between our work and the existing
literature on quenches in the transverse field Ising model, we provide
relevant details of this mapping in the following. In Sec. III.3, we exploit
this connection to formulate conjectures for the time dependence of the
subsystem parity.

Usually, the Jordan-Wigner transformation is applied directly to the
Hamiltonian of the transverse field Ising model, which is given by

$H_{\mathrm{Ising}}=-J\sum_{l=1}^{L}\left(\sigma_{l}^{x}\sigma_{l+1}^{x}+h\sigma_{l}^{z}\right).$
(76)

Spin operators can be expressed in terms of fermionic operators through the
Jordan-Wigner transformation [12]:

$\begin{split}\sigma_{l}^{+}&=e^{i\pi\sum_{l^{\prime}=1}^{l-1}c_{l^{\prime}}^{\dagger}c_{l^{\prime}}}c_{l},\\\
\sigma_{l}^{-}&=e^{i\pi\sum_{l^{\prime}=1}^{l-1}c_{l^{\prime}}^{\dagger}c_{l^{\prime}}}c_{l}^{\dagger},\\\
\sigma_{l}^{z}&=e^{i\pi c_{l}^{\dagger}c_{l}}.\end{split}$ (77)

By means of the Jordan-Wigner transformation, the Hamiltonian of the
transverse field Ising model is mapped to two copies of the Kitaev chain Eq.
(1) in decoupled sectors with even and odd fermion parity, and $J=\Delta$ and
$\mu=-2Jh$ [10]. In particular, when we restrict ourselves to positive values
of $h$ (and, therefore, negative values of $\mu$), the ferromagnetically
ordered phase of the Ising model with $0<h<1$ maps to the topological phase of
the Kitaev chain with $\mu<-2J$, and the paramagnetic phase with $h>1$ maps to
the trivial phase with $-2J<\mu<0$. Long-range order in the ground state in
the ferromagnetic phase is revealed as usual by adding a longitudinal field
$Jh_{x}\sum_{l=1}^{L}\sigma_{l}^{x}$ to the Hamiltonian Eq. (76) and
considering the following combination of limits:

$\lim_{\left\lvert l-l^{\prime}\right\rvert\to\infty}\lim_{h_{x}\to
0}\lim_{L\to\infty}\langle\psi_{0}\rvert
O_{l,l^{\prime}}\lvert\psi_{0}\rangle\neq 0,$ (78)

where, for $l^{\prime}>l$,

$O_{l,l^{\prime}}=\sigma_{l}^{x}\sigma_{l^{\prime}}^{x}=\left(c_{l}+c_{l}^{\dagger}\right)e^{i\pi\sum_{m=l}^{l^{\prime}-1}c_{m}^{\dagger}c_{m}}\left(c_{l^{\prime}}+c_{l^{\prime}}^{\dagger}\right).$
(79)

That is, longitudinal or order parameter correlations map to the string order
parameter $O_{l,l^{\prime}}$ [13] for the Kitaev chain.

We find it convenient to deviate from the route outlined above and employ the
Kramers-Wannier duality [14, 15] of the transverse field Ising model before
applying the Jordan-Wigner transformation. The Kramers-Wannier self-duality of
the Ising model is established by defining new spin operators through the
relations

$\tau^{z}_{l}=\sigma^{x}_{l}\sigma^{x}_{l+1},\qquad\tau^{x}_{l}\tau^{x}_{l+1}=\sigma^{z}_{l}.$
(80)

Then, the Hamiltonian maps onto itself with the global energy scale
transforming as $J\to Jh$, and the transverse field as $h\to 1/h$, i.e., the
ferromagnetic phase of $\sigma$ spins maps to the paramagnetic phase of $\tau$
spins. Consequently, the ferromagnetic and paramagnetic phases of the Ising
model for $\tau$ spins map to the trivial and topological phases of the Kitaev
chain, respectively. In particular, the topologically trivial vacuum state,
which we choose as the initial state for our studies of quench dynamics,
corresponds to the ground state of the Ising model in the strongly
ferromagnetic limit. Longitudinal correlations of $\tau$ spins,

$P_{l,l^{\prime}}=\tau_{l}^{x}\tau_{l^{\prime}+1}^{x}=\prod_{m=l}^{l^{\prime}}\tau_{m}^{x}\tau_{m+1}^{x}=\prod_{m=l}^{l^{\prime}}\sigma_{m}^{z},$
(81)

serve as a “disorder parameter” for $\sigma$ spins, and the expectation value
of $P_{l,l^{\prime}}$ is maximized in the paramagnetic phase of $\sigma$
spins. Through the Jordan-Wigner transformation Eq. (77), $P_{l,l^{\prime}}$
maps to the fermion parity of a subsystem that consists of the lattice sites
$\\{l,\dotsc,l^{\prime}\\}$,

$P_{l,l^{\prime}}=e^{i\pi\sum_{m=l}^{l^{\prime}}c_{m}^{\dagger}c_{m}}=\prod_{m=l}^{l^{\prime}}\left(iw_{2m-1}w_{2m}\right).$
(82)

When the expectation value of $P_{l,l^{\prime}}$ is taken in a translationally
invariant state, the result depends only on the difference $l^{\prime}-l$.
Then, without loss of generality, we set $l=1$ and $l^{\prime}=\ell$, and
define the subsystem parity as

$P_{\ell}=e^{i\pi\sum_{l=1}^{\ell}c_{l}^{\dagger}c_{l}}=\prod_{l=1}^{\ell}\left(iw_{2l-1}w_{2l}\right).$
(83)

As stated in the main text, the expectation value of the subsystem parity is
given by the Pfaffian of the reduced covariance matrix [16, 17,
*Barouch1971II, *Barouch1971III],

$\langle P_{\ell}\rangle=\mathop{\mathrm{Pf}}(\Gamma_{\ell}),$ (84)

where

$\Gamma_{\ell}=\left(\Gamma_{l,l^{\prime}}\right)_{l,l^{\prime}=1}^{2\ell}.$
(85)

According to the above discussion, the ground-state expectation value $\langle
P_{\ell}\rangle_{0}=\langle\psi_{0}|P_{\ell}|\psi_{0}\rangle$ can be regarded
as a “topological disorder parameter:” The ground-state expectation value
$\langle P_{\ell}\rangle_{0}$ vanishes in the topological phase of the Kitaev
chain and is nonzero in the trivial phase. Accordingly, the results of Refs.
[9, 10, 11] for the time evolution of order parameter correlations for
quenches from the ordered to the ordered and disordered phases of the
transverse field Ising model, when interpreted as applying to $\tau$ spins,
describe the behavior of $\langle P_{\ell}\rangle$ for quenches from the
trivial to the trivial and topological phases of the Kitaev chain,
respectively.

##  II Driven-dissipative Kitaev chain

We move on to consider quantum quenches in the driven-dissipative Kitaev
chain. As specified in the main text, we assume that the system is prepared in
the ground state $\lvert\psi_{0}\rangle$ of the Hamiltonian of the isolated
Kitaev chain in Eq. (1) with parameters $J=\Delta$ and $\mu_{0}$. At $t=0$,
the chemical potential is quenched to the value $\mu$, and, concurrently, the
system is coupled to Markovian reservoirs.

###  II.1 Master equation and Liouvillian

The time evolution of the driven-dissipative Kitaev chain after the quench is
described by a quantum master equation in Lindblad form [20, 21],

$i\frac{d}{dt}\rho=\mathcal{L}\rho,$ (86)

where we define the Liouvillian $\mathcal{L}$ as the sum of the Hamiltonian
contribution $\mathcal{H}$ and the dissipator $\mathcal{D}$,

$\mathcal{L}=\mathcal{H}+i\mathcal{D},$ (87)

with

$\mathcal{H}\rho=[H,\rho],$ (88)

and

$\mathcal{D}\rho=\sum_{l=1}^{L}\left(2L_{l}\rho
L_{l}^{\dagger}-\\{L_{l}^{\dagger}L_{l},\rho\\}\right).$ (89)

We consider local jump operators as given in the main text,

$L_{l}=\sqrt{\gamma_{l}}c_{l}+\sqrt{\gamma_{g}}c_{l}^{\dagger}.$ (90)

In terms of Majorana operators $w_{l}$, the jump operators can be represented
as

$L_{l}=\sum_{l^{\prime}=1}^{2L}B_{l,l^{\prime}}w_{l^{\prime}},$ (91)

where the coefficients $B_{l,l^{\prime}}$ form an $L\times 2L$ matrix. For
future reference, we define the bath matrix $M$ as the product

$M=B^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}B^{*}.$
(92)

According to its definition, the bath matrix is Hermitian; Hermiticity implies
the following symmetry conditions for the real and imaginary parts of $M$,
which we denote by $M_{R}$ and $M_{I}$, respectively:

$M_{R}=M_{R}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},\qquad
M_{I}=-M_{I}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}.$
(93)

In turn, these relations imply

$M+M^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}=2M_{R},\qquad
M-M^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}=i2M_{I}.$
(94)

Like the matrix $A$ in Eq. (18), the bath matrix $M$ is a block Toeplitz
matrix, with $2\times 2$ blocks $m_{l}$ given by

$m_{l}=m_{R,l}+im_{I,l}.$ (95)

The real and imaginary parts, $m_{R,l}$ and $m_{I,l}$, respectively, read

$m_{R,l}=\frac{\delta_{l,0}}{2}\left[\gamma\mathbbm{1}+\sqrt{\gamma_{l}\gamma_{g}}\sigma_{z}\right],\qquad
m_{I,l}=\frac{\delta_{l,0}}{4}\delta i\sigma_{y}.$ (96)

Here, as in the main text, we introduce the mean rate $\gamma$, which
represents a measure of the overall strength of dissipation, and the
difference of loss and gain rates $\delta$,

$\gamma=\frac{\gamma_{l}+\gamma_{g}}{2},\qquad\delta=\gamma_{l}-\gamma_{g}.$
(97)

###  II.2 Time evolution of the covariance matrix

As stated in Sec. I.2, for the quench protocol specified in the main text, the
system is in a Gaussian state at all times, and Gaussian states are fully
determined by the covariance matrix. Therefore, in the following, we specify
the equation of motion for the covariance matrix and its numerical and
analytical solution.

####  II.2.1 Numerical time evolution

The equation of motion for the covariance matrix reads

$\frac{d\Gamma}{dt}=-X\Gamma-\Gamma
X^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}-Y,$
(98)

where

$X=4\left(iH+M_{R}\right)=-A+4M_{R},\quad Y=-8M_{I},$ (99)

and its formal solution is given by [22]

$\Gamma(t)=\Gamma_{1}(t)+\Gamma_{2}(t),$ (100)

with

$\begin{split}\Gamma_{1}(t)&=e^{-Xt}\Gamma(0)e^{-X^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}t},\\\
\Gamma_{2}(t)&=-\int_{0}^{t}dt^{\prime}\,e^{-X\left(t-t^{\prime}\right)}Ye^{-X^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\left(t-t^{\prime}\right)}.\end{split}$
(101)

The initial conditions and the steady state are encoded, respectively, in
$\Gamma_{1}(t)$ and $\Gamma_{2}(t)$. In particular, $\Gamma_{1}(0)=\Gamma(0)$
while $\Gamma_{2}(0)=0$, and $\Gamma_{1}(t)\to 0$ while
$\Gamma_{2}(t)\to\Gamma_{\mathrm{SS}}$ for $t\to\infty$. For all quenches
considered in the main text, the initial state is chosen as the ground state
of the Kitaev chain for $\mu_{0}\to-\infty$, i.e., the vacuum state, and the
corresponding covariance matrix $\Gamma(0)=\Gamma_{0}$ is given in Eq. (27).
The integral in the expression for $\Gamma_{2}(t)$ can be performed explicitly
if we express the matrix $X$ as $X=V\Lambda V^{-1}$, where $\Lambda$ is a
diagonal matrix with the eigenvalues $\lambda_{l}$ of $X$ on its diagonal.
Then,

$\Gamma_{2}(t)=-V\left[\left(V^{-1}YV^{-{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\right)\circ
K(t)\right]V^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(102)

where we use the shorthand notation
$V^{-{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}$
to denote the inverse of the transpose, and where $A\circ B$ is the Hadamard
product (i.e., element-wise) product of $A$ and $B$,

$\left(A\circ B\right)_{l,l^{\prime}}=A_{l,l^{\prime}}B_{l,l^{\prime}}.$ (103)

Further, the matrix $K(t)$ is defined as

$K_{l,l^{\prime}}(t)=\frac{1-e^{-\left(\lambda_{l}+\lambda_{l^{\prime}}\right)t}}{\lambda_{l}+\lambda_{l^{\prime}}}.$
(104)

Numerically, we find that $\lambda_{l}+\lambda_{l^{\prime}}\neq 0$ for all
$l,l^{\prime}$ when $\gamma>0$. To obtain the numerical data shown in the
plots in the main text, we implement Eqs. (101) and (102).

####  II.2.2 Analytical time evolution for periodic boundary conditions

For periodic boundary conditions, due to translational invariance of the
Hamiltonian, the coupling to Markovian baths, and the initial state, all
matrices in Eq. (98) are block-circulant Toeplitz matrices. Therefore, Eq.
(98) can be rewritten as

$\frac{d\gamma_{l}}{dt}=-\sum_{l^{\prime}=1}^{L}\left(x_{l-l^{\prime}}\gamma_{l^{\prime}}+\gamma_{l-l^{\prime}}x_{-l^{\prime}}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\right)-y_{l},$
(105)

where the $2\times 2$ blocks $\gamma_{l}$ are defined in Eq. (21), and $x_{l}$
and $y_{l}$ are defined analogously in terms of the matrices $X$ and $Y$ in
Eq. (99), i.e.,

$x_{l}=-a_{l}+4m_{R,l},\qquad y_{l}=-8m_{I,l},$ (106)

with $a_{l}$ and the real and imaginary parts of the blocks of the bath matrix
given in Eqs. (12) and (96), respectively. In momentum space, the equation of
motion for the Fourier transform $\gamma_{k}$ in Eq. (50) takes the form

$\frac{d\gamma_{k}}{dt}=-i\left(x_{k}\gamma_{k}-\gamma_{k}x_{k}^{\dagger}\right)-y_{k},$
(107)

where we define the Fourier transforms of $x_{l}$ and $y_{l}$ with a prefactor
$-i$ as

$\begin{split}x_{k}&=-i\sum_{l=1}^{L}e^{-ikl}x_{l}=x_{\mathbbm{1}}\mathbbm{1}+\mathbf{x}_{k}\cdot\boldsymbol{\sigma},\\\
y_{k}&=-i\sum_{l=1}^{L}e^{-ikl}y_{l}=\mathbf{y}\cdot\boldsymbol{\sigma}.\end{split}$
(108)

The coefficients in the expansions of $x_{k}$ and $y_{k}$ in terms of Pauli
matrices read

$x_{\mathbbm{1}}=-i2\gamma,\quad\mathbf{x}_{k}=\mathbf{a}_{k}+i\mathbf{b},\quad\mathbf{b}=-2\sqrt{\gamma_{l}\gamma_{g}}\hat{\mathbf{e}}_{z},$
(109)

where $\mathbf{a}_{k}$ is given in Eq. (54), and

$\mathbf{y}=-2\delta\hat{\mathbf{e}}_{y}.$ (110)

In the above expressions, we denote unit vectors in the $x$, $y$, and $z$
direction by $\hat{\mathbf{e}}_{x,y,z}$. The solution to Eq. (107) is given by

$\gamma_{k}(t)=\gamma_{1,k}(t)+\gamma_{2,k}(t),$ (111)

where

$\begin{split}\gamma_{1,k}(t)&=e^{-ix_{k}t}\gamma_{k}(0)e^{ix_{k}^{\dagger}t},\\\
\gamma_{2,k}(t)&=-\int_{0}^{t}dt^{\prime}\,e^{-ix_{k}\left(t-t^{\prime}\right)}y_{k}e^{ix_{k}^{\dagger}\left(t-t^{\prime}\right)},\end{split}$
(112)

and the initial conditions are determined by Eq. (53),
$\gamma_{k}(0)=\hat{\mathbf{a}}_{0,k}\cdot\boldsymbol{\sigma}$. In particular,
for $\mu_{0}\to-\infty$, we find $\gamma_{k}(0)=\sigma_{y}$.

####  II.2.3 Spectrum of the covariance matrix

For future reference, we briefly discuss the spectrum of the time-evolved
covariance matrix. To that end, it is convenient to define the coefficients
$\gamma_{\mathbbm{1},k}(t)$ and $\boldsymbol{\gamma}_{k}(t)$ in an expansion
in terms of Pauli matrices through the relation

$\gamma_{k}(t)=\gamma_{\mathbbm{1},k}(t)\mathbbm{1}+\boldsymbol{\gamma}_{k}(t)\cdot\boldsymbol{\sigma}.$
(113)

The eigenvalues of $\boldsymbol{\gamma}_{k}\cdot\boldsymbol{\sigma}$ are
$\pm\xi_{k}=\pm\lvert\boldsymbol{\gamma}_{k}\rvert$, where we omit the
explicit time dependence to lighten the notation. Therefore, the spectrum
$\sigma(\gamma_{k})$ of the covariance matrix in Eq. (113) is given by

$\sigma(\gamma_{k})=\\{\gamma_{\mathbbm{1},k}\pm\xi_{k}\\}.$ (114)

We note that for an isolated system with unitary dynamics, time evolution
corresponds to precession of the vector $\boldsymbol{\gamma}_{k}(t)$ around
$\mathbf{a}_{k}$, with the explicit form of $\boldsymbol{\gamma}_{k}(t)$ given
in Eq. (164) below. Then,

$\gamma_{\mathbbm{1},k}=0,\quad\xi_{k}=\left\lvert\hat{\mathbf{a}}_{0,k}\right\rvert=1,$
(115)

and the spectrum is symmetric with respect to zero. This is, in general, not
the case for an open system. However, there are restrictions on the spectrum
due to PHS, which imply that the union of $\sigma(\gamma_{k})$ and
$\sigma(\gamma_{-k})$ is again symmetric with respect to zero. This can be
seen as follows: Equation (74) implies that $\xi_{k}=\xi_{-k}$. Then, with Eq.
(73) we obtain

$\sigma(\gamma_{-k})=\\{-\gamma_{\mathbbm{1},k}\pm\xi_{k}\\}.$ (116)

Therefore, when we define

$\xi_{+,k}=\gamma_{\mathbbm{1},k}+\xi_{k},\quad\xi_{-,k}=-\gamma_{\mathbbm{1},k}+\xi_{k},$
(117)

we can write

$\sigma(\gamma_{k})\cup\sigma(\gamma_{-k})=\\{\pm\xi_{+,k},\pm\xi_{-,k}\\}.$
(118)

###  II.3 Time evolution of the complex covariance matrix

The equation of motion for the covariance matrix in the basis of complex
fermions, which is given in the main text, can be obtained by transforming Eq.
(107) according to Eq. (50),

$\frac{dg_{k}}{dt}=-i\left(z_{k}g_{k}-g_{k}z_{k}^{\dagger}\right)-s_{k},$
(119)

where

$z_{k}=rx_{k}r^{\dagger}=z_{\mathbbm{1}}\mathbbm{1}+\mathbf{z}_{k}\cdot\boldsymbol{\sigma},\qquad
s_{k}=ry_{k}r^{\dagger}=\mathbf{s}\cdot\boldsymbol{\sigma}.$ (120)

By using Eqs. (15) and (16), we obtain

$z_{\mathbbm{1}}=-i2\gamma,\quad\mathbf{z}_{k}=\mathbf{h}_{k}-i2\sqrt{\gamma_{l}\gamma_{g}}\hat{\mathbf{e}}_{x},\quad\mathbf{s}=-2\delta\hat{\mathbf{e}}_{z}.$
(121)

Rather than solving Eq. (119) explicitly, we rotate Eq. (113) to the complex
basis to obtain the time-evolved covariance matrix in the form

$g_{k}(t)=r\gamma_{k}(t)r^{\dagger}=g_{\mathbbm{1},k}(t)\mathbbm{1}+\mathbf{g}_{k}(t)\cdot\boldsymbol{\sigma},$
(122)

where

$g_{\mathbbm{1},k}(t)=\gamma_{\mathbbm{1},k}(t),\quad\mathbf{g}_{k}(t)=P^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\boldsymbol{\gamma}_{k}(t).$
(123)

###  II.4 Symmetries of the Liouvillian

Before we go into details of quench dynamics which are encoded in the time
dependence of the covariance matrix, we discuss relevant symmetries of the
Liouvillian $\mathcal{L}$ that generates the dynamics. According to Eqs. (109)
and (121), the matrices $x_{k}$ and $z_{k}$ can be regarded as generalizations
of the matrices $a_{k}$ and the BdG Hamiltonian $h_{k}$, respectively, to open
systems. Therefore, by considering $x_{k}$ and $z_{k}$, we can check directly
how the symmetries of the isolated Kitaev chain, which are summarized in Sec.
I.4, are affected by dissipation. We find that of the two versions of, e.g.,
PHS in Eq. (56) for $\gamma=0$, only one remains valid when $\gamma>0$.
Further, IS and, therefore, PTS, apply only to the traceless parts of $x_{k}$
and $z_{k}$, given by $x_{k}^{\prime}=x_{k}+i2\gamma\mathbbm{1}$ and
$z_{k}^{\prime}=z_{k}+i2\gamma\mathbbm{1}$, respectively. This is known as
passive PTS [23, 24, 25]. In particular, using the terminology of Ref. [26],
the symmetries of $x_{k}$ are given by

PHS†: $\displaystyle x_{k}$ $\displaystyle=-x_{-k}^{*},$ (124) TRS†:
$\displaystyle x_{k}$
$\displaystyle=\sigma_{z}x_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{z},$
(125) CS: $\displaystyle x_{k}$
$\displaystyle=-\sigma_{z}x_{k}^{\dagger}\sigma_{z},$ (126) IS†:
$\displaystyle x_{k}^{\prime}$
$\displaystyle=\sigma_{y}x_{-k}^{\prime\dagger}\sigma_{y},$ (127) PTS:
$\displaystyle x_{k}^{\prime}$
$\displaystyle=\sigma_{x}x_{k}^{\prime*}\sigma_{x},$ (128)

where PTS is the combination of TRS† and IS†. For the sake of completeness, we
also list the symmetries of $z_{k}$, which result from rotating the above
relations to the complex basis:

PHS†: $\displaystyle z_{k}$ $\displaystyle=-\sigma_{x}z_{-k}^{*}\sigma_{x},$
(129) TRS†: $\displaystyle z_{k}$
$\displaystyle=z_{-k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(130) CS: $\displaystyle z_{k}$
$\displaystyle=-\sigma_{x}z_{k}^{\dagger}\sigma_{x},$ (131) IS†:
$\displaystyle z_{k}^{\prime}$
$\displaystyle=\sigma_{z}z_{-k}^{\prime\dagger}\sigma_{z},$ (132) PTS:
$\displaystyle z_{k}^{\prime}$
$\displaystyle=\sigma_{z}z_{k}^{\prime*}\sigma_{z}.$ (133)

We next discuss implications of the above symmetry conditions that are
relevant to our studies of quench dynamics. First, we show that IS† is not
sufficient to ensure that the time-evolved covariance matrix obeys IS as
specified in Eq. (75). To that end, we apply inversion to the time-evolved
covariance matrix Eq. (112), whereby we take into account that the initial
state in Eq. (53), as well as $y_{k}$ defined Eq. (108), obey inversion
symmetry in the sense of Eq. (75),

$\begin{split}\sigma_{y}\gamma_{1,-k}(t)\sigma_{y}&=e^{-i2x_{\mathbbm{1}}t}e^{-i\mathbf{x}_{k}^{\dagger}\cdot\boldsymbol{\sigma}t}\gamma_{k}(0)e^{i\mathbf{x}_{k}\cdot\boldsymbol{\sigma}t},\\\
\sigma_{y}\gamma_{2,-k}(t)\sigma_{y}&=-\int_{0}^{t}dt^{\prime}\,e^{-i2x_{\mathbbm{1}}\left(t-t^{\prime}\right)}e^{-i\mathbf{x}_{k}^{\dagger}\cdot\boldsymbol{\sigma}\left(t-t^{\prime}\right)}y_{k}e^{i\mathbf{x}_{k}\cdot\boldsymbol{\sigma}\left(t-t^{\prime}\right)}.\end{split}$
(134)

By comparing these expressions to Eq. (112), we see that inversion amounts to
reversing the sign of the vector $\mathbf{b}$ in Eq. (109), which, in general,
does not leave the state invariant. For quenches to the PT-symmetric phase,
this can be checked explicitly in Eqs. (157), (158), (161), and (162) below.

Further, as stated in the main text, and as we explain in detail in the
following, PTS has important consequences for the spectrum of $z_{k}^{\prime}$
[27, *Bender2008]. Since
$z_{k}^{\prime}=z_{k}+i2\gamma\mathbbm{1}=\mathbf{z}_{k}\cdot\boldsymbol{\sigma}$,
the eigenvalues of $z_{k}^{\prime}$ are given by
$\lambda_{\pm,k}^{\prime}=\pm\sqrt{\mathbf{z}_{k}\cdot\mathbf{z}_{k}}$. We
denote the corresponding eigenvectors, which are also eigenvectors of $z_{k}$
with eigenvalues $\lambda_{\pm,k}=\lambda_{\pm,k}^{\prime}-i2\gamma$, by
$\lvert\psi_{\pm,k}\rangle$. Then, the condition of PTS leads to

$z_{k}^{\prime}\sigma_{z}\lvert\psi_{+,k}\rangle^{*}=\lambda_{+,k}^{\prime*}\sigma_{z}\lvert\psi_{+,k}\rangle^{*}.$
(135)

That is, $\sigma_{z}\lvert\psi_{+,k}\rangle^{*}$ is an eigenvector of
$z_{k}^{\prime}$, and the corresponding eigenvalue is given by
$\lambda_{+,k}^{\prime*}$. Now, there are two possibilities: (i)
$\lambda_{+,k}^{\prime*}=\lambda_{+,k}^{\prime}\in\mathbb{R}$. Then,
$\sigma_{z}\lvert\psi_{+,k}\rangle^{*}=\lvert\psi_{+,k}\rangle$ (potentially,
up to a phase that can be absorbed in a redefiniton of
$\lvert\psi_{+,k}\rangle$). Further, also
$\lambda_{-,k}^{\prime}=-\lambda_{+,k}^{\prime}\in\mathbb{R}$, and
$\sigma_{z}\lvert\psi_{-,k}\rangle^{*}=\lvert\psi_{-,k}\rangle$. For the
eigenvalues of $z_{k}$, this implies
$\mathop{\mathrm{Re}}(\lambda_{+,k})=-\mathop{\mathrm{Re}}(\lambda_{-,k})$ and
$\mathop{\mathrm{Im}}(\lambda_{\pm,k})=-2\gamma$. (ii)
$\lambda_{+,k}^{\prime*}=-\lambda_{+,k}^{\prime}=\lambda_{-,k}^{\prime}\in
i\mathbb{R}$, and
$\sigma_{z}\lvert\psi_{\pm,k}\rangle^{*}=\lvert\psi_{\mp,k}\rangle$. For the
eigenvalues $z_{k}$, this implies $\mathop{\mathrm{Re}}(\lambda_{\pm,k})=0$
and
$\mathop{\mathrm{Im}}(\lambda_{+,k})=-\mathop{\mathrm{Im}}(\lambda_{-,k})-4\gamma$.
We are, therefore, lead to distinguish three phases: In the PT-symmetric
phase, (i) applies to all values of $k\in\mathrm{BZ}$, while (ii) applies to
all $k$ in the PT-broken phase. Finally, in the PT-mixed phase, both PT-
symmetric and PT-breaking momentum modes exist.

Before we study how these phases are realized in the driven-dissipative Kitaev
chain, we briefly describe the form of PTS in position space, which is
relevant, in particular, for the analysis of systems with open boundary
conditions. By taking the inverse Fourier transform of Eqs. (125), (127), and
(128) according to Eq. (108), we obtain

TRS†: $\displaystyle x_{l}$
$\displaystyle=\sigma_{z}x_{-l}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{z},$
(136) IS†: $\displaystyle x_{l}^{\prime}$
$\displaystyle=-\sigma_{y}x_{l}^{\prime{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\sigma_{y},$
(137) PTS: $\displaystyle x_{l}^{\prime}$
$\displaystyle=-\sigma_{x}x_{-l}^{\prime}\sigma_{x},$ (138)

where we define $x_{l}^{\prime}=x_{l}-2\gamma\mathbbm{1}$. For the matrix $X$,
which is constructed from the $2\times 2$ blocks $x_{l}$ like the matrix $A$
in Eq. (18) is constructed from $a_{l}$, and the shifted matrix
$X^{\prime}=X-2\gamma\mathbbm{1}$, we find

TRS†: $\displaystyle X$
$\displaystyle=\Sigma_{z,L}X^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\Sigma_{z,L},$
(139) IS†: $\displaystyle X^{\prime}$
$\displaystyle=-\Sigma_{y,L}\Xi_{L}X^{\prime{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\Xi_{L}\Sigma_{y,L},$
(140) PTS: $\displaystyle X^{\prime}$
$\displaystyle=-\Sigma_{x,L}\Xi_{L}X^{\prime}\Xi_{L}\Sigma_{x,L}.$ (141)

where, for $\mu\in\\{x,y,z\\}$,

$\Sigma_{\mu,\ell}=\bigoplus_{l=1}^{\ell}\sigma_{\mu},$ (142)

and $\Xi_{\ell}$ is a $2\ell\times 2\ell$ block anti-diagonal matrix, with
$2\times 2$ identity matrices on the anti-diagonal,

$\Xi_{\ell}=\begin{pmatrix}&&\mathbbm{1}\\\ &\iddots&\\\
\mathbbm{1}&&\end{pmatrix}.$ (143)

According to our convention for the Fourier transform in Eq. (108), for a
system with periodic boundary conditions, the union of the spectra of the
matrices $x_{k}$ for all $k\in\mathrm{BZ}$ yields the spectrum of $-iX$.
Therefore, we focus on the consequences of PTS for the spectrum of the latter.
PTS in the form given in Eq. (141) implies that for each nonzero eigenvalue
$\lambda^{\prime}$ of $-iX^{\prime}=-i\left(X-2\gamma\mathbbm{1}\right)$,
which obeys the eigenvalue equation
$-iX^{\prime}\lvert\psi\rangle=\lambda^{\prime}\lvert\psi\rangle$, also
$-\lambda^{\prime}$ is an eigenvalue, and the corresponding eigenvector is
given by $\Sigma_{x,L}\Xi_{L}\lvert\psi\rangle$. Further, since $-iX^{\prime}$
is purely imaginary, its eigenvalues $\lambda^{\prime}$ are either purely
imaginary or come in pairs $\lambda^{\prime}$ and $-\lambda^{\prime*}$.
Therefore, for each eigenvalue $\lambda^{\prime}$, there are three
possibilities: (i) There is a pair of real eigenvalues given by
$\pm\lambda^{\prime}\in\mathbb{R}$. (ii) There is a pair of purely imaginary
eigenvalues given by $\pm\lambda^{\prime}\in i\mathbb{R}$. (iii) If
$\mathop{\mathrm{Re}}(\lambda^{\prime})\neq 0$ and
$\mathop{\mathrm{Im}}(\lambda^{\prime})\neq 0$, there are four eigenvalues
given by $\pm\lambda^{\prime}$ and $\pm\lambda^{\prime*}$. Numerically, we
find that only cases (i) and (ii) are realized for the driven-dissipative
Kitaev chain. That is, the eigenmodes of
$-iX=-i\left(X^{\prime}+2\gamma\mathbbm{1}\right)$ are either (i) PT-
symmetric, with eigenvalues
$\lambda=\lambda^{\prime}-i2\gamma\in\mathbb{R}-i2\gamma$, or (ii) PT-
breaking, with $\lambda\in i\mathbb{R}$.

###  II.5 Spectrum of the Liouvillian

As can be shown by employing the formalism of third quantization [29, 30], for
a system with periodic boundary conditions, the spectrum of the Liouvillian is
determined by the spectrum of the matrix $x_{k}$ in Eq. (108), or,
equivalently, the matrix $z_{k}$ defined in Eq. (120), similarly to the
spectrum of the Hamiltonian Eq. (30) that can be obtained by populating the
single-particle energy levels determined by the BdG Hamiltonian Eq. (31). The
eigenvalues of $z_{k}$ are given by

$\lambda_{\pm,k}=-i2\gamma\pm\omega_{k},$ (144)

where

$\omega_{k}=\sqrt{\varepsilon_{k}^{2}-4\gamma_{l}\gamma_{g}}.$ (145)

As stated in the main text, the PT-symmetric phase, in which
$\omega_{k}\in\mathbb{R}_{>0}$, is realized for

$2\sqrt{\gamma_{l}\gamma_{g}}<\left\lvert
2J-\left\lvert\mu\right\rvert\right\rvert.$ (146)

Then, the bands $\lambda_{\pm,k}$ are separated by a real line gap [26, 31].
We note that for a system with open boundary conditions, we find numerically
that the spectrum of $-iX$ obeys $\sigma(-iX)\in-i2\gamma+\mathbb{R}$ as
suggested by Eq. (144) only for $\left\lvert\mu\right\rvert>2J$. In contrast,
for $\left\lvert\mu\right\rvert<2J$ there are PT-breaking edge modes with
purely imaginary eigenvalues in addition to fully PT-symmetric bulk of the
system [31]. For $\delta=0$, the eigenvalues corresponding to the edge modes
are $0$ and $-i4\gamma$. Increasing the strength of dissipation leads to the
occurrence of PT-breaking modes also in the bulk. In particular, in the PT-
broken phase for

$2\sqrt{\gamma_{l}\gamma_{g}}>2J+\left\lvert\mu\right\rvert,$ (147)

we write $\omega_{k}=i\kappa_{k}\in i\mathbb{R}_{>0}$ with

$\kappa_{k}=\sqrt{4\gamma_{l}\gamma_{g}-\varepsilon_{k}^{2}},$ (148)

such that Eq. (144) becomes

$\lambda_{\pm,k}=-i\left(2\gamma\mp\kappa_{k}\right).$ (149)

In this phase, the bands $\lambda_{\pm,k}$ are separated by an imaginary line
gap. Finally, for intermediate strengths of dissipation, the spectrum is
gapless. In the gapless or PT-mixed phase, the bands $\lambda_{\pm,k}$ cross
at exceptional points at $k=\pm k_{*}$, at which $\omega_{k}$ defined in Eq.
(145) vanishes, and which separate PT-symmetric from PT-breaking modes.

For future reference, we note that $z_{k}$ can be diagonalized in the PT-
symmetric phase by combining an ordinary rotation with a hyperbolic rotation,
i.e., a rotation by an imaginary angle. In particular, with

$U_{k}=U_{x,\theta_{k}}U_{y,-i\beta_{k}}=e^{i\left(\theta_{k}/2\right)\sigma_{x}}e^{\left(\beta_{k}/2\right)\sigma_{y}},$
(150)

where

$\beta_{k}=-\mathop{\mathrm{arctanh}}(2\sqrt{\gamma_{l}\gamma_{g}}/\varepsilon_{k}),$
(151)

we obtain

$z_{k}=-i2\gamma\mathbbm{1}+\omega_{k}U_{k}\sigma_{z}U_{k}^{-1}.$ (152)

By rotating this relation to the real basis, we find that $x_{k}$ can be
diagonalized as (note the appearance of $\sigma_{y}$ instead of $\sigma_{z}$)

$x_{k}=-i2\gamma\mathbbm{1}+\omega_{k}S_{k}\sigma_{y}S_{k}^{-1},$ (153)

where

$S_{k}=r^{\dagger}U_{k}r=U_{z,\theta_{k}}U_{x,-i\beta_{k}}=e^{i\left(\theta_{k}/2\right)\sigma_{z}}e^{\left(\beta_{k}/2\right)\sigma_{x}}.$
(154)

The symmetries of $z_{k}$ impose restrictions on the transformation $U_{k}$.
Analogously to the derivation of Eq. (70), on can show that inversion symmetry
leads to the relation

$U_{k}=\sigma_{z}U_{-k}^{-\dagger}\sigma_{z},$ (155)

which implies $\theta_{k}=-\theta_{-k}$ and $\beta_{k}=\beta_{-k}$.

##  III Quench dynamics in the PT-symmetric phase and relaxation to the PTGGE

We now focus on quench dynamics in the PT-symmetric phase and, in particular,
on relaxation to the PT-symmetric generalized Gibbs ensemble (PTGGE). In the
main text, we derive the PTGGE, which describes the late-time dynamics in the
PT-symmetric phase, as the maximum-entropy ensemble that is compatible with
given dephased expectation values of commutators of eigenmodes of the
Liouvillian, and with the modified statistics of these modes as encoded in
their anticommutators. Below, we present details of this calculation. But
first, we present a derivation of the PTGGE that is based on the explicit form
of the time-evolved covariance matrix Eq. (111).

###  III.1 Dephasing and relaxation to the PTGGE

The most direct way to derive the PTGGE is to determine simplifications of the
time-evolved covariance matrix Eq. (111) in the late-time limit due to
dephasing. Technical details of the time evolution generated by the matrix
$x_{k}$, which can be interpreted as a non-Hermitian Hamiltonian, are
presented in the Appendix.

####  III.1.1 The covariance matrix in the PT-symmetric phase

First, we provide explicit expressions for the components $\gamma_{1,k}(t)$
and $\gamma_{2,k}(t)$ of the time-evolved covariance matrix defined in Eq.
(112). With the initial condition
$\gamma_{k}(0)=\hat{\mathbf{a}}_{0,k}\cdot\boldsymbol{\sigma}$ specified in
Eq. (53), we write $\gamma_{1,k}(t)$ in the form

$\gamma_{1,k}(t)=e^{-ix_{k}t}\hat{\mathbf{a}}_{0,k}\cdot\boldsymbol{\sigma}e^{ix_{k}^{\dagger}t}=\gamma_{1,\mathbbm{1},k}(t)\mathbbm{1}+\boldsymbol{\gamma}_{1,k}(t)\cdot\boldsymbol{\sigma}.$
(156)

Then, Eq. (336) leads to

$\gamma_{1,\mathbbm{1},k}(t)=-\frac{e^{-4\gamma
t}}{\omega_{k}^{2}}\left(1-\cos(2\omega_{k}t)\right)\hat{\mathbf{a}}_{0,k}\cdot\left(\mathbf{a}_{k}\times\mathbf{b}\right),$
(157)

and from Eq. (337) we obtain

$\boldsymbol{\gamma}_{1,k}(t)=e^{-4\gamma
t}\left[\hat{\mathbf{a}}_{0,k}-\frac{\varepsilon_{k}^{2}}{\omega_{k}^{2}}\left(1-\cos(2\omega_{k}t)\right)\mathbf{a}_{\perp,k}\right.\\\
\left.+\frac{\varepsilon_{k}}{\omega_{k}}\sin(2\omega_{k}t)\mathbf{a}_{o,k}\right],$
(158)

where, with $\hat{\mathbf{a}}_{k}=\mathbf{a}_{k}/\varepsilon_{k}$, the
parallel, perpendicular, and out-of-plane components are defined as

$\begin{split}\mathbf{a}_{\parallel,k}&=\left(\hat{\mathbf{a}}_{0,k}\cdot\hat{\mathbf{a}}_{k}\right)\hat{\mathbf{a}}_{k},\\\
\mathbf{a}_{\perp,k}&=\hat{\mathbf{a}}_{0,k}-\mathbf{a}_{\parallel,k},\\\
\mathbf{a}_{o,k}&=-\hat{\mathbf{a}}_{0,k}\times\hat{\mathbf{a}}_{k}.\end{split}$
(159)

Analogously, for

$\gamma_{2,k}(t)=\gamma_{2,\mathbbm{1},k}(t)\mathbbm{1}+\boldsymbol{\gamma}_{2,k}(t)\cdot\boldsymbol{\sigma},$
(160)

we find

$\gamma_{2,\mathbbm{1},k}(t)=\int_{0}^{t}dt^{\prime}\,\frac{e^{-4\gamma\left(t-t^{\prime}\right)}}{\omega_{k}^{2}}\\\
\times\left(1-\cos(2\omega_{k}\left(t-t^{\prime}\right))\right)\mathbf{y}\cdot\left(\mathbf{a}_{k}\times\mathbf{b}\right),$
(161)

and

$\boldsymbol{\gamma}_{2,k}(t)=-\int_{0}^{t}dt^{\prime}\,e^{-4\gamma\left(t-t^{\prime}\right)}\\\
\times\left[\mathbf{y}-\frac{\varepsilon_{k}^{2}}{\omega_{k}^{2}}\left(1-\cos(2\omega_{k}\left(t-t^{\prime}\right))\right)\mathbf{y}_{\perp,k}\right.\\\
\left.+\frac{\varepsilon_{k}}{\omega_{k}}\sin(2\omega_{k}\left(t-t^{\prime}\right))\mathbf{y}_{o,k}\right],$
(162)

where

$\begin{split}\mathbf{y}_{\parallel,k}&=\left(\mathbf{y}\cdot\hat{\mathbf{a}}_{k}\right)\hat{\mathbf{a}}_{k},\\\
\mathbf{y}_{\perp,k}&=\mathbf{y}-\mathbf{y}_{\parallel,k},\\\
\mathbf{y}_{o,k}&=-\mathbf{y}\times\hat{\mathbf{a}}_{k}.\end{split}$ (163)

The integration in Eqs. (161) and (162) is elementary and we omit the lengthy
results.

####  III.1.2 Late-time limit and dephasing

In isolated systems, the defining signature of generalized thermalization
after a quantum quench is local observables assuming stationary expectation
values that are determined by the GGE. For noninteracting and translationally
invariant systems, generalized thermalization occurs through dephasing of
momentum modes that oscillate at different frequencies
$\varepsilon_{k}\neq\varepsilon_{k^{\prime}}$ for $k\neq k^{\prime}$ [32, 33].
In particular, for the isolated Kitaev chain, the component of the time-
evolved covariance matrix with momentum $k$ can be obtained by setting
$\gamma=0$ in Eqs. (157), (158), (161), and (162), which yields
$\gamma_{1,\mathbbm{1},k}(t)=\gamma_{2,k}(t)=0$, such that
$\gamma_{k}(t)=\boldsymbol{\gamma}_{k}(t)\cdot\boldsymbol{\sigma}$ with

$\boldsymbol{\gamma}_{k}(t)=\mathbf{a}_{\parallel,k}+\cos(2\varepsilon_{k}t)\mathbf{a}_{\perp,k}+\sin(2\varepsilon_{k}t)\mathbf{a}_{o,k}.$
(164)

The long-time behavior that is established through dephasing is captured by
setting the oscillatory components to zero. Thus, with subscript “d” for
dephased,

$\boldsymbol{\gamma}_{d,k}=\mathbf{a}_{\parallel,k}=\left(\hat{\mathbf{a}}_{0,k}\cdot\hat{\mathbf{a}}_{k}\right)\hat{\mathbf{a}}_{k},$
(165)

and the generalized Gibbs ensemble (GGE) is the Gaussian state that is
uniquely determined by the dephased covariance matrix
$\gamma_{d,k}=\boldsymbol{\gamma}_{d,k}\cdot\boldsymbol{\sigma}$. Since
$\boldsymbol{\gamma}_{d,k}\parallel\hat{\mathbf{a}}_{k}$, the dephased
covariance matrix, and, therefore, the GGE, is diagonal in the eigenbasis of
the postquench Hamiltonian. The diagonal elements are determined by the
conserved mode occupation numbers of the eigenmodes $d_{k}$ in Eq. (38) of the
postquench Hamiltonian in the initial state $\lvert\psi_{0}\rangle$,

$\hat{\mathbf{a}}_{0,k}\cdot\hat{\mathbf{a}}_{k}=\langle\psi_{0}|[d_{k},d_{k}^{\dagger}]|\psi_{0}\rangle=1-2\langle\psi_{0}|d_{k}^{\dagger}d_{k}|\psi_{0}\rangle,$
(166)

as follows from setting $\gamma=0$ in Eq. (222) below. For quenches in the
driven-dissipative Kitaev chain, this picture is modified in two fundamental
ways:

(i) The time dependence of the covariance matrix is described by renormalized
oscillation frequencies $\omega_{k}$ given in Eq. (145) and, crucially, there
is additional exponential decay. However, due to PTS, there is only a single
relaxation rate, which, for the driven-dissipative Kitaev chain, is determined
by the imaginary part $-\mathop{\mathrm{Im}}(\lambda_{\pm,k})=2\gamma$ of the
eigenvalues of $z_{k}$ in Eq. (144). After factoring out the resulting trivial
exponential time dependence, the dispersiveness of the real part
$\mathop{\mathrm{Re}}(\lambda_{\pm,k})=\pm\omega_{k}$ induces dephasing and,
consequently, relaxation dynamics in analogy to isolated systems. In contrast,
the dynamics of open systems without PTS is typically governed by a range of
exponential relaxation rates, and the slowest of these rates determines the
late-time behavior, while dephasing does not play an important role. For the
PT-symmetric driven-dissipative Kitaev chain, as in Eq. (165), dephasing can
be accounted for by setting oscillatory contributions to zero. Therefore, the
late-time behavior is described by

$\gamma_{k}(t)\sim\gamma_{d,k}(t)=\gamma_{d,1,k}(t)+\gamma_{d,2,k}(t)+\gamma_{\mathrm{SS},k},$
(167)

with

$\gamma_{d,1,k}(t)=e^{-4\gamma
t}\gamma_{d,1,k}^{\prime},\quad\gamma_{d,1,k}^{\prime}=\gamma_{d,1,\mathbbm{1},k}^{\prime}\mathbbm{1}+\boldsymbol{\gamma}_{d,1,k}^{\prime}\cdot\boldsymbol{\sigma},$
(168)

where

$\gamma_{d,1,\mathbbm{1},k}^{\prime}=-\frac{1}{\omega_{k}^{2}}\hat{\mathbf{a}}_{0,k}\cdot\left(\mathbf{a}_{k}\times\mathbf{b}\right)=\frac{2\varepsilon_{k}\sqrt{\gamma_{l}\gamma_{g}}}{\omega_{k}^{2}}\sin(\Delta\theta_{k}),$
(169)

and

$\begin{split}\boldsymbol{\gamma}_{d,1,k}^{\prime}&=\hat{\mathbf{a}}_{0,k}-\frac{\varepsilon_{k}^{2}}{\omega_{k}^{2}}\mathbf{a}_{\perp,k}\\\
&=\cos(\Delta\theta_{k})\hat{\mathbf{a}}_{k}+\frac{4\gamma_{l}\gamma_{g}}{\omega_{k}^{2}}\sin(\Delta\theta_{k})\hat{\mathbf{e}}_{z}\times\hat{\mathbf{a}}_{k}.\end{split}$
(170)

In these expressions, $\Delta\theta_{k}$ is the angle between
$\hat{\mathbf{a}}_{0,k}$ and $\hat{\mathbf{a}}_{k}$,

$\cos(\Delta\theta_{k})=\hat{\mathbf{a}}_{0,k}\cdot\hat{\mathbf{a}}_{k}.$
(171)

Analogously,

$\gamma_{d,2,k}(t)=e^{-4\gamma t}\gamma_{d,2,k}^{\prime},$ (172)

where

$\gamma_{d,2,k}^{\prime}=-\frac{1}{4\gamma}\left[\frac{1}{\omega_{k}^{2}}\mathbf{y}\cdot\left(\mathbf{a}_{k}\times\mathbf{b}\right)\mathbbm{1}-\left(\mathbf{y}-\frac{\varepsilon_{k}^{2}}{\omega_{k}^{2}}\mathbf{y}_{\perp,k}\right)\cdot\boldsymbol{\sigma}\right].$
(173)

Finally,

$\gamma_{\mathrm{SS},k}=\frac{1}{4\gamma}\left\\{\frac{1}{\varepsilon_{k}^{2}+\delta^{2}}\mathbf{y}\cdot\left(\mathbf{a}_{k}\times\mathbf{b}\right)\mathbbm{1}\right.\\\
\left.-\left[\mathbf{y}-\frac{\varepsilon_{k}}{\varepsilon_{k}^{2}+\delta^{2}}\left(\varepsilon_{k}\mathbf{y}_{\perp,k}-2\gamma\mathbf{y}_{o,k}\right)\right]\cdot\boldsymbol{\sigma}\right).$
(174)

(ii) The contribution $\gamma_{2,k}(t)$ in Eq. (111), which describes the
approach to the steady state $\rho_{\mathrm{SS}}$ defined by
$\mathcal{L}(\rho_{\mathrm{SS}})=0$ and given explicitly in terms of the
covariance matrix in Eq. (174), does not have a counterpart in isolated
systems. To understand the physical significance of $\gamma_{2,k}(t)$, it is
instructive to consider once more the isolated Kitaev chain for comparison.
There, due to the conservation of mode occupation numbers Eq. (166), a memory
of the initial state is kept in Eq. (165), and, therefore, in the GGE. In the
case of the driven-dissipative Kitaev chain, this memory of the initial state
is contained in $\gamma_{1,k}(t)$, and fades away due to the exponential
damping factor $e^{-4\gamma t}$ in Eqs. (157) and (158). At the same time, as
time progresses, the contribution $\gamma_{1,k}(t)$ is being overwritten by
$\gamma_{2,k}(t)$, which, for $t\to\infty$, approaches the stationary form
given in Eq. (174). This picture shows that the PTGGE, as a generalization of
the GGE, is captured by the dephased form of $\gamma_{1,k}(t)$ given in Eq.
(168),

$\gamma_{\mathrm{PTGGE},k}(t)=e^{-4\gamma
t}\gamma_{\mathrm{PTGGE},k}^{\prime},$ (175)

where $\gamma_{\mathrm{PTGGE},k}^{\prime}=\gamma_{d,1,k}^{\prime}$, while the
growth and eventual dominance of $\gamma_{2,k}(t)$ restricts the time frame
during which relaxation to the PTGGE can be observed. As pointed out in the
main text, $\gamma_{2,k}(t)$ vanishes for $\delta=\gamma_{l}-\gamma_{g}=0$.
Then, there is no temporal constraint on the validity of the PTGGE. In
contrast, when $\delta\neq 0$, the PTGGE captures the late-time relaxation of
local observables up to a crossover time scale $t_{\times}$. As stated in the
main text and detailed below, the precise value of $t_{\times}$ depends on the
observable under consideration.

####  III.1.3 From the covariance matrix to the density matrix

According to the above discussion, the PTGGE is uniquely and fully determined
by the dephased form of the contribution to the covariance matrix
$\gamma_{1,k}(t)$ given in Eq. (168). For the sake of completeness, we
describe in the following how the explicit form of the density matrix
corresponding to a given covariance matrix $\gamma_{k}$ can be found.

The dynamics generated by the Liouvillian $\mathcal{L}$ in Eq. (86) couples
only modes with momenta $\pm k$, and we can consider the subspaces
corresponding to these modes separately. In one such subspace, for
$k\notin\\{0,\pi\\}$ such that $k\neq-k$, a complete set of fermionic
operators is given by the doubled spinors of complex $\mathsf{C}_{k}$ and
Majorana $\mathsf{W}_{k}$ fermions:

$\mathsf{C}_{k}=\begin{pmatrix}c_{k}\\\ c_{k}^{\dagger}\\\ c_{-k}\\\
c_{-k}^{\dagger}\end{pmatrix},\qquad\mathsf{W}_{k}=\begin{pmatrix}w_{1,k}\\\
w_{2,k}\\\ w_{3,k}\\\ w_{4,k}\end{pmatrix},$ (176)

which are related by

$\mathsf{C}_{k}=\frac{1}{\sqrt{2}}\mathsf{R}\mathsf{W}_{k},\qquad\mathsf{R}=r\oplus
r.$ (177)

We note that the operators $\mathsf{W}_{k}$ describe true Majorana fermions,
i.e., they are their own antiparticles, in contrast to the Fourier transform
of the Majorana modes defined in Eq. (51). In analogy to Eq. (20), we define a
covariance matrix of $\mathsf{W}_{k}$ fermions as

$\mathsf{y}_{k}=i\left(\langle\mathsf{W}_{k}\mathsf{W}_{k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\rangle-\mathbbm{1}\right)=i\mathsf{R}^{\dagger}\left(2\langle\mathsf{C}_{k}\mathsf{C}_{k}^{\dagger}\rangle-\mathbbm{1}\right)\mathsf{R}.$
(178)

To establish the relation between $\mathsf{y}_{k}$ and the covariance matrix
$\gamma_{k}$ in Eq. (50) as well as the complex covariance matrix $g_{k}$ in
Eq. (47), we note that the spinors in Eq. (29) can be expressed in terms of
the doubled spinors $\mathsf{C}_{k}$ as

$\begin{pmatrix}C_{k}\\\ C_{-k}\end{pmatrix}=\begin{pmatrix}c_{k}\\\
c_{-k}^{\dagger}\\\ c_{-k}\\\
c_{k}^{\dagger}\end{pmatrix}=\mathsf{P}\mathsf{C}_{k},\quad\mathsf{P}=\begin{pmatrix}1&0&0&0\\\
0&0&0&1\\\ 0&0&1&0\\\ 0&1&0&0\end{pmatrix}.$ (179)

Therefore,

$\begin{split}\mathsf{y}_{k}&=i\mathsf{R}^{\dagger}\mathsf{P}\left(2\begin{pmatrix}\langle
C_{k}C_{k}^{\dagger}\rangle&\langle C_{k}C_{-k}^{\dagger}\rangle\\\ \langle
C_{-k}C_{k}^{\dagger}\rangle&\langle
C_{-k}C_{-k}^{\dagger}\rangle\end{pmatrix}-\mathbbm{1}\right)\mathsf{P}\mathsf{R}\\\
&=i\mathsf{R}^{\dagger}\mathsf{P}\begin{pmatrix}g_{k}&0\\\
0&g_{-k}\end{pmatrix}\mathsf{P}\mathsf{R}=i\mathsf{Q}\begin{pmatrix}\gamma_{k}&0\\\
0&\gamma_{-k}\end{pmatrix}\mathsf{Q},\end{split}$ (180)

where

$\mathsf{Q}=\mathsf{R}^{\dagger}\mathsf{P}\mathsf{R}.$ (181)

The corresponding density matrix is given by [34, 35, 36]

$\rho_{k}=\frac{1}{Z_{k}}e^{\frac{i}{2}\mathsf{W}_{k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\mathop{\mathrm{arctanh}}(\mathsf{y}_{k})\mathsf{W}_{k}}=\frac{1}{Z_{k}}e^{-2C_{k}^{\dagger}\mathop{\mathrm{arctanh}}(g_{k})C_{k}+\mathop{\mathrm{tr}}(g_{k})},$
(182)

where the normalization $Z_{k}$ is to be chosen such that
$\mathop{\mathrm{tr}}(\rho_{k})=1$, which leads to

$\rho_{k}=\frac{1}{\det\\!\left(\mathbbm{1}+e^{-2\mathop{\mathrm{arctanh}}(g_{k})}\right)}e^{-2C_{k}^{\dagger}\mathop{\mathrm{arctanh}}(g_{k})C_{k}}.$
(183)

This form applies also to $k\in\\{0,\pi\\}$ when we define, for these momenta,
$C_{k}=c_{k}$ and $g_{k}=\langle[c_{k},c_{k}^{\dagger}]\rangle$. Then, the
full density matrix is given by

$\rho=\prod_{k>0}\rho_{k}.$ (184)

According to the discussion in Sec. III.1.2, we obtain the PTGGE by inserting
$g_{k}(t)=r\gamma_{\mathrm{PTGGE},k}(t)r^{\dagger}$, with
$\gamma_{\mathrm{PTGGE},k}(t)$ given in Eq. (175). In the next section, we
establish the equivalence between the resulting expression for
$\rho_{\mathrm{PTGGE}}(t)$ and the form given in Eq. (8) in the main text in
terms of eigenmodes of the adjoint Liouvillian.

###  III.2 Derivation of the PTGGE from the principle of maximum entropy

As described in the main text, the PTGGE can be obtained as the maximum
entropy ensemble that is compatible with the modified statistics and dephased
expectation values of commutators of eigenmodes of the Liouvillian. In the
following, we present details of this derivation.

####  III.2.1 Eigenmodes of the adjoint Liouvillian

Our goal is to find eigenmodes of the Liouvillian, which generalize the
Bogoliubov quasiparticle operators defined in Eq. (38) to the driven-
dissipative setting. Specifically, for the isolated Kitaev chain, an operator
$O$ evolves in the Heisenberg picture according to $dO/dt=i[H,O]$, i.e., the
dynamics is generated by the superoperator $\mathcal{H}O=[H,O]$. The
Bogoliubov quasiparticle operators $d_{k}$ are eigenmodes of $\mathcal{H}$ in
the sense that they obey the eigenvalue equation
$\mathcal{H}d_{k}=-\varepsilon_{k}d_{k}$. As we explain below, this definition
of eigenmodes of the generator of dynamics cannot directly be generalized to
fermionic open systems. Instead of fermionic mode operators $d_{k}$
themselves, we consider bilinear forms of operators,

$\eta_{k}=[d_{k},d_{k}^{\dagger}],\qquad\chi_{k}=[d_{k},d_{-k}],$ (185)

which in the isolated Kitaev chain obey $\mathcal{H}\eta_{k}=0$ and
$\mathcal{H}\chi_{k}=-2\varepsilon_{k}\chi_{k}$ such that
$\langle\eta_{k}(t)\rangle=\langle\eta_{k}\rangle_{0}$ is conserved, while
$\langle\chi_{k}(t)\rangle$ is purely oscillatory and vanishes via dephasing.
We generalize this structure to open systems in three steps: (i) We specify
the generator of dynamics for operators, and thereby clarify which eigenvalue
equation we actually want to solve. (ii) To find a solution to the eigenvalue
equation, we make a bilinear ansatz in the form of $\eta_{k}$, which leads to
an explicit expression for the mode operators $d_{k}$. We find that the
expectation value $\langle\eta_{k}(t)\rangle$ is not conserved but
nonoscillatory, i.e., not subject to dephasing. (iii) The form of $\eta_{k}$
and, in particular, the mode operators $d_{k}$, completely determines
$\chi_{k}$. As we show, the expectation value $\langle\chi_{k}(t)\rangle$ is
the sum of constant and purely oscillatory contributions. The constant part
vanishes when $\delta=0$, such that, as in the isolated Kitaev chain,
$\langle\chi_{k}(t)\rangle$ does not contribute to the PTGGE.

##### (i) Adjoint Liouvillian.

For any operator $O$, the expectation value $\langle
O\rangle=\mathop{\mathrm{tr}}(O\rho)$ obeys the equation of motion

$\frac{d}{dt}\langle
O\rangle=\mathop{\mathrm{tr}}\\!\left(O\frac{d\rho}{dt}\right)=i\mathop{\mathrm{tr}}\\!\left(O\mathcal{L}\rho\right)=i\mathop{\mathrm{tr}}\\!\left(\rho\mathcal{L}^{\dagger}O\right)=i\langle\mathcal{L}^{\dagger}O\rangle,$
(186)

with the adjoint Liouvillian,

$\mathcal{L}^{\dagger}=\mathcal{H}-i\mathcal{D}^{\dagger},$ (187)

where $\mathcal{H}$ is given in Eq. (88), and

$\begin{split}\mathcal{D}^{\dagger}O&=\sum_{l=1}^{L}\left(2L_{l}^{\dagger}OL_{l}-\\{L_{l}^{\dagger}L_{l},O\\}\right)\\\
&=\sum_{l=1}^{L}\left(L_{l}^{\dagger}[O,L_{l}]+[L_{l}^{\dagger},O]L_{l}\right).\end{split}$
(188)

As a general remark, we note that the adjoint Liouvillian generates the
dynamics of operators only in the sense of expectation values, while the
Heisenberg picture for open systems is defined in terms of quantum Langevin
equations [37, 38].

Now, if $O$ is an eigenmode of $\mathcal{L}^{\dagger}$ in the following sense:

$\mathcal{L}^{\dagger}O=\lambda^{*}\left(O-O_{\mathrm{SS}}\right),$ (189)

where $O_{\mathrm{SS}}\in\mathbb{R}$, then,

$\frac{d}{dt}\langle O\rangle=i\lambda^{*}\left(\langle O\rangle-
O_{\mathrm{SS}}\right),$ (190)

which is solved by

$\langle O(t)\rangle=O_{\mathrm{SS}}+e^{i\lambda^{*}t}\left(\langle
O(0)\rangle-O_{\mathrm{SS}}\right).$ (191)

For $\mathop{\mathrm{Im}}(\lambda^{*})>0$, we obtain $\langle O(t)\rangle\to
O_{\mathrm{SS}}=\langle
O_{\mathrm{SS}}\rangle=\mathop{\mathrm{tr}}(O\rho_{\mathrm{SS}})$ for
$t\to\infty$. As we show in the following, nonoscillatory eigenmodes, for
which $\lambda$ is purely imaginary, $\lambda^{*}=i\kappa$, can be found for
operators $O$ in the form of $\eta_{k}$ in Eq. (185), while an ansatz in the
form of $\chi_{k}$ leads to oscillatory eigenmodes. Here and in the following,
we use the term eigenmodes both for the bilinears $\eta_{k}$ and $\chi_{k}$
and the fermionic operators $d_{k}$.

##### (ii) Nonoscillatory normal commutators.

To find nonoscillatory eigenmodes of the adjoint Liouvillian, we make an
ansatz in the form of $\eta_{k}$ in Eq. (185). In fact, for the time being, we
do not assume translational invariance, i.e., we do not assume that the
eigenmodes can be labelled by the momentum $k$. Then, $\eta=[d,d^{\dagger}]$
can be written in the form

$\eta=\frac{i}{4}\sum_{l,l^{\prime}=1}^{2L}w_{l}Q_{l,l^{\prime}}w_{l^{\prime}},$
(192)

where, without loss of generality, $Q$ can be assumed to be antisymmetric,
$Q=-Q^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}$.
Further, $\eta$ is Hermitian if $Q$ is real, and $\eta$ is traceless by
construction. We proceed to insert our ansatz for $\eta$ in Eq. (189). The
contribution from the Hamiltonian is given by

$\mathcal{H}\eta=-\frac{1}{4}\sum_{l,l^{\prime}=1}^{2L}w_{l}[A,Q]_{l,l^{\prime}}w_{l^{\prime}},$
(193)

and applying the adjoint dissipator to $\eta$ yields

$\begin{split}\mathcal{D}^{\dagger}\eta&=\sum_{m,m^{\prime}=1}^{2L}M_{m^{\prime},m}\left(w_{m}[\eta,w_{m^{\prime}}]+[w_{m},\eta]w_{m^{\prime}}\right)\\\
&=-i\sum_{l,l^{\prime}=1}^{2L}w_{l}\\{M_{R},Q\\}_{l,l^{\prime}}w_{l^{\prime}}+2\mathop{\mathrm{tr}}(M_{I}Q).\end{split}$
(194)

where the bath matrix $M=M_{R}+iM_{I}$ is defined in Eq. (92). We identify the
last, constant term in the above expression with the steady-state expectation
value $\eta_{\mathrm{SS}}$,

$\eta_{\mathrm{SS}}=-\frac{2}{\kappa}\mathop{\mathrm{tr}}(M_{I}Q),$ (195)

where we also anticipate that $\lambda^{*}=i\kappa$ is purely imaginary. Then,
Eq. (189) takes the form

$\sum_{l,l^{\prime}=1}^{2L}w_{l}\left([A,Q]+4\\{M_{R},Q\\}\right)_{l,l^{\prime}}w_{l^{\prime}}=\kappa\sum_{l,l^{\prime}=1}^{2L}w_{l}Q_{l,l^{\prime}}w_{l^{\prime}},$
(196)

Next, we multiply this equation by $w_{l}w_{l^{\prime}}$ and use

$\frac{1}{2^{L+1}}\mathop{\mathrm{tr}}(w_{l^{\prime}}w_{l}w_{m}w_{m^{\prime}})=\frac{1}{2}\left(\delta_{l,l^{\prime}}\delta_{m,m^{\prime}}-\delta_{l,m^{\prime}}\delta_{l^{\prime},m}+\delta_{l,m}\delta_{l^{\prime},m^{\prime}}\right),$
(197)

to obtain

$X^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}Q+QX=\kappa
Q,$ (198)

where matrix $X$ is defined in Eq. (99). Now that we have illustrated the
derivation of the eigenvalue equation for the bilinear $\eta$, let us briefly
comment on the possibility of finding eigenmodes of $\mathcal{L}^{\dagger}$
that are linear in fermionic operators. In general, such linear eigenmodes
cannot be found because there is no closed eigenvalue equation in the space of
linear operators. As a simple example, consider a lattice model without
Hamiltonian dynamics, $H=0$, and incoherent loss,
$L_{l}=\sqrt{\gamma_{l}}c_{l}$. Then, by applying the adjoint dissipator Eq.
(188) to the fermionic operator $c_{l}$ yields

$\mathcal{D}^{\dagger}c_{l}=-\gamma_{l}\left(1+4\sum_{l^{\prime}=1}^{L}c_{l^{\prime}}^{\dagger}c_{l^{\prime}}\right)c_{l},$
(199)

i.e., a sum of linear and cubic terms. In contrast, the second equality in Eq.
(188) shows that, if $O$ is bilinear and $L_{l}$ is linear in fermionic
operators, also $\mathcal{D}^{\dagger}O$ is bilinear, and, therefore, a closed
eigenvalue equation can always be found for a bilinear ansatz. The situation
is very different for bosons: Then, the second equality in Eq. (188) shows
that, if both $O$ and $L_{l}$ are linear in bosonic operators, also
$\mathcal{D}^{\dagger}O$ is linear.

To proceed, we exploit translational invariance of the driven-dissipative
Kitaev chain, which implies that there is an eigenmode $\eta_{k}$, expressed
in terms of a matrix $Q_{k}$, for each momentum mode $k$. We refine our ansatz
in Eq. (192) as follows:

$\eta_{k}=\frac{i}{4}\sum_{l,l^{\prime}=1}^{L}W_{l}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}q_{k,l-l^{\prime}}W_{l^{\prime}},$
(200)

where the spinors $W_{l}$ of Majorana operators are defined in Eq. (9), and
where we assume that $Q_{k}$ is a block Toeplitz matrix with $2\times 2$
blocks

$q_{k,l-l^{\prime}}=\begin{pmatrix}Q_{k,2l-1,2l^{\prime}-1}&Q_{k,2l-1,2l^{\prime}}\\\
Q_{k,2l,2l^{\prime}-1}&Q_{k,2l,2l^{\prime}}\end{pmatrix}.$ (201)

We define the Fourier transform of $q_{k,l}$ as

$q_{k,k^{\prime}}=i\sum_{l=1}^{L}e^{-ik^{\prime}l}q_{k,l}.$ (202)

For future reference, we note that, since $Q$ is real and antisymmetric, the
Fourier transform $q_{k,k^{\prime}}$ of $q_{k,l}$ obeys the relations

$q_{k,k^{\prime}}=-q_{k,-k^{\prime}}^{*},\qquad
q_{k,k^{\prime}}=-q_{k,-k^{\prime}}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}},$
(203)

which imply $q_{k,k^{\prime}}=q_{k,k^{\prime}}^{\dagger}$. By inserting Eqs.
(201) and (202) in Eq. (198) we obtain

$x_{k^{\prime}}^{\dagger}q_{k,k^{\prime}}-q_{k,k^{\prime}}x_{k^{\prime}}=i\kappa
q_{k,k^{\prime}}.$ (204)

In the PT-symmetric phase, according to Eq. (144), the spectrum of $x_{k}$ is
given by $\sigma(x_{k})=\\{-i2\gamma\pm\omega_{k}\\}$, and $x_{k}$ can be
diagonalized as specified in Eq. (153). This leads to

$\omega_{k^{\prime}}\left(S_{k^{\prime}}^{-\dagger}\sigma_{y}S_{k^{\prime}}^{\dagger}q_{k,k^{\prime}}-q_{k,k^{\prime}}S_{k^{\prime}}\sigma_{y}S_{k^{\prime}}^{-1}\right)=i\left(\kappa-4\gamma\right)q_{k,k^{\prime}},$
(205)

where we use the shorthand notation
$S_{k}^{-\dagger}=\left(S_{k}^{\dagger}\right)^{-1}$. On the right-hand side,
we identify $\kappa=4\gamma$, and are thus led to the following commutation
relation:

$\left[\sigma_{y},S_{k^{\prime}}^{\dagger}q_{k,k^{\prime}}S_{k^{\prime}}\right]=0.$
(206)

The matrices that commute with $\sigma_{y}$ are $\mathbbm{1}$ and $\sigma_{y}$
itself. Therefore, the general solution for $q_{k,k^{\prime}}$, which obeys
the above commutation relation, reads

$q_{k,k^{\prime}}=S_{k^{\prime}}^{-\dagger}\left(\alpha_{+,k,k^{\prime}}\sigma_{y}+\alpha_{-,k,k^{\prime}}\mathbbm{1}\right)S_{k^{\prime}}^{-1},$
(207)

with undetermined coefficients $\alpha_{+,k,k^{\prime}}$ and
$\alpha_{-,k,k^{\prime}}$. This solution can be seen to obey the conditions
Eq. (203) if $\alpha_{+,k,k^{\prime}}$ and $\alpha_{-,k,k^{\prime}}$ are,
respectively, symmetric and antisymmetric functions of $k^{\prime}$, and by
using that PHS† of $x_{k}$ in Eq. (124) implies

$S_{k}=S_{-k}^{*},$ (208)

in analogy to Eq. (69). We proceed to fix the coefficients
$\alpha_{+,k,k^{\prime}}$ and $\alpha_{-,k,k^{\prime}}$ such that $\eta_{k}$
can be written in terms of mode operators $d_{k}$ as
$\eta_{k}=[d_{k},d_{k}^{\dagger}]$. To this end, we first express $\eta_{k}$
in terms of complex fermions. By inserting Eq. (51) in Eq. (200) we obtain

$\begin{split}\eta_{k}&=\frac{1}{4}\sum_{k^{\prime}\in\mathrm{BZ}}W_{k^{\prime}}^{\dagger}q_{k,k^{\prime}}W_{k^{\prime}}=\frac{1}{2}\sum_{k^{\prime}\in\mathrm{BZ}}C_{k^{\prime}}^{\dagger}rq_{k,k^{\prime}}r^{\dagger}C_{k^{\prime}}\\\
&=\frac{1}{2}\sum_{k^{\prime}\in\mathrm{BZ}}C_{k^{\prime}}^{\dagger}U_{k^{\prime}}^{-\dagger}\left(\alpha_{+,k,k^{\prime}}\sigma_{z}+\alpha_{-,k,k^{\prime}}\mathbbm{1}\right)U_{k^{\prime}}^{-1}C_{k^{\prime}},\end{split}$
(209)

where, in the last line, we use Eq. (154). For solutions $\eta_{k}$ with fixed
momentum $k$, the sum over $k^{\prime}$ should collapse to the two terms
$k^{\prime}=\pm k$. Symmetric and antisymmetric combinations of these
components are given by

$\alpha_{\pm,k,k^{\prime}}=-\frac{\alpha_{0,k}}{2}\left(\delta_{k^{\prime},k}\pm\delta_{k^{\prime},-k}\right),$
(210)

with the coefficient $\alpha_{0,k}$ to be determined below. Then, with

$P_{z,\pm}=\frac{1}{2}\left(\mathbbm{1}\pm\sigma_{z}\right),$ (211)

we find

$\begin{split}\eta_{k}&=\frac{\alpha_{0,k}}{2}\left(C_{-k}^{\dagger}U_{-k}^{-\dagger}P_{z,-}U_{-k}^{-1}C_{-k}-C_{k}^{\dagger}U_{k}^{-\dagger}P_{z,+}U_{k}^{-1}C_{k}\right)\\\
&=\frac{\alpha_{0,k}}{2}\left(C_{k}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}U_{k}^{-{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}P_{z,+}U_{k}^{-*}C_{k}^{\dagger{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}-C_{k}^{\dagger}U_{k}^{-\dagger}P_{z,+}U_{k}^{-1}C_{k}\right),\end{split}$
(212)

where in the second equality we use that $U_{k}$, which is defined in Eq.
(150), obeys the same PHS as $U_{x,\theta_{k}}$ in Eq. (69). $\eta_{k}$ takes
its final form, $\eta_{k}=[d_{k},d_{k}^{\dagger}]$, if we define the
eigenmodes of the Liouvillian as

$D_{k}=\begin{pmatrix}d_{k}\\\
d_{-k}^{\dagger}\end{pmatrix}=V_{k}^{\dagger}C_{k},$ (213)

where

$V_{k}=\sqrt{\frac{2}{\mathop{\mathrm{tr}}\\!\left(U_{k}^{-1}U_{k}^{-\dagger}\right)}}U_{k}^{-\dagger}=\begin{pmatrix}\cos\\!\left(\frac{\theta_{k}+\phi_{k}}{2}\right)&i\sin\\!\left(\frac{\theta_{k}-\phi_{k}}{2}\right)\\\
i\sin\\!\left(\frac{\theta_{k}+\phi_{k}}{2}\right)&\cos\\!\left(\frac{\theta_{k}-\phi_{k}}{2}\right)\end{pmatrix}.$
(214)

The angle $\theta_{k}$ is defined in Eq. (35), and $\phi_{k}$ is defined
through the equality

$\varepsilon_{k}e^{i\phi_{k}}=\omega_{k}+i2\sqrt{\gamma_{l}\gamma_{g}},$ (215)

and is related to $\beta_{k}$ defined in Eq. (151) via

$\tan(\phi_{k})=-\frac{\varepsilon_{k}}{\omega_{k}}\tanh(\beta_{k}).$ (216)

When we insert Eq. (213) in Eq. (212), the normalization in Eq. (214) cancels
if we set

$\alpha_{0,k}=\mathop{\mathrm{tr}}\\!\left(U_{k}^{-1}U_{k}^{-\dagger}\right)=2\cosh(\beta_{k}).$
(217)

The normalization in Eq. (214) is chosen such that
$\\{d_{k},d_{k}^{\dagger}\\}=1$. As in the main text, we collect the
anticommutators of the modes $d_{k}$ in a matrix,

$\begin{pmatrix}\\{d_{k},d_{k^{\prime}}^{\dagger}\\}&\\{d_{k},d_{-k^{\prime}}\\}\\\
\\{d_{-k}^{\dagger},d_{k^{\prime}}^{\dagger}\\}&\\{d_{-k}^{\dagger},d_{-k^{\prime}}\\}\end{pmatrix}=f_{k}\delta_{k,k^{\prime}},$
(218)

where

$f_{k}=\begin{pmatrix}\\{d_{k},d_{k}^{\dagger}\\}&\\{d_{k},d_{-k}\\}\\\
\\{d_{-k}^{\dagger},d_{k}^{\dagger}\\}&\\{d_{-k}^{\dagger},d_{-k}\\}\end{pmatrix}=V_{k}^{\dagger}V_{k}=\mathbbm{1}+\sin(\phi_{k})\sigma_{y}.$
(219)

Next, we evaluate Eq. (195). A short calculation yields

$\eta_{\mathrm{SS},k}=\langle[d_{k},d_{k}^{\dagger}]\rangle_{\mathrm{SS}}=-\frac{\delta\omega_{k}}{4\varepsilon_{k}\gamma}\left(\cos(\theta_{k})-\frac{2\sqrt{\gamma_{l}\gamma_{g}}}{\omega_{k}}\sin(\theta_{k})\right).$
(220)

Finally, to completely specify the time evolution of the expectation value of
$\eta_{k}$, we need to know the expectation value $\langle\eta_{k}\rangle_{0}$
in the initial state. Starting from Eq. (212), we use Eqs. (47) to find

$\langle[d_{k},d_{k}^{\dagger}]\rangle=\frac{\omega_{k}}{\varepsilon_{k}}\mathop{\mathrm{tr}}\\!\left(U_{k}^{-\dagger}P_{z,+}U_{k}^{-1}g_{k}\right).$
(221)

For the ground state, with covariance matrix given in Eq. (49), we obtain

$\langle[d_{k},d_{k}^{\dagger}]\rangle_{0}=\cos(\Delta\theta_{k}-\phi_{k}),$
(222)

where $\Delta\theta_{k}$ is defined in Eq. (171). Putting everything together,
we obtain the result stated in the main text:

$\langle[d_{k},d_{k}^{\dagger}](t)\rangle=e^{-4\gamma
t}\langle[d_{k},d_{k}^{\dagger}]\rangle_{0}+\left(1-e^{-4\gamma
t}\right)\langle[d_{k},d_{k}^{\dagger}]\rangle_{\mathrm{SS}},$ (223)

where $\langle[d_{k},d_{k}^{\dagger}]\rangle_{0}$ and
$\langle[d_{k},d_{k}^{\dagger}]\rangle_{\mathrm{SS}}$ are given in Eqs. (222)
and (220), respectively.

##### (iii) Oscillatory anomalous commutators.

We next show that also $\chi_{k}=[d_{k},d_{-k}]$ is an eigenmode of the
Liouvillian in the sense of Eq. (189), however, where
$\mathop{\mathrm{Re}}(\lambda^{*})\neq 0$, such that the nonconstant term in
Eq. (191) is oscillatory. To that end, in analogy to Eq. (209), we write
$\chi_{k}$ as

$\chi_{k}=\frac{1}{2}\sum_{k^{\prime}\in\mathrm{BZ}}C_{k^{\prime}}^{\dagger}rp_{k,k^{\prime}}r^{\dagger}C_{k^{\prime}},$
(224)

where we define

$p_{k,k^{\prime}}=-2\left(\delta_{k,k^{\prime}}-\delta_{k,-k^{\prime}}\right)r^{\dagger}V_{k^{\prime}}\sigma_{-}V_{k^{\prime}}^{\dagger}r.$
(225)

A straightforward calculation shows that $p_{k,k^{\prime}}$ satisfies an
eigenvalue equation similar to Eq. (204),

$x_{k^{\prime}}^{\dagger}p_{k,k^{\prime}}-p_{k,k^{\prime}}x_{k^{\prime}}=-2\left(\omega_{k}-i2\gamma\right)p_{k,k^{\prime}},$
(226)

which leads us to identify $\lambda^{*}=-2\left(\omega_{k}-i2\gamma\right)$.
Consequently, as stated in the main text, we find

$\langle[d_{k},d_{-k}](t)\rangle=e^{-i2\left(\omega_{k}-i2\gamma\right)t}\langle[d_{k},d_{-k}]\rangle_{0}\\\
+\left(1-e^{-i2\left(\omega_{k}-i2\gamma\right)t}\right)\langle[d_{k},d_{-k}]\rangle_{\mathrm{SS}},$
(227)

where the steady-state expectation value is given by

$\begin{split}\chi_{\mathrm{SS},k}&=\langle[d_{k},d_{-k}]\rangle_{\mathrm{SS}}\\\
&=\frac{i\delta}{4\left(\omega_{k}-i2\gamma\right)}\sum_{k^{\prime}\in\mathrm{BZ}}\mathop{\mathrm{tr}}(\sigma_{y}p_{k,k^{\prime}})=\frac{\delta\sin(\theta_{k})}{\omega_{k}-i2\gamma}.\end{split}$
(228)

####  III.2.2 Maximum entropy ensemble

The PTGGE describes the late-time behavior for $\delta=0$, such that, as can
be seen explicitly in Eqs. (220) and (228),

$\langle[d_{k},d_{k}^{\dagger}]\rangle_{\mathrm{SS}}=\langle[d_{k},d_{-k}]\rangle_{\mathrm{SS}}=0,$
(229)

and when the contribution from the anomalous commutators becomes negligible
due to dephasing,

$\langle[d_{k},d_{-k}](t)\rangle\to 0.$ (230)

To summarize, the PTGGE is characterized by the following expectation values:

$\begin{pmatrix}\langle[d_{k},d_{k}^{\dagger}](t)\rangle_{\mathrm{PTGGE}}&\langle[d_{k},d_{-k}](t)\rangle_{\mathrm{PTGGE}}\\\
\langle[d_{-k}^{\dagger},d_{k}^{\dagger}](t)\rangle_{\mathrm{PTGGE}}&\langle[d_{-k}^{\dagger},d_{-k}](t)\rangle_{\mathrm{PTGGE}}\end{pmatrix}=\zeta_{k}(t),$
(231)

where

$\zeta_{k}(t)=e^{-4\gamma t}\zeta_{k}^{\prime},$ (232)

with

$\begin{split}\zeta_{k}^{\prime}&=\begin{pmatrix}\langle[d_{k},d_{k}^{\dagger}]\rangle_{0}&0\\\
0&\langle[d_{-k}^{\dagger},d_{-k}]\rangle_{0}\end{pmatrix}\\\
&=\begin{pmatrix}\cos(\Delta\theta_{k}-\phi_{k})&0\\\
0&-\cos(\Delta\theta_{k}+\phi_{k})\end{pmatrix}.\end{split}$ (233)

Specifically, we define the PTGGE as the maximum entropy ensemble that is
compatible with the expectation values of commutators summarized in
$\zeta_{k}(t)$, and with the statistics of the modes $d_{k}$ as determined by
their anticommutators in Eq. (219). The same information is stored in the
covariance matrix in the complex basis, which, according to Eqs. (47) and
(213), is given by

$g_{\mathrm{PTGGE},k}(t)=V_{k}^{-\dagger}\zeta_{k}(t)V_{k}^{-1}=e^{-4\gamma
t}g_{\mathrm{PTGGE},k}^{\prime},$ (234)

where

$g_{\mathrm{PTGGE},k}^{\prime}=V_{k}^{-\dagger}\zeta_{k}^{\prime}V_{k}^{-1}.$
(235)

For given two-point functions, the entropy is maximized for a Gaussian state,
which is determined uniquely by Eq. (183). Therefore, with

$C_{k}^{\dagger}\mathop{\mathrm{arctanh}}\\!\left(g_{\mathrm{PTGGE},k}(t)\right)C_{k}=D_{k}f_{k}^{-1}\mathop{\mathrm{arctanh}}\\!\left(\zeta_{k}(t)f_{k}^{-1}\right)D_{k},$
(236)

we obtain the PTGGE as stated in the main text:

$\rho_{\mathrm{PTGGE}}(t)=\frac{1}{Z_{\mathrm{PTGGE}}(t)}e^{-2\sum_{k\geq
0}D_{k}^{\dagger}f_{k}^{-1}\mathop{\mathrm{arctanh}}\left(\zeta_{k}(t)f_{k}^{-1}\right)D_{k}}.$
(237)

Finally, a straightforward calculation shows that
$\gamma_{\mathrm{PTGGE},k}(t)=r^{\dagger}g_{\mathrm{PTGGE},k}(t)r$ for
$\gamma_{\mathrm{PTGGE},k}(t)$ given in Eq. (175) and that, therefore, the
derivations of the PTGGE in Secs. III.1 and III.2 yield the same result.
However, the key benefit of the derivation presented in the current section is
that it clarifies how the notions of eigenmodes of the generator of the
dynamics and of a maximum entropy ensemble generalize to PT-symmetric open
systems. In particular, this derivation shows that conservation of occupation
numbers of eigenmodes is not necessary for the PTGGE to capture the late-time
limit of quench dynamics. Instead, nonoscillatory decay is sufficient.

###  III.3 Subsystem parity

In the following, we first describe how to obtain the data for the subsystem
fermion parity in Eq. (84) that is shown in Figs. 1 and 2 in the main text,
and the PTGGE prediction that describes the late-time behavior of the
subsystem parity and is shown as horizontal lines in the figures. Then, we
fill in details of the derivation of the conjectured analytical expressions
for the space-time scaling limit of the subsystem parity given in Eqs. (9) and
(10) in the main text. Next, we derive the estimate for the crossover time
$t_{\times}$ that is shown in the inset of Fig. 2, and finally, we derive
explicit expressions for the soft modes $k_{s,\pm}$, and we argue that the
existence of two distinct time scales $t_{s,+}$ and $t_{s,-}$ is a unique
feature of driven-dissipative systems with broken inversion symmetry.

####  III.3.1 Numerical time evolution

In the PT-symmetric phase, it is convenient to define a rescaled covariance
matrix as

$\Gamma^{\prime}(t)=e^{4\gamma t}\Gamma(t),$ (238)

where $\Gamma(t)$ is given in Eqs. (101) and (102). Through this rescaling,
one can avoid having to perform numerical calculations with extremely small
quantities when $\gamma t\gg 1$. In particular, we obtain the rescaled
subsystem parity, which is shown in Figs. 1 and 2 in the main text, directly
as

$e^{4\ell\gamma t}\langle
P_{\ell}(t)\rangle=\mathop{\mathrm{Pf}}(\Gamma^{\prime}(t)),$ (239)

where we use the following property of the Pfaffian: for a $2\ell\times 2\ell$
matrix $A$ and $\alpha\in\mathbb{C}$, $\mathop{\mathrm{Pf}}(\alpha
A)=\alpha^{\ell}\mathop{\mathrm{Pf}}(A)$.

####  III.3.2 Long-time behavior

As we show in Sec. III.2, the PTGGE is equivalent to the dephased covariance
matrix $\gamma_{\mathrm{PTGGE},k}(t)$ in Eq. (175). The corresponding
covariance matrix in position space can be found through an inverse Fourier
transform according to Eq. (50) and by using Eq. (21). In analogy to Eq.
(175), we obtain the late-time behavior in the form

$\Gamma(t)\sim\Gamma_{\mathrm{PTGGE}}(t)=e^{-4\gamma
t}\Gamma_{\mathrm{PTGGE}}^{\prime},$ (240)

where $\Gamma_{\mathrm{PTGGE}}^{\prime}$ is time-independent. The late-time
behavior of the subsystem parity, which according to Eq. (84) is just the
Pfaffian of the reduced covariance matrix Eq. (85), is thus given by

$\langle P_{\ell}(t)\rangle\sim\langle
P_{\ell}(t)\rangle_{\mathrm{PTGGE}}=e^{-4\ell\gamma
t}\mathop{\mathrm{Pf}}\\!\left(\Gamma_{\mathrm{PTGGE},\ell}^{\prime}\right).$
(241)

We obtain the horizontal lines in Figs. 1 and 2 of the main text by evaluating
the right-hand side of the above equation numerically. In particular, we
reconstruct $\Gamma_{\mathrm{PTGGE},\ell}^{\prime}$ from
$\gamma_{\mathrm{PTGGE},k}^{\prime}$ defined in Eq. (175) with the aid of Eqs.
(21) and (50). In the latter, we impose the thermodynamic limit $L\to\infty$
through

$\gamma_{l}=\frac{i}{L}\sum_{k\in\mathrm{BZ}}e^{ikl}\gamma_{k}\sim
i\int_{-\pi}^{\pi}\frac{dk}{2\pi}\,e^{ikl}\gamma_{k}.$ (242)

####  III.3.3 Space-time scaling limit

As pointed out in Sec. I.5, for the isolated Kitaev chain, the subsystem
fermion parity can be mapped to order parameter correlations in the transverse
field Ising model through the combination of a Jordan-Wigner transformation
and the Kramers-Wannier duality. The relaxation of order parameter
correlations after a quantum quench is studied in Refs. [9, 10, 11]. In
particular, for a quench to the topologically trivial phase of the isolated
Kitaev chain with $\left\lvert\mu\right\rvert>2J$, which corresponds to a
quench to the ferromagnetic phase of the Ising model, Eq. (19) of Ref. [10]
yields the following prediction for the behavior of the subsystem parity in
the space-time scaling limit $\ell,t\to\infty$ with $\ell/t$ fixed:

$\langle P_{\ell}(t)\rangle\sim
P_{0}e^{\int_{0}^{\pi}\frac{dk}{\pi}\min(2\left\lvert v_{k}\right\rvert
t,\ell)\ln(\left\lvert\cos(\Delta\theta_{k})\right\rvert)},$ (243)

where $v_{k}=d\varepsilon_{k}/dk$. The first step to obtain our conjecture in
Eq. (9) in the main text for the subsystem parity in the driven-dissipative
Kitaev chain is to rewrite the above equation in terms of the covariance
matrix in the complex basis:

$\begin{split}\ln(\left\lvert\cos(\Delta\theta_{k})\right\rvert)&=\frac{1}{2}\left(\ln(\left\lvert\cos(\Delta\theta_{k})\right\rvert)+\ln(\left\lvert\cos(\Delta\theta_{-k})\right\rvert)\right)\\\
&=\frac{1}{2}\left(\ln\\!\left(\left\lvert\langle[d_{k}^{\dagger},d_{k}]\rangle_{0}\right\rvert\right)+\ln\\!\left(\left\lvert\langle[d_{-k},d_{-k}^{\dagger}]\rangle_{0}\right\rvert\right)\right)\\\
&=\frac{1}{2}\mathop{\mathrm{tr}}\\!\left(\ln\\!\left(\left\lvert
g_{\mathrm{GGE},k}\right\rvert\right)\right),\end{split}$ (244)

where by taking the limit $\gamma\to 0$ in Eqs. (231), (232), (233), and
(234), we find that $g_{\mathrm{GGE},k}$ is given by

$g_{\mathrm{GGE},k}=U_{x,\theta_{k}}\begin{pmatrix}\langle[d_{k},d_{k}^{\dagger}]\rangle_{0}&0\\\
0&\langle[d_{-k}^{\dagger},d_{-k}]\rangle_{0}\end{pmatrix}U_{x,\theta_{k}}^{\dagger}.$
(245)

Then, Eq. (243) can be written as

$\langle P_{\ell}(t)\rangle\sim
P_{0}e^{\int_{0}^{\pi}\frac{dk}{2\pi}\min(2\left\lvert v_{k}\right\rvert
t,\ell)\mathop{\mathrm{tr}}\left(\ln\left(\left\lvert
g_{\mathrm{GGE},k}\right\rvert\right)\right)}.$ (246)

To generalize this relation to the driven-dissipative Kitaev chain, we note
that as in Eq. (241), the exponential decay of $g_{\mathrm{PTGGE},k}(t)$ in
Eq. (234) leads to the appearance of a prefactor $e^{-4\ell\gamma t}$.
Further, we set $v_{k}=d\omega_{k}/dk$, and we replace $g_{\mathrm{GGE},k}$ by
$g_{\mathrm{PTGGE},k}^{\prime}$. Finally, by using the cyclic invariance of
the trace and Eq. (219), we obtain

$\mathop{\mathrm{tr}}\\!\left(\ln\\!\left(\left\lvert
g_{\mathrm{PTGGE},k}^{\prime}\right\rvert\right)\right)=\mathop{\mathrm{tr}}\\!\left(\ln\\!\left(\left\lvert\zeta_{k}^{\prime}f_{k}^{-1}\right\rvert\right)\right),$
(247)

and are thus lead to Eq. (9) of the main text. As stated there, the constant
prefactor $P_{0}$ can be determined by fitting the late-time limit,

$\lim_{t\to\infty}e^{4\ell\gamma t}\langle P_{\ell}\rangle\sim
P_{0}e^{\ell\int_{0}^{\pi}\frac{dk}{2\pi}\mathop{\mathrm{tr}}\left(\ln\left(\left\lvert\zeta_{k}^{\prime}f_{k}^{-1}\right\rvert\right)\right)},$
(248)

to the numerical result in Eq. (241).

Next, we consider quenches to the topological phase of the isolated Kitaev
chain for $\left\lvert\mu\right\rvert<2J$, corresponding to quenches to the
paramagnetic phase of the transverse field Ising model. According to Eq. (32)
of Ref. [10], for $t<t_{F}$, the oscillatory decay of the subsystem parity in
the space-time scaling limit is given by

$\langle P_{\ell}(t)\rangle\sim
P_{0}\left(1+\cos(2\varepsilon_{k_{s}}t+\alpha)\right)e^{2t\int_{0}^{\pi}\frac{dk}{\pi}\left\lvert
v_{k}\right\rvert\ln(\left\lvert\cos(\Delta\theta_{k})\right\rvert)},$ (249)

where in the isolated system $k_{s}=k_{s,+}=k_{s,-}$. To generalize this
prediction to the driven-dissipative Kitaev chain, we proceed as above, and,
additionally, we note that the oscillatory prefactor is a complete square due
to the identity $2\cos(x)^{2}=1+\cos(2x)$. Based on this simple observation,
we conjecture that the proper generalization of the prefactor reads as
follows:

$2\cos(\varepsilon_{k_{s}}t+\alpha/2)^{2}\to
2\cos(\omega_{k_{s,+}}t+\alpha_{+})\cos(\omega_{k_{s,-}}t+\alpha_{-}),$ (250)

where, for $\gamma\to 0$, we set $\alpha_{+}=\alpha_{-}=\alpha/2$. Thereby we
obtain Eq. (10) of the main text.

####  III.3.4 Condition to observe relaxation to the PTGGE

The condition for the existence of a finite time window during which
relaxation to the PTGGE can be observed depends on the observable under
consideration. Here, we focus on the subsystem parity given in Eq. (84), and
on quenches that originate from the trivial phase of the isolated Kitaev
chain. For $\delta=0$, when $\Gamma_{2}(t)=0$ in Eq. (100), the behavior of
the subsystem parity is given by $\langle
P_{\ell}(t)\rangle=\mathop{\mathrm{Pf}}(\Gamma_{1,\ell}(t))$, and, as shown in
Figs. 1 and 2 of the main text, $\langle P_{\ell}(t)\rangle$ relaxes on a time
scale $t_{F}$ to the PTGGE prediction given in Eq. (241) 111More precisely, at
$t=t_{F}$, the initial fast relaxation of $\langle P_{\ell}(t)\rangle$ turns
into much slow decay towards the PTGGE prediction, which is reached after a
time $t\sim\ell^{4/3}>t_{F}\sim\ell$. We note that for quenches from the
trivial phase, we expect the subsystem parity to relax to the PTGGE prediction
on an even longer time scale that is exponentially large in $\ell$ [11]. If
$\delta\neq 0$, for relaxation of $\langle P_{\ell}(t)\rangle$ to the PTGGE
prediction to be observable at all, the contribution to $\langle
P_{\ell}(t)\rangle$ due to $\Gamma_{2}(t)\neq 0$ has to be small at $t=t_{F}$.
Then, for $t>t_{F}$, $\langle P_{\ell}(t)\rangle$ obeys the PTGGE prediction
up to the crossover time scale $t_{\times}$. To estimate $t_{\times}$, we
expand the Pfaffian to lowest nontrivial order,

$\langle
P_{\ell}\rangle=\mathop{\mathrm{Pf}}\\!\left(\Gamma_{1,\ell}+\Gamma_{2,\ell}\right)\sim\mathop{\mathrm{Pf}}\\!\left(\Gamma_{1,\ell}\right)\left(1+\frac{1}{2}\mathop{\mathrm{tr}}\\!\left(\frac{\Gamma_{2,\ell}}{\Gamma_{1,\ell}}\right)\right),$
(251)

where we omit the explicit time dependence. Therefore, the condition for
relaxation to the PTGGE to be observable can be stated as

$\frac{1}{2}\left\lvert\mathop{\mathrm{tr}}\\!\left(\frac{\Gamma_{2,\ell}}{\Gamma_{1,\ell}}\right)\right\rvert\ll
1\quad\text{for }t=t_{F}.$ (252)

At $t=t_{F}$, we may estimate $\Gamma_{1,\ell}$ and $\Gamma_{2,\ell}$ by their
“dephased” expressions, i.e., by keeping only nonoscillatory contributions.
According to Eq. (167), and with Eqs. (50) and (21), we can write

$\Gamma_{1}(t)\sim\Gamma_{d,1}(t),\qquad\Gamma_{2}(t)\sim\Gamma_{d,2}(t)+\Gamma_{\mathrm{SS}},$
(253)

where

$\Gamma_{d,1}(t)\sim e^{-4\gamma
t}\Gamma_{d,1}^{\prime},\qquad\Gamma_{d,2}(t)\sim e^{-4\gamma
t}\Gamma_{d,2}^{\prime}.$ (254)

Then, Eq. (252) becomes

$\frac{1}{2}\left\lvert\mathop{\mathrm{tr}}\\!\left(\frac{\Gamma_{d,2,\ell}^{\prime}}{\Gamma_{d,1,\ell}^{\prime}}\right)+e^{4\gamma
t_{F}}\mathop{\mathrm{tr}}\\!\left(\frac{\Gamma_{\mathrm{SS},\ell}}{\Gamma_{d,1,\ell}^{\prime}}\right)\right\rvert\ll
1.$ (255)

A sufficient condition for the validity of this inequality can be obtained by
using the triangle inequality. We find

$\frac{1}{2}\left(\left\lvert\mathop{\mathrm{tr}}\\!\left(\frac{\Gamma_{d,2,\ell}^{\prime}}{\Gamma_{d,1,\ell}^{\prime}}\right)\right\rvert+e^{4\gamma
t_{F}}\left\lvert\mathop{\mathrm{tr}}\\!\left(\frac{\Gamma_{\mathrm{SS},\ell}}{\Gamma_{d,1,\ell}^{\prime}}\right)\right\rvert\right)\ll
1,$ (256)

which is satisfied if $t_{F}$ is small in comparison to the crossover time
scale $t_{\times}$, for which the left-hand side of the above inequality is
equal to one,

$t_{F}\ll
t_{\times}=\frac{1}{4\gamma}\ln\\!\left(\frac{2-\left\lvert\mathop{\mathrm{tr}}(\Gamma_{d,2,\ell}^{\prime}/\Gamma_{d,1,\ell}^{\prime})\right\rvert}{\left\lvert\mathop{\mathrm{tr}}(\Gamma_{\mathrm{SS},\ell}/\Gamma_{d,1,\ell}^{\prime})\right\rvert}\right).$
(257)

For $\delta\to 0$, by taking into account that $\Gamma_{\mathrm{SS},\ell}$ and
$\Gamma_{d,2,\ell}^{\prime}$ vanish linearly in $\delta$ according to Eqs.
(173), (174), and (110), while $\Gamma_{d,1,\ell}^{\prime}$ is finite, we find

$t_{\times}\sim\frac{1}{4\gamma}\left\lvert\ln(c_{\times}\left\lvert\delta\right\rvert)\right\rvert,\quad
c_{\times}=\lim_{\delta\to
0}\frac{1}{2\delta}\left\lvert\mathop{\mathrm{tr}}\\!\left(\frac{\Gamma_{\mathrm{SS},\ell}}{\Gamma_{d,1,\ell}^{\prime}}\right)\right\rvert.$
(258)

The limit $\delta\to 0$ can be taken explicitly in Eq. (174),

$\lim_{\delta\to
0}\frac{\gamma_{\mathrm{SS},k}}{\delta}=-\frac{1}{2\gamma}\left(\frac{1}{\varepsilon_{k}^{2}}\hat{\mathbf{e}}_{y}\cdot\left(\mathbf{a}_{k}\times\mathbf{b}\right)\mathbbm{1}\right.\\\
\left.-\left\\{\hat{\mathbf{e}}_{y}-\left[\varepsilon_{k}\mathbf{e}_{y,\perp,k}-\frac{2\gamma}{\varepsilon_{k}}\mathbf{e}_{y,o,k}\right]\right\\}\cdot\boldsymbol{\sigma}\right),$
(259)

where we set $\mathbf{b}\to-2\gamma\hat{\mathbf{e}}_{z}$ for $\delta\to 0$,
and where

$\begin{split}\mathbf{e}_{y,\parallel,k}&=\left(\hat{\mathbf{e}}_{y}\cdot\hat{\mathbf{a}}_{k}\right)\hat{\mathbf{a}}_{k},\\\
\mathbf{e}_{y,\perp,k}&=\hat{\mathbf{e}}_{y}-\mathbf{e}_{y,\parallel,k},\\\
\mathbf{e}_{o,k}&=-\hat{\mathbf{e}}_{y}\times\hat{\mathbf{a}}_{k}.\end{split}$
(260)

This allows us to reconstruct $\lim_{\delta\to
0}\Gamma_{\mathrm{SS},\ell}/\delta$ numerically by employing Eqs. (50) and
(21). Similarly, we can obtain $\Gamma_{d,1,\ell}^{\prime}$ from
$\gamma_{d,1,k}^{\prime}$, and, thereby, calculate $c_{\times}$. This leads to
the semianalytical estimate in the inset of Fig. 2.

####  III.3.5 Oscillatory dynamics from zero modes

According to Eq. (233), soft modes of the PTGGE are determined by the
condition

$\cos(\Delta\theta_{k}-\phi_{k})=0.$ (261)

For $\phi_{k}=0$, due to $\Delta\theta_{k}=-\Delta\theta_{-k}$, which
according to Eq. (70) is a consequence of IS, the existence of a soft mode
$k_{s,+}>0$ implies the existence of another inversion conjugate soft mode
$k_{s,-}=-k_{s,+}$. Due to $\varepsilon_{k}=\varepsilon_{-k}$, the associated
oscillation periods $t_{s,\pm}=\pi/(2\varepsilon_{k_{s},\pm})$ are identical.
Interestingly, IS† as stated in Eq. (132) implies, via Eqs. (155) and (216),
that $\phi_{k}=\phi_{-k}$. However, the sum $\Delta\theta_{k}-\phi_{k}$ that
appears in Eq. (261) does not have definite transformation behavior under
inversion, just as the sums and differences of $\theta_{k}$ and $\phi_{k}$
that appear in the Liouvillian eigenmode transformation matrix $V_{k}$ in Eq.
(214). Consequently, the solutions $k_{s,\pm}$ of Eq. (261) are not related by
inversion. In particular, for quenches with pre- and postquench parameters
given by, respectively, $J_{0}=\Delta_{0}$ and $\mu_{0}$, and $J=\Delta$ and
$\mu$, we find

$k_{s,\pm}=\mp\mathop{\mathrm{sgn}}(\mu)\arccos\\!\left(-\frac{1}{2\left(J_{0}\mu+J\mu_{0}\right)^{2}}\right.\\\
\times\left\\{4J^{2}J_{0}\mu_{0}+J_{0}\left(\mu^{2}-4\gamma_{l}\gamma_{g}\right)\mu_{0}+J\mu\left(4J_{0}^{2}+\mu_{0}^{2}\right)\mp
2\mathop{\mathrm{sgn}}(\mu)\vphantom{\sqrt{\gamma_{l}\gamma_{g}\left[4J_{0}^{4}\mu^{2}-J_{0}^{2}\left(4J^{2}+\mu^{2}-4\gamma_{l}\gamma_{g}\right)\mu_{0}^{2}+J^{2}\mu_{0}^{4}\right]}}\right.\\\
\left.\left.\times\sqrt{\gamma_{l}\gamma_{g}\left[4J_{0}^{4}\mu^{2}-J_{0}^{2}\left(4J^{2}+\mu^{2}-4\gamma_{l}\gamma_{g}\right)\mu_{0}^{2}+J^{2}\mu_{0}^{4}\right]}\right\\}\right).$
(262)

In the limit $\mu_{0}\to-\infty$ as in the examples considered in the main
text, the above general form reduces to

$k_{s,\pm}=\mp\mathop{\mathrm{sgn}}(\mu)\arccos\\!\left(-\frac{1}{2J}\left(\mu\mp\mathop{\mathrm{sgn}}(\mu)2\sqrt{\gamma_{l}\gamma_{g}}\right)\right).$
(263)

####  III.3.6 Necessary conditions for $t_{s,+}\neq t_{s,-}$

As claimed in the main text, different periods of parity pumping for
subsystems at the left and right ends of the chain require both mixedness of
the time-evolved state and breaking of inversion symmetry and are, therefore,
unique to driven-dissipative systems. To see that this is the case, consider,
for simplicity, an equal bipartition into subsystems
$L_{L/2}=\\{1,\dotsc,L/2\\}$ and $R_{L/2}=\\{L/2+1,\dotsc,L\\}$, corresponding
to the left and right halves of the chain. Then, for a pure state, the
eigenvalues of $\Gamma_{L}$ and $\Gamma_{R}$, which determine the entanglement
spectra for subsystems $L_{L/2}$ and $R_{L/2}$, are identical 222For a
bipartition into subsystems of different sizes, e.g., $L$ smaller than $R$,
the spectrum of $\Gamma_{R}$ is given by the spectrum of $\Gamma_{L}$ and
additional pairs of eigenvalues $\pm i$.. Consequently, due to
$\det(\Gamma_{\ell})=\mathop{\mathrm{Pf}}(\Gamma_{\ell})^{2}=\langle
P_{\ell}\rangle^{2}$, zero crossings of $\langle P_{L}\rangle$ and $\langle
P_{R}\rangle$ occur at the same time. In contrast, for a mixed state, the
spectra of $\Gamma_{L}$ and $\Gamma_{R}$ are, in general, different. However,
if the state is symmetric under inversion, with inversion center in the middle
of the chain, the matrices $\Gamma_{L}$ and $\Gamma_{R}$ are, as we show
below, unitarily equivalent. Consequently, even for mixed states, the
entanglement spectra and the subsystem parities, given by the Pfaffians, are
the same. Therefore, as stated above, both mixedness of the state and
inversion symmetry breaking are necessary for $t_{s,+}\neq t_{s,-}$.

To show that $\Gamma_{L}$ and $\Gamma_{R}$ are unitarily equivalent for an
inversion-symmetric state, we first note that Eq. (59) for open boundary
conditions implies the following transformation behavior of the spinors of
Majorana operators defined in Eq. (9):

$IW_{l}I^{\dagger}=i\sigma_{y}W_{L+1-l}.$ (264)

Therefore, with Eq. (21), we find that the $2\times 2$ blocks of the
covariance matrix of an inversion-symmetric state $\rho=I\rho I^{\dagger}$
obey the condition

$\gamma_{l,l^{\prime}}=\sigma_{y}\gamma_{L+1-l,L+1-l^{\prime}}\sigma_{y},$
(265)

which implies that the full covariance matrix behaves as

$\Gamma=\Sigma_{y,L}\Xi_{L}\Gamma\Xi_{L}\Sigma_{y,L},$ (266)

where $\Sigma_{y,\ell}$ and $\Xi_{\ell}$ are defined in Eqs. (142) and (143),
respectively. The product $\Sigma_{y,\ell}\Xi_{\ell}$ can be seen to be
unitary by noting that both $\Sigma_{y,\ell}$ and $\Xi_{\ell}$ are Hermitian
and involutory. Therefore, the restriction of Eq. (266) to the left half of
the chain,

$\Gamma_{L}=\Sigma_{y,L/2}\Xi_{L/2}\Gamma_{R}\Xi_{L/2}\Sigma_{y,L/2},$ (267)

confirms the claimed unitary equivalence between $\Gamma_{L}$ and
$\Gamma_{R}$.

###  III.4 Subsystem entropy

Figure 4 in the main text shows the quasiparticle-pair contribution to the
subsystem entropy, which is defined as

$S_{\mathrm{vN},\ell}^{\mathrm{QP}}=S_{\mathrm{vN},\ell}-\frac{\ell}{L}S_{\mathrm{vN}}^{\mathrm{stat}}.$
(268)

In the following, we describe how to calculate the subsystem entropy
$S_{\mathrm{vN},\ell}$ and the statistical entropy
$S_{\mathrm{vN}}^{\mathrm{stat}}=S_{\mathrm{vN},L}$, as well as the PTGGE
prediction that captures the late-time behavior of
$S_{\mathrm{vN},\ell}^{\mathrm{QP}}$. Further, we present a short derivation
that motivates our conjecture for the space-time scaling limit of
$S_{\mathrm{vN},\ell}^{\mathrm{QP}}$ given in Eq. (11) in the main text. In
this section, we restrict ourselves to the case $\delta=0$.

####  III.4.1 Subsystem entropy from the covariance matrix

The von Neumann entropy of a subsystem that consists of $\ell$ contiguous
lattice sites is given by [41]

$S_{\mathrm{vN},\ell}=-\mathop{\mathrm{tr}}(\rho_{\ell}\ln(\rho_{\ell}))=\sum_{l=1}^{\ell}S(\xi_{l}),$
(269)

where

$S(\xi)=-\frac{1+\xi}{2}\ln\\!\left(\frac{1+\xi}{2}\right)-\frac{1-\xi}{2}\ln\\!\left(\frac{1-\xi}{2}\right),$
(270)

and where $\pm i\xi_{l}$ with $0\leq\xi_{l}\leq 1$ and
$l\in\\{1,\dotsc,\ell\\}$ are the eigenvalues of the reduced covariance matrix
$\Gamma_{\ell}$ defined in Eq. (85). When $\delta=0$, the time-evolved
covariance matrix in Eq. (100) reduces to $\Gamma(t)=\Gamma_{1}(t)$, which,
according to Eq. (156), is exponentially small for $\gamma t\gg 1$. Then, to
calculate the subsystem entropy Eq. (269), we can restrict ourselves to the
lowest nontrivial order in the expansion of $S(\xi)$ in powers of $\xi$,

$S(\xi)=\ln(2)-\xi^{2}/2+O(\xi^{4}).$ (271)

In particular, from the rescaled covariance matrix defined in Eq. (238), we
obtain the spectrum of the reduced rescaled covariance matrix:

$\sigma(\Gamma_{\ell}^{\prime}(t))=\\{\pm i\xi_{l}^{\prime}(t)\\}=\\{\pm
ie^{4\gamma t}\xi_{l}(t)\\}.$ (272)

For $\gamma t\gg 1$, the subsystem entropy can be written in terms of the
eigenvalues $\xi_{l}^{\prime}(t)$ as

$S_{\mathrm{vN},\ell}(t)\sim\ell\ln(2)-\frac{1}{2}\sum_{l=1}^{\ell}\xi_{l}(t)^{2}=\ell\ln(2)-\frac{e^{-8\gamma
t}}{2}\sum_{l=1}^{\ell}\xi_{l}^{\prime}(t)^{2}.$ (273)

####  III.4.2 Statistical entropy

Rather than calculating the statistical entropy
$S_{\mathrm{vN}}^{\mathrm{stat}}=S_{\mathrm{vN},L}$ numerically for a finite
system, we perform an analytical calculation in the thermodynamic limit as
detailed in the following: The covariance matrix is block-diagonal in momentum
space, with $4\times 4$ blocks $\mathsf{y}_{k}$ given in Eq. (178). According
to Eqs. (180) and (118), the spectrum of $\mathsf{y}_{k}$ is given by
$\sigma(\mathsf{y}_{k})=\sigma(\gamma_{k})\cup\sigma(\gamma_{-k})=\\{\pm\xi_{+,k},\pm\xi_{-,k}\\}$.
Therefore, in the thermodynamic limit,

$S_{\mathrm{vN}}^{\mathrm{stat}}(t)=\sum_{k\geq
0}\sum_{\sigma=\pm}S(\xi_{\sigma,k}(t))\to
L\int_{0}^{\pi}\frac{dk}{2\pi}\sum_{\sigma=\pm}S(\xi_{\sigma,k}(t)).$ (274)

The eigenvalues $\xi_{\pm,k}(t)$ are given in Eq. (117), where, for
$\delta=0$, $\gamma_{\mathbbm{1},k}(t)=\gamma_{1,\mathbbm{1},k}(t)$ is given
in Eq. (157), and $\xi_{k}(t)=\lvert\boldsymbol{\gamma}_{1,k}(t)\rvert$ can be
obtained by using Eq. (342), which yields

$\xi_{k}(t)=e^{-4\gamma
t}\left\\{1-\frac{\varepsilon_{k}^{2}}{\omega_{k}^{2}}\sin(\Delta\theta_{k})^{2}\left[2\left(1-\cos(2\omega_{k}t)\right)\vphantom{\frac{\varepsilon_{k}^{2}}{\omega_{k}^{2}}}\right.\right.\\\
\left.\left.-\frac{\varepsilon_{k}^{2}}{\omega_{k}^{2}}\left(1-\cos(2\omega_{k}t)\right)^{2}-\sin(2\omega_{k}t)^{2}\right]\right\\}^{1/2},$
(275)

where $\Delta\theta_{k}$ is defined in Eq. (171). At late times, when $\gamma
t\gg 1$, the eigenvalues $\xi_{\pm,k}(t)\sim e^{-4\gamma t}$ are exponentially
small, and we can use the expansion in Eq. (271) to simplify Eq. (274). To
make the exponential time dependence explicit, we define the rescaled
quantities

$\gamma_{\mathbbm{1},k}^{\prime}(t)=e^{4\gamma
t}\gamma_{\mathbbm{1},k}(t),\qquad\xi_{k}^{\prime}(t)=e^{4\gamma
t}\xi_{k}(t),$ (276)

in terms of which we obtain the statistical entropy in the form

$S_{\mathrm{vN}}^{\mathrm{stat}}(t)\sim L\left[\ln(2)-e^{-8\gamma
t}\int_{0}^{\pi}\frac{dk}{2\pi}\left(\gamma_{\mathbbm{1},k}^{\prime}(t)^{2}+\xi_{k}^{\prime}(t)^{2}\right)\right].$
(277)

####  III.4.3 Numerical time evolution

At short times $\gamma t\lesssim 1$, we calculate the covariance matrix for a
finite system numerically according to Eq. (100), and we use Eqs. (269) and
(274) to obtain the data shown in Fig. 4 in the main text. For longer times
$\gamma t\gg 1$, a numerically stable procedure that avoids calculations with
extremely small quantities is to take the difference between Eqs. (273) and
(277), which yields

$S_{\mathrm{vN},\ell}^{\mathrm{QP}}(t)\sim-e^{-8\gamma
t}\left[\frac{1}{2}\sum_{l=1}^{\ell}\xi_{l}^{\prime}(t)^{2}-\ell\int_{0}^{\pi}\frac{dk}{2\pi}\left(\gamma_{\mathbbm{1},k}^{\prime}(t)^{2}+\xi_{k}^{\prime}(t)^{2}\right)\right],$
(278)

where the exponential time dependence is contained in the prefactor that is
determined analytically.

####  III.4.4 Late-time behavior

The late-time limit that is described by the PTGGE can be obtained from Eq.
(278) by properly taking into account dephasing of oscillatory contributions.
Namely, $\xi_{l}^{\prime}(t)$ should be replaced by $\xi_{d,l}^{\prime}$,
where $\pm i\xi_{d,l}^{\prime}$ are the eigenvalues of
$\Gamma_{\mathrm{PTGGE}}^{\prime}$ defined in Eq. (240). In contrast, in the
terms that describe the statistical entropy, oscillatory terms have to be
dropped only after taking the squares of $\gamma_{\mathbbm{1},k}^{\prime}(t)$
and $\xi_{k}^{\prime}(t)$,

$S_{\mathrm{vN},\ell}^{\mathrm{QP}}(t)\sim-e^{-8\gamma
t}\left[\frac{1}{2}\sum_{l=1}^{\ell}\xi_{d,l}^{\prime
2}-\ell\int_{0}^{\pi}\frac{dk}{2\pi}\left(\gamma_{\mathbbm{1},k}^{\prime
2}+\xi_{k}^{\prime 2}\right)_{d}\right],$ (279)

where, according to Eq. (157),

$\left(\gamma_{\mathbbm{1},k}^{\prime
2}\right)_{d}=\frac{3}{2\omega_{k}^{4}}\left[\hat{\mathbf{a}}_{\mathbbm{1},k}\cdot\left(\mathbf{a}_{k}\times\mathbf{b}\right)\right]^{2},$
(280)

and, with Eq. (275),

$\left(\xi_{k}^{\prime
2}\right)_{d}=1+\frac{6\varepsilon_{k}^{2}\gamma^{2}}{\omega_{k}^{4}}\sin(\Delta\theta_{k})^{2}.$
(281)

####  III.4.5 Space-time scaling limit

The quasiparticle picture [41, 42, 43] suggests the following functional form
for the time dependence of the subsystem entropy after a quench in an isolated
integrable system:

$S_{\mathrm{vN},\ell}(t)\sim\int_{-\pi}^{\pi}\frac{dk}{2\pi}\min(2\left\lvert
v_{k}\right\rvert t,\ell)s_{k},$ (282)

where $v_{k}$ is the quasiparticle velocity, and the weight $s_{k}$ determines
the contribution of quasiparticles with momentum $k$ to the entropy. The above
form holds in the space-time scaling limit, $\ell,t\to\infty$ with $\ell/t$
fixed, and does not capture $O(1)$ contributions such as area-law entanglement
of the initial state. To provide a quantitative description of the evolution
of the subsystem entropy for a specific model, concrete expressions for
$v_{k}$ and $s_{k}$ are required. According to Ref. [42], the weight $s_{k}$
can be fixed to reproduce the stationary value of the entropy that is reached
for $t\to\infty$. For free fermionic models, this leads to

$s_{k}=-\langle n_{k}\rangle_{0}\ln(\langle n_{k}\rangle_{0})-\left(1-\langle
n_{k}\rangle_{0}\right)\ln(1-\langle n_{k}\rangle_{0}),$ (283)

where the occupation numbers $\langle
n_{k}\rangle_{0}=\langle\psi_{0}|d_{k}^{\dagger}d_{k}|\psi_{0}\rangle$ of
eigenmodes of the postquench Hamiltonian are determined by the initial state.
Further, $v_{k}$ is given by the group velocity $v_{k}=d\varepsilon_{k}/dk$,
where $\varepsilon_{k}$ is the dispersion relation of quasiparticles.

To derive Eq. (11) of the main text, which we conjecture to generalize Eq.
(282) to PT-symmetric and potentially strongly dissipative open systems, we
first rewrite Eq. (282) in terms of the covariance matrix in the complex
basis. Under the assumption that $\left\lvert v_{k}\right\rvert=\left\lvert
v_{-k}\right\rvert$, which is satisfied for the Kitaev chain, the range of
integration Eq. (282) can be restricted as

$S_{\mathrm{vN},\ell}(t)\sim\int_{0}^{\pi}\frac{dk}{2\pi}\min(2\left\lvert
v_{k}\right\rvert t,\ell)\left(s_{k}+s_{-k}\right).$ (284)

Next, by noting that, for $S(\xi)$ given in Eq. (270),

$s_{k}=S(2\langle n_{k}\rangle_{0}-1),$ (285)

and with $2n_{k}-1=-[d_{k},d_{k}^{\dagger}]$ and $S(\xi)=S(-\xi)$, we obtain

$s_{k}+s_{-k}=\mathop{\mathrm{tr}}(S(g_{\mathrm{GGE},k})),$ (286)

where $g_{\mathrm{GGE},k}$ is defined in Eq. (245), and we use the cyclic
invariance of the trace. We are thus lead to the desired expression for the
subsystem entropy in the space-time scaling limit in terms of the complex
covariance matrix:

$S_{\mathrm{vN},\ell}\sim\int_{0}^{\pi}\frac{dk}{2\pi}\min(2\left\lvert
v_{k}\right\rvert t,\ell)\mathop{\mathrm{tr}}(S(g_{\mathrm{GGE},k})).$ (287)

Further, in analogy to our approach for the subsystem parity in Sec. III.3.3,
to account for modifications of quasiparticle dynamics and statistics due to
strong dissipation, we set $v_{k}=d\omega_{k}/dk$ and we replace
$g_{\mathrm{GGE},k}$ by $g_{\mathrm{PTGGE},k}(t)$ given in Eq. (234). In
particular, to show explicitly the effect of modified quasiparticle
statistics, which are described by the anticommutation relations in Eq. (219),
we use the cyclic invariance of the trace to write

$\mathop{\mathrm{tr}}(S(g_{\mathrm{PTGGE},k}(t)))=\mathop{\mathrm{tr}}\\!\left(S\\!\left(\zeta_{k}(t)f_{k}^{-1}\right)\right).$
(288)

Finally, we subtract the statistical entropy to single out the contribution
due the propagation of pairs of entangled quasiparticles [44, 45, 46]. To that
end, we rewrite the statistical entropy in Eq. (274) by using $S(\xi)=S(-\xi)$
as

$\sum_{\sigma=\pm}S(\xi_{\sigma,k}(t))=\mathop{\mathrm{tr}}(S(\gamma_{k}(t))=\mathop{\mathrm{tr}}(S(g_{k}(t)).$
(289)

In the space-time scaling limit, oscillatory contributions dephase in the
integration over momenta, and, therefore, we have to keep only the
nonoscillatory part, which we denote by
$\mathop{\mathrm{tr}}(S(g_{k}(t))_{d}$. Then, following Ref. [45], we obtain
the quasiparticle-pair contribution to the subsystem entropy in the form given
in Eq. (11) in the main text. For $\gamma t\gg 1$, we can simplify Eq. (11) by
using the expansion of $S(\xi)$ given in Eq. (271). First, for the subsystem
entropy in Eq. (287), and with Eq. (232), we obtain

$\mathop{\mathrm{tr}}\\!\left(S\\!\left(\zeta_{k}(t)f_{k}^{-1}\right)\right)\sim\ln(2)-\frac{e^{-8\gamma
t}}{2}\mathop{\mathrm{tr}}\\!\left(\left(\zeta_{k}^{\prime}f_{k}^{-1}\right)^{2}\right).$
(290)

Further, for the statistical entropy, as in Eq. (279), we find

$\mathop{\mathrm{tr}}(S(g_{k}(t)))_{d}\sim\ln(2)-e^{-8\gamma
t}\left(\gamma_{\mathbbm{1},k}^{\prime 2}+\xi_{k}^{\prime 2}\right)_{d}.$
(291)

By combining these results, we are lead to

$S_{\mathrm{vN},\ell}^{\mathrm{QP}}(t)\sim-e^{-8\gamma
t}\int_{0}^{\pi}\frac{dk}{2\pi}\min(2\left\lvert v_{k}\right\rvert t,\ell)\\\
\times\left(\frac{1}{2}\mathop{\mathrm{tr}}\\!\left(\left(\zeta_{k}^{\prime}f_{k}^{-1}\right)^{2}\right)-\left(\gamma_{\mathbbm{1},k}^{\prime
2}+\xi_{k}^{\prime 2}\right)_{d}\right).$ (292)

This is the form used to obtain the solid lines in Fig. 4 of the main text for
$t/t_{F}\gtrsim 0.5$.

##  IV Universality of PT-symmetric relaxation

Generalized thermalization, i.e., relaxation to the GGE, is expected to occur
universally in the quench dynamics of isolated integrable systems. Is
relaxation to the PTGGE, which we have discussed so far for the specific
example of a driven-dissipative Kitaev chain, an equally universal property
for PT-symmetric and integrable driven-dissipative systems? In the following,
we address this question by studying the quench dynamics of a selection of PT-
symmetric models. We demonstrate relaxation to the PTGGE for two broad classes
of noninteracting fermionic many-body systems, and—for a subset of
observables—also for an interacting spin chain that can be mapped to fermions
with quadratic jump operators. For noninteracting bosons, we show that the
relaxation dynamics is captured by a dephased ensemble that generalizes the
PTGGE for fermions while maintaining the key property of conserving an
extensive amount of information about the initial state. Crucially, the same
mechanism—dephasing within a subspace that is spanned by eigenmodes of the
Liouvillian with equal decay rates—underlies relaxation to both the PTGGE for
fermions and the dephased ensemble for bosons. Therefore, while the PTGGE in
the form presented in the main text is not as general for driven-dissipative