Quantum Liouvillian exceptional and diabolical points for bosonic fields with quadratic Hamiltonians: The Heisenberg-Langevin equation approach

Equivalent approaches to determine eigenfrequencies of the Liouvillians of open quantum systems are discussed using the solution of the Heisenberg-Langevin equations and the corresponding equations for operator moments. A simple damped two-level atom is analyzed to demonstrate the equivalence of both approaches. The suggested method is used to reveal the structure as well as eigenfrequencies of the dynamics matrices of the corresponding equations of motion and their degeneracies for interacting bosonic modes described by general quadratic Hamiltonians. Quantum Liouvillian exceptional and diabolical points and their degeneracies are explicitly discussed for the case of two modes. Quantum hybrid diabolical exceptional points (inherited, genuine, and induced) and hidden exceptional points, which are not recognized directly in amplitude spectra, are observed. The presented approach via the Heisenberg-Langevin equations paves the general way to a detailed analysis of quantum exceptional and diabolical points in infinitely dimensional open quantum systems.


Introduction
Non-Hermitian Hamiltonians, for systems with properly balanced dissipation and amplification, have real energy spectra if they exhibit the parity-time (PT ) symmetry, as shown by Bender and Boettcher [1]. That discovery has lead to the development of non-Hermitian quantum mechanics [2,3,4,5] and has triggered impressive research interest ranging from studying fundamental aspects of quantum physics [6,7,8,9,10,11] to proposing applications in quantum metrology, optics, optomechanics, and photonics [12,13,14]. Studies of non-Hermitian quantum mechanics include also finding quantum analogues of general relativity concepts (like Einstein's quantum elevator [15]), proving various no-go theorems for information processing with non-Hermitian systems [16,17], and even proposals to apply specific non-Hermitian systems, with energy spectra corresponding to the zeros of the Riemann zeta function, aimed at proving the Riemann hypothesis [18]. dissipative and/or amplified systems have been proposed (for reviews see [19,20,14] and references therein).
Most of these studies on EPs are limited to Hamiltonian EPs, which correspond to the degeneracies of the eigenvalues of non-Hermitian Hamiltonians associated with their coalescent eigenvectors. We note that these EPs may be called semiclassical, because they are not affected by quantum jumps, as explicitly discussed in [21] based on the quantum-trajectory theory [22], which is also referred to as the Monte Carlo wave-function method [23,24] or the quantum-jump approach [25]. Indeed, non-Hermitian Hamiltonians describe only a continuous nonunitary dissipation of an open system, but its full quantum description requires also the inclusion of quantum jumps in the system evolution. Recently, these semiclassical Hamiltonian EPs have been generalized to the quantum regime by analyzing the EPs of quantum Liouvillians instead of those of non-Hermitian Hamiltonians [21]. Specifically, quantum EPs (QEPs) can be defined as the degeneracies of the eigenvalues corresponding to coalescent eigenmatrices (eigenoperators) of the quantum Liouvillian superoperator for a Lindblad master equation. Clearly, such an approach to EPs includes the effect of quantum jumps, as it is based on the standard master-equation approach [26] with a trivial metric for arbitrary systems. In contrast to this, a complete description of the evolution of a non-Hermitian system within the Bender-Boettcher quantum mechanics requires also calculating the evolution of a nontrivial metric of the system. This is rather complicated and nonintuitive, but necessary to obtain physical results without violating any no-go theorems [16,17]. The effects of quantum jumps on QEPs have been recently confirmed experimentally in [27,28].
QEPs and Hamiltonian EPs are, in general, different, although they can be equivalent for classical-like systems (e.g., linearly coupled harmonic oscillators [21]) or specific finitedimensional systems [29]. Anyway, one can experimentally observe the transition of QEPs into Hamiltonian EPs by a proper postselection of quantum trajectories within the hybrid-Liouvillian formalism of Ref. [30]. We also note that if the eigenvectors corresponding to degenerate eigenvalues of Hamiltonians or Liouvillians do not coalesce, then such points are refereed to as diabolical points. These points have also wide applications in witnessing quantum effects including phase transitions. For example, quantum Liouvillian diabolical points can reveal dissipative phase transitions and a Liouvillian spectral collapse [31,32].
For the above reasons, the determination of QEPs and their properties represents an important task, which can lead to applications in quantum technologies, including quantum metrology. Unfortunately, the standard approach of finding QEPs via the eigenvalue problem of Liouvillians becomes quite inefficient for multi-qubit or multi-level quantum systems. For systems with infinitely-dimensional Hilbert spaces, the determination of EPs and QEPs is even more challenging. Here, we develop an efficient method based on the Heisenberg-Langevin equations for finding QEPs, and we show the equivalence of QEPs found by these two approaches. We apply the developed method for the determination of QEPs for a system of M bosonic modes mutually interacting via quadratic nonlinear Hamiltonians. Such Hamiltonians are appealing as they lead to the linear exactly-soluble Heisenberg operator equations and, simultaneously, allow to describe nonclassical optical fields including squeezed [33] and sub-Poissonian fields [34]. Such fields are analyzed in the framework of a generalized superposition of a signal and noise [35] that describes quantum Gaussian fields.
The applied method uses the analysis of the dynamical equations for the moments of field operators. It is based upon the equivalence of the eigenfrequency analysis of the Liouvillian and the analysis of an arbitrary complete set of operators of measurable quantities. In the case of quadratic Hamiltonians the field operators moments (FOMs) can be chosen as a suitable set of such measurable operators. These moments form specific structures [36,37,38] that exhibit QEPs of varying (exceptional) degeneracies [QEP degeneracies] and multiplicities [diabolical degeneracies (QDP degeneracies) for quantum hybrid diabolical exceptional points (QHPs)]. Contrary to other more general Hamiltonians, the first-order FOMs of quadratic Hamiltonians form a closed set of dynamical equations. Moreover, the eigenfrequencies obtained from these equations allow for the determination of those for higher-order FOMs. These FOMs then provide a complete structure of eigenfrequencies whose numbers increase with the increasing order of the moments. Among others, this results in the identification of eigenfrequencies with infinite QEP and QDP degeneracies forming QEPs and QHPs. There also occur cases in which degenerate eigenfrequencies remain the same at a QEP, i.e. this QEP is not identified by the spectrum. We may refer to hidden quantum exceptional points (hidden QEPs).
The paper is organized as follows. The equivalence of the system eigenfrequencies analyses based on the Liouvillian and Heisenberg equations for the operators of measurable quantities is discussed in general in Sec. II. Both approaches are compared in Sec. III by considering a simple system of a damped two-level atom. Section IV is devoted to the application of the Heisenberg-Langevin equations for an M -mode bosonic system with a quadratic Hamiltonian and the eigenfrequency analysis of the equations for FOMs. The eigenfrequencies and their general structure are discussed in Sec. V by considering a twomode bosonic system with a general quadratic Hamiltonian. Conclusions are drawn in Sec. VI. The correspondence between the generalized master equation and the Heisenberg-Langevin equations is discussed in Appendix A.

Equivalence of eigenfrequency analyses in the spaces of statistical operators and measurable operators
Before we address the analysis of QEPs for specific models, we show that the eigenfrequency analysis of a LiouvillianL in the space of statistical operators can be equivalently replaced by the eigenfrequency analysis for an arbitrary complete set of measurable operators. The equivalence of projection operator techniques, used for the derivation of generalized master equations with their Liouvillians, and the set of the Heisenberg equations is discussed in detail in [39,40,41]. The Liouvillian superoperatorL· ≡ −i/h[Ĥ, ·], where · stands for an arbitrary operator, [, ] denotes the commutator, andh means the reduced Planck constant, is defined in terms of the overall HamiltonianĤ involving the system and its reservoir in its complete description. It is used to evolve the statistical operatorρ(t) according to the Liouville equation We consider an N -dimensional Hilbert space with an arbitrary basis |j for j = 1, . . . , N . The corresponding basis in the Liouville space of a statistical operatorρ is formed by vectorsρ jk ≡ |j k|. They allow to expressρ in terms of coefficients ρ jk as follows: The Liouville equation in (1) transforms into the following set of equations for the coefficients ρ jk : and L mn jk = m|(L|j k|)|n . The diagonalization of the matrix L mn jk provides us complex eigenvalues −iΩ α that define eigenfrequencies Ω α . The corresponding eigenvectors l α jk form the eigenoperatorsl α = jk l α jk |j k|. It holds for an anti-Hermitian superoperatorL that for any eigenvalue −iΩ α with eigenoperatorl α there also exists the eigenvalue iΩ * α with the corresponding eigenoperatorl †α . Moreover the LiouvilliansL have one eigenfrequency Ω 0 = 0 with the Hermitian eigenoperatorl 0 =l †0 and Tr{l 0 } = 1 that describes the steady state. Contrary tol 0 , the remaining eigenoperatorsl α are non-Hermitian and obey Tr{l α } = 0.
In parallel to Eq. (2), an arbitrary statistical operatorρ can be decomposed in the basiŝ l α as followsρ where ρ α (t) = Tr{ρ(t)l α † }. The evolution of the statistical operatorρ(t) is then expressed along the formulaρ that 'separates' the evolution for different eigenfrequencies Ω α and −Ω * α . The coefficients ρ α (0) in Eq. (5) characterize the initial state. The mean value A (t) of an arbitrary Hermitian operatorÂ(0) of a measurable quantity is then determined in the Schrödinger picture as using Eq. (5) and the coefficients A α (0) = Tr{Â(0)l α † }. Now let us consider an arbitrary Hermitian operatorÂ(t) of a measurable quantity and its evolution in the Heisenberg picture according to the Heisenberg equation The eigenoperators of the superoperator −L coincide with those of L and the corresponding eigenvalues differ by sign. Decomposing the operatorÂ(0) at t = 0 into the basisl α , we may express the solution of the Heisenberg equation (7) as follows: Using Eq. (9), the mean value A (t) of the operatorÂ in the Heisenberg picture is expressed in terms of the coefficient ρ α (0) of the initial statistical operatorρ(0): Provided that the LiouvillianL is constructed using a Hermitian HamiltonianĤ the eigenfrequencies Ω α are real and the formulas for the mean values A (t) in Eqs. (6) and (10) coincide. This means that the time evolution of the operatorÂ(t) is described only by the eigenfrequencies Ω α and −Ω α known from the evolution of the statistical operatorρ fulfilling Eq. (1). Once we construct an arbitrary basis from the Hermitian operatorsÂ of measurable quantities and analyze the evolution of their mean values A (t), we reveal all the eigenfrequencies Ω α . The used basis can even be chosen more generally involving also non-Hermitian operatorsÂ. When deriving a master equation for the reduced statistical operatorρ s , by tracing out the part of the whole statistical operatorρ belonging to the reservoir, we apply the perturbation solution of the general Liouville equation (1) valid up to the second power of the interaction Hamiltonian between the system and its reservoir. This results in a new Liou-villianL s that has a more general form compared to that expressed as a commutator with the HamiltonianĤ. We may diagonalize this more general form of the system Liouvillian L s to reveal, in general, complex eigenvalues −Ω α and Ω * α and the accompanying eigenoperatorsl α andl α † . Alternatively, we convert the system Liouville equation in Eq. (1) into a coupled set of differential equations for the mean values A (t) of operatorsÂ that form a basis in the space of system operators of measurable quantities. According to Eq. (6), the complex eigenfrequencies of the dynamics matrix of this set of differential equations coincide with the eigenfrequencies revealed by a direct diagonalization of the system Liouvillian L s . For details, see Appendix A. The differential equations for the mean values A (t) can alternatively be derived from a closed set of the Heisenberg equations written for both the system operatorsÂ forming the basis and operators of the reservoir. The reservoir operators can suitably be eliminated keeping the validity of the solution up to the second power of the interaction Hamiltonian. This leads to the Heisenberg-Langevin equations. This method known as the Wigner-Weisskopf model of damping [35] for a bosonic mode interacting with bosonic reservoir provides equivalent derivation of the differential equations for the mean values A (t). This stems from the above discussed equivalence of both descriptions in the Schrödinger and Heisenberg pictures, and equivalent approximations when eliminating the reservoir degrees of freedom [42].
The use of the Heisenberg equations for suitable operators instead of the master equation may provide qualitative advantage compared to the eigenfrequency analysis of the Liouvillian. We demonstrate this by analysing two important examples: a damped twolevel atom in Sec. III and the system of mutually interacting bosonic modes with the quadratic interaction Hamiltonian in Secs. IV and V. Whereas the analysis of a damped two-level atom in a finite-dimensional Liouville space leads to a deeper insight into the method based on the Heisenberg equations by its comparison with a direct diagonalization of the Liouvillian, we reveal the structure and values of eigenfrequencies in the model of bosonic modes relying just on the Heisenberg equations.

Eigenfrequency analyses of a damped two-level atom
We demonstrate the above discussed equivalence of both approaches in the determination of the system eigenfrequencies by analyzing, probably, the simplest possible physical system -a damped two-level atom. Its Liouville space of statistical operators has dimension 4 and so its Liouvillian L has 4 eigenfrequencies and eigenoperators. Moreover, as the system is finite dimensional, powers of the operators of measurable quantities do not play an important role in the eigenfrequency analysis because they are just a linear superposition of operators from any chosen basis in the space of measurable operators. The eigenfrequency analysis performed by both approaches also results in their detailed comparison.
The considered two-level atom has the ground state |0 and the excited state |1 . It is described by the Pauli operatorsσ x =σ + +σ − ,σ y = (σ HamiltonianĤ s of the two-level atom with the excitation energyhω is written aŝ

Analysis via the Liouvillian superoperator
We begin with the eigenfrequency analysis of the corresponding Liouvillian. We assume that the atom is damped via the interaction with a reservoir based on theσ x operator. The corresponding LiouvillianL x is derived in the form where the dot · stands for an arbitrary operator. For a more general form of damping and the corresponding analysis of the Liouvillian, see [21]. Expressing the statistical operatorρ s of the two-level atom in a suitable basis, we transform the Liouville equation (1) into the following set of linear differential equations for the coefficients of the decomposition: The dynamics matrix M defined in Eq. (15) has 4 eigenfrequencies: where Ω s = ω 2 − γ 2 x . The eigenfrequencies Ω s 1,2 identify a possible QEP for We may also determine the corresponding eigenvectors that are conveniently written as columns of the transformation matrix P: We note that the first eigenvectorρ s 0 for Ω s 0 that describes the stationary state is normalized such that Tr{ρ s 0 } = 1. All other eigenvectors are traceless (Tr{ρ s j } = 0 for j = 1, 2, 3) and no specific norm is introduced to normalize them. We have Ω s = 0 for the condition in Eq. (17) and so the second and third eigenvectors (columns) in the matrix P in Eq. (18) coalesce confirming the presence of an QEP.

Analysis via the Heisenberg-Langevin equations
Now let us apply the second approach. To derive the appropriate Heisenberg-Langevin equations we need to specify the system interaction with the reservoir represented by a large group of two-level atoms. Their HamiltonianĤ r 0 and a suitable interaction Hamiltonian H r int are expressed as: where ω r j stands for the frequency of the jth reservoir two-level atom that is coupled with the analyzed two-level atom via the real coupling constant κ r j . The Pauli operators of the reservoir two-level atoms are introduced in analogy to those of the analyzed system.
We may write the Heisenberg equations for both analyzed two-level atom and reservoir two-level atoms described by the overall HamiltonianĤ s +Ĥ r int +Ĥ r 0 . A systematic elimination of the reservoir operators then results in the following Heisenberg-Langevin equation for an arbitrary operatorX [42]: In Eq. (20) the Langevin operator forceF represents the back-action of the reservoir two-level atoms to the analyzed atom. The damping constant γ x is derived from the properties of the reservoir Langevin operator forceF x (t) along the formula: Moreover, the reservoir properties imply that F x (t) = 0.
Using the general formula in Eq. (20), we write the Heisenberg-Langevin equations for the four operatorsσ 0 ,σ 1 ,σ + , andσ − that form a basis in the space of operators of the measurable quantities: Applying averaging over the reservoir operators in Eq. (23), we arrive at the following equations for the mean values of the system operators: where the dynamics matrix in Eq. (24) coincides with that in Eq. (14) because, using the representation of the statistical operatorρ in Eq. (13), we have σ 0 = ρ s 00 , σ + = ρ s 01 , σ − = ρ s 10 , and σ 1 = ρ s 11 . Thus, the eigenfrequencies obtained by both approaches are identical.
We may also exploit the transformation matrix P in Eq. (18) to express the evolution of the mean values of the above operators: where c(t) = cos(Ω s t), s(t) = sin(Ω s t), and Ω s = ω 2 − γ 2 x ; sh and ch denote the hyperbolic sine and cosine functions. We can see from the solution in Eq. (25) that there exist two subspaces in the space of the operators of measurable quantities whose dynamics differ. Whereas there occurs only damped dynamics in the subspace spanned by operatorŝ σ 0 andσ 1 , the oscillatory evolution at frequency Ω s in the subspace spanned by vectorsσ + andσ − allows to identify an QEP (Ω s = 0). We may conclude in general that to identify a QEP we need to follow the dynamics of an operator that has a nonzero overlap with the subspace spanned byσ + andσ − .
For a real atom, monitoring the dynamics of the mean values σ + (t) and σ − (t) via the observation of field polarization generated by the atom is a natural choice. It can be accomplished by measuring, e.g., optical nutation or free induction decay [43]. The presence of a QEP can also be confirmed by measuring the amplitude spectra of the field emitted from the atom or, in case of continuous-wave fields, by spectral analysis of field intensity correlation functions [44,45].
We note that the consideration of powers of operators in any basis in finite-dimensional spaces does not bring advantage into the eigenfrequency analysis, because these operators can be expressed as linear combinations of the operators from this basis. This, among others, also implies that only the first-and second-order correlation functions of the reservoir operator forcesF are needed. This contrasts with the behavior of systems defined in the infinitely dimensional Liouville spaces in which moments of operators are the most useful, as we can see below in Sec. IV.

Eigenfrequency analysis of an M -mode bosonic system
We apply the above discussed equivalence of eigenfrequency analyses in the Schrödinger and the Heisenberg pictures to discuss the system of the eigenfrequencies of M mutually interacting bosonic modes described by their annihilation (â j , j = 1, . . . M ) and creation (â † j ) operators. To avoid approximations, when solving nonlinear differential equations, we assume that the M -mode system is described by the general quadratic HamiltonianĤ 0 : where the elements of matrix ( * ij = ji ) describe the linear coupling between pairs of modes, whereas the elements of the matrix κ characterize the nonlinear interactions between pairs of modes (the annihilation and creation of photon pairs). Such Hamiltonian allows to describe both squeezed-light generation and production of entangled states. Symbol H.c. in Eq. (26) replaces the Hermitian conjugated term.
The modes may be either damped or amplified. The system LiouvillianL 0 ≡ −i/h[Ĥ 0 , ·] has to be extended by the termsL d j , provided that the jth mode is damped with damping constant γ d j . On the other hand, the amplification of the mode k is described by the following additional termsL a k , using the amplification constant γ a k . The master equation for the statistical operatorρ(t) of the M -mode bosonic system comprising the LiouvilliansL 0 ,L d j andL a k is equivalently replaced [42,35] by the following system of the Heisenberg-Langevin equations conveniently written in the matrix form: The dynamics matrix M Ω is derived from the HamiltonianĤ 0 in Eq. (26) using the canonical commutation relations. Moreover, for the jth damped mode it includes additional terms −γ d jâ j (t)/2 and −γ d jâ † j (t)/2 on the r.h.s. of equations for dâ j (t)/dt and dâ † j (t)/dt, respectively. The stochastic Langevin operator forcesL j andL † j obey the relations L j (t)L † j (t ) = γ d j δ(t − t ) and L † j (t)L j (t ) = 0 in this case. Similarly, concerning the kth amplified mode the dynamics matrix M Ω contains the terms γ a kâ k (t)/2 and γ a kâ † k (t)/2 on the r.h.s. of equations for dâ k (t)/dt and dâ † k (t)/dt, respectively. The relations L k (t)L † k (t ) = 0 and L † k (t)L k (t ) = γ a k δ(t − t ) are obtained in this case. The Dirac δ functions characterizing the temporal correlations of Langevin forces reflect the Markovian character of the interaction with spectrally broadband reservoirs in individual modes [43]. More details are given in Appendix A.
The frequency analysis of the dynamics matrices of the differential equations for the FOMs allows us to completely determine all the eigenfrequencies of the LiouvillianL and to reveal their structure. This is owing to the linear form of the corresponding Heisenberg-Langevin equations in Eq. (30).
We first transform the original field operatorsâ j andâ † j , grouped in the vector a, to field operatorsb j andb † j forming the vector b where The stochastic Langevin operator forces, which are grouped in the vectorK T = (K 1 ,K † 1 , . . . , K M ,K † M ), induce their Gaussian and Markovian character from the original stochastic operator Langevin forces, as written in the vectorL. The set of equations in Eq. (31) can be solved as follows: Now let us analyze the equations for FOMs step by step by increasing the order of FOMs. The analysis of equations for the first-order FOMs written for either the original operators, or the transformed operators, immediately gives us the basic set of eigenfrequencies contained on the diagonal of the matrix Ω.
Combining Eqs. (31) and (33), we arrive at the following differential equations for the second-order FOMs (for the case without the Langevin forces, see [36]): in which the matrixK contains time-independent second-order correlation functions of the Langevin forces written in vectorK: The solution of the 'diagonal' form of the equations for the second-order FOMs, written in Eq. (36), is obtained in the form similar to that of Eq. (33). This reveals that the emerging eigenfrequencies Ω kk + Ω ll can be expressed in terms of the sums of those from the basic set. The structure of differential equations for the third-order FOMs is more general, as it contains also the time-dependent first-order FOMs on their r.h.s.: The solution to the set of equations in Eq. (39) can again be expressed in the form of Eq. (33). Whereas the homogeneous part of the solution to Eq. (39) oscillates at frequencies Ω kk + Ω ll + Ω mm , the nonhomogeneous part of the solution contains additional frequencies Ω nn , as it can be checked by a direct calculation. The general structure of the differential equations for the fourth-and higher-order FOMs is the same as that for the third-order FOMs. To demonstrate this, we derive the following equations for the fourth-order FOMs: The solution to Eq. (41) in the form of Eq. (33) reveals the two types of eigenfrequencies: Ω kk + Ω ll + Ω mm + Ω nn and Ω kk + Ω ll . In general, the analysis of differential equations for pth-order FOMs (for p > 3) reveals the sums of p (for the homogeneous solution) and p − 2 (for the nonhomogeneous solution) eigenfrequencies from the basic set.
The analysis of differential equations for the FOMs of increasing order gradually reveals all the eigenfrequencies and their degeneracies. At a general level, we may determine the number n (p) Ω of different eigenfrequencies provided by the analysis of the homogeneous solution of the set of equations for pth-order FOMs. This number is given by the number of independent pth-order moments. The overall number of pth-order moments equals (2M ) p . However, the mean values of the products of commuting field operators are insensitive to their ordering [e.g., â 1â2 = â 2â1 for the field operatorsâ 1 andâ 2 ]. Also, if two field operators do not commute, the mean values of the products of p operators with different ordering of noncommuting operators differ by the mean value of the product of p − 2 operators [e.g., xââ †ŷ = xâ †âŷ + xŷ for the field operatorsâ andâ † and arbitrary operatorsx andŷ]. Such lower-order FOMs form the r.h.s. of the differential equations for the pth-order FOMs, together with the terms arising from the interaction with the reservoirs.
The numbers of independent FOMs and, thus, the numbers n (p) Ω of eigenfrequencies can be easily determined using the combinations of numbers. Considering the moments up to the fourth order, we arrive at the following formulas: where C(k, l) = k!/[(k − l)!l!] is the binomial coefficient. The eigenfrequencies arising from the equations for the FOMs of different orders form specific structures. Their analysis then allows to identify QEPs and QHPs and their degeneracies. In the following section, we study this structure for a system composed of M = 2 modes and exhibiting both damping and amplification.
The eigenvalue analysis of the dynamics matrices of different FOMs provides also the corresponding eigenvectors. Whereas the obtained eigenvalues coincide with those revealed by the eigenvalue analysis of the corresponding Liouvillian L, the obtained eigenvectors do not allow to construct the eigenvectors of the Liouvillian L. The reason is that the appropriate equations are of a different kind. Whereas the influence of the reservoir noise is involved in the set of equations for second-order FOMs, in Eq. (36), via its nonhomogeneous solution, it is directly embedded in the form of the Liouvillian L whose eigenvalue problem is solved. Also, the presence of reservoir noise leads to the coupling of equations for the odd [and similarly even] orders of FOMs, as shown in Eq. (39) [Eq. (41)]. However and most importantly, this coupling is specific and, as discussed above, keeps the eigenvalues obtained for the homogenenous solutions unchanged.  Table 1: Real and imaginary parts of the complex eigenfrequencies Ω r gen − iΩ i gen derived from the equations for the FOMs up to fourth order which reveal QEPs and QHPs for g = 0. The corresponding moments in the 'diagonalized' field operators are written together with their degeneracy coming from different orderings of field operators. QDP degeneracy of QHPs (partial QDP degeneracy) derived from the indicated FOMs and QEP degeneracy of the constituting QEPs are given.

Spectral eigenfrequencies of a two-mode system with damping and amplification
In this section, we analyze the system of two interacting modes: one exhibiting damping (γ d 1 > 0), the other being amplified (γ a 2 > 0). We consider both linear exchange of photons between the modes (real ), as well as emission and annihilation of photon pairs in these modes. Photon pairs can be annihilated and created either inside both modes (real g) or their photons belong to different modes (real κ). The corresponding HamiltonianĤ 0 is written as follows: which leads to the following Heisenberg-Langevin equations: The only nonzero second-order correlation functions of the stochastic Langevin operator The diagonalization of the dynamics matrixM in Eq. (46) reveals four eigenfrequencies from the basic set: , and γ + = (γ d 1 + γ a 2 )/4. According to Eq. (48), all the imaginary parts of four eigenfrequencies are equal. When damping in mode 1 is stronger (weaker) than amplification in mode 2 [γ d 1 > γ a 2 (γ d 1 < γ a 2 )], the overall system is damped (amplified).
We note that the constant g is assumed real without the loss of generality: It corresponds to a suitable choice of the phases of the field operatorsâ 1 andâ 2 in Eq. (44). On the other hand, the phases of possibly complex constants and κ have a good physical meaning and influence the system dynamics to a certain extent. The derivation of eigenfrequencies Ω in this most general case results in the formulas that are derived from those in Eq. (48) by the formal replacement | | → and |κ| → κ. This means that the spectrum of eigenfrequencies with its QEPs, QDPs, and QHPs discussed below, remains the same. However, the formulas for eigenvectors in Eq. (49), given below, have to be replaced by more general ones in this case. We also note that no QEP can be observed when different frequencies of the modes are considered. The unnormalized eigenvectors Y 1 ,Ȳ 1 , Y 2 , andȲ 2 belonging in turn to the eigenfrequencies written in Eq. (48) are derived as follows (assuming g ≥ 0): Assuming g = 0 and provided that the condition or the condition 2 − κ 2 + γ 2 for the system parameters is fulfilled, two of the eigenfrequencies coincide. This identifies a QEP for which the eigenvectors corresponding to both eigenfrequencies coalesce. Each of the above conditions forms a hypersurface of dimension 4 in the space of independent parameters (γ d 1 , γ a 2 , , κ, g). Replacing the parameters γ d 1 and γ a 2 by γ + and considering linearity of the dynamics equations, the positions of QEPs form two doubled concentric cones in the 3-dimensional space (γ + / , κ/ , g/ ) plotted in Fig. 1.
If g = 0, two-fold degeneracy in eigenfrequencies occurs. Real nonzero eigenfrequencies Ω r 1,2 exist only for 2 − κ 2 − γ 2 + > 0. When all the four eigenfrequencies coincide. They form doubly degenerate QEPs localized at a hypersurface of dimension 3 in the parameter space (γ d 1 , γ a 2 , , κ) defined by the condition in Eq. (52). The positions of the QEPs fulfilling Eq. (52) form the circle with radius 1 in the space of parameters γ + / and κ/ shown in Fig. 1 (the intersection of the yellow and blue cones). At these QEPs, there exist two different eigenvectors each arising from the original collapsing two dimensional spaces [46]. We have the QHPs in this case. We note that the nonclassical properties of optical fields generated at and around these QHPs were analyzed in [46] where even a more general Hamiltonian involving the additional Kerr nonlinear terms was considered. Now let us have a deeper look at the eigenfrequencies and their spectral bifurcations that identify QEPs. To correctly identify QEPs, we also need to know the eigenvectors that correspond to the analyzed eigenfrequencies. The eigenvectors Y 1 ,Ȳ 1 , Y 2 , andȲ 2 given in Eq. (49), arising in the diagonalization of the dynamics matrix of the Heisenberg-Langevin equations (45), and belonging in turn to the eigenfrequencies Ω 1 , −Ω * 1 , Ω 2 , and −Ω * 2 , may be used to form the eigenvectors of the dynamics matrices of FOMs with increasing order. They directly represent the eigenvectors of the first-order FOMs dynamics matrix and, when formed into the supervector Y T ≡ (Y 1 ,Ȳ 1 , Y 2 ,Ȳ 2 ), they allow to express the eigenvectors of the dynamics matrix of a pth-order FOMs via the tensor product Y ⊗ . . . ⊗ Y p [38]. This allows us, among other properties, to identify the number of eigenfrequencies for a given order of FOMs and their degeneracies occurring at QEPs.
The spectra of the eigenfrequencies and the numbers of eigenvectors differ in the abovediscussed two cases (g = 0 and g = 0). Whereas three independent eigenvectors of the dynamics matrix of the Heisenberg-Langevin equations occur at the QEPs for g = 0, only two of the eigenvectors are found at the QHPs when g = 0. We note that, in both cases, all the eigenvalues contribute to the dynamics of the original field operatorsâ 1 ,â † 1 ,â 2 , andâ † 2 and so the analysis of the evolution of any of them allows, in principle, to identify QEPs. In the following eigenfrequency analysis, we pay a detailed attention to the eigenfrequencies belonging to the FOMs up to the fourth order and draw some conclusions concerning general orders.

Spectra of eigenfrequencies for a single nondegenerate QEP
We first consider the case for g = 0, in which a single QEP with a double QEP degeneration occurs in the spectrum of the dynamics matrix of the Heisenberg-Langevin equations. At this QEP, three independent eigenvectors suffice in describing the system evolution. In general, the eigenfrequency analysis of the equations for the first-order FOMs provides four eigenfrequencies ±Ω r 1,2 − iΩ i . Two of them (say ±Ω r 1 − iΩ i ) reveal the position of a QEP identified by the condition Ω r 1 = 0. This position of the QEP is also indicated by the eigenfrequencies originating in the analysis of higher-order FOMs. We have in turn 4, 16, 64, and 256 moments of the first, second, third, and fourth orders. However, as discussed above, some of these moments differ just in the positions of field operators. This means that they are either equal or differ by the moments of lower orders if the involved field operators do not commute. Taking into account this moment degeneracy, we may expect at maximum 4, 10, 20, and 35 different eigenfrequencies from the analysis of moments of the first, second, third, and fourth order. Degeneracy of the moments is given in Tables 1 and 2 for this case. This degeneracy is either mapped into the multiplicity of the corresponding QEPs (forming QHPs from QEPs and being characterized by a QDP degeneracy) or results in higher QEP degeneracies of QEPs (higher-order QEPs [36]). Whereas the moments b 1b2 , b 2b1 , b † 1b 2 , and b 2b † 1 serve as examples in the former case, the moments b † 1b 1 and b 1b † 1 participate in forming a four-fold degenerated QEP (see Table 1). The eigenfrequencies attained by the analysis of the first-, second-, third-and fourth-order FOMs are summarized in Tables 1 and 2 depending on their ability to form QEPs. We note that some of the eigenfrequencies summarized in Table 2, that do not form QEPs, are degenerated, i.e. they exhibit QDP degeneracy.
Provided that Ω i = 0, we observe in the spectra of eigenfrequencies in turn 1, 3, 6, and 10 bifurcations coming from the behavior of the first-, second-, third-, and fourth-order FOMs, as listed in Table 1. Following Table 1, there occur QHPs with QDP degeneracies 2, 6, and 12 considering in turn the moments of the second, third, and fourth order. On the other hand, the maximal QEP degeneracy of QEPs reached for the first-, second-, third-, and fourth-order FOMs equals 2, 4, 8 and 16. Here, we conclude that, in general, the analysis of pth-order FOMs gives QEPs with a 2 p -fold degeneracy. We note that some eigenfrequencies remain single as seen in Table 2. In Tables 1 and 2, we mention, side by side with the eigenfrequencies, the corresponding moments of the 'diagonalized' field operators, as there is one-to-one mapping between these moments and the structure of the eigenfrequencies.
We note that, at the QEPs observed in the dynamics matrix of the Heisenberg-Langevin equations, the number of independent eigenvectors, given in Eq. (49), decreases from 4 to 3. This is so as two eigenvectors have to coalesce at a nondegenerate QEP. This results in the reduction of the complexity of the system dynamics and leads to all physical effects discussed in relation to QEPs. Considering higher-order FOMs, the number of independent eigenvectors arising from the dynamics matrix of pth-order FOMs decreases from 4 p to 3 p , which also gives the maximal number of possibly different eigenfrequencies. Thus, the dynamics of higher-order field-operator correlation functions is more simplified at QEPs than that of the field mean operator amplitudes.
The eigenfrequencies that are related to the moments containing the 'building block' b † 1b 1 (e.g., the moments b † 1b 1 and b † 1b 1b2 in Table 1) form hidden QEPs. The contribution to the overall eigenfrequency of a higher-order FOM from this 'building block' equals zero as Ω r 1 − Ω r 1 = 0 independently of whether there is a QEP or not. However, the presence of a QEP is identified by the reduction of the number of eigenvectors by one that happens in this case without any manifestation in the eigenfrequencies spectrum. Thus, no spectral bifurcation commonly used for identifying QEPs is observed. We note that, for the QEPs listed in Table 1, such eigenfrequencies form QEPs more than doubly degenerated, together with other eigenfrequencies, and so these QEPs are still identified in the spectrum by bifurcations.
In the usually discussed PT -symmetric systems, gain and loss are in balance giving Ω i = 0. In this case, multiplicties (i.e., QDP degeneracies) of QEPs may even be higher as the imaginary parts of eigenfrequencies coming from different orders of FOMs coincide. For example, the QEPs positioned at ±Ω r 1 arise in the eigenfrequencies of both the firstand third-order FOMs, in the latter case even with QDP degeneracy 6 (see Table 1). Such QEPs are called genuine QEPs, as discussed below.
We note that, when describing the evolution of FOMs of a given order, we may neglect the redundant moments and keep just the differential equations for the remaining ones. The eigenfrequencies, originating in the analysis of such equations, remain the same as those analyzed above and summarized in Tables 1 and 2, but their multiplicities are just 1. Also, the reduction in the number of involved moments results in the change of the structure of the space of eigenvectors. This reduction conceals the QHPs identified in the last column of Table 1. We may call such QHPs as the induced QHPs as they originate in the extended space of FOMs that includes also the redundant moments. However, some QHPs, which we refer to as the genuine QHPs, still remain. These occur for Ω i = 0 and are formed by identical eigenfrequencies with different eigenvectors arising in the analysis of dynamics matrices for different orders FOMs (e.g., b 1 and b † 1 versus b 1b2b †  We also note that when reducing the number of necessary FOMs of given order, we face the problem of non-commuting operators. However, the FOMs that contain noncommuting operators at different positions mutually differ by FOMs of orders lower by 2, 4, 6, . . .: Each application of nontrivial commutation relation reduced the moments order by two. When writing the differential equations for the set of FOMs of given order, we arrive at the equations similar to those found in Eqs. (36), (39) and (41) that, however, contain additional terms formed by lower-order FOMs at their r.h.s. Nevertheless, these terms modify the nonhomogeneous solution of the equations qualitatively in the same way as those arising in the fluctuating Langevin forces, i.e. the solution is enriched by the terms oscillating at the eigenfrequencies appropriate to these lower-order FOMs. Thus the terms arising from the non-commuting operators do not change the eigenfrequencies and we may apply the above results concerning the eigenfrequencies also in this case. On the other hand, the eigenvectors in both approaches naturally differ. Their mutual relations [47] were discussed in relation to the quantum versus classical descriptions in [38].

Spectra of eigenfrequencies for a doubly degenerated QEP -QHP
The squeezing-effect part of the HamiltonianĤ 0 in Eq. (44) is often not considered (g = 0). This leads to two doubly degenerated eigenfrequencies ±Ω r − iΩ i , when the dynamics of the first-order FOMs is investigated. They form a QHP occurring directly in the dynamics matrix of the Heisenberg-Langevin equations. This QDP degeneracy is then directly transformed into the diabolical degeneracies of eigenfrequencies arising from the analysis of higher-order FOMs. We may call such QHP the inherited QHP. This inherited QDP degeneracy considerably reduces the number of different spectral eigenfrequencies provided by higher-order FOMs, at the expense of their increasing degeneracies. We note that this behavior was observed in [21], where the spectrum of the Liouvillian of a simplified two-mode bosonic system with g = κ = 0 was numerically analyzed. The obtained eigenfrequencies, their QDP and QEP degeneracies and the corresponding QHPs are summarized in Table 3, together with the corresponding 'diagonalized' FOMs. According to Table 3, the QEPs found for the dynamics matrix of pth-order FOMs exhibit a 2 p -fold QEP degeneracy. Due to the inherited spectral degeneracy, the frequencies of all QEPs belonging to pth-order FOMs equal zero (Ω r = 0), which results in the occurrence of the QHP with a QDP degeneracy 2 p formed by 2 p QEPs with a 2 p -fold QEP degeneracy, as shown in Table 3. This results from the fact that the number of independent eigenvectors of the dynamics matrix of the Heisenberg-Langevin equations decreases from 4 to 2 at QHPs.
As seen in Table 3, we have the hidden QEPs/QHPs also in this case. They occur not only in relation to the moments of a single mode (e.g., b † 1b 1 and b 1b † 1 ), but also when cross moments of different modes are considered (e.g., b † 1b 2 and b 1b † 2 , or b 2b † 1 and b † 2b 1 ). Apart from the hidden QEPs, we observe spectral bifurcations at ±pΩ r , ±(p − 2)Ω r , ±(p − 4)Ω r , . . . from the pth-order FOMs.
Provided that we exclude the redundant FOMs from the description, the QDP degeneracy of QHPs identified in Table 3, in general, lowers. Whereas it remains 2 for the first-order FOMs, it decreases from 2 p to p + 1 for pth-order FOMs.
Finally, we discuss some properties of eigenfrequencies when the FOMs of all orders are considered. Provided that the gain and loss are in balance, we have Ω i = 0 and we can directly compare the eigenfrequencies arising in the dynamics matrices of FOMs of different orders. Then the QEPs are localized with the help of pairs of eigenfrequencies ±pΩ r 1 (g = 0) and ±pΩ r (g = 0) for p = 1, 2, . . . occurring infinitely-many times in the spectrum. For g = 0, the other pairs of eigenfrequencies, as explicitly written in Table 1, are, among others, also found in the spectrum with infinite degeneration.
At the end, we note that, when the eigenfrequency analysis is accompanied by the determination of the corresponding eigenvectors, we can discuss the modes behavior at a general level using the quantities based on the determined FOMs. For example, phase squeezing is revealed by the behavior of the second-order FOMs [33,48], whereas sub-Poissonian photon-number statistics [34] of the modes and their sub-shot-noise photonnumber correlations [49] are quantified by the fourth-order FOMs. Also different types of nonclassicalities can be discussed [50].
We also note that our approach relies on the linear Heisenberg-Langevin equations. Nevertheless, it may be successfully applied also in investigations of quantum systems described by the nonlinear Heisenberg-Langevin equations provided that a suitable operator linearization of the nonlinear operator equations is applied, e.g. around a stationary state or a classical time-dependent solution [46,51,52,53]. Table 3: Real and imaginary parts of the complex eigenfrequencies Ω r gen − iΩ i gen derived from the equations for the FOMs up to fourth order which indicate a QHP for g = 0. The corresponding moments in the 'diagonalized' field operators are written together with their degeneracies coming from different orderings of field operators. QDP degeneracies of QHPs (partial QDP degeneracies) derived from the indicated FOMs and QEP degeneracies of the constituting QEPs are given.

Conclusions
We have shown that the eigenfrequency analysis of the Liouvillians of open quantum systems can alternatively be performed in the space of operators of measurable quantities provided that they form a complete basis. This is especially important for systems defined in infinite-dimensional Liouville spaces including those formed by the interacting bosonic fields. Considering a damped two-level atom, we have demonstrated the equivalence of both approaches in obtaining the system eigenfrequencies and positions of quantum exceptional points. Analysing the dynamical equations of field-operator moments of general M -mode fields, we have revealed the structure of eigenfrequencies attainable from dynamical equations for a given order of the field-operator moments. This shed light to the occurrence of quantum exceptional points identified from the obtained eigenfrequencies: All quantum exceptional points are recognized already from the eigenfrequencies obtained from the first-order field-operator moments. The eigenfrequencies obtained from higher-order field-operator moments are important in revealing multiple (i.e., diabolical) degeneracies of these quantum exceptional points only. We have developed a general approach to analyze a two-mode bosonic system described by a general quadratic Hamiltonian. In its general configuration, two distinct sets of quantum exceptional points occur for nonzero mode squeezing that, however, collapse into a single set with quantum hybrid diabolical exceptional points, when mode squeezing is not considered. In the analysis, we have observed the inherited, genuine and induced quantum hybrid diabolical exceptional points. Moreover, the hidden quantum exceptional points, whose presence is not directly inferred from the behavior of eigenfrequency spectra, were identified.
The consideration of the Heisenberg-Langevin equations for the operators of measurable quantities and the derived dynamical equations for field-operator moments represent a convenient starting point for the system eigenfrequency analysis that allows to reveal the eigenfrequencies of open quantum systems encoded in their Liouvillians. We believe that this approach is qualitatively less demanding compared to a direct diagonalization of the Liouvillians, at least when the linear Heisenberg-Langevin equations describe the analyzed system. This approach paves the way to a general and detailed analysis of quantum exceptional points in open quantum infinitely-dimensional systems.