Realizing exceptional points of any order in the presence of symmetry

Exceptional points~(EPs) appear as degeneracies in the spectrum of non-Hermitian matrices at which the eigenvectors coalesce. In general, an EP of order $n$ may find room to emerge if $2(n-1)$ real constraints are imposed. Our results show that these constraints can be expressed in terms of the determinant and traces of the non-Hermitian matrix. Our findings further reveal that the total number of constraints may reduce in the presence of unitary and antiunitary symmetries. Additionally, we draw generic conclusions for the low-energy dispersion of the EPs. Based on our calculations, we show that in odd dimensions the presence of sublattice or pseudo-chiral symmetry enforces $n$th order EPs to disperse with the $(n-1)$th root. For two-, three- and four-band systems, we explicitly present the constraints needed for the occurrence of EPs in terms of system parameters and classify EPs based on their low-energy dispersion relations.


I. INTRODUCTION
The appearance of symmetry-protected degeneracies in the energy dispersion of various Hermitian topological systems has attracted much attention in the past decades [1][2][3][4][5][6][7][8]. These Hermitian topological systems, aside from their space-group symmetry, are classified using ten symmetry classes [9] identified based on three discrete symmetries, namely time-reversal symmetry, particle-hole (or charge conjugation) symmetry, and chiral (or sublattice) symmetry [7]. Topological semimetals [10][11][12] and multifold fermions [13][14][15] are excellent representatives of such systems in which two-or multiband crossings can be observed in the energy spectra. In the absence of symmetry, these band touchings are generally unstable in lower-dimensional models due to the hybridization of the bands resulting in the gaping out of degeneracies. However, this band repulsion mechanism is absent in topological systems in which crystalline symmetries and/or discrete symmetries, e.g., time-reversal symmetry, may protect band touching points [16].
It has further been shown that the commonly observed linear energy dispersion close to nontrivial degeneracies might be forbidden due to certain symmetry constraints present in some systems [17]. As a result, higher-order band dispersions, such as cubic or quadratic, may find room to arise close to band touching manifolds [18][19][20]. These distinct characters of energy spectra are considered as an additional tool to classify various nontrivial degeneracies in Hermitian systems [16].
The recent surge of theoretical and experimental interests in the field of non-Hermitian systems has advanced our understanding of the intrinsic properties of systems with no Hermitian counterparts. Some of these exotic properties are i) the piling up of bulk states on the boundaries known as the non-Hermitian skin effect [21], which goes hand in hand with a violation of the conventional (Hermitian) bulk-boundary correspondence [22][23][24], ii) the emergence of exceptional points (EPs) [25,26] as defective degeneracies at which the geometric multiplicity is smaller than the algebraic multiplicity, and iii) the observation of different non-Hermitian topological systems due to the closure of non-Hermitian (line or point) gaps [27][28][29].
The emergence of these unique properties of non-Hermitian systems is linked to the extended 38 symmetry classes [25], which are the non-Hermitian counterparts of the tenfold Altland-Zirnbauer classification [7,9] in Hermitian systems. As Hermiticity is not respected in non-Hermitian systems, particle-hole symmetry (PHS and PHS † ) and time-reversal symmetry (TRS and TRS † ) acquire two different flavors, and chiral symmetry (CS) is discerned from sublattice symmetry (SLS). These six symmetries combined with psuedo-Hermiticity (psH) [30] give rise to 38 symmetry classes as defined in Ref. [25]. Aside from these seven symmetries, pseudo-chiral symmetry (psCS) [31], inversion (I) symmetry [32], parity (P) symmetry and its combination with time-reversal (PT ) [33] and particle-hole (CP) symmetries [34] have been considered in exploring various properties of non-Hermitian systems. We summarize these symmetries in Table I. We note that these symmetries can also be written in terms of the classification of random matrices [35] introduced by Bernard and LeClair [36], see also Appendix A for details.
These investigations, which usually explore case stud-arXiv:2202.07009v1 [quant-ph] 14 Feb 2022 ies, mainly address a number of questions as follows: i) How many constraints need to be satisfied to find an nth order EP, dubbed as EPn? It has been argued that 2(n − 1) real constraints should be imposed to detect EPns in systems with no symmetry [26]. Even though a description for these constraints is discussed in Ref. 34, a generic recipe to generate and understand these constraints in the presence of any symmetry is absent in the literature. Nevertheless, it has been suggested that relating each of these constraints to a momentum coordinate implies that merely EP2s can be realized in three spatial dimensions [53,54]. ii) What role is played by symmetries in the appearance of EPns? Recent researches reveal that including symmetries may reduce the number of constraints to realize EPs. As a result, various case studies reported the occurrence of Ep2s [38,39], Ep3s [52], and Ep4s [41] in one, two, and three spatial dimensions, respectively. More extended studies have also explored the link between observing EPns and the presence of either PT [40] or antiunitary [34] symmetries. iii) Is it possible to distinguish EPns based on the low-energy dispersions close to them? Similar to Hermitian systems at which linear, cubic, and quadratic dispersions were reported close to nontrivial degeneracies, nth-root dispersion in the vicinity of EPns were numerously identified [39]. A recent study reports the square-root behavior of band spectra close to EP3s in the presence of SLS [52], where this possibility has also been studied in Ref. 55 without reference to symmetry.
In this work, we revisit these questions using a generic mathematical formulation to explore the appearance of EPns in Hamiltonians represented by n-dimensional matrices. Based on our formalism, we are able to count the number of constraints in the presence of any symmetry and evaluate each constraint based on the traces and the determinant of the Hamiltonian of our interests. In particular, we find that one needs to satisfy tr[H k ] = 0 with k = 2, . . . , n − 1 and det[H] = 0 to find an EP with order n arriving at a total of 2(n − 1) constraints in agreement with the literature. Imposing symmetry considerations, we show that in the presence of CS, psCS, SLS, psH, PT , and CP symmetries, some traces or the determinant of n-band systems generally disappear. Moreover, when psH, PT or CP symmetry is present, we find that the number of constraints is reduced to half, i.e., n − 1 constraints. When we instead consider psCS or SLS, we recover n constraints for n ∈ even and n − 1 constraints when n ∈ odd. CS is only defined in even dimensions, in which case we find n − 1 constraints. We summarize these results in Table II. We, furthermore, identify conditions to characterize various EPns based on their low-energy dispersions. To do so, we introduce an alternative approach based on the Frobenius companion matrix of the characteristic polynomial, which can be interpreted as representing a perturbation close to an EPn. With this matrix in mind, we rederive the above statement pertaining to the 2(n − 1) constraints as well as explicitly calculate the low-energy band dispersions around an EPn. Despite the common assumption that EPns disperse with the nth root, we find that in the presence of SLS or psCS with n ∈ odd, the leading order term of the dispersion around an EPn generically scales with the (n − 1)th root.
We emphasize that our formulation is not limited to any specific spatial dimension. For completeness purposes, we calculate explicit forms for the nonzero constraints for all twelve symmetries listed in Table I for two-, three-, and four-band systems and present their nonzero parameters.
The outline of this paper is as follows. In Sec. II, we present our generic mathematical formulation to describe EPns. We further draw generic symmetry-based arguments on the behavior of EPns when a specific symmetry is respected. Using the generic decomposition of two-band systems in terms of Pauli matrices, we discuss the properties of EP2s, explicit forms of constraints, and collections of nonzero parameters in the presence of each twelve symmetries in Sec. III. In Sections IV and V we pursue similar lines of thought for EP3s and EP4s, respectively. Using the Gell-Mann matrices and their generalization, we rewrite three-and four-band Hamiltonians and identify their nonzero components when a symmetry constraint is enforced. We also discuss various possibilities to observe different energy dispersions close to EP3s and EP4s in Sections IV and V, respectively. We conclude our paper in Sec. VI.

II. EPS IN N-BAND SYSTEMS
Given a generic n×n matrix H, the characteristic polynomial is defined by where and other σ k 's are the sum of kth order diagonal minors of H. Defining p k = (−1) k σ k and s k = tr[H k ], we have  ]] = 0 with k = 2, . . . , n − 1. We emphasize that EPns occur when all of these 2(n − 1) constraints are simultaneously enforced. In parameter regimes in which smaller number of constraints are satisfied, lower-order EPs, e.g., EPms with m ≤ n, may find room to emerge in the spectrum of n × n dimensional matrices.
An alternative approach to counting the number of constraints in matrices with EPns is based on perturbing H close to EPns [60]. Here, we introduce a Jordan block J n as a description for the EPns in H with dimension n and, without loss of generality, diagonal value λ n = 0. Introducing the perturbation matrix δS, one can find all insignificant, trivial perturbations using [δS, J n ] [60]. The remaining non-trivial perturbation is a n × n matrix δJ n . The matrix elements of δJ n are always zero except for n − 1 complex-valued elements, which are δJ n,j with j = 1, . . . n − 1. The summation of the Jordan block and its nontrivial perturbation, namely J n + δJ, describes the low-energy behavior of H close to EPns, which reads Note that when tr[H] is nonzero, the (n, n) matrix element of δJ n is also nonzero. The matrix elements δJ n,j are related to the coefficients σ k since the characteristic polynomial of J n + δJ is identical to Eq. (1).
In fact, H 0 constructs the (transpose) Frobenius companion matrix for the characteristic polynomial in Eq. (1) [61], and each of the δJ n,j is proportional to σ n+1−j , in particular, δJ n,j = (−1) n+j σ n+1−j . This result was also derived in Ref. 62, and further generalized to describe perturbations of any matrix written in the Jordan normal form. From this approach, we again realize that 2(n − 1) constraints are needed to determine the presence of EPns in matrix H, i.e., Re[δJ n,j ] = 0 and Im[δJ n,j ] = 0 with j = 1, . . . , n − 1.
From the characteristic polynomial in Eq. (1) as well as the perturbed Jordan block in Eq. (4), we can also deduce how the EPns disperse. While it is commonly assumed that the series expansion resulting from a perturbation with ω around an EPn, i.e., writing H 0 = J n + ωδJ n , results in the Puiseux series, λ = λ 0 + ∞ j=1 ω j/n λ j with λ 1 = δJ 1/n n,1 , this is not generally the case. Indeed, only when δJ n,1 = 0, the Puiseux series is recovered for the energy eigenvalues close to an EPn [55,63]. When δJ n,1 = 0, the perturbed eigenvalues generally split in different cycles of the form λ = λ 0 + ∞ j=1 ω j/p λ j with p < n, and the different values of p summing up to n [55,63,64].
Let us now see how this translates into our perturbed Jordan block in Eq. (4). In particular, when σ n = 0 and all other σ j = 0 (or equivalently, when δJ n,1 = 0 and all other δJ n,j = 0), we straightforwardly find that the characteristic polynomial reduces to λ n + (−1) n σ n = 0 (or λ n − δJ n,1 = 0). In this case, the EPn disperses with e 2π i r/n [(−1) n σ n ] 1/n (e 2π i r/n [δJ n,1 ] 1/n ) for r = 1, . . . , n [65,66]. When σ n−1 = 0 and all other σ j = 0, (or equivalently, when δJ n,2 = 0 and all other δJ n,j = 0,) we find λ(λ n−1 + (−1) n−1 σ n−1 ) = 0 (or λ(λ n−1 − δJ n,2 ) = 0) for the characteristic polynomial. Now, the EPn disperses with the n−1th root, i.e., ∼ [(−1) n−1 σ n−1 ] 1/(n−1) (∼ [δJ n,2 ] 1/(n−1) ), combined with a flat band with λ = 0. In general, we thus find that when σ j = 0 (δJ n,j = 0) and all other σ k = 0 (δJ n,k = 0), the low-energy approximation around the EPn reads ). When all σ k = 0 (or all δJ n,k = 0), it is no longer possible to find complete analytical solutions for the eigenvalues λ when n ≥ 5. Nevertheless, one can numerically compute explicit solutions for the leading terms [64,67]. So far, we discussed EPns in systems with no additional symmetries. Let us now see how the presence of symmetries affects the appearance of EPns. Writing the determinant and traces as det[H] = i i and tr[H k ] = i k i with i the eigenvalues of H allows us to make general statements when making use of the energy constraints listed in Table I. We immediately see that PHS, PHS † , TRS, TRS † , I, and P symmetry are nonlocal in parameter space as they relate eigenvalues with momentum k to eigenvalues with momentum −k. As such, the presence of these symmetries does not reduce the number of constraints but instead puts a constraint on whether the entries in the Hamiltonian are symmetric or antisymmetric. We thus find that the number of constraints for finding EPns in the presence of these symmetries remains at 2(n − 1). For the remaining symmetries listed in Table I, however, there is a reduction in the number of constraints.
In the presence of SLS and psCS, { (k)} = {− (k)} dictates that in the case of n ∈ odd, at least one of the eigenvalues is necessarily zero, such that det[H] = 0. For n ∈ even, there is no such argument and we thus generally find det[H] = 0. Turning to the traces, we see that tr[H k ] = 0 when k ∈ even, while tr[H k ] = 0 when k ∈ odd for all n. To find EPns, we thus need to satisfy n constraints when n ∈ even, and n − 1 constraints when n ∈ odd. The fact that { (k)} = {− (k)} also leads to an interesting consequence when considering the possibility of realizing lower-order EPs in n-band systems, namely, the addition of an extra band to an (n − 1)-band system immediately promotes a possibly existing EP(n − 1) to an EPn as long as this additional band is coupled to the other bands. As such, there is a notion of fragility in these systems, as also pointed out in Ref. 52 for the case of SLS. However, if a band is added that does not couple to any of the other bands, the EP(n − 1) survives even though the energy eigenvalues are n-fold degenerate at the EP.
If we instead consider PT and psH symmetry, we see that This means that we need to satisfy n − 1 constraints to find an EPn in agreement with what is found in Ref. 34. Lastly, considering CS and CP symme- This gives us again n − 1 constraints, which was also found in Ref. 34, see also Ref. 68. We summarize the results for SLS, psCS, CS, PT , psH and CP symmetry in Table II and refer to Appendix C for details on the derivation of these findings. Now turning back to our results for the dispersion around EPns, we see that in the case of SLS and psCS with n ∈ odd, where det[H] = 0 (i.e., σ n = 0 or δJ 1,n = 0), EPns disperse with O(ω 1/(n−1) ). Interestingly, this is the only instance of symmetries generically preventing the recovery of the nth root dispersion for EPns.
In the following, we explore EPns and the implications of symmetry in greater detail by deriving exact results. Galois theory [69] implies that characteristic polynomials with dimensions greater than four cannot be expressed as combinations of radicals of rational functions of the polynomial coefficients. Therefore, to present analytical results in terms of radicals, we explore the role of symmetries in modifying the structure and numbers of constraints to detect EPns with n = 2, 3, 4.

III. EPS IN TWO-BAND SYSTEMS
To study second order EPs, we perform a matrix decomposition in the Pauli basis. The most generic twoband Hamiltonian in this representation is given by where σ = (σ x , σ y , σ z ) is the vector of Pauli matrices (see Appendix D 1), 1 2 is the 2 × 2 identity matrix, k denotes the momentum with the appropriate dimensions, and d 0 and d = (d x , d y , d z ) are complex-valued momentum dependent variables. In the following, we drop the momentum dependence for the purpose of brevity, and reinstate it when needed.
The characteristic polynomial given in Eq. (1) in this case reads The degenerate points are then obtained by setting the discriminant of F λ (k) in Eq. (7) to zero, i.e., The defective degenerate points are the EP2s. Without loss of generality, we can set tr[H] = 2d 0 = 0. To find EP2s, we introduce two constraints based on the real (η) and imaginary (ν) parts of D[H] as Here η and ν describe N -spatial-dimensional surfaces. Note that in Hermitian systems, i.e., d I = 0, band touching points occur when all components of d R vanish amounting to at most three constraints.
In the vicinity of EP2s, the Hamiltonian casts the perturbed Jordan block with dimension n = 2 in Eq. (4) and reads Using the similarity transformation H 0 = S∆S −1 , the dispersion relation close to the EP2s yields ± − det[H(k)], which are the diagonal elements of ∆. This result is central in various studies on systems with det[H(k)] = −|k| due to the nonanalytical energy dispersion [54].
In two spatial systems, solutions to η = ν = 0 [cf. Eq. (10)] describe two closed curves in k space, such that EP2s appear when these curves intersect. In three spatial dimensions, the intersection between the twodimensional surfaces described by η = 0 and ν = 0 forms a closed exceptional curve, which can give rise to exceptional knots and result in exotic features such as open real/imaginary Fermi surfaces [54].
We now turn to the symmetries listed in Table I and see how the presence of one or the coexistence of multiple symmetries constraints the appearance of EP2s. This problem for EP2s was also studied in Ref. 39 for the symmetries defined by Bernard and LeClair (BLC) [36] (see also Appendix A). We cast it here in the form of the symmetries as given in Table I, which also includes additional symmetries to the BLC classification. To demonstrate our procedure, we treat two symmetries explicitly as well as their combination in the following.
As an example, we start by considering PHS † symmetry with T − T * − = −1. We choose T − = i σ y , such that the most generic form of a particle-hole (PH)-symmetric Hamiltonian from Eq. (5) casts Here we have introduced an additional label onto d, where each of the d parameters is represented as . The trace and determinant then read Setting the discriminant in Eq. (9) to zero (D[H PHS † ] = 0), we immediately find modified η, ν constraints, which are The presence of PHS † thus does not reduce the number of constraints for finding EP2s but merely restricts the momentum-dependency of parameters d.
If we instead consider P symmetry with P = σ x , the most generic form of a parity-symmetric Hamiltonian reads The trace and determinant then read Here dO = dOR + idOI with O ∈ {x, y, z}. Symmetric and antisymmetric components of dO with respect to k → −k are labelled by dOαs and dOαa with α ∈ {R, I}, respectively. η and ν are introduced in Eq. (10). Note that non-zero parameters might vary by changing the chosen Pauli matrix for each symmetry operator, an example of which is presented for the parity symmetry for which we include two representations. Nevertheless, the number of parameters and constraints remain intact. Here To find EP2s, we satisfy η and ν constraints for this system as Similar to PHS † , P symmetry puts restrictions on the momentum dependency of the d i 's, while not reducing the number of constraints for realizing EP2s.
If we now consider the presence of both PHS † with T − = i σ y and P symmetry imposed by σ x , we get The trace and determinant then read and we find Clearly one merely should satisfy η = 0 to find EP2s in this system. Therefore, even though PHS † and P symmetry individually do not reduce the number of constraints, the combination of these symmetries leaves only one constraint nonzero.
We summarize the results for these and the other symmetries in Table III. There we specify the symmetry generator and number of nonvanishing constraints and d parameters in the presence of each symmetry. Table IV summarizes various combinations of psH symmetry and other symmetries in the system. We note that one can simply compare the number of parameters in the fourth column of Table III to see which terms survive in the presence of multiple symmetries. As expected, the results in Table III are in agreement with our general findings  regarding EPns in Table II.
To demonstrate our findings in this section by a concrete example, we now look at an effective description of the driven-dissipative Kitaev model presented in Ref. 29.
Here the traceless Hamiltonian is given by where k stands for the momentum index, J is the nearestneighbor hopping amplitude, µ denotes the chemical potential, and γ l and γ g are, respectively, loss and gain coupling rates between the 1D system and the dissipative reservoir. This model displays TRS † with generator σ z , PHS † with generator 1, and CS with generator σ z [29].
The trace and determinant of H ddK read As a result, the η and ν constraints cast One can find nonzero d parameters by keeping common nonzero parameters given by each symmetry individually presented in Table III. While we only list the coexistence of psH symmetry with other non-Hermitian symmetries, the found recipe for determining the number of parameters is generic.

IV. EPS IN THREE-BAND SYSTEMS
To study EPs of order three, we perform a matrix decomposition in the Gell-Mann basis. Within this decomposition, the most generic three-band Hamiltonian is given by The three solutions λ 1 , λ 2 , λ 3 of F λ are eigenvalues of H in Eq. (29) and are given explicitly in Appendix E. The associated discriminant for Eq. (30) then casts Here, the complex-valued constraints read In the presence of symmetries, the number of nonzero constraints may reduce and different d's may vanish. Table V summarizes these constraints and the number of nonzero parameters in Hamiltonians with a specific symmetry, listed in Table I. As before, although we depict a particular symmetry generator for each symmetry, the number of constraints and nonzero parameters do not depend on our choice of generator. This can be explicitly seen for PT symmetry, for which we have presented two possible symmetry operators. Similar to the case of EP2s, three-band touchings occur in the Hermitian case (d I = 0) when d R = 0 amounting to at most eight constraints.  We summarize these three types of EP3s in Table VI. Without reference to symmetries, types I and III were also reported in Ref. 55.
2 (ηR, νR) 12 (d1I , d2R, d2I , d3R, d3I , d4R, d5R, d5I , d6R, d6I , d7I , d8R) Here dO = dOR + i dOI with O ∈ {x, y, z}. Symmetric and antisymmetric components of dO with respect to k → −k are labelled by dOs and dOs, respectively. Complex-valued η and ν are introduced in Eqs. (32,33). νR, ηR stand for the real components of ν, η. Note that nonzero parameters might vary by changing the depicted symmetry operators for each symmetry, see, e.g., the two different choices for PT symmetry. Nevertheless, the number of parameters and constraints remain intact. Aside from EP3s, three-band systems may also host EP2s. To explore the conditions in which these EPs can be realized, we introduce a subclass of traceless 3 × 3 Hamiltonians, which read where the second factor originates from h 2×2 . For this factor, we can write the companion matrix which explicitly shows the possibility of observing EP2s in this subsystem of three-band Hamiltonians.
To explore the effect of imposing symmetries on the behavior of EPs and the associated conditions for their appearance, we introduce an explicit three-band model in the following. Our model Hamiltonian reads where h x = α x + i sin(k x ), h y = α y + i sin(k y ), and h z = α z + i(−2 + cos(k x ) + cos(k y )). Our model is a non-Hermitian generalization of the effective Hamiltonian for three-fold fermions at k z = π/2 introduced in Ref. 13. This model Hamiltonian displays the pseudo-chiral symmetry with generator −1 3 and hosts a threefold degeneracy in its Hermitian spectrum at α x = α y = α z = 0, as shown in Fig. 1 (a). The traces and the determinant of this model read For the purpose of simplicity, we set α x = α, α y = α and α z = i √ 2α 2 with α a real-valued number. The real and imaginary parts of the eigenvalues with α = 0.3 are shown in Figs. 1(b,c) and Figs. 2(a,b), respectively, and reveal that our system exhibits an EP3 when (k x , k y ) → 0. Based on the low-energy dispersion of the spectrum with one flat (with energy zero) and two dispersive bands, -π -π/2 0 π/2 π The factor which is under the square root in 2 and 3 is − tr[H 2 ]/2. Thus, our model in Eq. (39) indeed gives rise to type I EP3s.
Imposing PT symmetry with generator 1 3 /3 + M 7 − M 8 / √ 3 = diag(1, −1, 1) on this model leads to a pseudochiral-PT -symmetric Hamiltonian, which reads where h α = √ 2α + cos(k x ) + cos(k y ) − 2 . The band structure of this system at α = 0.3 is plotted in Figs. 1 (d,e) and Fig. 2 (c,d). Even though we observe three-band crossings in the band structure of this system, we emphasize that EP2s, instead of EP3s, emerge at momenta slightly away from the origin. To demonstrate this statement, we look at the characteristic polynomial, which reads where 2α cos(k y ) − cos 2 (k y ) + 4 cos(k y ) − 4. The characteristic polynomial in Eq. (47) factorizes into a first order and a second-order polynomial, similar to Eq. (37). This means that it is possible to find a unitary transformation such that the Hamiltonian matrix in Eq. (47) features a zero row and zero column, or in other words, that the zero-energy flat band is not coupled to the other bands. We also note that Ω α (k x , k y ) = 0 delineates the region in which EP2s exist.
Lastly, we consider the case of finding EP2s in a threeband model in the presence of SLS. We start with the Hamiltonian in Eq. (35), which is SL-symmetric with The eigenvalues for this Hamiltonian read 0, ± √ b 2 − cd. Even though the eigenvalues are three-fold degenerate when b 2 = −cd, only two eigenvectors coalesce onto one at this point, and we thus find an EP2 in the system. There are two important things to note for this example. Firstly, it is only possible to find EP2s in threeband models with SLS as long as the zero-energy band is not coupled to the other bands, or in other words, as long as the three-band model can be described by a Hamiltonian like Eq. (35). Indeed, if the zero-energy band were to be coupled to the other bands such that the most generic three-band SL-symmetric Hamiltonian reads , any previously existing EP2 would immediately be promoted to an EP3. Secondly, to retrieve H 1,SLS , one has to tune d 1 = d 4 = 0. This means that to find an EP2 in a three-band SL-symmetric model, one has to satisfy six real constraints, namely the two constraints, Re[b 2 ] = − Re[cd], Im[b 2 ] = − Im[cd], that one needs to satisfy to find an EP2 in the presence of SLS (cf.
To exemplify the role of symmetries on the low-energy dispersion of EP4s, we present a case study in the following. We start with a traceless non-Hermitian four-band These constraints simultaneously vanish when k x = k z = 0. As a result, EP4s appear in this system as we have shown in Figs. 3(b,c) and Fig. 4(a,b). As close to this EP4, η, ν and κ are nonzero, we identify this EP4 as type 0, see Table (IX). Aside from EP4s, our model also exhibits EP2s close to k x = k z ≈ 0.47, as shown in Figs. 3(b,c) and Fig. 4(a,b).
Here we see that bands are doubly degenerate come in pairs as merely two bands are visible. The momenta at which EP2s occur are k x = ±α. This can be obtained from the associated constraints for H psH η and ν are zero when k x = ±α. Finally, in agreement with our findings, the number of constraints reduces when we impose psH symmetry to our non-Hermitian system in Eq. (64).

VI. DISCUSSION AND CONCLUSION
In this work, we have studied the appearance of exceptional points of any order in the presence of symmetries. In particular, we have addressed three questions pertaining to the number of constraints to find EPns, the implications of symmetries on the number of constraints for realizing EPs, and the low-energy behavior of these EPs. By expressing the characteristic polynomial of an n-dimensional non-Hermitian Hamiltonian in terms of the determinant and traces of the Hamiltonian, we have shown that one can identify 2n − 2 real constraints for finding EPns. We, furthermore, have discussed that in the presence of various symmetries, the number of constraints may reduce. Our results show that combining symmetries generally results in further decreasing the number of constraints. By interpreting the companion matrix as a perturbation close to an EPn, we have explicitly identified plausible low-energy dispersions of EPns. Besides these general considerations for EPs of any order, we have derived exact results for EPs of orders two, three, and four. Through looking at the companion matrix, we have also calculated explicit expressions for the dispersion around an EP, allowing us to characterize EP3s and EP4s based on their low-energy spectrum. In addition, we have presented the appearance of lowerorder EPs in n-dimensional models and find that EP2s can be realized in three-band systems, while both EP2s and EP3s can appear in four-band systems.
While we have focused on EPns in this work, we emphasize that our results can be straightforwardly generalized to exceptional structures of higher dimensions. Associating a parameter with each constraint, we have shown that EPns generally appear in n − 1-dimensional setups in the presence of, e.g., psH symmetry. Consequently, exceptional one-dimensional lines or two-dimensional surfaces of order n appear generically in n-and n + 1dimensional systems, respectively. In other words, the number of constraints is related to the codimension of the exceptional structure, i.e., the difference between the total dimension of the system and the dimension of the exceptional structure, and our results can thus be readily applied to study the realization of higher-dimensional exceptional structures in the presence of symmetries.
Besides exceptional degeneracies, ordinary (Hermitian) degeneracies may appear in non-Hermitian systems where the eigenvalues coalesce, but the eigenbasis is complete. As we briefly discussed in Sects. III-V, this requires setting d = 0 for the various Hamiltonians such that these Hamiltonians are proportional to an identity matrix. Generally, this results in having to satisfy a large number of constraints to find these degeneracies. Indeed, one needs to satisfy 2(n 2 − 1) constraints to find an ordinary n-fold degeneracy, where n 2 − 1 is the number of dimensions of the group SU(n). Clearly, d = 0 is a solution to the characteristic polynomial in Eq. (1). We note that one of the crucial differences between EPs and ordinary degeneracies on the level of polynomial equations sits in the relation between the characteristic and the minimal polynomials: For EPs, the characteristic polynomial equals the minimal polynomial, whereas, for ordinary degeneracies, the characteristic polynomial is a multiple of the minimal polynomial [72].
In Ref. 34, it is stated that symmetry-protected multifold exceptional points are points at which the symmetry is spontaneously broken. This is indeed the case for the symmetries the authors consider there (CS, psH, PT , and CP symmetry), which are antiunitary symmetries that are local in parameter space. We here show that unitary, local symmetries such as SLS and psCS can also stabilize higher-order EPs in lower dimensions. These EPs do not mark a transition between broken and unbroken symmetry, thus showing that not all symmetry-protected EPns necessarily correspond to spontaneous symmetrybreaking points.
While we have presented an extensive study here on the realization of exceptional points of any order in the presence of symmetry, we did not touch upon the possibility of defining topological invariants. Former studies proposed to define Z 2 index based on either sign(det[H]) (sign(det[i H])) in two-band models with PT (CP) symmetry [38,73] or the sign of the discriminant [34] for systems of any dimension with CS, psH, PT , and CP symmetry. It would be intriguing to investigate whether more generic invariants could be defined based on our rigorous mathematical framework. We leave this problem for later studies.
Acknowledgments.-We would like to thank Emil J. Bergholtz for pointing out Ref. 60. Bernard and LeClair define non-Hermitian symmetries as follows [36].
When n = 2j + 1, at least one of the eigenvalues needs to be zero, such that det[H] = 0. We also find tr[H k ] = 0, ∀k ∈ odd as before. We thus are left with n−1 constraints that need to be satisfied to find an EPn with n = 2j + 1.

Parity-time and pseudo-Hermitian symmetries
In the presence of PT or psH symmetry, eigenvalues satisfy { (k)} = { * (k)}. This implies that for n = 2j +1 at least one of the eigenvalues should be real. We save this real eigenvalue in j+1 for n = 2j +1 in the following.

Chiral and parity-particle-hole symmetries
In the presence of CS or CP symmetry, eigenvalues display { (k)} = {− * (k)}. We note that CS is not defined for n = 2j + 1 as a result for odd dimensions in the following are merely relevant for the CP symmetry. For n = 2j + 1, we infer from the relation between the sets of eigenvalues that at least one of the eigenvalues is imaginary. We save this eigenvalue in j+1 for n = 2j +1.