††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Topological zero modes and edge symmetries of metastable Markovian bosonic
systems

Vincent P. Flynn  Department of Physics and Astronomy, Dartmouth College, 6127
Wilder Laboratory, Hanover, New Hampshire 03755, USA Emilio Cobanera
Department of Mathematics and Physics, SUNY Polytechnic Institute, 100 Seymour
Rd, Utica, NY 13502, USA Department of Physics and Astronomy, Dartmouth
College, 6127 Wilder Laboratory, Hanover, New Hampshire 03755, USA Lorenza
Viola  Department of Physics and Astronomy, Dartmouth College, 6127 Wilder
Laboratory, Hanover, New Hampshire 03755, USA

###### Abstract

Tight bosonic analogs of free-fermionic symmetry-protected topological phases,
and their associated edge-localized excitations, have long evaded the grasp of
condensed-matter and AMO physics. In this work, building on our initial
exploration [Phys. Rev. Lett. 127, 245701 (2021)], we identify a broad class
of quadratic bosonic systems subject to Markovian dissipation that realize
tight bosonic analogs of the Majorana and Dirac edge modes characteristic of
topological superconductors and insulators, respectively. To this end, we
establish a general framework for topological metastability for these systems,
by leveraging pseudospectral theory as the appropriate mathematical tool for
capturing the non-normality of the Lindbladian generator. The resulting
dynamical paradigm, which is characterized by both a sharp separation between
transient and asymptotic dynamics and a non-trivial topological invariant, is
shown to host edge-localized modes, which we dub Majorana and Dirac bosons.
Generically, such modes consist of one conserved mode and a canonically
conjugate generator of an approximate phase-space translation symmetry of the
dynamics. The general theory is exemplified through several representative
models exhibiting the full range of exotic boundary physics that topologically
metastable systems can engender. In particular, we explore the extent to which
Noether’s theorem is violated in this dissipative setting and the way in which
certain symmetries can non-trivially modify the edge modes. Notably, we also
demonstrate the possibility of anomalous parity dynamics for a bosonic cat
state prepared in a topologically metastable system, whereby an equal
distribution between even and odd parity sectors is sustained over a long
transient. For both Majorana and Dirac bosons, observable multitime signatures
in the form of anomalously long-lived quantum correlations and divergent zero-
frequency power spectral peaks are proposed and discussed in detail. Our
results point to a new paradigm for symmetry-protected topological physics in
free bosons, embedded deeply in the long-lived transient regimes of metastable
dynamics.

##  I Introduction

###  I.1 Context and motivation

Indistinguishable quantum particles come in two flavors: fermions and bosons.
While the distinction is kinematical and, as such, unrelated to any
Hamiltonian specification, it can be explained most clearly when the particles
are independent, or “free”. Systems of free fermions (bosons) are described by
Hamiltonians that are quadratic in their respective canonical fermionic
(bosonic) operators, and have long played a paradigmatic role as tractable –
either genuinely non-interacting or mean-field – models for both equilibrium
and non-equilibrium many-body physics [1]. For a quadratic fermionic
Hamiltonian (QFH), there is always a state of lowest energy, the ground state,
which captures the statistical behavior of the system in equilibrium at (and
close to) zero temperature. A quantum phase transition is a phase transition
at zero temperature that occurs as some parameter of the Hamiltonian is
varied. Generically, the phases of free-fermion quantum matter are gapped,
display no local order parameter, and the energy gap closes at a phase
transition. Since there is no local order parameter, there is no general
Landau theory of the quantum phases of QFHs. Rather, the general theory of
such phases is based on a different set of notions: that of protecting global
symmetries, space dimension, and topological invariants. The phases of free
fermions are examples of symmetry-protected topological (SPT) phases of
quantum matter [2].

In the absence of local order parameters, how can one tell apart the different
SPT phases of free fermions? A compelling answer is provided by the bulk-
boundary correspondence. This powerful principle states that the topological
invariant that characterizes an SPT phase also mandates the emergence of
robust zero-energy boundary-localized modes [3, 4]. These zero modes (ZMs) are
regarded as the main experimental manifestation of the underlying SPT phase.
For example, the integer quantum Hall regimes in two (spatial) dimensions are
SPT phases. In this case, the protecting symmetry is particle number and the
topological invariant is the Chern number of the occupied single-particle
energy bands. A measurement of the quantized Hall conductance probes directly
the associated chiral surface modes. In one dimension, the Su-Schrieffer-
Heeger model of polyacetylene also displays topologically-mandated edge modes
[5]. The protecting symmetries are particle number, spin rotations, (spinful)
time reversal, and a many-body particle-hole symmetry that exchanges fermionic
creation and annihilation operators. Likewise, superconductors can exist in a
variety of SPT phases. While particle number cannot be one of the protecting
symmetries, the superconducting classes are protected by a combination of spin
symmetry, time reversal, and many-body particle-hole. The celebrated Majorana
chain of Kitaev provides a paradigmatic example for $p$-wave topological
superconductivity and can display edge ZMs [6], which are protected against
(weak) perturbations that do not change the symmetry of the model [4].
Altogether, SPT phases of free fermions are distinguished by the following key
features: (i) The translationally invariant (bulk) system is gapped; (ii) The
system displays certain combinations of protecting many-body symmetries; (iii)
The ground state has an associated topological number, which can only change
across a quantum phase transition as long the protecting symmetries are
preserved; (iv) When the topological number is non-trivial and the system is
terminated by imposing open (‘hard-wall’) boundary conditions, the ground
energy level is degenerate due to low-energy surface quasi-particles
111Importantly, 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..

Having said that, if one understands the bulk-boundary correspondence as the
relationship between a topological invariant and quasiparticle surface or edge
modes, then there is nothing particularly “fermionic,” or symmetry-protected,
about it. Rather, it is a general property of the Helmoltz wave equation in a
structured medium. The point was forcefully made in Ref. 8, where a bulk-
boundary correspondence for photonic crystals was identified and
experimentally confirmed only one year later [9]. These advances launched in
earnest the new field of topological photonics, which has since witnessed a
dramatic development [10]. From a fundamental perspective, one may ask whether
topological photonics truly represents the mirror image of topological
electronics. If the answer is in the affirmative, then quantum statistics
would appear to have very little to do with topological physics. More
generally, to what extent can systems of free bosons exhibit SPT physics
analogous to the above? The answer is complicated, and indeed depends on what
concepts are emphasized.

Consider, for example, the concept of a “quantum phase.” Unlike a QFH, a
quadratic bosonic Hamiltonian (QBH) may not have a ground state: take a simple
two-mode QBH like $H=\omega(a_{1}^{\dagger}a_{1}-a_{2}^{\dagger}a_{2})$, with
frequency $\omega>0$. For such a QBH (which can arise, e.g., in cavity QED
systems [11]), the energy eigenvalues $E_{n_{1},n_{2}}=\omega(n_{1}-n_{2})$,
with $n_{1,2}$ non-negative integers, are unbounded in both directions without
external constraints placed on the total particle number. We say such
Hamiltonians are thermodynamically unstable, in the sense that no well-defined
Gibbs state $\rho\sim e^{-\beta H}$ exists. Since the ground state plays such
a crucial role in fermionic SPT physics, the most conservative extension into
the bosonic realm is to consider only those QBHs that are thermodynamically
stable. As it turns out, this subclass fails to exhibit any of the
characteristics of topological free-fermionic matter. The situation is neatly
captured by three no-go theorems [12], which we summarize as follows. If $H$
is a thermodynamically stable, gapped, translationally-invariant system, then:
(i) $H$ can be adiabatically deformed into any other QBH within the same class
without closing the gap or breaking symmetries; (ii) no edge-localized ZMs
emerge when the system is terminated at a boundary; and (iii) the ground state
of $H$ has always even bosonic parity. 222No-go (iii) does not require
translation invariance. In other words, there are no SPT phases of free-boson
quantum matter. This is a puzzling conclusion, because the role of topology
for QFHs is precisely to classify their gapped quantum phases. Fundamentally
speaking, mean-field bosonic matter has very little in common with its
fermionic counterparts from the point of view of quantum many-body physics.
Nonetheless, thermodynamically stable bosonic systems, such as certain
photonic, magnonic, and phononic crystals, may exhibit topologically mandated
edge modes at higher, non-zero energies [14, 15, 16, 17, 10]. These
topological features, however, are completely disconnected from the low-
energy, low-temperature physics.

If one is willing to accept the loss of a many-body ground state and remove
the constraint of thermodynamic stability, a gap condition may still be
imposed by requiring that the quasiparticle energies (which are now
necessarily not strictly positive) are gapped. That is, one can require that
the quasiparticle energy bands are bounded away from zero. Thermodynamically
unstable, “gapped” QBHs can then display genuinely topological ZMs under open
boundary conditions. At face value, this (rather significant) compromise
allows for the possibility of obtaining tight bosonic analogs of the ZMs
characteristic of fermionic SPT phases. In Ref. 18, we presented a QBH model
hosting Hermitian edge modes that are canonically conjugate (once properly
normalized), commute with the Hamiltonian in the infinite-size limit, are
determined by a winding number, and exist in a well-defined (complex) spectral
gap. Nonetheless, we are confronted with another major issue: due to the
intrinsic non-Hermiticity of the dynamical matrix that govern the system’s
behavior, these bosonic ZMs are intrinsically unstable, in a dynamical sense
[19, 20, 18]. Specifically, tools from Krein stability theory [21] reveal that
any QBH hosting bosonic ZMs is either (i) dynamically unstable (characterized
by unbounded evolution of observables), or (ii) able to be destabilized by
arbitrarily small perturbations. While the instability of these modes may have
interesting implications in its own [22, 23], and the modes themselves can be
thought to as “shadows” of conventional Majorana ZMs [18], again we find
ourselves far removed from the world of topological many-body physics. What
(if anything) can be done to bring us closer to a bosonic analogue of the rich
SPT physics enjoyed by free-fermionic matter?

Our fundamental realization is the need to let go of a more subtle assumption
of the fermionic paradigm: unitarity. Specifically, in this work we will
provide extensive evidence that open systems of free bosons can display SPT-
like many-body physics that comes as close as possible to the SPT physics of
QFHs [24]. In this sense, our work fully embraces the idea of “topology by
dissipation” introduced in Ref. 25, 26. However, there is a crucial
difference. For fermions, as considered in these works, dissipation is a twist
that one can add to a very well-developed theory for Hamiltonian systems: with
dissipation or not, free fermions can be topological in the many-body sense.
As we have emphasized, bosons do not seem amenable to that. Topology by
dissipation could then very well be the only hope to bring about bosonic SPT
physics without strong interactions.

The program we undertake begins by appropriately adapting the key ingredients
of free fermionic SPTs to the open bosonic setting. Focusing on the simplest
case of Markovian dissipation in one spatial dimension, we retain the non-
interacting property by restricting to the class of dynamics described by
“quasi-free,” or quadratic (Gaussian) Lindblad generators [27, 28, 29]. Under
certain stability assumptions, the notion of a ground state naturally maps to
that of a steady state (SS), while the many-body gap condition maps to one
placed on the Lindblad, or spectral gap. With these identifications made, we
answer the question: to what extent can SPT-like physics emerge in quadratic
bosonic Lindbladians (QBLs) possessing a unique SS and a finite spectral gap?
Remarkably, signatures of SPT physics are found to emerge in a newly
identified dynamical phase that we deem topologically metastable. Topological
metastability is, in turn, a specific instance of the more general dynamical
metastabilty, a phenomenon in which the dynamical stability of the system
changes abruptly in the infinite-size limit. As we have pointed out in Ref.
24, dynamical metastability may or may not be topological in nature.
Topological metastability arises precisely from requiring non-trivial bulk
topology, on top of dynamical metastability. The key features of topologically
metastable systems may be summarized as follows [24]:

  * •

A unique SS and a finite spectral gap are maintained for all finite system
sizes. In particular, dynamical stability is present for all finite system
sizes.

  * •

Tight bosonic analogues of Majorana fermions, which we deem Majorana bosons
(MBs), emerge localized on opposite ends of the chain. They consist of an
approximate ZM and an approximate-symmetry generator (SG) and are canonically
conjugate, despite macroscopic spatial separation.

  * •

A manifold of degenerate quasi-steady states manifest in the finite-size
chains. Physically, they are displacements (generated by the SG) of the unique
steady state.

  * •

Both the ZM, and the quasi-steady states persist in a transient dynamical
regime whose duration diverges with system size. Further, their existence
elicits divergent zero-frequency peaks in certain power spectra.

Working within the aforementioned identification scheme for linking closed and
dissipative many-body phenomena, the first three of these features (save for
the split roles of SGs and ZMs) closely resemble the generic features of
topologically non-trivial free fermionic systems subject to open boundary
conditions. For instance, in addition to the mathematical similarities between
the relevant edge modes, the joint presence of a steady state and a manifold
of quasi-steady states of the third, are reminiscent of the (nearly)
degenerate ground states found in, e.g., the Kitaev chain. A less-obvious
similarity arises by analogizing the long, but finite lifetimes of the ZM and
the quasi-steady states with the (generically) exponentially small, but non-
zero, energies associated to Majorana fermions in finite systems. Both
phenomena are finite-size effects in nature, and may be related to one another
by identifying the exponentially small energies in the fermionic case to the
exponentially small decay rates in the bosonic case.

Perhaps the most dramatic conceptual difference between topological fermions
and topologically metastable bosons arises when one considers their bulk
physics. As previously noted, for free fermions topological transitions are
inextricably tied to a bulk quantum phase transition (without spontaneous
symmetry-breaking). In particular, topological phases retain a unique (bulk)
ground state. In sharp contrast, topologically metastable bosonic systems are
necessarily bulk-unstable, despite remaining stable for all finite sizes.
Consequently, they lack a bulk steady state altogether. This leads us to
conclude that fair comparisons can only be made by considering finite systems
with open boundary conditions. We summarize the main comparisons for systems
lacking total number symmetry in Table 1.

Topologically non-trivial QFH | Topologically metastable QBL
---|---
Edge modes (finite $N$)
Approximate ZM pairs | Approximate ZM and SG pairs
Exponentially small energies | Exponentially small decay rates
Canonically conjugate | Canonically conjugate
Hermitian | Hermitian
Robust to weak perturbations | Robust to weak perturbations
“Low-energy” boundary physics (finite $N$)
(Nearly) degenerate ground states | Manifold of quasi-steady states
“Low-energy” boundary physics ($N\to\infty$)
Degenerate ground states | No steady states, unstable
Bulk physics ($N\to\infty$)
Non-zero bulk invariant | Non-zero bulk invariant
Unique ground state | No steady state, unstable
Table 1: Comparisons between topologically non-trivial QFHs and topologically
metastable QBLs, both in the absence of total number symmetry. The first five
rows correspond to properties of the edge modes that arise in finite systems
with open boundary conditions. The next two rows correspond to “low-energy”
(ground/steady state) features for a system with 2 ($N$ finite) or 1
($N\to\infty$) boundaries. The last two rows correspond to bulk (i.e.,
boundary-less) properties.

###  I.2 Outline and summary of main results

This paper solidifies and expands greatly on the above core ideas. From a
technical standpoint, a main tool for investigation is pseudospectral theory
[30, 31], which provides the appropriate mathematical framework for
convergence and stability analysis in the presence of non-normal dynamical
generators. In terms of general results, we identify four major ones:

  * •

We provide a self-contained proof that there exists a canonical correspondence
between ZMs and a class of linear symmetry generators in QBLs, despite the
explicit breakdown of Noether’s theorem in the dissipative Markovian setting.
That is, a partial restoration of Noether’s theorem is afforded by QBLs.

  * •

We introduce and apply two design principles for engineering QBLs with desired
features. The first constitutes an explicit mapping from topologically non-
trivial QFHs to QBLs that host MBs. The second provides a reservoir-
engineering scheme for designing a QBL that relaxes to the quasiparticle
vacuum of a given QBH.

  * •

We uncover the existence of tight bosonic analogues of the Dirac-fermion edge
modes characteristic of topological insulators in QBLs that possess total-
number symmetry. These bosonic modes, which we call Dirac bosons, are tightly
connected to two pairs of MBs and, as such, possess many of the same notable
properties.

  * •

We demonstrate a strong connection between topological metastability and the
existence of arbitrarily long-lived two-time quantum correlation functions.
Specifically, the macroscopically separated MBs show such long-lived
correlations.

To illustrate the consequences of these general results, we provide several
representative models – beyond the dissipative bosonic Kitaev chain we
introduced in Ref. 24 and we also revisit here. The first new model is a
purely dissipative ($H=0$) chain that descends from the fermionic Kitaev chain
and most closely mirrors the purely dissipative setting considered for
quadratic fermionic Lindbladians (QFLs) [25, 26]. This model demonstrates that
a coherent, Hamiltonian contribution is not needed for topological
metastability. It also possesses a new type of MBs, whereby both members of
the pairs are approximate ZMs and generate an approximate symmetry. We call
these non-split MBs. The second new model consists of the bosonic Kitaev chain
Hamiltonian [32, 18] with a specially engineered dissipator that ensures
relaxation to a pure SS. Remarkably, while the dissipator explicitly breaks
translation invariance, a restricted translation invariance is retained,
allowing for relevant spectral properties to be computed. Purity of the SS
grants us analytical access to the exact dynamics of the quasi-steady states;
in particular, these may exhibit highly non-trivial bosonic parity dynamics,
with transient odd bosonic parity, in an appropriate sense. The third new
model is a number-symmetric chain that possesses the aforementioned Dirac
bosons. Being this the first example exhibiting such modes, we study their
algebraic and dynamical properties in great detail.

In more detail, the content is organized as follows. In Sec. II, we first
establish relevant notation and second-quantization formalism in the closed-
system setting of quadratic fermionic and bosonic Hamiltonians; we then move
to Markovian quantum systems, and provide the necessary background about the
Lindblad formalism and quadratic Lindbladians, with emphasis on symmetry
properties; lastly, we summarize basic notions and results about
pseudospectra. In Sec. III, we present three main foundational results on
QBLs: the correspondence between certain conserved quantities and symmetry
generators, and two design protocols for generating QBLs with certain non-
trivial features to be utilized in later examples, as mentioned above. In Sec.
IV, we introduce the notion of dynamical metastability in one-dimensional (1D)
bulk-translationally invariant QBLs. We then focus on our main goal of
exhibiting topologically nontrivial metastable dynamics, by synthesizing
dynamical metastability and non-trivial bulk topology. We do so by separately
discussing the case where number symmetry is broken due to the presence of
“bosonic pairing” of either Hamiltonian or dissipative nature (Sec. V), and
the case where number symmetry is unbroken (Sec. VI). In particular, we
present general results of, and several example models displaying Majorana and
Dirac bosons, respectively. Observable signatures are explored in detail in
Sec. VII, whereby multitime correlation functions and their associated power
spectra are computed and analyzed in detail for several of our models.

A number of additional results are included in separate appendixes. In
particular, SPT phases of free fermions are further discussed in App. A,
including a detailed analysis of the edge ZMs emerging in the topologically
non-trivial regimes of the two paradigmatic Su-Schrieffer-Heeger and Kitaev
chains. In App. B and App. C, we collect a number of proofs and explicit
calculations supporting claims made in the main text. Finally, App. D presents
a fourth new model which, unlike all models discussed in the main text,
support two MBs of the same type (e.g., ZMs) on one edge in the topological
regime and can thus be of independent interest.

###  I.3 Relation to existing work

With the program and our key results laid out, we wish to further place our
contributions in the context of existing work. Firstly, concepts of
metastability motivated by classical statistical physics have been extended to
Markovian systems, and studied in great detail [33, 34]. In essence, this form
of metastability is characterized by multi-step relaxation, mandated by the
presence of large gaps in the Lindblad spectra. While our dynamically
metastable systems possess no such spectral gaps, it turns out they do exhibit
pseudospectral gaps. That is, the nontrivial pseudospectral (determined, as we
will see, by the spectra under semi-infinite boundary conditions) remains
gapped away from the exact spectrum in such a way to mandate anomalous
transient dynamics (see Fig. Sec. IV below). In fact, pseudospectra has been
conjectured [35] to also play a role in the recent discoveries of anomalous
relaxation dynamics in Markovian systems exhibiting a non-Hermitian skin-
effect [36, 37]. They have further been successfully applied to study
anomalous dynamics in random quantum circuits [38, 39, 40] and explicitly non-
Hermitian Hamiltonians [41, 35].

A second branch of research which, while not motivated by the many-body
physics of topological free fermions, shares notable points of contact with
our analysis, pertains to topological amplification [42, 43, 44, 45, 46]. The
standard approach to topological amplifiers employs input-output theory.
Central to the input-output treatment is the susceptibility matrix, or Green’s
function, say, $\bm{\chi}(\omega)$, which connects incoming fields at
frequency $\omega$ to outgoing fields at frequency $\omega$. For photonic
systems with $N$ modes, whose coherent and incoherent dynamics can be cast in
the form of a Lindblad master equation, this susceptibility matrix takes the
form $\bm{\chi}(\omega)=i(\omega\mathds{1}_{2N}-\mathbf{G})^{-1}$, with
$\mathbf{G}$ being the dynamical matrix we will extensively discuss in Sec.
II.2.2. It is then known from the above works that topological amplification
can only take place if $\bm{\chi}^{-1}(\omega)$ winds around the origin, in a
suitable sense.

As it turns out, in our language this is, in the simplest case, equivalent to
nontrivial winding of a certain spectral, “rapidity band” about the point
$i\omega$. According to our pseudospectral approach, this implies that there
exist pseudonormal modes with pseudoeigenvalues $i\omega+\lambda$, with
$\lambda>0$. Such pseudonormal modes must necessarily amplify in the
transient, and hence contribute to gain in the output signal. Moreover, we see
that when these systems amplify zero-frequency ($\omega=0$) signals, they must
possess MBs and all of the associated non-trivial transient dynamics they
entail. To summarize, in our framework, general topological amplifiers are
classified as dynamically metastable, while those that amplify zero-frequency
input signals are classified as topologically metastable. We further observe
that the introduction and application of the “doubled matrix” approach in
Refs. 42, 43 can be connected naturally to pseudospectral theory (see App.
B.1).

##  II Background

###  II.1 Warm up: Quadratic bosonic and fermionic Hamiltonians

In the language of second quantization, a closed system of independent bosons
(fermions) in equilibrium is described by a time-independent QBH (QFH).
Despite the profound physical differences that exist between bosons and
fermions, these Hamiltonians are diagonalized in essentially the same way. The
key step is to solve a commutator equation involving the quadratic Hamiltonian
of interest, to find the normal modes of the system. In this section, we
summarize some machinery adapted to this task both for fermions and bosons,
following exactly the same logic and comparing the mathematical structures
that emerge. In doing so, we also introduce several concepts and notational
conventions we will use in the more general setting of open quadratic dynamics
in the next section and throughout this paper.

Let $c_{i}$, $i=1,\dots,N$, denote a set of fermionic annihilation operators
satisfying the canonical anticommutation relations,
$\\{c_{i},c_{j}^{\dagger}\\}=\delta_{ij}1_{\mathcal{F}}$ and
$\\{c_{i},c_{j}\\}=0$, with $1_{\mathcal{F}}$ the fermionic Fock space
identity. A QFH is given by

$H_{F}=\frac{1}{2}\sum_{i,j=1}^{N}\left(\mathbf{K}_{ij}c_{i}^{\dagger}c_{j}+\bm{\Delta}_{ij}c_{i}^{\dagger}c_{j}^{\dagger}+\text{H.c.}\right),$

with $\mathbf{K}=\mathbf{K}^{\dagger}$ the $N\times N$ hopping matrix and
$\bm{\Delta}^{T}=-\bm{\Delta}$ the $N\times N$ pairing matrix. Such
Hamiltonians are more compactly expressed in terms of the fermionic Nambu
array,

$\Psi=[c_{1},c_{1}^{\dagger},\ldots,c_{N},c_{N}^{\dagger}]^{T}.$

It then follows that

$H_{F}=\frac{1}{2}\Psi^{\dagger}\mathbf{H}_{F}\Psi,$

where the $2N\times 2N$ Bogoliubov-de Gennes (BdG) $\mathbf{H}_{F}$ is a block
matrix with the $(ij)$-th block given by [1]

$[\mathbf{H}]_{ij}=\begin{bmatrix}\mathbf{K}_{ij}&\bm{\Delta}_{ij}\\\
-\bm{\Delta}_{ij}^{*}&-\mathbf{K}_{ij}^{*}\end{bmatrix}.$

The stated conditions on $\mathbf{K}$ and $\bm{\Delta}$ imply that (i) the
matrix $\mathbf{H}_{F}$ Hermitian; and (ii)
$\mathbf{H}_{F}=-\bm{\tau}_{1}\mathbf{H}_{F}^{T}\bm{\tau}_{1}$ in terms of the
matrices $\bm{\tau}_{j}\equiv\mathds{1}_{N}\otimes\bm{\sigma}_{j}$, with
$\mathds{1}_{N}$ the $N\times N$ identity matrix and $\bm{\sigma}_{j}$,
$j=1,2,3$, the usual Pauli matrices.

Condition (ii) can be formalized in terms of a fermionic projector, which we
define according to

$\mathcal{F}(\mathbf{M})\equiv\frac{1}{2}\left(\mathbf{M}-\bm{\tau}_{1}\mathbf{M}^{T}\bm{\tau}_{1}\right),$
(1)

where $\mathbf{M}$ is any $2N\times 2N$ complex matrix. Property (ii) then
says that fermionic BdG matrices are fixed points of this projection, i.e.,
$\mathcal{F}(\mathbf{H}_{F})=\mathbf{H}_{F}$. While this explains the
“fermionic” moniker, the use of the word “projector” follows from the fact
that $\mathcal{F}$, viewed as a linear map on the space of $2N\times 2N$
complex matrices, is idempotent: $\mathcal{F}^{2}=\mathcal{F}$. We may also
define the complementary bosonic projector,

$\mathcal{B}(\mathbf{M})=\frac{1}{2}\left(\mathbf{M}+\bm{\tau}_{1}\mathbf{M}^{T}\bm{\tau}_{1}\right)=\mathbf{M}-\mathcal{F}(\mathbf{M}).$
(2)

It is complementary in the sense that
$\mathbf{M}=\mathcal{F}(\mathbf{M})+\mathcal{B}(\mathbf{M})$. Moreover, these
projectors are orthogonal in the sense that
$\mathcal{B}(\mathcal{F}(\mathbf{M}))=0=\mathcal{F}(\mathcal{B}(\mathbf{M}))$.
Given a matrix $\mathbf{M}$, we call $\mathcal{F}(\mathbf{M})$ and
$\mathcal{B}(\mathbf{M})$ its fermionic and bosonic projections. Fixed points
of these projectors will be called fermionic and bosonic matrices,
respectively. While the “bosonic” moniker will be made clear later, it is
natural to ask what happens if $\mathbf{H}_{F}$ has a non-zero bosonic
projection. If this is the case, then

$\displaystyle
H_{F}=\frac{1}{2}\Psi^{\dagger}\mathbf{H}_{F}\Psi=\frac{1}{2}\Psi^{\dagger}\mathcal{F}(\mathbf{H}_{F})\Phi-\frac{1}{2}\text{tr}[\mathcal{B}(\mathbf{H}_{F})]\,1_{\mathcal{F}}.$

Thus, each Hermitian fermionic matrix $\mathcal{F}(\mathbf{H}_{F})$ defines an
equivalence class of QFHs that differ by a constant shift. We eliminate this
redundancy by requiring $\mathbf{H}_{F}$ be fermionic.

To diagonalize a QFH, one seeks a Bogoliubov transformation mapping the
original degrees of freedom to a set of independent fermionic quasiparticles.
That is, we search for a set of transformed fermionic operators $d_{n}$,
$n=1,\ldots,N$ with

$[H_{F},{\color[rgb]{0,0,0}d_{n}}]=-\epsilon_{n}d_{n},$

with $\epsilon_{n}\geq 0$ the associated quasiparticle energies. If $\ket{E}$
is an eigenstate of $H_{F}$ with energy $E$, then, if $d_{n}\ket{E}$
($d_{n}^{\dagger}\ket{E}$) is non-zero, it is a state with energy
$E-\epsilon_{n}$ ($E+\epsilon_{n}$). That is, $d_{n}$ and $d_{i}^{\dagger}$
annihilate and create a quasiparticle of energy $\epsilon_{n}$. It follows
that the many-body ground state is the state annihilated by each $d_{n}$, and
that eigenstates of $H_{F}$ are built up from it by populating quasiparticles
states.

A simple way to to find these quasiparticles is to introduce the fermionic hat
map. To each numerical vector $\vec{\alpha}\in\mathbb{C}^{2N}$, we associate a
linear form

$\widehat{\vec{\alpha}}\equiv\vec{\alpha}^{\dagger}\Psi=\alpha_{1}^{*}c_{1}+\alpha_{{2}}^{{\color[rgb]{0,0,0}*}}c_{1}^{\dagger}+\cdots+\alpha_{{2N-1}}^{{\color[rgb]{0,0,0}*}}c_{N}+\alpha_{2N}^{*}c_{N}^{\dagger}.$

This map has two notable properties:

$\displaystyle{\widehat{\vec{\alpha}}}^{\dagger}$ $\displaystyle=$
$\displaystyle\widehat{\bm{\tau}_{1}\vec{\alpha}^{*}},$ (3)
$\displaystyle\\{\widehat{\vec{\alpha}},\widehat{\vec{\beta}}\,\mbox{}^{{\dagger}}\\}$
$\displaystyle=$
$\displaystyle\vec{\alpha}^{\dagger}\vec{\beta}\,1_{\mathcal{F}},$ (4)

with the right hand-side of the second equation being the linear form
associated to a numerical vector $\bm{\tau}_{1}\vec{\alpha}^{*}$. With some
algebra, one may further verify that

$[H_{F},\widehat{\vec{\alpha}}]=-\widehat{\mathbf{H}_{F}\vec{\alpha}}.$ (5)

If we then take $\vec{\alpha}=\vec{d}_{n}$, with $\vec{d}_{n}$ an eigenvector
of $\mathbf{H}_{F}$ with eigenvalue $\epsilon_{n}$, it follows that

$[H_{F},\widehat{\vec{d}_{n}}]=-\widehat{\mathbf{H}_{F}\vec{d}_{n}}=-\epsilon_{n}\widehat{\vec{d}_{n}}.$

The eigenvalue equation $\mathbf{H}_{F}\vec{d}_{n}=\epsilon_{n}\vec{d}_{n}$
provides the fermionic BdG equation and resembles the time-independent
Schrödinger equation on a $2N$-dimensional Hilbert space.

Now, since $\mathbf{H}_{F}$ is fermionic, it follows that if $\vec{d}_{n}$ is
an eigenvector with eigenvalue $\epsilon_{n}$, then
${\vec{d}}^{\prime}_{n}\equiv\bm{\tau}_{1}\vec{d}_{n}^{*}$ is an eigenvector
with eigenvalue $-\epsilon_{n}$. We can then construct the quasiparticles as
$d_{n}=\widehat{\vec{d}_{n}}$, with $\vec{d}_{n}$ being the $N$ orthonormal
eigenvectors of $\mathbf{H}_{F}$ with eigenvalues $\epsilon_{n}\geq 0$. From
Eq. (3), it follows that $\widehat{\vec{d}_{n}^{\prime}}=d_{n}^{\dagger}$,
whereby orthonormality and Eq. (4) ensure the quasiparticles satisfy the
canonical anti-commutation relations.

Since the dynamics generated by $H_{F}$ is unitary, determining the evolution
of $\Psi$ suffices to determine the evolution of any observable built up from
products and sums of fermionic creation and annihilation operators. The
Heisenberg equation of motion (in units where $\hbar=1$),

$\displaystyle\frac{d}{dt}\Psi(t)=i[H_{F},\Psi(t)]=-i\mathbf{H}_{F}\Psi(t),$
(6)

can be mapped to a linear time-invariant (LTI) dynamical system on
$\mathds{C}^{2N}$ by using the above hat map. Consider a general linear form
$\widehat{\vec{\alpha}}$ in the Heisenberg picture and the Ansatz [18]

$\widehat{\vec{\alpha}}(t)=\vec{\alpha}^{\dagger}\Psi(t)=\vec{\alpha}^{\dagger}(t)\Psi(0)\equiv\widehat{\vec{\alpha}(t)},$

where the right hand-side is the linear form associated to the now time-
dependent coefficient vector $\vec{\alpha}(t)$. From Eq. (5), it follows that
$\widehat{\vec{\alpha}}(t)$ is a solution to Eq. (6) if and only if
$\vec{\alpha}(t)$ is a solution to the LTI matrix equation

$\dot{\vec{\alpha}}(t)=i\mathbf{H}_{F}\vec{\alpha}(t).$

The quasiparticles are then interpreted dynamically as normal modes of this
LTI system, i.e., $\vec{d}_{n}(t)=e^{i\epsilon_{n}t}\vec{d}_{n}(0)$. Since
$\epsilon_{n}\geq 0$, this corresponds to bounded motion for all times.

Let us now run through the exactly same ideas for bosons. Let $a_{i}$,
$i=1,\dots,N$, denote a set of bosonic annihilation operators satisfying the
canonical commutation relations (CCRs),
$[a_{i},a_{j}^{\dagger}]=\delta_{ij}1_{\mathcal{F}}$ and $[a_{i},a_{j}]=0$,
with $1_{\mathcal{F}}$ now being the bosonic Fock-space identity. A QBH is
given by

$H_{B}=\frac{1}{2}\sum_{i,j=1}^{N}\left(\mathbf{K}_{ij}a_{i}^{\dagger}a_{j}+\bm{\Delta}_{ij}{\color[rgb]{0,0,0}a_{i}^{\dagger}a_{j}^{\dagger}}+\text{H.c.}\right),$
(7)

with $\mathbf{K}=\mathbf{K}^{\dagger}$ the $N\times N$ hopping matrix and
$\bm{\Delta}^{T}=\bm{\Delta}$ the $N\times N$ pairing matrix (note that
bosonic pairing matrices are symmetric, unlike fermionic ones which are
antisymmetric). As with fermions, we define the bosonic Nambu array,

$\Phi\equiv[a_{1},a_{1}^{\dagger},\ldots,a_{N},a_{N}^{\dagger}]^{T}.$ (8)

It then follows that

$H_{B}=\frac{1}{2}\Phi^{\dagger}\mathbf{H}_{B}\Phi,$ (9)

where $\mathbf{H}_{B}$ is the $2N\times 2N$ block matrix with the $(ij)$-th
block given by

$[\mathbf{H}]_{ij}=\begin{bmatrix}\mathbf{K}_{ij}&\bm{\Delta}_{ij}\\\
\bm{\Delta}_{ij}^{*}&\mathbf{K}_{ij}^{*}\end{bmatrix}.$ (10)

The stated conditions on $\mathbf{K}$ and $\bm{\Delta}$ imply that (i) the
matrix $\mathbf{H}_{B}$ is Hermitian; and (ii) bosonic, i.e.,
$\mathbf{H}_{B}=\mathcal{B}(\mathbf{H}_{B})$. Analogously to the fermionic
case, any non-zero fermionic projection of $\mathbf{H}_{B}$ simply shifts the
Hamiltonian according to

$H_{B}=\frac{1}{2}\Phi^{\dagger}\mathbf{H}_{B}\Psi=\frac{1}{2}\Phi^{\dagger}\mathcal{B}(\mathbf{H}_{B})\Phi-\frac{1}{2}\text{tr}[\bm{\tau}_{3}\mathcal{F}(\mathbf{H}_{B})]\,1_{\mathcal{F}}.$
(11)

An equivalence class of QBHs is thus defined uniquely by a Hermitian, bosonic
matrix $\mathbf{H}_{B}$.

Diagonalization proceeds by searching for a Bogoliubov transformation to a set
of bosonic creation and annihilation operators $b_{n}$, $n=1,\ldots,N$
satisfying the CCRs and

$[H_{B},b_{n}]=-\epsilon_{n}b_{n},$ (12)

for some $\epsilon_{n}\in\mathbb{R}$. In sharp contrast to their fermionic
counterparts and Eq.(II.1), this is not always possible, however. To
understand why, we define the bosonic hat map. To each numerical vector
$\vec{\alpha}\in\mathbb{C}^{2N}$, we associate a linear form

$\widehat{\vec{\alpha}}\equiv\vec{\alpha}^{\dagger}\bm{\tau}_{3}\Psi=\alpha_{1}^{*}a_{1}-\alpha_{{2}}^{{\color[rgb]{0,0,0}*}}a_{1}^{\dagger}+\cdots+\alpha_{{2N-1}}^{{\color[rgb]{0,0,0}*}}a_{N}-\alpha_{2N}^{*}a_{N}^{\dagger},$

where $\bm{\tau}_{3}$ is included for later convenience. This map has two
notable properties (cf. Eqs. (3) and (4)):

$\displaystyle\widehat{\vec{\alpha}}\,\mbox{}^{\dagger}$ $\displaystyle=$
$\displaystyle-\widehat{\bm{\tau}_{1}\vec{\alpha}^{*}},$ (13)
$\displaystyle{[}\widehat{\vec{\alpha}},\widehat{\vec{\beta}}\,\mbox{}^{\dagger}{]}$
$\displaystyle=$
$\displaystyle\vec{\alpha}^{\dagger}\bm{\tau}_{3}\vec{\beta}\,1_{\mathcal{F}},$
(14)

With some algebra, one may further verify that

$[H_{B},\widehat{\vec{\alpha}}]=-\widehat{\mathbf{G}\vec{\alpha}},\quad\mathbf{G}\equiv\bm{\tau}_{3}\mathbf{H}_{B},$
(15)

where we have introduced the bosonic BdG Hamiltonian $\mathbf{G}$. The hunt
for quasiparticle excitations then requires solving the bosonic BdG equation,
$\mathbf{G}\vec{b}_{n}=\epsilon_{n}\vec{b}_{n}$. However, unlike the fermionic
BdG Hamiltonian, $\mathbf{G}$ is non-Hermitian whenever bosonic pairing is
present, $\bm{\Delta}\neq 0$. Instead, it is generally pseudo-Hermitian [47],
that is, $\mathbf{G}^{\dagger}=\bm{\tau}_{3}\mathbf{G}\bm{\tau}_{3}$. Non-
Hermiticity eliminates the guarantee that $\mathbf{G}$ has real eigenvalues,
or even that it is diagonalizable at all.

It turns out that the desired Bogoliubov transformation exists if and only if
$\mathbf{G}$ is diagonalizable with an entirely real spectrum [18]. In this
case, there are a set of $N$ eigenvectors $\vec{b}_{n}$ of $\mathbf{G}$
satisfying (i) $\mathbf{G}\vec{b}_{n}=\epsilon_{n}\vec{b}_{n}$; and (ii)
$\vec{b}_{n}^{\dagger}\bm{\tau}_{3}\vec{b}_{m}=\delta_{nm}$. The linear forms
$b_{n}\equiv\widehat{\vec{b}_{n}}$ then satisfy the CCRs and Eq. (12). The
remaining $N$ eigenvectors $\vec{b}_{n}^{\prime}=\bm{\tau}_{1}\vec{b}_{n}^{*}$
correspond to the eigenvalues $-\epsilon_{n}$ and
$\widehat{\vec{b}_{n}^{\prime}}=-b_{n}^{\dagger}$. Constructing the
eigenstates of $H_{B}$ proceeds exactly as in the fermionic case with one
notable exception. If there exists both a positive and a negative
quasiparticle energy, then then $H_{B}$ is thermodynamically unstable
(unbounded in both directions). Thus, there is no well-defined ground state
(nor a well-defined Gibbs state, even if allowing for negative temperatures).
Instead, the state annihilated by each $b_{n}$ is simply a quasiparticle
vacuum.

As a consequence of the fact that $\mathbf{G}$ is, generically, no longer
normal, the dynamical perspective for bosons is much richer than it is for
fermions. While the dynamics generated by $H_{B}$ is still unitary at the
many-body level (allowing us again to build up the dynamics of any observable
from that of $\Phi$), the dynamics of the Nambu array is effectively non-
Hermitian, that is,

$\displaystyle\frac{d}{dt}\Phi(t)=i[{\color[rgb]{0,0,0}H_{B}},\Phi(t)]=-i\mathbf{G}\Phi(t).$
(16)

Adopting the Ansatz
$\widehat{\vec{\alpha}}(t)=\vec{\alpha}^{\dagger}(t)\Phi(0)$, as we did in the
fermionic case, now leads to the LTI dynamical system

$\dot{\vec{\alpha}}(t)=i\mathbf{G}\vec{\alpha}(t).$ (17)

If $\mathbf{G}$ is diagonalizable with entirely real spectrum, then the
eigenvectors associated to the quasiparticles are normal modes of this system
with $\vec{b}_{n}(t)=e^{i\epsilon_{n}t}\vec{b}_{n}(0)$. However, if
$\mathbf{G}$ fails to meet these requirements, there will be (possibly
generalized, in the sense of the Jordan canonical form) eigenvectors with
possibly non-real eigenvalues. In this case, $H_{B}$ is dynamically unstable,
i.e., there exists observable expectation values that exhibit unbounded motion
(see Ref. 18 for a detailed account of dynamical instabilities in QBHs). The
normal mode picture for bosons is thus more general than the quasiparticle
picture, as it persists even when the system is dynamically unstable. For this
reason, it is more appropriate to refer to $\mathbf{G}$ as the dynamical
matrix (as opposed to the BdG Hamiltonian) of the QBH in the general case.

###  II.2 Open Markovian bosonic dynamics

####  II.2.1 Basics of Lindblad formalism

Our focus in this paper will be on Markovian open quantum systems, whose
dynamics are governed by a semigroup Lindblad master equation (LME) of the
form $\dot{\rho}(t)=\mathcal{L}(\rho(t))$, with $\rho(t)$ a density operator
describing the state of the system at time $t>t_{0}\equiv 0$ and the (time-
independent) superoperator $\mathcal{L}$ – the Lindbladian – being the
generator of the dynamics. Physically, a LME provides the most general form of
continuous-time quantum dynamics that obeys complete positivity and trace
preservation. In its canonical (or diagonal) form, the Markovian generator may
be written as [48, 49]

$\displaystyle\mathcal{L}(\rho)$
$\displaystyle=-i[H,\rho]+\sum_{\alpha}\Big{(}L_{\alpha}\rho
L_{\alpha}^{\dagger}-\frac{1}{2}\\{L_{\alpha}^{\dagger}L_{\alpha},\rho\\}\Big{)}$
$\displaystyle\equiv-i[H,\rho]+\sum_{\alpha}\mathcal{D}[L_{\alpha}](\rho),$
(18)

where in second line we have defined the dissipator $\mathcal{D}[A]$. Here,
$H=H^{\dagger}$ represents the Hamiltonian, coherent contribution to the
dynamics, whereas $\\{L_{\alpha}\\}$, $\alpha=1,\ldots,d$, are the Lindblad
(or “jump”) operators, with each $\mathcal{D}[L_{\alpha}]$ phenomenologically
accounting for a distinct, irreversible “noise channel”. It is well known that
the representation ${\cal L}(H,\\{L_{\alpha}\\})$ of the generator in terms of
Hamiltonian and Lindblad operators is not unique, nor is the separation
between a Hamiltonian and dissipative component [49, 50, 51]. In applications,
the choice of a given representation is typically dictated by physical
requirements; for instance, Lindblad operators are often naturally traceless.
In a similar venue, it may be preferable to work in a non-diagonal
representation, whereby both the Hamiltonian and the Lindblad operators are
expressed in terms of a fixed set of operators, which may enjoy special
mathematical (e.g., completeness) or physical (e.g., locality) properties.
Specifically, let $\\{A_{j}\\}$, $j=1,\ldots,d^{\prime}$, such a set of
operators, with $L_{\alpha}\equiv\sum_{j}\ell_{j\alpha}A_{j}$,
$\ell_{j\alpha}\in{\mathbb{C}}$. The Gorini-Kossakowski-Lindblad-Sudarshan
(GKLS) representation for the Markovian generator $\mathcal{L}$ then reads
[52, 49]

$\displaystyle\mathcal{L}(\rho)$
$\displaystyle=-i[H,\rho]+\sum_{jk}\mathbf{M}_{jk}\Big{(}A_{k}\rho
A_{j}^{\dagger}-\frac{1}{2}\\{A_{j}^{\dagger}A_{k},\rho\\}\Big{)}$
$\displaystyle\equiv-i[H,\rho]+\sum_{jk}\mathbf{M}_{jk}\,\mathcal{D}[A_{k},A_{j}^{\dagger}](\rho),$
(19)

where the square, positive-semidefinite (relaxation or GKLS) matrix
$\mathbf{M}_{jk}=\sum_{\alpha=1}^{d}\ell_{j\alpha}^{*}\ell_{k\alpha}$ accounts
for the non-unitary contribution to the dynamics, and we have introduced the
shorthand notation $\mathcal{D}[A,B](\rho)\equiv A\rho B-\\{BA,\rho\\}/2$.
Clearly, $\mathcal{D}[A,A^{\dagger}]=\mathcal{D}[A]$, as previously defined.

In addition to the above Schrödinger dynamics for the (not necessarily pure)
state of the system, one can define an equivalent Heisenberg picture. As
usual, all states are stationary in this picture and it is the operators
associated to physical observables or auxiliary quantities (e.g., the
electromagnetic potential) which evolve in time. An observable $B=B^{\dagger}$
then obeys the equation of motion $\dot{B}(t)=\mathcal{L}^{\star}(B(t))$, with
the dual (Heisenberg) Markovian generator being given by

$\displaystyle\mathcal{L}^{\star}(B)$ $\displaystyle\equiv$ $\displaystyle
i[H,B]+\sum_{\alpha}\Big{(}L_{\alpha}^{\dagger}BL_{\alpha}-\frac{1}{2}\\{L_{\alpha}^{\dagger}L_{\alpha},B\\}\Big{)}$
$\displaystyle=$ $\displaystyle
i[H,B]+\sum_{jk}\mathbf{M}_{jk}\Big{(}A_{k}^{\dagger}BA_{j}-\frac{1}{2}\\{A_{k}^{\dagger}A_{j},B\\}\Big{)},$

in its canonical and GKLS form, respectively. As for unitary dynamics, the
mathematical relationship between $\mathcal{L}$ and $\mathcal{L}^{\star}$
follows from demanding that expectation values agree in the two pictures at
all times, that is, $\text{tr}(B(t)\rho)=\text{tr}(B\rho(t))$, for all
observables and density operators. Unlike the unitary case, however, Markovian
evolution is not multiplicative in general: that is, generically,
$(B_{1}B_{2})(t)\neq B_{1}(t)B_{2}(t)$.

If the underlying Hilbert space $\mathcal{H}$ is finite-dimensional, say, of
dimension $D$, the LME describes an LTI dynamical system whose generator
$\mathcal{L}$ may be thought of as a (generally) non-Hermitian linear operator
acting on a $D^{2}$-dimensional space. A number of general conclusions can
then be made about the spectral properties of $\mathcal{L}$ and, in turn, the
ensuing transient and asymptotic relaxation dynamics into the steady state
manifold [53, 54]. In particular, the SS manifold always comprises at least
one SS, say, $\rho_{\text{ss}}$, satisfying $\mathcal{L}(\rho_{\text{ss}})=0$.
Furthermore, the spectrum of the Lindbladian, $\sigma(\mathcal{L})$, is bound
to the closed left-half complex-plane, with any non-real eigenvalues existing
in complex conjugate pairs. In the important situation where the only element
in $\sigma(\mathcal{L})$ with vanishing real part is zero, the SS is unique,
and arbitrary initial states $\rho_{0}$ asymptotically undergo exponential
relaxation to $\rho_{\text{ss}}$, that is, $\rho_{\text{ss}}$ is globally
asymptotically stable [50]. Let the spectral gap (or dissipative gap) of
$\mathcal{L}$ be defined by

$\Delta_{\mathcal{L}}\equiv|\sup\text{Re}[\sigma(\mathcal{L})\setminus\\{0\\}]|\geq
0,$ (20)

that is, by the closest distance between the imaginary axis and the set of
non-zero eigenvalues of $\mathcal{L}$. Then, $\Delta_{\mathcal{L}}$ sets the
asymptotic decay rate through the inequality

$d_{\text{max}}(t)\equiv\sup_{\rho_{0}}\lVert
e^{t\mathcal{L}}(\rho_{0})-\rho_{\text{ss}}\rVert_{\text{tr}}\leq
Ke^{-\Delta_{\mathcal{L}}t},\quad K>0,$

where the trace distance $\lVert
A\rVert_{\text{tr}}\equiv\text{tr}[\sqrt{A^{\dagger}A}]$. The minimum time it
takes for the above $d_{\text{max}}(t)$ to fall below a predetermined accuracy
$\epsilon>0$ defines the mixing time of the Markovian semigroup, in direct
analogy to classical Markov processes. Since the prefactor $K$ may be a priori
very large, the mixing time plays a central role in characterizing the
convergence behavior, by also accounting for the non-exponential pre-
relaxation regime. Notably, dissipative “quasi-free” fermionic dynamics,
described by QFLs, are known to exhibit rapid convergence to stationarity
[55], including via the occurrence of cutoff phenomena [56].

If $\mathcal{H}$ is infinite-dimensional, as it is necessarily the case for
bosonic systems, a rigorous mathematical description becomes substantially
more challenging and major departures from the above picture may arise. On the
one hand, proving the existence of a generator of the form in Eq. (18)
requires appropriate boundedness assumptions to hold; on the other hand, even
when a Markovian generator can be shown or assumed to exist, the resulting
Lindbladian may lack a SS, and relating spectral properties to dynamical
evolution is far more involved in general [57]. Leaving mathematical rigor
aside, LMEs have been remarkably successful in describing a wide class of open
quantum systems, notably, across quantum optics and photonics [58, 10, 59].
Infinite dimensionality opens up the possibility for the Hamiltonian and other
observables to be unbounded. As a consequence, ${\cal L^{\star}}$ may now
obtain eigenvalues with strictly positive real part, e.g., in the case of
optical pumping with $H=\omega a^{\dagger}a$ and $L=\sqrt{\kappa}a^{\dagger}$.
With these caveats in mind, we now narrow our focus to the class of quadratic
dissipative bosonic systems we will be specifically interested in throughout
this work.

####  II.2.2 Quadratic bosonic Lindbladians

As before, let $a_{j}^{\dagger}$ and $a_{j}$ denote canonical creation and
annihilation operators for bosons, arranged in a bosonic Nambu array
$\Phi^{\dagger}$ as in Eq. (8). The bosons can be regarded as quasi-free if
the Lindblad generator is of the quadratic form

$\mathcal{L}(\rho)=-i[H,\rho]+\sum_{i,j=1}^{2N}\mathbf{M}_{ji}\mathcal{D}(\Phi_{i},\Phi_{j}^{\dagger})(\rho),$
(21)

with $H=\frac{1}{2}\Phi^{\dagger}\mathbf{H}\Phi$ being a QBH. We will call
$\mathcal{L}$ (and its adjoint $\mathcal{L}^{\star}$) quadratic bosonic
Lindbladians (QBLs). We note that a QBL so defined is not the most general
instance of quasi-free bosonic dynamics, as it is also possible for $H$ to
include a term linear in the creation and annihilation operators [28]. In the
context of continuous-variable quantum information, such quasi-free Lindblad
generators are also often referred to as Gaussian [60]. Notably, Gaussian
Markovian generators may be equivalently defined in terms of their
preservation of the set of Gaussian states [29]. Owing to the non-uniqueness
in representing the generator, linear Hamiltonian contributions can often be
removed by way of an appropriate constant shift of the Lindblad operators
[50], hence the creation and annihilation operators. In any case, our primary
focus here will be purely quadratic Gaussian generators, of the form given in
Eq. (21).

Gaussian LMEs are exactly solvable thanks to the fact that a quadratic
Lindblad generator maps any polynomial in the creation and annihilation
operators to another such polynomial of the same or lesser degree. As a rule,
the dynamics of a form of odd (even) order depends on all the forms of odd
(even) order of lesser degree. Hence, the linear forms satisfy a closed
equation of motion and so do quadratic forms, up to a contribution
proportional to the identity operator (the form of degree zero). To obtain a
closed system of equations for cubic forms, one needs to consider cubic and
and linear forms together, and so on. In short, one can map Gaussian LME to a
hierarchy of linear ordinary differential equations. In this sense, a Gaussian
LME is exactly solvable. However, fully taking advantage of these facts
requires some machinery.

Our discussion in Sec. II.1 showed that fermionic projections $\mathcal{F}$
play no role in the dynamics of closed bosonic systems. Remarkably,
dissipation brings them right back into the picture to play a very distinct
role. Let $\mathcal{D}^{\star}$ denote the dissipator of a general QBL. Then,
by writing $\mathbf{M}=\mathcal{B}(\mathbf{M})+\mathcal{F}(\mathbf{M})$ in Eq.
(21), and leveraging the CCRs, one can show that

$\displaystyle\mathcal{D}^{\star}(A)=\frac{1}{2}\sum_{i,j}\mathcal{B}(\mathbf{M})_{ij}\\!\left(\left[[\Phi_{i}^{\dagger},A],\Phi_{j}\right]+\left[\Phi_{i}^{\dagger},[A,\Phi_{j}]\right]\right)$
$\displaystyle+\,\mathcal{F}(\mathbf{M})_{ij}\\!\left(\left\\{[\Phi_{i}^{\dagger},A],\Phi_{j}\right\\}+\left\\{\Phi_{i}^{\dagger},[A,\Phi_{j}]\right\\}\right)$
,

with complete generality. Thinking of $A$ as a form of degree $d$, we see that
the first term reduces the degree of $A$ by 2, while the second preserves the
degree. Thus, the bosonic projection of $\mathbf{M}$ is responsible for
connecting the degree-$d$ and degree-$(d-2)$ sectors of the system’s operator
algebra, while the fermionic projections leaves the degree-$d$ sector
invariant.

Finally, we are ready to present the equations of motion for linear and
quadratic forms. Let $Q\equiv\Phi\Phi^{\dagger}$ denote the square array of
elementary quadratic forms $a_{i}a_{j}$ and so on. Then, the Heisenberg
equations of motion for arbitrary linear and quadratic form follow from the
compact formulas

$\displaystyle\left\\{\begin{array}[]{l}\\!\dot{\Phi}(t)=\mathbf{G}\Phi(t),\\\
\\!\dot{Q}(t)=-i(\mathbf{G}Q(t)-Q(t)\mathbf{G}^{\dagger})+\bm{\tau}_{3}\mathbf{M}\bm{\tau}_{3}\,1_{F},\end{array}\right.$
(24)

where

$\mathbf{G}\equiv\bm{\tau}_{3}\mathbf{H}-i\bm{\tau}_{3}\mathcal{F}(\mathbf{M})$
(25)

is the dynamical matrix associated to the QBL. Remarkably, the dynamics of a
linear bosonic form is entirely controlled by a fermionic matrix in the case
of pure dissipation ($\mathbf{H}=0)$. For future reference, we note the
relationships

$\mathcal{L}^{\star}(\widehat{\vec{\alpha}})=\widehat{i\widetilde{\mathbf{G}}\vec{\alpha}},\quad\widetilde{\mathbf{G}}\equiv\bm{\tau}_{3}\mathbf{G}^{\dagger}\bm{\tau}_{3}.$
(26)

It follows that the dynamics generated by $\mathcal{L}^{\star}$, restricted to
the subalgebra of linear forms, are encoded by an associated LTI system within
the BdG space $\mathbb{C}^{2N}$ generated by $i\widetilde{\mathbf{G}}$. That
is,

$\widehat{\vec{\alpha}}(t)=\widehat{\vec{\alpha}(t)}\Longleftrightarrow\dot{\vec{\alpha}}(t)=i\widetilde{\mathbf{G}}\vec{\alpha}(t).$
(27)

Notice that $\widetilde{\mathbf{G}}=\mathbf{G}$ if and only if $\mathcal{D}=0$
(cf. Eq. (17)).

Since the bosonic Fock space is infinite-dimensional, a QBL can become
dynamically unstable. Then, by definition, there exists an observable $A$ and
a state $\rho$ such that $\braket{A\,}\\!(t)$ is unbounded in time. As a rule,
dynamical instability can be diagnosed through the dynamical matrix by way of
the rapidity spectrum, that is, $\sigma(-i\mathbf{G})$. In addition to the
spectral gap $\Delta_{\cal L}$ of Eq. (20), the stability gap is defined by

$\Delta_{S}\equiv\max\text{Re}\,\sigma(-i\mathbf{G}).$ (28)

If $\Delta_{S}<0$ ($\Delta_{S}>0$, then the QBL is dynamically stable
(unstable). The crossover regime $\Delta_{S}=0$ can be either stable or
unstable contingent on the existence of a steady state [28]. If the QBL is
dynamically stable, then the stability gap determines the spectral gap of
$\mathcal{L}$ by way of the simple formula $\Delta_{\mathcal{L}}=-\Delta_{S}$
[61]. Moreover, there exists a unique, globally attractive, and Gaussian
$\rho_{ss}$. Its first and second moments are stationary solutions to Eqs.
(24), that is,

$\displaystyle\left\\{\begin{array}[]{l}\mathbf{G}\braket{\Phi}_{ss}=0,\\\
-i\mathbf{G}\\!\braket{Q}_{ss}+i\braket{Q}_{ss}\\!\mathbf{G}^{\dagger}+\bm{\tau}_{3}\bf{M}\bm{\tau}_{3}=0.\end{array}\right.$
(31)

The condition of a strictly negative stability gap identifies $-i\mathbf{G}$
as a Hurwitz matrix in the language of matrix analysis. Then, uniqueness of
the SS follows from the uniqueness of the solutions of the above two equations
when $-i\mathbf{G}$ is Hurwitz [61]. In particular, $\braket{\Phi}=0_{ss}$
necessarily. If $\Delta_{S}\geq 0$, there can be either zero or infinitely
many SSs [28]. Since $\braket{\mathcal{F}(Q)}=\bm{\tau}_{3}/2$ automatically
in any state, one can focus on $\braket{B(Q)}_{ss}$ which, after an
appropriate basis transformation, is unitarily equivalent to the SS covariance
matrix [62].

####  II.2.3 Symmetries, conserved quantities, and
bulk-translationally invariant QBLs

Just as in the closed-system case, symmetries and conserved quantities provide
crucial information about the system. While the notion of a conserved
quantity, i.e., an observable $Q$ such that $\dot{Q}(t)=0$, has a
straightforward characterization ($\mathcal{L}^{\star}(Q)=0$), two distinct
notions of symmetry emerge in the Markovian setting [63, 64, 54]. A unitary or
anti-unitary operator $S$ is a weak symmetry of the dynamics if, for all
$\rho$, $(S\rho S^{-1})(t)=S\rho(t)S^{-1}$ or, equivalently,
$[\mathcal{L},\mathcal{S}]=0$, with $\mathcal{S}(\rho)\equiv S\rho S^{-1}$.
Further, $S$ is a strong symmetry if $[H,S]=0$ and
$[L_{\color[rgb]{0,0,0}\alpha},S]=0$, for all ${\color[rgb]{0,0,0}\alpha}$.
For simplicity, unless otherwise noted, in this work we will generically refer
to a weak symmetry as just a symmetry.

If the dynamics possesses a continuous unitary family of strong symmetries
$S(\theta)$, the existence of a corresponding conserved quantity is
guaranteed. In fact, the conserved quantity $Q$ coincides with the generator
of the one-parameter group $S(\theta)=e^{i\theta Q}$. This fact generalizes
naturally to higher-dimensional Lie groups and is reminiscent of Noether’s
theorem, which establishes the existence of conserved quantities given a
continuous symmetry. However, in sharp contrast to the Hamiltonian case, (i)
the existence of a continuous weak symmetry does not imply the existence of a
conserved quantity; and (ii) the existence of a conserved quantity does not
imply the existence of either a weak or strong continuous symmetry. In this
sense, we say the quantum Hamiltonian form of Noether’s theorem breaks down in
the presence of Markovian dissipation. We further note that
$S(\theta)=e^{i\theta Q}$ is a weak symmetry of the dynamics if and only if

$\mathcal{L}^{\star}([Q,A])-[Q,\mathcal{L}^{\star}(A)]=0,\qquad\forall A,$
(32)

that is, $\mathcal{L}^{\star}$ commutes with the adjoint action of $Q$.

The class of QBLs we consider are further characterized by (weak) discrete
bulk-translation symmetry in 1D. Here, bulk-translation invariance refers to
the presence of translation invariance when suitably far from a well-defined
boundary. Specifically, bulk-translationally invariant QBLs may correspond to
one of four main configurations:

  * •

A finite number of lattice sites on a chain (open boundary conditions, OBCs);

  * •

A finite number of lattice sites on a ring (periodic boundary conditions,
PBCs);

  * •

An infinite number of lattice sites extending only in one direction (semi-
infinite boundary conditions, SIBCs);

  * •

An infinite number of lattice sites extending in two directions (bi-infinite
boundary conditions, BIBCs).

All the Hamiltonian and dissipative couplings are assumed to be independent of
lattice site, modulo BCs, and of a finite range $R$, with $R\ll N$ [65]. In
this way, the equation of motion for operators at site $j$ only involve
operators at most a distance $R$ away, e.g., $R=1$ for nearest-neighbor (NN)
hopping. To each lattice site $j$, let us associate $d_{\text{int}}$ bosonic
degrees of freedom arranged in a site-local Nambu array,
$\phi_{j}\equiv[a_{j,1},a_{j,1}^{\dagger},\ldots,a_{j,d_{\text{int}}},a_{j,d_{\text{int}}}^{\dagger}]^{T}$.
The Hamiltonian contribution is then constrained to a particular simple form,

$H=\frac{1}{2}\sum_{j}\sum_{r=-R}^{R}\phi_{j}^{\dagger}\mathbf{h}_{r}\phi_{j+r},$

with $\mathbf{h}_{r}$ a $2d_{\text{int}}\times 2d_{\text{int}}$ matrix that
encodes the coherent hopping and pairing mechanisms between sites $j$ and
$j+r$. Necessarily, $\mathbf{h}_{r}^{\dagger}=\mathbf{h}_{-r}$ and
$\mathbf{h}_{r}^{*}=\bm{\tau}_{1}\mathbf{h}_{r}\bm{\tau}_{1}$. In the above
expression, the sum over $j$ encodes the BCs. Likewise, the Heisenberg-picture
dissipator must take the form

$\mathcal{D}^{\star}(A)=\sum_{j}\sum_{r=-R}^{R}\Big{(}\phi_{j}^{\dagger}A\mathbf{m}_{r}\phi_{j+r}-\frac{1}{2}\Big{\\{}\phi_{j}^{\dagger}\mathbf{m}_{r}\phi_{j+r},A\Big{\\}}\Big{)},$

with $\mathbf{m}_{r}$ a $2d_{\text{int}}\times 2d_{\text{int}}$ matrix that
encodes the incoherent/dissipative damping, pumping, and pairing mechanisms
between sites $j$ and $j+r$. Necessarily,
$\mathbf{m}_{r}^{\dagger}=\mathbf{m}_{-r}$, and again, the $j$-sum encodes the
BCs.

Working with a fixed number of lattice sites $j=1,\ldots,N$, the dynamical
matrix and the GKLS matrix are then given by

$\displaystyle\mathbf{G}_{N}$ $\displaystyle=$
$\displaystyle\mathds{1}_{N}\otimes\mathbf{g}_{0}+\sum_{r=1}^{R}\mathbf{S}_{N}^{r}\otimes\mathbf{g}_{r}+\mathbf{S}_{N}^{\dagger}{}^{r}\otimes\mathbf{g}_{-r},$
$\displaystyle\mathbf{M}_{N}$ $\displaystyle=$
$\displaystyle\mathds{1}_{N}\otimes\mathbf{m}_{0}+\sum_{r=1}^{R}\mathbf{S}_{N}^{r}\otimes\mathbf{m}_{r}+\mathbf{S}_{N}^{\dagger}{}^{r}\otimes\mathbf{m}_{-r},$

in terms of the BC-dependent left-shift operator

$\mathbf{S}_{N}\equiv\sum_{j=1}^{N-1}\vec{e}_{j}\vec{e}_{j+1}^{\,{\dagger}}+\begin{cases}0,&\text{OBCs},\\\
\vec{e}_{N}\vec{e}_{1}^{\,{\dagger}},&\text{PBCs},\end{cases}$

with $\vec{e}_{j}$ the $j$’th canonical basis vector of $\mathbb{C}^{N}$. When
focusing exclusively on OBCs (PBCs), we will usually denote
$\mathbf{S}_{N}=\mathbf{T}_{N}$ ($\mathbf{V}_{N}$). By using Eq. (25), the
internal coupling matrices of the dynamical matrix are identified by letting

$\mathbf{g}_{r}\equiv\bm{\tau}_{3}\mathbf{h}_{r}-\frac{i}{2}\bm{\tau}_{3}(\mathbf{m}_{r}-\bm{\tau}_{1}\mathbf{m}^{*}_{r}\bm{\tau}_{1}),$

with $\bm{\tau}_{j}=\mathds{1}_{d_{\text{int}}}\otimes\bm{\sigma}_{j}$ in this
context. Notably, we have
$\mathbf{g}_{r}^{*}=-\bm{\tau}_{1}\mathbf{g}_{r}\bm{\tau}_{1}$. Under OBCs
(PBCs), the matrices $\mathbf{h}_{r},\mathbf{m}_{r}$ belong to the class of
banded block-Toeplitz (banded block-circulant) matrices [65]. We label semi-
infinite and bi-infinite lattices according to $j=1,2,\ldots,\infty$ and
$j=0,\pm 1,\pm 2,\ldots,\pm\infty$, respectively. In the case of SIBCs
(BIBCs), the dynamical matrices belong to the classes of banded block-Toeplitz
(banded block-Laurent) operators. For OBCs, PBCs, SIBCs, and BIBCs, we denote
the matrices $\mathbf{X}=\mathbf{G},\mathbf{M}$ by
$\mathbf{X}^{\text{OBC}}_{N}$, $\mathbf{X}^{\text{PBC}}_{N}$,
$\mathbf{X}^{\text{SIBC}}$, and $\mathbf{X}^{\text{BIBC}}$, respectively (see
Table 2). When the particular BC is unimportant, or when we wish to refer to
the full family of configurations, we simply write $\mathbf{X}$.

Boundary conditions | Dynamical/GKLS matrix type
---|---
Open | Banded block-Toeplitz matrix
Periodic | Banded block-circulant matrix
Semi-infinite | Banded block-Toeplitz operator
Bi-infinite | Banded block-Laurent operator
Table 2: The correspondence between BCs and the structure of the dynamical and
GKLS matrices $\mathbf{G}$ and $\mathbf{M}$ in bulk-translationally invariant
QBLs.

In the translationally invariant cases (PBCs and BIBCs), we may introduce the
conserved crystal momenta $k$ through Fourier modes,
$b_{k,m}\equiv\sum_{j}e^{-ijk}a_{j,m}$. Here, $k$ takes on either discrete
(PBCs) or continuous (BIBCs) values in the Brillouin zone $[-\pi,\pi]$. The
dynamics of the Fourier modes, organized again in the form of a Nambu array,
$\widetilde{\phi}_{k}\equiv[b_{k,1},b_{-k,1}^{\dagger},\ldots,b_{k,d_{\text{int}}},b_{-k,d_{\text{int}}}^{\dagger}]^{T}$,
are then governed by a $2d_{\text{int}}\times 2d_{\text{int}}$ $k$-dependent
dynamical matrix $\mathbf{g}(k)$, that is,

$\dot{{\widetilde{\phi}}}_{k}(t)=-i\mathbf{g}(k)\widetilde{\phi}_{k}(t),\quad\mathbf{g}(k)\equiv\sum_{r=-R}^{R}\mathbf{g}_{r}e^{ikr}.$
(33)

We call $\mathbf{g}(k)$ the Bloch dynamical matrix to draw analogy with the
Bloch Hamiltonian of condensed-matter systems. Similarly, we define the
rapidity bands, which are the dissipative generalization of energy bands, to
be the eigenvalues of $-i\mathbf{g}(k)$. These provide the rapidities of the
translationally invariant configurations PBCs and BIBCs. Consistent with the
known symmetries of Lindbladian eigenvalues, the property
$\mathbf{G}^{*}=-\bm{\tau}_{1}\mathbf{G}\bm{\tau}_{1}$ (or, equivalently,
$\mathbf{g}(k)^{*}=-\bm{\tau}_{1}\mathbf{g}(-k)\bm{\tau}_{1}$ with
$\bm{\tau}_{1}$ understood here to be $2d\times 2d$) ensures that if
$\lambda(k)$ is a rapidity band, then so is $\lambda(k)^{*}$.

Thus, by exploiting translational invariance, a simple description of the
rapidities for PBCs and BIBCs is possible, and solving for the dynamics of a
large class of operators of interest amounts to solving the finite-dimensional
LTI system in Eq. (33). The dynamics of QBLs for which translational symmetry
is broken by OBCs or SIBCs are considerably more difficult to describe. While
the spectral theory of block-Toeplitz matrices and operators is well-
established in the mathematical-physics literature [65, 30, 31], it turns out
that spectral properties are not, in general, the appropriate tool to
consider, due to the fact that the relevant dynamical matrices need not be
Hermitian or even normal.

###  II.3 The pseudospectrum

Pseudospectral theory is an extension of spectral theory adapted to non-normal
matrices and operators. As our work will illustrate, it is a powerful
alternative to the traditional spectral analysis in terms of eigenvalues and
invariant subspaces, because it can predict approximate dynamical transient
behavior that is hard to identify directly from the exact normal modes of the
dynamical system. In addition, the notion of the pseudospectrum can
accommodate the idea of an “approximate mode” in a mathematically sharp
framework.

Let $\mathbf{X}$ be an $n\times n$ complex matrix. A complex number $\lambda$
is in the spectrum of $\mathbf{X}$, $\lambda\in\sigma(\mathbf{X})$, if
$\mathbf{X}-\lambda\mathds{1}_{n}$ is not invertible. To motivate the notion
of the pseudospectrum, let us recast this definition in an equivalent way: if
$\\{\lambda_{i}\\}$ is a a sequence converging to $\lambda$, and
$\mathbf{X}-\lambda_{i}$ is invertible for all $i$, then
$\lambda\in\sigma(\mathbf{X})$ if
$\lVert(\mathbf{X}-\lambda_{i})^{-1}\rVert\to\infty$ for some matrix norm. The
choice of norm is of no consequence if ${\cal H}$ is finite-dimensional,
because all norms are equivalent (induce the same topology). This definition
of the spectrum suggests a natural generalization. Given some fixed matrix
norm $\lVert\cdot\rVert$ and $\epsilon>0$, the
$(\epsilon,\lVert\cdot\rVert)$-pseudospectrum of $\mathbf{X}$ is defined as

$\sigma_{\epsilon,\lVert\cdot\rVert}(\mathbf{X})\equiv\set{\lambda\in\mathbb{C}:\lVert(\mathbf{X}-\lambda\mathds{1}_{n})^{-1}\rVert>1/\epsilon},$
(34)

with the understanding that $\lVert(\mathbf{X}-\lambda)^{-1}\rVert=\infty$ for
$\lambda\in\sigma(\mathbf{X})$; that is, the spectrum is always a subset of
the pseudospectrum. The elements of the pseudospectrum are the
$(\epsilon,\lVert\cdot\rVert)$-eigenvalues, or pseudoeigenvalues of
$\mathbf{X}$. If $\lVert\cdot\rVert$ is induced by a vector norm, that is,
$\lVert\mathbf{X}\rVert\equiv\text{sup}_{\lVert\vec{v}\rVert=1}\lVert\mathbf{X}\vec{v}\rVert$,
we have the equivalent (and more useful for our applications) definition:

$\sigma_{\epsilon,\lVert\cdot\rVert}(\mathbf{X})=\set{\\!\lambda\in\mathbb{C}:\exists\vec{v}\in\mathbb{C}^{n},\,\lVert\vec{v}\rVert=1,\,\lVert(\mathbf{X}-\lambda)\vec{v}\rVert<\epsilon}.$
(35)

The normalized vectors $\vec{v}$, with
$\lVert(\mathbf{X}-\lambda)\vec{v}\rVert<\epsilon$, are the
$(\epsilon,\lVert\cdot\rVert)$-pseudoeigenvectors associated to the
$(\epsilon,\lVert\cdot\rVert)$-eigenvalue $\lambda$. When $\epsilon$,
$\lVert\cdot\rVert$ are either understood or inconsequential, we will drop one
or both of them in our notations.

The pseudospectrum becomes increasingly relevant as the matrix $\mathbf{X}$
becomes highly nonnormal. One reason is the way it relates to perturbations.
It follows from the definition that if
$\lambda\in\sigma_{\epsilon,\lVert\cdot\rVert}(\mathbf{X})$, then there exists
a perturbation $\mathbf{E}$ of size $\lVert\mathbf{E}\rVert<\epsilon$, such
that $\lambda\in\sigma(\mathbf{X}+\mathbf{E})$ [30]. Hence,

$\sigma_{\epsilon,\lVert\cdot\rVert}(\mathbf{X})=\bigcup_{\mathbf{E}:\,\lVert\mathbf{E}\rVert<\epsilon}\sigma(\mathbf{X}+\mathbf{E}).$
(36)

Another reason can be seen simply by considering the case where
$\lVert\cdot\rVert=\lVert\cdot\rVert_{2}$ is the matrix norm induced by the
standard (Euclidean) vector $2$-norm. From Eq. (36), a perturbation of size
$\epsilon$ can only shift the spectrum of a normal matrix by at most
$\epsilon$. However, if $\mathbf{X}$ is nonnormal, small perturbations can
drastically modify the spectrum – one manifestation of this being the non-
Hermitian skin-effect (NHSE) [66, 67, 68]. In this case, the pseudospectra can
dramatically influence the transient dynamics of an LTI system generated by
$\mathbf{X}$, e.g., $\dot{\vec{v}}=\mathbf{X}\vec{v}$. In particular, the
pseudospectra bounds the maximal dilation of norm in the sense of

$\displaystyle\sup_{t\geq 0}\lVert
e^{t\mathbf{X}}\rVert\geq\frac{\alpha_{\epsilon}(\mathbf{X})}{\epsilon},\quad\alpha_{\epsilon}(\mathbf{X})\equiv\sup\text{Re}\,\sigma_{\epsilon}(\mathbf{X}).$
(37)

The quantity $\alpha_{\epsilon}$ is called the pseudospectral abscissa and
measures the extent to which the $\epsilon$-pseudospectrum of $\mathbf{X}$
extends towards, or into, the right-half complex plane [30].

The bound of Eq. (37) is particularly relevant for Hurwitz matrices, that is,
matrices with spectrum bounded in the left-half complex plane. If
$\alpha_{\epsilon}(\mathbf{X})$ is positive and large compared to $\epsilon$,
$\lVert e^{t\mathbf{X}}\rVert$ will experience transient growth before
asymptotically decaying to zero. That is, highly non-normal, but
asymptotically stable, dynamical systems can appear unstable during a
transient period. In fact, given a particular $\epsilon$-pseudoeigenvector
$\vec{v}$ with $\epsilon$-pseudoeigenvalue $\lambda$, we have

$\\!\\!\lVert e^{t\mathbf{X}}\vec{v}-e^{\lambda
t}\vec{v}\rVert\leq\lVert(\mathbf{X}-\lambda)\vec{v}\rVert
t+\mathcal{O}(t^{2})<\epsilon t+\mathcal{O}(t^{2}),$ (38)

so that $\vec{v}$ evolves like a normal mode with eigenfrequency $\lambda$ for
sufficiently small (set by $\epsilon$) timescales. We call such modes
pseudonormal modes. We note that $2$-norm pseudospectra can be directly
related to singular values, as we briefly discuss in App. B.1.

The theory of pseudospectra is especially well-developed and useful for the
Toeplitz matrices and operators which, as we have seen, describe bulk-
translationally invariant QBLs. While this theory applies generally to block-
Toeplitz matrices [30, 31], we recount here only the non-block case for
simplicity. Let $\mathbf{X}_{N}^{\text{OBC}}$, denote a banded, $N\times N$
Toeplitz matrix, $\mathbf{X}^{\text{SIBC}}$ its associated Toeplitz operator,
and $\mathbf{X}^{\text{BIBC}}$ its associated Laurent operator, respectively.
Explicitly,

$\mathbf{X}_{N}^{\text{OBC}}=x_{0}\mathds{1}_{N}+\sum_{r=1}^{R}(x_{r}\mathbf{T}_{N}^{r}+x_{-r}\mathbf{T}_{N}^{\dagger}{}^{r}),\quad
x_{0},x_{\pm r}\in{\mathbb{C}},$

with $R$ the range. With this, we may define the complex-valued symbol
$x(z)\equiv\sum_{r=-R}^{R}x_{r}z^{r}$, with $z\in\mathbb{C}$. Upon comparison,
it follows that $x(e^{ik})$ is precisely the non-block analogue of the Bloch
dynamical matrix we introduced for translationally invariant QBLs. As in that
case 333To 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_{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}$., one
finds that $x(e^{ik})$ determines the spectrum of the Laurent operator
$\mathbf{X}^{\text{BIBC}}$. The spectrum of $\mathbf{X}^{\text{SIBC}}$ may
also be characterized as the spectrum of $\mathbf{X}^{\text{BIBC}}$, together
with all complex numbers about which the symbol winds. Explicitly, let the
integer-valued winding number of $x(e^{ik})$ about a complex number $\lambda$
not in the image of the unit circle under $x$, i.e.,
$\lambda\not\in\\{x(e^{ik}):k\in[-\pi,\pi]\\}$, be given by

$\displaystyle\nu(\lambda)$ $\displaystyle\equiv$ $\displaystyle\frac{1}{2\pi
i}\oint_{|z|=1}\\!\frac{x^{\prime}(z)}{x(z)-\lambda}dz$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi
i}\int_{-\pi}^{\pi}\\!\frac{d}{dk}\ln[x(e^{ik})-\lambda]\,dk.$

Then, if $\nu(\lambda)\neq 0$, $\lambda$ is in the spectrum of
$\mathbf{X}^{\text{SIBC}}$.

Figure 1: Representative behavior of rapidity spectra under different BCs, for
a QBL describing the dissipative bosonic Kitaev chain model of Sec. V.3. The
solid curves give the bulk BIBC spectra, whereas the points along the curves
are the rapidities for PBCs. The points on the vertical lines are the
rapidities for OBCs. The shaded regions denote the SIBC spectra. (a) The
doubly-degenerate rapidity spectrum when $\Gamma=0$. The filled (open) markers
represent the topologically metastable (anomalously relaxing) regime, with
$\kappa/\Delta=0.6$ $(1.4)$. (b) The doubly-degenerate rapidity spectrum when
$\Gamma=0.12$, so that the winding around $\lambda=0$ is zero. This represents
non-topological dynamical metastability. In all cases, $J=2$, $\Delta=0.5$,
$\mu=0$, and $N=25$. Adapted from Ref. [24].

Finding a general and sharp characterization of the spectrum of
$\mathbf{X}_{N}^{\text{OBC}}$ remains an open problem. At best, it is known
that the eigenvalues fall on curves contained in the SIBC spectrum as
$N\to\infty$. Shockingly, however, the finite-$N$ spectra need not converge to
the semi-infinite spectra, as illustrated in the Fig. 1. The spectra shown are
those of Toeplitz matrices associated to the dynamical matrix of a particular
QBL model to be considered in Sec. V.3. In this example, the Toeplitz (OBC)
spectra of the system with $N$ sites fall on a vertical line contained with in
the 2D shaded area corresponding to the SIBC spectrum. This suggests that the
spectral properties of the infinite-size limit are generically distinct from
those of the corresponding finite-size systems. In fact, a precise
mathematical relationship does exist when one zooms out to the pseudospectra.
Namely, we have [30]

$\lim_{\epsilon\to
0}\lim_{N\to\infty}\sigma_{\epsilon}(\mathbf{X}_{N}^{\text{OBC}})=\sigma(\mathbf{X}^{\text{SIBC}})\neq\lim_{N\to\infty}\lim_{\epsilon\to
0}\sigma_{\epsilon}(\mathbf{X}_{N}^{\text{OBC}}).$ (39)

That is, the two limits do _not_ commute in general: While the right hand-side
reflects the generic discontinuity in finite-size spectra as $N\to\infty$, the
left hand-side of the above equation can be understood as follows. Given
$\lambda$ in the SIBC spectrum and an arbitrary $\epsilon>0$, then $\lambda$
is in the $\epsilon$-pseudospectra of the finite-size OBC system for all
suitably large $N$. This can be understood intuitively by considering an edge-
localized (see Sec. IV for a brief discussion on localization) eigenvector of
the semi-infinite system corresponding to an eigenvalue $\lambda$ not in the
bulk bands, and then projecting this eigenvector onto a finite chain of length
$N$. Generically, this will not provide an exact eigenvector of the finite
chain but, rather, an approximate one with pseudoeigenvalue $\lambda$ and
accuracy $\epsilon$ set by the localization length of the original
eigenvector. In this sense, we can say that the SIBC spectrum “imprints
itself” into the $\epsilon$-pseudospectra of the finite-size OBC system, for
arbitrary $\epsilon$.

##  III Some fundamental results about QBLs

###  III.1 A correspondence between conserved quantities and symmetry
generators

The edge-localized Majorana ZMs of a free-fermion Hamiltonian are an
indication that the system is in an SPT phase. A mode is, by definition, a
linear combination of creation and annihilation operators and a ZM, in
particular, commutes with the Hamiltonian. There are two other important
classes of operators that commute with the Hamiltonian. They are the
generators of continuous symmetry groups and the observables associated to
conserved quantities. Hence, for Hamiltonian systems, three conceptually quite
different objects arise as the solutions of one and the same equation:
$[H,A]=0$.

The situation is different for Markovian systems. The symmetry generators and
conserved quantities of a LME are characterized by two very different
equations which need not share any solutions [64]. On the one hand, the
adjoint action of a symmetry generator commutes with $\mathcal{L}^{\star}$,
Eq. (32); on the other hand, a conserved quantity is represented by an
observable in the kernel of $\mathcal{L}^{\star}$. In preparation for
searching for SPT-like physics in QBLs, we will develop in this section a
theory of modes that are ZMs when regarded as symmetries or when regarded as
conserved quantities, but not necessarily both. In this sense, dissipation
“splits” the set of ZMs, and a ZM that is both a conserved quantity and a
symmetry generator occurs only in special circumstances. This suggests that if
one is going to rely on boundary physics as an indicator of SPT-like behavior
in QBLs, then one should consider both kinds of ZMs since there is no clear
reason as yet to privilege edge-localized conserved quantities over symmetry
generators, or the other way around.

Let $\mathcal{L}^{\star}$ be a Heisenberg-picture QBL. Referring to Eq. (26),
it follows that $\widehat{\vec{\alpha}}$ is a ZM if and only if
$\vec{\alpha}\in\ker\widetilde{\mathbf{G}}$. Note that we can, without loss of
generality, take $\widehat{\vec{\alpha}}$ to be Hermitian 444This 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.. We denote the real
vector space of Hermitian ZMs by $\mathcal{Z}$. Likewise, by virtue of Eq.
(26), we have

$\mathcal{L}^{\star}([\widehat{\vec{\alpha}},A])-[\widehat{\vec{\alpha}},\mathcal{L}^{\star}(A)]=[\widehat{-i\mathbf{G}\vec{\alpha}},A],\quad\forall
A.$

Referring to Eq. (32),
$\widehat{\vec{\alpha}}(=\widehat{\vec{\alpha}}^{\dagger})$ is a SG if and
only if $\widehat{-i\mathbf{G}\vec{\alpha}}$ is proportional to the identity.
Since it is linear in $a$ and $a^{\dagger}$, it must be zero. It follows that
$\widehat{\vec{\alpha}}$ generates a symmetry if and only if
$\vec{\alpha}\in\ker\mathbf{G}$. Since the operator
$e^{i\theta\widehat{\vec{\alpha}}}$ has the form of Weyl displacement
operator, we refer to this class of symmetries as Weyl symmetries and their
corresponding generators as Weyl SGs. The corresponding real vector space of
Weyl SGs will be denoted by $\mathcal{W}$. We refer to the both ZMs and Weyl
SGs as Noether modes.

Note that for closed-system evolution ($\mathcal{D}=0$), the dynamical matrix
is $\bm{\tau}_{3}$-pseudo-Hermitian, i.e., $\mathbf{G}=\widetilde{\mathbf{G}}$
555In 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$.. Hence, as expected,
$\mathcal{Z}=\mathcal{W}$. Interestingly, in the purely dissipative case
($H=0)$, the dynamical matrix is skew $\bm{\tau}_{3}$-pseudo-Hermitian, i.e.,
$\mathbf{G}=-\widetilde{\mathbf{G}}$. Again, $\mathcal{Z}=\mathcal{W}$. In
such cases, we say the Noether modes are non-split.

Generically, Noether-modes are split, in the sense that a given mode is either
a ZM or an SG, and not both. Remarkably, however, for QBLs we are able to
establish a direct one-to-one correspondence between ZMs and Weyl SGs. One may
show that, given a fixed ZM, there always exists a corresponding canonically
conjugate Weyl SG. This canonical isomorphism, which we anticipated in Ref.
24, is rather surprising given the general decoupling of ZMs and SGs in open
quantum systems [64]. Formally, we have the following:

Theorem 1. For an arbitrary QBL, we have
$\dim_{\mathbb{R}}\mathcal{Z}=\dim_{\mathbb{R}}\mathcal{W}$. If, in addition,
the zero rapidity hosts only length-one Jordan chains, then for each ZM there
is a canonically conjugate Weyl SG, and vice-versa.

###### Proof.

Consider the antilinear operator $\mathcal{C}$ defined by
$\mathcal{C}\vec{\alpha}=\bm{\tau}_{1}\vec{\alpha}^{*}$ so that
$\widehat{\vec{\alpha}}\,\mbox{}^{\dagger}=-\widehat{\mathcal{C}\vec{\alpha}}$.
In particular, $\widehat{\vec{\alpha}}$ is Hermitian if and only if
$\mathcal{C}\vec{\alpha}=-\vec{\alpha}$. Since
$-\mathbf{G}=\bm{\tau}_{1}\mathbf{G}^{\color[rgb]{0,0,0}*}\bm{\tau}_{1}=\mathcal{C}\mathbf{G}\mathcal{C}$,
it follows that $\ker\mathbf{G}$ is invariant under $\mathcal{C}$. Since
$\mathcal{C}^{2}=\mathds{1}_{2N}$, we have a real structure on
$\ker\mathbf{G}$. That is, $\ker{\mathbf{G}}\cong\mathcal{C}_{-}\\!\oplus
i\mathcal{C}_{-}$, with $\mathcal{C}_{-}$ the real vector space of kernel
vectors $\vec{\alpha}$, with $\mathcal{C}\vec{\alpha}=-\vec{\alpha}$.
“Quantizing” these vectors with the bosonic hat map yields the set of Weyl
SGs. Moreover,
$\dim_{\mathbb{R}}\mathcal{W}=\dim_{\mathbb{R}}\mathcal{C}_{-}=\dim_{\mathbb{C}}\ker\mathbf{G}$.
An identical analysis yields
$\dim_{\mathbb{R}}\mathcal{Z}=\dim_{\mathbb{C}}\ker\widetilde{\mathbf{G}}$.
Elementary linear algebra yields
$\dim_{\mathbb{C}}\ker\mathbf{G}=\dim_{\mathbb{C}}\ker\widetilde{\mathbf{G}}$,
establishing the first claim.

For the second statement, let us take a biorthogonal basis
$\\{\vec{\gamma}^{z}_{j},\vec{\eta}_{j}\\}$ for $\ker\widetilde{\mathbf{G}}$.
By definition, $\\{\vec{\gamma}^{z}_{j}\\}$ spans
$\ker\widetilde{\mathbf{G}}$, $\\{\vec{\eta}_{j}\\}$ spans
$\ker\widetilde{\mathbf{G}}^{\dagger}$, and
$\vec{\eta}_{j}^{\dagger}\vec{\gamma}^{z}_{k}=\delta_{jk}$. The existence of
such a basis hinges upon the length-one Jordan chain assumption. Per the
invariance under $\mathcal{C}$, we can take
$\mathcal{C}\vec{\gamma}_{j}^{z}=-\vec{\gamma}_{j}^{z}$ and
$\mathcal{C}\vec{\eta}_{j}=-\vec{\eta}_{j}$. Now, let
$\vec{\gamma}_{j}^{s}=-i\bm{\tau}_{3}\vec{\eta}_{j}$. These vectors span
$\ker{\mathbf{G}}$ and are odd under $\mathcal{C}$. Finally,

$[\widehat{\vec{\gamma}}_{j}^{s},\widehat{\vec{\gamma}}_{k}^{z}]=\vec{\gamma}_{j}^{s}{}^{\dagger}\bm{\tau}_{3}\vec{\gamma}_{k}^{z}\,1_{\mathcal{F}}=i\vec{\eta}_{j}^{\dagger}\vec{\gamma}^{z}_{k}\,1_{\mathcal{F}}=i\delta_{jk}\,1_{\mathcal{F}}.$

Thus, we have constructed canonically conjugate bases of $\mathcal{Z}$ and
$\mathcal{W}$, as desired. ∎

A relevant generalization of this result establishes a correspondence between
the set of approximate ZMs and the generators of approximate symmetries. We
say that a Hermitian operator $\widehat{\vec{\gamma}}^{z}$ is an approximate
ZM if $\lVert\widetilde{\mathbf{G}}\vec{\gamma}^{z}\rVert<\epsilon$ for a
prescribed accuracy $\epsilon>0$. These arise from
$\epsilon$-pseudoeigenvectors of $\widetilde{\mathbf{G}}$ corresponding to the
$\epsilon$-pseudoeigenvalue $0$. Since, in our work, the most useful norm is
the $2$-norm, $\vec{\gamma}^{z}$ can equivalently characterized as a singular
vector of $\widetilde{\mathbf{G}}$ corresponding to a singular value
$s<\epsilon$ (see also App. B.1). The second-quantized operator then satisfies
$\mathcal{L}^{\star}(\widehat{\vec{\gamma}}{}^{z})=K\alpha$, with
$|K|<\epsilon$ and $\alpha$ some linear form with bounded coefficients.
Generators of approximate Weyl symmetries are defined analogously. See App.
B.2 for a detailed account of this generalization, including an extension of
Theorem 1 to the approximate setting.

###  III.2 Constructive procedures for QBL design

If they are robust against an appropriate set of perturbations, the edge-
localized ZMs of a free-fermion Hamiltonian are likely to be an indication
that the system is in a SPT phase. An SPT phase is a gapped phase of matter
characterized by a degenerate ground state and the absence of a local order
parameter. The ZMs connect the various ground states essentially by changing
the number of the quasiparticles they create or annihilate, with no associated
energy cost. In order to search for this kind of physics in QBLs, one would
greatly benefit from two key capabilities:

  1. 1.

The ability to synthesize QBLs with desirable physical properties (e.g.,
locality) and a rich array of possibilities in terms of their ZMs.

  2. 2.

The ability to synthesize QBLs with pure SSs.

There is a systematic way to meet these needs and we describe them in this
section.

####  III.2.1 A zero-mode-preserving map of Hamiltonian free fermions
to Markovian free bosons

Suppose we have a QFH $H_{F}$ with a midgap state. In the BdG formalism, this
means there is a vector $\vec{\gamma}_{+}$ such that
$\mathbf{H}_{F}\vec{\gamma}_{+}=\epsilon_{N}\vec{\gamma}_{+}$, with
$\epsilon_{N}\sim\mathcal{O}(e^{-N})$ exponentially small in system size and
$\mathbf{H}_{F}$ the BdG Hamiltonian corresponding to $H_{F}$. Because
$\mathcal{F}(\mathbf{H}_{F})=\mathbf{H}_{F}$, there is also a vector
$\vec{\gamma}_{-}$ satisfying
$\mathbf{H}_{F}\vec{\gamma}_{-}=-\epsilon_{N}\vec{\gamma}_{-}$. Moreover, we
can ensure that $\vec{\gamma}_{-}=\mathcal{C}\vec{\gamma}_{+}$ and
$\vec{\gamma}_{+}^{\dagger}\vec{\gamma}_{-}=0$. Upon defining the vector
$\vec{\gamma}\equiv(\vec{\gamma}_{+}-\vec{\gamma}_{-})/\sqrt{2}$, note that

$\lVert\mathbf{H}_{F}\vec{\gamma}\rVert=\epsilon_{N}\frac{\lVert\vec{\gamma}_{+}+\vec{\gamma}_{-}\rVert}{\sqrt{2}}=\epsilon_{N}\sim\mathcal{O}(e^{-N})$
(40)

and, by construction, $\mathcal{C}\vec{\gamma}=-\vec{\gamma}$. If $\Psi$ is
the Nambu array of fermionic operators, then $i\vec{\gamma}{}^{\dagger}\Psi$
is a Majorana fermion whose commutator with $H_{F}$ is exponentially small in
system size.

Now consider the QBL defined by

$\mathbf{H}=0,\qquad\mathbf{M}=\mathbf{H}_{F}+\mathbf{B},$

with $\mathbf{B}=\mathbf{B}^{\dagger}$ any bosonic matrix such that
$\mathbf{H}_{F}+\mathbf{B}\geq 0$. This QBL is purely dissipative (has zero
Hamiltonian) and its dynamical matrix is
$\mathbf{G}=-i\bm{\tau}_{3}\mathbf{H}_{F}$. The main virtue of this
“transmutation” of fermions into bosons is that it preserves ZMs. If
$\mathbf{H}_{F}$ has a midgap state, then we can construct a Hermitian
operator $\gamma\equiv\vec{\gamma}{}^{\dagger}\tau_{3}\Phi$ that (i) is
approximately conserved; and (ii) generates an approximate symmetry. Since
Noether modes are non-split in purely dissipative systems, (i) and (ii) follow
simply from

$\lVert\mathbf{G}\vec{\gamma}\rVert=\lVert\bm{\tau}_{3}\mathbf{H}_{F}\vec{\gamma}\rVert\leq\lVert\mathbf{H}_{F}\vec{\gamma}\rVert\sim\mathcal{O}(e^{-N}).$

Importantly, this operator has the same spatial profile as the Majorana
fermion, and hence is localized on one edge of the chain. The existence of a
second, canonically conjugate Hermitian mode localized on the opposite edge
follows from Theorem 2 in App. B.2.

A potential drawback to our procedure is that the constructed QBL may fail to
be dynamically stable. The reason is that there is no simple relationship
between the spectrum of $\mathbf{H}_{F}$ and $\bm{\tau}_{3}\mathbf{H}_{F}$, in
general. If there is one because $\bm{\tau}_{3}$ commutes with
$\mathbf{H}_{F}$, it means that the QFH commutes with the fermion number
operator. A closely related but nonetheless surprising fact is that a pair of
unitarily equivalent fermionic matrices may lead to QBLs with different
stability properties. Nonetheless, we will also see how useful our map turns
out to be in practice.

####  III.2.2 Designing QBLs with pure steady states via dualities

We now provide a constructive procedure for engineering a QBL that is
guaranteed to relax to the quasiparticle vacuum (QPV) of an arbitrary, but
fixed, dynamically stable QBH. This procedure hinges on the existence of a
particular duality transformation that maps any number-non-conserving QBH to a
number-conserving one [72].

Specifically, consider a dynamically stable QBH $H$, defined by its dynamical
matrix $\mathbf{G}_{0}$. Dynamical stability ensures $\mathbf{G}_{0}$ is
diagonalizable and has an entirely real spectrum. In turn, these features
ensure that there exists a positive-definite matrix $\mathbf{S}$, such that
$\mathbf{G}_{0}^{\dagger}=\mathbf{S}\mathbf{G}_{0}\mathbf{S}^{-1}$. Moreover,
$\mathbf{S}$ and, more importantly, its unique positive-definite square root
$\mathbf{R}=\mathbf{S}^{1/2}$ provide a Bogoliubov transformation
$\Phi\mapsto\mathbf{R}^{-1}\Phi$ that maps $H$ to a number-conserving
Hamiltonian $H^{D}$. The dynamical matrix of $H^{D}$ is
$\mathbf{G}^{D}=\mathbf{R}\mathbf{G}_{0}\mathbf{R}^{-1}$. This transformation
can be understood physically by noting that the bosonic covariance matrix of
the QPV $\ket{\widetilde{0}}$ of $H$ is

$\displaystyle\mathcal{B}(\mathbf{Q}_{\text{QPV}})=\braket{\widetilde{0}}{\mathcal{B}(Q)}{\widetilde{0}}=\frac{1}{2}\mathbf{S}^{-1}=\frac{1}{2}\mathbf{R}^{-2}.$

We claim that the QBL $\mathcal{L}$ defined via the Hamiltonian $H$ and a GKLS
matrix $\mathbf{M}\equiv\kappa(\mathbf{S}+\bm{\tau}_{3})$, with $\kappa>0$ is
(i) dynamically stable; and (ii) relaxes to the unique SS
$\rho_{\text{ss}}=\ket{\widetilde{0}}\bra{\widetilde{0}}$. Firstly, one may
show that this QBL is well-defined, in that $\mathbf{M}\geq 0$. This follows
from the explicit form of $\mathbf{S}$ given in Ref. 72. Dynamical stability
follows from the observation that the dynamical matrix of $\mathcal{L}$ is
$\mathbf{G}_{\kappa}=\mathbf{G}_{0}-i\kappa\mathds{1}_{2N}$. Hence, the
rapidites are simply $-\kappa\pm i\omega_{n}$, with $\omega_{n}$ the
quasiparticle energies of $H$. Dynamically stability ensures that the SS is
unique. Moreover, Gaussianity allows us to characterize it entirely by its
covariance matrix $\mathcal{B}(\mathbf{Q}_{\text{ss}})$. We compute this by
solving the Lyapunov equation Eq. (24) directly. Thus,

$\displaystyle\mathcal{B}(\mathbf{Q}_{\text{ss}})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}e^{-i\mathbf{G}_{\kappa}t}\bm{\tau}_{3}\mathcal{B}(\mathbf{M})\bm{\tau}_{3}e^{i\mathbf{G}_{\kappa}^{\dagger}t}\,{\color[rgb]{0,0,0}dt}.$

By construction, $\mathcal{B}(\mathbf{M})=\kappa\mathbf{S}$. Recalling that
$\mathbf{S}$ provides a Bogoliubov transformation, we have
$\bm{\tau}_{3}\mathbf{S}\bm{\tau_{3}}=\mathbf{S}^{-1}$. Also,

$\displaystyle e^{-i\mathbf{G}_{\kappa}t}=e^{-\kappa t}e^{-i\mathbf{G}_{0}t}.$

With this, we can simplify the integrand by computing

$\displaystyle
e^{-i\mathbf{G}_{0}t}\mathbf{S}^{-1}e^{i\mathbf{G}_{0}^{\dagger}t}$
$\displaystyle=$
$\displaystyle\mathbf{S}^{-1}e^{-i\mathbf{S}\mathbf{G}_{0}\mathbf{S}^{-1}t}e^{i\mathbf{G}_{0}^{\dagger}t}$
$\displaystyle=$
$\displaystyle\mathbf{S}^{-1}e^{-i\mathbf{G}_{0}^{\dagger}t}e^{i\mathbf{G}_{0}^{\dagger}t}=\mathbf{S}^{-1}.$

Altogether, this yields the solution

$\displaystyle\mathcal{B}(\mathbf{Q}_{\text{ss}})$ $\displaystyle=$
$\displaystyle\kappa\mathbf{S}^{-1}\int_{0}^{\infty}e^{-2\kappa
t}\,dt=\frac{1}{2}\mathbf{S}^{-1}.$

Since both the QPV and the SS are Gaussian, have the same mean vector (zero),
and the same covariance matrix, they are the same state. Alternatively, this
may be shown by noting that the dissipator becomes diagonal in a particular
Hamiltonian normal-mode basis (as we will see in the concrete example of Sec.
V.4).

##  IV Dynamical metastability

We have seen thus far that (i) the dynamical matrices of 1D bulk-
translationally invariant QBLs under OBCs are block-Toeplitz matrices; and
(ii) such matrices have non-trivial pseudospectra influenced by the bulk
topology. What are, then, the physical consequences of (ii) given (i)?

The spectrum of a translationally invariant QBL can be described in terms of
curves in the complex plane labeled by the crystal momentum $k$, namely, the
rapidity bands of the system. The pseudospectrum of a bulk-translationally
invariant QBL is, instead, best described in terms of the band winding number.
Consider a QBL with distinguishable rapidity bands, that is, the bands are
described by $M$ complex-valued functions
$\\{\lambda_{m}(k),\lambda_{m}(k)^{*}\\}_{M=1}^{d}$, and let
$\lambda\in{\mathbb{C}}$ not belong to any of the rapidity bands. In the
language of Ref. 73, $\lambda$ is a point-gap of the complex spectrum. For
each band $m=1,\ldots,M$, the winding number of the band about $\lambda$ is
given by

$\nu_{m}(\lambda)=\frac{1}{2\pi
i}\int_{-\pi}^{\pi}\frac{d}{dk}\ln[\lambda_{m}(k)-\lambda]\,dk.$

As we mentioned in Sec. II.3, the SIBC rapidity spectrum consists of the
rapidity bands and all the point-gaps for which $\nu_{m}(\lambda)\neq 0$ for
at least one $m$. Due to the Bloch-like structure of the dynamical matrix
${\mathbf{g}}(k)$ [Eq. (33)], and following the same logic that leads one to
conclude that mid-gap modes are localized in Hermitian systems, the modes
associated to a point-gap are necessarily boundary-localized [67, 65].
Generically, a positive winding number elicits left-localized modes, while a
negative one elicits right-localized modes [30]. Since we adopt the convention
where the semi-infinite chain retains its left boundary, the right localized
modes are lost 666One could, alternatively, consider a chain extending
infinitely to the left, which will retain all right-localized modes and lose
the left-localized ones.. As we also mentioned, the OBC rapidities are more
complicated to describe and less predictable. At a minimum, they are a
discrete subset of the SIBC spectrum and lie on continuous arcs within the
SIBC spectrum. This has especially notable consequences when the OBC spectra
consists entirely of point-gapped rapidities. Such systems turn out to exhibit
the NHSE [66, 67, 68, 35], whereby all finite-size normal modes are edge-
localized.

Strikingly, it is even possible for the OBC chain to be dynamically stable for
all finite $N$, despite having an unstable semi-infinite limit. Physical
intuition stipulates that the dynamics of increasingly large-$N$ truncations
should well-approximate those of the semi-infinite limit. How can we possibly
meld this with the spontaneous loss of dynamical stability, or the spectral
discontinuity in general? The pseudospectrum is precisely the tool for
answering this question. Mathematically, we have seen that the SIBC spectrum
imprints itself into the pseudospectrum of the finite OBC chain, in the
precise sense of Eq. (39). To understand this intuitively, consider a
localized normal mode of the semi-infinite chain, corresponding to a rapidity
far separated from the OBC rapidities for any finite $N$. Truncating this mode
to fit within the finite chain will not yield a normal mode, but rather a
pseudonormal modes with pseudorapidity (pseudoeigenvalue of $-i\mathbf{G}$)
equal to the semi-infinite rapidity. Edge localization implies that early-time
dynamics of this truncated mode cannot sense the existence of the two
boundaries and hence will behave as though the system were semi-infinite in
extent. As $N$ grows, the pseudospectral accuracy parameter $\epsilon$ goes to
zero, and so these pseudonormal modes will be indistinguishable from exact
normal modes over appreciably long timescales.

Making the discussion more concrete, let $\Delta_{S,N}^{\text{OBC}}$ and
$\Delta_{S}^{\text{SIBC}}$ denote the stability gaps of a length-$N$ OBC chain
and the corresponding semi-infinite chain, respectively, as defined in Eq.
(28). Note that $\Delta_{S}^{\text{SIBC}}$ is also the stability gap for
BIBCs. In general, we always have

$\Delta_{S,\infty}^{\text{OBC}}\equiv\lim_{N\to\infty}\Delta_{S,N}^{\text{OBC}}\leq\Delta_{S}^{\text{SIBC}}.$
(41)

However, if the OBC and SIBC spectra differ dramatically, a strict inequality
$\Delta_{S,\infty}^{\text{OBC}}<\Delta_{S}^{\text{SIBC}}$ is possible. We
distinguish two notable cases:

  * (i)

If both $\Delta_{S,\infty}^{\text{OBC}}$ and $\Delta^{\text{SIBC}}_{S}$ are
negative, we say that the OBC chain is in an anomalously relaxing dynamical
phase. This phase, which is dynamically stable both in the finite case and the
infinite-size limit, is characterized by an increasingly long transient time
whereby the system possesses pseudonormal modes decaying at rates much slower
than the rate set by the finite-size stability gap. This transient is followed
by asymptotic decay, whose rate is set by the finite-size stability gap. We
conjecture that a similar mechanism is responsible for the anomalously long
relaxation dynamics found in dissipative systems exhibiting the so-called
Liouvillian skin effect [36, 37, 75]. This dynamical behavior also bears some
resemblance to metastability associated with spectral separation [33], with
the spectral separation here being between the finite-size stability gap, and
the infinite-size (or pseudospectral) gap .

  * (ii)

If $\Delta_{S,\infty}^{\text{OBC}}<0<\Delta^{\text{SIBC}}_{S}$, we say that
the OBC chain is dynamically metastable (or just metastable, in our context,
for short). This phase is characterized by an increasingly long transient time
whereby the pseudo-normal modes whose pseudoeigenvalues are in the right half
plane amplify exponentially. Once this transient concludes, all modes decay
asymptotically with rate set by the finite-size stability gap. The terminology
“dynamical metastability” refers to the metastable amplifying transient
dynamical phase that eventually gives way to the necessarily stable asymptotic
dynamics. Note that, since our notion aims to capture metastable behavior of
the structural (stability) properties of the dynamical generator, it is, at
the current stage of understanding, a distinct notion than metastability for
states of open quantum systems that is known to stem from large intra-spectral
gaps [33, 76, 34]. However, both cases are characterized by delayed relaxation
to the true SS manifold.

Figure 2: Anomalously relaxing vs. dynamically metastable dynamics in the DBKC
(Sec. V.3). The solid lines are the trajectory of $|\braket{x_{N}}(t)|$
averaged over 250 random initial conditions. The filled regions are $\pm$ one
standard deviation from the mean. The black dashed lines are the dynamics
predicted from the SIBC stability $\Delta_{S}^{\text{SIBC}}$ gap, while the
gray dashed lines are the dynamics expected from the finite-size stability
gaps,
$\Delta_{S,\infty}^{\text{OBC}}=\lim_{N\to\infty}\Delta_{S,N}^{\text{OBC}}$.
For this example,
$\Delta_{S,N}^{\text{OBC}}=\Delta_{S,\infty}^{\text{OBC}}=-\kappa$. In all
cases, $J=2$, $\Delta=0.5$, $\mu=\Gamma=0$, $N=25$. The dynamically metastable
(anomalously relaxing) curve corresponds to $\kappa/\Delta=0.6$ (1.4).

Examples of the above dynamical behavior are found in the dissipative bosonic
Kitaev chain (DBKC) model whose rapidity spectra are shown in Fig. 1 (see Sec.
V.3 for a detailed discussion). In particular, the rapidities corresponding to
open marker in Fig. 1(a) show an anomalously relaxing phase when the
dissipation rate $\kappa$ exceeds the two-photon pumping, $\Delta$, while the
filled (darker) rapidities show a metastable phase. To explicitly illustrate
the two-step nature of the relaxation dynamics that classes (i) and (ii) above
support, in Fig. 2 we plot the trajectory of $|\braket{x_{N}}(t)|$ averaged
over an ensemble of random initial conditions. In the anomalously relaxing
phase, exponential relaxation indeed proceeds in two steps. The transient
decay rate is set by $\Delta_{S}^{\text{SIBC}}<0$, while the asymptotic decay
rate is the true one, i.e., $\Delta_{S,N}^{\text{OBC}}<0$. The dynamically
metastable phase shows a drastically different transient behavior, with
amplification engendered by $\Delta_{S}^{\text{SIBC}}>0$. Dynamical stability
nonetheless guarantees that the asymptotic dynamics coincide with those of the
anomalously relaxing phase.

It is natural to ask how long these anomalous transient regimes persist. For
this, consider the normal mode of the SIBC system corresponding to the
rapidity with the largest real part. In the anomalously relaxing case, this is
the slowest relaxing mode, while in the dynamically metastable case, it is the
one that amplifies the fastest. Truncating this mode to fit into a chain of
length $N$ will provide a pseudonormal mode of the OBC system of accuracy
$\epsilon(N)$. As before, the well-behaved nature of the pseudospectrum
infinite-size limit implies that $\epsilon(N)\to 0$ as $N\to\infty$. Referring
to Eq. (38), the lifetime of this pseudonormal mode will necessarily increase
as $\epsilon(N)\to 0$. Thus, the transient timescale diverges as $N\to\infty$.
One consequence of this fact is that, in the dynamically metastable case, the
linear mixing time of the QBL (roughly, the time it takes the first moments of
arbitrary states to come suitably close to their SS value) will diverge [24].

In what follows, we shall focus on showing how dynamically metastable phases
can be further categorized into those that are topological in an appropriate
sense (as in Fig. 1(a)), and those that are not (as in Fig. 1(b)). We further
remark that a more comprehensive analysis of dynamical metastability will be
presented in Ref. 77.

##  V Topological dynamical metastability with broken number symmetry

Model | Hamiltonian | Lindblad dissipator | U(1) symmetry | Noether modes | Steady state
---|---|---|---|---|---
DBKC | BKC + DPA | Uniform onsite plus next-NN damping | No | Split | Mixed
Pure-SS | BKC | Uniform damping in BKC normal-mode basis; | No | Split | Pure
DBKC |  | site-local; restricted translation symmetry |  |  |
PDMC | $H=0$ | Uniform onsite and NN damping, pumping, pairing | No | Non-split | Mixed
DDWC | DPA+ NDPA | Uniform onsite damping and pumping, NN pairing | No | Split, doubled | Mixed
DNSC | (i)TB | Uniform onsite and NN damping, onsite pumping | Yes | Split | Mixed
Table 3: Summary of the key features of the five topologically metastable
models analyzed in this work. Abbreviations: DBKC = dissipative bosonic Kitaev
chain; PDMC = purely dissipative Majorana chain; DDWC = dissipative double-
winding chain; DNSC = dissipative number-symmetric chain; DPA = degenerate
parametric amplifier; NDPA = non-degenerate parametric amplifier; NN =
nearest-neighbor; (i)TB = (imaginary) tight-binding. In all models (except the
pure SS DBKC), the purity of the SS is assessed numerically throughout the
relevant regions of parameter space (data not shown). The DDWC, which is found
to support two ZMs (SGs) on each edge, is discussed in App. D, while the rest
are in the main text.

Let us now restrict our focus on dynamically metastable QBLs, and further
assume that they (i) are point-gapped at zero (i.e., their Bloch dynamical
matrix $\mathbf{g}(k)$ is invertible for all $k$); and (ii) have at least one
rapidity band winding about zero. We shall call QBLs belonging to this class
topologically dynamically metastable (or just topologically metastable, for
short). Physically, an additional important restriction stems from the
possible presence of a (weak) U(1) number symmetry, generated by the total
bosonic number operator. While we defer the discussion of number-symmetric
QBLs to Sec.VI, we examine first QBLs which, in analogy to fermionic
topological superconducting phases, will feature some bosonic pairing
mechanisms, at either the Hamiltonian or dissipative level. For reference, a
summary of all the illustrative models we will analyze in this work is
provided in Table 3.

###  V.1 Majorana edge bosons: General properties

From the definition of topological metastability, in addition to the spectral
and pseudospectral properties of Toeplitz and Laurent matrices and operators
discussed in Sec. II.3, we may immediately conclude that
$\mathbf{G}_{N}^{\text{OBC}}$ possesses at least one pseudonormal mode with
zero pseudorapidity. Similarly, we can conclude that
$\widetilde{\mathbf{G}}_{N}^{\text{OBC}}$ also possesses a pseudonormal mode
with zero pseudorapidity. Following the conventions of Sec. III.1, we call
these two pseudomodes $\vec{\gamma}^{s}$ and $\vec{\gamma}^{z}$, respectively,
and ensure that the associated linear forms
$\gamma^{s}\equiv\widehat{\vec{\gamma}}{}^{s}$ and
$\gamma{}^{z}\equiv\widehat{\vec{\gamma}}^{z}$ are Hermitian, and have
commutator777Note that the canonical correspondence assumes the technical
assumptions of Theorem 2 in App. B.2 are met. equal to $i$. We deem the pair
$(\gamma^{z},\gamma^{s})$ Majorana bosons. MBs enjoy a number of remarkable
properties derived from their pseudospectral and topological origins:

  * (i)

An MB pair consists of one approximate ZM $\gamma^{z}$ and one generator of an
approximate Weyl symmetry $\gamma^{s}$. Both are necessarily Hermitian. That
is, MBs are a particular instance of Noether modes.

  * (ii)

The pair can always be normalized to satisfy canonical commutation relations,
while the roles of $\gamma^{z}$ and $\gamma^{s}$ as approximate ZM and SG are
maintained.

  * (iii)

One member of the pair is exponentially localized on the left half of the
chain, while the other is localized on the right. This follows because, due to
the adjoint relationship between $\widetilde{\mathbf{G}}$ and $\mathbf{G}$,
the winding associated to $\gamma^{z}$ and $\gamma^{s}$ necessarily have
opposite sign. This engenders the stated localization properties [30, 31, 79].

  * (iv)

Combining (i)-(iii) allows us to construct a spatially split bosonic degree of
freedom $\gamma^{z}+i\gamma^{s}$ whose quadrature components are the MBs. In
the case where the MBs are non-split (in the sense of Sec. III.1), this
creates a long-lived bosonic excitation in the system.

It is worth to zoom in specifically on point (ii). While a more general and
rigorous discussion will be provided in a mathematics-focused companion paper
[79], a typical instance of MBs takes on the form:

$\gamma^{z}\equiv\mathcal{M}_{z}(N)\sum_{j=1}^{N}\delta^{j}x_{j},\quad\gamma^{s}\equiv\mathcal{M}_{s}(N)\sum_{j=1}^{N}\delta^{N-j}p_{j},$
(42)

with $\mathcal{M}_{z,s}(N)$ size-dependent normalization constants to be
determined and $\delta$ real, with $|\delta|<1$. The above operators satisfy

$\displaystyle\mathcal{L}^{\star}(\gamma^{z})$ $\displaystyle=$
$\displaystyle\mathcal{M}_{z}(N)\delta^{N-1}\chi,$ (43)
$\displaystyle\mathcal{L}^{\star}([\gamma^{s},A])-[\gamma^{s},\mathcal{L}^{\star}(A)]$
$\displaystyle=$ $\displaystyle\mathcal{M}_{s}(N)\delta^{N-1}[\xi,A],\quad$
(44)

along with the algebraic relationship

$[\gamma^{z},\gamma^{s}]=i\mathcal{M}_{z}(N)\mathcal{M}_{s}(N)N\delta^{N-1}.$
(45)

Here, $\chi$ and $\xi$ represent (typically localized) Hermitian linear forms
whose coefficients in the Nambu basis are system-size independent. The goal of
normalization is to choose $\mathcal{M}_{z}(N)$ and $\mathcal{M}_{s}(N)$ such
that: (i) the right hand-sides of Eqs. (43) and (44) go to zero as
$N\to\infty$; and (ii) the right hand-side of Eq. (45) is $i$. That is, we
want canonically conjugate modes that provide asymptotically exact ZMs and
SGs. Mathematically, these conditions are

(i)
$\displaystyle\lim_{N\to\infty}\mathcal{M}_{z}(N)\delta^{N-1}=\lim_{N\to\infty}\mathcal{M}_{s}(N)\delta^{N-1}=0,$
(ii) $\displaystyle\mathcal{M}_{z}(N)\mathcal{M}_{s}(N)N\delta^{N-1}=1.$

Remarkably, a scheme satisfying both (i) and (ii) always exists [79]. The
simplest, and most natural choice is

$\mathcal{M}_{z}(N)=\mathcal{M}_{s}(N)=\frac{\delta^{-(N-1)/2}}{\sqrt{N}}.$

With this, condition (ii) is clearly satisfied. As for (i),

$\lim_{N\to\infty}\mathcal{M}_{z}(N)\delta^{N-1}=\lim_{N\to\infty}\frac{\delta^{(N-1)/2}}{\sqrt{N}}=0,$

and similarly for $\mathcal{M}_{s}(N)$. This scheme, which we call symmetric
normalization, will be employed throughout the text. Clearly, this choice is
non-unique, and another may be adopted if certain features are desired over
others. For example, if we wish for $\gamma^{z}$ to have bounded coefficients
in the $x_{j}$ basis, then we may take $\mathcal{M}_{z}(N)=1$ and
$\mathcal{M}_{s}(N)=\delta^{-(N-1)}/N$. Crucially, the exponential profile of
the MBs is independent of the choice of normalization. We also remark that
this approach to MB normalization shares many features with the fermionic
case, as detailed in App. A.

While the dynamical features of the approximate ZM MB $\gamma^{z}$ is
relatively easy to describe, those of $\gamma^{s}$ are initially more opaque.
We illuminate one particular implication of the existence of this approximate
symmetry by leveraging the uniqueness of the SS under OBCs. Denoting said SS
by $\rho_{\text{ss}}$, we introduce the family of Weyl-displaced Gaussian
states

$\rho_{\theta}=e^{i\theta\gamma^{s}}\rho_{\text{ss}}e^{-i\theta\gamma^{s}},\quad\theta\in\mathbb{R}.$

Because $\gamma^{s}$ generates an approximate symmetry, the states
$\rho(\theta)$ are quasi-steady, in the sense that $\dot{\rho}_{\theta}(0)\sim
0+\mathcal{O}(\theta\epsilon)$, with $\epsilon$ being an exponentially small
accuracy parameter such that $0$ is in the $\epsilon$-pseudospectrum of
$\mathbf{G}$. Unlike the SSs of a QBL, these quasi-SSs possess non-zero mean
vectors:

$\vec{m}_{\theta}(0)\equiv\text{tr}[\Phi\rho_{\theta}(0)]=i\theta\vec{\gamma}^{s}.$

Since $\vec{\gamma}^{s}$ is edge-localized and a pseudonormal mode of the
dynamical matrix, the quasi-steady mean vectors are exponentially localized on
one edge of the chain and are long-lived. More explicitly, we may write

$\frac{\lVert\vec{m}_{\theta}(t)-\vec{m}_{\theta}(0)\rVert}{\lVert\vec{m}_{\theta}(0)\rVert}\leq\epsilon\,t\,e^{\Omega
t},$ (46)

with $\Omega>0$ a system-size independent constant and $\epsilon\rightarrow 0$
as $N\rightarrow\infty$ (see Ref. 24 for an explicit derivation). On the other
hand, being only a Weyl displacement of the SS, it follows that
$\rho_{\theta}$ shares a covariance matrix (and hence all even moments) with
$\rho_{\text{ss}}$. Interestingly, these states can be used to construct long-
lived classical non-Gaussian states. Namely, any convex linear combination of
the $\rho_{\theta}$’s will be long-lived and exhibiting a non-Gaussian, albeit
still non-negative, Wigner function.

###  V.2 A purely dissipative bosonic Majorana chain

As the first concrete example of topological metastability, we focus on a
purely dissipative ($H=0$) setting. The motivation for this is twofold.
Firstly, the existence of a purely dissipative model displaying topological
metastability demonstrates that the phenomena does not require any Hamiltonian
contribution. In particular, this establishes that topological metastability
is not a property of the Hamiltonian that is stabilized through dissipation
but, rather, a fundamentally dissipative phenomenon. In fact, a similar
philosophy underpins early analysis of topological physics in quasi-free open
fermionic systems [26]. Secondly, as detailed in Sec. III.1, the MBs of purely
dissipative chain are necessarily non-split: Each MB is both a ZM and an SG.
Thus, they provide the most natural bosonic generalization of the Majorana
fermions of topological superconductivity.

Figure 3: Comparison between the topological phase diagram of the FKC and the
topological stability phase diagram of the corresponding PDMC.

To accomplish this, let us apply the construction of Sec. III.2.1 to the
paradigmatic fermionic Kitaev chain (FKC) (see also App. A), whose phase
diagram is shown in Fig. 3. Within the topological phase and under OBCs,
$H_{\text{FKC}}$ has two approximate Majorana ZMs (in the special case
$\mu=0$, these approximate ZMs are exact for odd $N$). The corresponding BdG
Hamiltonian is
$\mathbf{H}_{\text{FKC}}=\mathds{1}_{N}\otimes\mathbf{h}_{0}+\mathbf{S}\otimes\mathbf{h}_{1}+\mathbf{S}^{\dagger}\otimes\mathbf{h}_{1}^{\dagger}$,
with the internal matrices given by

$\displaystyle\mathbf{h}_{0}=\begin{bmatrix}-\mu&0\\\
0&\mu\end{bmatrix},\quad\mathbf{h}_{1}=\begin{bmatrix}-J&\Delta\\\
-\Delta&J\end{bmatrix}.$

Moving to the corresponding QBL, we set
$\mathbf{M}\equiv\mathbf{H}_{\text{FKC}}+\alpha\mathds{1}_{2N}$, with
$\alpha\geq|\min\sigma(\mathbf{H}_{\text{FKC}})|$, so that we ensure
$\mathbf{M}\geq 0$. Physically, the Lindblad dissipator describing this purely
dissipative Majorana chain (PDMC) has five distinct incoherent contributions:
$\mathcal{D}\equiv\mathcal{D}_{-,0}+\mathcal{D}_{+,0}+\mathcal{D}_{-,1}+\mathcal{D}_{+,1}+\mathcal{D}_{p,1}$,
with

$\displaystyle\mathcal{D}_{-,0}$ $\displaystyle=$
$\displaystyle(\alpha-\mu)\sum_{j=1}^{N}\mathcal{D}[a_{j},a_{j}^{\dagger}],$
$\displaystyle\mathcal{D}_{+,0}$ $\displaystyle=$
$\displaystyle(\alpha+\mu)\sum_{j=1}^{N}\mathcal{D}[a_{j}^{\dagger},a_{j}],$
$\displaystyle\mathcal{D}_{-,1}$ $\displaystyle=$
$\displaystyle-J\sum_{j=1}^{N}\mathcal{D}[a_{j},a_{j+1}^{\dagger}]+\mathcal{D}[a_{j+1},a_{j}^{\dagger}],$
$\displaystyle\mathcal{D}_{+,1}$ $\displaystyle=$ $\displaystyle
J\sum_{j=1}^{N}\mathcal{D}[a_{j}^{\dagger},a_{j+1}]+\mathcal{D}[a_{j+1}^{\dagger},a_{j}],$
$\displaystyle\mathcal{D}_{p,1}$ $\displaystyle=$
$\displaystyle\frac{\Delta}{2}\sum_{j=1}^{N}\mathcal{D}[a_{j},a_{j+1}]-\mathcal{D}[a_{j+1},a_{j}]-(a\leftrightarrow
a^{\dagger}).$

From top to bottom, these mechanisms can be described as uniform onsite
damping of strength $\alpha-\mu$, uniform onsite pumping of strength
$\alpha+\mu$, uniform NN damping of strength $J$, uniform NN pumping of
strength $J$, and incoherent pairing of strength $\Delta$, respectively.

Diagonalizing the GKLS matrix leads to a particularly interesting set of
Lindblad operators. Let $\\{\vec{\psi}_{n,\pm}\\}$ be the eigenvectors of
$\mathbf{H}_{\text{FKC}}$ (and hence, $\mathbf{M}$) with eigenvalues
$\\{\pm\epsilon_{n}\\}$. These represent the FKC quasiparticles ($+$) and
quasiholes ($-$). The Lindblad operators are then
$L_{n}^{\pm}=\sqrt{\alpha\pm\epsilon_{n}}\,\vec{\psi}^{\,{\dagger}}_{n,\pm}\Phi$.
Under PBCs, these correspond to delocalized plane-waves with $\pm\epsilon_{n}$
in the FKC bands. The OBC case depends on the topological phase. In the
trivial phase, the $L_{n}^{\pm}$’s are delocalized standing waves with
$\pm\epsilon_{n}$ in the bands. In the non-trivial phase, there are $2N-2$
standing waves with $\pm\epsilon_{n}$ confined to the bands and 2 edge modes
with $\pm\epsilon_{n}\sim\pm e^{-N}$ in the gap. These two modes provide
boundary dissipation and pumping for the corresponding bosonic chain.

The Bloch dynamical matrix of the PDMC is founds as

$-i\mathbf{g}(k)=(\mu+2J\cos(k))\mathds{1}_{2}-2i\Delta\sin(k)\bm{\sigma}_{1}.$

Thus, the rapidity bands for this model are given by

$\lambda_{\pm}(k)=\mu+2J\cos(k)\pm 2i\Delta\sin(k).$

Both dynamical instability and non-zero winding are guaranteed whenever
$|\mu/2J|<1$, which corresponds exactly to the non-trivial phase of the FKC.
However, the OBC chain is unstable whenever $\mu>0$, see Fig. 3. Explicitly,
the OBC rapidities for the PDMC are given by

$\lambda_{m}=\mu+2i\sqrt{\Delta^{2}-J^{2}}\cos\left(\frac{m\pi}{N+1}\right),$
(47)

with $m=0,\ldots,2N-1$ (see App. B.3 for a more detailed spectral analysis).
For simplicity, we take $\Delta\geq J>0$ and restrict $\mu<0$. The relevant
stability gaps are then

$\Delta_{S,N}^{\text{OBC}}=\mu,\quad\Delta_{S}^{\text{SIBC}}=2J+\mu.$

Thus, we conclude that the PDMC (under OBCs) is topologically metastable
whenever $\mu<0$ and $-\mu/(2|J|)<1$.

For general parameter values, the calculation of MBs is performed numerically,
as described in more detail in App. B.4. However, we may construct the
(unnormalized) MBs analytically in the special case $J=\Delta$:

$\gamma_{L}=\sum_{j=1}^{N}\delta^{j-1}x_{j},\quad\gamma_{R}=\sum_{j=1}^{N}\delta^{N-j}p_{j},\quad\delta\equiv-\frac{\mu}{2J}.$
(48)

We immediately see that (i) both $\gamma_{L}$ and $\gamma_{R}$ are approximate
ZMs and generate approximate symmetries; and (ii) they have the exact same
localization lengths as the MF edge modes of the FKC. While (ii) can be
inferred directly from Eq. (48), (i) is verified by the relations

$\displaystyle\mathcal{L}^{\star}(\gamma_{L})$ $\displaystyle=$
$\displaystyle-2J\delta^{N}x_{N},$
$\displaystyle\mathcal{L}^{\star}(\gamma_{R})$ $\displaystyle=$
$\displaystyle-2J\delta^{N}p_{1},$
$\displaystyle\mathcal{L}^{\star}([\gamma_{L},A])-[\gamma_{L},\mathcal{L}^{\star}(A)]$
$\displaystyle=$ $\displaystyle-2iJ\delta^{N}[x_{N},A],\,\,\forall A,$
$\displaystyle\mathcal{L}^{\star}([\gamma_{R},A])-[\gamma_{R},\mathcal{L}^{\star}(A)]$
$\displaystyle=$ $\displaystyle-2iJ\delta^{N}[p_{1},A],\,\,\forall A.$

Moreover, the unnormalized MB algebra is
$[\gamma_{L},\gamma_{R}]=2iN\delta^{N-1}$. Normalization proceeds exactly as
detailed in Sec. V.1. Finally, as anticipated, the MBs are non-split and have
the same spatial profile as the Majorana fermions of the FKC [6]. In Sec. VII,
we will show how non-split MBs may be distinguished from the more typical
split MBs via two-time correlation functions.

###  V.3 A dissipative bosonic Kitaev chain

Figure 4: (a) The topological phase diagram of the DBKC with $J\geq\Delta>0$
and $\mu=0$. The “Non-zero winding” region indicates the parameter regime
where the rapidity bands wind around the origin. The “Zero winding” region
indicates the parameter regime where neither band winds around the origin. If
the OBC chain is dynamically stable in either of these two regions, then it is
dynamically metastable, additionally being topologically metastable in the
former case. The “Dynamically stable bulk” region occurs when the rapidity
bands lie in the left-half plane. If the OBC chain is dynamically stable here,
then it is anomalously relaxing. In the “Ill-defined” region, $\mathbf{M}$ is
no longer positive-semidefinite. The black dots at
$(\Gamma/\Delta,\kappa/\Delta)=(0,0.6)$ and $(0.24,0.6)$ indicate the position
of two representative parameter choices used in later figures. (b) For later
reference, the SS purity under OBCs in the same phase space as (a) with $J=2$,
$\Delta=0.5$, and $N=25$. (c) The topological stability phase diagram of the
DBKC under OBCs with $\Gamma=0$. The inset figures show representative
rapidity band structure in each region. The number of bands winding around the
origin determines the number of MB pairs. (a) and (c) adapted from Ref. [24].

The MBs of the previous model represent the simplest bosonic extension of
Majorana fermions. However, the model itself is atypical in the sense that
QBLs generically possess both coherent (Hamiltonian) and incoherent processes.
As a quintessential example of topological metastability furnished by the
interplay between these two dynamical contributions, we consider the
dissipative BKC (DBKC) we introduced in Ref. 24. The Hamiltonian component now
consists of the BKC Hamiltonian [32], modified with an isotropic degenerate
parametric amplifier (DPA) term at each site. Explicitly,

$H_{\text{BKC}}=\frac{i}{2}\sum_{j=1}^{N}\\!\left(Ja_{j+1}^{\dagger}a_{j}+\Delta
a_{j+1}^{\dagger}a_{j}^{\dagger}+\mu
a_{j}^{\dagger}{}^{2}\right)+\text{H.c.},$ (49)

with $J\geq\Delta\geq 0$ being the hopping and pairing amplitudes,
$\mu\in\mathbb{R}$ the uniform DPA strength, and $a_{N+j}\equiv a_{j}$ for
PBCs, or $a_{N+j}\equiv 0$ for OBCs. By making the system open, we introduce
two damping mechanisms:

$\mathcal{D}=\sum_{j=1}^{N}\bigg{(}\kappa-\frac{(2\Gamma)^{2}}{\kappa}\bigg{)}\mathcal{D}[a_{j}]+\mathcal{D}\Big{[}\sqrt{\kappa}a_{j}+\frac{2\Gamma}{\sqrt{\kappa}}a_{j+2}\Big{]},$

where $\kappa\geq 0$ is a uniform onsite damping rate and $\Gamma\geq 0$ is a
next-nearest neighbor (NNN) damping rate. BCs are enforced in the usual way,
whereas we ensure the GKLS matrix $\mathbf{M}$ is positive-semidefinite by
requiring $\kappa\geq 2\Gamma$.

The Bloch dynamical matrix is given by

$\mathbf{g}(k)\\!=\\!-i[\kappa-
iJ\sin(k)+2\Gamma\cos(2k)]\mathds{1}_{2}+i[\mu+\Delta\cos(k)]\bm{\sigma}_{1},$

from which the two rapidity bands are determined to be

$\lambda_{\pm}(k)=-(\kappa\pm\mu)\mp\Delta\cos(k)-iJ\sin(k)-2\Gamma\cos(2k),$
(50)

where $k\in[-\pi,\pi]$ is the crystal momentum. The stability phase diagram
for PBCs/BIBCs is obtained by determining when the bands cross the imaginary
axis. In the unstable phase, two possibilities arise: (i) at least one band
winds around the origin or (ii) neither band wind around the origin. A
representative stability phase diagram, including the partition of the
unstable phase into topological and non-topological regions, is shown in Fig.
4(a).

####  V.3.1 The parameter regime $\Gamma=0$

The parameter regime where the NNN damping vanishes is especially interesting
from the perspective of topological dynamical metastability. In this case, the
OBC rapidities can be computed as

$\lambda_{\pm,m}=-\kappa\pm\mu+i\sqrt{J^{2}-\Delta^{2}}\cos\left(\frac{m\pi}{N+1}\right),$

with $m=0,\ldots,2N-1$. In the notation of Sec. IV, the relevant stability
gaps in this limit are

$\Delta_{S,N}^{\text{OBC}}=-\kappa+|\mu|,\quad\Delta_{S}^{\text{SIBC}}=\Delta-\kappa+|\mu|.$

Example rapidities for $\Gamma=\mu=0$ are shown in Fig. 1(a). Stability is
guaranteed whenever $\kappa\geq|\mu|$, independent of $N$. Under this
constraint, the winding numbers $\nu_{\pm}$ of the bands $\lambda_{\pm}(k)$
can be tuned to be either $(\nu_{+},\nu_{-})=(1,-1)$, $(1,0)$, $(0,-1)$, or
$(0,0)$, see Fig. 4(c). The first three regimes exhibit topological
metastability while the fourth is anomalously relaxing. Each non-zero winding
number corresponds to an MB pair. Thus, in a topologically metastable regime,
the number of MB pairs can either be 1 or 2: in the former case, a particular
edge hosts either a boundary-localized ZM or a SG, whereas in the latter case
both a boundary-localized ZM and a SG are present on each edge. We remark that
the existence of phases where $\sum\nu\neq 0$ explicitly requires dissipation,
and is interesting from the perspective of Fredholm theory [79].

For the special case $J=\Delta$, we compute the MBs exactly:

$\displaystyle\gamma_{L}^{z}$
$\displaystyle=\sum_{j=1}^{N}\delta_{-}^{j-1}x_{j},$
$\displaystyle\gamma_{R}^{s}$
$\displaystyle=\sum_{j=1}^{N}\delta_{-}^{N-j}p_{j},$
$\displaystyle\gamma_{L}^{s}$
$\displaystyle=\sum_{j=1}^{N}\delta_{+}^{j-1}x_{j},$
$\displaystyle\gamma_{R}^{z}$
$\displaystyle=\sum_{j=1}^{N}\delta_{+}^{N-j}p_{j},$

with $\delta_{\pm}=-(\mu\pm\kappa)/J$. The relevant non-vanishing commutators
are $[\gamma_{L}^{z},\gamma_{R}^{s}]=iN\delta_{-}^{N-1}$ and
$[\gamma_{L}^{s},\gamma_{R}^{z}]=iN\delta_{+}^{N-1}$. The pseudospectral
origin of each mode and its role as approximate ZM / SG are evident in the
following identities:

$\displaystyle\mathcal{L}^{\star}(\gamma_{L}^{z})$
$\displaystyle=-J\delta_{-}^{N}x_{N},$
$\displaystyle\mathcal{L}^{\star}(\gamma_{R}^{z})$
$\displaystyle=J\delta_{+}^{N}p_{1},$
$\displaystyle\mathcal{L}^{\star}([\gamma_{L}^{s},A])-[\gamma_{L}^{s},\mathcal{L}^{\star}(A)]$
$\displaystyle=-J\delta_{+}^{N}[x_{N},A],\,\,\forall A,$
$\displaystyle\mathcal{L}^{\star}([\gamma_{R}^{s},A])-[\gamma_{R}^{s},\mathcal{L}^{\star}(A)]$
$\displaystyle=J\delta_{-}^{N}[p_{1},A],\,\,\forall A.$

Whenever $\delta_{+}$ or $\delta_{-}$ has modulus less than one, the
corresponding $c/s$ mode is an approximate ZM / SG. Independence of
$\delta_{+}$ and $\delta_{-}$ leads to the phase diagram of Fig. 4(c). It is
interesting to note that, unlike the exact normal-modes of the system (which
are pinned to the edges via the NHSE), the localization lengths of the MBs,
$\xi_{\pm}\equiv[\ln|\delta_{\pm}|]^{-1}$, diverge at the topological phase
transition. In fact, divergence of the localization length explicitly define
the phase boundaries, reminiscent of Majorana fermions in topological
superconductors.

####  V.3.2 The parameter regime $\Gamma\neq 0$

When the NNN damping is turned on and exceeds a certain threshold, the DBKC
can exhibit non-topological metastability. The case where $\mu=0$ is
illustrated in both Fig. 4(a) and Fig. 1(b). In this case, Eq. (50) reveals
that the two rapidity bands are degenerate; with respect to the elliptical
bands of the BKC Hamiltonian, the onsite damping introduces an overall shift,
while the NNN damping introduces more complex curvature. The bulk stability
gap is found as

$\displaystyle\Delta_{S}^{\text{BIBC}}=\Delta_{S}^{\text{SIBC}}=\begin{cases}-\kappa+\Delta-2\Gamma,&\Gamma/\Delta<1/8,\\\
-\kappa+\Delta^{2}/16\Gamma+2\Gamma,&\Gamma/\Delta>1/8,\end{cases}$

with the region where $\Delta_{S}^{\text{BIBC}}>0$ being divided into a sector
where both bands wind around the origin and one where neither do. The latter
corresponds to a non-topological dynamically metastable phase, with no MBs.
Nonetheless, transient amplification is observed [24]. As explained in Sec.
IV, this follows due to the fact that non-trivial pseudospectra (equivalently,
the infinite-size limit spectra) emerges in the right-half complex plane. In
Sec. VII we will see that this non-topological metastability can be
distinguished from its topological counterpart using two-time correlation
functions.

###  V.4 A dissipative bosonic Kitaev chain with a pure steady state

In all models considered thus far, we have confirmed by numerically
determining the purity that the SS of the OBC configuration is always a mixed
state, see Fig. 4(b). This prevents us from exploring the potential interplay
between topological metastability and pure SSs (so-called “dark states” in
quantum optics). Following the procedure detailed in Sec. III.2.2, however, it
is always possible, given a dynamically stable Hamiltonian with QPV
$\ket{\widetilde{0}}$, to engineer a dissipator that relaxes any initial
condition to $\ket{\widetilde{0}}$. In the spirit of dissipative pure-state
stabilization [50, 80], such a procedure is especially interesting when such a
vacuum state possesses certain nontrivial properties such as non-zero
squeezing, as is the case with the BKC. Can topological metastability arise in
such a system? The answer is Yes, as we will now show by applying the
procedure to the BKC Hamiltonian under OBCs.

Consider the Lindbladian defined via the BKC Hamiltonian (with $\mu=0$ for
simplicity) and the site-local dissipator

$\mathcal{D}=2\kappa\sum_{j=1}^{N}\mathcal{D}[\alpha_{j}(r)],\quad\kappa>0,$
(51)

where we have introduced squeezed bosonic degrees of freedom according to

$\alpha_{j}(r)\equiv\cosh(jr)a_{j}-\sinh(jr)a_{j}^{\dagger},$

and the squeezing parameter $r$ is fixed according to $\tanh(r)=\Delta/J$.
This represents a local dissipative process of the squeezed degrees of freedom
with constant loss rate $2\kappa$. A number of remarks are in order.

First, the dissipator is not bulk-translationally invariant. Explicitly, the
squeezing of the mode $\alpha_{j}(r)$ increases with $j$. Remarkably, however,
the dynamical matrix
$\mathbf{G}=\mathbf{G}\big{|}_{\kappa=0}-i\kappa\mathds{1}_{2N}$ is
translationally invariant, up to BCs. In fact, it is precisely equivalent to
that of the DBKC of Sec. V.3. Mathematically, the translation invariance-
breaking terms are entirely contained in the bosonic projection of the GKLS
matrix $\mathbf{M}$, while $\mathscr{F}(\mathbf{M})=\kappa\tau_{3}$. Thus, the
Heisenberg equations of motion for the Nambu array, and thus all linear
observables, are translationally invariant, up to BCs. In this sense, we say
this model has a restricted translational symmetry. This allows for a
straightforward computation of rapidities and pseudospectra. Since the
dynamical matrix is equivalent to that of the DBKC, the chain is topologically
metastable for $|\kappa/\Delta|<1$ and anomalously relaxing for
$|\kappa/\Delta|>1$. In particular, there are two MB pairs in the
topologically metastable regime.

Second, the dissipator has a particularly simple representation in terms of a
set of normal modes of the BKC Hamiltonian. Recall that the set of normal
modes of the BKC is non-unique (as a consequence of anti-commutation with