Topological semimetal phase with exceptional points in 1D non-Hermitian systems

To describe eigenstates in non-Hermitian crystalline systems, the non-Bloch band theory has recently been established, and the generalized Brillouin zone (GBZ) has unique features which are absent in Hermitian systems. In this Letter, we show that in one-dimensional non-Hermitian systems with both sublattice symmetry and time-reversal symmetry, a topological semimetal phase with exceptional points is stable. This stems from the unique features of the GBZ. We can relate the motion of the exceptional points with the change of the value of a topological invariant characterizing a topological insulator phase. It is also shown that the energy bands are divided into three region, depending on the symmetry of the eigenstates.

To describe eigenstates in non-Hermitian crystalline systems, the non-Bloch band theory has recently been established, and the generalized Brillouin zone (GBZ) has unique features which are absent in Hermitian systems. In this Letter, we show that in one-dimensional non-Hermitian systems with both sublattice symmetry and time-reversal symmetry, a topological semimetal phase with exceptional points is stable. This stems from the unique features of the GBZ. We can relate the motion of the exceptional points with the change of the value of a topological invariant characterizing a topological insulator phase. It is also shown that the energy bands are divided into three region, depending on the symmetry of the eigenstates.
Non-Hermitian quantum mechanics has been attracting much attention in many fields of physics in the past decades. Many experimental studies have realized various physical systems with non-Hermitian effects 1- 19 . Among these experimental studies, appearance of exceptional points and rings where some energy eigenvalues become degenerate and the corresponding eigenstates coalesce 20,21 , and intriguing phenomena have been observed [22][23][24][25][26][27][28][29][30][31][32][33] . At such non-Hermitian degeneracy, since the Hamiltonian is non-diagonalizable, these degeneracies are unique to non-Hermitian systems. As is motivated by these experimental studies, existence of exceptional points, rings, and surfaces, and phenomena induced by them have been theoretically predicted in various physical systems  , and under some symmetries, they are classified in terms of topology 62-66 .
Recent theoretical studies have been focusing on topological systems in solid state physics 67-79 . The bulkedge correspondence has been notably under debate  since it seems to be violated in contrast to Hermitian systems. The main issue of the bulk-edge correspondence in non-Hermitian systems is that there is a difference between the energy spectrum in a periodic chain and that in an open chain. This difference is caused by the non-Hermitian skin effect 83 . In Refs. 83 and 102, it was shown that while the Bloch wave number k takes a real value in a periodic chain, it becomes complex in an open chain, and that the value of β ≡ e ik is confined on a loop on the complex plane so that continuum bands are reproduced in a large open chain. Then the loop of β is called the generalized Brillouin zone (GBZ) C β , which is a generalization of the Brillouin zone in Hermitian systems. We note that C β is deformed as the system parameters change and that it can have cusps 102 . As a result, one can establish bulk-edge correspondence between the topological invariant defined in terms of C β and existence of the topological edge states in some cases.
In this Letter, we show that the non-Hermitian modification of the GBZ C β gives qualitative changes to physics of topological semimetals (TSMs). Because of this modification, the TSM phase with exceptional points is stable in one-dimensional (1D) non-Hermitian systems. The phase is topologically protected by both sublattice symmetry (SLS) and time-reversal symmetry (TRS), and therefore, it disappears if either the SLS or the TRS is broken. As system parameters change continuously, C β is deformed so that gapless points always lie on C β , and the system remains in the TSM phase. This TSM phase can be regarded as an intermediate phase between a normal insulator (NI) phase and a topological insulator (TI) phase characterized by a topological invariant. Then the creation and annihilation of the exceptional points can be related to the change of this topological invariant. We also find that the continuum bands are divided into three regions, where the energies become real, pure imaginary, and complex, depending on the symmetry of the eigenstates.
For illustration of the stabilization of the TSM phase, in Fig. 1, we show an evolution of our model, whose details are explained later, upon a change of a system parameter µ. In Fig. 1(c-1)-(h-1), a gap in the system closes if the GBZ C β goes through one of gap-closing points represented by the red or blue dots and squares. In particular, one can see from Figs. 1(e-1) and (f-1) that as µ changes, C β is deformed so that it keeps the gap closed. We also find that the phase transitions between the NI and TI phases are qualitatively different from Hermitian systems. For example, in Fig. 3(b), at the phase transition, the exceptional points come to the cusps, and in total, the gap closes at three points on C β . These novel behaviors never occur when the Bloch wave number takes real values. Since in such a case, the Brillouin zone becomes a unit circle regardless of the values of the system parameters, the Brillouin zone cannot hold the gap-closing point (Figs. 2(c) and (d)), and the gap closes only at isolated points along the µ-axis (Figs. 2(b) and (e)).
In the following, we study a 1D non-Hermitian tightbinding system with both the SLS and the TRS. For a real-space Hamiltonian H, these symmetries are defined as ΓHΓ −1 = −H, T H * T −1 = H, where Γ and T are unitary matrices. Here, due to the SLS, one can write where R ± (β) are N × N matrices. Then the eigenvalue equation det [H (β) − E] = 0 yields the bands symmetric with respect to E = 0. Next we mention the main point of the non-Bloch band theory in Ref. 102. In non-Hermitian systems, the Bloch wave number k becomes complex so that it describes continuum bands in a long open chain. The set of k forms the GBZ C β . We note that the eigenvalue equation det [H (β) − E] = 0 is an algebraic equation for β with an even degree 2M in general, and we number the 2M solutions of det [H (β) − E] = 0 so as to satisfy |β 1 | ≤ · · · ≤ |β 2M |. Then the condition for continuum bands is given by and the trajectories of β M and β M+1 give C β . Now we describe the reason why the TSM phase is stable under the SLS and TRS. From the Bloch Hamiltonian (1), a condition for a gap closing at E = 0 is decomposed into two equations det R + (β) = 0 and det R − (β) = 0. Thanks to the TRS, det R ± (β) are polynomials of β and β −1 with real coefficients, and therefore, it follows that any complex solutions of det R ± (β) = 0 appear in complex conjugate pairs (β, β * ). Then, if we suppose β M and β M+1 form a pair of the complex conjugate solutions of det R + (β) = 0, they satisfy Eq. (2), meaning that β M and β M+1 are on the GBZ, and the gap is zero. Even when system parameters change, the gap remains zero as long as this pair gives M th and (M + 1)th largest absolute values among the 2M solutions. In conclusion, the gapless region in the non-Hermitian system with the SLS and the TRS is robust against the change of system parameters.
To study how the TSM phase appears, we investigate the non-Hermitian Kitaev chain, which has been studied in some previous works 104,111-119 . Its real-space Hamiltonian is written as This Hamiltonian reduces to that of the conventional Kitaev chain 120 when We assume that the real-space Hamiltonian (3) satisfies both the SLS and the TRS, which sets all the parameters to be real. Then the Bloch Hamiltonian H (β) is expressed as the off-diagonal form (1) with Im(β) Since it is a quartic equation for β, the condition for continuum bands can be written as |β 2 | = |β 3 | when the solutions satisfy |β 1 | ≤ |β 2 | ≤ |β 3 | ≤ |β 4 | 102 . We note that the GBZ is always a single loop encircling the origin on the complex plane 106,121 . The insulator phases of this system are classified in terms of a Z topological invariant called winding number w because it has the SLS, defined as 102 where [arg R ± (β)] C β means the change of the phase of the functions R ± (β) as β goes along the GBZ C β in a counterclockwise way. As long as there is a gap at E = 0, R ± (β) never vanish along C β , and w is well defined. Figure 1(a) shows the phase diagram of our model, with the NI phase with w = 0 (white region), the TI phase with w = 1 (blue region), and the TSM phase (orange region). In the TSM phase, the gap closes at E = 0, which means that the equation R + (β) = 0 or R − (β) = 0 holds somewhere on the GBZ C β . At such a point on C β , the Hamiltonian cannot be diagonalizable, and such point is called exceptional point. In other words, the orange region represents that the system has the exceptional points. We note that the TSM phase appears as an intermediate phase between the NI and TI phases. Now we explain the mechanism of the appearance of this TSM phase. In our model, the solutions of the equa- tions R ± (β) = 0 are gap-closing points, shown as the red and blue dots and squares in Figs. 1(c-1)-(h-1), and they become the exceptional points when the GBZ C β goes through them. Let β = β a i (i = 1, 2, a = +, −) denote the gap-closing points of R a (β) = 0, with |β a 1 | ≤ |β a 2 |. In the regions A and B in Fig. 1(a), R − (β) = 0 has two complex-conjugate gap-closing points, satisfying β − 1 = β − 2 * , and their common absolute value is between the values of β + 1 and β + 2 : Thus the condition (2) is satisfied, and therefore, β − 1 and β − 2 are on C β and are the exceptional points. The condition β − 1 = β − 2 * remains satisfied even when the system parameter changes because R − (β) = 0 is an algebraic equation with real coefficients. Thus the exceptional points move along C β as shown in Figs. 1(c-2)-(h-2), and the system remains in the TSM phase. In the region C and D in Fig. 1(a), a similar scenario holds true, by exchanging R + (β) and R − (β). Because of this mechanism for the stabilization of the exceptional points, annihilations (and likewise creations) of them are limited to two patterns as shown in Figs. 3(a) and (b). Figure 3(a) represents a coalescence of two exceptional points, and Fig. 3 (b) represents an encounter between the gap-closing point and the cusp. In Fig. 3(a), the two exceptional points meet and become two real gap-closing points. It occurs on the real axis. This can be seen in Fig. 1(g). On the other hand, the case of Fig. 3(b) occurs when two complex-conjugate exceptional points and one gap-closing point share the same absolute value. Therefore the gap closes at three points on the GBZ, for example, as shown in Fig. 1(d). At this point, Finally we show that in non-Hermitian systems with the SLS and TRS, the continuum bands are divided into three regions in terms of the symmetry of the eigenstates. In the first one, the energies are real, and eigenenergies of a time-reversal pair |ψ , T |ψ * are degenerate, i.e., E = E * , where |ψ is an eigenstate of the Hamiltonian. In the second region, the energies are pure imaginary, and a pair of states |ψ and ΓT |ψ * related by the sublattice-time-reversal symmetry (STS) is degenerate, i.e., E = −E * . In the third region, the energies are complex, and neither the time-reversal pair nor the sublattice-time-reversal pair is degenerate. We call these three regions TRS-unbroken region, STS-unbroken region, and TRS/STS-broken region, respectively. For example, in the non-Hermitian Kitaev chain, we show the above three regions (Figs. 4(a-2)-(b-2)).
Importantly, these regions are connected to each other at the cusps or the exceptional points. In Fig.4(a-2) and (b-2), three curves meet at one point, where the GBZ C β have cusps represented by A 1 and B 1 (or A 2 and B 2 ). It is consistent with the property of the cusp, where three points on C β share the same energy. Furthermore, in Fig. 4(b), the green and blue lines are connected at the exceptional point with E = 0. Thus the exceptional point connecting the real and pure-imaginary energies becomes stable because such structure is topologically protected by the symmetries.
In summary, we show that in 1D non-Hermitian systems with both the SLS and the TRS, the TSM phase with exceptional points is stable, unlike Hermitian systems. The appearance of the TSM phase is attributed to the unique features of the GBZ. It is shown that this TSM phase can be regarded as an intermediate phase between the NI and TI phases. We also find that the continuum bands are divided into three regions in terms of the symmetry of the eigenstates, and the regions change only at the cusps and the exceptional points. Thus non-Hermiticity brings about qualitative changes to the topological phase transition.
So far, we have treated the parameters ∆ b and ∆ f as real. When ∆ b , ∆ f ∈ C, both the SLS and the TRS are broken. In this case, this system has only pseudo particle-hole symmetry (PHS) 110 . We show that the non-Hermitian systems with the pseudo PHS is classified in terms of a Z 2 topological invariant, and importantly, we find that the pseudo PHS cannot stabilize the TSM phase with exceptional points. We show some analyses in the Supplemental Material 122 .
We note that whether such TSM phase appears or not depends on symmetries of systems. When the systems have either chiral symmetry, pseudo TRS, or PHS 110 , the Bloch wave number becomes real, independent of any boundary conditions 101,122 . Therefore the TSM phase cannot stably exist under either of these symmetries.
We are grateful for Ryo Okugawa and Ryo Taka In this section, we show that in one-dimensional (1D) non-Hermitian systems with both sublattice symmetry (SLS) and time-reversal symmetry (TRS), the value of the winding number changes through the creation and annihilation of the exceptional points as shown in Fig. S1. Here, for a real-space Hamiltonian H, the SLS and the TRS are defined as where Γ and T are unitary matrices representing the SLS and the TRS, respectively.

A. General cases
First of all, we focus on a two-band model. Due to the SLS, we can write down the Bloch Hamiltonian H (β) of this system as where β ≡ e ik , k ∈ C, and R ± (β) are holomorphic functions for β. In the following, without loss of generality, we can write these functions as where C ± are real constants due to the TRS. Then the eigenvalue equation det [H (β) − E] = R + (β) R − (β) − E 2 = 0 can be explicitly written as which is an algebraic equation for β with a degree 4m.
Here, by numbering the solutions of Eq. (S4) so as to satisfy |β 1 | ≤ · · · ≤ |β 2M |, the condition for continuum bands 1 is given by and one can get the generalized Brillouin zone (GBZ) C β from the trajectories of β M and β M+1 . We note that C β always encircles the origin on the complex plane 2,3 . In this case, this system is classified in terms of a Z topological invariant called winding number w, and it can be defined as 1 where [arg R ± (β)] C β means the change of the phase of the functions R ± (β) as β goes along C β in a counterclockwise way. We note that as long as there is a gap at E = 0, R ± (β) never vanish along C β , and w is well defined. On the other hand, when the gap closes at E = 0, the gap-closing condition can be decoupled into two equations R + (β) = 0 and R − (β) = 0 being satisfied somewhere on C β . At such a point, the Bloch Hamiltonian (S2) cannot be diagonalizable, and such point is called exceptional point. In this case, w is not well defined. From Eq. (S3), we can rewrite this form as where N ± zeros expresses the number of the solutions β ± i of R ± (β) = 0 inside C β , respectively. Furthermore it is worth noting that from Ref. 2, when the system has a gap around E = 0, we can get where N zeros (= N + zero + N − zero ) expresses the number of solutions of the equation det H (β) = 0 inside C β . Equation (S8) tells us that the total number of the solutions of det H (β) = 0 inside C β , namely the number of the gap-closing points inside C β , is unchanged as long as the system has a gap. Now we focus on two insulator phases separated by the topological semimetal (TSM) phase with exceptional points. When the system enters the TSM phase from one of the insulator phases, the exceptional points are created by the inverse process as shown in Fig. S1(a) or  (b). Then, after a further change of the system parameters, the system becomes the other insulator phase from the TSM phase, and here, the exceptional points are annihilated by the process as shown in Fig. S1(a) or (b). If the creation is by the inverse process of Fig. S1(a) and the annihilation is by the process of Fig. S1(b) (or vice versa), while the total number N zeros (= N + zero + N − zero ) of the solutions of det H (β) = 0 inside C β is unchanged, such a motion of the exceptional points change the value of N + zeros and N − zeros by 1 and −1 (or by −1 and 1), respectively, resulting in the change of w by 1 (or by −1). Therefore we conclude that these insulator phases have different values of w. We note that this scenario can be extended to multi-band systems.

B. Non-Hermitian Kitaev chain
We show an example of the non-Hermitian Kitaev chain with the SLS and TRS introduced in Eq. (4) in the main text. This model has the Bloch Hamiltonian H (β) as and all the parameters are real. We note that the gap closes at E = 0 due to the SLS when either of the two quantities R ± (β) is zero. We show the phase diagram of this model in Fig. S2(a) and its GBZ C β , the solutions of the equations R ± (β) = 0, and the motion of the exceptional points in Fig. S2(c)-(h). In the topological phase transition along the black arrow in Fig. S2(a), we show the position of the solutions β + 1 and β + 2 of R + (β) = 0 (and β − 1 and β − 2 of R − (β) = 0) as the red (and blue) dots and squares, respectively, as shown in Figs. S2(c-1)-(h-1). When we decrease the value of the parameter µ, at µ = 2.683 ( Fig. S2(g-1)), β + 1 and β + 2 change from real values to complex values via coalescence, and it corresponds to a pair creation of the exceptional points. After the coalescence, β + 1 and β + 2 become complex, with β + 1 = β + 2 * , and their common absolute value is between the values of β − 1 and β − 2 , meaning that β + 1 and β + 2 stay on C β . Then, at µ = 0.5813 (Fig. S2(d-1)), the value of β − 2 becomes equal to that of β + 1 = β + 2 , and after passing that point (i.e., µ becomes less than 0.5813), β − meaning that β + 1 and β + 2 are no longer on C β , and the gap opens (Fig. S2(c-1)). At the phase transition point µ = 0.5813 (Fig. S2(d-1)), three solutions β + 1 , β + 2 , β − 2 share the same absolute value, and the gap closes at three points on C β , and at such points, C β has cusps. The change of the winding number w defined as Eq. (S6) readily follows from the following argument. In the normal insulator (NI) phase with w = 0, C β surrounds one solution of R + (β) = 0 and one solution of R − (β) = 0 (Fig. S2(h-1)). In this case, On the other hand, in the topological insulator (TI) phase with w = 1, there exist two solutions of R − (β) = 0 and no solutions of R + (β) = 0 inside C β (Fig. S2(c-1)), which leads to w − = 1, w + = −1, and w = 1 as expected. Therefore the creation and annihilation of the exceptional points changes the number of the solutions of R ± (β) = 0 inside C β , and the value of the topological invariant also changes.
We note that in addition to the topological phase transition between two insulator phases via the TSM phase, a direct phase transition from the NI phase to the TI phase is also possible as shown in the inset of Fig. S2(a).
Here the gap closes on the real axis, where R + (β) and R − (β) simultaneously become zero at this value of β.

SII. 1D NON-HERMITIAN SYSTEM WITH PSEUDO PARTICLE-HOLE SYMMETRY
In this section, we investigate a 1D non-Hermitian system with pseudo particle-hole symmetry (PHS) and show that it is classified in terms of a Z 2 topological invariant. Here, for a real-space Hamiltonian H, this symmetry is defined as where C is the unitary matrix. In the following, we focus on a two-band model.

A. Two-band model
The Bloch Hamiltonian H (β) can be written as where σ 0 is a 2 × 2 identity matrix, σ i (i = x, y, z) are the Pauli matrices, and the complex Bloch wave number is defined as β ≡ e ik , k ∈ C. We assume that it satisfies and then, the coefficients H i (β) (i = 0, x, y, z) in Eq. (S12) satisfy (S14) It is worth noting that H i (β) (i = 0, x, y) are a pure imaginary and H z (β) is real when arg β = 0 and arg β = π.
B. Z2 topological invariant Such a system is classified in terms of the Z 2 topological invariant ν (= 0, 1), and for the Bloch Hamiltonian (S12), we can define it as where β 0 (or β π ) is the value of β at arg β = 0 (or arg β = π) on the GBZ C β (for example, see Fig. S3). In Eq. (S15), the integral contour β goes along C β , and we select the branch cut of the square root so that the both functions R ± (β) becomes continuous on C β . In the following, we assume that a system has a gap. Here we note that two continuum bands are separated by a line which determines a complex gap on the complex energy plane. This gap is called a line gap 5 . In this case, one can show that ν takes only 0 or 1. To this end, we calculate the value of exp (2πiν). As mentioned in Sec. SII A, since the functions R ± (β 0 ) and R ± (β π ) take real values, we can rewrite the expression of exp (2πiν) as (S16) Here we note that the quantities i H 2 i (β 0 ) and i H 2 i (β π ) are real. Since we assume presence of a line gap, we conclude that i H 2 i (β 0 ) and i H 2 i (β π ) have the same sign, as we prove by contradiction in the following.
Suppose i H 2 i (β 0 ) and i H 2 i (β π ) have different signs. We can set i=x,y,z without loss of generality. First of all, we assume H 0 (β) = 0 for simplicity. At β = β 0 , the energies are , and we change β along C β in a counterclockwise way from β = β 0 to β = β π . Here, let C + denote this path on the complex plane. Then, at β = β π , the energy is given by E = ε π = i H 2 i (β π ) ∈ iR, where the branch of the square root is chosen in such a way that E = i H 2 i (β) is continuous along C + . Next we consider a path C − along C β in a clockwise way from β = β 0 to β = β π . Because , the energy at β and that at β * are complex conjugate. Therefore, because (C + ) * = C − , the energy at β = β π along C − is E = ε * π = −ε π . Thus, by encircling C β (= −C − + C + ) once from β π to β π , the branch changes from E = −ε π to E = ε π , meaning that the two energies E = ± i H 2 i (β) are continuously connected, and there is no line gap, contradicting the assumption. Thus we can conclude that i H 2 i (β 0 ) and i H 2 i (β π ) have the same sign. So far, we assume H 0 (β) = 0, but even in the case of H 0 (β) = 0, because this term does not affect above argument, the above proof remains valid.
Therefore, from Eq. (S16), we can get exp (2πiν) = 1 and can conclude that the value of ν is an integer. As a result, ν can take only 0 or 1 (mod 2), and we conclude that ν can be interpreted as the Z 2 topological invariant in this system.
In particular, in Hermitian cases, we can greatly simplify the formula of ν in Eq. (S15). The Bloch wave number k becomes real, and the all coefficients included in Eq. (S12) become real functions. In the following, we replace H i (β) by H i (k) (i = 0, x, y, z) , k ∈ [−π, π].
Here the system has the conventional PHS represented as and H i (k) (i = 0, x, y) become zero at k = 0 and k = π. Furthermore we can get and then, Eq. (S15) can be rewritten as Since Eq. (S18) tells us that both argH z (0) and argH z (π) take 0 or π (mod 2π), Eq. (S20) can be further rewritten as and we obtain the known formula 6 (−1) Since Eq. (S25) is a quartic equation for β, the condition for continuum bands can be written as |β 2 | = |β 3 | when the solutions of Eq. (S25) satisfy |β 1 | ≤ |β 2 | ≤ |β 3 | ≤ |β 4 |. We note that the trajectories of β 2 and β 3 give the GBZ C β . The examples of C β and the continuum bands are given in Fig. S3 and in Fig. S4(a-3), respectively. Let ℓ ± denote the loops drawn by the functions R ± (β) on the complex plane when β goes along C β in a counterclockwise way. Equation (S15) tells us how to determine the value of ν; when neither ℓ + nor ℓ − surrounds the origin O on the complex plane, ν is equal to 0, and when two loops simultaneously surrounds O, ν is equal to 1. For example, in the case of Fig. S4(a-4), ν becomes 1. We note that the system has the exceptional points when either ℓ + or ℓ − passes O, and ν is not well defined in this case.
We can get the phase diagram in the generalized non-Hermitian Kitaev chain as shown in Fig. S4(a) and can confirm that the topological edge states appear when ν takes the nonzero value as shown in Fig. S4(b). Therefore we can establish the bulk-edge correspondence between the Z 2 topological invariant ν and existence of the topological edge states.
In this model, we can see that the pseudo PHS cannot protect a TSM phase with exceptional points. For simplicity, let the parameter γ be zero. We note that by treating the parameters t b and t f as complex, one can add H 0 (β) term to the Bloch Hamiltonian H (β) in the case of Fig. S4(a). In fact, for ℑ (t b ) = ℑ (t f ) = 0 being infinitesimal values as an example, we can obtain the phase diagram as shown in Fig. S4(b). We can confirm that the TSM phase with exceptional points in Fig. S4(a) disappears by adding H 0 (β). We note that the exceptional points appear on the orange lines on the phase diagram. However, this exceptional point can be removed by adding other perturbation terms. Therefore we conclude that the pseudo PHS cannot protect the TSM phase with exceptional points.

SIII. RESTRICTION ON THE BLOCH WAVE NUMBER BY EITHER PARTICLE-HOLE SYMMETRY OR CHIRAL SYMMETRY
In this section, we show that in 1D non-Hermitian tight-binding systems, some symmetries make the Bloch wave number real even in an open chain. The previous work 5 showed that it takes real values even in non-Hermitian cases under pseudo time-reversal symmetry (TRS) represented as where H is a real-space Hamiltonian, and T is a unitary matrix. Here we show that PHS or chiral symmetry (CS) also restricts the Bloch wave number to real values. In the following, we focus on a 1D real-space tight-binding where the unit cell is composed of q degrees of freedom, and the range of hopping is N . We note that this Hamiltonian can be non-Hermitian, meaning that t i,µν is not necessarily equal to t * −i,νµ . Then a Bloch Hamiltonian H (β) is given by [H (β)] µν = N i=−N t i,µν β i , (µ, ν = 1, · · · , q) , where the complex Bloch wave number is defined as β = e ik , k ∈ C. We assume that the eigenvalue equation det [H (β) − E] = 0 is an algebraic equation for β with an even degree 2M , and 2M solutions of this equation satisfy |β 1 | ≤ · · · ≤ |β 2M |. Here the condition for continuum bands is written as |β M | = |β M+1 |, of which the detail is given in the main text and Ref.