Multiple quantum exceptional, diabolical, and hybrid points in multimode bosonic systems: II. Nonconventional PT -symmetric dynamics and unidirectional coupling

Multiple quantum exceptional, diabolical, and hybrid points in multimode bosonic systems: II. Nonconventional PT-symmetric dynamics and unidirectional coupling Jan Peˇ rina Jr.1, Kishore Thapliyal1, 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 We analyze the existence and degeneracies of quantum exceptional, diabolical, and hybrid points in simple bosonic systems - comprising up to six modes with damping and/or amplification - under two complementary scenarios to those of Ref. [1]: (i) nonconventional PT-symmetric dynamics confined to a subspace of the full Liouville space, and (ii) systems featuring unidirectional coupling. The system dynamics described by quadratic non-Hermitian Hamiltonians is governed by the Heisenberg-Langevin equations. Conditions for the observation of inherited quantum hybrid points with up to sixth-order exceptional and second-order diabolical degeneracies are revealed, though relevant only for short-time dynamics. This raises the question of whether higherorder inherited singularities exist in bosonic systems under general conditions. Nevertheless, for short times, unidirectional coupling of various types enables the concatenation of simple bosonic systems with secondand third-order exceptional degeneracies such that arbitrarily high exceptional degeneracies are reached. Methods for numerical identifying the quantum exceptional and hybrid points together with their degeneracies are addressed. Following Ref. [1] rich dynamics of second-order field-operator moments is analyzed from the point of view of the presence of exceptional and diabolical points and their degeneracies. 1 Introduction Non-Hermitian bosonic parity-time (PT) -symmetric systems exhibit a range of remarkable phenomena at and near quantum exceptional and hybrid points (QEPs and QHPs, respectively). Several studies [2,3,4,5,6,7] have shown that these singularities can enhance measurement precision beyond the classical limit and amplify effective system nonlinearities, leading to the generation of highly nonclassical and entangled states with tailored Jan Peˇ rina Jr.: [email protected] Kishore Thapliyal: [email protected] Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 1 arXiv:2405.01667v2 [quant-ph] 4 Nov 2025 properties [8,9,10,11]. We note that several studies (e.g., [12,13,14,15,16,17]) have demonstrated that in the linear regimes of quantum systems the measurement sensitivity boost may not occur. At present, investigations of QEPs in nonlinear quantum systems from the point of view of the measurement precision proceed to demonstrate the advantage of PT-symmetric dynamics in connection with the system nonlinearity that allows on its own to beat the classical limit. However, the intrinsic noise accompanying damping and amplification in PT-symmetric quantum systems must also be taken into account, as it typically degrades quantumness over extended timescales [18,19]. Despite this, highly squeezed and sub-Poissonian states remain achievable. Under various configurations of passive and active PT-symmetric bosonic systems, one can also generate entangled states exhibiting asymmetric steering and Bell nonlocality [20]. Moreover, these systems — and their field-operator moments (FOMs) [21] — can simulate complex many-body dynamics, including boundary effects [21,22]. On the applied side, optical switching via encirclement of a QHP has been demonstrated [23], and unidirectional light propagation — along with invisibility cloaking — has been realized in such platforms [24,25,26,27]. For a comprehensive overview of the rich physics around the spectral singularities of PT -symmetric systems, see the reviews in [28,29]. The strength of these effects depends in many cases on the order of the degeneracy of exceptional points (EPs): The higher-order the degeneracy, the more enhanced the processes become. This led us in Part I of Ref. [1] to the analysis of QEPs and QHPs of PT-symmetric bosonic systems with up to five modes considering different configurations under very general conditions. However, this analysis revealed only the systems with QHPs with the secondand third-order exceptional degeneracies (EDs) and second-order diabolical degeneracies (DDs), despite the fact that the bosonic systems with five modes were considered with the promise of observation of QHPs with the fifth-order ED. For this reason, we consider more general bosonic system compared to those analyzed in [1] to seek for the observation of higher-order EDs for inherited QEPs and QHPs. First, we weaken our requirements for the observation of PT-symmetric dynamics by considering only subspace(s) of the whole Liouville space of the statistical operators. We note here that, owing to the linearity of quantum mechanics, we can equivalently describe the system dynamics [30,31] in the Liouville space of the statistical operators and the complete space spanned by the operators of measurable quantities. We also note that the linearity gives the one-to-one correspondence between the subspaces of the above mentioned spaces. When such PT -symmetric-like behavior is restricted to only a subspace we refer to nonconventional PT-symmetric system behavior. Second, we admit in our analysis more general non-Hermitian PT -symmetric Hamiltonians. We recall here that, in Ref. [1], the non-Hermiticity of the investigated Hamiltonians originated only in the presence of damping and amplification terms, whose non-Hermiticity was ‘remedied’ by the presence of the Langevin fluctuating operator forces [32,31]. This guarantees the system evolution preserving the bosonic canonical commutation relations. Here, we consider also the Hamiltonians that describe bosonic systems with unidirectional coupling between the modes. The reason is that unidirectional coupling allows to concatenate two bosonic subsystems such that they keep their original eigenvalues. Moreover the original subspaces belonging to the same eigenvalues merge together which results in the increased EDs. This property gives rise to the method suggested and elaborated in Refs. [33,34] that provides QEPs with higher-order EDs. Realization of unidirectional propagation is based upon ring resonators [33,34]. However, we note that, apart from a complex experimental realization, unidirectional coupling is highly non-Hermitian and violates reciprocity of physical processes. Nevertheless, up to our best knowledge, this Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 2 is the only method for reaching QEPs and QHPs with high-order EDs for open bosonic systems under general conditions. However, as we show here, the quantum-mechanical consistency of the models with unidirectional coupling is guaranteed only for short times. Higher-order QEPs using unidirectional coupling were already realized in [35]. We note that if quantum field-amplitude fluctuations are completely neglected and pure states are considered, higher-order EPs in the Hamiltonian spectra can be observed relatively easily. For example, higher-order EPs were predicted in optomechanical [36] and cavity magnonic systems [37], those described by the Bose-Hubbard model [38], or photonic structures [39,40]. Higher-order EDs were also studied in Refs. [41,42]. They play significant role in amplification [43] and sensing [5] as well as speeding up entanglement generation [44]. However, the dynamics of open quantum systems brings new features and qualitatively modifies the systems’ behavior [45,46]. Increasing complexity of the bosonic systems also poses the question about identification of inherited QEPs and QHPs and the determination of their degeneracies. Whereas simple bosonic systems allow for analytical derivation of the eigenvalues and eigenvectors of their dynamical matrices, more complex bosonic systems admit only numerical treatment. In this case, we may numerically decompose a given dynamical matrix into its Jordan form that directly reveals QEPs and QHPs with their degeneracies. Alternatively, we may add a little ϵperturbation to any element of the dynamical matrix that can remove both EDs and DDs and even allow for distinguishing ED and DD. Some of the effects related to the presence of QEPs (quantum exceptional points) and QHPs (quantum hybrid points) with higher-order EDs (exceptional degeneracies) and DDs (diabolical degeneracies) are observed also in the behavior of higher-order FOMs (fieldoperator moments) [32,30,47]. We note that we refer to genuine QEPs and QHPs in the case of higher-order FOMs. This originates from the fact that the dynamics of nthorder FOMs is built, in certain sense, as a ‘multiplied dynamics’ of first-order FOMs with their inherited QEPs and QHPs. This ‘multiplied dynamics’ then naturally contains the genuine QEPs and QHPs with higher-order EDs and DDs, as it was discussed in [47]. In general, genuine QEPs with ED orders up to nth power of ED orders of inherited QEPs of the first-order FOMs are expected in the dynamics of nth-order FOMs (for examples, see Ref. [47]). The structure of genuine QEPs and QHPs in the dynamics of higher-order FOMs intimately depends on that of the inherited QEPs and QHPs found for the first-order FOMs. For this reason, we have explicitly revealed the structure of genuine QEPs and QHPs belonging to the second-order FOMs for the bosonic systems in the tables of Appendix B of Ref. [1]. They explicitly elucidate the relation between EDs and DDs of these genuine QEPs and QHPs and degeneracies of inherited QEPs and QHPs. We note that the induced QEPs and QHPs, which were introduced in Ref. [47], have their origin in the existence of identical or similar (related by the field commutation relations) FOMs in the formal construction of higher-order FOMs spaces and they further increase the multiplicity of spectral degeneracies. Nevertheless, they do not enrich the system dynamics. There exist some general properties of EDs and DDs of genuine QEPs and QHPs of bosonic systems independent of their configuration that we address here to complete the analysis of specific bosonic systems. The paper is organized as follows. Section 2is devoted to bosonic systems exhibiting nonconventional PT-symmetric dynamics, as an example, four-mode systems are analyzed. Section 3contains the analysis of bosonic systems with unidirectional coupling in linear configurations involving in turn from two to six modes. The dynamics of the two-mode bosonic system with unidirectional coupling and its applicability are analyzed in Sec. 4. A general analysis of genuine and induced QHPs in the dynamics of arbitrary-order FOMs Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 3 ˆa1 ˆa2 ˆa3 ˆa4 , κ , κ , κ , κ γ13 γ24 γ13 γ24 (a) ˆa1 ˆa2 ˆa3 ˆa4 , κ , κ , κ , κ , κ , κ γ13 γ24 γ13 γ24 (b) Figure 1: Schematic diagrams of the four-mode bosonic systems in (a) circular and (b) tetrahedral configurations that exhibit quantum exceptional points (QEPs) and quantum hybrid points (QHPs) with various exceptional degeneracies (EDs) and diabolical degeneracies (DDs) observed in their nonconventional PT-symmetric dynamics of field-operator moments (FOMs) of different orders. Strengths ϵand κcharacterize, respectively, the linear and nonlinear coupling between the modes, while γ, with subscripts indicating the mode number(s), are the damping or amplification rates, and the annihilation operators ˆaidentify the mode number via their subscripts. is given in Sec. 5. Section 6brings conclusions. Tables summarizing QEPs and QHPs— along with their eigenvalue EDs and DDs—for firstand second-order FOM dynamics are given in Appendix A. Appendix Bpresents the numerical methods used to identify these singularities and their degeneracies. In Appendix C, we analyze the properties of the Langevin operator forces in the unidirectional-coupling model, while Appendix D details the statistical characteristics of the two-mode bosonic system under unidirectional coupling. 2 Bosonic systems with bidirectional coupling and nonconventional PTsymmetric dynamics When seeking for QEPs and QHPs in simple bosonic systems, we have observed the situations in which the behavior typical for PT -symmetric systems occurs only in certain subspaces of the whole space spanned by the field operators and their moments. We speak about nonconventional PT -symmetric dynamics in these subspaces. We note that we may alternatively specify the corresponding subspaces in the Liouville space of statistical operators [30]. As the conditions for the observation of nonconventional PT-symmetric dynamics are less restrictive than those required for the usual PT-symmetric dynamics, we analyze here simple bosonic systems exhibiting this form of dynamics from the point of view of the occurrence of higher-order QEPs and QHPs. In the following, we consider in turn four-mode bosonic systems in their circular and tetrahedral configurations (see Fig. 1). Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 4 2.1 Circular configuration The Hamiltonian ˆ H4,cof a four-mode bosonic system in the circular configuration depicted in Fig. 1(a) takes the following form : ˆ H4,c=h¯hϵˆa† 1ˆa2+ ¯hϵˆa† 2ˆa3+ ¯hϵˆa† 3ˆa4+ ¯hϵˆa† 4ˆa1+ ¯hκˆa1ˆa2 +¯hκˆa2ˆa3+ ¯hκˆa3ˆa4+ ¯hκˆa4ˆa1] + H.c. (1) where ˆaj(ˆa† j) for j= 1,...,4denotes the annihilation (creation) operator of the jth mode, ϵ(κ) is the linear (nonlinear) coupling strength between the modes [48]. Symbol H.c. replaces the Hermitian-conjugated terms. Damping or amplification of mode jis described by the damping (amplification) rate γjand the corresponding Langevin stochastic operator forces, ˆ Ljand ˆ L† jthat occur in the dynamical Heisenberg-Langevin equations written below in Eq. (2). The Langevin stochastic operator forces are assumed to have the Markovian and Gaussian properties specific to the damping and amplification processes [47]. In Eq. (49), we present an example of firstand second-order correlation functions of the Langevin forces associated with damping (in mode 1) and amplification (in mode 2). Their presence in the Heisenberg-Langevin equations guarantees the fulfillment of the field-operator commutation relations. Moreover the properties of the Langevin stochastic operator forces are related to the damping (amplification) rates via the fluctuation-dissipation theorems [31,32]. The Heisenberg-Langevin equations corresponding to the Hamiltonian ˆ H4,cin Eq. (1) are written in the form: dˆ a dt =−iM(4) cˆ a+ˆ L,(2) where the vectors ˆ aof field operators and ˆ Lof the Langevin operator forces are given as ˆ a= [ˆ a1,ˆ a2,ˆ a3,ˆ a4]T≡hˆa1,ˆa† 1,ˆa2,ˆa† 2,ˆa3,ˆa† 3,ˆa4,ˆa† 4iTand ˆ L=hˆ L1,ˆ L† 1,ˆ L2,ˆ L† 2,ˆ L3,ˆ L† 3,ˆ L4,ˆ L† 4iT. The dynamical 8×8matrix M(4) cintroduced in Eq. (2) is derived in the form M(4) c=     −i˜γ1ξ0ξ ξ−i˜γ2ξ0 0ξ−i˜γ3ξ ξ0ξ−i˜γ4      .(3) The 2×2submatrices ˜γj,j= 1,...,4, and ξoccurring in Eq. (3) are defined as: ˜γj="γj/2 0 0γj/2#,ξ="ϵ κ −κ−ϵ#,(4) and γjstands for the damping or amplification rate of the mode j. The submatrices ˜γj,j= 1,...,4, and ξin Eq. (4) can simultaneously be diagonalized. As the submatrices ˜γjare linearly proportional to the identity matrix, the appropriate transformation does not change their form. On the other hand, it yields two eigenvalues λξ 1,2of the matrix ξwhich we identify by a new variable ξ. In this diagonalized form of Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 5 the above submatrices, the original 8×8dynamical matrix M(4) cattains the form M(4) c=                 −iγ10λξ 1000λξ 10 0−iγ10λξ 2000λξ 2 λξ 10−iγ20λξ 1000 0λξ 20−iγ20λξ 200 00λξ 10−iγ30λξ 10 000λξ 20−iγ30λξ 2 λξ 1000λξ 10−iγ40 0λξ 2000λξ 20−iγ4                 . In this form, the matrix M(4) cis decomposed into the direct sum of two 4×4matrices belonging to ξ=λξ 1and ξ=λξ 2. This allows us to effectively consider the dynamical matrix M(4) cin its 4×4block structure as if there are just numbers instead of 2×2submatrices. This effectively halves the dimension of the diagonalization procedure, allowing us to derive important analytical results. To reveal QEPs and QHPs, inspired by PT -symmetry and the results of Part I of Ref. [1], we apply the conditions γ1=γ3≡γ13, γ2=γ4≡γ24,(5) in the dynamical 4×4matrix M(4) cin Eq. (3). Then, we reveal its eigenvalues λM(4) c j: λM(4) c 1=−iγ13, λM(4) c 2=−iγ24, λM(4) c 3,4=−iγ+∓β. (6) The corresponding eigenvectors are derived as follows: yM(4) c 1= [−1,0,1,0]T, yM(4) c 2= [0,−1,0,1]T, yM(4) c 3=−2ξ χ∗,1,−2ξ χ∗,1T , yM(4) c 4=2ξ χ,1,2ξ χ,1T ,(7) where χ=iγ++β,β2= 4ξ2−γ2 −, and 4γ±=γ13 ±γ24. If β= 0 then λM(4) c 3=λM(4) c 4and also yM(4) c 3=yM(4) c 4. We note that the eigenvalues λM(4) c 3and λM(4) c 4share also their imaginary parts, which is important for the observation of PT -symmetric-like dynamics in a suitable interaction frame [49]. Provided that the system initial conditions are chosen such that only the eigenvalues λM(4) c 3and λM(4) c 4determine its dynamics, we observe a second-order QEP. This QEP changes into a QHP with secondorder ED and DD when the 8×8matrix M(4) cis analyzed (see below). The condition β= 0 assuming ξ2=ζ2=ϵ2−κ2[see below Eq. (10)] transforms into the formula κ2 ϵ2+γ2 − 4ϵ2= 1 (8) Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 6 (a) (b) Figure 2: Real parts λrof the eigenvalues (a) λM(4) c 3,4of the matrix M(4) c, given in Eq. (6), for the four-mode bosonic system in the circular configuration with different damping and/or amplification rates of neighbor modes and (b) λM(4) t 3,4of the matrix M(4) t, given in Eq. (17), for ξ=±ζfor the four-mode bosonic system in the tetrahedral configuration with the same damping and/or amplification rates of neighbor modes are drawn in the parameter space (κ/ϵ, γ−/ϵ). The dashed red curves indicate the values at the positions of the QHPs of the four-mode bosonic systems as given by Eq. (8). for an ellipse in the parameter space (κ/ϵ, γ−/ϵ)that identifies the positions of QHPs. Real parts of two eigenvalues λM(4) c j,j= 3,4, that form QHPs are plotted in this space in Fig. 2(a). The diagonalized 8×8dynamical matrix M(4) cis obtained with the help of the eigenvalues and eigenvectors given in Eqs. (6) and (7), respectively, and the eigenvalues λξ 1,2 and the eigenvectors yξ 1,2of the matrix ξ: λξ 1,2=∓ζ, (9) yξ 1,2=−ϵ∓ζ κ,1T ,(10) where ζ=√ϵ2−κ2. In detail, the eigenvalues ΛM(4) c jfor j= 1,...,8and the corresponding eigenvectors YM(4) c jof the full 8×8matrix M(4) care obtained using the general scheme applicable to the fields composed of, in general, nmodes. Relying on the results applicable to a general 2n×2ndimensional matrix M(n)presented in Part I of Ref. [1] (Appendix A), for which we have ΛM(n) 2j−1=λM(n) j(ξ=λξ 1), ΛM(n) 2j=λM(n) j(ξ=λξ 2),(11) 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       ,(12) j= 1, . . . , n, Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 7 the appropriate results are obtained assuming n= 4 and M(n)=M(4) c. In the basis in which the 8×8dynamical matrix M(4) cattains its diagonal form, the system dynamics is described by the new field operators ˆ b=hˆ b1,ˆ b† 1,ˆ b2,ˆ b† 2,ˆ b3,ˆ b4,ˆ b† 3,ˆ b† 4iT arising in this diagonalization. We note that the diagonalization and introduction of the new field operators in the two-mode bosonic system is explicitly given in Eq. (14) of Sec. II of Part I of Ref. [1]. We also note that the order of elements in the operator vector ˆ bis given by the numbering of the eigenvalues and the corresponding diagonalization transform. If the initial conditions allow to describe the complete system dynamics in terms of the field operators ˆ b3,ˆ b4,ˆ b† 3, and ˆ b† 4with equal damping or amplification rate γ+the firstand second-order FOMs exhibit in their dynamics QEPs and QHPs summarized in Tab. 2in Appendix A. 2.2 Tetrahedral configuration In the tetrahedral configuration depicted in Fig. 1(b), the Hamiltonian ˆ H4,tof four-mode bosonic system attains the form: ˆ H4,t=h¯hϵˆa† 1ˆa2+ ¯hϵˆa† 1ˆa3+ ¯hϵˆa† 1ˆa4+ ¯hϵˆa† 2ˆa3+ ¯hϵˆa† 2ˆa4 +¯hϵˆa† 3ˆa4+ ¯hκˆa1ˆa2+ ¯hκˆa1ˆa3+ ¯hκˆa1ˆa4+ ¯hκˆa2ˆa3 +¯hκˆa2ˆa4+ ¯hκˆa3ˆa4] + H.c. (13) The Heisenberg-Langevin equations corresponding to the Hamiltonian ˆ H4,tare derived in the form: dˆ a dt =−iM(4) tˆ a+ˆ L(14) using the following dynamical matrix M(4) t: M(4) t=     −i˜γ1ξ ξ ξ ξ−i˜γ2ξ ξ ξ ξ −i˜γ3ξ ξ ξ ξ −i˜γ4      .(15) In seeking QEPs, we assume equal damping and/or amplification rates of modes 1 and 2, and also of modes 3 and 4: γ1=γ2≡γ12, γ3=γ4≡γ34.(16) We note that, due to the symmetry, identical results are obtained when assuming equal damping and/or amplification rates of modes 1 and 3 and also of modes 2 and 4 [compare Eq. (5)]. Under these conditions, diagonalization of the 4×4dynamical matrix M(4) tin Eq. (15) leaves us with the following eigenvalues: λM(4) t 1=−iγ12 −ξ, λM(4) t 2=−iγ34 −ξ, λM(4) t 3,4=−iγ++ξ∓β. (17) Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 8 The corresponding eigenvectors are written as: yM(4) t 1= [−1,1,0,0]T, yM(4) t 2= [0,0,−1,1]T, yM(4) t 3=1−2iγ− χ− ,1−2iγ− χ− ,1,1T , yM(4) t 4=1−2iγ− χ+ ,1−2iγ− χ+ ,1,1T ,(18) and χ±=iγ−±β+ 2ξ,β2= 4ξ2−γ2 −, and 4γ±=γ12 ±γ34. Provided that β= 0, we have λM(4) t 3=λM(4) t 4and yM(4) t 3=yM(4) t 4as χ−=χ+. As the imaginary parts of eigenvalues λM(4) t 1,2differ from those of λM(4) t 3,4, the system can exhibit only the non-conventional PT-symmetric dynamics: If the system initial conditions are such that only the eigenvalues λM(4) t 3and λM(4) t 4suffice in describing its dynamics, we observe a second-order QEP for the 4×4dynamical matrix M(4) t. As the eigenvalues λM(4) t jin Eq. (17) show the linear dependence on ξ, the diabolical second-order degeneracy of the 8×8dynamical matrix M(4) c, originating in the form of the eigenvalues λM(4) c jin Eq. (6), is not observed in the tetrahedral configuration. Instead, for β= 0, we find one secondorder QEP for ξ=ζand another second-order QEP for ξ=−ζ. These QEPs occur at the positions described in Eq. (8) in the parameter space (κ/ϵ, γ−/ϵ). Real parts of four eigenvalues λM(4) t 3,4for ξ=±ζthat build two QEPs are drawn in this space in Fig. 2(b). In the basis with the diagonal 8×8dynamical matrix M(4) t, the system dynamics is described by the new field operators ˆ b=hˆ b1,ˆ b† 1,ˆ b2,ˆ b† 2,ˆ b3,ˆ b4,ˆ b† 4,ˆ b† 3iT. The corresponding eigenvalues and eigenvectors are discussed in general in Appendix A of Ref. [1]. The firstand second-order FOMs exhibit in their dynamics QEPs and QHPs provided in Tab. 2of Appendix A. 3 Concatenated bosonic systems with unidirectional coupling: Higherorder quantum exceptional points on demand The analysis of PT -symmetric bosonic systems with up to five modes in their linear, circular, tetrahedron, and pyramid configurations presented above and in Part I of Ref. [1] revealed only the inherited QEPs and QHPs with secondand third-order EDs. That is why, we extend our analysis to more general non-Hermitian Hamiltonians that involve unidirectional coupling between the modes. Whereas the non-Hermiticity of the above discussed systems is given solely by the presence of damping and/or amplification, the Hamiltonians with unidirectional coupling are non-Hermitian per se. Despite their non-Hermiticity, they have direct effective physical implementations based upon counter-directional field propagation and mutual scattering [33,34] or nonlinear Kerr interaction [50]. Moreover, it has been shown in Ref. [51] that exponential improvement of measurement precision can be reached in QEPs in systems with unidirectional coupling. It was shown in Refs. [33,34] that using a specific unidirectional coupling of two field modes belonging to different PT -symmetric systems with QEPs, the combined system exhibits a QEP with ED given as the sum of those of the constituting systems. This opens the door for observing inherited QEPs with EDs of orders higher than three. Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 9 The solution in Eq. (50) can be recast into a simpler form written for the annihilation operators ˆa1and ˆa2: "ˆa1(t) ˆa2(t)#=U(t)"ˆa1(0) ˆa2(0) #+V(t)"ˆa† 1(0) ˆa† 2(0) #+"ˆ f1(t) ˆ f2(t)#. (54) The elements of the matrices Uand Vare defined as Uj,k(t) = P2j−1,2k−1(t, 0) and Vj,k(t)=P2j−1,2k(t, 0) for j, k = 1,2, and we also have ˆ fj(t) = ˆ F2j−1(t)for j= 1,2. Using the eigenvalues and eigenvectors written in Eqs. (9), (10), (21), and (22), we arrive at the formulas specific to our model: U(t) = "µ(t) 0 −iϵs(t) γ 1 µ(t)#,V(t) = "0 0 −iκs(t) γ0#,(55) ⟨ˆ F(t)ˆ F†T(t)⟩="F1(t)F12(t) F∗T 12 (t)F2(t)#,(56) F1(t) = "1−µ2(t) 0 0 0 #,F12(t) = iσ(t) 2γ"ϵ−κ 0 0 #, F2(t) = s(t)−2γt 2γ2"ϵ2−ϵκ −ϵκ −κ2#+1 µ2(t)−1"0 0 0 1 #, (57) where µ(t) = exp(−γt),σ(t) = exp(−2γt)−1 + 2γt, and s(t) = sinh(2γt)using the hyperbolic sinus function. To check consistency of the model with unidirectional coupling, we determine the mean values of equal time field-operators commutation relations. Compared to the usual canonical commutation relations ⟨[ˆaj(t),ˆak(t)]⟩= 0,⟨[ˆaj(t),ˆa† k(t)]⟩=δjk, and ⟨[ˆa† j(t),ˆa† k(t)]⟩= 0 for j, k = 1,2, the following two relations are found: ⟨[ˆa2(t),ˆa† 2(t)]⟩= 1 + (ϵ2−κ2)ϕ(t) 2γ2, ⟨[ˆa1(t),ˆa† 2(t)]⟩=−iϵψ(t) 2γ,(58) where ψ(t) = exp(−2γt)−1−2γt and ϕ(t) = exp(2γt)−1−2γt. For short times tassuming t≪1/γ we have ⟨[ˆa2(t),ˆa† 2(t)]⟩= 1 + (ϵ2−κ2)t2and ⟨[ˆa1(t),ˆa† 2(t)]⟩=−2iϵt. Thus, we additionally require t≪1/ϵ and t≪1/√ϵ2−κ2. As we usually assume in PT -symmetric systems that κ≤ϵ, we are left with the following conditions for applicability of the model with unidirectional coupling: t≪min 1 γ,1 ϵ.(59) The form of the commutation relations given in Eq. (58) and the ensuing restricted validity of the model poses the question about possible corrections of the model using suitable properties of the reservoir Langevin operator forces. Similarly as it is done when damping and amplification are introduced into the Heisenberg equations (the Wigner– Weisskopf model of damping, see Ref. [32]). However, as discussed in detail in Appendix C, this approach is not successful. Also, in Appendix Dthe properties of the modes are analyzed in the framework of the Gaussian states, their nonclassicality depths and logarithmic negativity are determined. Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 16 k z }| { Λ1. . . Λ1 | {z } k1 . . . . . . . . . ΛΣΛ. . . ΛΣΛ | {z } kΣΛ ˆ B1,1. . . ˆ B1,1......... ˆ BΣΛ,1. . . ˆ BΣΛ,1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ˆ B1,n1. . . ˆ B1,n1......... ˆ BΣΛ,nΣΛ. . . ˆ BΣΛ,nΣΛ Figure 4: Scheme for constructing general kth-order FOMs for an n-mode bosonic system with ΣΛ eigenvalues Λ1,...,ΛΣΛhaving n1, . . . , nΣΛcoalescing eigenvectors, i.e. PΣΛ j=1 nj= 2n. The field operators ˆ bland ˆ b† lfor l= 1, . . . , n in the diagonalized basis form elements of the field-operator vectors ˆ B1,..., ˆ BΣΛcomposed of in turn n1, . . . , nΣΛelements. The kth-order FOMs are divided into the groups identified with the vectors k≡[k1, . . . , kΣΛ],PΣΛ j=1 kj=k, containing FOMs of the form ⟨ˆ Bk1 1··· ˆ BkΣΛ ΣΛ⟩. These results point out at specific properties of the two-mode bosonic system with unidirectional coupling applicable only under the conditions given in Eq. (59). These results lead us to the conclusion that this method for concatenating simpler bosonic systems with QEPs based on unidirectional coupling to arrive at QEPs with higher-order ED has limitations. The question how to obtain higher-order inherited QEPs in bosonic systems without these limitations and under general conditions is open. 5 Higher-order hybrid points revealed by field-operator moments In the last section, we derive general formulas that give us the orders of EDs and DDs of QHPs that occur in the dynamics of higher-order FOMs. We note that the structure of higher-order FOM spaces mapped onto suitable lattices and its influence to system’s properties has been analyzed in detail in Ref. [21] for twoand three-mode bosonic systems. Let us consider a bosonic system composed of nmodes and described by an appropriate quadratic non-Hermitian Hamiltonian. Such a system is described by the 2nannihilation and creation operators, and the 2n×2ndynamical matrix M(n)of its Heisenberg-Langevin equations has 2neigenvalues. Let us fix a position in the system parameter space. After diagonalization of the 2n×2ndynamical matrix M(n), we identify the eigenvalues Λj with different eigenvectors and count their degeneration numbers njaccording to the number of coalescing eigenvectors (from the original maximal Hilbert space). We note that the eigenvalues Λjcan coincide, leading to diabolical degeneracies. The corresponding eigenvalue structure is illustrated in Fig. 4for reference. Denoting the number of such eigenvalues by ΣΛ, we have ΣΛ X j=1 nj= 2n. (60) We assign, to any eigenvalue Λj, an operator vector ˆ Bjthat encompasses the diagonalized field operators ˆ bland ˆ b† l(l= 1, . . . , n) that are associated with the coalescing vectors of this eigenvalue, as shown in Fig. 4. Now we consider the dynamics of kth-order FOMs. These kth-order FOMs can formally be expressed in their general form using the elements of the above-defined operator vectors ˆ Bjas ⟨QΣΛ j=1 ˆ Bkj j⟩, where the nonnegative integers kjobey PΣΛ j=1 kj=k. The Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 17 vector kdefined as [k1, k2, . . . , kΣΛ]then points out at a possible QHP with the complex eigenfrequency Λgiven as Λ = ΣΛ X j=1 kjΛj.(61) The degrees dk eof ED and dk dof DD of this QHP are given by the following formulas (compare the scheme in Fig. 4): dk e= ΣΛ Y j=1 nkj j,(62) dk d=k! QΣΛ j=1 kj!.(63) The number N(k) Bof different FOMs contributing to DD observed for the kth-order FOMs (see the right-parts of the third main columns of Tabs. II—IV in Appendix A and Tabs. I—IV in Appendix B of Ref. [1]) is determined using the combination number (permutations with repetition, kballs in ΣΛdrawers) as: N(k) B= ΣΛ+k−1 k!= ΣΛ+k−1 ΣΛ−1!.(64) The number N(k) bof different kth-order FOMs written in the operators ˆ bland ˆ b† l, that are the elements of the vectors ˆ Bj, and belonging to a fixed vector kis expressed as: N(k) b= ΣΛ Y l=1 nl+kl−1 kl!= ΣΛ Y l=1 nl+kl−1 nl−1!.(65) We note that a formula similar to Eq. (64) gives the number N(˜ k) bof different ˜ kth-order FOMs written in the operators ˆ bland ˆ b† lfor l= 1 . . . n,⟨ˆ b˜ k1 1···ˆ b˜ kl l···ˆ b†˜ kl+1 1···ˆ b†˜ k2l l⟩for ˜ k= [˜ k1,...,˜ k2l]and P2l j=1 ˜ kj=˜ k(permutations with repetition, ˜ kballs in 2ndrawers): N(˜ k) b= 2n+˜ k−1 ˜ k!= 2n+˜ k−1 2n−1!.(66) These FOMs are explicitly written in the left-parts of the third main columns of Tabs. II— IV in Appendix Aand Tabs. I—IV in Appendix B of Ref. [1]. To demonstrate the general approach and formulas, we apply them to the four-mode bosonic system in the circular configuration analyzed in Sec. 2in the regime of nonconventional PT-symmetric dynamics. Its firstand second-order FOMs and the revealed QHPs with their degeneracies are given in Tab. 2in Appendix A. The general approach gives us Tab. 1that allows to easily derive the content of Tab. 2in Appendix A: Two lines belonging to Λi j=γ+are obtained assuming k= (0,0,0,0,1,0) and (0,0,0,0,0,1). For Λi j= 2γ+the entries to Tab. 2are in turn generated assuming k= (0,0,0,0,1,1), (0,0,0,0,2,0), and (0,0,0,0,0,2). The orders of the corresponding EDs and DDs are derived using Eqs. (62) and (63), respectively. These results can be used for detailed discussions of genuine and induced QHPs and their degeneracies for FOMs of arbitrary orders and considering different systems with their specific inherited QHPs observed in the dynamics of field operators governed by the Heisenberg-Langevin equations. The occurrence of a QHP with nkth ED and second-order DD in the dynamics of kth-order FOMs in a bosonic system with an inherited QHP with nth-order ED and second-order DD is probably the most valuable result. Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 18 −iΛjnjˆ Bjkj γ13 1ˆ B1≡[ˆ b1]k1= 0 γ13 1ˆ B2≡[ˆ b† 1]k2= 0 γ24 1ˆ B3≡[ˆ b2]k3= 0 γ24 1ˆ B4≡[ˆ b† 2]k4= 0 γ+2ˆ B5≡[ˆ b3,ˆ b† 3]k5 γ+2ˆ B6≡[ˆ b4,ˆ b† 4]k6 PΣΛ j=1 Λj= Λ PΣΛ j=1 nj= 2n⟨QΣΛ j=1 ˆ Bkj j⟩PΣΛ j=1 kj=k Table 1: Eigenvalues Λj, their degeneracies nj, the corresponding field-operator vectors ˆ Bjdefined in the third column of the table, and their varying powers kjfor the four-mode bosonic system in the circular configuration analyzed in Tab. 2in Appendix Aunder the condition in Eq. (8). QHPs are uniquely identified by powers of k5and k6. 6 Conclusions Quantum exceptional, diabolical, and hybrid points were analyzed in simple bosonic systems with unusual properties. The bosonic systems exhibiting nonconventional PT - symmetric dynamics characterized by the observation of quantum exceptional and hybrid points only in certain subspace(s) of the whole system Liouville space were revealed and their properties investigated. Nonconventional second-order inherited quantum exceptional and hybrid points were identified in four-mode bosonic systems. Applying the method of concatenating simple bosonic systems via unidirectional coupling, the conditions for the observation of up to sixth-order inherited quantum exceptional and hybrid points were found. However, by analyzing the behavior of the two-mode bosonic system with unidirectional coupling we have shown that the bosonic systems with unidirectional coupling are applicable only for short times in which they ensure the physically consistent behavior. In short times, more complex bosonic systems with diverse structures and arbitrary-order exceptional degeneracies can be built by concatenating simple bosonic systems via unidirectional coupling of several types. Nevertheless, tailoring the properties of the Langevin operator stochastic forces in systems with unidirectional coupling does not enable extending their applicability to arbitrary times, in contrast to how the damping and amplification are consistently described. The operation of two numerical methods for the identification of quantum exceptional points and their degeneracies, namely (1) the transformation of a dynamical matrix into its Jordan form and (2) the introduction of a suitable perturbation δinto the dynamical matrix and its subsequent eigenvalue analysis, was demonstrated analytically in two-mode bosonic systems. The exceptional and diabolical degeneracies of inherited quantum hybrid points were used to derive higher-order degeneracies observed in the dynamics of higher-order fieldoperator moments. The quantum exceptional and hybrid points of second-order fieldoperator moments were summarized in tables that evidence a rich dynamics of the fieldoperator moments. Numbers of genuine and induced quantum hybrid points and their exceptional and diabolical degeneracies were expressed as they depend on the order of the field-operator moments in the general form using parameters of the inherited quantum exceptional and hybrid points and their degeneracies. The presented analysis considerably broadens the knowledge of bosonic systems with PT -symmetry exhibiting well-accepted physical behavior, though it does not reveal bosonic Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 19 systems with higher-order exceptional and hybrid singularities found under general conditions (e.g. long times). The present analysis, together with the results of Ref. [1], significantly advances our understanding of bosonic systems with PT -symmetry and their spectral singularities. It demonstrates that observing higher-order exceptional and hybrid singularities in physically well-behaved (i.e., real) bosonic systems at arbitrary times remains a challenging task. 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 QEPs and QHPs in firstand second-order FOMs spaces for uniand bidirectional coupling Considering both unidirectional and bidirectional coupling schemes discussed in the main text, we provide tables listing the QEPs and QHPs, along with their associated degeneracies, as found in firstand second-order FOM spaces. Based on the inherited QEPs and QHPs and their degeneracies observed in the dynamics of first-order FOMs, we construct the genuine and induced QEPs and QHPs in the second-order FOM dynamics, following the scheme outlined in Ref. [47]. This analysis complements the results presented in Part I of Ref. [1] and relies on the general formulas for spectral degeneracies and their multiplicities provided in Sec. V. The tables below illustrate the diversity and richness of spectral degeneracies that can emerge in simple bosonic models. The circular four-mode bosonic system described by the dynamical matrix M(4) cgiven in Eq. (3) with different damping and/or amplification rates of neighbor modes [see Eq. (5)] exhibits nonconventional PT-symmetric dynamics. Its QEPs and QHPs found in the dynamics of the firstand second-order FOMs are given in Tab. 2. Considering another system giving the QEPs and QHPs with nonconventional PT - symmetric dynamics — the four-mode bosonic system in tetrahedral configuration with the dynamical matrix M(4) tgiven in Eq. (15) and equal damping and/or amplification rates of neighbor modes [see Eq. (16)], we reveal the QEPs and QHPs belonging to the dynamics of the firstand second-order FOMs as summarized in Tab. 3. Table 4presents the QEPs and QHPs observed in the dynamics of firstand secondorder FOMs for the n-mode bosonic models with n= 3, ..., 6and unidirectional coupling. This includes systems characterized by the dynamical matrices M(1+2) uin Eq. (24) under the condition in Eq. (25), M(2+2) uin Eq. (30) under the condition in Eq. (32), M(2+3) u in Eq. (35) under the condition in Eqs. (37) and (40), and M(3+3) uin Eq. (42) under the condition in Eq. (43). In Tabs. 2—4, we can see the most typical feature of QEPs and QHPs in higher-order FOMs: The second (n-th) -order ED observed in first-order FOMs gives raise in Tabs. 2 and 3(4) the fourth (n2-th) -order ED in second-order FOMs. These are examples of the general rule that assigns nkED in k-th-order FOMs provided that there occurs n-th-order ED in first-order FOMs [47]. Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 20 Λ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. γ+±β⟨ˆ b3⟩,⟨ˆ b† 3⟩ ⟨ ˆ B5⟩1 1x2 2x2 1x2 2x2 ⟨ˆ b4⟩,⟨ˆ b† 4⟩ ⟨ ˆ B6⟩1 1x2 1x2 2γ+±2β⟨ˆ b3ˆ b4⟩, ⟨ˆ b† 3ˆ b† 4⟩ ⟨ ˆ B5 ˆ B6⟩2 2x4 4x4 1x4 1x4 β−β⟨ˆ b† 3ˆ b4⟩2 + β−β⟨ˆ b3ˆ b† 4⟩2 ±2β⟨ˆ b2 3⟩, ⟨ˆ b†2 3⟩ ⟨ ˆ B2 5⟩1 1x4 1x3 2x3 β−β⟨ˆ b† 3ˆ b3⟩2 ±2β⟨ˆ b2 4⟩, ⟨ˆ b†2 4⟩ ⟨ ˆ B2 6⟩1 1x4 1x3 β−β⟨ˆ b† 4ˆ b4⟩2 Table 2: Real and imaginary parts of the complex eigenfrequencies Λr j−iΛi jof the matrix M(4) c, given in Eq. (3) for the four-mode bosonic system in the circular configuration with different damping and/or amplification rates of neighbor modes valid for the regime of nonconventional PT -symmetric dynamics and derived from the equations for the FOMs up to second order. The corresponding moments written in the ‘diagonalized’ field operators involving the operators ˆ b3,ˆ b† 3,ˆ b4, and ˆ b† 4are given together with their degeneracies (deg.) coming from different possible relative 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= 5,6are defined in the rows written for Λi j=γ+devoted to the first-order FOMs, i.e. ˆ B5≡[ˆ b3,ˆ b† 3],ˆ B6≡[ˆ b4,ˆ b† 4]. Symbol ˆ Bjˆ Bk,j, k = 5,6, stands for the tensor product that gives 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. B Numerical identification of exceptional points and their degeneracies The performance of two different numerical approaches to reveal exceptional points and their orders of exceptional degeneracies is demonstrated analytically considering the simplest two-mode bosonic systems, one with the usual bidirectional coupling, the other with unidirectional coupling. B.1 Jordan canonical form of a general matrix The Jordan form JMof a matrix Mcontains nonzero elements on the diagonal and also the nearest upper diagonal. It has the dimension of the matrix M. Some elements at the nearest upper diagonal equal one if the matrix Mis non-diagonalizable. The neighbor elements equal to one form groups. The number of elements in a given group gives the order of ED (equal to the number of elements + 1) related to the corresponding eigenvalue. For example and considering the dynamical matrix M(2), given in Eq. (29), belonging to the two-mode bosonic system with usual coupling under the condition β= 0, i.e., where a second-order QEP occurs, the Jordan form JM(2) and the corresponding similarity transformation SM(2) such that M(2) =SM(2) JM(2) S−1 M(2) (67) take the form: JM(2) ="−iγ+1 0−iγ+#,SM(2) ="i−1/γ− 1 0 #.(68) Accepted in Quantum 2025-11-04, 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. γ+−ζ±β⟨ˆ b3⟩,⟨ˆ b† 4⟩ ⟨ ˆ B5⟩1 1x2 1x2 1x2 1x2 ζ±β⟨ˆ b4⟩,⟨ˆ b† 3⟩ ⟨ ˆ B6⟩1 1x2 1x2 1x2 1x2 2γ+±2β⟨ˆ b3ˆ b4⟩, ⟨ˆ b† 4ˆ b† 3⟩ ⟨ ˆ B5 ˆ B6⟩2 2x4 2x4 1x4 1x4 ±2ζ⟨ˆ b3ˆ b† 3⟩, ⟨ˆ b† 4ˆ b4⟩ 2 −2ζ±2β⟨ˆ b2 3⟩, ⟨ˆ b†2 4⟩ ⟨ ˆ B2 5⟩1 1x4 1x4 1x3 1x3 −2ζ⟨ˆ b† 4ˆ b3⟩2 2ζ±2β⟨ˆ b2 4⟩, ⟨ˆ b†2 3⟩ ⟨ ˆ B2 6⟩1 1x4 1x4 1x3 1x3 2ζ⟨ˆ b† 3ˆ b4⟩2 Table 3: Real and imaginary parts of the complex eigenfrequencies Λr j−iΛi jof the matrix M(4) t, given in Eq. (15), for the four-mode bosonic system in the tetrahedral configuration valid for the regime of nonconventional PT -symmetric dynamics and derived from the equations for the FOMs up to second order. We have ˆ B5≡[ˆ b3,ˆ b† 4],ˆ B6≡[ˆ b4,ˆ b† 3], and more details are given in the caption to Tab. 2. In Eq. (68), the element (1,2) of the matrix JM(2) , equal to 1, identifies a second-order QEP. Similarly, when analyzing the two-mode bosonic system with unidirectional coupling and the dynamical matrix M(1+1) u, written in Eq. (20), under the condition γ1=γ2 guaranteeing the existence of a second-order QEP, we arrive at: JM(1+1) ="−iγ2/2 1 0−iγ2/2#,SM(1+1) ="0 1/ξ 1 0 #.(69) B.2 Perturbation of a dynamical matrix The second approach is based upon introducing a suitable small perturbation on a dynamical matrix Mthat can remove both EDs and DDs. In this approach, the eigenvalues of the perturbed dynamical matrix Mδare determined and the degeneracies are removed as the degenerated eigenvalues split and then gradually diverge with the increasing perturbation. Perturbation can also be used for characterizing ED [57,58,59]. Identifying EDs, we demonstrate different kinds of the influence of perturbation δon the eigenvalues of the matrix Mconsidering the two-mode systems with the usual and unidirectional couplings and different positions of the perturbation δinside the matrix M. Specifically, we demonstrate splitting of the eigenvalues in their real and/or imaginary parts and splitting proportional to √δand δwhen a QEP with second-order ED is disturbed. We quantify the strength of the perturbation by the overlap Fof the normalized eigenvectors y1and y2that become gradually distinguishable as the perturbation δ increases F=|⟨y2|y1⟩| p⟨y1|y1⟩⟨y2|y2⟩(70) and symbol ⟨|⟩ stands for the scalar product of complex vectors. 1. A two-mode system with bidirectional coupling and β= 0 described by the perturbed Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 22 dynamical matrix M(2) δ,1=M(2) +"δ0 0 0 #.(71) The eigenvalues and the corresponding eigenvectors take, respectively, the form: λM(2) δ,1 1,2=−iγ++δ 2∓sδδ 4−iγ−(72) and yM(2) δ,1 1,2="i−δ∓pδ(δ−4iγ−) 2γ− ,1#T .(73) 2. A two-mode system with bidirectional coupling and β= 0 described by the perturbed dynamical matrix M(2) δ,2=M(2) +"0−δ 0 0 #.(74) The eigenvalues and the corresponding eigenvectors are derived, respectively, as follows: λM(2) δ,2 1,2=−iγ+∓qδγ−(75) and yM(2) δ,2 1,2="i±sδ γ− ,1#T .(76) 3. A two-mode system with unidirectional coupling and γ1=γ2described by the perturbed dynamical matrix M(1+1) u,δ =M(1+1) u+"δ0 0 0 #.(77) The eigenvalues and the corresponding eigenvectors are obtained, respectively, in the form: λM(1+1) u,δ 1=−iγ2/2, λM(1+1) u,δ 2=−iγ2/2+δ(78) and yM(1+1) u,δ 1= [0,1]T, yM(1+1) u,δ 2=δ ξ,1T .(79) The real and imaginary parts of the eigenvalues λj,j= 1,2, from Eqs. (72), (75), and (78) are plotted in Figs. 5(a,b) as they depend on the perturbation δ. The perturbation δ in general disturbs more strongly the system with the usual coupling, as it is apparent both from the graphs of the eigenvalues and the overlap Fof eigenvectors shown in Fig. 5(c). Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 23 (a) 0.00 0.02 0.04 0.06 0.08 0.10 δ 0.2 0.1 0.0 0.1 0.2 λr 1 2 1 2 1 2 (b) 0.00 0.02 0.04 0.06 0.08 0.10 δ 0.2 0.1 0.0 0.1 0.2 λi 2 1 1,2 1,2 (c) 0.02 0.04 0.06 0.08 0.10 δ 0.92 0.94 0.96 0.98 1.00 F Figure 5: (a) Real λrand (b) imaginary λiparts of eigenvalues λ1,2of the dynamical matrices M(2) δ,1 (γ+= 0,γ−= 0.5, blue solid curves), M(2) δ,2(γ+= 0,γ−= 0.5, red dashed curves), and M(1+1) u,δ (γ2= 0, green dot-dashed curves) given in turn in Eqs. (71), (74), and (77) as they depend on perturbation parameter δ. In (c) the overlap Fof the corresponding eigenvectors is plotted. In (a) and (b), the numbers denote the curves of the corresponding eigenvalues. In the above-discussed cases, the second-order DD present in both two-mode systems was not modified by the perturbation δbecause of the structure of these systems. However, suitable positioning of the perturbation δinside a dynamical matrix Mmay also result in revealing the DDs. As the DD is embedded in the 2×2matrix ξgiven in Eq. (4), the perturbation δhas to affect this matrix. We consider two kinds of perturbation in the two-mode system with the dynamical matrix M(2): The first one splits all eigenvalues Λj,j= 1,...,4, in their real parts, whereas the second one distinguishes two eigenvalues in their real parts and the remaining two eigenvalues in their imaginary parts. We note that the perturbation δprimarily modifies the eigenvalues of the matrix ξ, which removes the DD. Secondarily, as the eigenvalues of ξdetermined for nonzero δdiffer from those valid for δ= 0, the conditions for having a QEP of the matrix Mchange and the original setting of the system parameters for a QEP is lost and so also the corresponding ED is lost. 1. A two-mode system with bidirectional coupling and β= 0 described by the dynamical matrix M(2) where ξδ,1=ξ+"δ0 0 0 #.(80) The eigenvalues and the corresponding eigenvectors of matrix ξδ,1are written, respectively, as: λξδ,1 1,2=δ 2∓sζ2+δϵ+δ 4(81) and yξδ,1 1,2= −ϵ+λξδ,1 1,2 κ,1  T .(82) 2. A two-mode system with bidirectinal coupling and β= 0 described by the dynamical matrix M(2) where ξδ,2=ξ+"δ δ 0 0 #.(83) The eigenvalues and the corresponding eigenvectors of matrix ξδ,2are obtained, respectively, as: λξδ,2 1,2=δ 2∓sζ2+δϵ−κ+δ 4(84) Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 24 (a) 0.00 0.02 0.04 0.06 0.08 0.10 δ -0.4 -0.2 0.0 0.2 0.4 Λr 3 1 2 4 1 3,4 2 (b) 0.00 0.02 0.04 0.06 0.08 0.10 δ -0.2 -0.1 0.0 0.1 0.2 Λi 1,2,3,4 3 1,2 4 Figure 6: (a) Real Λrand (b) imaginary Λiparts of the eigenvalues Λ1,...,4of the dynamical matrix M(2) involving ξδ,1(blue solid curves) and ξδ,2(red dashed curves) given in Eqs. (80) and (83), respectively, as they depend on perturbation parameter δ;γ+= 0,β= 0,κ/ϵ = 1/2. In (a) and (b), the numbers denote the curves of the corresponding eigenvalues. and yξδ,2 1,2= −ϵ+λξδ,2+δ 1,2 κ,1  T .(85) The looked for eigenvalues and eigenvectors are then reached using the eigenvalues and eigenvectors of the 2×2matrix M(2) written in Eq. (29): λM(2) 1,2=−iγ+∓β(86) and yM(2) 1,2=−iγ−±β ξ,1T ,(87) where 4γ±=γ1±γ2and β2=ξ2−γ2 −. The real and imaginary parts of the eigenvalues Λj,j= 1,...,4, of the 4×4matrix M(2) are plotted in Fig. 6. Finally, we mention a specific case of the perturbation δthat removes DD but keeps one from the two originally diabolically-degenerated QEPs present in the system. The dynamical matrix M(2) of the two-mode bosonic system with bidirectional coupling is perturbed in the following way: ˜γδ,3 1=˜γ1+"δ0 0 0 #.(88) The eigenvalues ΛM(2) δ,3of the dynamical matrix in Eq. (88) are derived in the form: ΛM(2) δ,3 1,2=−iγ+∓β, ΛM(2) δ,3 3,4=−iγ++δ 2∓sβ2+δδ 4−iγ−.(89) The corresponding eigenvectors are expressed as follows: YM(2) δ,3 1,2=h0,iγ−±β ϵ,−κ ϵ,1iT, YM(2) δ,3 3,4="iγ−−Λ M(2) δ,3 3,4 κ,0,−ϵ κ,1#T .(90) Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 25 Λ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. γ+±β1,..., ⟨ˆ b1⟩,⟨ˆ b† 1⟩,...,⟨ˆ bn/2⟩,⟨ˆ b† n/2⟩ ⟨ ˆ B1⟩1 1xn 2xn 1xn 2xn ±βn⟨ˆ bn/2+1⟩,⟨ˆ b† n/2+1⟩,...,⟨ˆ bn⟩,⟨ˆ b† n⟩ ⟨ ˆ B2⟩1 1xn 1xn 2γ+±(βk+βl)⟨ˆ bkˆ bn/2+l⟩,⟨ˆ b† kˆ b† n/2+l⟩ ⟨ ˆ B1 ˆ B2⟩2 2xn24xn21xn21xn2 βl−βk⟨ˆ b† kˆ bn/2+l⟩2 + βk−βl⟨ˆ bkˆ b† n/2+l⟩2 k, l = 1, . . . n/2 ±2βk⟨ˆ b2 k⟩,⟨ˆ b†2 k⟩ ⟨ ˆ B2 1⟩1 1xn21x (n+ 1)n/2 2x (n+ 1)n/2 k= 1, . . . n/2 ±(βk+βl)⟨ˆ bkˆ bl⟩,⟨ˆ b† kˆ b† l⟩2 k, l = 1, . . . n/2, l < k βl−βk⟨ˆ b† kˆ bl⟩2 k, l = 1, . . . n/2 ±2βk⟨ˆ b2 n/2+k⟩,⟨ˆ b†2 n/2+k⟩ ⟨ ˆ B2 2⟩1 1xn21x (n+ 1)n/2 k= 1, . . . n/2 ±(βk+βl)⟨ˆ bn/2+kˆ bn/2+l⟩,⟨ˆ b† n/2+kˆ b† n/2+l⟩2 k, l = 1, . . . n/2, l < k βl−βk⟨ˆ b† n/2+kˆ bn/2+l⟩2 k, l = 1, . . . n/2 Table 4: Real and imaginary parts of the complex eigenfrequencies Λr j−iΛi jof the matrix M(n) ufor n-mode bosonic system (n > 2) with unidirectional coupling having a QHP with nth-order ED and second-order DD derived from the equations for the FOMs up to second order. Table is valid for even n, where the vector ˆ bof the diagonalized field operators is written as ˆ b= [ˆ b1,ˆ b2,ˆ b† 1,ˆ b† 2, . . .ˆ bn−1,ˆ bn,ˆ b† n−1,ˆ b† n]. The vector ˆ battains the form ˆ b= [ˆ b1,ˆ b† 1,...,ˆ bm,ˆ b† m,...,ˆ bm+1,ˆ bm+2,ˆ b† m+1,ˆ b† m+2, . . .ˆ bn−1,ˆ bn,ˆ b† n−1,ˆ b† n]if there exist mun-paired eigenvalues in the subsystems that compose the analyzed bosonic system [see, e.g., Eqs. (23), (47), and (54) in Ref. [1]]. In such cases, the columns entitled Moments have to be modified accordingly, but all other columns remain valid. We have ˆ B1≡[ˆ b1,ˆ b† 1,...,ˆ bn/2,ˆ b† n/2],ˆ B2≡[ˆ bn/2+1,ˆ b† n/2+1,...,ˆ bn,ˆ b† n], and more details are given in the caption to Tab. 2. For m > 0un-paired eigenvalues in the subsystems, additional operators ˆ B3,..., ˆ B2m+2 have to be introduced and the table has to be extended. Accepted in Quantum 2025-11-04, click title to verify. Published under CC-BY 4.0. 32