The discovery of phases of matter exhibiting topological order revolutionized our understanding of many-body quantum systems, by challenging the Ginzburg-Landau theory of local order parameters. Symmetry-protected topological (SPT) phases of free (mean-field) bulk-gapped fermionic matter in equilibrium have been fully classified [1]. A pillar of the ensuing bulk-boundary correspondence is the emergence of zeroenergy boundary-localized Majorana fermions (MFs) that are robust against symmetry-preserving local perturbations and, thanks to the non-Abelian braiding statistics they can engender, could enable topological quantum computation [2,3]. In contrast, a series of no-go theorems rules out the existence of SPT phases and bosonic analogues of MFs -Majorana bosons (MBs) -in systems described by a gapped, stable quadratic Hamiltonian [4]. Non-trivial topology may emerge if one allows for gapless phases [5], strong interactions [6][7][8], or instabilities [9][10][11]. Does this then altogether rule out Majorana physics for non-interacting bosons and the possibility to build insight starting from a single-particle picture?

The search for topological order in non-interacting matter has extended into the realm of open quantum dynamics -both in a semiclassical limit described by non-Hermitian effective Hamiltonians [12][13][14][15][16], and a fully quantum setting described by a quadratic Markovian semigroup [17]. For fermionic matter, the use of engineered dissipation has been shown to provide a compelling paradigm for accessing topologically non-trivial steady states, and significant advances have been made towards uncovering a dissipative bulk-boundary correspondence which links bulk invariants to dissipative Majorana edge modes [18,19]. By furnishing almost conserved operators, these long-lived modes can dramatically alter the transient dynamics. They elicit correlation and relaxation times that may diverge with system size [20,21], reminiscent of metastability [22]. For bosons, the situation remains far more opaque. On the one hand, connections have been made between certain bulk topology and directional amplification of an input driving field [23][24][25][26]; on the other hand, there has yet to be any identification of the analogous MBs that are so fundamentally tied to SPT phases.

In this Letter we establish the existence of MBs for systems described by quadratic bosonic Lindbladians (QBLs). MBs are linked to non-trivial bulk topology and to a form of topo-logical metastability that we argue is unique to bosons. By leveraging tools from pseudospectral analysis [27], we show that, despite the system being dynamically stable for all finite size, unstable behavior may arise in the thermodynamic limit, enforcing anomalously long prethermalization and transient amplification. Metastability is necessary but not sufficient for MBs, however: topological metastability, featuring MBs, additionally requires that a winding number around zero be nonvanishing. Consistent with the fact that symmetries and conserved quantities are generally independent in open dynamics [28][29][30][31], we find that a MB pair generically comprises a conserved zero mode (ZM) localized on one edge, and a symmetry generator, localized on the other -together forming a split bosonic degrees of freedom. The symmetry is responsible for a continuum of quasi-steady states that can lead to unexpected features, like persistent non-Gaussianity. We identify simple models and show how MBs lead to distinctive signatures in experimentally accessible steady-state power spectra. Thanks to the robustness of pseudospectra against perturbations, our results carry through in the presence of weak disorder.

Quadratic bosonic Lindbladians.-Consider the Lindblad master equation ρ(t) = L(ρ(t)) for a density operator ρ(t) on an N -mode bosonic Fock space. The observables A carry the time dependence in the adjoint ( ⋆ ) Heisenberg picture, so that

where the Hamiltonian H = H † need not be bounded below and the dissipator D ⋆ is bilinear in the Lindblad operators

T be a Nambu array of bosonic annihilation, creation operators. Then, H ≡ 1 2 Φ † HΦ is a quadratic bosonic Hamiltonian associated to the matrix H = τ 1 H * τ 1 . Hereafter, τ j ≡ 1 N ⊗ σ j , in terms of the Pauli matrices. As the Lindblad operators

The equations of motion of the linear forms Φ(t) and the quadratic forms

with the dynamical matrix

. The rapidity spectrum, σ(-iG), is the set of eigenvalues of -iG. Since G = -τ 1 G * τ 1 , we have σ(-iG) = σ(-iG) * .

A quantum system is dynamically stable if the expectation value of an arbitrary observable in any state is bounded for all t ≥ t 0 . For a QBL, a sufficient condition for stability is that -iG is a Hurwitz matrix, that is, σ(-iG) is bound to the open left-half plane. In this case, there exists a unique, globally attractive, Gaussian steady state ρ ss , L(ρ ss ) = 0, which is completely determined by the expectation values ⟨Φ⟩ ss = 0 and ⟨Q⟩ ss = Q ss [32,33]. When, in addition, G is diagonalizable, L is as well and its spectrum follows from σ(-iG) [34][35][36]. The convergence to the steady state is exponential, with an asymptotic rate determined by the spectral gap,

The minimum time it takes for d max (t) to fall below a pre-determined accuracy δ > 0 is the mixing time t mix (δ) of the semigroup generated by L [17].

We focus on one-dimensional, bulk translation-invariant QBLs. There are three basic configurations: periodic (PBCs), open BCs (OBCs, two terminations), and semi-open BCs (infinite system with one termination). With hindsight, an infinite system without terminations is well described as the limiting case of PBCs. The dynamical matrices of these QBLs are known as block-Toeplitz or circulant matrices for OBCs and PBCs, respectively [37,38]. Translation invariance implies the rapidity spectrum for PBCs form closed curves in C as N → ∞. We will rely on two model QBLs for illustration. The Hamiltonian is the bosonic Kitaev chain (BKC) [10,11],

) where J ≥ ∆ ≥ 0 are nearest-neighbor hopping and nondegenerate parametric amplification amplitudes, and µ ∈ R is a uniform degenerate parametric amplification strength [39]. It is the dissipation mechanisms that separates the two models:

• Model 1: Uniform onsite dissipation, L j = √ 2κ a j for j = 1, . . . , N , with κ > 0 the uniform damping strength.

• Model 2: Add to the previous model next-nearest-neighbor damping of strength Γ > 0.

Metastable QBLs and ϵ-pseudospectrum.-Model 1 is dynamically unstable for PBCs whenever ∆ + |µ| ≥ κ, whereas it is stable for OBCs when κ -|µ| > 0 (see in Fig. 1 for µ = 0 and [40] for details). The dependence on BCs is dramatic for both models and cannot be weakened by increasing N . As we will show, any QBL that is (i) dynamically stable for OBCs, and (ii) unstable for PBCs is metastable, in the sense that the relaxation proceeds in two steps: a long-lived transient where the system behaves as if it were in an unstable phase, followed by stable asymptotic dynamics. To contrast, we shall call any system that is dynamically stable but not metastable unconditionally stable. This metastable behavior is not captured, in general, by the spectral properties of the QBL. The problem is that G is non-normal and so the spectral decomposition is fragile; one must drop it in favor of the more robust notion of the ϵ-pseudospectrum [27], as also recently used for non-Hermitian Hamiltonians [41,42].

For a matrix X and ϵ > 0, the ϵ-pseudospectrum of X is

The vectors ⃗ v are the ϵ-pseudoeigenvectors. When X = -iG is a Toeplitz matrix, the ϵ-pseudospectrum can be computed by analyzing the winding number of the bulk rapidity bands for PBCs [43]: if a rapidity band winds around λ ∈ C, then λ eventually joins the ϵ-pseudospectrum for any ϵ > 0 as the system size increases [40]. Thus, the semi-open limit of a metastable system is dynamically unstable. Unlike the spectrum, the ϵ-pseudospectrum is robust against perturbations, in that it scales linearly in the perturbation strength [27]. Pseudospectra control the transient behavior of a linear dynamical system. Notably, for all ϵ > 0, we have that sup t≥0 ∥e Xt ∥ 2 ≥ sup Re(σ ϵ (X))/ϵ [27].

The above result may be used to lower-bound the mixing time of an arbitrary QBL. A simple proxy for t mix (δ) is the linear mixing time t lin (δ), determined by the time it takes for the worst-case distance d lin (t) ≡ sup ⟨Φ⟩ 0 ∥⟨Φ⟩ t -⟨Φ⟩ ss ∥ ∥⟨Φ⟩ 0 -⟨Φ⟩ ss ∥ to drop below δ > 0. From Eq. ( 2), one finds d lin (t) = ∥e -iGt ∥ 2 . If the QBL is metastable, one can show [40] that for any fixed δ < sup t≥0 d lin (t) and r ≥ 0, there exists N max such that

where Ω > 0 is finite and determined by the dynamical matrix G B of the bi-infinite system. This is consistent with the observation that systems exhibiting the non-Hermitian skin effect can experience abnormally long relaxation without a closing of the spectral gap in the thermodynamic limit and suggests that pseudospectra may explain this generally [44][45][46]. Another manifestation of metastability is a time of transient amplification that increases with system size. In Fig. 2(a)-(b) we plot the expectation values of the boundary observable x N over uniformly sampled initial conditions, in the metastable and stable regimes for the two models. In the former regime, x N is exponentially amplified (while x 1 decays, not shown). In the unconditionally stable regime, exponential decay dominates, regardless of spatial location. We also plot in (d) the mean peak location for ⟨x N ⟩ t in the metastable regime against system size. This quantity, which lower-bounds a generic relaxation time-scale, diverges linearly with N .

Topological metastability and MBs.-The MFs of condensed matter physics are the Hermitian zero-energy boundary modes of topologically non-trivial superconductors: e.g., the fermionic Kitaev chain (FKC) [47] can display two MFs localized on opposite ends of the chain. These modes can be combined into a single (Dirac) fermionic ZM with support only on the boundary. Tight bosonic analogues of MFs do not exist for thermodynamically stable quadratic bosonic Hamiltonians [4]. Hints of them appear if stability is lifted: e.g., the BKC of Eq. ( 4) with J = ∆ hosts approximate ZMs [11],

which mimic closely the MFs of the FKC when |µ| < J: they are Hermitian, exponentially boundary-localized, are (almost) conserved, [H,

0 , so they can be normalized and combined to form a split bosonic mode. However, the system is either dynamically unstable, or on the cusp of instability.

One can circumvent the dynamical instability of the BKC by allowing for dissipation. For Model 1 (J = ∆), the modes

L ) provides a right-(left-)localized approximate ZM whenever |δ -| < 1 (|δ + | < 1), that is, precisely when the x(p)-quadrature rapidity band winds around the origin. These operators have the physical meaning of (approximately) conserved quantities, and coincide with Eq. ( 6) as κ → 0.

Noting that the system is stable whenever |µ| -κ < 0, it appears that we achieved our goal of finding MBs. However, there is a glaring discrepancy with the FKC: whenever κ ̸ = 0, it is not generally possible to combine γ c L and γ c R into a split bosonic degree of freedom. This follows because (i) they can exist independently of each other (see also Fig. 3); and (ii) they need not be canonically conjugate. For example, when µ = 0, [γ c L , γ c R ] = 0 for N even. Instead, the conjugate modes required to form the split degree of freedom are

which satisfy

We propose that the pairs (γ c L , γ s R ) and (γ s L , γ c R ) are the bosonic analogues of the MFs. There is a doubling of MB pairs with respect to the FKC. This is consistent with the breakdown of Noether's theorem for open systems [28,30], in the sense that conserved quantities (c) and symmetry generators (s) need no longer coincide (although a correspondence still exists for linear forms [40]). Explicitly, while in the Hamiltonian case the edge modes can each be thought of as supplying both a conserved quantity and a generator of a U(1) symmetry, they now generically split into pairs of distinct ZMs and generators. The latter, given in Eq. ( 7) obey

∀A, and similarly for γ s R , with

L,R , θ ∈ R, the dynamics are invariant under the unitary action (or "weak" symmetry [29,40]) L(U (θ)ρU (θ) † ) = U (θ)L(ρ)U (θ) † , up to exponentially small corrections in system size.

To see the connection between MBs, topology, and pseudospectra more generally, let us associate to each linear form

From the definition of ϵ-pseudospectra, one can verify that σ ϵ (G) = σ ϵ ( G). So, when the rapidity bands wind around λ = 0, it follows that 0 ∈ σ ϵ ( G), with ϵ exponentially small in N . The conserved MB is then associated with an ϵ-pseudoeigenvector of G; the symmetry generator MB is formed from the corresponding ϵ-pseudoeigenvector of G [40]. It may be possible for these pseudoeigenvectors to coincide, yielding 'non-split' MBs. Importantly, MBs survive in the presence of disorder, thanks to the robustness of pseudospectra [40]. We call this phenomenon topological metastability to contrast it with systems like Model 2, which can have vanishing winding number at zero (hence no MBs) when metastable. This topological metastability is forbidden in quadratic fermionic Lindbladians, due to the rapidity spectra being bounded to the left-half plane [35].

The approximate symmetries arising from topological metastability imply a degenerate quasi-steady state manifold. Let γ s be any linear approximate symmetry and consider the Weyl displacements ρ θ ≡ e iθγ s ρ ss e -iθγ s , with ρ ss the unique, Gaussian steady state under OBCs. In general, these quasisteady states are arbitrarily long-lived. To see this, note that ρ θ is also Gaussian, hence its lifetime is fully determined by those of the first and second moments, ⃗ m θ (t) ≡ tr[Φρ θ (t)] and Q θ (t) ≡ tr[Qρ θ (t)]. One can show that

with Ω the size-independent constant of Eq. ( 5) and a similar bound holding for Q θ (t). Thus, the expectation of an arbitrary linear form ⟨α⟩ t /[∥⃗ α∥∥ ⃗ m θ (0)∥] remains within an accuracy δ from ⟨α⟩ 0 /[∥⃗ α∥∥ ⃗ m θ (0)∥] for all times t with te Ωt < δ/ϵ [40]. Since, for these systems, ϵ decays exponentially with N , this condition can be met for arbitrarily long t (see also Fig. 2(c)-(d)). The quasi-steady states ρ θ can be used to construct long-lived initial states with properties uncharacteristic of stable QBLs: e.g., convex combinations of the ρ θ 's are generically non-Gaussian. Moreover, the existence of longlived states with non-zero first moments is unique to systems exhibiting topological metastability [40].

Towards experimental signatures of MBs.-Consider the two-time correlators of two linear forms, C α,β (t, τ ) ≡ ⟨α(t + τ )β † (t)⟩ = tr[α(τ )β † (0)ρ(t)], where the second equality relies on the quantum regression theorem [32]. Because of Eq. ( 2 the (one-sided) Fourier transform of this correlator as t → ∞ yields the steady-state power spectrum,

Within input-output theory, χ N (ω) is related to the susceptibility matrix that determines the transformation of an input field to an outgoing one, and a correspondence exists between topological properties of G and exponentially enhanced endto-end amplification in driven-dissipative photonic lattices [23][24][25][26]. Utilizing pseudospectra, a non-zero winding of the bulk rapidity spectra around iω mandates that ∥χ N (ω)∥ diverges with N , signaling amplification [48]. While this amplification by itself is a manifestation of metastability, a necessary condition for the existence of MBs is that ∥χ N (0)∥ = ∥G -1 ∥ diverges with N . To isolate the influence of MBs, we must account for possible divergent behavior of steady-state correlations, by computing the normalized power spectra [49] Sα,β (ω) ≡ S α,β (ω)/C α,β (∞, 0). Hence, we expect MBs to be diagnosed by a divergent zerofrequency peak of certain steady-state normalized power spectra. We assess this conjecture in Fig. 4: a clear zero-frequency peak is seen for the BKC, in contrast to the MB-free model. Moreover, this peak appears to diverge as a power law in N .

Outlook.-We have shown that Majorana bosons, despite being strictly forbidden for gapped free-boson Hamiltonians, can emerge in metastable QBLs. The appropriate phase diagram characterizes dynamical rather than thermodynamical stability. While there are natural ways to extend these ideas to higher dimension, a key question for characterizing SPT phases of free bosons is whether dynamical stability phase diagrams are dictated by some topological classification with an associated, dissipative bulk-boundary correspondence. The nature of our Majorana bosons suggests that the breakdown of Noether's theorem may play a role in answering this question. Experimental realizations of metastable QBLs offer another exciting venue for future research. Since the basic nonlinearities and dissipative couplings are readily available [50][51][52], superconducting arrays appear especially well-positioned to possibly uncover Majorana bosons and probe their physics.