Multiple quantum exceptional, diabolical, and hybrid points in multimode bosonic systems: I. Inherited and genuine singularities

Multiple quantum exceptional, diabolical, and hybrid points in multimode bosonic systems: I. Inherited and genuine singularities Kishore Thapliyal1, Jan Peˇ rina Jr.1, Grzegorz Chimczak2, Anna Kowalewska-Kud laszyk2, and Adam Miranowicz2 1Joint Laboratory of Optics, Faculty of Science, Palack´ y University, Czech Republic, 17. listopadu 12, 771 46 Olomouc, Czech Republic 2Institute of Spintronics and Quantum Information, Faculty of Physics, Adam Mickiewicz University, 61-614 Pozna´ n, Poland The existence and degeneracies of quantum exceptional, diabolical, and hybrid (i.e., diabolically degenerated exceptional) singularities of simple bosonic systems composed of up to five modes with damping and/or amplification are analyzed. Their dynamics governed by quadratic non-Hermitian Hamiltonians is followed using the Heisenberg-Langevin equations. Their dynamical matrices generally exhibit specific structures that allow for an effective reduction of their dimension by half. This facilitates analytical treatment and enables efficient spectral analysis based on characteristic second-order diabolical degeneracies. Conditions for the observation of inherited quantum hybrid points, observed directly in the dynamics of field operators, having up to third-order exceptional and second-order diabolical degeneracies are revealed. Surprisingly, exceptional degeneracies of only second and third orders are revealed, even though the systems with up to five modes are considered. Exceptional and diabolical genuine points and their degeneracies observed in the dynamics of second-order field-operator moments are also analyzed. Each analyzed bosonic system exhibits its own unique and complex dynamical behavior. 1 Introduction Non-Hermitian Hamiltonians had been for a long time considered not being suitable for describing real physical systems. This opinion has changed after the seminal work by Bender and Boettcher [1] who showed that the non-Hermitian Hamiltonians endowed with a parity and time symmetry (PT-symmetry) exhibit real spectra in certain areas of the system parameter space. This leads to the formulation of new area of physics, i.e., non-Hermitian quantum mechanics [2–5], which has already provided numerous models [6–11] suitable for describing real quantum systems in many areas of physics. Moreover, non-Hermitian Hamiltonians exhibit new algebraic structures: It has been shown that, for certain values of parameters, there occur spectral degeneracies accompanied by degeneracies of the eigenvectors of a given Hamiltonian. Such points in a system parameter Kishore Thapliyal: [email protected] Jan Peˇ rina Jr.: [email protected] Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 1 arXiv:2405.01666v2 [quant-ph] 8 Dec 2025 space with exceptional degeneracies (EDs) are called exceptional points (EPs) for which the dimension of the corresponding Hilbert space is reduced. This has interesting physical consequences and leads to new unexpected physical effects. It allows to enhance the precision of measurements of suitable physical quantities [12–16]. It also leads to enhanced nonlinear interactions [17,18]. For this reason, PT-symmetric systems of different kinds have been found appealing in many areas of physics including: optical waveguides [19,20], optical lattices [21–24], spin lasers [25], optical coupled structures [26–30], coupled optical microresonators [13,31–35], quantum-electrodynamics circuits (QED) [36], systems with complex potentials [37], optomechanical systems [38,39], photonics molecules [40], among others. Moreover, schemes for engineering properties of EPs have been suggested (see, e.g., [41] and references therein). Subsequent studies have revealed also other features observed in PT -symmetric systems in addition to EPs. For example, the spectral degeneracies which are not accompanied by the corresponding eigenvector degeneracies, were observed [42]. Such degeneracies do not usually lead to the above discussed physical effects and so the corresponding points in the system parameter space were named diabolical points (DPs). Note that DPs, contrary to EPs, can also be observed in Hermitian systems [42]. As pointed out in [43], it may happen that this DP with its diabolical degeneracy (DD) occurs at an EP. In this case, ‘independent’ (i.e. with different eigenvectors) multiple spectral and eigenvector degeneracies are found in systems and we refer to them as hybrid diabolic exceptional points (HPs). Such systems then naturally exhibit a physical behavior similar to that observed at EPs. Moreover, new effects originating in diabolical degeneracy may arise. For example, the system behavior when encircling an HP has been used to construct a multi-mode optical switch [44]. Non-Hermitian PT -symmetric optical bosonic systems are especially interesting from the point of view of their behavior at EPs. Their infinite-dimensional Hilbert space leads to numerous manifestations of modifications of their dynamics at quantum EPs (QEPs), i.e. EPs observed in quantum systems [18], including the effect of quantum jumps [45,46]. The system dynamics may be followed either directly in the Hilbert space [the Liouville space of statistical operators] or in the complementary space of field-operator moments (FOMs) of all orders [47,48]. The studies performed in the moment space of field operators in multimode bosonic systems described by non-Hermitian quadratic Hamiltonians revealed different types of degeneracies related to QEPs [43,49]. Moreover, they allowed to sort QEPs into three classes: inherited QEPs, genuine QEPs and induced QEPs. The dynamical equations for the mean values of field operators indicated the presence of inherited QEPs and quantum HPs (QHPs) [43] that represent the core of the studied unusual behavior. The presence of such inherited QEPs and QHPs then implies the existence of genuine QEPs and QHPs [43] observed in the dynamics of higher-order FOMs. With the increasing FOM order, the degeneracies of genuine QEPs and QHPs increase. Moreover, similar or identical FOMs, as being related by the field commutation relations, arise in the formal construction of higher-order FOMs. Thus, we can also define induced QEPs and QHPs [43] that further enlarge the multiplicity of the spectral degeneracies. However, as these redundant FOMs share their time evolution with the FOMs contributing to genuine QEPs and QHPs, they do not lead to additional diversity of the system evolution. Thus, they are not interesting from the point of view of the dynamics of FOMs of a given order. This dynamics is fully characterized by the corresponding genuine QEPs and QHPs. As the properties of genuine QEPs and QHPs originate in those of the inherited QEPs and QHPs, the analysis of the latter is crucial for the understanding of a system evolution. For this reason, it is important to identify the inherited QEPs and QHPs and their deAccepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 2 generacies in simple bosonic systems formed by smaller numbers of bosonic modes. This analysis may then be exploited for further studies of physical effects in such systems. The system composed of two mutually interacting modes, one being damped and the other amplified, was already analyzed from this point of view in Ref. [43]. This analysis may be considered as the simplest building block useful for investigations of more complex bosonic systems in which much richer and diverse physical behavior is expected. Here, we consider systems composed of up to five modes in different configurations promising for the observation of QEPs and QHPs. Looking for inherited QEPs with higher-order EDs is the main goal of our investigations. It is motivated by the fact that the higher is the ED order, the more modified is the system dynamics at a QEP. This then enhances the physical effects specific to QEPs like improvement in measurement precision or enhancement of nonlinear effects. As, surprisingly, we have been able to reveal only up to the third-order ED of QEPs in the analyzed systems, we continue our analysis in Ref. [50], being the second part of this paper, in which we pay attention to the spectral degeneracies of non-Hermitian Hamiltonians observed only in specific subspaces of the systems’ Liouville spaces as well as non-Hermitian Hamiltonians with unidirectional coupling. In particular, systems with unidirectional coupling exhibit significantly stronger non-Hermitian features compared to those analyzed here, enabling the observation of QEPs and QHPs of arbitrary order. In Ref. [50], we also address numerical identification of QEPs and QHPs, useful in the bosonic systems with higher dimensions. We also extend the analysis of the genuine and induced QEPs and QHPs by considering the FOMs of a general order. At the end of Introduction, we would like to note that while the existence of higherorder EPs has been discussed previously in the literature, those analyses typically rely on semiclassical models that neglect quantum fluctuations, that are necessary for quantum consistency of the models. Moreover, many prior works investigate the spectra of semiclassical non-Hermitian Hamiltonians, which exhibit infinite (though countable) sets of eigenfrequencies. In contrast, the physically relevant spectra are those of the corresponding dynamical matrices governing the evolution equations (e.g., the Heisenberg ones in the Heisenberg picture), which involve only a small number of frequencies. It is precisely these frequencies that determine the observable behavior of the system. In semiclassical models, higher-order EPs can often be identified with relative ease. However, the validity of such models is generally restricted to specific conditions, such as short evolution times or weak damping/amplification. This raises a critical and largely unexplored question: What are the spectral degeneracies in PT-symmetric bosonic systems under general conditions, where fluctuating quantum forces - governed by fluctuation-dissipation theorems - play a non-negligible role? The paper is organized as follows. A two-mode bosonic system with unequal damping and/or amplification rates, as the simplest considered model, is analyzed in Sec. II. Section III brings the analysis of a related three-mode system in the linear configuration. The corresponding generalized four-mode systems in their linear and circular configurations are investigated in Sec. IV whereas the analysis of the five-mode systems in their linear and pyramid configurations is found in Sec. V. Section VI brings conclusions. In Appendix A, the eigenvalues and eigenvectors of the general 2n×2ndynamical matrices are summarized. The structure of the second-order FOMs for the analyzed systems is detailed in the tables provided in Appendix B. Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 3 2 Two-mode bosonic system: Basic building blocks We begin with the consideration of one of the simplest PT-symmetric bosonic systems that is composed of two modes: one being attenuated and the other amplified [for the scheme, see Fig. 1(a)]. We note that even a one-mode bosonic system may exhibit PT -symmetric behavior, as shown in Ref. [51]. We also pay attention only to the systems described by quadratic Hamiltonians that lead to linear exactly-solvable Heisenberg equations. Though these Hamiltonians lead to linear dynamical equations, they allow for describing the nonlinear effect of photon-pair generation and annihilation. This effect is commonly used in various quantum optical systems to generate entangled [52] and squeezed [53,54] states of light. PT -symmetry restricts the form of the studied non-Hermitian Hamiltonian such that the underlying dynamics is described by the dynamical matrix composed of characteristic 2×2submatrices. They are also used to build the matrices of more complex PT -symmetric bosonic systems. The characteristic block structure — typical of bosonic systems with quadratic Hamiltonians — forms a critical ingredient of the analytical framework presented below, allowing us to uncover the complex structures of eigenvalues and their associated eigenvectors of the corresponding dynamical matrices in the space of system’s parameters. Moreover, we also consider bosonic systems in which damping and amplification are not in balance, which is a typical situation of PT-symmetric systems. In this unbalanced case, relying on the results presented in Ref. [55], we may introduce a specific interaction picture in which the average damping or amplification dynamics is projected out and the remaining dynamics exhibits the features found in PT -symmetric systems. Under these conditions, we may write the quadratic Hamiltonian of the considered two-mode bosonic system in the interaction picture as follows: ˆ H2=h¯hϵˆa† 1ˆa2+ ¯hκˆa1ˆa2i+ H.c.,(1) where ˆaj(ˆa† j) stands for the annihilation (creation) operator of the jth mode, j= 1,2,ϵis the linear coupling strength between the modes, and κis the nonlinear coupling strength between the modes. Symbol H.c. replaces the Hermitian-conjugated terms. Whereas the linear coupling originates in the spatial overlap of the mode electric-field amplitudes, the nonlinear coupling arises in the three-mode parametric process with strong pumping [56]. The damping (amplification) of modes occurs as a consequence of the interaction with the reservoir whose two-level atoms are in the ground (excited) state [43]. Projecting out the reservoir two-level atoms, we are left with the damping (amplification) rate γjof jth damped (amplified) mode and the corresponding Langevin stochastic operator forces, ˆ Ljand ˆ L† j, in the dynamical Heisenberg-Langevin equations. We note that properties of the Langevin operator forces differ for the damping and amplification processes and they have to be chosen such that the field-operator commutation relations are fulfilled. This results in the corresponding fluctuation-dissipation theorems [47,57] formulated within the Heisenberg-Langevin formalism. Using the two-mode Hamiltonian in Eq. (1), we derive the Heisenberg-Langevin equations written for the vector ˆ a= [ˆ a1,ˆ a2]T≡hˆa1,ˆa† 1,ˆa2,ˆa† 2iTof field operators and vector ˆ L=hˆ L1,ˆ L† 1,ˆ L2,ˆ L† 2iTof the Langevin operator forces as follows [58]: dˆ a dt =−iM(2) ˆ a+ˆ L.(2) Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 4 ˆa1ˆa2 , κ γ1γ2 (a) ˆa1ˆa2ˆa3 , κ , κ γ1γ1+γ3 2γ3 (b) ˆa2 ˆa1ˆa3ˆa4 , κ , κ, κ γ12 γ12 γ34 γ34 (c) ˆa1 ˆa2 ˆa3 ˆa4 , κ , κ , κ , κ γ12 γ12 γ34 γ34 (d) ˆa1ˆa2ˆa3ˆa4 , κ , κ, κ γ13 γ24 γ13 γ24 (e) ˆa2 ˆa1ˆa3ˆa4ˆa5 , κ , κ , κ , κ γ12 γ12 γ12+γ45 2γ45 γ45 (f) ˆa1 ˆa2 ˆa3 ˆa4 ˆa5, κ , κ , κ , κ , κ , κ , κ , κ γ12 γ12 γ34 γ34 γ12+γ34 2 (g) Figure 1: Schematic diagrams of the bosonic systems composed of (a) two, (b) three, (c—e) four, and (f,g) five modes with typical linear, circular, and pyramid configurations that exhibit quantum exceptional points (QEPs) and quantum hybrid points (QHPs) with various exceptional degeneracies (EDs) and diabolical degeneracies (DDs) observed in the dynamics of field-operator moments (FOMs) of different orders. Strengths ϵand κcharacterize, respectively, the linear and nonlinear coupling between the modes, γ, with subscripts indicating the mode number(s), are the damping or amplification rates, and annihilation operators ˆaidentify the mode number via their subscripts, and γjk indicates that γj=γk. In Eq. (2), the dynamical matrix M(2), M(2) ="−i˜γ1ξ ξ−i˜γ2 ,#(3) is expressed in terms of 2×2submatrices ˜γj,j= 1,2, and ξis defined as: ˜γj="γj/2 0 0γj/2#,ξ="ϵ κ −κ−ϵ#,(4) where γjis the damping or amplification rate of the mode jthat is accompanied by the corresponding Langevin operator forces. The 2×2matrices, ˜γj(j= 1,2) and ξ, given in Eq. (4) can be simultaneously diagonalized using the diagonalization transformation appropriate to the matrix ξas the remaining two matrices are linearly proportional to the unity matrix and, thus, are not modified by the transformation. This diagonalization transformation then decomposes the 4×4matrix M(2) in Eq. (3) into the direct sum of two independent 2×2matrices belonging to the eigenvalues λξ 1and λξ 2of the matrix ξ[see Eq. (7) below] and having the Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 5 structure of the original 2×2matrix M(2): M(2) =      −iγ1/20λξ 10 0−iγ1/20λξ 2 λξ 10−iγ2/20 0λξ 20−iγ2/2       . This step is critical to our analytical approach, and it is applied analogously to systems with a larger number of bosonic modes. It enables an effective reduction in the dimensionality of the diagonalization problem, thereby allowing the derivation of comprehensive analytical expressions essential for analyzing QEPs and QDPs. Denoting a general eigenvalue of the 2×2matrix ξas ξ, we may express the eigenfrequencies λM(2) 1,2and eigenvectors yM(2) 1,2of the decomposing 2×2matrices of the 4×4matrix M(2) in Eq. (2) in the following common form: λM(2) 1,2=−iγ+∓β(5) and yM(2) 1,2=−iγ−±β ξ,1T ,(6) where 4γ±=γ1±γ2and β2=ξ2−γ2 −. The eigenvalues λξ 1,2and the corresponding eigenvectors yξ 1,2of the matrix ξare simply derived in the form: λξ 1,2=∓ζ(7) and yξ 1,2=−ϵ∓ζ κ,1T ,(8) where ζ=√ϵ2−κ2. Combing the above two results for the matrix diagonalization, we can express the eigenvalues ΛM(2) of the 4×4matrix M(2) in Eq. (2)asΛM(2) 1,3=λM(2) 1,2for ξ=λξ 1and ΛM(2) 2,4=λM(2) 1,2for ξ=λξ 2, i.e.: ΛM(2) 1= ΛM(2) 2=−iγ+−β, ΛM(2) 3= ΛM(2) 4=−iγ++β, (9) and β=qζ2−γ2 −. We note that, in general, we use capital Greek letters Λto denote the eigenvalues of the dynamical matrices Min their full dimensions, and lowercase Greek letters λfor the eigenvalues in the reduced (half) dimensions. We also note that the average damping or amplification rate γ+= 0 in the usual PT -symmetric systems, that are, however, only specific cases in our general calculations. Similarly as the eigenvalues, we obtain the eigenvectors along the formulas YM(2) j=hyM(2) 1,1(ξ=λξ j)yξ j, yM(2) 1,2(ξ=λξ j)yξ jiT, for j= 1,2, YM(2) j=hyM(2) 2,1(ξ=λξ j−2)yξ j−2, yM(2) 2,2(ξ=λξ j−2)yξ j−2iT, for j= 3,4,(10) Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 6 ΛM(n) 2j−1=λM(n) j(ξ=λξ 1) ΛM(n) 2j=λM(n) j(ξ=λξ 2) YM(n) 2j−1=              yM(n) j,1(ξ=λξ 1)"yξ 1,1 yξ 1,2# yM(n) j,2(ξ=λξ 1)"yξ 1,1 yξ 1,2# ... yM(n) j,n (ξ=λξ 1)"yξ 1,1 yξ 1,2#              YM(n) 2j=              yM(n) j,1(ξ=λξ 2)"yξ 2,1 yξ 2,2# yM(n) j,2(ξ=λξ 2)"yξ 2,1 yξ 2,2# ... yM(n) j,n (ξ=λξ 2)"yξ 2,1 yξ 2,2#              j= 1,...,n. Figure 2: Schematic diagram of the structure of eigenvalues ΛM(n) jand the corresponding eigenvectors YM(n) j,j= 1,...,2n, of a general dynamical 2n×2nmatrix M(n)built from the eigenvalues λM(n) k and the corresponding eigenvectors yM(n) k,k= 1, . . . , n, of the matrix M(n)considered as an n×n matrix composed of 2×2submatrices at the positions of its elements; n= 2,3, . . .. The eigenvalues λξ 1,2and the corresponding eigenvectors yξ 1,2belong to the submatrix ξ. in the form: YM(2) 1,2=h(∓ϵ+ζ)χ κζ ,±χ ζ,−(ϵ∓ζ) κ,1iT, YM(2) 3,4=h(±ϵ−ζ)χ∗ κζ ,∓χ∗ ζ,−(ϵ∓ζ) κ,1iT,(11) where χ=iγ−+β. The structure of eigenvalues and eigenvectors of a general dynamical 2n×2nmatrix M(n)formed by the eigenvalues λM(n) kand eigenvectors yM(n) kof the n×n matrix M(n)expressed using 2×2submatrices and the eigenvalues λξ 1,2and eigenvectors of the submatrix ξis shown in Fig. 2. The formula in Eq. (9) provides the two pairs of coinciding eigenvalues. However, Eq. (11) for the corresponding eigenvectors reveal no eigenvector degeneracy in general. On the other hand, χis purely imaginary for β= 0, which leads to YM(2) 1=YM(2) 3and YM(2) 2=YM(2) 4, while having all the eigenvalues the same. So we observe the secondorder degeneration in the Hilbert space that is formed by two second-order QEPs with identical eigenvalues. We, thus, have a QHP with second-order DD and ED in this case. The condition β= 0 implies that the eigenvalues and eigenvectors of the 2×2matrix M(2) given in Eqs. (5) and (6) coincide. The second-order ED, thus, originates in the 2×2form of matrix M(2). We note that this degeneracy can be verified by transforming the matrix M(2) for β= 0 into its Jordan form JM(2) ="−iγ+1 0−iγ+#.(12) On the other hand, the eigenvalues of 2×2submatrix ξwritten in Eq. (7) point out at the origin of DD: They differ just by the sign, but they lead to the same value of β, i.e., to the same eigenvalue λM(2) 1,2of the 2×2matrix M(2). DD is then implied by the fact that the eigenvectors yξ 1,2of 2×2submatrix ξin Eq. (8) differ for ζ= 0. These findings about the structure of the matrices describing the dynamical equations have their counterpart in the structure of the analyzed two-mode PT -symmetric system. The 2×2matrix ξ, defined in Eq. (4), connects pairs of the annihilation and creation operators of different modes. As such it forms the basic building block of more complex PT -symmetric systems, together with the damping and/or amplification matrices ˜γjin Eq. (4). In more complex PT-symmetric systems, the dependencies of the eigenvalues of the dynamical n×nmatrices M(n)(n= 2,3, . . .) on the eigenvalues λξ 1,2of the ξ Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 7 matrix are typically quadratic. This, thus, results in the observation of second-order DD in the dynamical features of the matrices M(n). In this case, the eigenvalue analysis of the considered systems considerably simplifies and we may restrict our attention to only the n×nmatrices M(n)with their elements in the form of 2×2submatrices, when a detailed eigenvalue analysis is performed. The consequences of the eigenvalue analysis are then combined with the above second-order DD. The condition β= 0 for observing QHPs can be analyzed in the space of system parameters (ϵ, κ, γ1, γ2)as follows. The quadratic Hamiltonian ˆ H(2) in Eq. (1) provides the linear Heisenberg-Langevin equations that allow for the temporal rescaling ϵt, i.e., only the relative parameters (κ/ϵ, γ1/ϵ, γ2/ϵ)suffice in characterizing the system dynamics. Moreover, the structure of the analyzed systems is such that the average damping or amplification rate γ+influences equally only the eigenvalues, but it does not modify the eigenvectors. This makes the space of the system parameters effectively two-dimensional with the spanning parameters (κ/ϵ, γ−/ϵ). Using these parameters, the condition β= 0 is expressed as follows: κ2 ϵ2+γ2 − ϵ2= 1.(13) Thus, the QHPs form a circle with the unit radius in the space (κ/ϵ, γ−/ϵ), as shown by the red dashed curve in Fig. 3(a). We note that, at this circle, there is a specific point at γ−= 0 in which all four eigenvectors, given in Eq. (10), are the same which gives raise to the fourth-order QEP. However, this corresponds to the system in which both modes are equally damped or amplified. It is worth noting that the eigenvectors of both matrices M(2), given in Eq. (3), and ξ, in Eq. (4), are degenerated at this point. Deeper insight into the structure of higher-order FOMs, as well as the system dynamics, can be obtained once we transform the Heisenberg-Langevin equations (2) into the form in which the dynamical matrix has the diagonal form. Using the transformation matrix P=hYM(2) 1,YM(2) 2,YM(2) 3,YM(2) 4iformed from the eigenvectors in Eq. (10), we define the corresponding field operators ˆ b=hˆ b1,ˆ b2,ˆ b† 1,ˆ b† 2iT, the Langevin operator forces ˆ K= hˆ K1,ˆ K2,ˆ K† 1,ˆ K† 2iT, and the diagonal dynamical matrix Λ(2): Λ(2) =P−1M(2)P,ˆ b=P−1ˆ a,ˆ K=P−1ˆ L.(14) We note that the order of elements in the operator vector ˆ b(and similarly in the vector ˆ K of the accompanying Langevin operator forces) is given by the numbering of the eigenvalues in Eq. (9) and the corresponding diagonalization transform. In the transformed basis, the Heisenberg-Langevin equations take the form: dˆ b dt =−iΛ(2)ˆ b+ˆ K.(15) We note that the positions of the newly-defined annihilation and creation operators in the vector ˆ b, as well as the positions of the accompanying Langevin operator forces in the vector ˆ K, differ from those in the original vectors ˆ aand ˆ Ldefined above Eq. (2). The solution to Eq. (15) is expressed as: ˆ b(t) = exp(−iΛ(2)t)ˆ b(0) +Zt 0 dt′exp[−iΛ(2)(t−t′)] ˆ K(t′).(16) Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 8 (a) (b) (c) (d) (e) (f) (g) (h) Figure 3: Real parts λrof the eigenvalues (a) λM(2) 1,2of the matrix M(2), given in Eq. (5), for the two-mode bosonic system, (b) λM(3) 1,2,3of the matrix M(3), given in Eq. (23), for the three-mode linear bosonic system, (c) [(d)] λM(4) l1 1,...,4[λM(4) l2 1,...,4] of the matrix M(4) l1 [M(4) l2 ], given in Eq. (30) [(34)] for the four-mode linear bosonic system with neighbor modes having equal [different] damping and/or amplification rates, (e) λM(4) c1 1,...,4of the matrix M(4) c1 , given in Eq. (40), for the four-mode circular bosonic system with neighbor modes having equal damping and/or amplification rates, (f) λM(5) l 2,...,5of the matrix M(5) l, given by Eq. (47), for the five-mode linear bosonic system, and (g,h) λM(5) p 2,...,5of the matrix M(5) p, given in Eq. (54), for the five-mode pyramid bosonic system assuming (g) β1= 0 and (h) β2= 0. The eigenvalues are drawn in the parameter space (κ/ϵ, γ−/ϵ), where ϵ(κ) is the linear (nonlinear) coupling strength and γ−the difference of the damping/amplification rates in individual models. Dashed red curves indicate the positions of the QHPs given by (a,e,g) Eq. (13), (b) Eq. (25), (c) Eq. (32), (d) Eq. (36), (f) Eq. (49), and (h) Eq. (56). Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 9 Defining the vectors ˆ a= [ˆ a1,ˆ a2,ˆ a3,ˆ a4,ˆ a5]T≡hˆa1,ˆa† 1,ˆa2,ˆa† 2,ˆa3,ˆa† 3,ˆa4,ˆa† 4,ˆa5,ˆa† 5iTof field operators and ˆ L=hˆ L1,ˆ L† 1,ˆ L2,ˆ L† 2,ˆ L3,ˆ L† 3,ˆ L4,ˆ L† 4,ˆ L5,ˆ L† 5iTof the Langevin operator forces, we obtain the following Heisenberg-Langevin equations: dˆ a dt =−iM(5) lˆ a+ˆ L.(43) The dynamical 10 ×10 matrix M(5) lintroduced in Eq. (43) is written with the help of 2×2submatrices introduced in Eq. (4) as M(5) l=       −i˜γ1ξ000 ξ−i˜γ2ξ0 0 0ξ−i˜γ3ξ0 0 0 ξ−i˜γ4ξ 0 0 0 ξ−i˜γ5        .(44) Motivated by PT symmetry we assume: γ1=γ2≡γ12, γ4=γ5≡γ45.(45) Moreover, inspired by the condition in Eq. (22) derived for the three-mode linear system, we additionally assume: 2γ3=γ12 +γ45.(46) Then, diagonalization of the 5×5matrix M(5) lresults in the following eigenvalues λM(5) l j: λM(5) l 1=−iγ+, λM(5) l 2,3=−iγ+±α−, λM(5) l 4,5=−iγ+±α+.(47) The corresponding eigenvectors are derived as follows: yM(5) l 1="1,iγ− ξ,−γ2 − ξ2−1,−iγ− ξ,1#T , yM(5) l 2,4="−1∓n∓ 4ξ3,−m∓ ξ3∓2βχ∓ ξ2,χ2 ∓ ξ2−1,χ∓ ξ,1#T , yM(5) l 3,5="−1∓n∗ ∓ 4ξ3,m∗ ∓ ξ3±2βχ∗ ∓ ξ2,χ∗2 ∓ ξ2−1,−χ∗ ∓ ξ,1#T , (48) where n∓= (2β∓ξ)[(2β∓ξ)2−4iγ−α∓],m∓=iγ−(χ2 ∓+ξ2),χ±=−iγ++α±, α2 ±=β2+ 7ξ2/4±2βξ,β2=ξ2/4−γ2 −, and 4γ±=γ12 ±γ45. If β= 0 then α+=α−and χ+=χ−. Under these conditions, we have λM(5) l 2=λM(5) l 4 and yM(5) l 2=yM(5) l 4. We also have λM(5) l 3=λM(5) l 5and yM(5) l 3=yM(5) l 5. Thus, we observe two QEPs with second-order EDs that give rise to two QHPs with second-order EDs and Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 16 DDs when the 10×10 matrix M(5) lis analyzed. In the parameter space (κ/ϵ, γ−/ϵ), these QHPs are localized at the positions obeying the condition κ2 ϵ2+4γ2 − ϵ2= 1.(49) The real parts of the eigenvalues λM(5) l jforming QHPs (j= 2,...,5) are plotted in Fig. 3(f) where their coinciding values for the condition in Eq. (49) are indicated by the red dashed curves. The eigenvalues and eigenvectors of the 10 ×10 matrix M(5) lconstructed from the formulas given in Eqs. (47) and (48) together with those in Eqs. (7) and (8) (for details, see Appendix A) allow us to describe the system dynamics via a diagonal dynamical matrix. The appropriate operators ˆ b=hˆ b1,ˆ b† 1,ˆ b2,ˆ b3,ˆ b† 2,ˆ b† 3,ˆ b4,ˆ b5,ˆ b† 4,ˆ b† 5iTthen reveal QEPs and QHPs found in the dynamics of FOMs of different orders. For the firstand second-order FOMs they are written in Tab. 5of Appendix B. 5.2 Pyramid configuration Second-order QEPs and QHPs can also be identified in the pyramid configuration [the only considered non-planar configuration, see Fig. 1(g)] described by the following Hamiltonian ˆ H5,p: ˆ H5,p=h¯hϵˆa† 1ˆa2+ ¯hϵˆa† 1ˆa4+ ¯hϵˆa† 1ˆa5+ ¯hϵˆa† 2ˆa3+ ¯hϵˆa† 2ˆa5+ ¯hϵˆa† 3ˆa4+ ¯hϵˆa† 3ˆa5 +¯hϵˆa† 4ˆa5+ ¯hκˆa1ˆa2+ ¯hκˆa1ˆa4+ ¯hκˆa1ˆa5+ ¯hκˆa2ˆa3+ ¯hκˆa2ˆa5 +¯hκˆa3ˆa4+ ¯hκˆa3ˆa5+ ¯hκˆa4ˆa5] + H.c. (50) The Hamiltonian ˆ H5,pgiven in Eq. (50) leads to the Heisenberg-Langevin equations (43) in which the dynamical 10 ×10 matrix M(5) pattains the form using the 2×2submatrices introduced in Eq. (4): M(5) p=       −i˜γ1ξ0ξ ξ ξ−i˜γ2ξ0ξ 0ξ−i˜γ3ξ ξ ξ0ξ−i˜γ4ξ ξ ξ ξ ξ −i˜γ5        .(51) Motivated by the PT symmetry we assume: γ1=γ2≡γ12, γ3=γ4≡γ34.(52) Moreover, inspired by the condition (22) derived for the three-mode linear bosonic system, we additionally assume: 2γ5=γ12 +γ34.(53) Under these conditions, diagonalization of the 5×5matrix M(5) pleaves us the following five eigenvalues λM(5) p j: λM(5)p 1=−iγ+, λM(5) p 2,3=−iγ+−ξ±β1, λM(5) p 4,5=−iγ++ξ±β2.(54) Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 17 The corresponding eigenvectors are obtained as follows: yM(5) p 1=−iξ γ− ,−iξ γ− ,iξ γ− ,iξ γ− ,1T , yM(5) p 2=χ1 ξ,−χ1 ξ,−1,1,0T , yM(5) p 3=−χ∗ 1 ξ,χ∗ 1 ξ,−1,1,0T , yM(5) p 4=σ+, σ+, σ∗ +, σ∗ +,1T, yM(5) p 5=σ∗ −, σ∗ −, σ−, σ−,1T, (55) where σ±= (ξ±χ2)/(4ξ),β2 1=ξ2−γ2 −,β2 2= 5ξ2−γ2 −,χ1,2=β1,2−iγ−, and 4γ±= γ12 ±γ34. Contrary to the eigenvalues of the models discussed above, the eigenvalues λM(5) p jfor j= 2,...,5in Eq. (55) exhibit the linear dependence on ξ. This means that the DD inherited to all the above-discussed models is removed in this model. It remains only for the eigenvalue λM(5) p 1that, however, has no ability to form EDs. Provided that β1= 0, we have λM(5) p 2=λM(5) p 3and yM(5) p 2=yM(5) p 3. Thus, we have the QEP with second-order ED. When the 10 ×10 matrix M(5) pis analyzed, there occur one second-order QEP for ξ=ζand one second-order QEP for ξ=−ζ. These QEPs occur in the parameter space (κ/ϵ, γ−/ϵ)under the condition written in Eq. (13). Moreover, in parallel, if β2= 0, it holds that λM(5) p 4=λM(5) p 5and yM(5) p 4=yM(5) p 5. Thus, similarly as above, we observe the QEP with second-order ED for the 5×5matrix M(5) p. There occur one second-order QEP for ξ=ζand one second-order QEP for ξ=−ζfor the 10 ×10 matrix M(5) p. These QEPs are localized in the parameter space (κ/ϵ, γ−/ϵ) at the points fulfilling the condition: κ2 ϵ2+γ2 − 5ϵ2= 1.(56) The real parts of the eigenvalues λM(5) p jforming QEPs (j= 2,...,5) are drawn in Fig. 3(g,h); the coinciding values occurring for the condition in Eq. (13) [(56)] are plotted in Fig. 3(g) [(h)] as red dashed curves. We note that for ξ= 0 (i.e., γ−= 0,κ=ϵ) we have a specific QHP with ten-fold frequency degeneracy that belongs to four doubly-degenerate eigenvectors and another two eigenvectors. The eigenvalues and eigenvectors of the 10 ×10 matrix M(5) pformed from the expressions given in Eqs. (54) and (55) and also in Eqs. (7) and (8) (for details, see Appendix A) using the scheme in Fig. 2allow us to analyze the system dynamics via a diagonal dynamical matrix. The appropriate operators ˆ b=hˆ b1,ˆ b† 1,ˆ b2,ˆ b3,ˆ b† 3,ˆ b† 2,ˆ b4,ˆ b5,ˆ b† 5,ˆ b† 4iTthen allow to identify the QEPs and QHPs observed in the dynamics of FOMs of different orders. For the firstand second-order FOMs, they are explicitly given in Tab. 6in Appendix B. Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 18 6 Conclusions We have analyzed the dynamics of simple bosonic systems described by quadratic nonHermitian Hamiltonians from the point of view of the occurrence of quantum exceptional, diabolical, and hybrid points. Non-Hermiticity of the considered systems, composed of from two to five coupled modes, originated in their damping and/or amplification, that are accompanied by the corresponding Langevin fluctuating forces to assure the physically consistent behavior. We have identified specific configurations defined by two-mode couplings and conditions for the damping and amplification rates of the modes at which the inherited quantum exceptional and diabolical points occur. Surprisingly, in these physically consistent models, we have found only secondand third-order inherited quantum exceptional points including their doubling due to second-order diabolical degeneracies. On the other hand, doubled second-order inherited quantum exceptional points have been observed. We have shown that the analyzed bosonic systems naturally exhibit the secondorder diabolical degeneracies. Nevertheless, we have found an exception from this behavior in which no diabolical degeneracy occurs. The exceptional and diabolical degeneracies of inherited quantum hybrid points have then been used to construct higher-order degeneracies observed in the dynamics of secondorder field-operator moments. The corresponding quantum exceptional and hybrid points are summarized in tables that demonstrate richness of the evolution of the general-order field-operator moments and serve for future studies of the related physical effects. The investigations have revealed the need for further looking for the bosonic systems exhibiting higher-order inherited quantum exceptional and hybrid points by considering more general bosonic systems. In Ref. [50] we extend our investigations to the systems with partial PT -symmetry like dynamics (nonconventional PT-symmetry) as well as nonHermitian bosonic systems with unidirectional coupling. We may conclude in general that the performed analysis opens the door for further detailed investigations of the role of exceptional and diabolical degeneracies responsible for inducing unusual physical effects observed in physically well-behaved systems at exceptional, diabolical, and hybrid points. 7 Acknowledgements The authors thank Ievgen I. Arkhipov for useful discussions. J.P. and K.T. acknowledge support by the project OP JAC CZ.02.01.01/00/22 008/0004596 of the Ministry of Education, Youth, and Sports of the Czech Republic. J.P. acknowledges support by the project No. 25-15775S of the Czech Science Foundation. A.K.-K., G.Ch., and A.M. were supported by the Polish National Science Centre (NCN) under the Maestro Grant No. DEC-2019/34/A/ST2/00081. A Eigenvalues and eigenvectors of three-, fourand five-mode bosonic systems In this Appendix, we present the eigenvalues and eigenvectors of the dynamical matrices of the Heisenberg-Langevin equations for the three-, four-, and five-mode bosonic systems in configurations depicted in Fig. 1. Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 19 A.1 General n-mode bosonic systems The eigenvalues of the 2n×2ndynamical matrix M(n)belonging to a system with n modes (n= 2,3, . . .) are constructed from the eigenvalues λM(n) j,j= 1,...,n, derived for the form of the dynamical matrix containing 2×2submatrices and the eigenvalues λξ j, j= 1,2, of the matrix ξgiven in Eq. (7) as follows: ΛM(n) 2j−1=λM(n) j(ξ=λξ 1), ΛM(n) 2j=λM(n) j(ξ=λξ 2), j = 1,2, . . . , n. (57) Similarly, using the corresponding eigenvectors yM(n) jof the n×nmatrix M(n)formed by submatrices and the eigenvectors yξ j,j= 1,2, of the matrix ξin Eq. (8), we arrive at the following eigenvectors of the 2n×2nmatrix M(n)associated with the eigenvalues given in Eq. (57): YM(n) 2j−1=      yM(n) j,1(ξ=λξ 1)yξ 1 yM(n) j,2(ξ=λξ 1)yξ 1 . . . yM(n) j,n (ξ=λξ 1)yξ 1       , YM(n) 2j=      yM(n) j,1(ξ=λξ 2)yξ 2 yM(n) j,2(ξ=λξ 2)yξ 2 . . . yM(n) j,n (ξ=λξ 2)yξ 2       , j = 1, . . . , n. (58) A.2 Three-mode bosonic system In the three-mode system (n= 3) assuming β=q2ζ2−γ2 −= 0, we find a QHP with third-order ED and second-order DD. All the eigenvalues in Eq. (57) are the same in this case and the eigenvectors in Eq. (58) obey the relations YM(3) 1=YM(3) 3=YM(3) 5and YM(3) 2=YM(3) 4=YM(3) 6. A.3 Four-mode bosonic systems In the four-mode system (n= 4) in the linear configuration with equal damping and/or amplification rates of neighbor modes and assuming µ= 2qζ2(5ζ2/16 −γ2 −) = 0, we observe two QHPs with second-order ED and second-order DD. The eigenvalues ΛM(4) l1 1,2,5,6 and also the eigenvalues ΛM(4) l1 3,4,7,8in Eq. (57) coincide. The eigenvectors in Eq. (58) obey the relations YM(4) l1 j=YM(4) l1 j+4 for j= 1,...,4. In the four-mode system (n= 4) in the linear configuration with different damping and/or amplification rates in neighbor modes and assuming α+= 0 [α−= 0], α±= q(3 ±√5)ζ2/2−γ2 −)=0, we find a single QHP with second-order ED and DD. The eigenvalues ΛM(4) l2 1,2,3,4[ΛM(4) l2 5,6,7,8] in Eq. (57) are equal. The eigenvectors in Eq. (58) fulfil: YM(4) l2 1=YM(4) l3 3and YM(4) l2 2=YM(4) l2 4[YM(4) l2 5=YM(4) l2 7and YM(4) l2 6=YM(4) l2 8]. In the four-mode system (n= 4) in the circular configuration with equal damping and/or amplification rates of neighbor modes and assuming β=qζ2−γ2 −= 0, we have Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 20 two QHPs with second-order EDs and second-order DDs. The eigenvalues and eigenvectors in this case behave analogous to those given above for the four-mode system in the linear configuration with equal damping and/or amplification rates of neighbor modes. A.4 Five-mode bosonic systems In the five-mode system (n= 5) in the linear configuration and assuming β=qζ2/4−γ2 −= 0, we have two QHPs with second-order EDs and second-order DDs. The eigenvalues ΛM(5) l 3,4,7,8and also the eigenvalues ΛM(5) l 5,6,9,10 in Eq. (57) are the same. The eigenvectors in Eq. (58) fulfill the relations: YM(5) l 3=YM(5) l 7,YM(5) l 4=YM(5) l 8,YM(5) l 5=YM(5) l 9, and YM(5) l 6=YM(5) l 10 . In the five-mode system (n= 5) in the pyramid configuration and assuming β1= qζ2−γ2 −= 0, we have two QEPs with second-order ED. The eigenvalues ΛM(5) p 3and ΛM(5) p 5together with their accompanying eigenvectors YM(5) p 3=YM(5) p 5coincide. The same is true for eigenvalues ΛM(5) p 4and ΛM(5) p 6and eigenvectors YM(5) p 4and YM(5) p 6. Provided that β2=q5ζ2−γ2 −= 0, we reveal two QEPs with second-order ED. The eigenvalues ΛM(5) p 7 and ΛM(5) p 9together with their accompanying eigenvectors YM(5) p 7and YM(5) p 9are the same. Similarly, the eigenvalues ΛM(5) p 8and ΛM(5) p 10 and eigenvectors YM(5) p 8and YM(5) p 10 equal. B QEPs and QHPs in firstand second-order field-operator-moments spaces: Comparison Here, we present tables that summarize the QEPs and QHPs, along with their degeneracies, as observed in the firstand second-order FOM spaces for the bosonic systems analyzed in the main text. This analysis follows the methodology developed in Ref. [43], which derives the genuine and induced QEPs and QHPs — together with their degeneracies — from the inherited QEPs and QHPs identified in the dynamics of first-order FOMs. The genuine QEPs can be viewed as analogues of the inherited QEPs in the dynamical matrix of first-order FOMs, now extended to the dynamical matrices of higher-order FOMs. However, a crucial distinction arises due to the structure of higher-order FOMs, which are constructed as products of a fixed number of field operators, incorporating all combinations of annihilation and creation operators. These combinations include terms that are mutually related by commutation relations. Since such terms contribute to the system dynamics only once, only one representative term should be considered when analyzing spectral degeneracies — this yields the genuine QEPs. In contrast, if all such terms are retained in the analysis, without accounting for redundancy due to commutation relations, the resulting degeneracies are referred to as induced QEPs. The same distinction applies to QHPs, leading to the identification of both genuine and induced QHPs [43]. Different configurations of QEPs and QHPs emerge across different models, highlighting the rich variety of spectral degeneracy structures and their implications for diverse physical phenomena. For the two-mode bosonic system with its dynamical matrix M(2) given in Eq. (3), we summarize the corresponding degeneracies up to second-order FOMs in Tab. 1. This provides a direct comparison of the QEP and QHP degeneracies observed in the dynamics Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 21 Λi jΛr jMoments Moment Genuine and induced QHPs Genuine QHPs deg. Partial Partial Partial Partial QDP x QDP x QDP x QDP x QEP deg. QEP deg. QEP deg. QEP deg. γ+∓β⟨ˆ b1⟩,⟨ˆ b† 1⟩ ⟨ ˆ B1⟩1 1x2 2x2 1x2 2x2 ⟨ˆ b2⟩,⟨ˆ b† 2⟩ ⟨ ˆ B2⟩1 1x2 1x2 2γ+∓2β⟨ˆ b1ˆ b2⟩,⟨ˆ b† 1ˆ b† 2⟩ ⟨ ˆ B1ˆ B2⟩2 2x4 4x4 1x4 1x4 β−β⟨ˆ b† 1ˆ b2⟩2 + β−β⟨ˆ b1ˆ b† 2⟩2 ∓2β⟨ˆ b2 1⟩,⟨ˆ b†2 1⟩ ⟨ ˆ B2 1⟩1 1x4 1x3 2x3 β−β⟨ˆ b† 1ˆ b1⟩2 ∓2β⟨ˆ b2 2⟩,⟨ˆ b†2 2⟩ ⟨ ˆ B2 2⟩1 1x4 1x3 β−β⟨ˆ b† 2ˆ b2⟩2 Table 1: Real and imaginary parts of the complex eigenfrequencies Λr j−iΛi jof the matrix M(2) given in Eq. (3) for the two-mode bosonic system derived from the equations for the FOMs up to second order. The corresponding moments in the ‘diagonalized’ field operators are written together with their degeneracies (deg.) coming from different positions of the field operators. The DDs of QHPs (partial DDs) derived from the indicated FOMs and the EDs of the constituting QEPs are given. Both genuine and induced QEPs and QHPs are considered. The operator vectors ˆ Bjfor j= 1,2are defined in the rows written for Λi j=γ+devoted to the first-order FOMs, i.e. ˆ Bj≡[ˆ bj,ˆ b† j]for j= 1,2. Symbol ˆ Bjˆ Bk,j, k = 1,2, stands for the tensorial product that provides four terms explicitly written in the rows for Λi j= 2γ+; the terms derived from those explicitly written by using the commutation relations are omitted. of firstand second-order FOMs. Notably, we identify genuine QEPs with a fourth-order ED in the second-order FOM dynamics, in contrast to the inherited QEPs exhibiting only second-order EDs in the first-order FOM dynamics. It is instructive to compare the QEP and QHP degeneracies listed in Tab. 1with those observed in systems featuring larger numbers of modes, as presented below. This comparison highlights how the spectral degeneracy structures in first-order FOMs evolve when extended to second-order FOMs. Notably, characteristic features emerge that can be generalized to spectral degeneracies of arbitrary-order FOMs, as discussed in Ref. [43]. For instance, one can infer the maximal EDs and maximal DDs occurring in the dynamics of FOMs of a given order. Considering the linear three-mode bosonic system with its dynamical matrix M(3) given in Eq. (21) under the condition given in Eq. (22), the QEPs and QHPs predicted in the dynamics of firstand second-order FOMs are summarized in Tab. 2. As shown there, the inherited QEPs with a third-order ED give rise to genuine QEPs with a ninth-order ED when analyzing the dynamics of second-order FOMs. In the linear four-mode bosonic system described by the dynamical matrix M(4) lgiven in Eq. (28), which features equal damping and/or amplification rates for neighboring modes [see Eq. (29)], the QEPs and QHPs in the dynamics of firstand second-order FOMs are summarized in Tab. 3. These results also apply to the circular four-mode bosonic system with dynamical matrix M(4) cshown in Eq. (38), assuming equal damping and/or amplification rates of neighboring modes [see Eq. (39]. According to Tab. 3, the inherited Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 22 Λi jΛr jMoments Moment Genuine and induced QHPs Genuine QHPs deg. Partial Partial Partial Partial QDP x QDP x QDP x QDP x QEP deg. QEP deg. QEP deg. QEP deg. γ+0, ∓β⟨ˆ b1⟩,⟨ˆ b2⟩,⟨ˆ b† 2⟩ ⟨ ˆ B1⟩1 1x3 2x3 1x3 2x3 0, ∓β⟨ˆ b† 1⟩,⟨ˆ b3⟩,⟨ˆ b† 3⟩ ⟨ ˆ B2⟩1 1x3 1x3 2γ+0⟨ˆ b2 1⟩ ⟨ ˆ B2 1⟩1 1x9 4x9 1x6 2x6 ∓β⟨ˆ b1ˆ b2⟩,⟨ˆ b1ˆ b† 2⟩2 + ∓2β⟨ˆ b2 2⟩,⟨ˆ b†2 2⟩1 β−β⟨ˆ b† 2ˆ b2⟩2 0⟨ˆ b† 1ˆ b1⟩ ⟨ ˆ B2ˆ B1⟩2 2x9 1x9 1x9 ∓β⟨ˆ b† 1ˆ b2⟩,⟨ˆ b† 1ˆ b† 2⟩2 ∓β⟨ˆ b3ˆ b1⟩,⟨ˆ b† 3ˆ b1⟩2 ∓2β⟨ˆ b3ˆ b2⟩,⟨ˆ b† 3ˆ b† 2⟩2 β−β⟨ˆ b3ˆ b† 2⟩,⟨ˆ b† 3ˆ b2⟩2 0⟨ˆ b†2 1⟩ ⟨ ˆ B2 2⟩1 1x9 1x6 ∓β⟨ˆ b† 1ˆ b3⟩,⟨ˆ b† 1ˆ b† 3⟩2 ∓2β⟨ˆ b2 3⟩,⟨ˆ b†3 3⟩1 β−β⟨ˆ b† 3ˆ b3⟩2 Table 2: Real and imaginary parts of the complex eigenfrequencies Λr j−iΛi jof the matrix M(3), given in Eq. (21), for the linear three-mode bosonic system derived from the equations for the FOMs up to second order. We have ˆ B1≡[ˆ b1,ˆ b2,ˆ b† 2]and ˆ B2≡[ˆ b† 1,ˆ b3,ˆ b† 3]and more details are given in the caption to Tab. 1. QEPs with doubled second-order ED give rise to genuine QEPs exhibiting fourth-order ED and sixth-order DD when analyzing second-order FOM dynamics. On the other hand, for the linear four-mode bosonic system described by the dynamical matrix M(4) lgiven in Eq. (28), with varying damping and/or amplification rates of neighboring modes [see Eq. (33)], the QEPs and QHPs observed in the dynamics of firstand second-order FOMs are summarized in Tab. 4. The spectral structure of first-order FOMs, featuring two inherited QEPs with second-order EDs, evolves into a more complex second-order FOM structure that includes genuine QEPs with second-, third-, and fourthorder EDs. Notably, the genuine QHPs with second-order ED exhibit a sixteenth-order DD. In the linear five-mode bosonic system described by the dynamical matrix M(5) lgiven in Eq. (44), with damping and/or amplification rates satisfying the conditions in Eqs. (45) and (46), the QEPs and QHPs observed in the dynamics of firstand second-order FOMs are summarized in Tab. 5. Notably, two inherited QEPs with second-order EDs give rise to genuine QEPs exhibiting second-, third-, and fourth-order EDs in various configurations within the spectral structure of second-order FOMs. Finally, for the five-mode bosonic system in the pyramid configuration, described by the dynamical matrix M(5) pgiven in Eq. (51) with damping and/or amplification rates satisfying the conditions in Eqs. (52) and (53), the QEPs and QHPs identified in the dynamics of firstand second-order FOMs are systematically summarized in Tab. 6. GenAccepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 23 Λi jΛr jMoments Moment Genuine and induced QHPs Genuine QHPs deg. Partial Partial Partial Partial QDP x QDP x QDP x QDP x QEP deg. QEP deg. QEP deg. QEP deg. γ+α±⟨ˆ b1⟩,⟨ˆ b3⟩ ⟨ ˆ B1⟩1 2x2 4x2 2x2 4x2 ⟨ˆ b4⟩,⟨ˆ b2⟩ ⟨ ˆ B2⟩1 −α±⟨ˆ b† 1⟩,⟨ˆ b† 3⟩ ⟨ ˆ B3⟩1 2x2 2x2 ⟨ˆ b† 4⟩,⟨ˆ b† 2⟩ ⟨ ˆ B4⟩1 2γ+2α±⟨ˆ b2 1⟩,⟨ˆ b2 3⟩ ⟨ ˆ B2 1⟩1 1x4 16x4 1x3 4x3 α−+α+⟨ˆ b1ˆ b3⟩2 + 2α∓⟨ˆ b2 2⟩,⟨ˆ b2 4⟩ ⟨ ˆ B2 2⟩1 1x4 1x3 α−+α+⟨ˆ b2ˆ b4⟩2 −2α±⟨ˆ b†2 1⟩,⟨ˆ b†2 3⟩ ⟨ ˆ B2 3⟩1 1x4 1x3 −α−−α+⟨ˆ b† 1ˆ b† 3⟩2 −2α∓⟨ˆ b†2 2⟩,⟨ˆ b†2 4⟩ ⟨ ˆ B2 4⟩1 1x4 1x3 −α−−α+⟨ˆ b† 2ˆ b† 4⟩2 −α±+α±⟨ˆ b† 1ˆ b1⟩,⟨ˆ b† 3ˆ b3⟩ ⟨ ˆ B3ˆ B1⟩2 2x4 1x4 6x4 ±α−∓α+⟨ˆ b† 1ˆ b3⟩,⟨ˆ b1ˆ b† 3⟩2 −α∓+α∓⟨ˆ b† 2ˆ b2⟩,⟨ˆ b† 4ˆ b4⟩ ⟨ ˆ B4ˆ B2⟩2 2x4 1x4 ∓α−±α+⟨ˆ b† 2ˆ b4⟩,⟨ˆ b† 4ˆ b2⟩2 α++α−⟨ˆ b1ˆ b2⟩,⟨ˆ b3ˆ b4⟩ ⟨ ˆ B1ˆ B2⟩2 2x4 1x4 2α±⟨ˆ b1ˆ b4⟩,⟨ˆ b3ˆ b2⟩2 α∓−α±⟨ˆ b† 1ˆ b2⟩,⟨ˆ b† 3ˆ b4⟩ ⟨ ˆ B3ˆ B2⟩2 2x4 1x4 α±−α±⟨ˆ b† 1ˆ b4⟩,⟨ˆ b† 3ˆ b2⟩2 α±−α∓⟨ˆ b† 2ˆ b1⟩,⟨ˆ b† 4ˆ b3⟩ ⟨ ˆ B4ˆ B1⟩2 2x4 1x4 α±−α±⟨ˆ b† 4ˆ b1⟩,⟨ˆ b† 2ˆ b3⟩2 −α+−α−⟨ˆ b† 1ˆ b† 2⟩,⟨ˆ b† 3ˆ b† 4⟩ ⟨ ˆ B3ˆ B4⟩2 2x4 1x4 −2α±⟨ˆ b† 1ˆ b† 4⟩,⟨ˆ b† 3ˆ b† 2⟩2 Table 3: Real and imaginary parts of the complex eigenfrequencies Λr j−iΛi jof the matrix M(4) l1 , given in Eq. (28) with Eq. (29), for the linear four-mode bosonic system with equal damping and/or amplification rates for neighbor modes derived from the equations for the FOMs up to second order. We note that α±(ζ)=α∓(−ζ)is used here. We have ˆ B1≡[ˆ b1,ˆ b3],ˆ B2≡[ˆ b2,ˆ b4],ˆ B3=ˆ B† 1, ˆ B4=ˆ B† 2. More details are given in the caption to Tab. 1. Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 24 Λi jΛr jMoments Moment Genuine and induced QHPs Genuine QHPs deg. Partial Partial Partial Partial QDP x QDP x QDP x QDP x QEP deg. QEP deg. QEP deg. QEP deg. γ+±α+⟨ˆ b1⟩,⟨ˆ b† 1⟩ ⟨ ˆ B1⟩1 1x2 2x2 1x2 2x2 ⟨ˆ b2⟩,⟨ˆ b† 2⟩ ⟨ ˆ B2⟩1 1x2 1x2 α−⟨ˆ b3⟩ ⟨ ˆ B3⟩1 2x1 2x1 2x1 2x1 ⟨ˆ b4⟩ ⟨ ˆ B4⟩1 −α−⟨ˆ b† 3⟩ ⟨ ˆ B5⟩1 2x1 2x1 2x1 2x1 ⟨ˆ b† 4⟩ ⟨ ˆ B6⟩1 2γ+±2α+⟨ˆ b2 1⟩,⟨ˆ b†2 1⟩ ⟨ ˆ B2 1⟩1 1x4 4x4 1x3 2x3 α+−α+⟨ˆ b† 1ˆ b1⟩2 + ±2α+⟨ˆ b2 2⟩,⟨ˆ b†2 2⟩ ⟨ ˆ B2 2⟩1 1x4 1x3 α+−α+⟨ˆ b† 2ˆ b2⟩2 ±α+⟨ˆ b1ˆ b2⟩,⟨ˆ b† 1ˆ b† 2⟩ ⟨ ˆ B1ˆ B2⟩2 2x4 1x4 1x4 α+−α+⟨ˆ b† 1ˆ b2⟩,⟨ˆ b1ˆ b† 2⟩2 α−±α+⟨ˆ b1ˆ b3⟩,⟨ˆ b† 1ˆ b3⟩ ⟨ ˆ B1ˆ B3⟩2 2x2 16x2 2x2 16x2 ⟨ˆ b1ˆ b4⟩,⟨ˆ b† 1ˆ b4⟩ ⟨ ˆ B1ˆ B4⟩2 2x2 2x2 α−±α+⟨ˆ b2ˆ b3⟩,⟨ˆ b† 2ˆ b3⟩ ⟨ ˆ B2ˆ B3⟩2 2x2 2x2 ⟨ˆ b2ˆ b4⟩,⟨ˆ b† 2ˆ b4⟩ ⟨ ˆ B2ˆ B4⟩2 2x2 2x2 −α−±α+⟨ˆ b1ˆ b† 3⟩,⟨ˆ b† 1ˆ b† 3⟩ ⟨ ˆ B1ˆ B5⟩2 2x2 2x2 ⟨ˆ b1ˆ b† 4⟩,⟨ˆ b† 1ˆ b† 4⟩ ⟨ ˆ B1ˆ B6⟩2 2x2 2x2 −α−±α+⟨ˆ b2ˆ b† 3⟩,⟨ˆ b† 2ˆ b† 3⟩ ⟨ ˆ B2ˆ B5⟩2 2x2 2x2 ⟨ˆ b2ˆ b† 4⟩,⟨ˆ b† 2ˆ b† 4⟩ ⟨ ˆ B2ˆ B6⟩2 2x2 2x2 2α−⟨ˆ b2 3⟩,⟨ˆ b2 4⟩ ⟨ ˆ B2 3⟩, ⟨ˆ B2 4⟩ 1 4x1 4x1 3x1 3x1 ⟨ˆ b3ˆ b4⟩ ⟨ ˆ B3ˆ B4⟩2 −2α−⟨ˆ b†2 3⟩,⟨ˆ b†2 4⟩ ⟨ ˆ B2 5⟩, ⟨ˆ B2 6⟩ 1 4x1 4x1 3x1 3x1 ⟨ˆ b† 3ˆ b† 4⟩ ⟨ ˆ B5ˆ B6⟩2 α−−α−⟨ˆ b† 3ˆ b3⟩,⟨ˆ b† 4ˆ b4⟩ ⟨ ˆ B5ˆ B3⟩, ⟨ˆ B6ˆ B4⟩ 2 8x1 8x1 4x1 4x1 ⟨ˆ b† 3ˆ b4⟩,⟨ˆ b3ˆ b† 4⟩ ⟨ ˆ B5ˆ B4⟩, ⟨ˆ B3ˆ B6⟩ 2 Table 4: Real and imaginary parts of the complex eigenfrequencies Λr j−iΛi jof the matrix M(4) l2 , given in Eq. (28) with Eq. (33) for the linear four-mode bosonic system with different damping and/or amplification rates of neighbor modes, as derived from the equations for the FOMs up to second order assuming α+= 0. We have ˆ Bj≡[ˆ bj,ˆ b† j]for j= 1,2,ˆ Bj≡[ˆ bj]for j= 3,4,ˆ B5=ˆ B† 3,ˆ B6=ˆ B† 4 and more details are given in the caption to Tab. 1. Accepted in Quantum 2025-12-07, click title to verify. Published under CC-BY 4.0. 25