time-reversal and a resulting $\pm$ symmetry in the many-body energies), with
each choice of normal modes corresponding to a particular QPV [18]. We may
choose one particular set and, under OBCs, diagonalize the BKC Hamiltonian as
$H_{\text{BKC}}=\sum_{\mu}\omega_{\mu}(\psi_{\mu}^{\dagger}\psi_{\mu}+1/2)$,
for bosonic quasi-particles $\\{\psi_{\mu}\\}$ that satisfy the CCRs [32, 18]
and are ‘standing-waves’ in the squeezed basis $\alpha_{j}(r)$,

$\psi_{\mu}=\sqrt{\frac{2}{N+1}}\,\sum_{j=1}^{N}i^{j}\sin\left(\frac{\mu\pi
j}{N+1}\right)\alpha_{j}(r).$

From here, the dissipator is given by

$\mathcal{D}=2\kappa\sum_{\mu=1}^{N}\mathcal{D}[\psi_{\mu}].$ (52)

That is, in the normal mode basis, the chain is a set of $N$ decoupled
dissipative harmonic oscillators with mode-varying frequencies $\omega_{\mu}$
and uniform damping $\kappa$ 888The fact that the dissipator has more than one
diagonal representation, Eqs. (51) and (52), results from a large degeneracy
in the GKLS matrix spectrum. In particular, $\mathbf{M}$ has an $N$-fold
degenerate $0$ eigenvalue. Different diagonal representation of $\mathcal{D}$
are derived from different choices of bases for the $0$ eigenvalue eigenspace.
. Equivalently, $L_{\mu}=\sqrt{2\kappa}\,\psi_{\mu}$. We further note that
dissipators diagonal in the normal mode basis of the system Hamiltonian often
arise in self-consistent derivations of quadratic master equations [82].

####  V.4.1 Relaxation dynamics

Eq. (52) makes the identification of the SS trivial. If $\ket{\widetilde{0}}$
is the QPV of the BKC, that is, $\psi_{\mu}\ket{\widetilde{0}}=0$ for all
$\mu$, then $\mathcal{L}(\ket{\widetilde{0}}\bra{\widetilde{0}})=0$. This
state is a squeezed Gaussian state that we have constructed explicitly in Ref.
18. Since it is defined as a QPV of $H_{\text{BKC}}$, it is necessarily a
function of the Hamiltonian parameters $J$ and $\Delta$ only. In particular,
it is insensitive to $\kappa$, and hence to the topology of the rapidity
bands. This demonstrates that topological metastability need not be reflected
in the structure of the SS.

Given a system with both a pure SS and topological metastability, the manifold
of quasi-SSs admits a simple characterization. Let $\gamma^{s}$ be one of two
edge-symmetries $\gamma_{L}^{s}$ or $\gamma_{R}^{s}$. Following Sec. V.1, we
construct the Weyl-displaced quasi-SSs, $\rho_{\theta}\equiv
e^{i\theta\gamma^{s}}\ket{\widetilde{0}}\bra{\widetilde{0}}e^{-i\theta\gamma^{s}}$,
which are also pure. In fact, the states
$\ket{\vec{\alpha}(\theta),\\{\psi_{\mu}\\}}\equiv
e^{i\theta\gamma^{s}}\ket{\widetilde{0}}$, which are generically squeezed
coherent states with respect to the physical degrees of freedom $\\{a_{j}\\}$,
can be interpreted as coherent states with respect to the squeezed bosonic
normal modes $\\{\psi_{\mu}\\}$. Explicitly,

$\psi_{\mu}\ket{\vec{\alpha}(\theta),\\{\psi_{\mu}\\}}=\alpha_{\mu}(\theta)\ket{\vec{\alpha}(\theta),\\{\psi_{\mu}\\}},\quad\alpha_{\mu}(\theta)\in{\mathbb{C}}.$

The amplitudes $\alpha_{\mu}(\theta)$ encode the extent to which the normal
modes commute with the SG, i.e.,
$[\psi_{\mu},\gamma^{s}]=i\theta\vec{\psi}_{\mu}^{\dagger}\bm{\tau}_{3}\vec{\gamma}^{s}\,1_{\mathcal{F}}$
$=\alpha_{\mu}(\theta)1_{\mathcal{F}}$, where $\vec{\psi}_{\mu}$ and
$\vec{\gamma}^{s}$ are the single-particle vectors associated to $\psi_{\mu}$
and $\gamma^{s}$, respectively. As in the general case, the mean vector
$\vec{m}(t)=\braket{\vec{\alpha}(\theta),\\{\psi_{\mu}\\}}{\Phi(t)}{\vec{\alpha}(\theta),\\{\psi_{\mu}\\}}=i\theta\vec{\gamma}^{s}(t)$
is long-lived (recall Eq. (46)). Due to the correspondence between dynamical
matrices, this is the same mean vector as the mixed quasi-SSs in the DBKC.

Identifying the quasi-SSs as coherent states evolving under decoupled
dissipative harmonic dynamics (in the normal-mode basis) affords us the
ability to compute their dynamics exactly [83, 84] and analyze their
relaxation dynamics more precisely. As they evolve, the states remain
coherent, i.e.,
$\rho_{\theta}(t)=\ket{\vec{\alpha}(\theta,t),\\{\psi_{\mu}\\}}\bra{\vec{\alpha}(\theta,t),\\{\psi_{\mu}\\}}$,
with amplitudes

$\alpha_{\mu}(\theta,t)=\alpha_{\mu}(\theta)e^{-(\kappa+i\omega_{\mu})t}=\braket{\psi_{\mu}}(t),$

which relax to equilibrium exponentially fast in time. This apparent paradox
is resolved by noting that these amplitudes correspond to the expectation
values of $\psi_{\mu}$ which, due to squeezing, have exponentially large
coefficients when expressed in the basis of the physically relevant degrees of
freedom $\\{a_{j}\\}$. This explains why the normal-mode amplitudes decay
exponentially, despite the mean vector being long-lived, with a lifetime
increasing in system size. While this property holds independently of which SG
is chosen (left- or right-localized), the exact nature of the relaxation
dynamics do.

Figure 5: (a) Relaxation dynamics of the left-localized quasi-SS as $N$
increases. (b) Same as in (a), but for the right-localized quasi-SS. (c)
Relaxation times $t_{\text{rel}}(\delta)$ for the left- (red) and right-
(blue) localized quasi-SSs for accuracies $\delta=0.75$ (solid) and
$\delta=0.90$ (dashed). (d) Scaling of the initial amplitude vector norms for
the left- (red) and right- (blue) localized quasi-SSs with system size. In all
cases, we take $\theta=1$, symmetrically normalize the MBs, and set $J=2$,
$\Delta=0.5$, and $\kappa=0.3$.

Consider the trace distance between two quantum states,
$T(\rho,\sigma)=(1/2)\lVert\rho-\sigma\rVert_{\text{tr}}$, which provides a
measure of their distinguishability. Thanks to our knowledge of the exact
dynamics, we can study the relaxation time of the quasi-SSs via the distance
from equilibrium,

$\displaystyle T(\rho_{\theta}(t),\rho_{\text{ss}})$ $\displaystyle=$
$\displaystyle\sqrt{1-|\braket{\vec{\alpha}(\theta,t),\\{\psi_{\mu}\\}}{\widetilde{0}}|^{2}}$
$\displaystyle=$
$\displaystyle\sqrt{1-\exp\left(-\lVert\vec{\alpha}(\theta,0)\rVert^{2}e^{-2\kappa
t}\right)},$

where $\vec{\alpha}(\theta,t)$ is either $\vec{\alpha}_{L}(\theta,0)$ or
$\vec{\alpha}_{R}(\theta,0)$ depending which SG is used to generate the quasi-
SSs. We then define the relaxation time to be the time
$t_{\text{rel}}(\delta)$ such that the relative distance from equilibrium
$T(\rho_{\theta}(t),\rho_{\text{ss}})/T(\rho(0),\rho_{\text{ss}})$ falls (and
remains) below a prescribed accuracy $\delta>0$. It follows that

$2\kappa\,t_{\text{rel}}(\delta)=\ln\left[\frac{\lVert\vec{\alpha}(\theta,0)\rVert^{2}}{\ln\left[\left(1-\delta(1-e^{-\lVert\vec{\alpha(\theta,0)}\rVert^{2}})\right)^{-1}\right]}\right].$

The system-size scaling of $t_{\text{rel}}(\delta)$ is thus explicitly tied to
the system-size scaling of the norm of the initial amplitude vector,
$\lVert\vec{\alpha}(\theta,0)\rVert$. Remarkably, the two manifolds of quasi-
SSs display dramatically different relaxation dynamics, as inferred from their
relaxation times. We find, numerically, that
$\lVert\vec{\alpha}_{L}(\theta,0)\rVert$ increases exponentially with $N$,
whereas $\lVert\vec{\alpha}_{R}(\theta,0)\rVert$ decreases exponentially. The
consequences of this observation are displayed in Fig. 5: the left-localized
quasi-SS exhibits a increasingly long relaxation time, while the right-
localized shows just the opposite. The asymmetry in the dynamics may be
understood by noting that, while the dynamical matrix is translationally
invariant (up to boundaries), the full generator is not. Specifically, the
Lindblad operators $\alpha_{j}(r)$ are right-localized.

The boundedness of the relaxation times for the right-localized state does not
contradict our claim of long-livedness. Per the general theory, the lifetime
of the physically accessible degrees of freedom, as captured by the
macroscopic mean vector $\vec{m}_{\theta}^{R}(t)$, remains indistinguishable
from its initial value $\vec{m}_{\theta}^{R}(0)$, as exemplified by Eq. (46).
On the other hand, the quasi-SS and the true SS have overlap exponentially
close to one as $N$ increases. The resolution comes by noting that
$\vec{m}^{R}_{\theta}(t)$ is obtained from the exponentially small normal mode
amplitudes $\vec{\alpha}_{R}(\theta,t)$ via a squeezing transformation. This
transformation dilates the exponentially small amplitudes into the physical
amplitudes. Similarly, the same transformation squeezes the exponentially
large normal mode amplitudes $\vec{\alpha}_{L}(\theta,t)$ into the physical
amplitudes $\vec{m}_{\theta}^{L}(\theta,t)$.

####  V.4.2 Transient odd-parity behavior

To conclude our analysis of this model, we pose the question: Do the pure
quasi-SSs span a subspace of long-lived states? That is, do linear
combinations of these pure states provide long-lived states themselves? This
question is interesting from two perspectives. Firstly, from the perspective
of continuous-variable quantum information, linear combinations of coherent
states provide the simplest realization of bosonic cat-codes [85, 86]. Thus,
the ability to potentially dissipatively prepare and sustain such states could
have unforeseen implications, for instance in the context of robust encodings.
Secondly, from the perspectives of condensed-matter physics and quantum
optics, ground states and SSs are intrinsically averse to odd bosonic parity.
Specifically, if a QBH (QBL) has a unique ground (steady) state, then it must
have even bosonic parity. In the Hamiltonian case, this eliminates the
possibility for the parity switching behavior characteristic of topological
superconductors. Moreover, negative expectation values of parity signify
quantum non-Gaussianity.

In general, the question is difficult to answer. However, we can make
interesting statements about particularly relevant linear combinations.
Concretely, let us consider bosonic “cat states” of the form

$\ket{\mathcal{C}_{\phi}(\vec{\alpha})}=\mathcal{N}_{\phi}(\vec{\alpha})\Big{(}\ket{\vec{\alpha}(\theta),\\{\psi_{\mu}\\}}+e^{i\phi}\ket{-\vec{\alpha}(\theta),\\{\psi_{\mu}\\}}\Big{)},$

where $\mathcal{N}_{\phi}(\vec{\alpha})$ is a normalization constant, and
$\phi\in[0,2\pi]$. Being generally non-Gaussian, these states are not simply
characterized by their first and second cumulants. Moreover, computation of
overlaps becomes substantially more difficult than in the previous case.
Instead of a direct analysis of relaxation times, we focus explicitly on
parity dynamics,

$\displaystyle\braket{P}_{\phi,\vec{\alpha}}(t)=\braket{\mathcal{C}_{\phi}(\vec{\alpha})}{P(t)}{\mathcal{C}_{\phi}(\vec{\alpha})},$

with $P$ the bosonic parity operator introduce in Sec. II.2.2. It is then
possible to derive the following exact analytical result (see App. C for
detail):

$\braket{P}_{\phi,\vec{\alpha}}(t)=\frac{e^{-2\lVert\vec{\alpha}\rVert^{2}}+\cos(\phi)e^{-2\lVert\vec{\alpha}\rVert^{2}(1-e^{-2\kappa
t})}}{1+\cos(\phi)e^{-2\lVert\vec{\alpha}\rVert^{2}}}.$ (53)

Again, which quasi-SS manifold is chosen strongly affects the ensuing
behavior. If $\vec{\alpha}=\vec{\alpha}_{L}(\theta,0)$, the parity shows a
dramatic dependence on both $N$ and $\phi$, as seen in Fig. 6. Unless
$\phi=\pm\pi/2$, the parity drops incredibly fast to zero, remains zero for a
transient time that scales linearly with $N$, and then eventually rises up to
the asymptotic value of $1$. The extremely fast initial drop to zero
corresponds to a singularity in the derivative

$\braket{\dot{P}}_{\phi,\vec{\alpha}}(0)=-4\lVert\vec{\alpha}_{L}(\theta,0)\rVert^{2}\bigg{(}\frac{e^{-2\lVert\vec{\alpha}\rVert^{2}}-\cos(\phi)}{1+e^{-2\lVert\vec{\alpha}\rVert^{2}}\cos(\phi)}\bigg{)},$

as $N$ (and hence $\lVert\vec{\alpha}\rVert$) goes to infinity. The transient
state of zero parity that follows is interesting, as it represents a long-
lived period where the measurement statistics of parity are split evenly
between the $+1$ and $-1$ outcomes. While we do not sustain a state of odd
parity explicitly, we do, in fact, sustain an even mixture between even and
odd parity sectors, for a transient that grows linearly with system size.

In sharp contrast, if $\vec{\alpha}=\vec{\alpha}_{R}(\theta,0)$, the parity
dynamics are much more well-behaved. The exponentially small norms of
$\vec{\alpha}_{R}(\theta,0)$ ensure that, for sufficiently large, $N$,

$\lim_{N\to\infty}\braket{P}_{\phi,\vec{\alpha}_{R}}(t)=\begin{cases}1&\phi\neq\pi,\\\
1-2e^{-2\kappa t}&\phi=\pi.\end{cases}$

Accordingly, unless $\phi=\pi$, the parity approaches $1$ for all $t$ as $N$
increases. Moreover, when $\phi=\pi$, the parity is indistinguishable from
$1-e^{-2\kappa t}$ for sufficiently large $N$. Regardless, unlike their left-
localized partners, these right-localized cat states fail to support any
semblance of an odd-parity state for any meaningful amount of time.

Figure 6: (a) Parity dynamics of the initial odd-parity cat state formed from
the left-localized quasi-SS for increasing $N$. Inset: Short time behavior.
(b) Parity dynamics of the initial cat state formed from the left-localized
quasi-SS for $N=25$ and varying $\phi$. (c) Parity dynamics of the initial
odd-parity cat state formed from the right-localized quasi-SS for increasing
$N$. The dashed line indicates the $N\to\infty$ limit. (d) Parity dynamics of
the initial cat state formed from the right-localized quasi-SS for $N=25$ and
varying $\phi$. In all cases, we take $\theta=1$, symmetrically normalize the
MBs, and set $J=2$, $\Delta=0.5$, and $\kappa=0.3$.

##  VI Topological dynamical metastability in number-symmetric QBLs

###  VI.1 Dirac edge bosons: General properties

The topological modes in all the model systems of Sec. V bear the closest
resemblance to the Majorana fermions of topological superconductivity. Thus,
one could suspect that our analysis may fail to produce topological modes
analogous to the ones hosted by topological insulators. To see why this is
false, we consider QBLs with number symmetry to provide the natural
dissipative-bosonic analogues of insulators.

Number-symmetric QBLs possess a weak U(1) symmetry of the form $a_{j}\mapsto
e^{i\phi}a_{j}$, generated by the total number operator
$\sum_{j}a_{j}^{\dagger}a_{j}$. The corresponding Bogoliubov transformation is
$\Phi\mapsto e^{i\bm{\tau}_{3}\theta}\Phi$. Thus, number symmetry implies
$[\mathbf{G},\bm{\tau}_{3}]=0$, and the dynamical matrix may be written as

$\mathbf{G}=\mathbf{K}\otimes\begin{bmatrix}1&0\\\
0&0\end{bmatrix}+(-\mathbf{K}^{*})\otimes\begin{bmatrix}0&0\\\
0&1\end{bmatrix},$ (54)

with $\mathbf{K}$ an arbitrary $N\times N$ complex matrix. Just as the spectra
of $\mathbf{G}$ is determined by that of $\mathbf{K}$, so is the
pseudospectra. In particular, if the system is topologically metastable, both
$\mathbf{G}$ and $\mathbf{K}$ possess approximate kernel vectors.

To be more concrete, let us suppose the system is topologically metastable.
This affords us at least one MB pair $(\gamma^{z}_{1},\gamma^{s}_{1})$.
Expanding in the bosonic basis yields

$\displaystyle\gamma^{z}_{1}=\sum_{j=1}^{N}(u_{j}^{*}a_{j}+u_{j}a^{\dagger}_{j}),\quad\gamma^{s}_{1}=\sum_{j=1}^{N}(v_{j}^{*}a_{j}+v_{j}a^{\dagger}_{j}),$

where $u_{j}$ and $v_{j}$ are related to the Nambu vectors via
$\vec{\gamma}^{z}_{1}=[u_{1},-u_{1}^{*},\ldots,u_{N},-u_{N}^{*}]^{T},$
$\vec{\gamma}^{s}_{1}=[v_{1},-v_{1}^{*},\ldots,v_{N},-v_{N}^{*}]^{T}.$ The
HWRs imply that $\text{Im}\sum_{j=1}^{N}u_{j}^{*}v_{j}=1/2$, since

$\displaystyle
i=[\gamma_{1}^{z},\gamma_{1}^{s}]=\sum_{j=1}^{N}(u_{j}^{*}v_{j}-v_{j}^{*}u_{j})=2i\,\text{Im}\sum_{j=1}^{N}u_{j}^{*}v_{j}.$

Number symmetry additionally implies that an approximate (or exact) ZM or SG
remains so after a phase rotation $a_{j}\mapsto e^{i\phi}a_{j}$. In
particular, we may rotate each of our MBs to construct a second, linearly
independent MB pair. Specifically, if we fix $\phi=\pi/2$, then we immediately
find another linearly independent MB pair given by

$\displaystyle\gamma^{z}_{2}=i\sum_{j=1}^{N}(u_{j}^{*}a_{j}-u_{j}a^{\dagger}_{j}),\quad\gamma^{s}_{2}=i\sum_{j=1}^{N}(v_{j}^{*}a_{j}-v_{j}a^{\dagger}_{j}).$

We have not yet encountered anything resembling the edge modes of topological
insulators. For that, we need to utilize both MB pairs. Explicitly, consider
the operators

$\displaystyle\alpha$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2\sqrt{C^{z}}}(\gamma^{z}_{1}-i\gamma^{z}_{2})=\frac{1}{\sqrt{C^{z}}}\sum_{j=1}^{N}u_{j}^{*}a_{j},$
$\displaystyle\beta$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2\sqrt{C^{s}}}(\gamma^{s}_{1}-i\gamma^{s}_{2})=\frac{1}{\sqrt{C^{s}}}\sum_{j=1}^{N}v_{j}^{*}a_{j},$

where $C^{z}=\sum_{j}|u_{j}|^{2}$ and $C^{s}=\sum_{j}|v_{j}|^{2}$ are positive
real numbers chosen to ensure that the Dirac boson modes $\alpha$ and $\beta$
satisfy the following properties:

  * •

They are bosonic:
$[\alpha,\alpha^{\dagger}]=[\beta,\beta^{\dagger}]=1_{\mathcal{F}}$.

  * •

They are algebraically related via $[\alpha,\beta]=0$ and
$[\alpha,\beta^{\dagger}]\neq 0$.

  * •

They are edge-localized according to the localization of the constituent MBs.

  * •

The bosonic mode $\alpha$ is approximately conserved, while the real and
imaginary quadratures of $\beta$ generate two (non-commuting) approximate
symmetries.

Algebraically speaking, we have constructed edge modes on which the number
symmetry acts trivially. The necessary and sufficient ingredients for this
construction are topological metastability and number symmetry, and so it
applies broadly.

The normalization of these operators is far more straightforward than in the
number non-symmetric case. To see this, consider a left-localized,
approximately conserved, bosonic operator $\alpha$, i.e.,

$\alpha=\mathcal{M}(N)\sum_{j=1}^{N}\delta^{j-1}a_{j},\quad|\delta|<1.$ (55)

Taking $\delta>0$ without loss of generality, canonical commutation relations
require

$\displaystyle\mathcal{M}(N)=\sqrt{\frac{1-\delta^{2}}{1-\delta^{2N}}}$ (56)

which, unlike in the case of the MBs, converges to a finite value as
$N\to\infty$. Namely, $\lim_{N\to\infty}\mathcal{M}(N)=\sqrt{1-\delta^{2}}$.
With this, the exact ZM in the $N\to\infty$ limit is

$\alpha=\sqrt{1-\delta^{2}}\,\sum_{j=1}^{\infty}\delta^{j-1}a_{j}.$

While, in general, the exact expression for $\alpha$ (and, similarly, $\beta$)
could be more complicated than Eq. (55), this argument demonstrates the
unambiguous normalizability of the modes in the infinite size limit.

It is interesting to ask whether a number-symmetric analogue to the PDMC may
exist, namely, a purely dissipative, number-symmetric chain that exhibits
dynamical metastability. In the purely dissipative case, we have
$\mathbf{G}=-i\bm{\tau}_{3}\mathcal{F}(\mathbf{M})$, and the number symmetry
property $[\mathbf{G},\bm{\tau}_{3}]=0$ manifests as
$[\mathcal{F}(\mathbf{M}),\bm{\tau}_{3}]=0$. Thus,
$\mathbf{G}^{\dagger}=(-i\bm{\tau}_{3}\mathcal{F}(\mathbf{M}))^{\dagger}=i\bm{\tau}_{3}\mathcal{F}(\mathbf{M})=-\mathbf{G}.$
That is, a purely dissipative number-symmetric chain must have an anti-
Hermitian (therefore, normal) dynamical matrix $\mathbf{G}$. We conclude that
no such model can exhibit dynamical metastability.

###  VI.2 A number-symmetric dissipative chain

Let us explore the interplay between number symmetry and topological
metastability in a concrete example. The dynamical matrix of a number-
symmetric QBL has the general form Eq. (54), where $\mathbf{K}$ an arbitrary
$N\times N$ complex matrix. Topological metastability can then be engineered
through an appropriate choice of $\mathbf{K}$. Specifically, let us consider

$\mathbf{K}=-i\kappa\mathds{1}_{N}+J_{L}\mathbf{S}+J_{R}\mathbf{S}^{\dagger},$
(57)

with $\kappa\geq 0$, $J_{L},J_{R}\in\mathbb{R}$, and $\mathbf{S}$ the usual
BC-dependent shift operator. We identify this matrix as that of the Hatano-
Nelson (HN) asymmetric hopping model [87], with an identity shift that will
ultimately serve to stabilize the QBL. For convenience, we define
$J_{\pm}\equiv(J_{L}\pm J_{R})/2$. The Hamiltonian can be unambiguously
determined from $\mathbf{G}$ and is given by

$H_{\text{NS}}=\frac{J_{+}}{2}\sum_{j=1}^{N}\left(a_{j}^{\dagger}a_{j+1}+a_{j+1}^{\dagger}a_{j}\right).$

Per usual, the QBL is not fully determined until we specify
$\mathcal{B}(\mathbf{M})$. Moreover, the second necessary condition for the
U(1) symmetry is that $\bm{\tau}_{3}\mathcal{B}(\mathbf{M})$ commutes with
$\bm{\tau}_{3}$. Together, $[\mathbf{G},\bm{\tau}_{3}]=0$ and
$[\bm{\tau}_{3}\mathcal{B}(\mathbf{M}),\bm{\tau}_{3}]=0$ are necessary and
sufficient for the requisite U(1) symmetry. We specify
$\mathcal{B}(\mathbf{M})$ implicitly, by defining the dissipator
$\mathcal{D}_{\text{NS}}\equiv\mathcal{D}_{-,0}+{\color[rgb]{0,0,0}\mathcal{D}_{+,0}+}\mathcal{\mathcal{missing}}{D}_{-,1}$,
with

$\displaystyle\mathcal{D}_{-,0}$ $\displaystyle=$ $\displaystyle
2\kappa_{-}\sum_{j=1}^{N}\mathcal{D}[a_{j}],$ $\displaystyle\mathcal{D}_{+,0}$
$\displaystyle=$ $\displaystyle
2\kappa_{+}\sum_{j=1}^{N}\mathcal{D}[a_{j}^{\dagger}],$
$\displaystyle\mathcal{D}_{-,1}$ $\displaystyle=$ $\displaystyle
2iJ_{-}\sum_{j=1}^{N}\mathcal{D}[a_{j},a_{j+1}^{\dagger}]-\mathcal{D}[a_{j+1},a_{j}^{\dagger}].$

Here, $2\kappa_{-}\geq 0$ and $2\kappa_{+}\geq 0$ are the onsite loss and gain
rates, respectively, while $2J_{-}$ takes the role of the NN loss rate. The
GKLS matrix is positive-semidefinite for OBCs and PBCs, and for all $N$, if
$\kappa_{-}\geq 2|J_{-}|$. This QBL has a dynamical matrix specified by Eq.
(57) if we further identify $\kappa\equiv\kappa_{-}-\kappa_{+}$. We refer to
this model as the dissipative number-symmetric (DNS) chain. The SS behavior of
a related model has been considered in Ref. 88.

The rapidities can be easily obtained from the well-known HN spectrum, and
closely resemble that of the DBKC. For BIBCs, the bands are given by
$\\{\lambda(k),\lambda(k)^{*}\\}$, with

$\lambda(k)=-\kappa+2J_{-}\sin(k)+i2J_{+}\cos(k).$

These bands trace out an ellipse centered at $-\kappa$ in the complex plane.
Winding about the origin requires $2|J_{-}|>\kappa$. For (finite-size) OBCs,
the eigenvalues are given by $\lambda_{m}$, $m=0,\ldots,2N-1$, where

$\lambda_{m}=-\kappa+2i\sqrt{J_{+}^{2}-J_{-}^{2}}\cos\left(\frac{m\pi}{N+1}\right).$

To simplify the discussion (and fix the OBC Lindblad gap
$\Delta_{\mathcal{L}}=\kappa$, for all $N$), we focus on $|J_{+}|\geq|J_{-}|$.
Combining the GKLS-matrix positivity condition and the rapidity band-winding
condition, we identify a topologically metastable regime whenever
$\kappa_{-}/|J_{-}|\geq 2$ and $0\leq\kappa/|J_{-}|\leq 2$. The OBC stability
phase diagram is shown in Fig. 7.

Figure 7: The stability phase diagram for the DNS under OBCs with
$|J_{+}|>|J_{-}|$ under OBCs. The “Ill-defined” region corresponds to the
parameter regime where $\mathbf{M}$ is no-longer positive semidefinite.

In the special case $J_{R}=0$, the pseudoeigenvectors of $\mathbf{K}$ (and
$\mathbf{K}^{\dagger}$) with zero pseudoeigenvalue may be computed
analytically. Such pseudoeigenvectors can be used to build approximate kernel
vectors of $\mathbf{G}$ and $\widetilde{\mathbf{G}}=\mathbf{G}^{\dagger}$,
which correspond to approximate SGs and ZMs, respectively. Specifically,
consider the bosonic modes

$\displaystyle\alpha\equiv\mathcal{M}(N)\sum_{j=1}^{N}(i\delta)^{N-j}a_{j},\quad\beta\equiv\mathcal{M}(N)\sum_{j=1}^{N}(-i\delta)^{j-1}a_{j},$

where $\delta\equiv-\kappa/2J_{-}$ and $\mathcal{M}(N)$ is given in Eq. (56).
First, we have
$[\alpha,\alpha^{\dagger}]=[\beta,\beta^{\dagger}]=1_{\mathcal{F}}$ and
$[\alpha,\beta]=0$. Second, these operators satisfy

$\displaystyle\mathcal{L}^{\star}(\alpha)$
$\displaystyle=-\kappa\mathcal{M}(N)\,(-i\delta)^{N-1}a_{1},$ (58)
$\displaystyle\mathcal{L}^{\star}([\beta,A])-[\beta,\mathcal{L}^{\star}(A)]$
$\displaystyle=-\kappa\mathcal{M}(N)\,(i\delta)^{N-1}[a_{N},A],\;\forall A.$
(59)

Utilizing
$\mathcal{L}^{\star}(A^{\dagger})=[\mathcal{L}^{\star}(A)]^{\dagger}$ yields
similar expressions for $\alpha^{\dagger}$ and $\beta^{\dagger}$. For
$|\delta|<1$, the right hand-side of each equation goes to zero as
$N\to\infty$. Physically, this means $\alpha$ is a bosonic approximate ZM
while the real and imaginary quadratures of $\beta$ generate approximate
symmetries. That is, $\alpha$ and $\beta$ are DBs.

We may trace these two bosonic modes back to two pairs of MBs using Eqs.
(VI.1). Specifically, we have

$\displaystyle\gamma_{1}^{z}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{z}(N)\sum_{j=1}^{N}\left((i\delta)^{N-j}a_{j}+(-i\delta)^{N-j}a_{j}^{\dagger}\right),$
$\displaystyle\gamma_{1}^{s}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{s}(N)\sum_{j=1}^{N}\left((-i\delta)^{j-1}a_{j}+(i\delta)^{j-1}a_{j}^{\dagger}\right),$
$\displaystyle\gamma_{2}^{z}$ $\displaystyle=$ $\displaystyle
i\mathcal{M}_{z}(N)\sum_{j=1}^{N}\left((i\delta)^{N-j}a_{j}-(-i\delta)^{N-j}a_{j}^{\dagger}\right),$
$\displaystyle\gamma_{2}^{s}$ $\displaystyle=$ $\displaystyle
i\mathcal{M}_{s}(N)\sum_{j=1}^{N}\left((-i\delta)^{j-1}a_{j}-(i\delta)^{j-1}a_{j}^{\dagger}\right),$

for normalization constants $\mathcal{M}_{z}(N)$ and $\mathcal{M}_{s}(N)$
chosen to ensure that $[\gamma_{j}^{z},\gamma_{j}^{s}]=i1_{\mathcal{F}}$ for
$j=1,2$. The relevant equations of motion follow from Eqs. (58)-(59), that is,

$\displaystyle\mathcal{L}^{\star}(\gamma_{1}^{z})$
$\displaystyle=-\kappa\mathcal{M}_{z}(N)\delta^{N-1}Q_{L},$
$\displaystyle\mathcal{L}^{\star}(\gamma_{2}^{z})$
$\displaystyle=\kappa\mathcal{M}_{z}(N)\delta^{N-1}P_{L},$
$\displaystyle\mathcal{L}^{\star}([\gamma_{1}^{s},A])-[\gamma_{1}^{s},\mathcal{L}^{\star}(A)]$
$\displaystyle=-\kappa\mathcal{M}_{s}(N)\delta^{N-1}[Q_{R},A],$
$\displaystyle\mathcal{L}^{\star}([\gamma_{2}^{s},A])-[\gamma_{2}^{s},\mathcal{L}^{\star}(A)]$
$\displaystyle=\kappa\mathcal{M}_{s}(N)\delta^{N-1}[P_{R},A],$

with $Q_{L}\equiv\chi_{N}^{*}a_{1}+\chi_{N}a_{j}^{\dagger}$,
$P_{L}\equiv-i(\chi_{N}^{*}a_{1}-\chi_{N}a_{1}^{\dagger})$,
$Q_{R}\equiv\chi_{N}a_{N}+\chi_{N}^{*}a_{N}^{\dagger}$,
$P_{R}\equiv-i(\chi_{N}a_{N}-\chi_{N}^{*}a_{N}^{\dagger})$, and
$\chi_{N}\equiv i^{N-1}$. Following the discussion of Sec. V.1, there exists
normalization schemes ensuring that the right hand-sides of these four
equations vanish as $N\to\infty$, while keeping each pair canonically
conjugate (for instance, the symmetric normalization scheme we have previously
employed).

We remark that these pairs of MBs are not simply the real and imaginary
quadratures of $\alpha$ and $\beta$, i.e., the operators
$x_{\alpha}\equiv(\alpha+\alpha^{\dagger})/\sqrt{2}$, $p_{\alpha}\equiv
i(\alpha-\alpha^{\dagger})/\sqrt{2}$ (and similarly for $\beta$). Instead, the
MBs are proportional to these quadratures, according to the expressions above.
The system-size-dependent proportionality constants are key to ensuring that
the macroscopically separated pairs $(\gamma_{1}^{z},\gamma_{1}^{s})$ and
$(\gamma_{2}^{z},\gamma_{2}^{s})$ are canonically conjugate. We further remark
that each edge supports two Noether modes of the same type: the left edge
supports two SGs, whereas the right edge supports two ZMs. This is to be
contrasted with the previous models, which all featured at most one of each
type on each edge. A similar phenomena is uncovered in a model without the
additional constraint of number symmetry in App. D.

##  VII Observable signatures of topological metastability

Just like spatial, equal-time correlation functions of physical observables
play a central role in characterizing equilibrium phases of matter and phase
transitions, multi-time correlation functions provide an essential tool for
probing dynamical response properties of many-body systems away from
equilibrium [89, 90]. In particular, standard material characterization
techniques ranging from neutron scattering to reflectivity and photoemission
spectroscopy rely on access to two-time correlation functions of the form
$\langle A(t)B(t^{\prime})\rangle$ [91, 92]; likewise, in the context of
quantum optics, statistical properties of the photon field are directly
related to first- and second-order coherence functions,
$g_{ij}^{(1)}(t,t^{\prime})\propto\braket{a_{i}^{\dagger}(t)a_{j}(t^{\prime})}$,
and
$g_{ij}^{(2)}(t,t^{\prime})\propto\braket{a_{i}^{\dagger}(t)a^{\dagger}_{j}(t^{\prime})a_{i}(t^{\prime})a_{j}(t)}$,
which can be inferred through photon counting and quantum interference
experiments [58, 93].

Two-time correlation functions have been suggested as a means for
operationally characterizing multi-step relaxation dynamics in Markovian
quantum systems [33, 34]. Since, for the QBL systems we are interested in,
dynamical metastability is inherently a transient phenomena, it is natural to
expect that they may also serve as an experimentally accessible diagnostic
tool, in principle. As we will show next, it is indeed possible to identify
response functions whose behavior can not only distinguish topologically
trivial from non-trivial dynamical metastability but, in the non-trivial case,
further provide distinctive signatures of the different types of topological
bosonic edge modes the system may host.

###  VII.1 Two-time correlation functions and power spectra

Since systems in equilibrium are time-translation invariant, any two-time
response function can only depend upon the relative time, say, $\tau\equiv
t^{\prime}-t$. In non-equilibrium situations, in contrast, the full dependence
on both $t$ and $\tau$ becomes important in general. For Markovian dynamics,
two-time averages may be computed under the assumption that the quantum
regression theorem (QRT) holds [58]. Specifically, given two operators $A$ and
$B$ and a state $\rho$, their two-point correlation function may then be
expressed as

$\displaystyle C^{+}_{A,B}(t,\tau)$ $\displaystyle\equiv$
$\displaystyle\braket{A(t+\tau)B(t)}{\color[rgb]{0,0,0}\overset{\text{QRT}}{=}\text{tr}{\\{A(0)e^{\tau{\mathcal{L}}}[B(0)\rho(t)]\\}}}$
$\displaystyle=$ $\displaystyle\text{tr}[A(\tau)B(0)\rho(t)],\quad t,\tau\geq
0,$

where $A(\tau)=e^{\tau\mathcal{L}^{\star}}(A)$ and
$\rho(t)=e^{t\mathcal{L}}(\rho)$ and the superscript $+$ signifies that $B$ is
measured first. In this way, under the same assumptions, we also have

$\displaystyle C^{-}_{A,B}(t,\tau)$ $\displaystyle\equiv$
$\displaystyle\braket{A(t)B(t+\tau)}\overset{\text{QRT}}{=}\text{tr}[A(0)B(\tau)\rho(t)]$
$\displaystyle=$ $\displaystyle[C^{+}_{B,A}(t,\tau)]^{*},\quad t,\tau\geq 0.$

In the case of a unique SS $\rho_{\text{ss}}$, we consider the SS two-time
correlation functions $\lim_{t\to\infty}C^{\pm}_{A,B}(t,\tau)$, which we can
express compactly as

$C^{\text{ss}}_{A,B}(\tau)=\begin{cases}\text{tr}[A(\tau)B(0)\rho_{\text{ss}}],&\tau\geq
0,\\\ \text{tr}[A(0)B(|\tau|)\rho_{\text{ss}}],&\tau<0.\end{cases}$ (60)

As it may be more practically significant in certain situations, we
additionally define the (two-sided) SS power spectrum

$S^{\text{ss}}_{A,B}(\omega)=\int_{-\infty}^{\infty}e^{i\omega\tau}C^{\text{ss}}_{A,B}(\tau)\,d\tau.$

Long-lived correlations are then revealed through large power-spectral peaks
at zero frequency. To more appropriately captures the relative decay of
correlations, we define the normalized correlation functions and power spectra
[94],

$\widetilde{C}^{\text{ss}}_{A,B}(\tau)\equiv\frac{C^{\text{ss}}_{A,B}(\tau)}{C^{\text{ss}}_{A,B}(0)},\quad\widetilde{S}^{\text{ss}}_{A,B}(\omega)\equiv\frac{S^{\text{ss}}_{A,B}(\omega)}{C^{\text{ss}}_{A,B}(0)}.$
(61)

For QBLs, we will primarily focus on the case where $A$ and $B$ are linear
forms, i.e., $A=\widehat{\vec{\alpha}}$ and
$B=\widehat{\vec{\beta}}\mbox{}^{\dagger}$, with
$\vec{\alpha},\vec{\beta}\in\mathbb{C}^{2N}$. In this case, the unnormalized
SS correlation functions and power spectra take on a simple closed form:

$\displaystyle
C^{\text{ss}}_{\alpha,\beta^{\dagger}}(\tau)=\begin{cases}\vec{\alpha}{}^{\dagger}\bm{\tau}_{3}e^{-i\mathbf{G}\tau}\mathbf{Q}_{\text{ss}}\bm{\tau}_{3}\vec{\beta}&\tau\geq
0,\\\
\vec{\alpha}{}^{\dagger}\bm{\tau}_{3}\mathbf{Q}_{\text{ss}}e^{i\mathbf{G}^{\dagger}\tau}\bm{\tau}_{3}\vec{\beta}&\tau<0,\end{cases}$
$\displaystyle
S^{\text{ss}}_{\alpha,\beta^{\dagger}}(\omega)=\vec{\alpha}\bm{\tau}_{3}[\bm{\chi}(\omega)\mathbf{Q}_{\text{ss}}+\mathbf{Q}_{\text{ss}}\bm{\chi}^{\dagger}(\omega)]\bm{\tau}_{3}\vec{\beta},$

in terms of a susceptibility matrix [44, 46]

$\bm{\chi}(\omega)\equiv i(\omega\mathds{1}_{2N}-\mathbf{G})^{-1}.$

Mathematically, $\bm{\chi}(\omega)$ is resolvent of $-i\mathbf{G}$ evaluated
at $-i\omega$.

The restriction to linear forms further yields a state-independent notion of
correlation functions. Note that

$\displaystyle\alpha(\tau)\beta^{\dagger}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\\{\alpha(\tau),\beta^{\dagger}(0)\\}+\frac{1}{2}[\alpha(\tau),\beta^{\dagger}(0)]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\\{\alpha(\tau),\beta^{\dagger}(0)\\}+\frac{1}{2}\vec{\alpha}^{\dagger}e^{-i\widetilde{\mathbf{G}}^{\dagger}\tau}\bm{\tau}_{3}\vec{\beta}\,1_{\mathcal{F}}.$

Since quantum states have unit trace, we have

$C^{+}_{\alpha,\beta^{\dagger}}(t,\tau)=C^{+,\text{cl}}_{\alpha,\beta^{\dagger}}(t,\tau)+\frac{1}{2}\vec{\alpha}^{\dagger}e^{-i\widetilde{\mathbf{G}}^{\dagger}\tau}\bm{\tau}_{3}\vec{\beta},$

where
$C^{+,\text{cl}}_{\alpha,\beta^{\dagger}}(t,\tau)=\frac{1}{2}\text{tr}[\\{\alpha(\tau),\beta^{\dagger}(0)\\}\rho(t)]$
is the classical (symmetrized) correlation function. We identify the quantum
(antisymmetrized) correlation functions as

$C^{+,\text{qu}}_{\alpha,\beta^{\dagger}}(\tau)=\frac{1}{2}\vec{\alpha}^{\dagger}e^{-i\widetilde{\mathbf{G}}^{\dagger}\tau}\bm{\tau}_{3}\vec{\beta},$

whose state- and $t$-independence follows from the state independent nature of
$\braket{[\alpha(\tau),\beta^{\dagger}(0)]}$ when $\alpha$ and $\beta$ are
linear forms. We define $C^{-,\text{cl}}_{\alpha,\beta^{\dagger}}(t,\tau)$ and
$C^{-,\text{qu}}_{\alpha,\beta^{\dagger}}(\tau)$ analogously. In the case
where $\alpha$ and $\beta$ are observables (self-adjoint), we have

$i\text{Im}\,C_{\alpha,\beta^{\dagger}}^{\pm}(t,\tau)=C^{\pm,\text{qu}}_{\alpha,\beta^{\dagger}}(\tau),\quad\forall
t,$

which directly relates the quantum correlation function to the imaginary part
of the full one, independently of the state. Similar to Eq. (60), we can drop
the $\pm$ by defining

$C^{\text{qu}}_{\alpha,\beta^{\dagger}}(\tau)\equiv\begin{cases}C^{+,\text{qu}}_{\alpha,\beta^{\dagger}}(\tau)&\tau\geq
0,\\\ C^{-,\text{qu}}_{\alpha,\beta^{\dagger}}(|\tau|)&\tau<0,\end{cases}$

The corresponding frequency-space quantum power spectrum may be defined and
expressed as

$\displaystyle S_{\alpha,\beta^{\dagger}}^{\text{qu}}(\omega)$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}e^{i\omega\tau}C^{\text{qu}}_{\alpha,\beta^{\dagger}}(\tau)\,d\tau$
$\displaystyle=$
$\displaystyle\vec{\alpha}^{\dagger}\left(\bm{\tau}_{3}\bm{\chi}(\omega)+\bm{\chi}^{\dagger}(\omega)\bm{\tau}_{3}\right)\vec{\beta},$

which is, again, a state-independent quantity.

###  VII.2 Distinguishing topologically non-trivial from trivial dynamical
metastability

Dynamically metastable systems can be either topologically trivial or non-
trivial. Starting from the susceptibility matrix, one can predict certain
properties of the SS power spectrum that can distinguish the two regimes. We
will focus on the normalized power spectra [Eq. (61)] in order to capture the
relative behavior of correlations. This eliminates the influence of
exponentially large steady state second moments $\mathbf{Q}_{\text{ss}}$
(e.g., occupation numbers) that may arise in systems displaying transient
amplification.

Figure 8: (a-d) The 2-norm of the DBKC’s OBC susceptibility matrix in the
anomalously relaxing, topologically metastable, and non-topological
dynamically metastable regimes, respectively. (a) and (b) corresponds to the
regimes whose rapidities are the open and filled markers in Fig. 1(a),
respectively. (c) corresponds to the regime whose rapidities are shown in Fig.
1(b). (d) shows the zero-frequency susceptibility norm in the three
aforementioned regimes. The dashed line is a linear fit.

Let $\bm{\chi}_{N}(\omega)$ denote the susceptibility matrix/resolvent of the
dynamical matrix $-i\mathbf{G}_{N}$ for an open chain of length $N$. On the
one hand, if the chain is anomalously relaxing, then necessarily
$\bm{\chi}_{N}(\omega)$ is bounded (in norm) for all $\omega$. The reason is
that the rapidity bands of anomalously relaxing systems are strictly bound to
the left-half of the complex plane and so $-i\omega$ is not in the SIBC
spectrum. Thus, we have a system-size independent upper bound on
$\lVert{\color[rgb]{0,0,0}\bm{\chi}}_{N}(\omega)\rVert$. On the other hand, if
the chain is dynamically metastable, then there is necessarily a subset of the
imaginary axis contained within the SIBC rapidity spectrum. Equivalently,
there are intervals on the imaginary axis about which the rapidity bands wind.
Ultimately, the restriction of $\bm{\chi}_{N}(\omega)$ to these intervals will
necessarily grow without bound as $N\to\infty$. Since topological
metastability is characterized by the presence of zero in these non-trivial
intervals, we conjecture that it generically elicits a peak of the power
spectrum at zero frequency, which grows without bound with system size. There
should be no such peak in a dynamically metastable system that is
topologically trivial.

The distinctive behavior of $\bm{\chi}_{N}(\omega)$ in these regimes is
exemplified in Fig. 8. In (a), the $2$-norm of the susceptibility matrix
converges for all values of $\omega$ considered. In (b) and (c), we see a
divergence of the norm at frequencies $\omega$ such that $i\omega$ is
contained in the non-trivial interior of the rapidity bands. The zero
frequency behavior of $\lVert\bm{\chi}_{N}(\omega)\rVert_{2}$ is shown in (d).
In particular, the topological regime is distinguished from both the non-
topological dynamically metastable regime and the anomalously relaxing regime
by an exponential divergence in system size. We will further see that this
manifests directly in certain quantum power spectra.

####  VII.2.1 Signatures of Majorana edge bosons

Within the class of topological metastable QBLs with no U(1) symmetry, models
may be distinguished based on whether their MBs are split or non-split. Let
$(\gamma^{z},\gamma^{s})$ denote a split MB pair whose constituent modes
satisfy $[\gamma^{z},\gamma^{s}]=i$ and are localized on opposite sides of a
chain. During the transient timescale ($t<t_{N}$, for some $t_{N}$ increasing
with system size $N$), we have $\gamma^{z}(t)\simeq\gamma^{z}(0)$. However,
because the MBs are split, $\gamma^{s}(t)$ deviates significantly from
$\gamma^{s}(0)$ over the same timescale. Consider the associated correlation
function $C^{\text{ss}}_{\gamma^{z},\gamma^{s}}(\tau)$, obtained from Eq.
(60). Firstly, we note that canonical commutation implies the existence of
non-zero equal-time quantum correlations, at $\tau=0$. Explicitly,
$C^{\text{qu}}_{\gamma^{z},\gamma^{s}}(0)=i/2$. Remarkably, this persists in
spite of the macroscopic spatial separation of the two modes.

However, non-split MBs satisfy the same identity. To distinguish between split
and non-split pairs, it is necessary to go beyond $\tau=0$. For
$0<\tau<t_{N}$, we have

$\displaystyle C^{\text{ss}}_{\gamma^{z},\gamma^{s}}(\tau)$ $\displaystyle=$
$\displaystyle\text{tr}[\gamma^{z}(\tau)\gamma^{s}(0)\rho_{\text{ss}}]$
$\displaystyle\simeq$
$\displaystyle\text{tr}[\gamma^{z}(0)\gamma^{s}(0)\rho_{\text{ss}}]=C^{\text{ss}}_{\gamma^{z},\gamma^{s}}(0).$

On the other hand, for $\tau<0$ we have

$C^{\text{ss}}_{\gamma^{z},\gamma^{s}}(\tau)=\text{tr}[\gamma^{z}(0)\gamma^{s}(|\tau|)\rho_{\text{ss}}]\not\simeq
C^{\text{ss}}_{\gamma^{z},\gamma^{s}}(0).$

If instead the MBs were non-split, we would additionally have
$\gamma^{s}(t)\simeq\gamma^{s}(0)$ for $t<t_{N}$. So,
$C^{\text{ss}}_{\gamma^{z},\gamma^{s}}(\tau)\simeq
C^{\text{ss}}_{\gamma^{z},\gamma^{s}}(0)$. Therefore, split and non-split MBs
may be distinguished by asymmetries in the associated steady state correlation
function around $\tau=0$, reflecting the fact that split MBs are approximately
stationary for $0<\tau<t_{N}$, whereas non-split MBs are approximately
stationary for $-t_{N}<\tau<t_{N}$. Two remarks are in order: (i) While we
have treated the full steady state correlation function explicitly, the same
conclusions hold for both the classical and the (state-independent) quantum
contributions of the two-time correlation function. (ii) Per the properties of
MBs, the stationarity of the correlation functions in each case will become
more pronounced as system size is increased.

To exemplify the distinctions, we focus on the DBKC and the PDMC as
representative examples of these two classes. Moreover, by leveraging the
mapping described in App. B.3 between these two models [Eq. (66)] allows us to
directly compare topological phases. First, to distinguish the two models, we
relabel the quantities $J$, $\Delta$, and $\mu$ in the PDMC by $J_{F}$,
$\Delta_{F}$, and $\mu_{F}$, respectively; then, we identify $J_{F}=\Delta/2,$
$\Delta_{F}=J/2$, $\mu_{F}=-\kappa$, and take $\mu=\Gamma=0$. Let
$(\gamma_{L}^{c},\gamma_{R}^{s})$ and $(\gamma_{L}^{s},\gamma_{R}^{z})$ denote
the (split) MB pairs of the DBKC, and $(\gamma_{L},\gamma_{R})$ the (non-
split) MB pair of the PDMC. The above parameter identification yields

$\gamma_{L}=\gamma_{L}^{z},\quad\quad\gamma_{R}=\gamma_{R}^{s},$

which may be directly verified in the case $J=\Delta$. In particular, the
second MB pair of the DBKC $(\gamma_{L}^{s},\gamma_{R}^{z})$ are not
approximate ZMs, nor Weyl SGs, in the PDMC. This has several implications for
certain two-time correlation functions. Since the steady states of these two
models may differ in meaningful ways, we can directly compare the state-
independent quantum correlations of the MBs. Our general analysis above
predicts that the DBKC correlation function
$C^{\text{qu}}_{\gamma_{L}^{z},\gamma_{R}^{s}}(\tau)$ will be asymmetric about
zero and increasingly stationary in the positive $\tau$ direction as $N$
increases. On the contrary, the correlation function
$C^{\text{qu}}_{\gamma_{c},\gamma_{R}}(\tau)$ for the PDMC should be symmetric
and increasingly stationary in both the positive and negative $\tau$ direction
as $N$ increases. All of these predictions are numerically verified in Fig. 9.

Figure 9: (a) A MB correlation function for the DBKC. (b) A MB correlation
function for the PDMC. (c) A different MB correlation function for the DBKC.
(d) A correlation function for the same operators in (c) but for the PDMC. In
all cases, the modes are normalized so that canonical commutation relations
hold at $\tau=0$. Note that $\gamma_{L}^{z}=\gamma_{L}$,
$\gamma_{R}^{s}=\gamma_{R}$, $\gamma_{R}^{z}=\chi_{R}$, and
$\gamma_{L}^{s}=\chi_{L}$ at $\tau=0$, when the parameter mapping discussed in
the main text is applied.

A further distinction between these two models is identified by noting that
the second MB pair in the BKC has no analogue in the PDMC. Explicitly, the
operators $\gamma_{L}^{s}$ and $\gamma_{R}^{z}$, when mapped to the PDMC, are
neither approximate ZMs nor Weyl SGs. To distinguish which model we are
working in, let $\gamma_{L}^{s}\mapsto\chi_{L}$ and
$\gamma_{R}^{z}\mapsto\chi_{R}$ denote the image of the BKC’s second MB pair
in the PDMC. As argued above, the correlation function
$C_{\gamma_{R}^{z},\gamma_{L}^{s}}(\tau)$ will become more and more stationary
for $\tau\geq 0$ as $N\to\infty$. On the contrary, no such argument applies to
$C_{\chi_{R},\chi_{L}}(\tau)$ and so we generally expect exponentially
decaying correlations.

####  VII.2.2 Signatures of Dirac edge bosons

We have thus far focused on observable signatures of MBs. How, if at all, do
the correlation signatures we identified above change in the presence of
number symmetry? That is, can we detect the DBs of Sec. VI.2? Consider the
elementary steady state correlator $C_{i,j}^{\text{ss}}(\tau)\equiv
C_{\Phi_{i},\Phi_{j}^{\dagger}}^{\text{ss}}(\tau).$ This correlator is
elementary in the sense that any correlation function of linear observables
can be written as a linear combination of the above:

$C_{\alpha,\beta^{\dagger}}^{\text{ss}}(\tau)=\sum_{i,j}c^{i,j}_{\alpha,\beta^{\dagger}}C_{i,j}^{\text{ss}}(\tau),$

with $c^{i,j}_{\alpha,\beta^{\dagger}}$ determined by the coefficients of
$\Phi_{i}$ and $\Phi_{j}^{\dagger}$ in the definitions of $\alpha$ and
$\beta^{\dagger}$. Explicitly,
$c^{i,j}_{\alpha,\beta^{\dagger}}=\alpha_{i}\beta_{j}^{*}$, with $\alpha_{i}$
and $\beta_{j}$ elements of $\vec{\alpha}$ and $\vec{\beta}$.

Number symmetry ensures that $\mathcal{L}$ commutes with the superoperator
defined by the action $\rho\mapsto e^{i\theta\widehat{N}}\rho
e^{-i\theta\widehat{N}}$, with $\widehat{N}$ the total boson number operator
and $\theta\in\mathbb{R}$. Combining this with uniqueness of the steady state
immediately yields $[\rho_{\text{ss}},\widehat{N}]=0$. This guarantees that
the “off-diagonal” elementary correlators, i.e., correlators of the form
$\braket{a_{i}^{\dagger}(\tau)a_{j}^{\dagger}(0)}_{\text{ss}}$ and their
Hermitian conjugate counterparts, vanish. Mathematically,
$C_{i,j}^{\text{ss}}(\tau)=0$ for $i\leq N<j$ and $j\leq N<i$. In fact, the
off-diagonal, state-independent quantum correlations always vanish. This
follows because number symmetry guarantees

$\displaystyle a_{i}(\tau)=\sum_{j=1}^{N}d_{ij}(\tau)a_{j}(0),$

for some time-dependent coefficients
$d_{ij}(\tau)=2C_{i,j}^{\text{qu}}(\tau)$. This observation, combined with
canonical commutation relations, ensures that $[a_{i}(\tau),a_{j}(0)]=0$, for
all $\tau$. The equivalent statement in the quadrature basis is
$C^{\text{qu}}_{x_{j},p_{i}}(\tau)=C_{x_{i},p_{j}}^{\text{qu}}(-\tau)$. Since
this restatement is in terms of correlations between observables, it is
directly accessible in principle.

With this, we can characterize topologically metastable, number-symmetric
chains through the vanishing of their off-diagonal two-time correlation
functions, in conjunction with long-lived correlations / divergent zero
frequency power-spectral peaks. This behavior is reflected in Fig. 10. The
quantum power spectra $S_{a_{1}^{\dagger},a_{N}}(\omega)$ display exponential
divergence at zero frequency in both the DNS chain and the DBKC. However, the
off-diagonal spectra $S_{a_{1}^{\dagger},a_{N}^{\dagger}}(0)$ are exactly zero
for the DNS chain and diverging exponentially for the DBKC.

Figure 10: (Color online) Various quantum power spectra for the DNS chain and
the DBKC in their respective topologically-metastable regimes. The DNS shows
exponential divergence of the zero frequency peak
$S^{\text{qu}}_{a_{1}^{\dagger},a_{N}}(0)$ (green circles) and the vanishing
of the off-diagonal spectra $S^{\text{qu}}_{a_{1}^{\dagger},a_{N}^{\dagger}}$
(black squares). To contrast, the DBKC exhibits exponential growth of the off-
diagonal spectra (black diamonds). The parameters for the DNS are $J_{+}=1$,
$J_{-}=0.25$, and $\kappa=0.3$, while $J=2$, $\Delta=0.5$, $\mu=0$,
$\kappa=0.3$, and $\Gamma=0$ for the DBKC. These choices ensure an isospectral
relationship between the two models. An exponential fit is shown as a green
dashed line while $0$ is emphasized with a black dashed line.

##  VIII Conclusions and outlook

In this work, we have uncovered several convincing signatures of SPT-like
physics in a class of quadratic bosonic systems undergoing Markovian
dissipation. In this sense, we have progressed the search for SPT phases in
free bosonic systems by a great deal. The phenomenon of topologically non-
trivial dynamical metastability provides clear evidence that SPT phases may be
a transient, dynamical phenomena for non-interacting bosons. Specifically,
topologically metastable systems (i) host pairs of edge modes that bear an
extremely close resemblance to the fermionic edge modes of topological
superconductors and insulators, in a sense we have made mathematically
precise; and (ii) exhibit a sort-of “ground state degeneracy” in the form of a
manifold of Weyl-displaced quasi-steady states. Remarkably, a key ingredient
for topological metastability – namely, bulk-instability – is explicitly
forbidden in quadratic fermionic systems. Utilizing various new techniques, we
have produced and characterized several non-trivial models that sample the
expansive realm of possible instances of topological metastability.

With a robust theory of topological metastability at our fingertips, there are
a number of natural directions for future research. The most compelling next
step is to address the experimental accessibility of the phenomena we predict.
Potentially, there exists a number of promising proposals and experimental
implementations of related systems in cavity- and circuit-QED platforms [23,
32, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106], microlasers and
ring resonators [107, 108], optomechanical systems [109, 110], and vibronic
lattices [15]. While most of the contributions entering the driven-dissipative
models we considered have been individually realized, the dissipative hopping
and the (coherent and dissipative) bosonic pairing mechanisms are, arguably,
more exotic and challenging from a practical standpoint. Notably, Ref. 42
proposes a method for realizing dissipative hopping, while realizations of
both coherent and dissipative pairing have been proposed via three-wave mixing
with suitably tuned couplings with auxiliary modes [111, 112, 113, 114]. Once
a topologically metastable chain is successfully implemented, the two-time
correlation signatures we identified can, in principle, be accessed via their
connection to coherence functions. While significant challenges are likely to
arise, it is our hope that the possibility of validating our predictions of
tight analogues of fermionic zero modes in free bosons will spur a focused
effort across the community.

Beyond experimental detection of topological metastability, several
outstanding theoretical questions are also well-worth of further
investigation. One avenue involves extending the dynamical theory of
metastability into the realm of fermionic systems. While topological features
of said systems are fairly well-understood, we believe that our
pseudospectral-based metastability framework can provide major utility in
explaining anomalous relaxation [36] and cutoff phenomena [56], among other
things. Specifically, while dynamically metastable phases are explicitly
forbidden in fermions (since dynamical instabilities are [115, 28],
anomalously relaxing phases are not. Further exploration into the splitting of
zero modes and symmetry generators in fermionic systems are also of interest.
Still concerning bosons, addressing higher dimensions is a major necessary
extension of our framework. For example, what is the dissipative bosonic
analogue of the surface bands so prevalent in two-dimensional topological
insulators and superconductors, such as quantum Hall systems or the $p+ip$
superconductor? More pressingly, to what extent are the signatures we have
uncovered actual consequences of a genuine (one-dimensional) SPT phase? That
is, how do many-body symmetries factor into the theory of topological
metastability and to what extent do transitions into a topological metastable
regimes constitute genuine (dissipative) phase transition? This final question
is particularly interesting in light of recent work [116], which appears to
exclude criticality for stable 1D QBLs. Moreover, answering it would further
elucidate precisely how closely the parallel we are drawing between
topologically metastable QBLs and topologically non-trivial QFHs holds.

Finally, drawing analogy with Majorana-based quantum computing proposals, it
is natural to ask whether Majorana bosons could find utility in the context of
continuous-variable quantum information processing. In particular, central to
quantum computation schemes employing Gottesman-Kitaev-Preskill [117, 118,
119] is a pair of canonically conjugate quadratures. From here, the Gottesman-
Kitaev-Preskill code is built from fixed displacement operators within the
associated phase space and the logical states are built from the (ideally)
infinitely squeezed eigenstates of these operators. In our topological
metastability paradigm, we are provided with two canonically conjugate
quadratures that have an additional non-trivial property: they can be
macroscopically separated in space. This fact, which arises due to their
topological origin, affords an extra degree of robustness that may prove
beneficial for such applications. Notably, in the non-split case, the MBs
generate (orthogonal) phase-space displacements that leaves the overall
dynamics invariant to an arbitrarily high degree of precision (as set by the
system size). We look forward to address some of these questions in future
studies.

## Acknowledgements

It is a pleasure to thank Mariam Ughrelidze, Michiel Burgelman, Mattias
Fitzpatrick, and Joshuah Heath for valuable discussions. Work at Dartmouth was
partially supported by the US National Science Foundation through Grants No.
PHY-2013974 and No. OIA-1921199, the US Department of Energy, Office of
Science, Office of Advanced Scientific Computing Research, under the
Accelerated Research in Quantum Computing (ARQC) program, and the Constance
and Walter Burke Special Projects Fund in Quantum Information Science.

##  Appendix A Symmetry-protected topological phases
of free fermions

Given their core relevance in motivating our work, in this appendix we review
the basics of SPT phases of fermionic matter, with emphasis on the distinctive
properties that topological zero-energy edge states enjoy in topological
superconductor and insulator models.

Focusing on the most basic topological classification, the tenfold way [2],
there are four protecting many-body symmetries (electron number, spin
rotations, the physical time-reversal operation for electrons, and a particle-
hole anti-unitary transformation that exchanges creation and annihilation
operators) and ten (out of sixteen) symmetry classes that can be topologically
non-trivial depending on the space dimension. In any space dimension, five out
of the ten classes are topologically non-trivial. What changes with the
dimension is which ones are, and the change happens in a predictable, periodic
pattern.

Let us focus on the non-trivial classes in 1D. From the point of view of the
bulk, two classes break up into two disconnected components characterized by a
topological Pfaffian invariant in one case, and an invariant that mixes the
Pfaffian and the winding number in the other case; three other classes break
up into a countable infinite of components labeled by a winding number. From
the point of view of the edge, the Pfaffian invariant is associated to a
boundary invariant that counts the number of Majorana ZMs per edge modulo two.
This is the class of the Majorana chain of Kitaev. The “mixture” invariant is
associated to a boundary invariant that counts the number of Dirac ZMs per
edge modulo two or, equivalently, the number of pairs of Majorana modes per
edge, modulo two. And finally, the winding number invariant is associated to a
boundary invariant that counts the number of Dirac ZMs per edge. Loosely
speaking, a Dirac ZM is a fermionic degree of freedom that is not Hermitian,
that is, described by a creation and conjugate annihilation pair of operators.
In contrast, a Majorana mode is a Hermitian fermionic degree of freedom.

The general picture is well-captured by two models: The Majorana chain of
Kitaev [6] (FKC), and the mean-field model of the electronic structure of
trans-polyacetylene introduced by Su, Schrieffer, and Heeger [5] (SSH). Let us
consider the Majorana chain first from the point of view of the boundary. The
model is characterized by the Hamiltonian

$H_{\text{FKC}}=-\sum_{j=1}^{N}\mu
c_{j}^{\dagger}c_{j}-\sum_{j=1}^{N-1}\left(Jc_{j}^{\dagger}c_{j+1}-\Delta
c_{j}^{\dagger}c_{j+1}^{\dagger}+\text{H.c.}\right),$

with $\mu,J,\Delta\in\mathbb{R}$. In the language of the previous paragraph,
there is one Dirac fermion per site. However, since the symmetry of fermion
number is broken anyways, it is more natural to analyze the model in terms of
Majorana fermions $\gamma_{l}$, $l=1,\dots,2N$, such that
$2c_{j}\equiv\gamma_{2j-1}+i\gamma_{2j}.$ In this language, the Hamiltonian
becomes (up to an additive constant)

$\displaystyle
H_{\text{FKC}}=-\frac{\mu}{2}\sum_{j=1}^{N}i\gamma_{2j-1}\gamma_{2j}$
$\displaystyle-\frac{J+\Delta}{2}\sum_{j=1}^{N-1}i\gamma_{2j-1}\gamma_{2j+2}-i\frac{\Delta-J}{2}\sum_{j=1}^{N-1}\gamma_{2j}\gamma_{2j+1}.$

Looking individually at the three terms, the first term is topologically
trivial and the other two are non-trivial in the sense that they support
topologically mandated edge modes. Whether their sum is trivial or not depends
on the competition among these three terms. Let us focus for concreteness on
the case $J=\Delta$. Then, the model is non-trivial for $|\mu|<|2J|$ and
trivial otherwise. The lines $|\mu|=|2J|$ are critical.

Let $\delta=\mu/2J$. In the non-trivial regime, $|\delta|<1$ and the operators

$\displaystyle\gamma_{L}\equiv{\cal
M}(N)\sum_{j=1}^{N}\delta^{j-1}\gamma_{2j-1},\quad\gamma_{R}\equiv{\cal
M}(N)\sum_{j=1}^{N}\delta^{N-j}\gamma_{2j},$

are candidate edge ZMs because

$\displaystyle i[H_{\text{FKC}},\gamma_{L}]$ $\displaystyle=\delta^{N}{\cal
M}(N)2J\,\gamma_{2N},$ (62) $\displaystyle i[H_{\text{FKC}},\gamma_{R}]$
$\displaystyle=\delta^{N}{\cal M}(N)2J\,\gamma_{1}.$ (63)

may potentially vanish exponentially fast with $N$ However, the normalization
factor ${\cal M}(N)$ depends on the size of the chain and so its behavior must
be ascertained. To determine it, notice that

$\displaystyle\\{\gamma_{L},\gamma_{R}\\}=0,\quad\gamma_{L}^{2}=\gamma_{R}^{2}={\cal
M}(N)^{2}\sum_{j=1}^{N}\delta^{2(j-1)}.$

Hence, this pair of operators will satisfy the algebra of Majorana fermions
provided that

$\displaystyle{\cal M}(N)^{-1}=\sqrt{\frac{1-\delta^{2(N-1)}}{1-\delta^{2}}}.$
(64)

In particular, ${\cal M}(N)$ is bounded as a function of system size if
$|\mu|<2J$. It follows that the commutators of Eqs. (62)-(63) do indeed vanish
exponentially fast in the size of the system.

The two Majorana fermions at zero energy can be combined to form a single
Dirac fermion at the expense of locality. The resulting Dirac fermion is very
peculiar: it is neither localized nor uniformly spread over the length of the
chain. The degeneracy of the ground energy level is two, depending on whether
this Dirac fermion is or not present in the ground state. As a consequence,
the fermion number operator is distributed over the odd or even integers only
in one or the other ground state: the two ground states differ by their
fermionic parity.

Let us consider next the SSH model. The Hamiltonian is

$\displaystyle H_{\text{SSH}}=$
$\displaystyle\,t_{2}\sum_{\sigma}\sum_{j=1}^{M}(c_{2j-1,\sigma}^{\dagger}c_{2j,\sigma}+\text{H.c.})$
$\displaystyle+$
$\displaystyle\,t_{1}\sum_{\sigma}\sum_{j=1}^{M-1}(c_{2j,\sigma}^{\dagger}c_{2j+1,\sigma}+\text{H.c.}).$

Here, $\sigma$ is the label associated to the spin degree of freedom of the
electron. This Hamiltonian belongs to the most symmetric symmetry class, BDI.
The structure of the edge is identical for spin up and spin down electrons.
Hence, let us relabel $c_{2j-1,\sigma}$ as $a_{j}$ and $c_{2j,\sigma}$ as
$b_{j}$, write $N$ for $M/2$, and focus on any one spin direction. The
Hamiltonian of interest is then

$H=\
t_{2}\sum_{j=1}^{N}(a_{j}^{\dagger}b_{j}+\text{H.c.})+t_{1}\sum_{j=1}^{N-1}(b_{j}^{\dagger}a_{j+1}+\text{H.c.}).$

This Hamiltonian is sometimes also called the Peierls model. The first term is
topologically trivial and the second one is non-trivial, in the same sense as
before. Their sum is non-trivial if $|t_{1}|>|t_{2}|$. The line
$|t_{1}|=|t_{2}|$ is critical, as one expects from the fact that it marks the
transition from one topological phase to another within a symmetry class.

Let $\delta=-t_{2}/t_{1}$. In the nontrivial regime, $\delta<1$ and the
operator

$\displaystyle a_{L}={\cal M}(N)\sum_{j=1}^{N}\delta^{j-1}a_{j}$

and its Hermitian conjugate are candidate ZMs, because

$[H,a_{L}]=\delta^{N}{\cal M}(N)t_{1}b_{N}$

may vanish exponentially fast as a function of the size of the system. Let us
confirm that the normalization factor ${\cal M}(N)$ is bounded as a function
of $N$. Notice that

$\displaystyle
a_{L}^{2}=a_{L}^{\dagger\,2}=0,\quad\\{a_{L},a_{L}^{\dagger}\\}={\cal
M}(N)^{2}\sum_{j=1}^{N}\delta^{2(j-1)}.$

Hence, the pair $a_{L},a_{L}^{\dagger}$ constitutes a Dirac ZM localized on
the left edge provided that ${\cal M}(N)$ is given by Eq. (64). A similar
calculation shows that the operator

$\displaystyle b_{R}={\cal M}(N)\sum_{j=1}^{N}\delta^{j-1}b_{N-j+1}$

and its Hermitian conjugate constitutes a Dirac ZM localized on the right
edge. Since, in addition,

$\displaystyle\\{a_{L},b_{R}\\}=0,\quad\\{a_{L},b_{R}^{\dagger}\\}=0,$

the left and right ZMs taken together are two Dirac fermions. Hence, the
degeneracy of the ground level is four. Notice that, unlike the zero-energy
Dirac fermions of the FKC, those of the SSH chain are localized on the edges
of the chain.

##  Appendix B Additional technical results

###  B.1 Pseudospectra, singular values, and the
doubled-matrix approach

In this appendix, we will recount the connection between the 2-norm
pseudospectra and singular values, and establish the bridge to the so-called
doubled-matrix approach utilized prominently in topological amplification
contexts [42, 43].

Let $s_{j}(\mathbf{X})$, $j=1,\ldots,n$, denote the $n$ singular values of an
$n\times n$ matrix $X$. If $\bf{X}$ is invertible, it is well-known that
$\lVert\mathbf{X}^{-1}\rVert_{2}=1/s_{\text{min}}(\mathbf{X})$ with
$s_{\text{min}}(\mathbf{X})$ the minimal singular value of $\mathbf{X}$.
Immediately we can see that the 2-norm $\epsilon$-pseudospectrum of a matrix
$\mathbf{X}$ is precisely the set of complex numbers $\lambda$ such that
$s_{\text{min}}(\lambda\mathds{1}-\mathbf{X})<\epsilon$.

The doubled matrix approach starts by defining the “doubled-Hamiltonian”

$\overline{\mathbf{X}}(\omega)=\begin{bmatrix}0&\omega\mathds{1}-\mathbf{X}\\\
\omega\mathds{1}-\mathbf{X}^{\dagger}&0\end{bmatrix},$

with $\omega$ a real frequency and $\mathbf{X}$ typically taken to be the
dynamical matrix of interest. One may show that the built-in chiral structure
of $\overline{X}(\omega)$ implies a $\pm$ symmetry of the spectrum, and
moreover the eigenvalues are precisely $\pm
s_{j}(\omega\mathds{1}-\mathbf{X})$. Often, “midgap eigenvalues”, i.e.
exponentially-small-in-$N$ eigenvalues sandwiched between the positive and
negative bulk-bands, are of interest. Importantly, these correspond to the
minimal singular values of $\omega\mathds{1}_{n}-\mathbf{X}$. In particular,
if $\overline{\mathbf{X}}(\omega)$ has an exponentially small midgap
eigenvalue, then $\omega\mathds{1}-\mathbf{X}$ has a vanishing minimal
singular value indicating that $\omega$ exists within the pseudospectrum of
$\mathbf{X}$ at sufficiently large system size. In short, if
$\overline{\mathbf{X}}(\omega)$ has midgap eigenvalues, then $\omega$ is in
the pseudospectrum of $\mathbf{X}$.

###  B.2 A canonical correspondence
between approximate ZMs and SGs

Here, we prove that, under generic conditions, for each approximate SG there
is a canonically conjugate ZM, thereby generalizing Theorem 1 to the
approximate case.

Theorem 2. Given an arbitrary QBL with dynamical matrix $\mathbf{G}$, let
$\gamma^{s}=\widehat{\vec{\gamma}{}^{s}}$ be an approximate symmetry generator
of accuracy $\epsilon$ (in the sense that $\vec{\gamma}^{s}$ satisfies
$\lVert\mathbf{G}\vec{\gamma}^{s}\rVert<\epsilon$ for some $\epsilon>0$).
Then, if the matrix
$\mathbf{G}^{\prime}\equiv\mathbf{G}-\vec{\alpha}\vec{\gamma}^{s}{}^{\dagger}/\lVert\vec{\gamma}^{z}\rVert$
hosts only Jordan chains of length one at $0$, there exists a canonically
conjugate approximate zero mode $\gamma^{z}=\widehat{\vec{\gamma}^{z}}$ of
accuracy $\epsilon^{\prime}=\epsilon\lVert\vec{\gamma}^{z}\rVert^{2}$. The
converse holds as well.

###### Proof.

We will assume that the approximate symmetry generator is not exact. If it is,
apply Theorem 1. By the assumptions of the theorem, we have that
$\mathbf{G}\vec{\gamma}^{s}=\vec{\alpha}$ for some vector $\vec{\alpha}$ with
$\lVert\vec{\alpha}\rVert<\epsilon$. Hermiticity of $\gamma^{s}$ implies that
the Nambu representation satisfies
$\vec{\gamma}^{s}=-\bm{\tau}_{1}\vec{\gamma}^{s}{}^{*}$. Using this, and the
fact that $\bm{\tau}_{1}\mathbf{G}^{*}\bm{\tau}_{1}=-\mathbf{G}$, we have

$\displaystyle\bm{\tau}_{1}\vec{\alpha}^{*}=\bm{\tau}_{1}\mathbf{G}^{*}\vec{\gamma}^{s}{}^{*}-\mathbf{G}\bm{\tau}_{1}\vec{\gamma}^{s}{}^{*}=\mathbf{G}\vec{\gamma}^{s}=\vec{\alpha}.$
(65)

Now, consider the matrix
$\mathbf{G}^{\prime}\equiv\mathbf{G}-\vec{\alpha}\vec{\gamma}^{s}{}^{\dagger}/\lVert\vec{\gamma}^{s}\rVert^{2}$.
We claim that $\mathbf{G}^{\prime}$ may be interpreted as a dynamical matrix
of some other QBL and that $\mathbf{G}^{\prime}\vec{\gamma}^{s}=0$. The first
claim follows $\mathbf{G}^{\prime}$ obeys the only constraint set on dynamical
matrices:

$\displaystyle\bm{\tau}_{1}\mathbf{G}^{\prime}{}^{*}\bm{\tau}_{1}$
$\displaystyle=\bm{\tau}_{1}\mathbf{G}^{*}\bm{\tau}_{1}-\bm{\tau}_{1}\vec{\alpha}^{*}\vec{\gamma}^{s}{}^{T}\bm{\tau}_{1}$
$\displaystyle=-\mathbf{G}-(\bm{\tau}_{1}\vec{\alpha})^{*}(\bm{\tau}_{1}\vec{\gamma}^{s}{}^{*})^{\dagger}$
$\displaystyle=-\mathbf{G}+\vec{\alpha}\vec{\gamma}^{s}{}^{\dagger}=-\mathbf{G}^{\prime},$

where we have used Eq. (65) in the third equality. The second claim is
verified directly

$\displaystyle\mathbf{G}^{\prime}\vec{\gamma}^{s}=\mathbf{G}\vec{\gamma}^{z}-\vec{\alpha}\frac{\vec{\gamma}^{s}{}^{\dagger}\vec{\gamma}^{s}}{\lVert\vec{\gamma}^{s}\rVert^{2}}=\vec{\alpha}-\vec{\alpha}=0.$

Now, any QBL with dynamical matrix $\mathbf{G}^{\prime}$ has $\gamma^{s}$ as
an exact SG. We may then directly apply Theorem 1 to find a canonically
conjugate exact ZM $\gamma^{z}$. In Nambu space, this means there is a vector
$\vec{\gamma}^{z}=-\bm{\tau}_{1}\vec{\gamma}^{z}{}^{*}$ such that
$\widetilde{\mathbf{G}}^{\prime}\vec{\gamma}^{z}=0$, with

$\displaystyle\widetilde{\mathbf{G}}^{\prime}=\bm{\tau}_{3}\mathbf{G}^{\prime}{}^{\dagger}\bm{\tau}_{3}=\widetilde{\mathbf{G}}-\bm{\tau}_{3}\vec{\gamma}^{s}\vec{\alpha}^{\dagger}\bm{\tau}_{3},$

and $\vec{\gamma}^{s}{}^{\dagger}\bm{\tau}_{3}\vec{\gamma}^{z}=i$, so that
$[\gamma^{s},\gamma^{z}]=i$. It follows that

$\displaystyle\lVert\widetilde{\mathbf{G}}\vec{\gamma}^{z}\rVert$
$\displaystyle=\lVert\widetilde{\mathbf{G}}^{\prime}\vec{\gamma}^{z}+\frac{\bm{\tau}_{3}\vec{\gamma}^{s}\vec{\alpha}^{\dagger}\bm{\tau}_{3}}{\lVert\vec{\gamma}^{s}\rVert^{2}}\vec{\gamma}^{z}\rVert$
$\displaystyle=\frac{1}{\lVert\vec{\gamma}^{s}\rVert^{2}}\lVert\bm{\tau}_{3}\vec{\gamma}^{s}\vec{\alpha}^{\dagger}\bm{\tau}_{3}\vec{\gamma}^{z}\rVert=\frac{|\vec{\alpha}^{\dagger}\bm{\tau}_{3}\vec{\gamma}^{z}|}{\lVert\vec{\gamma}^{s}\rVert}.$

Using the Cauchy-Schwarz inequality, we may upper-bound the numerator of the
right hand-side by
$\lVert\vec{\alpha}\rVert\lVert\vec{\gamma}^{z}\rVert<\epsilon\lVert\vec{\gamma}^{z}\rVert$.
Thus,
$\lVert\widetilde{\mathbf{G}}\vec{\gamma}^{z}\rVert<\epsilon{\lVert\vec{\gamma}^{z}\rVert}/{\lVert\vec{\gamma}^{s}\rVert}.$
An additional application of the Cauchy-Schwarz inequality to the identity
$\vec{\gamma}^{s}{}^{\dagger}\bm{\tau}_{3}\vec{\gamma}^{z}=i$ yields
$1/\lVert\vec{\gamma}^{s}\rVert<\lVert\vec{\gamma}^{z}\rVert$, so that

$\displaystyle\lVert\widetilde{\mathbf{G}}\vec{\gamma}^{z}\rVert<\epsilon\lVert\vec{\gamma}^{z}\rVert^{2},$

as claimed. The converse holds by simply replacing $\mathbf{G}$ with
$\widetilde{\mathbf{G}}$ and $\vec{\gamma}^{s}$ with $\vec{\gamma}^{z}$. ∎

We remark that canonically conjugate modes can always be rescaled in a natural
way to make the accuracies $\epsilon$ and $\epsilon^{\prime}$ equal.
Generally, if $[\gamma^{s},\gamma^{z}]=i$, then
$[\gamma^{s}{}^{\prime},\gamma^{z}{}^{\prime}]=i$, with
$\gamma^{s}{}^{\prime}=\mathcal{M}\gamma^{s}$ and
$\gamma^{z}{}^{\prime}=\gamma^{z}/\mathcal{M}$ for any $\mathcal{M}>0$. Taking
$M=\lVert\vec{\gamma}^{z}\rVert$ provides approximate SGs and ZMs of accuracy
$\mathcal{M}\epsilon$. The most interesting examples are then those in which
$\epsilon\ll 1/\mathcal{M}$.

###  B.3 An isospectral mapping between the purely dissipative and the
dissipative bosonic Kitaev chains

Here, we establish an isospectral mapping between the PDMC and the DBKC. Let
$\mathbf{G}^{\text{DBKC}}(\kappa,J,\Delta)$ and
$\mathbf{G}^{\text{PDC}}(\mu_{F},J_{F},\Delta_{F})$ denote the dynamical
matrices of the DBKC and the PDMC under OBCs, respectively. Note that we have
distinguished the FKC hopping, pairing, and onsite potential with a subscript
$F$. Define the momentum-space translation operator

$\bm{\Lambda}(\delta k)\equiv\text{diag}(e^{-i\delta
k},\ldots,e^{-iN{\color[rgb]{0,0,0}\delta k}})\otimes\mathds{1}_{2}.$

With this, one may verify that

$\bm{\Lambda}\left(\frac{\pi}{2}\right)\mathbf{G}^{\text{PDC}}\left(0,-\frac{J}{2},-\frac{\Delta}{2}\right)\bm{\Lambda}^{-1}\left(\frac{\pi}{2}\right)=i\mathbf{G}^{\text{DBKC}}(0,J,\Delta).$
(66)

That is, the dynamical matrices are unitarily equivalent up to a phase when
$\mu_{F}=\kappa=0$. In particular, the $\mu_{F}=0$ OBC rapidity spectrum for
the PDMC is equal to $i$ times that of the DBKC, with the identifications
$\kappa=0$, $J_{F}\leftrightarrow-J/2$, and
$\Delta_{F}\leftrightarrow-\Delta/2$. This allows us to establish Eq. (47)
from the known spectral properties of the DBKC [24].

###  B.4 On the numerical determination of Majorana bosons

In this appendix, we present a procedure for numerically computing MBs in
topologically metastable QBLs. The basic inputs are (i) a family of dynamical
matrices $\mathbf{G}_{N}$, one for each system size $N$, corresponding to a
topologically metastable QBL under OBCs; and (ii) a sufficiently small (in a
sense to be explained) accuracy parameter $\epsilon$. This procedure outputs
pairs of canonically conjugate $\epsilon$-pseudoeigenvectors corresponding to
$\epsilon$-pseudoeigenvalue 0 for sufficiently large $N$ (to be explained).
The pseudospectra considered are $2$-norm pseudospectra.

To begin, we define a positive-semidefinite matrix
$\mathbf{R}_{N}\equiv\mathbf{G}_{N}^{\dagger}\mathbf{G}_{N}$. The eigenvalues
of $\mathbf{R}$ are the squares of the singular values of $\mathbf{G}_{N}$. In
the language of App. B.1, topological metastability ensures that
$\mathbf{R}_{N}$ possesses a finite number (say, $m$) of midgap modes, i.e.,
eigenvectors corresponding to eigenvalues that decay to $0$ as $N$ increases.
The remaining $2N-m$ eigenvalues are bounded below by a gap parameter
$\lambda_{\text{gap}}$. Fix $0<\epsilon<\lambda_{\text{gap}}^{1/2}$ and let
$N_{0}$ be the smallest $N$, such that the $m$ midgap eigenvalues are less
than $\lambda_{\text{gap}}$. By construction, the $m$ normalized midgap
eigenvectors, which we will label $\vec{u}_{j}$, $j=1,\ldots,m$, will be
$\epsilon$-pseudoeigenvectors of $\mathbf{G}_{N}$ corresponding to
$\epsilon$-pseudoeigenvalue zero for all $N>N_{0}$. Explicitly, we have

$\lVert\mathbf{G}_{N}\vec{u}\rVert_{2}=\left(\vec{u}^{\dagger}\mathbf{R}_{N}\vec{u}\right)^{1/2}=s_{j}<\epsilon,$

with $s_{j}$ the singular value corresponding to $\vec{u}_{j}$. Equivalently,
$s_{j}^{2}=\lambda_{j}$ the corresponding midgap eigenvalue of
$\mathbf{R}_{N}$.

Despite being pseudoeigenvectors of $\mathbf{G}_{N}$, these need not yet
provide approximate SGs. In particular, they need not generate a Hermitian
linear form. For this, we leverage the property
$\mathcal{C}\mathbf{G}_{N}\mathcal{C}^{-1}=-\mathbf{G}_{N}$. This implies
$\mathcal{C}\mathbf{R}_{N}\mathcal{C}^{-1}=\mathbf{R}_{N}$. Thus, if
$\vec{u_{j}}$ is an eigenvector of $\mathbf{R}$, so is
$\mathcal{C}\vec{u}_{j}$ with the same eigenvalues. In particular, the set
$\\{\vec{u}_{j}\\}$ is invariant under $\mathcal{C}$. If
$\mathcal{C}\vec{u}_{j}=e^{i\theta}\vec{u}_{j}$, then we can normalize to
obtain a new vector $\vec{\gamma}_{j}^{s}=-e^{i\theta}\vec{u}_{j}$ that
defines a Hermitian approximate SG $\widehat{\vec{\gamma}}^{s}$. If there are
degenerate singular values, then the restriction of $\mathcal{C}$ to this
subspace can be diagonalized to obtain a basis of vectors
$\vec{\gamma}_{j}^{s}$ that each define Hermitian SGs.

Constructing the approximate ZMs follows analogously, but with $\mathbf{G}$
replaced with $\widetilde{\mathbf{G}}$. Let
$\\{\vec{\gamma}_{j}^{s}{}^{\prime}\\}$ and
$\\{\vec{\gamma}_{j}^{z}{}^{\prime}\\}$ be the outputs of these two
procedures. It follows that canonical commutation relations can be ensured if
the matrix
$\mathbf{F}_{jk}=\vec{\gamma}_{j}^{s}{}^{\prime}{}^{\dagger}\bm{\tau}_{3}\vec{\gamma}_{k}^{z}{}^{\prime}$
is diagonalizable. In all examples considered in this paper, this is the case.

##  Appendix C The time-evolution of bosonic parity in a Markovian system
prepared in a cat state

We present here the explicit calculation of the cat-state parity dynamics
under the pure SS DBKC. As we noted in the main text, this model can be seen
as a set of independent damped quantum harmonic oscillators in the normal-mode
basis. This decoupling allows us to reduce the problem to that of computing
the parity dynamics of a single-mode cat-state under damped harmonic motion.
The multimode generalization then follows naturally.

Let $\ket{\alpha}$, $\alpha\in\mathbb{C}$, denote a single-mode coherent state
and define the single-mode cat state
$\ket{\mathcal{C}_{\phi}(\alpha)}\equiv\mathcal{N}_{\phi}(\alpha)\left(\ket{\alpha}+e^{i\phi}\ket{-\alpha}\right)$,
with $\mathcal{N}_{\phi}(\alpha)$ a normalization constant. We wish to compute
the expectation value of parity $P=e^{i\pi a^{\dagger}a}$ in the time-evolved
state $\rho_{\alpha,\phi}(t)$, resulting from the initial condition
$\rho(0)=\ket{\mathcal{C}_{\phi}(\alpha)}\bra{\mathcal{C}_{\phi}(\alpha)}$
under the QBL

$\mathcal{L}(\rho)=-i[\omega a^{\dagger}a,\rho]+2\kappa\Big{(}a\rho
a^{\dagger}-\frac{1}{2}\\{a^{\dagger}a,\rho\\}\Big{)}.$

Here, $\omega$ is the oscillator frequency and $\kappa$ is the damping rate.
We may compute $\rho_{\alpha,\phi}(t)$ exactly utilizing known results [83,
84]. We may write

$\displaystyle\rho_{\alpha,\phi}(t)=e^{t\mathcal{L}}\left(\rho_{\alpha,\phi}(0)\right)$
$\displaystyle=\mathcal{N}_{\phi}(\alpha)^{2}\left(\sigma_{\alpha}(t)+\sigma_{-\alpha}(t)+e^{i\phi}\chi_{\alpha}(t)+e^{-i\phi}\chi_{-\alpha}(t)\right),$

where

$\sigma_{\alpha}(t)\equiv
e^{t\mathcal{L}}(\ket{\alpha}\bra{\alpha}),\quad\chi_{\alpha}(t)\equiv
e^{t\mathcal{L}}(\ket{\alpha}\bra{-\alpha}).$

The terms $\sigma_{\pm\alpha}(t)$ can be quoted directly as

$\sigma_{\pm\alpha}(t)=\ket{\pm\alpha(t)}\bra{\pm\alpha(t)},\quad\alpha(t)=e^{-(\kappa+i\omega)t}\alpha(0).$

The terms $\chi_{\pm\alpha}(t)$ are slightly more complicated,

$\chi_{\alpha}(t)=D(t)\bigg{(}\sum_{k=0}^{\infty}\frac{(2e^{-\kappa
t}\sinh(\kappa
t))^{k}}{k!}a^{k}\chi_{\alpha}(0)a^{\dagger}{}^{k}\bigg{)}D(t)^{\dagger},$

with $D(t)\equiv e^{-(\kappa+i\omega)ta^{\dagger}a}$. Now,

$a^{k}\chi_{\alpha}(0)a^{\dagger}{}^{k}=\alpha^{k}\chi_{\alpha}(0)(-\alpha^{*})^{k}=(-|\alpha|^{2})^{k}\chi_{\alpha}(0),$

which leads us to

$\chi_{\alpha}(t)=\exp[-2|\alpha|^{2}e^{-\kappa t}\sinh(\kappa
t)]\,D(t)\chi_{\alpha}(0)D(t)^{\dagger}.$

The remaining time-dependence may be computed as

$D(t)\chi_{\alpha}(0)D(t)^{\dagger}\\!=\\!\exp[-2|\alpha|^{2}e^{-\kappa
t}\sinh(\kappa t)]\ket{\alpha(t)}\bra{-\alpha(t)}.$

Finally, we have

$\displaystyle\chi_{\alpha}(t)$ $\displaystyle=\exp[-4|\alpha|^{2}e^{-\kappa
t}\sinh(\kappa t)]\ket{\alpha(t)}\bra{-\alpha(t)}$ $\displaystyle\equiv
f_{\alpha}(t)\ket{\alpha(t)}\bra{-\alpha(t)},$

with $\chi_{-\alpha}(t)$ following accordingly.

With the exact time dependence of $\rho(t)$ computed, we can now evaluate the
expectation value of parity. A particularly useful identity is
$P\ket{\alpha}=\ket{-\alpha}$. Using this, we can compute
$P_{1}^{\alpha}(t)=\text{tr}[P\sigma_{\alpha}(t)]$ and
$P_{2}^{\alpha}(t)=\text{tr}[P\chi_{\alpha}(t)]$ to obtain

$\braket{P}\\!(t)=\mathcal{N}_{\phi}(\alpha)^{2}[P_{1}^{\alpha}(t)+P_{1}^{-\alpha}(t)+e^{i\phi}P_{2}^{\alpha}(t)+e^{-i\phi}P_{2}^{-\alpha}(t)].$

Proceeding, we find

$\displaystyle
P_{1}^{\alpha}(t)\\!=\\!\text{tr}[P\ket{\alpha(t)}\bra{\alpha(t)}]$
$\displaystyle=$ $\displaystyle\braket{\alpha(t)}{-\alpha(t)}$
$\displaystyle=$ $\displaystyle\exp\left(-2|\alpha|^{2}e^{-2\kappa t}\right),$
$\displaystyle
P_{2}^{\alpha}(t)\\!=\\!f_{\alpha}(t)\text{tr}[P\ket{\alpha(t)}\bra{-\alpha(t)}]$
$\displaystyle=$ $\displaystyle f_{\alpha}(t)\braket{-\alpha(t)}{-\alpha(t)}$
$\displaystyle=$ $\displaystyle f_{\alpha}(t).$

Putting this all together yields the single-mode evolution,

$\braket{P}\\!(t)=\frac{e^{-2|\alpha|^{2}}+\cos(\phi)e^{-2|\alpha|^{2}(1-e^{-2\kappa
t})}}{1+\cos(\phi)e^{-2|\alpha|^{2}}},$

from which the multimode generalization given by Eq. (53) in the main text
follows.

##  Appendix D A dissipative double-winding bosonic chain

It is natural to wonder to what degree we can control the number of ZMs and
SGs on each edge of a topologically metastable chains. The properties of
Toeplitz matrices ensure that each band with a positive winding number yield a
left-localized ZM and a right-localized SG. Conversely, negative winding
numbers yield a right-localized ZM and left-localized SGs. In all the examples
discussed in the main text, the number of ZMs (SGs) on a particular edge has
been either zero or one. Surprisingly, we can write down a range-one model
with a single internal degree of freedom that has two ZMs (SGs) on a
particular edge. The Hamiltonian is given by

$\displaystyle
H_{\text{DW}}=\frac{i}{2}\sum_{j=1}^{N}\left(\Delta_{h}a_{j}^{\dagger}a_{j+1}^{\dagger}+\mu
a_{j}^{\dagger}{}^{2}-\text{H.c.}\right),$

with the usual BC conventions. Here, $\Delta_{h},\mu\in\mathbb{R}$ denote the
non-degenerate and degenerate parametric amplification strength, respectively.
The dissipator consists of three distinct contributions
$\mathcal{D}_{\text{DW}}\equiv\mathcal{D}_{-,0}+\mathcal{D}_{+,0}+\mathcal{D}_{p,1}$.
Explicitly,

$\displaystyle\mathcal{D}_{-,0}$ $\displaystyle=$ $\displaystyle
2\kappa_{-}\sum_{j=1}^{N}\mathcal{D}[a_{j},a_{j}^{\dagger}],\quad$
$\displaystyle\mathcal{D}_{+,0}$ $\displaystyle=$ $\displaystyle
2\kappa_{+}\sum_{j=1}^{N}\mathcal{D}[a_{j}^{\dagger},a_{j}],$
$\displaystyle\mathcal{D}_{p,1}$ $\displaystyle=$
$\displaystyle\frac{\Delta_{d}}{2}\sum_{j=1}^{N}\mathcal{D}[a_{j}^{\dagger},a_{j+1}^{\dagger}]-\mathcal{D}[a_{j+1}^{\dagger},a_{j}^{\dagger}]-(a_{j}^{\dagger}\leftrightarrow
a_{j}).$

In words, we have uniform onsite loss with strength $\kappa_{-}\geq 0$,
uniform onsite gain with strength $\kappa_{+}\geq 0$, and NN incoherent
pairing with strength $\Delta_{d}$. Positivity is the associated GKLS matrix
demands that $4\kappa_{+}\kappa_{-}\geq\Delta_{d}^{2}$. For simplicity, we
additionally take $\Delta_{d}\geq\Delta_{h}\geq 0$. If
$2\kappa_{\pm}\geq\Delta_{d}$, we can write a simple diagonal representation:

$\displaystyle\mathcal{D}_{\text{DW}}$
$\displaystyle=(2\kappa_{-}-\Delta_{d})\sum_{j=1}^{N}\mathcal{D}[a_{j}]+(2\kappa_{+}-\Delta_{d})\sum_{j=1}^{N}\mathcal{D}[a_{j}^{\dagger}]$
$\displaystyle+\frac{\Delta_{d}}{2}\sum_{j=1}^{N}\mathcal{D}[a_{j}-a_{j+1}^{\dagger}]+\frac{\Delta_{d}}{2}\sum_{j=1}^{N}\mathcal{D}[a_{j}+a_{j+1}^{\dagger}].$

The symbol of the dynamical matrix is given by
$-i\mathbf{g}_{\text{DW}}(k)=-\kappa\mathds{1}_{2}+(\mu+\Delta_{h}\cos(k)+i\Delta_{d}\sin(k))\bm{\sigma}_{1}$
with $\kappa\equiv\kappa_{-}-\kappa_{+}$. The rapidity bands are

$\lambda_{\pm}(k)=-\kappa\pm(\mu+\Delta_{h}\cos(k)+i\Delta_{d}\sin(k)).$

The winding numbers $\nu_{\pm}$ of the bands $\lambda_{\pm}(k)$ are given by

$\nu_{\pm}=\begin{cases}1&-1<\tilde{\kappa}\mp\tilde{\mu}<1,\\\
0&\text{otherwise,}\end{cases}$

with $\tilde{\kappa}\equiv\kappa/\Delta_{h}$ and
$\tilde{\mu}\equiv\mu/\Delta_{h}$. By tuning $\tilde{\kappa}$ and
$\tilde{\mu}$ it is possible to obtain a total winding $\nu_{+}+\nu_{-}$ equal
to $0$, $1$, or $2$.

## References

  * Blaizot and Ripka [1986] J.-P. Blaizot and G. Ripka,  _Quantum Theory of Finite Systems_ (Cambridge: MIT Press Cambridge, 1986) Chap. 2-3..
  * Chiu  _et al._ [2016] C. K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Alldridge  _et al._ [2020] A. Alldridge, C. Max, and M. R. Zirnbauer, Bulk-boundary correspondence for disordered free-fermion topological phases, Commun. Math. Phys. 377, 1761 (2020).
  * Alase  _et al._ [2023] A. Alase, E. Cobanera, G. Ortiz, and L. Viola, Wiener-Hopf matrix factorization approach to a bulk-boundary correspondence and stability conditions for topological ground-state degeneracy, Ann. Phys. (NY) 458, 169457 (2023).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * Kitaev [2001] A. Kitaev, Unpaired Majorana fermions in quantum wires, Phys. Usp. 44, 131 (2001).
  * Note [1] Importantly, for a number-conserving system, its phase is determined by spectral properties of the Hamiltonian and the filling fraction. In addition, while in one dimension the topological surface modes are isolated in the gap, in higher dimensions the surface modes are organized into surface bands that either connect the lowest bulk band to zero energy (for superconductors), or cross the Fermi energy.
  * Raghu and Haldane [2008] S. Raghu and F. D. M. Haldane, Analogs of quantum-Hall-effect edge states in photonic crystals, Phys. Rev. A 78, 033834 (2008).
  * Wang  _et al._ [2009] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Solijacić, Observation of unidirectional backscattering-immune topological electromagnetic states, Nature 461, 772 (2009).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 015006 (2019).
  * Lu  _et al._ [2018] Y.-K. Lu, P. Peng, Q.-T. Cao, D. Xu, J. Wiersig, Q. Gong, and Y.-F. Xiao, Spontaneous $\mathcal{T}$-symmetry breaking and exceptional points in cavity quantum electrodynamics systems, Science Bull. 63, 1096 (2018).
  * Xu  _et al._ [2020] Q.-R. Xu, V. P. Flynn, A. Alase, E. Cobanera, L. Viola, and G. Ortiz, Squaring the fermion: The threefold way and the fate of zero modes, Phys. Rev. B 102, 125127 (2020).
  * Note [2] No-go (iii) does not require translation invariance.
  * Peano and Schulz-Baldes [2018] V. Peano and H. Schulz-Baldes, Topological edge states for disordered bosonic systems, J. Math. Phys. 59, 031901 (2018).
  * Peano  _et al._ [2015] V. Peano, C. Brendel, M. Schmidt, and F. Marquardt, Topological Phases of Sound and Light, Phys. Rev. X 5, 031011 (2015).
  * Shindou  _et al._ [2013] R. Shindou, R. Matsumoto, S. Murakami, and J.I. Ohe, Topological chiral magnonic edge mode in a magnonic crystal, Phys. Rev. B 87, 174427 (2013).
  * Laurell and Fiete [2017] P. Laurell and G. A. Fiete, Topological Magnon Bands and Unconventional Superconductivity in Pyrochlore Iridate Thin Films, Phys. Rev. Lett. 118, 177201 (2017).
  * Flynn  _et al._ [2020a] V. P. Flynn, E. Cobanera, and L. Viola, Deconstructing effective non-Hermitian dynamics in quadratic bosonic Hamiltonians, New J. Phys. 22, 083004 (2020a).
  * Barnett [2013] R. Barnett, Edge-state instabilities of bosons in a topological band, Phys. Rev. A 88, 063631 (2013).
  * Schulz-Baldes [2014] H. Schulz-Baldes, Signature and spectral flow of $J$-unitary $\mathbb{S}^{1}$-Fredholm operators, Integr. Equ. Oper. Theory 78, 323 (2014).
  * Yakubovich and Starzhinskiĭ [1975] V. A. Yakubovich and V. M. Starzhinskiĭ,  _Linear Differential Equations with Periodic Coefficients_ (New York: Wiley, 1975).
  * Galilo  _et al._ [2015] B. Galilo, D. K. K. Lee, and R. Barnett, Selective Population of Edge States in a 2d Topological Band System, Phys. Rev. Lett. 115, 245302 (2015).
  * Peano  _et al._ [2016] V. Peano, M. Houde, C. Brendel, F. Marquardt, and A. A. Clerk, Topological phase transitions and chiral inelastic transport induced by the squeezing of light, Nat. Comm. 7, 1 (2016).
  * Flynn  _et al._ [2021] V. P. Flynn, E. Cobanera, and L. Viola, Topology by Dissipation: Majorana Bosons in Metastable Quadratic Markovian Dynamics, Phys. Rev. Lett. 127, 245701 (2021).
  * Diehl  _et al._ [2011] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Topology by dissipation in atomic quantum wires, Nat. Phys. 7, 971 (2011).
  * Bardyn  _et al._ [2013] C.-E. Bardyn, M. A. Baranov, C. V. Kraus, E. Rico, A. İmamoglu, P. Zoller, and S. Diehl, Topology by dissipation, New J. Phys. 15, 085001 (2013).
  * Teretenkov [2019] A. E. Teretenkov, Irreversible quantum evolution with quadratic generator: Review, Infin. Dimens. Anal. Quantum Probab. Relat. Top. 22, 1930001 (2019).
  * Barthel and Zhang [2022a] T. Barthel and Y. Zhang, Solving quasi-free and quadratic Lindblad master equations for open fermionic and bosonic systems, J. Stat. Mech. 2022, 113101 (2022a).
  * Poletti [2022] D. Poletti, Characterization of Gaussian quantum Markov semigroups, Infin. Dimens. Anal. Quantum Probab. Relat. Top. 25, 2250014 (2022).
  * Trefethen and Embree [2005] L. N. Trefethen and M. Embree, in  _Spectra and pseudospectra: the behavior of nonormal matrices and operators_ (Princeton University Press, 2005).
  * Böttcher and Silbermann [2012] A. Böttcher and B. Silbermann, in  _Introduction to Large Truncated Toeplitz Matrices_ (Springer Science & Business Media, 2012).
  * McDonald  _et al._ [2018] A. McDonald, T. Pereg-Barnea, and A. A. Clerk, Phase-dependent chiral transport and effective non-Hermitian dynamics in a bosonic Kitaev-Majorana chain, Phys. Rev. X 8, 041031 (2018).
  * Macieszczak  _et al._ [2016] K. Macieszczak, M. Guţă, I. Lesanovsky, and J. P. Garrahan, Towards a theory of metastability in open quantum dynamics, Phys. Rev. Lett. 116, 240404 (2016).
  * Macieszczak  _et al._ [2021] K. Macieszczak, D. C. Rose, I. Lesanovsky, and J. P. Garrahan, Theory of classical metastability in open quantum systems, Phys. Rev. Res. 3, 033047 (2021).
  * Okuma and Sato [2021] N. Okuma and M. Sato, Quantum anomaly, non-Hermitian skin effects, and entanglement entropy in open systems, Phys. Rev. B 103, 085428 (2021).
  * Haga  _et al._ [2021] T. Haga, M. Nakagawa, R. Hamazaki, and M. Ueda, Liouvillian Skin Effect: Slowing Down of Relaxation Processes without Gap Closing, Phys. Rev. Lett. 127, 070402 (2021).
  * Mori and Shirai [2020] T. Mori and T. Shirai, Resolving a Discrepancy between Liouvillian Gap and Relaxation Time in Boundary-Dissipated Quantum Many-Body Systems, Phys. Rev. Lett. 125, 230604 (2020).
  * Bensa and Z̆nidaric̆ [2021] J. Bensa and M. Z̆nidaric̆, Fastest Local Entanglement Scrambler, Multistage Thermalization, and a Non-Hermitian Phantom, Phys. Rev. X 11, 031019 (2021).
  * Z̆nidaric̆ [2022] M. Z̆nidaric̆, Solvable non-Hermitian skin effect in many-body unitary dynamics, Phys. Rev. Res. 4, 033041 (2022).
  * Bensa and Z̆nidaric̆ [2023] J. Bensa and M. Z̆nidaric̆, Purity decay rate in random circuits with different configurations of gates, Phys. Rev. A 107, 022604 (2023).
  * [41] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, “Topological Phases of Non-Hermitian systems,” Phys. Rev. X 8, 031079 (2018).
  * Porras and Fernández-Lorenzo [2019] D. Porras and S. Fernández-Lorenzo, Topological amplification in photonic lattices, Phys. Rev. Lett. 122, 143901 (2019).
  * Ramos  _et al._ [2021] T. Ramos, J. J. García-Ripoll, and D. Porras, Topological input-output theory for directional amplification, Phys. Rev. A 103, 033513 (2021).
  * Wanjura  _et al._ [2020] C. C. Wanjura, M. Brunelli, and A. Nunnenkamp, Topological framework for directional amplification in driven-dissipative cavity arrays, Nat. Commun. 11, 3149 (2020).
  * Wanjura  _et al._ [2021] C. C. Wanjura, M. Brunelli, and A. Nunnenkamp, Correspondence between Non-Hermitian Topology and Directional Amplification in the Presence of Disorder, Phys. Rev. Lett. 127, 213601 (2021).
  * Brunelli  _et al._ [2022] M. Brunelli, C. C. Wanjura, and A. Nunnenkamp, Restoration of the non-Hermitian bulk-boundary correspondence via topological amplification, arXiv:2207.12427v2 [quant-ph] (2022).
  * Mostafazadeh [2002] A. Mostafazadeh, Pseudo-Hermiticity versus PT-symmetry I-III, J. Math. Phys. 43, 205; 2814; 3944 (2002).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * Alicki and Lendi [1987] R. Alicki and K. Lendi,  _Quantum Dynamical Semigroups and Applications_ (Springer-Verlag, Berlin, 1987).
  * Ticozzi and Viola [2008] F. Ticozzi and L. Viola, Quantum Markovian subsystems: Invariance, attractivity, and control, IEEE Trans. Aut. Control 53, 2048 (2008).
  * Gough and James [2009] J. Gough and M. R. James, The series product and its application to quantum feedforward and feedback networks, IEEE Trans. Aut. Control 54, 2530 (2009).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of $N$-level systems, J. Math. Phys. 17, 821 (1976).
  * Schirmer and Wang [2010] S. G. Schirmer and X. Wang, Stabilizing open quantum systems by Markovian reservoir engineering, Phys. Rev. A 81, 062306 (2010).
  * Baumgartner and Narnhofer [2008a] B. Baumgartner and H. Narnhofer, Analysis of quantum semigroups with GKS-Lindblad generators: II. General, J. Phys. A 41, 395303 (2008a).
  * Temme  _et al._ [2014] K. Temme, F. Pastawski, and M. J. Kastoryano, Hypercontractivity of quasi-free quantum semigroups, J. Phys. A 47, 405303 (2014).
  * Vernier [2020] E. Vernier, Mixing times and cutoffs in open quadratic fermionic systems, SciPost Phys. 9, 049 (2020).
  * Klausen and Werner [2022] F. R. Klausen and A. H. Werner, Spectra of Lindbladians on the infinite line: From non-Hermitian to full evolution via tridiagonal Laurent matrices, arXiv:2206.09879 [math-ph] (2022).
  * Gardiner and Zoller [2000] C. W. Gardiner and P. Zoller,  _Quantum Noise_ (Springer, Heidelberg, 2000).
  * Arkhipov  _et al._ [2020a] I. I. Arkhipov, A. Miranowicz, F. Minganti, and F. Nori, Quantum and semiclassical exceptional points of a linear system of coupled cavities with losses and gain within the Scully-Lamb laser theory, Phys. Rev. A 101, 013812 (2020a).
  * Genoni  _et al._ [2016] M. G. Genoni, L. Lami, and A. Serafini, Conditional and unconditional Gaussian quantum dynamics, Contemp. Phys. 57, 331 (2016).
  * Prosen and Seligman [2010] T. Prosen and T. H. Seligman, Quantization over boson operator spaces, J. Phys. A 43, 392004 (2010).
  * Krüger [2006] O. Krüger,  _Quantum Information Theory with Gaussian Systems_ , Ph.D. thesis, Technical University of Braunschweigs (2006).
  * Berislav and Prosen [2012] B. Berislav and T. Prosen, A note on symmetry reductions of the Lindblad equation: transport in constrained open spin chains, New J. Phys. 14, 073007 (2012).
  * Albert and Jiang [2014] V. V. Albert and L. Jiang, Symmetries and conserved quantities in Lindblad master equations, Phys. Rev. A 89, 022118 (2014).
  * Cobanera  _et al._ [2017] E. Cobanera, A. Alase, G. Ortiz, and L. Viola, Exact solution of corner-modified banded block-Toeplitz eigensystems, J. Phys. A 50, 195204 (2017).
  * Yao and Wang [2018] S. Yao and Z. Wang, Edge States and Topological Invariants of Non-Hermitian Systems, Phys. Rev. Lett. 121, 086803 (2018).
  * Borgnia  _et al._ [2020] D. S. Borgnia, A. J. Kruchkov, and R.-J. Slager, Non-Hermitian Boundary Modes and Topology, Phys. Rev. Lett. 124, 056802 (2020).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Topological Origin of Non-Hermitian Skin Effects, Phys. Rev. Lett. 124, 086801 (2020).
  * Note [3] To spell out the analogy further, one arrives at the Bloch dynamical matrix following a Fourier transform. Similarly, the action of $\mathbf{X}^{\text{BIBC}}$ on plane wave modes $\vec{f_{k}}=\sum@\slimits@_{j}e^{ijk}\vec{e}_{j}$ is characterized by $x(e^{ik})$ according to $\mathbf{X}^{\text{BIBC}}\vec{f}_{k}=x(e^{ik})\vec{f}_{k}$.
  * Note [4] This follows from the property ${\mathcal{L}}^{\star}(A^{\dagger})=[{\mathcal{L}}^{\star}(A)]^{\dagger}$, for all $A$. Therefore, if $\hat{\vec{\alpha}}$ is not Hermitian to begin with, we may form the linear combination $\hat{\vec{\alpha}}+\hat{\vec{\alpha}}^{\dagger}$, which is Hermitian and still belongs to ker ${\mathcal{L}^{\star}}$, as desired.
  * Note [5] In fact, $\mathbf{G}$ is $\bm{\tau}_{3}$-pseudo-Hermitian whenever $\mathcal{F}(\mathbf{M})=0$, which ensures unitality of the dynamics, $\mathcal{L}(1_{\mathcal{F}})=0$.
  * Flynn  _et al._ [2020b] V. P. Flynn, E. Cobanera, and L. Viola, Restoring number conservation in quadratic bosonic Hamiltonians with dualities, Europhys. Lett. 131, 40006 (2020b).
  * 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).
  * Note [6] One could, alternatively, consider a chain extending infinitely to the left, which will retain all right-localized modes and lose the left-localized ones.
  * Mori [2021] T. Mori, Metastability associated with many-body explosion of eigenmode expansion coefficients, Phys. Rev. Res. 3, 043137 (2021).
  * Vasiloiu  _et al._ [2018] L. M. Vasiloiu, F. Carollo, and J. P. Garrahan, Enhancing correlation times for edge spins through dissipation, Phys. Rev. B 98, 094308 (2018).
  * Ughrelidze  _et al._ [2024] M. Ughrelidze, V. P. Flynn, E. Cobanera, and L. Viola, Upcoming paper covering the interplay between finite- and infinite-size stability in bosonic systems,  (2024).
  * Note [7] Note that the canonical correspondence assumes the technical assumptions of Theorem 2 in App. B.2 are met.
  * Flynn  _et al._ [2024] V. P. Flynn, E. Cobanera, and L. Viola, Upcoming paper covering more mathematically sophisticated instances of topological metastability,  (2024).
  * Kraus  _et al._ [2008] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Preparation of entangled states by quantum Markov processes, Phys. Rev. A 78, 042307 (2008).
  * Note [8] The fact that the dissipator has more than one diagonal representation, Eqs. (51) and (52), results from a large degeneracy in the GKLS matrix spectrum. In particular, $\mathbf{M}$ has an $N$-fold degenerate $0$ eigenvalue. Different diagonal representation of $\mathcal{D}$ are derived from different choices of bases for the $0$ eigenvalue eigenspace.
  * D’Abbruzzo and Rossini [2021] A. D’Abbruzzo and D. Rossini, Self-consistent microscopic derivation of Markovian master equations for open quadratic quantum systems, Phys. Rev. A 103, 052209 (2021).
  * Fujii [2012] K. Fujii, Quantum damped harmonic oscillator, in  _Advances in Quantum Mechanics_ , edited by P. Bracken (InTech, 2012) Chap. 7, arXiv:1209.1437 [quant-ph].
  * Korsch [2019] H. J. Korsch, Lindblad dynamics of the damped and forced quantum harmonic oscillator, arXiv:1908.01187v2 (quant-ph) (2019).
  * Albert  _et al._ [2018] V. V. Albert, K. Noh, K. Duivenvoorden, D. J. Young, R. T. Brierley, P. Reinhold, C. Vuillot, L. Li, C. Shen, S. M. Girvin, B. M. Terhal, and L. Jiang, Performance and structure of single-mode bosonic codes, Phys. Rev. A 97, 032346 (2018).
  * Guillaud  _et al._ [2023] J. Guillaud, J. Cohen, and M. Mirrahimi, Quantum computation with cat qubits, SciPost Phys. Lect. Notes 72 (2023).
  * Hatano and Nelson [1996] N. Hatano and D. R. Nelson, Localization Transitions in Non-Hermitian Quantum Mechanics, Phys. Rev. Lett. 77, 570 (1996).
  * McDonald  _et al._ [2022] A. McDonald, R. Hanai, and A. A. Clerk, Non-equilibrium stationary states of quantum non-Hermitian lattice models, Phys. Rev. B 105, 064302 (2022).
  * Livi and Politi [2017] R. Livi and P. Politi,  _Nonequilibrium Statistical Physics: A Modern Perspective_ (Cambridge University Press, Cambridge, 2017).
  * Sciolla  _et al._ [2015] B. Sciolla, D. Poletti, and C. Kollath, Two-time correlations probing the dynamics of dissipative many-body quantum systems: Aging and fast relaxation, Phys. Rev. Lett. 114, 170401 (2015).
  * Brown  _et al._ [1997] G. Brown, P. A. Rikvold, M. Sutton, and M. Grant, Speckle from phase-ordering systems, Phys. Rev. E 56, 6601 (1997).
  * Freericks and Kemper [2021] J. K. Freericks and A. F. Kemper, What do the two times in two-time correlation functions mean for interpreting tr-ARPES?, J. El. Spectr. Rel. Phenom. 251, 147104 (2021).
  * Walls and Milburn [2007] D. F. Walls and G. J. Milburn,  _Quantum Optics_ (Springer, Berlin, Heidelberg,, 2007).
  * Arkhipov  _et al._ [2020b] I. I. Arkhipov, A. Miranowicz, F. Minganti, and F. Nori, Liouvillian exceptional points of any order in dissipative linear bosonic systems: Coherence functions and switching between $\mathcal{PT}$ and anti-$\mathcal{PT}$ symmetries, Phys. Rev. A 102, 033715 (2020b).
  * Koch  _et al._ [2010] J. Koch, A. A. Houck, K. Le Hur, and S. M. Girvin, Time-reversal-symmetry breaking in circuit-QED-based photon lattices, Phys. Rev. A 82, 043811 (2010).
  * Kamal  _et al._ [2011] A. Kamal, J. Clarke, and M. H. Devoret, Noiseless non-reciprocity in a parametric active device, Nat. Phys. 7, 311 (2011).
  * Abdo  _et al._ [2013a] B. Abdo, K. Sliwa, L. Frunzio, and M. Devoret, Directional Amplification with a Josephson Circuit, Phys. Rev. X 3, 031001 (2013a).
  * Sliwa  _et al._ [2015] K. M. Sliwa, M. Hatridge, A. Narla, S. Shankar, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Reconfigurable Josephson Circulator/Directional Amplifier, Phys. Rev. X 5, 041020 (2015).
  * Metelmann and Clerk [2015] A. Metelmann and A. A. Clerk, Nonreciprocal Photon Transmission and Amplification via Reservoir Engineering, Phys. Rev. X 5, 021025 (2015).
  * Anderson  _et al._ [2016] B. M. Anderson, R. Ma, C. Owens, D. I. Schuster, and J. Simon, Engineering Topological Many-Body Materials in Microwave Cavity Arrays, Phys. Rev. X 6, 041043 (2016).
  * Roushan et al. [2017] P. Roushan et al., Chiral ground-state currents of interacting photons in a synthetic magnetic field, Nat. Phys. 13, 146 (2017).
  * Fitzpatrick  _et al._ [2017] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Observation of a Dissipative Phase Transition in a One-Dimensional Circuit QED Lattice, Phys. Rev. X 7, 011016 (2017).
  * Owens  _et al._ [2018] C. Owens, A. LaChapelle, B. Saxberg, B. M. Anderson, R. Ma, J. Simon, and D. I. Schuster, Quarter-flux Hofstadter lattice in a qubit-compatible microwave cavity array, Phys. Rev. A 97, 013818 (2018).
  * Zhang  _et al._ [2019] Y. Zhang, B. J. Lester, Y. Y. Gao, L. Jiang, R. J. Schoelkopf, and S. M. Girvin, Engineering bilinear mode coupling in circuit QED: Theory and experiment, Phys. Rev. A 99, 012314 (2019).
  * Kim et al. [2021] E. Kim et al., Quantum Electrodynamics in a Topological Waveguide, Phys. Rev. X 11, 011015 (2021).
  * Gong  _et al._ [2021] Y. Gong, L. Guo, S. Wong, A. J. Bennett, and S. Soon Oh, Tailoring topological edge states with photonic crystal nanobeam cavities, Sci. Rep. 11, 1055 (2021).
  * Zhao  _et al._ [2018] H. Zhao, P. Miao, M. H. Teimourpour, S. Malzard, R. El-Ganainy, H. Schomerus, and L. Feng, Topological hybrid silicon microlasers, Nat. Comm. 9, 981 (2018).
  * Parto et al. [2018] M. Parto et al., Edge-Mode Lasing in 1D Topological Active Arrays, Phys. Rev. Lett. 120, 113901 (2018).
  * Malz  _et al._ [2018] D. Malz, L. D. Tóth, N. R. Bernier, A. K. Feofanov, T. J. Kippenberg, and A. Nunnenkamp, Quantum-Limited Directional Amplifiers with Optomechanics, Phys. Rev. Lett. 120, 023601 (2018).
  * Youssefi  _et al._ [2022] A. Youssefi, S. Kono, A. Bancora, M. Chegnizadeh, J. Pan, T. Vovk, and T. J. Kippenberg, Topological lattices realized in superconducting circuit optomechanics, Nature 612, 666 (2022).
  * Abdo  _et al._ [2013b] B. Abdo, A. Kamal, and M. Devoret, Nondegenerate three-wave mixing with the Josephson ring modulator, Phys. Rev. B 87, 014508 (2013b).
  * Frattini  _et al._ [2017] N. E. Frattini, U. Vool, S. Shankar, A. Narla, K. M. Sliwa, and M. H. Devoret, 3-wave mixing Josephson dipole element, Appl. Phys. Lett. 110, 222603 (2017).
  * Sivak  _et al._ [2019] V. V. Sivak, N. E. Frattini, V. R. Joshi, A. Lingenfelter, S. Shankar, and M. H. Devoret, Kerr-Free Three-Wave Mixing in Superconducting Quantum Circuits, Phys. Rev. Appl. 11, 054060 (2019).
  * Pocklington  _et al._ [2023] A. Pocklington, Y. X. Wang, and A. A. Clerk, Dissipative Pairing Interactions: Quantum Instabilities, Topological Light, and Volume-Law Entanglement, Phys. Rev. Lett. 130, 123602 (2023).
  * Prosen [2008] T. Prosen, Third quantization: a general method to solve master equations for quadratic open Fermi systems, New J. Phys. 10, 043026 (2008).
  * Zhang and Barthel [2022b] Y. Zhang and T. Barthel, Criticality and Phase Classification for Quadratic Open Quantum Many-Body Systems, Phys. Rev. Lett. 129, 120401 (2022b).
  * Gottesman  _et al._ [2001] D. Gottesman, A. Kitaev, and J. Preskill, Encoding a qubit in an oscillator, Phys. Rev. A 64, 012310 (2001).
  * Flühmann  _et al._ [2019] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Encoding a qubit in a trapped-ion mechanical oscillator, Nature 566, 513 (2019).
  * Campagne-Ibarcq  _et al._ [2020] P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zayls-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Quantum error correction of a qubit encoded in grid states of an oscillator, Nature 584, 368 (2020).