integrable systems as the GGE for isolated integrable systems, our results
indicate that the underlying mechanism is indeed very general.

Before we enter the detailed discussion of other models, let us point out that
relaxation to the PTGGE as described in the main text is not tied to the
particular form of the Kitaev chain and the jump operators discussed there.
Instead, our specific choice of Hamiltonian and jump operators is to be
understood as one particular representative of a broad class of models. To see
why this is the case, we recall the necessary conditions for relaxation to the
PTGGE of the driven-dissipative Kitaev chain: (i) PT symmetry of the
Liouvillian; (ii) the existence of a fully PT-symmetric phase; (iii) for
relaxation to the PTGGE to persist to arbitrary times (for simplicity, we
restrict our discussion of other models to this case), the steady state must
be an infinite-temperature state. But these requirements leave considerable
freedom to modify the model. On the one hand, this concerns the Hamiltonian:
For example, relaxation to the PTGGE persists for long-range hopping and
pairing [47]. On the other hand, also the jump operators can be modified: An
infinite-temperature steady state is reached whenever the jump operators are
Hermitian, i.e., when they take the form given in Eq. (91) with real
coefficients $B_{l,l^{\prime}}$. However, PT symmetry poses restrictions on
the coefficients $B_{l,l^{\prime}}$. For example, for purely local
dissipation, the choices $L_{l}=\sqrt{\gamma}w_{2l-1}$ and
$L_{l}=\sqrt{\gamma}w_{2l}$ are compatible with PT symmetry. Moreover, for
both choices, a fully PT-symmetric phase exists up to a finite critical value
of $\gamma$. These considerations extend straightforwardly to systems in
higher dimensions.

###  IV.1 Fermionic SSH model with loss and gain

Noninteracting fermionic systems can broadly be classified into models with
and without pairing terms in the Hamiltonian, and the Kitaev chain, on which
our discussion has focused so far, is a particular representative of the
former class. In the following, we consider the SSH model [48, 49] as a
paradigmatic example of a one-dimensional tight-binding chain without pairing.
We present a short summary of results that demonstrate relaxation to the
PTGGE. A more detailed discussion will be provided elsewhere [47].

The SSH model describes hopping along a one-dimensional bipartite lattice with
sublattices $A$ and $B$. For a system that consists of $L$ unit cells, the
Hamiltonian reads

$H_{\mathrm{SSH}}=\sum_{l=1}^{L}\left(J_{1}c_{A,l}^{\dagger}c_{B,l}+J_{2}c_{A,l+1}^{\dagger}c_{B,l}+\mathrm{H.c.}\right),$
(293)

where $J_{1},J_{2}\in\mathbb{R}$. We find it convenient to parameterize
$J_{1}$ and $J_{2}$ in terms of the total and relative hopping amplitudes,

$J=\frac{1}{2}\left(J_{1}+J_{2}\right),\qquad\Delta
J=\frac{1}{2}\left(J_{1}-J_{2}\right).$ (294)

The topologically nontrivial phase of the SSH model is realized for $\Delta
J<0$ [5, 6, 7].

Isolated quadratic fermionic systems without pairing terms in the Hamiltonian
do not exhibit anomalous correlations, and, therefore, their dynamics is
captured by a closed equation of motion for the covariance matrix which is
defined as

$G_{l,l^{\prime}}^{s,s^{\prime}}=\langle[c_{s,l},c_{s^{\prime},l^{\prime}}^{\dagger}]\rangle.$
(295)

For such systems, a natural choice of dissipation which results again in a
closed equation of motion for the covariance matrix and can lead to an
infinite-temperature steady state, is loss and gain with jump operators

$L_{l,s,l}=\sqrt{\gamma_{l,s}}c_{s,l},\qquad
L_{g,s,l}=\sqrt{\gamma_{g,s}}c_{s,l}^{\dagger},$ (296)

where $\gamma_{l,s}$ and $\gamma_{g,s}$ are the loss and gain rates on
sublattice $s\in\\{A,B\\}$, respectively. In particular, the steady state is
at infinite temperature if the loss and gain rates are equal,
$\gamma_{l,s}=\gamma_{g,s}$ [47]. Focusing on this case, we set
$\gamma_{l,A}=\gamma_{g,A}=\gamma+\Delta\gamma$, and
$\gamma_{l,B}=\gamma_{g,B}=\gamma-\Delta\gamma$. The time evolution of the
density matrix is then given by a master equation in Lindblad form,

$i\frac{d\rho}{dt}=[H_{\mathrm{SSH}},\rho]+i\left(\mathcal{D}_{l}\rho+\mathcal{D}_{g}\rho\right),$
(297)

with dissipators $\mathcal{D}_{l}$ and $\mathcal{D}_{g}$ that describe loss
and gain, respectively. To state the resulting equation for the covariance
matrix, it is convenient to introduce $2L$ fermionic operators $c_{l}$ such
that $c_{A,l}=c_{2l-1}$ and $c_{B,l}=c_{2l}$, and to write the Hamiltonian in
the form

$H_{\mathrm{SSH}}=\sum_{l,l^{\prime}=1}^{2L}c_{l}^{\dagger}H_{l,l^{\prime}}c_{l^{\prime}},$
(298)

and the jump operators as

$L_{l,l}=\sum_{l^{\prime}=1}^{2L}B_{l,l,l^{\prime}}c_{l^{\prime}},\qquad
L_{g,l}=\sum_{l^{\prime}=1}^{2L}B_{g,l,l^{\prime}}c_{l^{\prime}}^{\dagger}.$
(299)

In terms of the coefficient matrices $B_{l}$ and $B_{g}$, we define the bath
matrices

$M_{l}=B_{l}^{\dagger}B_{l},\qquad
M_{g}=B_{g}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}B_{g}^{*}.$
(300)

Finally, we write the covariance matrix as

$G_{l,l^{\prime}}=\langle[c_{l},c_{l^{\prime}}^{\dagger}]\rangle.$ (301)

Then, the equation of motion for the covariance matrix reads

$\frac{dG}{dt}=-i\left(ZG-GZ^{\dagger}\right)+2\left(M_{l}-M_{g}\right),$
(302)

where $Z=H-i\left(M_{l}+M_{g}\right)$. The matrix $Z$ has PT symmetry, and a
fully PT-symmetric phase is realized for sufficiently small values of $\gamma$
and $\Delta\gamma$ [47]. For equal loss and gain rates, also the bath matrices
are identical, $M_{l}=M_{g}$, and the inhomogeneity in Eq. (302) vanishes.
Consequently, the covariance matrix decays to zero for $t\to\infty$, which
reflects the fact that the steady state has infinite temperature.

Figure 1: Quasiparticle-pair contribution to the entanglement entropy after
quenches to the trivial (green, $\Delta J=0.5J$) and topological (blue,
$\Delta J=-0.5J$) PT-symmetric phases of the SSH model for
$\gamma=\Delta\gamma=0.1J$, and $\ell=20$. The numerical data (dashed lines)
is close to the analytical prediction (solid lines), which is based on a
conjecture in analogy to Eq. (11) of the main text [47]. Inset: For the
trivial quench and at time $t=2t_{F}$, the difference between the numerical
data and the analytical prediction (blue dots) vanishes as $1/\ell$ (orange
line).

In Fig. 1, we illustrate relaxation to the PTGGE with the time evolution of
the quasiparticle-pair contribution to the subsystem entropy following a
quench from the ground state with $\Delta J_{0}=J$. We consider a subsystem
that consist of $\ell$ unit cells. As shown in the inset of the figure, the
numerical data for the quasiparticle-pair entropy converges with subsystem
size as $1/\ell$ to our analytical prediction. Details of the underlying
analytical and numerical computations will be presented in a forthcoming
publication [47].

Finally, let us emphasize that, as we have already discussed for the Kitaev
chain, relaxation to the PTGGE does not require the specific forms of the
Hamiltonian and the jump operators given in Eqs. (293) and (296),
respectively. Instead, both can be modified as long as the requirements
(i)–(iii) stated above are satisfied. Therefore, the behavior shown in Fig. 1
is representative for a whole class of models.

###  IV.2 Bosonic SSH model with loss and gain

Next, we consider a bosonic version of the SSH model with loss and gain. As
anticipated above, the PT-symmetric relaxation dynamics of noninteracting
bosons is described by a dephased ensemble that differs from the PTGGE for
fermions. The nature and origin of the modifications for bosons can be
understood as follows: The Gaussian states which describe the noninteracting
systems we consider are fully determined by two-point correlations. However,
the decomposition of a two-point function such as $\langle
c_{l}c_{l^{\prime}}^{\dagger}\rangle$ into correlation- and statistical
contributions differs for fermions and bosons. In both cases we can write
$c_{l}c_{l^{\prime}}^{\dagger}=\left([c_{l},c_{l^{\prime}}^{\dagger}]+\\{c_{l},c_{l^{\prime}}^{\dagger}\\}\right)\\!/2$.
Fermionic statistics are encoded in the anticommutator,
$\\{c_{l},c_{l^{\prime}}^{\dagger}\\}=\delta_{l,l^{\prime}}$, and correlations
in the expectation value of the commutator which, therefore, features in the
covariance matrix Eq. (301). For the fermionic SSH model, we have chosen equal
loss and gain rates such that the steady state is the infinite-temperature
state, $\rho_{\mathrm{SS}}=\rho_{\infty}=1/2^{2L}$. Therefore, taking the
expectation value of any operator in the steady state reduces to taking the
trace of that operator. But the trace of a commutator vanishes, and,
consequently, the covariance matrix vanishes in the steady state,
$\langle[c_{l},c_{l^{\prime}}^{\dagger}]\rangle_{\mathrm{SS}}=\mathop{\mathrm{tr}}\\!\left([c_{l},c_{l^{\prime}}^{\dagger}]\right)=0$.
In the PT-symmetric phase, the decay of the covariance matrix is described by
a single exponential factor, $G(t)=e^{-4\gamma t}G^{\prime}(t)$, and dephasing
among PT-symmetric modes induces relaxation of $G^{\prime}(t)$ according to
$G^{\prime}(t)\to G^{\prime}_{\mathrm{PTGGE}}$. In contrast, for bosons,
statistics are contained in the commutator,
$[a_{l},a_{l^{\prime}}^{\dagger}]=\delta_{l,l^{\prime}}$, and correlations are
described by the covariance matrix defined as the expectation value of the
anticommutator,

$G_{l,l^{\prime}}=\langle\\{a_{l},a_{l^{\prime}}^{\dagger}\\}\rangle.$ (303)

However, the anticommutator
$\\{a_{l},a_{l}^{\dagger}\\}=1+2a_{l}^{\dagger}a_{l}$ is positive-definite,
and its expectation value is generically nonzero. This leads to a covariance
matrix of the general form

$G(t)=e^{-2\delta t}G^{\prime}(t)+G_{\mathrm{SS}},$ (304)

with a nonvanishing steady-state contribution $G_{\mathrm{SS}}$. While the
fact that this constant term is nonzero for any choice of dissipation
represents a fundamental difference to fermionic systems, the structure of the
time-dependent term is the same for both fermions and bosons: PT symmetry
leads to decay with a single exponential factor, and, also for bosons,
relaxation of $G^{\prime}(t)\to G_{d}^{\prime}$ is induced by dephasing.
Therefore, the late-time relaxation dynamics of a PT-symmetric bosonic system
is described by a dephased ensemble,

$G(t)\sim G_{d}(t)=e^{-2\delta t}G_{d}^{\prime}+G_{\mathrm{SS}}.$ (305)

In analogy to the fermionic PTGGE, information about the initial state is
contained in the bosonic dephased ensemble in the term $G_{d}^{\prime}$.
Additionally, information about the final state enters through both
$G_{d}^{\prime}$ and $G_{\mathrm{SS}}$.

To illustrate relaxation to this dephased ensemble, we study quench dynamics
of a bosonic SSH model that is described by the Hamiltonian

$H_{\mathrm{BSSH}}=\sum_{l=1}^{L}\left(J_{1}a_{A,l}^{\dagger}a_{B,l}+J_{2}a_{A,l+1}^{\dagger}a_{B,l}+\mathrm{H.c.}\right),$
(306)

where $a_{s,l}$ and $a_{s,l}^{\dagger}$ with $s\in\\{A,B\\}$ are bosonic
annihilation and creation operators, respectively, and obey the canonical
commutation relations
$[a_{s,l},a_{s^{\prime},l^{\prime}}^{\dagger}]=\delta_{s,s^{\prime}}\delta_{l,l^{\prime}}$
and $[a_{s,l},a_{s^{\prime},l^{\prime}}]=0$. As in the fermionic case, we
consider local loss and gain, described by the jump operators

$L_{l,s,l}=\sqrt{\gamma_{l,s}}a_{s,l},\qquad
L_{g,s,l}=\sqrt{\gamma_{g,s}}a_{s,l}^{\dagger}.$ (307)

For fermions, we have focused exclusively on the case of equal loss and gain
rates, $\gamma_{l,s}=\gamma_{g,s}$. Here we consider the more general case of
unequal loss and gain rates on both sublattices,
$\gamma_{l/g,A}=\gamma_{l/g}+\Delta\gamma_{l/g}$ and
$\gamma_{l/g,B}=\gamma_{l/g}-\Delta\gamma_{l/g}$. We find it convenient to
express the loss and gain rates in terms of the four parameters $\gamma$,
$\delta$, $\Delta\gamma$, and $\Delta$, which are defined as follows:

$\displaystyle\gamma$ $\displaystyle=\frac{\gamma_{l}+\gamma_{g}}{2},$
$\displaystyle\delta$ $\displaystyle=\gamma_{l}-\gamma_{g}$ (308)
$\displaystyle\Delta\gamma$
$\displaystyle=\frac{\Delta\gamma_{l}+\Delta\gamma_{g}}{2},$
$\displaystyle\Delta$ $\displaystyle=\Delta\gamma_{l}-\Delta\gamma_{g}.$

The equation of motion for the bosonic covariance matrix defined in Eq. (303)
reads [47]

$\frac{dG}{dt}=-i\left(ZG-GZ^{\dagger}\right)+2\left(M_{l}+M_{g}\right),$
(309)

with $Z=H-i\left(M_{l}-M_{g}\right)$ and where the Hamiltonian matrix $H$ and
the bath matrices $M_{l}$ and $M_{g}$ are defined in Eqs. (298) and (300),
respectively. The only but important difference from the fermionic version in
Eq. (302) is the sign of the gain matrix $M_{g}$. Due to this sign change, the
inhomogeneity in Eq. (309), which is the sum of two positive-semidefinite
matrices, is always nonzero for nonvanishing dissipation. In particular, it is
nonzero for equal loss and gain rates, $\delta=\Delta=0$, when the
inhomogeneity in the fermionic equation of motion (302) vanishes. This
reflects the fact that, as explained above, the covariance matix is
generically nonzero in the steady state. Nevertheless, in the fully PT-
symmetric phase, which is realized for sufficiently small values of loss and
gain rates, all eigenvalues $\lambda_{l}$ of $Z$ have constant imaginary part,
which leads to the general form of the covariance matrix given in Eq. (304).
Dephasing among modes with different real parts
$\mathop{\mathrm{Re}}(\lambda_{l})$ induces relaxation to an ensemble that is
determined by $G_{d}(t)$ defined in Eq. (305). As we proceed to show,
relaxation of $G^{\prime}(t)\to G_{d}^{\prime}$ is clearly visible in
quasiparticle-pair contribution to the subsystem entropy.

Figure 2: Quasiparticle-pair contribution to the entanglement entropy after
quenches of the bosonic SSH model with loss and gain to trivial (green,
$\Delta J=0.5J$) and topological (blue $\Delta J=-0.5J$) PT-symmetric phases
for $\gamma=\Delta\gamma=0.15J$, $\delta=\Delta=0.12J$ and $\ell=10$. In the
initial state in Eq. (311), we have set $\mu_{A}=1$, $\mu_{B}=10$. The red
horizontal lines indicate the asymptotic prediction given in Eq. (314).

We thus consider quench dynamics and choose a mixed initial state $\rho_{0}$
as follows: For $J_{2}=0$, the Hamiltonian in Eq. (306) is diagonal when it is
expressed in terms of bosonic modes $b_{s,l}$ that are defined as

$\begin{pmatrix}b_{A,l}\\\
b_{B,l}\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix}\begin{pmatrix}a_{A,l}\\\ a_{B,l}\end{pmatrix}.$ (310)

The mixed initial state $\rho_{0}$ is determined by two parameters $\mu_{A}$
and $\mu_{B}$,

$\rho_{0}=\frac{1}{Z_{0}}e^{-\sum_{s=A,B}\mu_{s}\sum_{l=1}^{L}b_{s,l}^{\dagger}b_{s,l}},$
(311)

where the normalization $Z_{0}$ ensures that
$\mathop{\mathrm{tr}}(\rho_{0})=1$.

Further, the loss and gain rates are chosen so as to lead to a particularly
simple steady state. The steady-state occupation of a single bosonic lattice
site is determined by the ratio between loss and gain rates. Therefore, to
achieve a spatially uniform distribution of particles, we choose equal ratios
on both sublattices, $\gamma_{l,A}/\gamma_{g,A}=\gamma_{l,B}/\gamma_{g,B}$ or,
equivalently, $\gamma/\Delta\gamma=\delta/\Delta$. The resulting uniform
distribution is also stationary under unitary evolution generated by the SSH
Hamiltonian in Eq. (306). The covariance in the steady state is thus given by
$G_{\mathrm{SS}}=\left(\delta/\gamma\right)\mathbbm{1}$.

As for the fermionic case, the entanglement entropy can be calculated from the
covariance matrix [50, 51]. In particular, for a system without anomalous
correlations, the von Neumann entanglement entropy of a subsystem that
consists of $\ell$ unit cells is given by

$S_{\mathrm{vN},\ell}=\sum_{l=1}^{2\ell}S(\xi_{l}),$ (312)

where $\xi_{l}\geq 1$ are the $2\ell$ eigenvalues of the covariance matrix
$G_{\ell}$ restricted to a subsystem that consists of $\ell$ unit cells, and

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

The function $S(\xi)$ for bosons differs from the fermionic version given in
Eq. (270) only by the sign of the first term. A procedure similar to that
described in Sec. III.4 for fermions leads to the following result for the
late-time asymptotics of the quasiparticle-pair contribution to the subsystem
entropy:

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

where $S_{2}=d^{2}S/d\xi^{2}$ is the second derivative of Eq. (313), evaluated
at $\xi=\delta/\gamma$, and $\xi_{d,l}^{\prime}$ are the eigenvalues of
$G_{d}^{\prime}$ in Eq. (305). Further, $G^{\prime}(t)$ in Eq. (304) is a
block Toeplitz matrix with $2\times 2$ blocks $g_{l}^{\prime}(t)$, which we
decompose in momentum space as
$g_{k}^{\prime}(t)=g_{\mathbbm{1},k}^{\prime}(t)\mathbbm{1}+\mathbf{g}_{k}^{\prime}(t)\cdot\boldsymbol{\sigma}$.
This decomposition defines both $g_{\mathbbm{1},k}^{\prime}(t)$ and
$\xi_{k}^{\prime}(t)=\left\lvert\mathbf{g}_{k}^{\prime}(t)\right\rvert$. As in
Eq. (279), the subscript “d” in the last term in Eq. (314) indicates that
oscillating terms should be dropped after taking the square.

In Fig. 2, the quasiparticle-pair contribution to the entanglement entropy is
shown for quenches to the trivial and topological PT-symmetric phases. The
characteristic signatures of the quasiparticle picture in the fermionic model
transfer well to the bosonic case, consequently reinforcing particular
universal aspects of a PT-symmetric ensemble description. Accordingly, after
rescaling with $e^{4\delta t}$, the bosonic entropy shows initial ballistic
growth, followed by saturation roughly at the characteristic Fermi time
$t_{F}$, which depends on the modified quasiparticle velocity. The late-time
asymptotic saturation of the entropy is well predicted by the dephased
ensemble through Eq. (314).

###  IV.3 Dimerized XY model with dephasing

In the study of quench dynamics of isolated integrable systems, an important
role is played by one-dimensional spin chains, or, equivalently, lattice
models of hard-core bosons. Paradigmatic examples are the XX model [52, 9, 10,
11, 53] and the transverse field Ising model [52, 9, 10, 11, 53]. For these
models of interacting spins or bosons, analytical progress is facilitated by
mapping them to noninteracting fermions through the Jordan-Wigner
transformation. For example, pertinent GGEs can be expressed in terms of
occupation numbers of free fermionic modes. In contrast, for interacting
systems that cannot be mapped to noninteracting systems, the construction of
the GGE is much more complicated [32]. Therefore, to address the question
whether our results for PT-symmetric relaxation of driven-dissipative
integrable systems also apply to interacting systems, attempting to generalize
these results to interacting systems that can be mapped to noninteracting
fermions [53] presents itself as a natural first step. However, applying the
Jordan-Wigner transformation to driven-dissipative spin chains with jump
operators that are linear in spin operators does not lead to noninteracting
fermionic models with jump operators that are linear in fermionic operators.
Instead, the Jordan-Wigner transformation in Eq. (77) maps spin flips
corresponding to the jump operators $\sigma_{l}^{\pm}$ to nonlocal products of
fermionic operators. Dephasing, which is described by local Hermitian jump
operators,

$L_{l}=\frac{\sqrt{\gamma}}{2}\sigma^{z}_{l},$ (315)

takes a much simpler but still quadratic form in fermionic operators. Thus, in
the following, we will study quench dynamics of a dimerized XY spin chain [54]
with dephasing as a canonical and experimentally relevant form of dissipation,
which occurs naturally and can also be induced deliberately, e.g., in quantum
simulations of the XY model with trapped ions [55]. We emphasize that there is
a fundamental difference between master equations with linear and quadratic
Hermitian jump operators such as the driven-dissipative Kitaev chain and the
XY model with dephasing, respectively: In the former case, the Gaussianty of
an initially Gaussian state is preserved; in the latter case, an initially
Gaussian state evolves into a mixture of Gaussian states [56], and, therefore,
nontrivial higher-order connected correlations are generated. Nevertheless, a
master equation with quadratic Hamiltonians and quadratic Hermitian jump
operators still leads to a closed equation of motion for the covariance
matrix. As we will show, the relaxation of the covariance matrix—and,
therefore, of all spin observables that can be expressed in terms of the
covariance matrix—is described by a PTGGE even though the state of the system
is not Gaussian.

The dimerized XY spin chain [54] is described by the following Hamiltonian:

$H_{\mathrm{XY}}=\frac{1}{2}\sum_{l=1}^{2L}\left(J-\left(-1\right)^{l}\Delta
J\right)\left(\sigma^{x}_{l}\sigma^{x}_{l+1}+\sigma^{y}_{l}\sigma^{y}_{l+1}\right),$
(316)

where $\sigma_{l}^{x,y,z}$ are Pauli matrices that act on the spin at site $l$
with $l\in\\{1,\dotsc,2L\\}$. We choose the spatially uniform component of the
XY couplings as antiferromagnetic, $J>0$. The spin-Peierls term with
dimerization strength $\Delta J$ describes a modulation of the XY couplings
across successive bonds between the values $J\pm\Delta J$ [54]. This
modulation results in dimer order in the ground state, which is detected by a
nonvanishing expectation value of the operator [57, *Haldane1982a]

$d_{l}=\sigma^{+}_{l}\sigma^{-}_{l+1}-\sigma^{+}_{l-1}\sigma^{-}_{l}+\mathrm{H.c.}$
(317)

In particular, for $\Delta J=J$, the spin chain is fully dimerized and the
ground state is given by a product of singlets across the odd bonds, which
leads to the dimerization order parameter taking the value $\langle
d_{2l-1}\rangle_{0}=-1=-\langle d_{2l}\rangle_{0}$.

With dephasing as given in Eq. (315), the master equation for the XY chain
reads

$i\frac{d\rho}{dt}=\mathcal{L}_{\mathrm{XY}}\rho=[H_{\mathrm{XY}},\rho]-i\sum_{l=1}^{2L}[L_{l},[L_{l},\rho]].$
(318)

Since the jump operators are Hermitian, the steady state is again at infinite
temperature. Further, the Liouvillian $\mathcal{L}_{\mathrm{XY}}$ of the XY
model with dephasing has PT symmetry [59].

To numerically study quench dynamics, we switch to a representation of the XY
model in terms of fermions. This representation can be obtained through the
Jordan-Wigner transformation in Eq. (77), which maps the XY Hamiltonian in Eq.
(316) to

$H_{\mathrm{XY}}=-\sum_{l=1}^{2L}\left(J-\left(-1\right)^{l}\Delta
J\right)\left(c_{l}^{\dagger}c_{l+1}+\mathrm{H.c.}\right).$ (319)

That is, the Hamiltonian takes the form of the SSH model in Eq. (293), with
hopping amplitudes $J_{1,2}=-\left(J\pm\Delta J\right).$ Further, the jump
operators for dephasing in Eq. (315) are quadratic in fermionic operators,

$L_{l}=\sqrt{\gamma}c_{l}^{\dagger}c_{l},$ (320)

where we omit a constant contribution that cancels due to the commutator
structure in Eq. (318).

Figure 3: Dimerization order parameter $\langle d_{2l-1}\rangle$ across odd
bonds. The initial value for a fully dimerized chain is $\langle
d_{2l-1}\rangle_{0}=-1$. After rescaling with $e^{2\gamma t}$, the
dimerization order parameter approaches a stationary value that is determined
by the PTGGE and indicated by horizontal lines in the figure.

In analogy to Eq. (301) for the SSH model with linear dissipation, we define
the covariance matrix as

$G_{l,l^{\prime}}=\langle[c_{l},c_{l^{\prime}}^{\dagger}]\rangle,$ (321)

where $l,l^{\prime}\in\\{1,\dotsc,2L\\}$. To formulate the equation of motion
for the covariance matrix, we write the Hamiltonian and the jump operators in
the form

$H_{\mathrm{XY}}=\sum_{l,l^{\prime}=1}^{2L}c_{l}^{\dagger}H_{l,l^{\prime}}c_{l^{\prime}},\quad
L_{l}=\sum_{m,m^{\prime}=1}^{2L}c_{m}^{\dagger}B_{l,m,m^{\prime}}c_{m^{\prime}},$
(322)

where, for dephasing as described in Eq. (320), the elements of the matrix
$B_{l}$ are given by
$B_{l,m,m^{\prime}}=\sqrt{\gamma}\delta_{l,m}\delta_{l,m^{\prime}}$. The
covariance matrix evolves in time according to

$\frac{dG}{dt}=-i[H,G]-\sum_{l=1}^{L}[B_{l},[B_{l},G]].$ (323)

To illustrate the effect of dephasing most transparently, we decompose the
covariance matrix into diagonal and offdiagonal elements as
$G=G_{\mathrm{OD}}+G_{\mathrm{OD}}$. Then, Eq. (323) reduces to

$\frac{dG}{dt}=-i[H,G]-2\gamma G_{\mathrm{OD}},$ (324)

which shows that dephasing affects only offdiagonal elements. This equation
can be simplified further by taking advantage of the inversion symmetry of the
SSH model: One can show that for translationally invariant and inversion-
symmetric initial states, the diagonal component $G_{\mathrm{D}}$ vanishes at
all times [47]. This applies, in particular, to the ground state of the SSH
Hamiltonian. Then, with $G=G_{\mathrm{OD}}$, we obtain

$\frac{dG}{dt}=-i\left(ZG-GZ^{\dagger}\right),$ (325)

where $Z=H-i\gamma\mathbbm{1}$; that is, time evolution with dephasing differs
from the dynamics of isolated system only by an overall decay of $G$ at the
rate $2\gamma$. Consequently, the PTGGE which describes relaxation of the
covariance matrix and, therefore, of all expectation values of spin
observables that can be expressed in terms of the covariance matrix, is
obtained simply by multiplying the Lagrange multipliers in the GGE by
$e^{-2\gamma t}$. In contrast to the combination of loss and gain discussed in
Sec. IV.1, the dynamics and statistics of elementary excitations are not
affected by dephasing.

Figure 4: Quasiparticle-pair contribution to the subsystem entropy, calculated
from the covariance matrix and thus neglecting deviations from Gaussianity,
after quenches of the dimerized XY model with dephasing for $\Delta J=0.3J$
and and $\ell=20$. For $\gamma=0$, the quasiparticle-pair entropy reduces to
the entanglement entropy, and its stationary value is determined by the
conventional GGE; for any $\gamma>0$, the rescaled quasiparticle-pair entropy
approaches a stationary value that is independent of $\gamma$ and determined
by the PTGGE. The PTGGE prediction is reached at $\gamma t\gg 1$ for very
small values of $\gamma$ with $\gamma<1/t_{F}\approx 0.07J$. In contrast, for
$\gamma>1/t_{F}$, the curves collapse, and converge towards the PTGGE value at
$t\gtrsim t_{F}$. The value of $t_{F}$ does not depend on $\gamma$, which
reflects the fact that dephasing does not modify the dynamics of quasiparticle
excitations.

To illustrate the resulting dynamics, we consider quenches from the ground
state for $\Delta J_{0}=J$ with covariance matrix

$G_{0}=\bigoplus_{l=1}^{L}\sigma_{x},$ (326)

to postquench values $\Delta J\neq\Delta J_{0}$. Relaxation of the covariance
matrix to the PTGGE is shown in Fig. 3 with the time evolution of the
dimerization order parameter Eq. (317), which can be expressed as $\langle
d_{l}\rangle=\mathop{\mathrm{Re}}(G_{l-1,l}-G_{l,l+1})$. As stated below Eq.
(317), the dimerization order parameter takes the value $\langle
d_{2l-1}\rangle_{0}=-1$ in the initial ground state. Due to dephasing, it
decays exponentially, but after rescaling with $e^{2\gamma t}$, it approaches
a stationary value as predicted by the PTGGE.

Since the state of the system after the quench is not Gaussian, the
entanglement entropy is not fully determined by the covariance matrix.
However, it is still informative to neglect deviations from Gaussianity and to
study the time evolution of the quasiparticle-pair entropy computed from the
covariance matrix alone [44]. In Fig. 4, we show the time evolution of the
rescaled quasiparticle-pair entropy, which is defined here as $e^{4\gamma
t}S_{\mathrm{vN},\ell}^{\mathrm{QP}}(t)$ and can be calculated from the
covariance matrix as outlined in Sec. IV.1. The fact that dephasing does not
affect the dynamics and statistics of quasiparticle excitations has two
notable consequences that can be seen in the figure: First, the PTGGE
prediction for the rescaled quasiparticle-pair entropy is independent of the
value of $\gamma$; and second, also the Fermi time, which marks the transition
from the initial linear increase of the rescaled entropy to the much slower
approach to the stationary value, does not depend on $\gamma$.

*

##  Appendix A Time evolution induced by a non-Hermitian $2\times 2$
Hamiltonian

In this short technical appendix, we summarize some relations that allow us to
determine the time evolution of the components of the covariance matrix in Eq.
(112). We consider a non-Hermitian $2\times 2$ Hamiltonian $x$ that is given
by

$x=x_{\mathbbm{1}}\mathbbm{1}+\mathbf{x}\cdot\boldsymbol{\sigma},\qquad\mathbf{x}=\mathbf{a}+i\mathbf{b},$
(327)

where $\mathbf{a},\mathbf{b}\in\mathbb{R}^{3}$. Further, we define

$\omega=\sqrt{\mathbf{x}\cdot\mathbf{x}}=\sqrt{\varepsilon^{2}+i2\mathbf{a}\cdot\mathbf{b}-b^{2}},$
(328)

where $\varepsilon=\left\lvert\mathbf{a}\right\rvert$, and
$b=\left\lvert\mathbf{b}\right\rvert$. In general, $\omega\in\mathbb{C}$.
Finally, we define a vector $\hat{\mathbf{x}}=\mathbf{x}/\omega$. In terms of
these quantities, the exponential of $x$ an be written as

$e^{-ixt}=e^{-ix_{\mathbbm{1}}t}\left(\cos(\omega t)\mathbbm{1}-i\sin(\omega
t)\hat{\mathbf{x}}\cdot\boldsymbol{\sigma}\right).$ (329)

Using this relation, and with

$\boldsymbol{\sigma}\left(\mathbf{x}\cdot\boldsymbol{\sigma}\right)=\mathbf{x}\mathbbm{1}+i\mathbf{x}\times\boldsymbol{\sigma},$
(330)

it is straightforward to show that

$\boldsymbol{\tau}=e^{-ixt}\boldsymbol{\sigma}e^{ix^{\dagger}t}=e^{2\mathop{\mathrm{Im}}(x_{\mathbbm{1}})t}\left\\{\left\lvert\cos(\omega
t)\right\rvert^{2}\boldsymbol{\sigma}\right.\\\ -i\sin(\omega
t)\cos(\omega^{*}t)\left(\hat{\mathbf{x}}-i\hat{\mathbf{x}}\times\boldsymbol{\sigma}\right)\\\
+i\sin(\omega^{*}t)\cos(\omega
t)\left(\hat{\mathbf{x}}^{*}+i\hat{\mathbf{x}}^{*}\times\boldsymbol{\sigma}\right)\\\
\left.+\left\lvert\sin(\omega
t)\right\rvert^{2}\left[\left(\hat{\mathbf{x}}\cdot\boldsymbol{\sigma}\right)\hat{\mathbf{x}}^{*}+\hat{\mathbf{x}}\left(\hat{\mathbf{x}}^{*}\cdot\boldsymbol{\sigma}\right)-i\hat{\mathbf{x}}\times\hat{\mathbf{x}}^{*}-\left(\hat{\mathbf{x}}\cdot\hat{\mathbf{x}}^{*}\right)\boldsymbol{\sigma}\right]\right\\}.$
(331)

Next, we consider the scalar product of $\boldsymbol{\tau}$ and
$\mathbf{a}_{0}\in\mathbb{R}^{3}$, which we denote by $\gamma$:

$\gamma=\mathbf{a}_{0}\cdot\boldsymbol{\tau}=\gamma_{\mathbbm{1}}\mathbbm{1}+\boldsymbol{\gamma}\cdot\boldsymbol{\sigma}.$
(332)

We obtain

$\gamma_{\mathbbm{1}}=2e^{2\mathop{\mathrm{Im}}(x_{\mathbbm{1}})t}\left[\vphantom{\left\lvert\frac{\sin(\omega
t)}{\omega}\right\rvert^{2}}\mathop{\mathrm{Im}}\\!\left(\sin(\omega
t)\cos(\omega^{*}t)\mathbf{a}_{0}\cdot\hat{\mathbf{x}}\right)\right.\\\
\left.-\left\lvert\frac{\sin(\omega
t)}{\omega}\right\rvert^{2}\mathbf{a}_{0}\cdot\left(\mathbf{a}\times\mathbf{b}\right)\right],$
(333)

and

$\boldsymbol{\gamma}=e^{2\mathop{\mathrm{Im}}(x_{\mathbbm{1}})t}\left[\left(\left\lvert\cos(\omega
t)\right\rvert^{2}-\left\lvert\sin(\omega
t)\right\rvert^{2}\hat{\mathbf{x}}\cdot\hat{\mathbf{x}}^{*}\right)\mathbf{a}_{0}\right.\\\
-2\mathop{\mathrm{Re}}\\!\left(\sin(\omega
t)\cos(\omega^{*}t)\mathbf{a}_{0}\times\hat{\mathbf{x}}\right)\\\
\left.+2\left\lvert\sin(\omega
t)\right\rvert^{2}\mathop{\mathrm{Re}}\\!\left(\left(\mathbf{a}_{0}\cdot\hat{\mathbf{x}}\right)\hat{\mathbf{x}}^{*}\right)\right].$
(334)

In the main text, we are primarily concerned with the case
$\mathbf{a}\cdot\mathbf{b}=0$ and $\varepsilon>b$, such that

$\omega=\sqrt{\varepsilon^{2}-b^{2}}\in\mathbb{R}_{>0},$ (335)

and $\mathbf{a}_{0}\cdot\mathbf{b}=0$. These conditions lead to the following
simplifications of the above expressions for $\gamma_{\mathbbm{1}}$ and
$\boldsymbol{\gamma}$:

$\gamma_{\mathbbm{1}}=-\frac{e^{2\mathop{\mathrm{Im}}(x_{\mathbbm{1}})t}}{\omega^{2}}\left(1-\cos(2\omega
t)\right)\mathbf{a}_{0}\cdot\left(\mathbf{a}\times\mathbf{b}\right),$ (336)

and

$\boldsymbol{\gamma}=e^{2\mathop{\mathrm{Im}}(x_{\mathbbm{1}})t}\\\
\times\left[\mathbf{a}_{0}-\frac{\varepsilon^{2}}{\omega^{2}}\left(1-\cos(2\omega
t)\right)\mathbf{a}_{\perp}+\frac{\varepsilon}{\omega}\sin(2\omega
t)\mathbf{a}_{o}\right],$ (337)

where, with $\hat{\mathbf{a}}=\mathbf{a}/\varepsilon$, the parallel,
perpendicular, and out-of-plane components are given by

$\mathbf{a}_{\parallel}=\left(\mathbf{a}_{0}\cdot\hat{\mathbf{a}}\right)\hat{\mathbf{a}},\qquad\mathbf{a}_{\perp}=\mathbf{a}_{0}-\mathbf{a}_{\parallel},\qquad\mathbf{a}_{o}=-\mathbf{a}_{0}\times\hat{\mathbf{a}}.$
(338)

Finally, we calculate the norm of the vector $\boldsymbol{\gamma}$, which we
denote by $\xi=\left\lvert\boldsymbol{\gamma}\right\rvert$. To that end, we
use that $\mathbf{a}_{\parallel}$, $\mathbf{a}_{\perp}$, and $\mathbf{a}_{o}$
are mutually orthogonal, and that

$\mathbf{a}_{0}\cdot\mathbf{a}_{\perp}=\left\lvert\mathbf{a}_{\perp}\right\rvert^{2}=\left\lvert\mathbf{a}_{o}\right\rvert^{2}=a_{0}^{2}-\left(\mathbf{a}_{0}\cdot\hat{\mathbf{a}}\right)^{2}.$
(339)

We further assume that $\mathbf{a}_{0}=\hat{\mathbf{a}}_{0}$ is a unit vector,
such that $a_{0}=\left\lvert\mathbf{a}_{0}\right\rvert=1$. Then,

$\hat{\mathbf{a}}_{0}\cdot\mathbf{a}_{\perp}=\left\lvert\mathbf{a}_{\perp}\right\rvert^{2}=\left\lvert\mathbf{a}_{o}\right\rvert^{2}=1-\left(\hat{\mathbf{a}}_{0}\cdot\hat{\mathbf{a}}\right)^{2}=\sin(\Delta\theta)^{2},$
(340)

where $\Delta\theta$ is the angle between $\hat{\mathbf{a}}_{0}$ and
$\hat{\mathbf{a}}$,

$\cos(\Delta\theta)=\hat{\mathbf{a}}_{0}\cdot\hat{\mathbf{a}},$ (341)

and we obtain

$\xi^{2}=e^{4\mathop{\mathrm{Im}}(x_{\mathbbm{1}})t}\left\\{1-\frac{\varepsilon^{2}}{\omega^{2}}\sin(\Delta\theta)^{2}\left[2\left(1-\cos(2\omega
t)\right)\vphantom{\frac{\varepsilon^{2}}{\omega^{2}}}\right.\right.\\\
\left.\left.-\frac{\varepsilon^{2}}{\omega^{2}}\left(1-\cos(2\omega
t)\right)^{2}-\sin(2\omega t)^{2}\right]\right\\}.$ (342)

## References

  * Wimmer [2012] M. Wimmer, Algorithm 923, ACM Trans. Math. Softw. 38, 1 (2012).
  * Surace and Tagliacozzo [2022] J. Surace and L. Tagliacozzo, Fermionic Gaussian states: An introduction to numerical approaches, SciPost Phys. Lect. Notes 54, (2022).
  * Kraus  _et al._ [2009] C. V. Kraus, M. M. Wolf, J. I. Cirac, and G. Giedke, Pairing in fermionic systems: A quantum-information perspective, Phys. Rev. A 79, 012306 (2009).
  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Phys. Rev. B 55, 1142 (1997).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Leumer  _et al._ [2021] N. Leumer, M. Grifoni, B. Muralidharan, and M. Marganska, Linear and nonlinear transport across a finite Kitaev chain: An exact analytical study, Phys. Rev. B 103, 165432 (2021).
  * Calabrese  _et al._ [2011] P. Calabrese, F. H. L. Essler, and M. Fagotti, Quantum Quench in the Transverse-Field Ising Chain, Phys. Rev. Lett. 106, 227203 (2011).
  * Calabrese  _et al._ [2012a] P. Calabrese, F. H. L. Essler, and M. Fagotti, Quantum quench in the transverse field Ising chain: I. Time evolution of order parameter correlators, J. Stat. Mech. Theory Exp. 2012, P07016 (2012a).
  * Calabrese  _et al._ [2012b] P. Calabrese, F. H. L. Essler, and M. Fagotti, Quantum quenches in the transverse field Ising chain: II. Stationary state properties, J. Stat. Mech. Theory Exp. 2012, P07022 (2012b).
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Über das Paulische Äquivalenzverbot, Zeitschrift für Phys. 47, 631 (1928).
  * Budich and Heyl [2016] J. C. Budich and M. Heyl, Dynamical topological order parameters far from equilibrium, Phys. Rev. B 93, 085416 (2016).
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Statistics of the two-dimensional ferromagnet. part i, Phys. Rev. 60, 252 (1941).
  * Fisher [1995] D. S. Fisher, Critical behavior of random transverse-field ising spin chains, Phys. Rev. B 51, 6411 (1995).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Ann. Phys. (N. Y). 16, 407 (1961).
  * Barouch  _et al._ [1970] E. Barouch, B. M. McCoy, and M. Dresden, Statistical Mechanics of the XY Model. I, Phys. Rev. A 2, 1075 (1970).
  * Barouch and McCoy [1971a] E. Barouch and B. M. McCoy, Statistical Mechanics of the $XY$ Model. II. Spin-Correlation Functions, Phys. Rev. A 3, 786 (1971a).
  * Barouch and McCoy [1971b] E. Barouch and B. M. McCoy, Statistical Mechanics of the XY Model. III, Phys. Rev. A 3, 2137 (1971b).
  * Gorini [1976] V. Gorini, Completely positive dynamical semigroups of N-level systems, J. Math. Phys. 17, 821 (1976).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * Prosen and Ilievski [2011] T. Prosen and E. Ilievski, Nonequilibrium Phase Transition in a Periodically Driven $XY$ Spin Chain, Phys. Rev. Lett. 107, 060403 (2011).
  * Ornigotti and Szameit [2014] M. Ornigotti and A. Szameit, Quasi $\mathcal{PT}$-symmetry in passive photonic lattices, J. Opt. 16, 065501 (2014).
  * Joglekar and Harter [2018] Y. N. Joglekar and A. K. Harter, Passive parity-time-symmetry-breaking transitions without exceptional points in dissipative photonic systems, Photonics Res. 6, A51 (2018).
  * Roccati  _et al._ [2022] F. Roccati, G. M. Palma, F. Ciccarello, and F. Bagarello, Non-Hermitian Physics and Master Equations, Open Syst. Inf. Dyn. 29, 2250004 (2022).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and Topology in Non-Hermitian Physics, Phys. Rev. X 9, 041015 (2019).
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real Spectra in Non-Hermitian Hamiltonians Having $\mathcal{PT}$ Symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * Bender [2007] C. M. Bender, Making sense of non-Hermitian Hamiltonians, Reports Prog. Phys. 70, 947 (2007).
  * Prosen [2008] T. Prosen, Third quantization: a general method to solve master equations for quadratic open Fermi systems, New J. Phys. 10, 043026 (2008).
  * Prosen [2010] T. Prosen, Spectral theorem for the Lindblad equation for quadratic open fermionic systems, J. Stat. Mech. Theory Exp. 2010, P07020 (2010).
  * Sayyad  _et al._ [2021] S. Sayyad, J. Yu, A. G. Grushin, and L. M. Sieberer, Entanglement spectrum crossings reveal non-Hermitian dynamical topology, Phys. Rev. Res. 3, 033022 (2021).
  * Essler and Fagotti [2016] F. H. L. Essler and M. Fagotti, Quench dynamics and relaxation in isolated integrable quantum spin chains, J. Stat. Mech. Theory Exp. 2016, 064002 (2016).
  * Barthel and Schollwöck [2008] T. Barthel and U. Schollwöck, Dephasing and the Steady State in Quantum Many-Particle Systems, Phys. Rev. Lett. 100, 100601 (2008).
  * Peschel [2003] I. Peschel, Calculation of reduced density matrices from correlation functions, Journal of Physics A: Mathematical and General 36, L205 (2003).
  * Peschel and Eisler [2009] I. Peschel and V. Eisler, Reduced density matrices and entanglement entropy in free lattice models, J. Phys. A Math. Theor. 42, 504003 (2009).
  * Vidal  _et al._ [2003] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Entanglement in Quantum Critical Phenomena, Phys. Rev. Lett. 90, 227902 (2003).
  * Gardiner and Zoller [2004] C. W. Gardiner and P. Zoller,  _Quantum Noise_, 3rd ed., Springer series in synergetics, Vol. 56 (Springer, Berlin Heidelberg, 2004).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn,  _Quantum Measurement and Control_, 1st ed. (Cambridge University Press, Cambridge, 2009).
  * Note [1] More 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].
  * Note [2] For 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$.
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems, J. Stat. Mech. Theory Exp. 2005, P04010 (2005).
  * Alba and Calabrese [2017] V. Alba and P. Calabrese, Entanglement and thermodynamics after a quantum quench in integrable systems, Proc. Natl. Acad. Sci. 114, 7947 (2017).
  * Alba and Calabrese [2018] V. Alba and P. Calabrese, Entanglement dynamics after quantum quenches in generic integrable systems, SciPost Phys. 4, 017 (2018).
  * Alba and Carollo [2021] V. Alba and F. Carollo, Spreading of correlations in Markovian open quantum systems, Phys. Rev. B 103, L020302 (2021).
  * Carollo and Alba [2022] F. Carollo and V. Alba, Dissipative quasiparticle picture for quadratic Markovian open quantum systems, Phys. Rev. B 105, 144305 (2022).
  * Alba and Carollo [2022] V. Alba and F. Carollo, Hydrodynamics of quantum entropies in Ising chains with linear dissipation, J. Phys. A Math. Theor. 55, 074002 (2022).
  * [47] E. Starchl and L. M. Sieberer, in preparation.
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * Heeger  _et al._ [1988] A. J. Heeger, S. Kivelson, J. R. Schrieffer, and W. P. Su, Solitons in conducting polymers, Rev. Mod. Phys. 60, 781 (1988).
  * Simon  _et al._ [1994] R. Simon, N. Mukunda, and B. Dutta, Quantum-noise matrix for multimode systems: U(n) invariance, squeezing, and normal forms, Phys. Rev. A 49, 1567 (1994).
  * Adesso and Illuminati [2007] G. Adesso and F. Illuminati, Entanglement in continuous-variable systems: recent advances and current perspectives, J. Phys. A Math. Theor. 40, 7821 (2007).
  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a Completely Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics of the Highly Excited States of 1D Lattice Hard-Core Bosons, Phys. Rev. Lett. 98, 050405 (2007).
  * Vidmar and Rigol [2016] L. Vidmar and M. Rigol, Generalized Gibbs ensemble in integrable lattice models, J. Stat. Mech. Theory Exp. 2016, 064007 (2016).
  * Giamarchi [2003] T. Giamarchi,  _Quantum Physics in One Dimension_, International Series of Monogr (Oxford University Press, Oxford, 2003).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Horstmann  _et al._ [2013] B. Horstmann, J. I. Cirac, and G. Giedke, Noise-driven dynamics and phase transitions in fermionic systems, Phys. Rev. A 87, 012108 (2013).
  * Haldane [1982a] F. D. M. Haldane, Spontaneous dimerization in the $S=\tfrac{1}{2}$ Heisenberg antiferromagnetic chain with competing interactions, Phys. Rev. B 25, 4925 (1982a).
  * Haldane [1982b] F. D. M. Haldane, Erratum: Spontaneous dimerization in the $S=\tfrac{1}{2}$ Heisenberg antiferromagnetic chain with competing interactions, Phys. Rev. B 26, 5257 (1982b).
  * Prosen [2012] T. Prosen, Generic examples of $\mathbb{PT}$-symmetric qubit (spin-1/2) Liouvillian dynamics, Phys. Rev. A 86, 044103 (2012).


* Hasso-Plattner Institute (2020) Hasso-Plattner Institute (2020), Hasso-Plattner Institute.

## About the authors

Kiev Gama is an associate professor of computer science at the Universidade
Federal de Pernambuco (UFPE), Brazil. He received his PhD in computer science
from the University of Grenoble, France. His research interests include
several aspects of time-bounded collaborative events (e.g., hackathons, game
jams) within contexts such as open innovation and education.

George Valença is a professor of software engineering at Universidade Federal
Rural de Pernambuco (UFRPE), Brazil. He received an MSc and a PhD in Computer
Science from Universidade Federal de Pernambuco (UFPE), Brazil. He researches
on open innovation and business process management approaching topics such as
power relationships, data protection and corporate hackathons.

Pedro Alessio is a professor at the Department of Graphic Expression at
Universidade Federal de Pernambuco (UFPE), Brazil. He holds a bachelor’s
degree in Industrial Design from the Federal University of Pernambuco, a
master’s degree and a PhD in Informatics from the Conservatoire National des
arts et métiers (CNAM), France.

Rafael Formiga holds Master’s degree in Design from the Universidade Federal
de Pernambuco (UFPE), Brazil. He is a product designer at VTEX. His main
research focuses on design methods, product design and UX research. He works
with digital design education, building bridges between academia and industry.

André Neves is an associate professor at the Universidade Federal de
Pernambuco (UFPE), Brazil. He has a bachelor’s degree in Industrial Design
from the Universidade Federal da Paraíba (UFPB), and a Master’s and PhD in
Computer Science from UFPE. He is experienced with the Design of computer
systems, working mainly in the investigation, development and application of
methods and design techniques as an instrument of innovation.

Nycolas Lacerda graduated in computer science from the Universidade Federal
Rural de Pernambuco (UFRPE), Brazil. His main research focus is on corporate
hackathon organizations, having participated in a research group ”Innovation
for Business and Software Evolution - IBSE” as a student researcher. He
currently works as an IT consultant.

##  Appendix A Interview Guide

  1. 1.

Demographic Information

     * •

Age

     * •

Educational background

     * •

Time of experience in the area

     * •

Experience with hackathons (total participations, prizes won)

  2. 2.

Hackathon Dynamics

    1. (a)

Groups

       * •

How groups were formed?

       * •

Was there interdisciplinaraity?

       * •

What was you role in the group?

    2. (b)

Creative process

       * •

Did the organizers / mentors suggest the use of any creative technique for
idea generation?

       * •

Did you know / use any specific Design Thinking methods?

       * •

What was the step-by-step process that your team followed to design and
develop applications in hackathons?

       * •

Did you work with someone from design domain? How did he/she influence the
project?

       * •

How were the tasks managed?

       * •

Did you already have any idea for the hackathon project? If yes, how was it
generated?

       * •

Could you identify or have access to the target audience or domain experts?

       * •

Were there any guidelines to generate and select the application solution
alternatives?

       * •

Did you research about potential competitors of your solution?

       * •

Did you think of a differentiation strategy in relation to competitors?

       * •

Did you do any low fidelity (lo-fi) prototype?

       * •

Could you do any validation of the prototype or application with users or
experts?

       * •

How tests were handled?

       * •

How the pitch was made?