Deconstructing effective non-Hermitian dynamics in quadratic bosonic Hamiltonians

Unlike their fermionic counterparts, the dynamics of Hermitian quadratic bosonic Hamiltonians are governed by a generally non-Hermitian Bogoliubov-de Gennes effective Hamiltonian. This underlying non-Hermiticity gives rise to a dynamically stable regime, whereby all observables undergo bounded evolution in time, and a dynamically unstable one, whereby evolution is unbounded for at least some observables. We show that stability-to-instability transitions may be classified in terms of a suitably generalized $\mathcal{P}\mathcal{T}$ symmetry, which can be broken when diagonalizability is lost at exceptional points in parameter space, but also when degenerate real eigenvalues split off the real axis while the system remains diagonalizable. By leveraging tools from Krein stability theory in indefinite inner-product spaces, we introduce an indicator of stability phase transitions, which naturally extends the notion of phase rigidity from non-Hermitian quantum mechanics to the bosonic setting. As a paradigmatic example, we fully characterize the stability phase diagram of a bosonic analogue to the Kitaev-Majorana chain under a wide class of boundary conditions. In particular, we establish a connection between phase-dependent transport properties and the onset of instability, and argue that stable regions in parameter space become of measure zero in the thermodynamic limit. Our analysis also reveals that boundary conditions that support Majorana zero modes in the fermionic Kitaev chain are precisely the same that support stability in the bosonic chain.


Introduction
Systems of non-interacting ("free") fermions or bosons have long provided a simple yet paradigmatic setting for investigating the differences that their statistical properties elicit in both equilibrium and non-equilibrium many-body physics [1]. The relevant model Hamiltonians, which are quadratic in the respective canonical fermionic or bosonic operators, either describe truly non-interacting degrees of freedom or, more commonly, they may arise from simplified treatments of interaction effects, for instance via random-phase or meanfield approximations [2]. At equilibrium, a striking manifestation of the different underlying statistics stems from the ability of bosons to exhibit Bose-Einstein condensation into a macroscopically occupied quantum state [3], as opposed to the formation of a sharp Fermi surface in systems non-interacting fermions at zero temperature. Away from equilibrium, the dynamical behavior of free fermions and bosons can also be dramatically different, however [1,4,5]. The Heisenberg equations of motion for the canonical creation and annihilation operators can be described in terms of an effective single-particle Hamiltonian (SPH), whose associated Schrödinger-like equation is precisely the Bogoliubov de-Gennes equation. Solving this equation yields the normal modes, which determine the elementary excitations of the model, as well as the dynamics of any physical observable of interest.
For both fermions and bosons, the effective SPH inherits a complex structure, the charge conjugation operation, arising from the fact that canonical creation and annihilation operators are mutually adjoint. Similarly, effective SPHs inherit a type of Hermiticity structure from underlying statistics. However, while fermionic effective SPHs are Hermitian with respect to the canonical inner product, bosonic effective SPHs are Hermitian with respect to an indefinite inner product -namely, they are "pseudo-Hermitian" [6,7] (see also [8] for an early physical application). As a consequence, fermionic effective SPHs are always diagonalizable by a unitary Bogoliubov transformation and possess a purely real spectrum, with the corresponding normal modes always exhibiting periodic (bounded) time evolution and obeying canonical anti-commutation relations. From a dynamical-system standpoint, quadratic fermionic Hamiltonians can only be dynamically stable. In contrast, a bosonic effective SPH can possess non-real eigenvalues and lose diagonalizability at exceptional points (EPs) in parameter space, where both eigenvalues and eigenvectors coalesce [9]; the corresponding normal modes can display both oscillatory and non-oscillatory dynamics, leading to unbounded growth (or decay) of physical observables in time, and can satisfy a wide range of algebraic relationships beyond the canonical commutation relations. In other words, both dynamically stable and dynamically unstable behavior is possible in general for a many-body system described by a quadratic bosonic Hamiltonian (QBH).
Physically, the emergence of effective non-Hermitian dynamics at the single-particle level is a direct manifestation of bosonic statistics in systems where particle-number conservation is broken at the manybody level. Such a scenario is thus typical for systems of massless bosons, conspicuous examples of which could include photons [10], phonons [11], and magnons [12]. Even for systems of massive bosons such as cold atoms, the mean-field quasiparticles and the small fluctuations of Bose-Einstein condensates are well described by non-particle-conserving QBHs [3]. As a result, the onset of dynamical instabilities has been linked to a wide range of phenomena, including certain Goldstone modes [1], parametric amplification by coherent driving in photonic systems [13], decay mechanisms for atoms in optical lattices [14], robustness properties of topological edge modes [15,16], and the formation of spin domains in spinor Bose-Einstein condensates [17].
Interest in the dynamical behavior of QBHs has heightened in recent years due to several reasons. On the one hand, topology plays an extraordinarily important role for fermions. Even in the absence of strong interactions, topological features of fermionic SPHs explain remarkable, robust physical effects like the integer quantum Hall that informs the current metrological standard for resistance. Along with the many successes of the topological classification of mean-field fermionic matter [18], this is prompting researchers to look for roles that topology may play for free bosons. Although no conclusive notion of a topological classification has been established as yet (see however [19,20] for up-to-date discussions on this fast-evolving subject), topologically non-trivial bands are by now well-documented for QBHs and a rigorous bulk-boundary correspondence has been identified, relating these bands to surface bands away from zero-frequency [16,20]. Likewise, topology has been conjectured to constrain dynamical rather than thermodynamical properties of free bosons [21], bringing to the fore the possibility that topological features might be most apparent in systems that need not be thermodynamically stable. On the other hand, the emergence of effective non-Hermitian dynamics makes QBHs natural candidates for the wealth of distinctive phenomena that non-Hermitian open dynamical systems are known to exhibit [22], including asymmetric mode switching, enhanced EP sensing, and non-Hermitian skin effect [9,23,24]. Meanwhile, alongside theoretical progress, a variety of photonic setups are becoming experimentally available as powerful platforms for simulating topological states of light and matter, and probing their statistical properties away from equilibrium [25].
Our goal in this paper is to obtain a deeper, unified understanding of the possible dynamical regimes that QBHs can support, and the ways in which transitions between different regimes may occur as parameters in the Hamiltonian and system size are changed. That is, for a given QBH, we aim to characterize its dynamical stability phase diagram, for any given system size. We tackle this problem by combining tools from non-Hermitian quantum mechanics [26,7,9] and fundamental results from linear algebra in indefinite inner-product spaces [27] with a particular type of stability theory of linear time-invariant dynamical systems, known as "Krein stability theory" [28]. In essence, our work leads to the following general conclusions: (i) All bosonic effective SPHs are PT -symmetric in a suitable sense; (ii) Stability-to-instability transitions are associated with breaking of this generalized PT symmetry; (iii) The transitions between dynamical phases can be detected by a new type of "phase rigidity" indicator that we call the Krein phase rigidity.
With reference to Fig. 1 for a pictorial summary of the different elements that enter our analysis and their interconnections, let us put our main results above in context.
(i) & (ii) All bosonic effective SPHs possess a generalized PT symmetry.-There are many investigations into the connections between pseudo-Hermiticity, PT symmetry, and more general antilinear symmetries in the mathematical physics literature. Motivated by the necessary conditions for a given diagonalizable non-Hermitian matrix to possess a real spectrum, it has been proposed that pseudo-Hermiticity should be regarded as the natural generalization of PT symmetry [6,7]. As we emphasize in our work, the condition of pseudo-Hermiticity is geometric in origin: it is equivalent to self-adjointness with respect to a possibly indefinite inner product. A very recent investigation shows that PT symmetry implies pseudo-Hermiticity, Figure 1. Pictorial summary of the main concepts relevant to our analysis, along with their logical interconnections. An arrow from box A to box B indicates that A leads naturally to B in some particular way. Black arrows (boxes) indicate connections (concepts) previously established in literature, while red arrows (boxes) indicate new connections (concepts) we established (introduced) in this paper. The red dashed arrows (and box) indicate connections and concepts we further elaborate on in a companion paper [29]. regardless of diagonalizability [30]. In Sec. 3.1, we complete this line of reasoning, by proving that the converse is also true if one allows for a suitably modified, yet mathematically and physically reasonable, notion of PT symmetry that we dub generalized PT (GPT ) symmetry. It follows that all bosonic effective SPHs, which are inherently pseudo-Hermitian as a reflection of bosonic canonical commutation relations, are GPT -symmetric. In particular, one can always classify bosonic dynamical phases according to whether this GPT symmetry is broken or not. Physically, as we also show, unbroken GPT symmetry is intimately connected to the fact that all normal modes may be chosen to obey canonical commutation rules and a normalizable Fock vacuum may be constructed, even for thermodynamically unstable systems. In fact, dynamical stability allows one to identify a dual, number-conserving Hamiltonian that is unitarily equivalent to the original, as we explicitly show in a separate study [29].
(iii) Krein Phase Rigidity.-One of the hallmarks of non-Hermitian quantum mechanics are EPs. The standard indicator of an EPs for complex-symmetric non-Hermitian Hamiltonians is a scalar quantity called the "phase rigidity" [31,32,33]. However, transitions between dynamical phases can occur without EPs. In Sec. 3.2, we introduce a KPR indicator which, in addition to detecting EPs, also leverages the indefinite inner-product geometry specific to bosons to detect transitions caused by the splitting of degenerate real eigenvalues into non-real ones, that entail no loss of diagonalizability -a so-called "Krein collision". In this sense, our KPR is a powerful extension of the phase rigidity to bosonic many-body systems, reducing to it whenever the relevant effective SPH happens to be complex-symmetric. Since there exist protocols for simulating non-Hermitian Hamiltonians of varied physical origin using photonic arrays [22], this establishes one of many possible routes for extending the applicability of our techniques beyond systems of free bosons.
To exemplify the application of our general theory of stability phase transitions for QBHs in a concrete setting, in the second half of the paper (Sec. 4) we present a detailed investigation of the dynamical stability phase diagram of the bosonic Kitaev-Majorana chain (BKC) introduced in [21]. Aiming to a fuller comparison between the fermionic and bosonic models and the interplay between bulk and boundary, we explore a broader family of BCs that interpolates smoothly between open and periodic and allows for an arbitrary twisting angle. Most of our work on the BKC is analytical or supported by analytical results, which is important in several respects. For example, it can be difficult to distinguish numerically whether a stability transition point is associated to loss of diagonalizability of the effective SPH (an EP), or a Krein collision, or both. Another difficult but crucial problem is that of calculating the dynamical stability phase diagram in the limit of infinite system size. While dynamical stability phase diagrams and thermodynamic phase diagrams are a priori very different objects, it is sensible to ask whether the limit of infinite system size can highlight gross regularities within the class of all possible dynamical phase diagrams. Our analysis for the BKC suggests that the answer is likely in the affirmative. For both finite and infinite size, the key to accessing spectral properties analytically is the exact diagonalization algorithm for corner-modified block-banded Toeplitz matrices of [34]. In the context of fermions, this algorithm yields a generalized Bloch's theorem for clean systems under arbitrary boundary conditions (BCs) [35]. Here, we explicitly apply this algorithm to bosonic SPHs for the first time, which in itself constitutes a new result of independent interest.
Finally, as we remarked above, in [21] the authors conjecture that certain topological properties of the BKC may be responsible for its peculiar dynamical features -including the strong sensitivity to BCs and the ability to propagate excitations in a "chiral" fashion, depending on the phase of the excitation. The BKC model itself is constructed by applying a particular mapping from free fermions to free bosons to the fermionic Kitaev chain at zero chemical potential. In Sec. 5, we characterize this mapping in full generality and show that can it only preserve certain topological invariants at the cost of mapping fermions to either dynamically unstable bosons or bosons that are at the cusp of instability. For the BKC, this is precisely the difference between closed (unstable) and open (at the cusp of instability) BCs. Hence, Krein stability theory explains the fragility of the dynamically stable phase of the open BKC against bulk disorder observed in [21]. Altogether, our analytical solutions point to a remarkable correspondence between (open and π/2-twisted) BCs and parameter regimes that support Majorana zero modes in the fermionic Kitaev-Majorana chain versus, respectively, dynamical stability in its bosonic counterpart. Further to that, as we outline in Sec. 5, bosonic models may also be constructed, which host boundary-localized, Hermitian zero-frequency analogs to Majorana modes -although, ascertaining a clear connection with topology remains a fascinating topic for further research.

Background
QBHs have been extensively considered in the literature as integrable models for both thermodynamically stable systems [4,1,2] and beyond [5,36,37,38,39]. In this section, after introducing the basic concepts, we provide a self-contained review of the mathematical techniques needed for casting these Hamiltonians in normal form, to a level of detail appropriate for subsequent analysis.

From quadratic bosonic Hamiltonians to effective non-Hermitian dynamics
Let a † i and a i denote canonical bosonic creation and annihilation operators for a mode i, satisfying The requirement that H be Hermitian imposes the condition K ij = K * ji . Furthermore, the bosonic commutation relations allow us to take ∆ ij = ∆ ji . In terms of the Nambu (column) arrayΦ ≡ [a 1 , a † 1 . . . , a N , a † N ] T , the above QBH can be re-written as H = 1 Note that our definition of Nambu arrays differs from the standard orderingΦ ≡ [a 1 , . . . , a N , a † 1 , . . . , a † N ] T . While this rearrangement clearly does not affect the physics, it will be instrumental to bring the role of translation invariance to the fore and enable a more straightforward application of the diagonalization techniques to be employed in Sec. 2.3. Beside being Hermitian, the matrix H in Eq. (2) obeys the condition We will also denote τ 2 ≡ 1 N ⊗ σ 2 and τ 3 ≡ 1 N ⊗ σ 3 , with σ j , j = 1, 2, 3, being the usual Pauli matrices. While formally H plays the role of a SPH, a key difference between bosonic and fermionic quadratic forms arises from the fact that diagonalizing H does not, in general, imply diagonalization of H [1], nor does it characterize the dynamics that H generates. The latter is determined by the solution of the Heisenberg equations which, in units where = 1, may be compactly written as A more transparent way to interpret this equation may be obtained by considering an arbitrary (column) vector |α ∈ C 2N and define α ≡ α| τ 3Φ (e.g., for a single mode, α = α * 1 a 1 − α * 2 a † 1 ). It then follows that (i) many-body commutators induce a geometric structure on single-particle space, in the sense that [ α, β † ] = α|τ 3 |β 1 F ; and (ii) the many-body adjoint operation similarly induces a charge-conjugation operation, α † = − Cα, with C ≡ τ 1 K = C −1 and K being complex conjugation. In order to exemplify the single-particle structure of the dynamics in Eq. (4), we further define the convention α(t) ≡ α(t)| τ 3Φ (0), that is, we absorb the time-dependence of the bosonic operators a j , a † j into the vector of coefficients |α . Then, by taking advantage of the identity −[ H, α(t)] = Gα, one finds that α(t) satisfies the Heisenberg equation of motion if and only if d dt Some symmetries of this single-particle equation are inherited from, and/or can be can be lifted back, to the original many-body problem. For example, CG = −GC in terms of the charge conjugation operation identified above. As for ordinary, commuting symmetries, suppose U is a "τ 3 -unitary" (or para-unitary) matrix, that is,  (5), the dynamics of a free bosonic system are governed by the effective SPH G which need not be Hermitian or even normal (in fact, G is Hermitian if and only if the "pairing" contribution ∆ ij vanishes). As a consequence, G can have non-real eigenfrequencies as well as non-trivial Jordan chains. Since, from a dynamical-system standpoint, Eq. (5) defines a linear time-invariant system with state matrix −iG, the normal modes of the dynamics generated by H are built up from the (generalized) eigenvectors of G. Two distinct notions of stability are then relevant for the QBHs in question: (i) The system H is dynamically stable if G is diagonalizable and all of its eigenvalues are real.
(ii) The system H is thermodynamically stable if there exists a finite lower bound on the expectations of H.
On the one hand, the normal modes of a dynamically stable QBH exhibit strictly bounded motion [16]: the onset of dynamical instability is signaled by the appearance of a non-trivial Jordan chain or a complex eigenvalue in the eigensystem of G ‡. On the other hand, the notion of thermodynamic stability matters because a mean-field ground state need not exist in the bosonic Fock space: for example, if there is a bosonic excitation of negative energy, in the thermodynamic limit one may occupy it with an arbitrary number of bosons and lower the energy of the system without ever hitting a ground state. Such an instability is sometimes called a Landau instability [17] and we will see an example in Sec. 4. Thermodynamic stability can be diagnosed in terms of the the Hermitian matrix H = τ 3 G: If H is positive semi-definite, then H is thermodynamically stable [1,4]. The converse implication, while likely to hold, is more subtle [38]. In particular, if H ≥ 0 and G is diagonalizable, then the associated QBH is both thermodynamically and dynamically stable. Nonetheless, the two notions of stability are independent, as one can see from the following examples: (1) A system of N decoupled harmonic oscillators is thermodynamically and dynamically stable. (2) A quantum harmonic chain with nearest-neighbor (NN) couplings (one-dimensional acoustic phonons) is thermodynamically stable and displays a non-trivial Jordan chain associated to the conserved total momentum and center of mass operators. Hence, it is dynamically unstable. (3) As we will see in Sec. 4, there are parameters regimes of the bosonic analogue of the Kitaev chain introduced in [21], in which the model is dynamically yet not thermodynamically stable.

Normal form of a quadratic bosonic Hamiltonian
In order to elucidate how the eigensystem of G determines the normal modes of Eq. (4) and a normal form of the corresponding QBH, we start from identifying some intrinsic symmetry properties that the spectrum of G enjoys. Firstly, Hermiticity of H implies that G † = τ 3 Gτ 3 . As a consequence, if ω n is an eigenvalue of G, then so is ω * n . Secondly, the charge conjugation symmetry CG = −GC is equivalent to G * = −τ 1 Gτ 1 . Hence, if ω n is an eigenvalue of G, then so is −ω * n . Accordingly, the eigenvalues of G come in quartets {ω n , ω * n , −ω n , −ω * n }. These two properties also have consequences on the eigenvectors of G and G † . Specifically, let |ψ n be an eigenvector of G associated to the eigenvalue ω n . Then, (i) τ 3 |ψ n is an eigenvector of G † with eigenvalue ω * n , and (ii) C |ψ n is an eigenvector of G with eigenvalue −ω * n .

Diagonalizable case.
Suppose that G is diagonalizable, in which case the eigenvectors of G form a complete basis of C 2N . Then, aside from minor notational differences, we may repeat the analysis in [5]. Let |ψ n * denote the eigenvector of G with eigenvalue ω * n if Im(ω n ) = 0, and |ψ n * = |ψ n otherwise. There are 2N eigenvectors |ψ n of G, which correspond to eigenvalues ω n , n = 1, . . . , 2N , and satisfy ψ n * |τ 3 |ψ m = δ nm , that is, there exists a τ 3 -orthonormal basis. In terms of this basis, This spectral decomposition of G = τ 3 H leads to the desired normal form of the QBH, namely: n=1 ω n ψ † n ψ n * − 1 2 tr K, (6a) ‡ Note that the condition of all eigenvalues of G being real is stronger than the standard (Hurwitz) stability condition for the linear system in Eq. (5), which would only require every eigenvalue to have a strictly positive imaginary part. This stems from the symmetry structure that the bosonic nature of the many-body problem imposes on the eigenvalue spectrum, see Sec. 2.2.
where | ψ n ≡ C |ψ n , [ ψ n * , ψ † m ] = δ nm 1 F , and [ ψ n * , ψ m * ] = 0 = [ ψ n , ψ m ]. At this point the analysis splits depending on whether ω n is or is not real. If ω n ∈ R, the pair ( ψ n , ψ † n * ) is called a canonical pseudo-bosonic normal mode [5] and there is no room for further simplification. By contrast, if ω n ∈ R, then one can choose the eigenvectors so that |ψ n * ∝ |ψ n and | ψ n is an eigenvector with eigenvalue −ω n . Furthermore, ψ n |τ 3 |ψ n = − ψ n |τ 3 | ψ n = 0 §. Hence, one can renormalize these states so that ψ n |τ 3 |ψ n = 1. This procedure yields a term of the form ω n ψ † n ψ n in the normal form of H. The pair ( ψ n , ψ † n ) is a bosonic normal mode, that is, it satisfies the canonical commutation relation [ ψ n , ψ † n ] = 1 F . The transformationΦ →Ψ ≡ [ ψ 1 * , ψ † 1 , . . . , ψ N * , ψ † N ] T from the physical bosonic modesΦ to the normal modes of the system is furnished by the modal matrix M , whose columns are the eigenvectors of G, that is, Ψ = M −1Φ . Inverting this transformation allows us to find the desired time dependence ofΦ. Explicitly, ψ n * (t) = e −iω * n t ψ n * (0), ψ n (t) = e −iωnt ψ n (0). Hence, the system is dynamically stable only if all the eigenvalues are real (we have excluded non-trivial Jordan chains by assumption). Otherwise, some normal modes are amplified exponentially with time, which can happen for example in parametrically driven optical systems. Each amplified mode is paired with an exponentially decaying (de-amplified) mode. Such a system is thermodynamically stable if it is dynamically stable and all the creation operators are associated to positive eigenfrequencies.

General case.
If G is not diagonalizable, at least one eigenvalue of G, say, ω 0 , must be associated to a non-trivial Jordan chain. Allowing for degeneracy, let |χ j1 for j = 1, . . . , N 0 denote a complete set of independent eigenvectors corresponding to ω 0 . Moreover, for each j, let r j > 1 be the length of the corresponding Jordan chain, satisfying We can then see that the normal modes defined by χ jk ≡ χ jk |τ 3Φ manifest as Jordan chains in the adjoint action of H, namely, and exhibit exponentially modulated polynomial time evolution, Hence, dynamical stability is excluded by non-trivial Jordan chains but, as it turns out, thermodynamic stability is not. Let us call the modes χ jk , k > 1, generalized normal modes of rank k. As we saw, one can enforce bosonic or pseudo-bosonic commutation relations for the diagonalizable sectors [37] (i.e., r j = 1). By contrast, the commutation relations of the generalized normal modes are more difficult to pin down. To gain physical insight into this problem, let us first focus on thermodynamically stable systems. When H ≥ 0, it is known that there can only be Jordan chains at ω 0 = 0 and their lengths are at most two (see Theorem 5.7.2 in [27]). The eigenvector and generalized eigenvector of a given Jordan chain can then be used to create pair of Hermitian normal modes satisfying Heisenberg-Weyl commutations relations. Physically, these Jordan chains often represent Goldstone modes [1]. Mathematically, the situation is neatly exemplified by a single mode with Hamiltonian H = p 2 /2m: then, p is the zero mode and any combination c 1 x + c 2 p, with c 1 = 0, can officiate as the generalized zero mode. Naturally, with hindsight, one chooses c 1 = 1, c 2 = 0, but there is nothing in the problem that makes this choice canonical. At best, one can require invariance under the charge conjugation operation and restrict c 1 , c 2 ∈ R. More generally (regardless of thermodynamic stability), since a generalized eigenvector can always be shifted by a constant multiple of an eigenvector, the commutation relations of the corresponding normal modes may obey a wide range of commutation relations. Notably, this ambiguity is employed in [5] to construct bosonic generalized normal modes at a non-zero frequency. § More precisely, if ωn = 0, we have ψn|τ 3 |ψn = ω −1 n ψn|τ 3 G|ψn = ω −1 n ψn|H|ψn = 0. The case ωn = 0 requires separate consideration; see [4,20] for a self-contained study of the zero subspace.
While the commutation relations of the normal modes are malleable, it is important to note that the evolution of observables is independent of the particular generalized eigenbasis used to determine their time dependence by way of Eq. (9). For concreteness, one can fix the commutation relations of a given Jordan chain by utilizing a particular Jordan normal form available for matrices satisfying G † = τ 3 Gτ 3 . By Theorem 5.1.1 in [27], we then know there exists an invertible matrix T , such that J = T GT −1 is a Jordan normal form for G and τ 3 = T † P T , with P a block-diagonal matrix whose j-th block is an r j × r j matrix, and either all 1's or −1's on the anti-diagonal and 0's elsewhere. By applying T −1 to the canonical basis of C 2N , we can construct Jordan chains at eigenfrequency ω j satisfying the orthonormality condition χ jk |τ 3 |χ p = ε j δ j * δ k,rj +1−p , with ε j ∈ {−1, 1} and * the index labeling a Jordan chain at ω * . The commutation relations of the corresponding normal modes then satisfy [ χ jk , χ † p ] = ε j δ j * δ k,rj +1−p 1 F . The commutators [ χ jk , χ p ] = − χ jk |τ 2 K|χ p 1 F , can then be determined on a case-by-case basis.

Diagonalization of corner-modified, banded block-Toeplitz matrices
As mentioned in the Introduction, our main goal is to develop a general theory of the possible dynamical behaviors of QBHs, by shedding light, in particular, on the extreme sensitivity of their dynamical response to BCs. A key tool for our analysis is an exact diagonalization procedure for G (loosely speaking, since G need not be diagonalizable) which, while originally developed with fermions in mind [34,35,40], works more generally for clean (disorder-free) quadratic Hamiltonians subject to arbitrary BCs. The key property is that G belongs to a class of structured matrices known as corner-modified, banded block-Toeplitz (BBT). In what follows, we summarize the essential steps of this diagonalization procedure both to explicitly show-case its first application to bosons and to make the presentation as self-contained as possible. For clarity, technical details are deferred to Appendix A.
Corner-modified BBT matrices arise naturally from QBHs whose couplings have finite range and possess translation invariance "up to a boundary". While the framework developed in [34,35,40] is more general (e.g., it allows for multiple internal modes), we specialize it here to a one-dimensional lattice of N sites with bosonic modes (a j , a † j ) attached to each site, and defineφ j = [a j , a † j ] T , 1 ≤ j ≤ N . The Nambu array of Sec. 2 then readsΦ = [φ 1 , . . . ,φ N ] T , and the relevant class of QBHs has the form where b, b ∈ {1, . . . , R, N − R + 1, . . . , N }, and h r , W bb are 2 × 2 matrices that couple the bosonic modes at different sites in the "bulk" and "boundary" respectively, with range R < N/2. In practice, one often has R N (e.g., R = 1 for NN couplings). As one can see from the structure of H in Eq. (10), translation invariance is maintained up to a boundary slab of thickness R. Hence, we speak of mildly broken translation invariance in these systems. The effective SPH associated to Eq. (10) can be split as and where we have . Mathematically, G O is an example of a BBT matrix, that is, a block matrix whose entries are constant along diagonals. It acts naturally on the tensor-product space C N ⊗ C 2 ≡ H L ⊗ H I , where the first (second) factor carries the lattice (internal) degrees of freedom. The bandwidth of G O , that is, the number of non-zero diagonals, is 2R + 1. Physically, one recognizes G O as the effective SPH of the system subject to open BCs.
The role of (mildly broken) translation invariance becomes more apparent if we define the left-shift operator, T ≡ N −1 j=1 |j j + 1|, which acts only on H L and is the finite-lattice truncation of the translation operator T ≡ j∈Z |j j + 1|. Importantly, T is taken to act on the infinite lattice space span {|j } j∈Z , without the corresponding 2 -inner product; thus, while remaining invertible, T is no longer unitary and may possess generalized eigenvectors in general. Using shift operators, the above G O rewrites as The effective SPH G of the translation-invariant system is recovered by replacing T with T . The matrix V is, in turn, a corner modification that encodes BCs other than open. Periodic BCs, for example, can be recovered by a suitable choice of V [35]. In the following, we will also need the matrix polynomial , which is the analytic continuation of the effective singleparticle Bloch Hamiltonian off the Brillouin zone.
Next, let us define the bulk and boundary projectors as The goal is to solve the eigenvalue equation G O |ψ = ω |ψ . Since P B +P ∂ = 1 2N and P B V = 0, the eigenproblem is equivalent to the following "bulk-boundary system of equations" [40]: The diagonalization proceeds by first solving the bulk equation, Eq. (13a), parametrically in ω, and then employing the resulting solutions as an Ansatz for the boundary equation, Eq. (13b). One can show that, generically, such a strategy yields all of the eigenvectors of G O + V and can also be applied for computing generalized eigenvectors [34]. For fixed ω ∈ C, the complete set of solutions to the bulk equation (13a) breaks up into three different types of solutions (for a derivation, see Appendix A). Solutions of the first type are obtained by restricting to the finite-lattice solutions of the translation-invariant equation (G − ω) n Ψ = 0, for some suitable n, and thus arise from eigenvectors and generalized eigenvectors of G. Specifically, these solutions take the form where the z are the roots of the equation det(G(z) − ω) = 0 with algebraic multiplicity s , and the vectors |z , ν are as follows: for ν = 1, |z , 1 = N j=1 z j |j represents a generalized Bloch wave, with possibly complex momentum k = −i log(z ); for ν > 1 the |z , ν are proportional to ∂ ν−1 z |z , ν , and hence contain amplitudes with a power-law pre-factor to the exponential weight z j . The other two types of solutions that can arise are localized on the boundary of the system and are no longer controlled by G and the corresponding (non-unitary) translation symmetry. Rather, they emerge entirely due to the truncation from the bi-infinite lattice to a finite one. We will denote these left (−) and right (+) localized emergent solutions by |ψ ± , with = 1, . . . , s 0 ≡ 2R − 1 2 n =1 s . Here, s 0 is the multiplicity of z = 0 as a root of det(G(z) − ω) = 0. Finally, we remark that there may exist exceptional, isolated values of ω, which physically correspond to dispersion-less "flat bands" and whose associated eigenvectors are not included among the previous three types of solutions. While we refer to [34] for more discussion, flat bands will not be encountered in the models under consideration in this paper.
The complete set of solutions to the bulk equation may thus be parameterized as follows: . Using |ω, α as an Ansatz for the boundary equation, Eq. (13b), leads to the identity Here, the boundary matrix B(ω) has elements B bs (ω) that, in our case, consist of 2×1 blocks and are given by us that if B(ω)α = 0, then |ω, α solves both the bulk and boundary equations and hence is an eigenvector of G O with eigenvalue ω, as desired. For diagonalizable matrices, the above procedure yields a Bloch-like diagonal basis. However, G may fail to be diagonalizable, in which case the generalized eigenvectors of G are needed in addition to the eigenvectors to complete a basis. One can calculate some generalized eigenvectors in Bloch-like form by repeating the above procedure (see also Appendix A) to determine ker (G − ω1 2N ) p for various powers p and each eigenvalue ω. However, there is a constraint p < (N − 1)/R ≡ p max on how high p can be because, for p ≥ p max , (G − ω1 2N ) p need not be a corner-modified BBT matrix. If there are any, generalized eigenvectors of rank greater than p max − 1 may have to be determined by means other than the above bulk-boundary separation. As it turns out, the models we consider offer examples of this peculiar phenomenon.

Manifestations of effective non-Hermitian dynamics
As we have just seen, the effective single-particle description of a QBH is based on mapping the Heisenberg dynamics of Fock space operators to the the dynamics of finite-dimensional vectors by way of Eq. (4). The resulting dynamical system, Eq. (5), features a structured matrix G that need not be Hermitian. The only two constraints on G are G † = τ 3 Gτ 3 , and G * = −τ 1 Gτ 1 . These constraints are the direct mathematical manifestation of the quantum statistics of bosons at the effective single-particle level. In this section, we develop a general framework for analyzing the dynamical system Eq. (5) subject to the two above constraints, with special emphasis on the consequences at the many-body level.

Pseudo-Hermiticity and generalized PT symmetry
The condition G † = τ 3 Gτ 3 is an instance of a more general mathematical concept known as pseudo-Hermiticity. A matrix M is pseudo-Hermitian if there exists a Hermitian, invertible matrix η such that M † = ηM η −1 [8,6,7]. A pseudo-Hermitian matrix need not be diagonalizable unless the spectrum of η is positive-or negative-definite, and that is precisely not the case for bosons since η = τ 3 . As we are now going to show, the conditions of pseudo-Hermiticity can nonetheless be recast as a particular form of anti-linear symmetry. Specifically, such an anti-linear symmetry turns out to be closely related to the PT symmetry of some non-Hermitian quantum systems [26,41]. In the current usage of the term, an n × n complex matrix M is said to be PT -symmetric if it commutes with an anti-linear operator of the form PT , with P linear and involutory, that is, P 2 = 1 n , and T the operation of complex conjugation with respect to the canonical basis. The physical origin or meaning of the PT symmetry is not important from a mathematical perspective. There are a series of results in the literature establishing a web of relationships between pseudo-Hermiticity and PT symmetry in the strict sense just described. In particular: • Every PT -symmetric matrix is pseudo-Hermitian [30].
• A diagonalizable matrix is pseudo-Hermitian if and only if it commutes with an anti-linear invertible mapping [42]. If, in addition, the spectrum the matrix is real and discrete, then there is a basis in which its anti-linear symmetry can be decomposed as a product of an involutory linear operator and complex conjugation with respect to this basis [7].
• Weak pseudo-Hermiticity (a condition that is more general that pseudo-Hermiticity but coincides with it if the spectrum is discrete) is equivalent to the existence of an involutory anti-linear symmetry [43,44].
In the following proposition, we show that a pseudo-Hermitian matrix is always PT -symmetric provided that one relaxes the notion of PT symmetry slightly: A PT -symmetric matrix in the usual sense is automatically GPT -symmetric with respect the canonical basis. We then have: Proof: If M commutes with an invertible anti-linear transformation, like a GPT symmetry, then it is necessarily pseudo-Hermitian by Theorem 3 in [44]. To establish the opposite implication, suppose M is pseudo-Hermitian. Let λ 1 , . . . , λ α denote the real eigenvalues of M (if any) and µ 1 , . . . , µ β denote the non-real eigenvalues of M in the upper half-plane (i.e., Im(µ j ) > 0) for j = 1, . . . , β). Furthermore, let r j and p j denote the lengths of the Jordan chains corresponding to λ j and µ j , respectively, in the Jordan normal form of M . Pseudo-Hermiticity implies that for each Jordan chain of length p j corresponding to eigenvalue µ * j , there is a Jordan chain of length p j corresponding to the eigenvalue µ * j (see e.g. Proposition 4.2.3 in [27]). Hence, we can construct a basis of H consisting of generalized eigenvectors of M , say, B ≡ {|v jk , |w jk , |w jk }, where |v jk , |w jk , and |w jk denote rank-k generalized eigenvectors of M at eigenvalues λ j , µ j , and µ * j , respectively, with Next, define an involutory linear transformation P on B by P |v jk ≡ |v jk , P |w jk ≡ |w jk , and P |w jk ≡ |w jk . Furthermore, let T denote complex conjugation with respect to the basis B and let Θ ≡ PT . It follows immediately that [M, Θ] |v jk = 0. Furthermore, we also have There is a shorter but non-constructive proof of necessity. Theorem 3 in [44] asserts that every pseudo-Hermitian operator possesses an anti-linear involutory symmetry. Furthermore, every anti-linear involutory operator coincides with complex conjugation in some basis [45]. Since conjugation is of the type PT , every pseudo-Hermitian operator possesses a GPT symmetry.
The following theorem is simply the instantiation of Proposition 3.1 to effective SPHs G of QBHs. Recall that these matrices are pseudo-Hermitian with respect to τ 3 , G † = τ 3 Gτ 3 .
The theorem is also true for fermions. The SPHs of fermionic systems are Hermitian, and hence pseudo-Hermitian with respect to η = 1 2N . In this case, the GPT symmetry of Proposition 3.1 is simply complex conjugation with respect to the eigenbasis of the SPH. The range of possibilities is richer for bosons because non-trivial Jordan chains are possible even for thermodynamically stable systems and so, even in the most elementary systems, this symmetry can be broken -a situation with no counterpart in fermions.
We say that the GPT symmetry Θ of a bosonic matrix G is unbroken if there exists a basis of simultaneous eigenvectors of G and Θ. In the GPT -unbroken phase, G is diagonalizable (by definition) and one can check that the spectrum of G is necessarily real. Hence, G is dynamically stable and, by combining these spectral features with results from Sec. 2.2, we see that the normal modes of the system may be chosen to satisfy canonical commutation relations. In the GPT -broken phase, in contrast, either G fails to be diagonalizable or the spectrum of G includes complex eigenvalues, and so G is dynamically unstable and at least some of the normal modes need not satisfy canonical commutation relations. Therefore, for bosons, the emergence of dynamical instabilities and non-canonical normal modes is rooted in symmetry breaking. With this mechanism firmly established, we call a transition between a dynamically stable and an unstable regime a dynamical phase transition (not to be confused with the notion of a phase transition occurring in the thermodynamic limit away from equilibrium).
Physically, a many-body feature distinctive of the GPT -unbroken phase is the existence of a convergent quasi-particle vacuum. Since all the normal modes of the system can be chosen to satisfy bosonic commutation relations, it follows that G can be diagonalized by a τ 3 -unitary matrix L satisfying By contruction, this transformation maps the physical bosonic modesΦ to the bosonic normal modesΨ = The quasi-particle vacuum |0 , that is, the state such that ψ j |0 = 0 for all j, is formally given by [1,37,2] where a j |0 = 0 for all j. One needs to check that |0 is normalizable just as |0 is. The normalizability of |0 is equivalent to the condition that all the singular values of X −1 Y be less than one. It is convenient to recast this condition as the requirement that all the eigenvalues of Z † Z, with Z ≡ X −1 Y , be strictly less than one. Since L is a canonical transformation, we know that XX † = 1 N + Y Y † . It follows that | det Y | 2 < | det X| 2 and this inequality in turn implies that det Z † Z < 1 (note that this also implies X is always invertible). Hence, the quasi-particle vacuum |0 is always normalizable, as claimed.
Knowing that |0 is normalizable, one may compute a basis of eigenvectors of H as where the eigenfrequencies ω m are precisely the eigenvalues of G. This familiar diagonalization procedure for H is only possible in the GPT -unbroken phase. Note that the quasi-particle vacuum |0 need not be the ground state unless the system also happens to be thermodynamically stable.

Characterizing GPT -symmetry-breaking phase transitions
The GPT symmetry of QBHs breaks if the effective SPH matrix G approaches a point in parameter space where it loses diagonalizablility -that is, an EP -or, when the spectrum splits off the real axis into the complex plane, regardless of diagonalizability. For the subclass of non-Hermitian matrices that are symmetric, so-called phase rigidity (PR) has been extensively used to detect EPs [31,32,33]. Let |ψ be an eigenvector of M = M T with eigenvalue λ. Then, |ψ * ≡ (|ψ ) * is an eigenvector of M † with eigenvalue λ * . If one imposes on |ψ the bi-orthonormalization condition ψ * |ψ = 1, then one can show that the PR, defined by vanishes smoothly as M approaches an EP. While useful, such an indicator is not sufficient to our purpose of detecting dynamical stability-to-instability transitions in bosons, because bosonic matrices G need not be symmetric to begin with and, in addition, such dynamical transitions can occur without loss of diagonalizability. In what follows, we introduce a new PR indicator and argue that it successfully characterizes transitions between GPT -broken and unbroken phases of QBHs. As it turns out, to do so it is necessary to leverage mathematical results from stability theory of linear dynamical systems governed by pseudo-Hermitian matrices, also known as Krein stability theory [28], as we recall next.
3.2.1. Tools from Krein stability theory. Let η denote an invertible, Hermitian linear transformation of C 2N featuring both positive and negative eigenvalues. The space C 2N endowed with the indefinite inner product (·|·) = ·|η|· is called a Krein space. An η-Hermitian linear transformation satisfies (M ψ|φ) = (ψ|M φ), for all ψ, φ. This condition can be seen to be equivalent to M being pseudo-Hermitian in the sense of Sec. 3.1, namely, M † = ηM η −1 . A vector |α ∈ C 2N is η-positive if α|η|α > 0 and η-negative if α|η|α < 0. If α|η|α = 0, then |α is η-null. If |α is either η-positive or η-negative, the sign of α|η|α is the Krein signature of |α . Given an η-Hermitian matrix M , the eigenspace E λ for the eigenvalue λ can be classified as follows: If all the eigenvectors of M with eigenvalue λ are η-positive (negative), we say that E λ , or λ itself, is η-positive (η-negative) definite. If we do not wish to specify the sign, we just call E λ , or λ, η-definite. If E λ is not definite, then we call it, or λ, η-indefinite. The following definition and results (adapted from [28], Chapter III) are then relevant in our present context:  This lemma has direct relevance to QBHs, as noticed in [16,46] (see also [47] for related, albeit less general, results). When the bosonic matrix G is diagonalizable and all of its eigenvalues are definite, then small perturbations of the system cannot drive it away from dynamic stability. By contrast, if some eigenvalue of G hosts a Krein collision, then it can be de-stabilized by an arbitrarily small perturbation. These results have important consequences for the stability of topological edge modes of QBHs [16]. Another interesting consequence of Lemma 3.1 is that the zero-frequency bosonic normal modes of dynamically stable QBHs are fragile. If G is dynamically stable, then these zero modes are constructed out of pairs of kernel vectors of G of opposite Krein signature. Put differently, QBHs that host zero modes are either dynamically unstable or at the cusp of dynamical instability by Lemma 3.1(ii). We investigate in detail bosonic zero-modes of thermodynamically stable systems in [20].
Thus, the KPR in Eq. (17) may be calculated by a formula seemingly identical to the PR of Eq. (16), once the input states are suitably normalized. In fact, the KPR coincides with the PR when G happens to be symmetric. To see why this is the case, consider for simplicity a non-degenerate, real eigenvalue ω of G = G T , with |ψ , |ψ the corresponding eigenvectors satisfying ψ|τ 3 |ψ = κ ∈ {1, −1} and ψ * |ψ = 1, respectively. Since ω is not degenerate, |ψ = µ |ψ , for some µ ∈ C. Comparing Eqs. (16) and (17), we see that KPR and PR coincide if we can set |µ| = 1. Since G is symmetric, both τ 3 |ψ and |ψ are eigenvectors of G † and pseudo-Hermiticity further implies that ω is also non-degenerate as an eigenvalue of G † . It follows that τ 3 |ψ = ν |ψ * for some ν. Equating norms and noting that τ 3 is unitary, we concludes that |µ| = |ν|. Then κ = ψ|τ 3 |ψ = µν * and taking the modulus of this equation one finally obtains that |µ| = 1, as desired. While the corresponding arguments for a non-degenerate, complex eigenvalue are more involved, the equality still holds.
Despite the above formal similarities, and the fact that the KPR consistently reduces to the PR for matrices that are pseudo-Hermitian and symmetric, the KPR is not the same quantity as the PR. Crucially, for the PR, the mapping from the eigenvector ψ to its bi-orthogonal partner |ψ * in Eq. (17) is antilinear (complex conjugation), whereas for the KPR of bosonic matrices, the corresponding mapping is linear and given by τ 3 . Thus, the bi-orthonormalization condition associated to the PR translates into a τ 3normalization condition for the KPR. This difference is ultimately responsible for the KPR to serve as a useful indicator of GPT symmetry breaking for arbitrary pseudo-Hermitian matrices -in particular, bosonic ones, which are our main focus in this paper. Specifically, our claim is as follows: Claim 3.1 If a bosonic matrix G undergoes a dynamical phase transition as a function of some parameter or parameters, then there exists at least one eigenvector of G such that its KPR vanishes at the phase boundary.
Heuristic justification: We first argue that the KPR vanishes at a dynamical phase boundary that hosts a Krein collision. Consider a bosonic matrix G(p) that depends on some parameters p. Suppose that at some point p c , there is an eigenvalue ω c of G(p c ) that hosts a Krein collision. Then there exist two eigenvectors |ψ c j , j = 1, 2, such that ψ c i |τ 3 |ψ c j = δ ij κ j , with κ j the Krein signature of |ψ c j and κ 1 = −κ 2 . Consider a smooth path p(s), with p(s c ) = p c and, at this point in the argument, let is add the requirement that p c lies on a phase boundary: that is, the system is dynamically stable for, say, s ≤ s c and unstable for s > s c (see also Fig. 2). From Lemma 3.1(ii), one expects generically that ω c should split into a pair of eigenvalues ω(s), ω(s) * related by complex conjugation, for some s > s c , and and that ω c should split into a pair of real eigenvalues ω 1 (s), ω 2 (s) for some s < s c . Now, let |ψ 1 (s) (|ψ 2 (s) ) be the eigenvector of G(p(s)) corresponding to the eigenvalue ω(s) (ω(s) * ) for s > s c and ω 1 (s) (ω 2 (s)) for s < s c . Furthermore, suppose these eigenvectors obey the pseudo-bosonic normalization condition ψ 1 (s)|τ 3 |ψ 2 (s) = 1 for s > s c and the bosonic normalization condition ψ i (s)|τ 3 |ψ j (s) = δ ij κ j for s < s c . Since |ψ 1 (s) and |ψ 2 (s) depend smoothly on s and must satisfy the bosonic normalization condition for s < s c , we can (without loss of generality) assume that |ψ j (s) smoothly connects to |ψ c j at s = s c . The KPR along the path is then where we again enforce the normalization ψ 1 (s) = ψ 2 (s) . Equivalently, r(s) = ψ 1 (s)|τ 3 |ψ 2 (s) , where |ψ j (s) are the normalized eigenvector of G(p(s)). Since |ψ 1 (s) and |ψ 2 (s) smoothly connect to τ 3orthogonal degenerate eigenvectors at s = s c , we have that r(s c ) = 0 at this point. Hence, the KPR should vanish at a phase boundary that hosts a Krein collision. Next, we argue that the KPR should also detect a dynamical phase boundary that hosts an EP, even if G = G T (recall that we already know the claim to be true for G = G T , since r = ρ). The setup is the same as above, but instead of having two linearly independent eigenvectors |ψ c j at s = s c , we now have one normalized eigenvector |χ a and one normalized generalized eigenvector |χ b . Imposing a smooth evolution in parameter space implies that |ψ 1 (s) and |ψ 2 (s) will both approach the eigenvector |χ a as s → s c from above and below. Since |χ a splits into two τ 3 -orthogonal eigenvectors for s > s c , we must have that χ a |τ 3 |χ a = 0 and so r(s) = ψ 1 (s)|τ 3 |ψ 2 (s) will vanish as s → s c from above.

Example.
It is instructive to check our claim for the simplest non-trivial model that may support possible dynamical features. Consider the quadratic single-mode QBH given by The corresponding effective SPH reads with eigenvalues ω ± = ± √ αβ. If sgn (α) = sgn (β), the system can be described a a dynamically but not necessarily thermodynamically stable quantum harmonic oscillator (QHO). For sgn (α) = sgn (β), the system Note that if pc was an isolated point of stability, then it would be possible for |ψ j (s) , with s > sc, to smoothly connect to eigenvectors of null Krein signature at s = sc.  18), with the dynamical phase diagram overlayed. Notice that r(α, β) vanishes between the phase boundaries separating the stable quantum harmonic oscillator (QHO) phase and the unstable parametric amplifier (PA) phases. (b) A plot of the KPR r(α, β) evaluated on the contours β = α n for n = 1, . . . , 6. Notice that the KPR vanishes at α = β = 0 even though the system is diagonalizable along these contours.
is dynamically unstable and equivalent to a degenerate parametric amplifier (PA) (see [48] for an in-depth analysis). The phase boundaries correspond to α = 0 or β = 0. When only one of the parameters is zero, G fails to be diagonalizable and the system can be described as a free particle. By contrast, G is trivially diagonalizable for α = 0 = β because it is the zero matrix. The system can be interpreted as a zero-frequency QHO. This dynamical phase diagram is summarized in Fig. 3. The isolated point α = β = 0 is an example of a dynamical transition point at which the system remains diagonalizable. In Sec. 4 we will investigate a model where these types of phase boundaries extend beyond isolated points.
Let us now investigate the phase boundaries in terms of the KPR. For G = 0, the KPR is r(α, β) = 2 |α||β| |α| + |β| , for both eigenvectors (see again Fig. 3(a)). The KPR vanishes at (α = 0, β = 0) and (α = 0, β = 0) as expected, since these points separate a dynamically stable phase from an unstable one and correspond to EPs of G. The situation at α = 0 = β is more delicate because four dynamical phases meet at the origin and a path through the origin could connect two dynamically (un)stable phases. Our heuristic argument for the behavior of the KPR suggests that the KP R may not be defined at such a point and indeed, the limit lim (α,β)→0 r does not exist. However, we do expect the KPR to vanish on any path through the origin that connects a dynamically stable phase to an unstable one. For concreteness, let β = f (α), with f (α) real analytic at α = 0 and f (0) = 0. Since f (α) = c 1 α + c 2 α 2 + · · ·, the KPR for |α| 1 is well approximated by which can take any value in [0, 1] as long as c 1 = 0. This behavior is exactly as expected according to our heuristic arguments because paths that behave linearly to lowest order near the origin cross from one (un)stable phase to another (un)stable phase. By contrast, paths with c 1 = 0 and even leading order cross from a stable to an unstable phase or, for odd leading order, flatten along the boundary of EPs of G. Either way, the KPR vanishes as it should, signaling that the point α = 0 = β = is indeed contained in a dynamical phase boundary. In Fig. 3(b) we demonstrate this fact by evaluating r(α, β) along the contours β = α n for n = 1, . . . , 6. For n = 1 the KPR is constant while for n > 1 the KPR vanishes, as expected.We will encounter this contour-dependent behavior of the KPR around Krein collisions again in the multi-mode model studied in Sec. 4. It is also instructive to investigate the behavior of the KPR at α, β = 0 from the point of view of the eigenvectors of G(α, β). Since G(0, 0) = 0, any choice of τ 3 -normalized, charge-conjugate vectors can act as the bosonic normal modes and so there is a Krein collision at 0 frequency. The coefficient c 1 in Eq. (20) can be traced back to the normalized eigenvectors of G(α, β = c 1 α), that read Clearly, different choices of c 1 = 0 yield different limiting eigenvectors at α = 0 and, for c 1 → 0, the eigenvectors converge to the τ 3 -null vector associated with the momentum operator p.

Case study: The bosonic Kitaev-Majorana chain
The foundation of our theory of dynamical stability for QBHs lies on the two results from the previous section: (i) Regimes of dynamical stability and instability can be, respectively, understood as unbroken or broken phases of an anti-linear GPT symmetry; and (ii) transitions between two distinct GPT phases can be detected by the KPR introduced in Eq. (17). In this section, we put our theory to work on an interesting model, the BKC of [21]. While this one-dimensional QBH was introduced as a bosonic analogue of the fermionic Kitaev chain, the conventional notion of a quantum phase diagram is not applicable because the system is thermodynamically unstable. The expectation in [21] was that the topology of the effective SPH would control, or impact significantly, the dynamical properties of the model. Our work was partly motivated by a desire to ground this conjecture. We will come back to this issue in Sec. 5.

The model
The QBH we investigate may be written in the form H(s, ϕ) ≡ H O + s W (ϕ), with s ∈ [0, 1] and Specifically, the QBH H O is the bosonic Kitaev-Majorana chain of [21] subject to open BCs. With respect to the general form in Eq. (1), we thus have K ij = it 2 (δ i,j+1 − δ i+1,j ), ∆ ij = i∆ 2 (δ i,j+1 + δ i+1,j ) in the bulk, and K 1N = ist 2 e iϕ = K * N 1 , ∆ 1N = is∆ 2 e iϕ = ∆ N 1 on the boundary. The boundary modification W (ϕ) imposes twisted BCs with a twisting angle ϕ, and the parameter s allows us to smoothly interpolate between open s = 0 and twisted s = 1 BCs (see Fig. 4). The relevant bosonic matrix is G(s, ϕ) ≡ G O + V (s, ϕ), with  with the matrices and, in addition, It is easy to see that H O is not thermodynamically stable. A quick check shows that odd under time reversal, that is, In [21], diagonalization of H(s, ϕ) for t = ∆ was achieved for open BCs (s = 0) thanks to a suitably devised position-dependent local squeezing transformation and for periodic BCs (s = 1, ϕ = 0) by standard momentum-space techniques. Here, we employ the general techniques from Sec. 2.3 to both recover these solutions and analyze different exactly solvable parameter regimes (see also Fig. 5). While we include in Appendix B full detail about the diagonalization of G(s, ϕ) by our approach, we highlight in the following section the salient features of the resulting solutions, as relevant to the subsequent stability analysis.
with N + = N/2 for N even and N + = (N − 1)/2 for N odd, respectively. For N odd there is an additional zero mode, at m = (N + 1)/2, which commutes with the Hamiltonian. The quasi-particle vacuum is The normal form of H O reveals a symmetry of the model embodied in the Bogoliubov transformation, with s m ∈ R arbitrary. The new normal modes are Hence, s m is a free parameter that determines the localization properties of each mode. For the particular choice s m = −j 0 r for all m, we recover the same parametric freedom identified in [21]. At the transition point t = ∆, the normal modes consist of two Jordan chains of the adjoint action of length N and, from the discussion in Sec. 2.2.2, we know that their algebra is fairly arbitrary. The corresponding Jordan chains of G O (t = ∆) can be chosen to be |χ 1k = (it) −k |N + 1 − k |+ and |χ 2k = (−it) −ik |k |− , with k = 1, . . . , N , and map to (multiples of) the Hermitian quadratures These normal modes cannot be combined to yield bosonic normal modes for k = (N +1)/2, as [x k , p N +1−k ] = 0. For N odd and k = (N + 1)/2, one can construct a single bosonic normal mode at zero frequency.

Periodic and anti-periodic boundary conditions. Periodic and anti-periodic BCs (PBCs and APBCs)
correspond to (s = 1, ϕ = 0) and (s = 1, ϕ = π), respectively. Owing to translational invariance, the Hamiltonian can be block-diagonalized the Fourier transform b k ≡ 1 N , x = A, P , labeling the set of wave vectors appropriate for each type of BC, that is, explicitly, Note that the blocks H k=±π/2 , when ±π/2 is in the Brillouin zone (which happens only if N is a multiple of 4 for PBCs), are already in normal form in terms of the bosonic Fourier modes b ±π/2 . The associated eigenfrequencies ω ±π/2 = ±t are the only real ones. Also, recall that the values k ∈ {0, ±π} (when in the spectrum) are precisely the momentum modes that are unpaired in the fermionic Kitaev-Majorana chain. Interestingly, since H k ∝ x k p k + p k x k for k ∈ {0, ±π}, these blocks can be put in a pseudo-bosonic form by letting ψ † k = −ip k and ψ k * = x k . For k ∈ {0, ±π/2, ±π}, the normal form of the blocks is in terms of pseudo-bosons satisfying [ ψ k , ψ † q * ] = δ kq and [ ψ k , ψ q ] = [ ψ k * , ψ q * ] = 0.
For N even, each of the above eigenvalues is doubly degenerate: ω m = ω N +m , for m = 1, . . . , N . Thus, we can take 1 ≤ m ≤ N . Up to a normalization constant, the two eigenvectors for each m are given by m , 1 are generalized Bloch waves (see also Eq. (A.2)), and |1 denoted the first canonical basis vector of C N . For N odd, the spectrum is non-degenerate and the eigenvectors are where again N m is a normalization constant that is chosen so that ψ m |τ 3 |ψ = δ m * , . Given these normalization conditions, one can construct the pseudo-bosonic normal modes just as for periodic and antiperiodic BCs. This exact solution is remarkable because it makes it possible to investigate analytically the flow of eigenvectors as the system approaches the EP that now exists at ϕ = π/2. Let us consider the way in which the pseudo-bosonic modes parametrically evolve into the Jordan chains of generalized normal modes at ϕ = π/2. Focusing for simplicity on N odd, let |χ 11 ≡ |1 |− be one of the two linearly dependent eigenvectors of G T (π/2; t = ∆), and let O m ≡ | χ11|ψm | χ11 ψm be the corresponding fidelity overlap with the eigenvector in Eq. (29). After determining the proper normalization constants N m (Appendix B.2.3), we find the following closed-form expression: In both cases, we have that O m → 1 as ϕ → π/2. Thus, all the eigenvectors coalesce to a single eigenvector at the EP. In particular, they all become perfectly localized at site j = 1 as the system approaches the EP. Before moving to an in-depth analysis of the stability properties of the chain, we briefly comment on the induced many-body action of the GPT symmetry as a function of BCs. The stable BCs (e.g., open and π/2twisted) support an unbroken GPT symmetry that (anti-linearly) maps each of the bosonic quasi-particle creation and annihilation operator to itself. For the unstable cases (e.g., periodic and anti-periodic), the GPT symmetry (anti-linearly) maps the normal modes to their pseudo-bosonic partners (e.g., z ψ k → z * ψ k * , for any z ∈ C). In particular, we observe that this symmetry is a function of the BCs.

Stability analysis
4.3.1. Stability phase diagram as a function of BCs. As we saw in terms of closed-form solutions, the bosonic Kitaev-Majorana chain of Eq. (21) displays both stable and unstable dynamical phases. There is a GPT symmetry that breaks at any stability-to-instability transition in one of two non-exclusive ways: either because the effective SPH loses diagonalizability, or because it develops Krein collisions that will generically cause a real eigenvalue to split into a pair of complex-conjugate eigenvalues. For for both open or π/2-twisted BCs, the onset of instability at t = ∆, is of the first type, due to loss of diagonalizability. There are also, however, dynamical phase transitions induced by changes in the the boundary parameters s and ϕ for t > ∆. In this section, we investigate both analytically and numerically these transitions. Figures 5(b) and 6(a)-(d) show the dynamical phase diagrams for various choices of parameters, obtained by numerically determining the eigenvalue spectrum for different system size N . From Lemma 3.1, we know that the open chain can undergo a dynamical phase transition for arbitrarily small perturbations, due to Krein collisions. While Fig. 5(b) suggests that the region of dynamical stability which connects OBCs to π/2-twisted BCs may be very nearly a line (hence, of zero measure), additional analysis for a smaller ∆/t ratio reveals that the thickness of the region surrounding the line ϕ = π/2 does not vanish for finite lattice size. That is, for each s, the chain is stable for ϕ ∈ [π/2 − δϕ N (s)/2, π/2 + δϕ N (s)/2] ≡ I ϕ N (s). For finite N , δϕ N (s) > 0 strictly, and the minimum width of the stability region is given by δϕ N ≡ δϕ N (s = 1). Analogously, for each ϕ, N , we let δs N (ϕ) be such that the system is dynamically stable for s ∈ [0, δs N (ϕ)] ≡ I s N (ϕ) and define the minimum height of the stability region by δs N ≡ δs N (ϕ = 0).
We can characterize both the way in which phase boundaries depend upon parameters and their nature in terms of Krein collisions versus EPs by combining analytical and numerical techniques. To this end, it is useful to analyze the spectral flow through the transitions (shown in Fig. 7 for the same t, ∆ used in Fig. 6) as well as the presence of the zero eigenvalue in the spectrum, which we track by computing minimum-modulus eigenvalue of G(s, ϕ) (see Figs. 6(e)-(h)). Four main points are worth noticing: (ii) For N odd, Fig. 7(a2) further indicates that the dynamical phase transition between the π/2-twisted chain and the periodic chain at s = 1 occurs when eigenvalues of opposite Krein signature become arbitrarily close.
(iii) For N even, Fig. 7(b) reveals that the transition away from stability happens when pairs of eigenvalues with opposite Krein signature become arbitrarily close at finite s > 0.
(iv) The value zero enters the spectrum of G(s, ϕ) at the left dynamical phase boundary for arbitrary N , (in fact, at both boundaries if N is odd), and the phase diagram is symmetric about ϕ = π/2. Interestingly, further numerics (data not shown) indicate that, for N a multiple of 4, the right phase boundary is still defined as the locus in parameter space where the largest (in modulus) eigenvalue of G(s, ϕ) is closest to zero as a function of s, ϕ.
Observation (i) demonstrates explicitly the behavior predicted by Lemma 3.1, namely, stability-toinstability transitions are generically mediated by the coalescence of eigenvalues with opposite Krein Going down the left column (a1), we show how the eigenvalues evolve within the stable phase, as s goes from 0 to 1 at ϕ = π/2. Recall that, for N odd, ϕ = π/2, and s = 1, each eigenvalue except the extremal ones are twofold degenerate. Stability is seen to be preserved along this flow. Going down the right column (a2), we show how that eigenvalues evolves around the transition between π/2-twisted and periodic, as ϕ goes from π/2 to 1.03 · (π/2) at s = 1. Note how eigenvalues are each split as ϕ increases, and eventually move symmetrically off the real axis. (b) N = 10. The spectral flow around the transition between open and periodic, as s goes from 0 to 0.1 with ϕ = 0. In this case, stability is retained for sufficiently small strength of the boundary perturbation due to non-zero s. signature. However, while this behavior is generic, it is not universal. Observations (ii) and (iii) together highlight an interesting even-odd effect, namely, the minimum stability height is only non-zero when N is even. Finally, observation (iv) provides evidence supporting the following conjecture: Conjecture 1: Generically, bosonic zero modes indicate dynamical phase boundaries.
While, as Fig. 6 illustrates, the right phase boundary for even N provides an example where no zero modes are hosted, we also see that if the two boundaries coalesce as N grows, then, in the thermodynamic limit, every dynamical phase transition of the BKC is defined by zero modes. We now show that this does, in fact, happen. We do so by utilizing the techniques of Sec. 2.3 to determine the values of s and ϕ that support zero modes. Combining calculations that are detailed in Appendix B.2.2 with the symmetry of the dynamical phase diagram about ϕ = π/2, we find that the phase boundaries are parameterized by the curves cos(ϕ) = ± 1 2 (s + s −1 )sech(N r), N even, 2 sech(N r), N odd, which are shown in Fig. 6(a)-(d) are white, dashed lines. In particular, in term of the parameter 2r = ln[(t + ∆)/(t − ∆)] previously introduced, the minimum height and width, which we introduced earlier on for characterizing the stable phase, are respectively determined by the following expressions: Both quantities decrease rapidly with system size, confirming that that stable regions become of measure zero in parameter space in the thermodynamic limit. Thus, the transition lines in this case are precisely defined by the occurrence of zero modes, suggesting that Conjecture 1 may hold for all QBHs in the thermodynamic limit. Similarly, for fixed N , both δs N and δϕ N decrease as ∆ approaches t, signaling the transition to the unstable ∆ > t phase. We note that in a recent paper [50], a quantity analogous to δs N was calculated for the non-Hermitian fermionic Hatano-Nelson chain. As shown in [21], the effective SPH of the BKC (without twisting) is unitarily equivalent to two copies of the Hatano-Nelson chain. Since twisting the BCs manifestly destroys this unitary equivalence, it is remarkable that such strong similarities still exist. As we explicitly show in Appendix B.2.2, the zero eigenvalue hosts one (two) Jordan chain(s) of length two on the (left) phase boundaries for N odd (even). This information allows us to better assess whether the rest of the spectrum hosts non-trivial Jordan chains or simple Krein collisions. For N odd, we can do this numerically, by calculating the distance between the two eigenvectors (after correcting for arbitrary phases) that coalesce at each eigenvalue. We find that for odd system sizes between N = 5 up to N = 55, various choices of t/∆ ∈ (0, 1), and at various points along the phase boundary, these distances vanish, indicating a loss of linear independence along phase boundaries at s > 0. Thus, each eigenvalue hosts a Jordan chain of length two. For N even, things are more complicated. As seen from Fig. 7(b), the spectrum is always at least doubly degenerate, and splittings occur at different points in parameter space (the splitting at zero defining the phase boundary). In contrast, each eigenvalue splits simultaneously, at the same values of s and ϕ, for N odd. Due to these complications, we do not further assess the nature of these splittings; we nonetheless conjecture that, like for zero frequency, they are induced by the formation of two length-two Jordan chains.
Aside from the specific BKC example and the single-mode model considered in Sec. 3, Conjecture 1 is further motivated from a physical point of view. To see this, consider perturbing a QBH of the form Eq. (1) by adding a small term linear in the creation and annihilation operators (e.g., arising from a constant force). Intuition suggests that such a term may only elicit dynamical instability from a dynamically stable system if the latter possesses zero modes (e.g., a harmonic oscillator exhibits bounded motion even in a constant gravitational field, whereas a free particle or zero-frequency oscillator does not). Further, it can be shown that a QBH can "absorb" any such linear perturbation, meaning that the total Hamiltonian can be shown to be unitarily equivalent to a purely QBH, unless it possesses zero modes. In this sense, zero modes are generically expected to define dynamical phase boundaries.
There is an intriguing point of contact between the BKC for t > ∆ and the single-mode model of Sec. 3. It appears that the only phase boundary of the BKC that hosts (a macroscopic number of) Krein collisions is the line s = 0, which physically represents a single (open) BC, hence a single point in parameter space. Similarly, there is only one point in the phase boundary of the single-mode model that hosts a Krein collision, the origin. In contrast, the one-dimensional phase boundaries in both models are dominated by EPs. These observations suggest a second conjecture: Conjecture 2: Generically, the (d − 1)-dimensional phase boundaries of the d-dimensional dynamical phase diagram of a QBH are characterized completely by EPs, whereas the (d − 2)-dimensional boundaries are characterized completely by Krein collisions.
We conclude this section by investigating the spectral behavior near the dynamical phase transition as a function of system size. In Fig. 8(a), we show the spectrum of G(1, ϕ) in the vicinity of the phase boundary around ϕ = π/2. We see that as N increases, the spectrum clings more strongly to the ellipse defined by the periodic/anti-periodic spectrum. Thus, the "speed" at which the spectrum splits from the real axis increases dramatically as N increases. This may be quantified by examining the spectral speed, d|ω m |/dϕ, for which we conjecture that that lim N →∞ d|ω m |/dϕ ∝ δ(ϕ − π/2). Although we cannot evaluate this quantity analytically for t = ∆, we can use the exact analytical solution at t = ∆, s = 1 and ϕ ∈ [0, π] as a point of Figure 8. (a) The spectrum with t = 1, ∆ = 0.5, s = 1, and ϕ = 0.99 · (π/2) for system sizes, N = 10, 20, and 100, from inner to outer. The solid outermost ellipse traces out the periodic spectrum for N → ∞. (b) The spectral speed for same t and ∆ for various N , averaged over all eigenvalues.
In particular, we see in this case that lim N →∞ d|ω m |/dϕ ∝ δ(ϕ − π/2), as conjectured. Thus, in the thermodynamic limit, this quantity contains an extreme non-analyticity, as we also illustrate in Fig. 8(b).

4.3.2.
Krein phase rigidity. In Sec. 3.2, we argued that the Krein phase rigidity of Eq. (17) should be able to detect both EPs and Krein collisions at a boundary between a stable and an unstable dynamical phase.
Here we demonstrate this capability in the context of the BKC by evaluating numerically the KPR of a representative eigenvector as a function of s and ϕ, see Fig. 9. We confirm that the KPR does, in fact, vanish at the phase boundaries for s > 0. In order to understand the behavior of the KPR at s = 0, we again point out that this corresponds to precisely one BC (open). Thus, the limiting value of the KPR at s = 0 is contour-dependent as it was in the single-mode model of Sec. 3.2.3. In particular, we emphasize that the KPR evaluated along any contour in parameter space whose transition from instability (ϕ ∈ I ϕ N ) to stability (ϕ ∈ I ϕ N ) is mediated by the point s = 0 will vanish at this point. This can be seen in Fig. 9(c), where the KPR is evaluated along a parabolic contour which passes through the point s = 0 at precisely the twisting angles defining the dynamical phase boundaries and along which the effective SPH G(s, ϕ) remains diagonalizable. As predicted in Sec. 3, the KPR detects these Krein-collision-dominated dynamical phase transitions, despite the lack of EPs. To illustrate the response of the KPR to system size, let us focus on twisted BCs, s = 1 and ϕ arbitrary, see Fig. 10(a). As N increases, the KPR r(ϕ) approaches a continuous curve, which however sharply detects the transition point ϕ = π/2, akin to a "stability order parameter". While this relatively tame behavior contrasts with the extreme sensitivity of the spectrum to system size, it is interesting to note that a similar, "less extreme" response to increase in system size has been reported for the eigenvectors of non-Hermitian asymmetric hopping models, through studies of fidelity decay and Loschmidt echo [51]. We can compare the numerical results of Fig. 10(a) for t = 1 and ∆ = 0.5 with analytical results available for t = ∆. In particular, our exact analytical solution for G(1, ϕ; t = ∆) allows us to investigate the KPR in the vicinity of the EP at ϕ = π/2. By taking N to be odd and using the eigenvectors from Eq. (29), we find A plot of r m (ϕ) for various system sizes N is given in Fig. 10(b). In the limit as N → ∞, lim N →∞ r m (ϕ) = | cos ϕ| ln(| cos ϕ|) | cos ϕ| 2 − 1 .
Since r m (ϕ) vanishes as | cos ϕ| → 0, the EP at s = 1, ϕ = π/2 is indeed detected by the KPR. d dt This feature is referred to as "phase-dependent chiral transport" in [21]. In this context, the word "chiral" is meant to further highlight the asymmetric way in which each quadrature is influenced by couplings between adjacent lattice sites, with maximal asymmetry and uni-directional transport being approached as t → ∆. Considering more general BCs, our analysis reveals that although the features characteristic of phasedependent chiral transport survive for arbitrary s ∈ [0, 1] and ϕ = 0, they are fragile from a dynamical perspective, since the system is either dynamically stable with Krein collisions in the spectrum or dynamically unstable. Is it possible to have phase-dependent chiral transport in chains with more favorable stability properties? Unfortunately, the answer is in the negative. The decoupling of the equations of motion is the key prerequisite for both phase-dependent and chiral transport, and this prerequisite condition forces the spectrum of G to host Krein collisions, assuming it is not already unstable. More formally: Proposition 5.1 Let H be a dynamically stable QBH that supports phase-dependent transport. Then, the spectrum of the effective SPH G necessarily hosts Krein collisions. In other words, the system is at the cusp of instability.
From a many-body perspective, if the matrices K and ∆ have purely imaginary entries (the "decoupling condition"), then the Hamiltonian H is, like the GKC H O , necessarily odd under time reversal. Thus, H cannot be thermodynamically stable and even if it is dynamically stable, this feature is necessarily fragile.

Interplay with topology
Looking for hints of topological physics in systems of free bosons, several researchers have developed mappings from free fermions to free bosons that preserve specific features of interest. For example, in our work in [20] we explore a mapping from fermions to bosons that preserve zero-energy modes. The BKC of the previous section is the result of another such mappings; specifically, by design, one that preserves a topological invariant, the winding number. To put things in context, we state here the mapping of [21] in greater generality. Let denote a general quadratic fermionic Hamiltonian [1]. The mapping of interest is restricted to the subclass of Hamiltonians obeying K T = K and ∆ † = ∆. Altogether, these additional conditions imply that K is purely real and ∆ is purely imaginary. The output of the mapping is the QBH that is, with respect to Eq. (1), we get Hence, according to Proposition 5.1, these bosonic counterparts of QFHs are either dynamically unstable or at the cusp of instability. This particular mapping from free fermions to free bosons may preserve topological invariants at the expense of sacrificing robust dynamical stability.
This conclusion is tightly linked the factor of i that appears in front of K in Eq. (33). Interestingly, the deleterious factor of i is introduced precisely so that the winding number of the fermionic Kitaev chain survives the passage to bosons. For periodic BCs and zero chemical potential, the Hamiltonian of the fermionic Kitaev chain may be written in momentum space as where BZ denotes the appropriate Brillouin zone,Ψ k = [c k , c † −k ] T , and d f (k) = [0, −∆ sin(k), t cos(k)] T . The SPHs H = d f (k)·σ has a well-defined topological invariant -namely, the number of times d f (k) wraps around the origin in the yz-plane -which is protected by a chiral symmetry. Implementing the above mapping, the bulk Hamiltonian for the corresponding QBH is given by d b (k) · σ where d b (k) = [0, −∆ cos(k), t sin(k)] T In real space, this is precisely the Hamiltonian of Eq. (21) (for PBCs). The vector d b (k) · σ winds around the origin just like the original fermionic one does, but only because of the factor of i in front of K in Eq. (33).
Taking all of this into consideration, it seems both peculiar and suggestive that the BKC is dynamically unstable for periodic BCs and stable, but at the cusp of instability, for OBCs. Is this feature somehow related to the winding number? The following analysis suggest the the answer is in the negative. Consider again the Kitaev chain, but now including a non-zero chemical potential, where µ, t, ∆ > 0 and now d f (k) = [0, −∆ sin(k), µ + t cos(k)] T . The winding number is non-zero for µ < t in the topologically non-trivial phase of the superconductor. In real space and for open BCs the associated bosonic Hamiltonian is now where H O is the BKC of Eq. (21). Therefore, the chemical potential term in the fermionic chain maps to a sum of degenerate parametric amplifier terms in the bosonic chain. These additional terms modify the bosonic effective SPH as G O → G O + iµτ 1 , causing the originally doubly degenerate eigenvalues ω m of G O to split into ω m ± iµ, with ω m = √ t 2 − ∆ 2 cos(mπ/(N + 1)), m = 1, . . . , N . Thus, for µ = 0, the system is always dynamically unstable. And, this conclusion holds irrespective of the condition µ < t required for a non-zero winding invariant.
At this point it would seem that we have managed to discount a direct topological origin to stability properties of QBH. However, as we noted in the Introduction, the role of topology for QBHs remains, at best, only partially understood as yet. The highly non-trivial nature of this interplay is nicely illustrated by two further observations pertaining to the comparison between the Kitaev chain and its bosonic counterpart. First, our analysis establishes that the BCs that render the BKC dynamically stable for all N and t > ∆ (open and π/2-twisted BCs) are precisely the same that host Majorana zero modes in the corresponding fermionic chain [49]. Second, the Hamiltonian H b of Eq. (34) displays localized approximate zero modes for t = ∆ and every finite N . Specifically, one can verify that the following left-and right-localized modes These Hermitian "Majorana bosons" at zero frequency are normalizable, that is, exponentially localized precisely if µ < t -in perfect correspondence with the topologically-sourced Majorana zero modes of the fermionic chain. Captivating as they are, we do not know as yet whether these "bosonic shadows" of Majorana physics are topological in any suitable sense of the word.

Conclusion and outlook
We have systematically investigated the landscape of free-boson dynamical phase diagrams from a general dynamical stability perspective and by way of two paradigmatic examples: a single bosonic mode and a bosonic version of the Kitaev-Majorana chain, as introduced in [21]. Our general framework for QBHs combines tools from pseudo-Hermitian, PT -symmetric, and non-Hermitian quantum mechanics with the Krein stability theory of dynamical systems in indefinite inner-product spaces. Two key new results emerge from this analysis: First, all free-boson systems are PT -symmetric in a suitable sense and their dynamical phase diagrams are controlled by the fate of this symmetry. From a many-body standpoint, this symmetry is broken precisely when the QBH can no longer be diagonalized in terms of canonically bosonic Bogoliubov quasi-particles. Second, we argued that dynamical phase boundaries can be detected by a KPR indicator, which naturally extends the notion of phase rigidity widely employed within semiclassical (non-Hermitian) treatments of open quantum systems. Bosonic dynamical phase boundaries can consist of loci of exceptional points, where diagonalizability of the effective SPH is lost, but also, remarkably, or Krein collisions where degenerate real eigenvalues split into the complex plane without loss of diagonalizability.
To illustrate and validate our framework, we obtained a complete characterization of the dynamical stability phase diagram of the BKC for a two-parameter family of BCs that interpolates between open and periodic BCs and includes twisted BCs as notable case. In particular, for both open and π/2-twisted BCs we were able to diagonalize in closed form the BKC, by employing for the first time an exact diagonalization procedure developed in the context of free fermions as a generalization of Bloch's theorem to systems where translational symmetry is broken by BCs. In particular, the use of this procedure proved instrumental to access to study the stability phase diagram as a function of system size, supporting the emergence of extreme non-analyticity of the spectral response in the thermodynamic limit. We confirmed explicitly that the KPR vanishes at all the phase boundaries, as expected from our general arguments, and sharply detects GPT symmetry-breaking phase transitions in the thermodynamic limit. Remarkably, our analytical solutions prove that the BCs that host Majorana zero modes in the fermionic Kitaev-Majorana chain are precisely the same that allow for dynamical stability in its bosonic counterpart.
One of the most interesting physical properties of the BKC is that it can support phase-dependent chiral transport, stemming from the decoupling of the evolution of the real and imaginary parts of coherent excitations. Using tools from Krein stability theory, we showed that any QBH exhibiting a similar "decoupling condition" describes a thermodynamically unstable system, which is either dynamically unstable or sitting at the cusp of dynamical instability. As a consequence, unless some additional protection mechanism is in place, stable phase-dependent transport is fragile against perturbations. Finally, the BKC is the result of applying a certain mapping to the fermionic chain. By exploring in more generality the idea of mappings between fermionic and bosonic systems, so that a specified set of topological invariants is preserved, we have shown that while topology may not directly influence the dynamical phase of a QBH, bosonic analogues to Majorana zero modes exist in the topological phase of a further generalized BKC model. Achieving a clearer picture of what role topology may play (if any) in informing the dynamical properties of bosonic systems, along with an explanation of whether the seemingly special status of certain BCs is simply a coincidence or rather has a deeper significance, are natural next questions we leave to future research.
Our analysis also points to a number of additional directions for investigation. On the one hand, for closed systems of bosons, it is worth to explore the role of KPR and GPT symmetry beyond the mean-field approximation. For example, a question of fundamental relevance would be to determine whether the KPR may be used to assess the validity of the underlying quadratic approximations for interacting systems, both at equilibrium or possibly under a time-dependent driving. Furthermore, is there a role for Krein stability theory or GPT symmetry-breaking to play in the full Fock space of systems of interacting bosons? On the other hand, our framework and tools can be extended to a large class of quadratic (fermionic or bosonic) open systems -either described, semi-classically, in terms of non-Hermitian many-body effective Hamiltonians or, within a fully quantum formalism, by Lindblad (Markovian) master equations [52]. In particular, we expect that the KPR may find natural applications in the context of exploring topological phenomena related to PT -symmetric quantum quenches [53] or dynamically encircling EPs [54], or in the context of PT symmetry-breaking enhanced quantum metrology [55]. Likewise, while for quadratic Lindbladians there exist (generically) non-Hermitian matrix analogues to the effective SPH [56,57], important differences are also to be expected and, indeed, have been recently pointed out for instance in the nature of the underlying EPs [58]. One of our next steps will thus be to understand the extent to which the KPR (or some suitable modification of it) may still provide a useful diagnostic tool in open quantum dynamical settings and, if so, be employed to characterize steady-state stability phase diagram in driven-dissipative many-body quantum systems, including the possible emergence of topological features or exotic transport phenomena [59,60].
simultaneous eigenvectors of T and T −1 are given by Φ z,1 ≡ j∈Z z j |j , where z is an arbitrary, non-zero complex number. Explicitly, T Φ z,1 = zΦ z,1 and T −1 Φ z,1 = z −1 Φ z,1 , which immediately lead to the identity where |u ∈ C 2 is arbitrary. We call G(z) the reduced bulk effective Hamiltonian and note that G(z = e ik ) is the usual Bloch Hamiltonian that arises in 1D systems under Born-von-Karman (periodic) BCs. Thus, G(z) is the analytic continuation of G(e ik ) off the unit circle. Furthermore, we see that for any z = 0 such that G(z) |u = ω |u , Φ z,1 |u is an eigenvector of G with eigenvalue ω.
To continue, we define the complex characteristic polynomial P (ω, z) ≡ z 4R det(H(z) − ω1 2 ). We call an eigenvalue ω regular if P (ω, z) is not the zero polynomial. Otherwise, we say ω is singular. For the applications in this paper, it suffices to restrict to the eigenvalues ω that are regular. For a fixed ω, let {z } n =1 denote the n distinct roots of P (ω, z) and {s } n =1 denote their corresponding multiplicities. Generically, G(z ) will have s eigenvectors {|u s } s s=1 satisfying G(z ) |u s = ω |u s , in which case, the vectors are solutions to the bulk equation, and akin to Bloch waves with complex momentum. When the reduced bulk Hamiltonian G(z ) has less than s eigenvectors, the remaining solutions are constructed from the generalized eigenvectors of the left and right translation operators. The sequences span the kernel of (T − z) s for ν = 1, . . . , s. Furthermore, One can then show that the sequence Ψ ≡ ν n=1 Φ z,n |u n satisfy where G ν (z) is an upper-triangular block-Toeplitz matrix with non-zero blocks It can then be shown that the eigenspace of G corresponding to eigenvalue ω is a direct sum of n vector spaces spanned by generalized eigenvectors of T ±1 of the form where the linearly independent vectors {u sν } are chosen in such a way that G s (z ) |u s = ω |u s with |u s = [|u s1 , . . . , |u ss ] T . With these, we obtain n =1 s solutions to the bulk equation given by |z , ν |u sν , |z , ν = P 1,N Φ z,ν .
If g ±R are not invertible, then there exists 2s 0 ≡ 4R− n =1 s additional boundary localized solutions to the bulk equation, where s 0 is the multiplicity of z = 0 as a root of the characteristic polynomial P (ω, z) for a given regular eigenvalue ω. We will now demonstrate how to construct the left (j = 1) localized solutions. Since these solutions emerge due to the truncation of the bi-infinite lattice to a finite one, we consider the half-infinite auxiliary effective Hamiltonian and unilateral shift operators |j + 1 j| .
Appendix B.2.2. Dynamical phase boundaries. In this section, we determine analytically the dynamical phase boundaries in boundary parameter space. An important assumption of this derivation is that certain phase boundaries are characterized by the emergence of zero modes and that the phase diagram is symmetric about ϕ = π/2. Thus, we will uncover the conditions on s and ϕ for G(s, ϕ) to possesses zero as an eigenvalue. As in the preceding Appendix, we will rotate via the unitary U and study the unitarily equivalent matrix G (s, ϕ). In contrast to the preceding section, however, we keep ϕ arbitrary and restrict to the non-open case s ∈ (0, 1]. Since the bulk (G O ) is unchanged, and the roots of the characteristic polynomial P (ω = 0, z) are distinct, we have the same four bulk solutions |ψ j , j = 1, 2, 3, 4, given in Eqns. (B.1)-(B.6). On the other hand, the boundary modification is now given by f cos(ϕ) f sin(ϕ) −J sin(ϕ) J cos(ϕ) = σ y v † −1 σ y .
Demanding that the determinant vanishes, we obtain the conditions cos(ϕ ± ) = 1 2 (s + s −1 )sech(N r), N even, ±2 sech(N r), N odd. (B.8) For N even, this specifies one angle ϕ + = ϕ − in the interval [0, π], in fact, smaller than π/2. On the other hand, for N odd, there are two distinct angles ϕ ± symmetric about each side of π/2. Thus, both phase boundaries host zero modes for N odd and just the left boundary for N even. When Eq. (B.8) is satisfied, the kernel of B(0) can be determined analytically. The cases s = 1 and s = 1 must be handled separately. We begin by taking s = 1. For N even, ker B(0) is two-dimensional and spanned by the vectors The important thing to note is that these calculations reveal that the dimension of the zero-mode subspace is one (two) for N odd (even). The four-fold symmetry of the spectra of bosonic effective SPHs implies that the algebraic multiplicity of the zero eigenvalue must always be even. This confirms that for N odd, there must be a Jordan chain of length two at zero, along the phase boundaries (s > 0). An additional symmetry of the even chain implies that all non-zero eigenvalues of G(s, ϕ) are at least doubly degenerate, implying that the zero eigenvalue has algebraic multiplicity four. Thus, the even chain possesses two lengthtwo Jordan chains at zero, along the left phase boundary. Alternatively, this can be concluded by checking that the dimension of kernel of the boundary matrix of G 2 at zero frequency is four.