Topological Phase Transitions of Generalized Brillouin Zone

It has been known that the bulk-boundary correspondence (BBC) of the non-Hermitian skin effect is characterized by the topology of the complex eigenvalue spectra, while the topology of the wave function gives rise to Hermitian BBC with conventional boundary modes. In this work, we go beyond the known description of the non-Hermitian topological phase by discovering a new type of BBC that appears in generalized boundary conditions. The generalized Brillouin zone (GBZ) possesses non-trivial topological structures in the intermediate boundary condition between open and periodic boundary conditions. Unlike the conventional BBC, the topological phase transition is characterized by the generalized momentum touching of GBZ, which manifests as exceptional points. As a realization of our proposal, we suggest the non-reciprocal Kuramoto oscillator lattice, where the phase slips accompany the exceptional points as a signature of such topological phase transition. Our work establishes an understanding of non-Hermitian topological matter by complementing the non-Hermitian BBC as a general foundation of the non-Hermitian topological systems.

The rapid progress of non-Hermitian mechanics has paved a new direction to obtain unconventional phases of matter that have no Hermitian analog [1][2][3].The representative example is the non-Hermitian skin effect (NHSE), which exhibits the macroscopic localization of bulk states on the boundary.There have been theoretical efforts to understand non-Hermitian BBC using the bulk topology in the periodic boundary condition (PBC).It has been shown that the winding of complex eigenvalue in one-dimensional systems signifies the presence of the NHSE [4][5][6][7].In a more general sense, the macroscopic spectral collapse of the NHSE can be understood within the theory of generalized Brillouin zone (GBZ) defined in the complex plane.
So far, it has been shown that the BBC of non-Hermitian systems is characterized by two distinct types of topology: wave function topology in GBZ and complex eigenvalue spectra topology in conventional BZ.While the wave function topology gives rise to the Hermitian-like boundary modes, the winding number of the eigenvalue spectra in the complex plane signals the presence of the NHSE in OBC.Here, for the first time, we present the third type of non-Hermitian BBC that emerges in generalized boundary condition (GBC).Unlike the case of PBC or OBC, the single wave functions in GBC possess distinct momenta with multiple GBZs.The existence of multiple GBZ can provoke the intrinsic topological structure of GBZ [See Fig .1],where the topological phase transitions are accompanied by exceptional points with touchings of GBZs.The information of the complex eigenvalue spectra itself is ignorant of such concomitant topological phase transition.As a result, our result falls in neither of the two previously known categories of non-Hermitian BBC.
In this work, we formulate our theory by presenting the example of the Hatano-Nelson model in GBC.After * moonjippark@ibs.re.krFIG. 1.
Schematic illustration of non-Hermitian bulkboundary correspondence in generalized boundary condition(GBC).In the GBC, the GBZ can possess the nontrivial topological invariant, which manifests as the topological boundary mode.These topological boundary modes are not captured in the conventional band topology.The topological phase transition is characterized by the touchings of the GBZs, which manifests as exceptional points.
constructing the single-particle theory, we propose the non-reciprocal Kuramoto lattice as an example of manybody systems.The topological phase transitions between different time-dependent phases are accompanied by the touching of two generalized momenta.Unlike the energy band touching, the touchings of GBZs generally manifest as the coalescence of the two eigenstates, known as the exceptional point (EP).The emergence of the exceptional point in the dynamical spectral results in the novel class of non-equilibrium phase transition, known as exceptional transitions (ET) [8].Our result establishes the important connection between the non-Hermitian BBC and the ET in the non-equilibrium phases of many-body interacting collective behaviors.
Non-Bloch band theory in GBC-We start our discussion by considering the Hatano-Nelson (HN) model [9] arXiv:2305.08584v1[cond-mat.mes-hall]15 May 2023 with the GBCs.The Hamiltonian is given as, , where t R(L) represents the hopping toward the right (left) nearest neighbor site.The second line represents the general deformation of the boundary terms, which include the cases of OBC (t R = t L = 1 = N = 0), PBC (t R = t R , t L = t L , 1 = N = 0) and any intermediate boundary conditions between them (See I A for full classifications of the boundary conditions).
We represent the eigenstates as, |ψ = 1/ √ N N i=1 ψ i |i , where ψ i = α c α z i α and z α ∈ C is the generalized complex momenta.The eigenstates in the GBC are required to simultaneously satisfy bulk and boundary equations.In the case of the Hatano-Nelson model having the only nearest neighbor couplings, the bulk equation is given by the second-order complex polynomial as, which admits the pair of the complex momenta solutions, (z 1 , z 2 ).We call the two momenta paired if the two momenta (z 1 , z 2 ) satisfy the bulk equation with the same eigenvalue, E, and the paired momenta (z 1 , z 2 ) are related by the condition z 1 z 2 = r 2 , where r 2 = t R t L .In addition to the bulk equation, the eigenstates are required to satisfy the boundary equation, H B (c 1 , c 2 ) T = 0, where The coefficients of the boundary matrix H B are explicitly given as, A(z In practice, we determine the solutions for (z 1 , z 2 ) by solving Eq. ( 3), and subsequently solve for the eigenvalue with Eq. (2).
Except for PBC, where the translational symmetry is intact, a wave function is generally described by the two different paired momenta (z 1 , z 2 ) ∈ C 2 , which map to the same eigenvalue.The contour of different momenta constitutes the generalized Brillouin zones L 1 × L 2 [10].In the case of the Hatano-Nelson model, the GBZs form two circles with the radii |z 1 |, |z 2 | in the complex plane.For example, in the case of the OBC, the two GBZs overlap each other (|z 1 | = |z 2 | = r), which forms a single GBZ as discussed in the previous studies [4,5].
In contrast, in the case of the GBCs, the two GBZs can form two disjoint contours, which is the focus of our study.Generalized momenta z i outside (inside) the unit circle (|z i | = 1) or GBZs describes the evanescent wave localized on the left (right) boundary, which corresponds to the boundary states.Touchings of the unpaired momenta of the two GBZs correspond to the coalescence of the wave functions, which give rise to the EP (See the representative examples of GBZs and eigenvalue spectra for different GBCs in Sec.I C).
Topological boundary matrix-The presence of the multiple GBZs promotes the spinor representations of the eigenstates.For example, even though the system consists of a single site unit cell, the Bloch wave function can be represented as the expanded spinor (c 1 , c 2 ) T , where each component represents the coefficients of the paired momenta z 1 and z 2 respectively.Correspondingly, the boundary equation can be transformed as the eigenvalue equation, HB (c 1 , c 2 ) T = (c 1 , c 2 ) T , with the effective chiral symmetry as, where h + B (z 1 , z 2 ) = A(z 2 )/A(z 1 ), h − B (z 1 , z 2 ) = B(z 1 )/B(z 2 ).The transformed boundary matrix HB preserves the chiral symmetry by satisfying the following conditions, { HB , σ z } = 0.The presence of the chiral symmetry indicates that for any right (left) eigenstates |χ R(L) with the complex eigenvalue E(E * ), we can define the right (left) chiral partner states | χR(L) ≡ σ z |χ R(L) with the eigenvalue −E(−E * ), which satisfy the bi-orthogonality relations: There is a suggestive similarity between the boundary equation [Eq.( 4)] and the Bloch Hamiltonian of the non-Hermitian SSH model [4,5].We can assign the topological winding number W to classify the boundary matrix as, where [arg h ± B ] L1×L2 is the change of phase of h ± B as (z 1 , z 2 ) goes along the GBZ.In general, another nonzero topological invariant, W = (W + +W − )/2, can exist.However, for (z 1 , z 2 ) ∈ L 1 × L 2 , the GBZ ensures that W + = −W − due to the boundary equation det(H B ) = 0.As a result, W can only have non-trivial topological numbers.The physical manifestation of the non-trivial winding number manifests as the topological bound state localized on the boundary of the chain [See Fig. 2 (c) and also Sec.I D for more details on the symmetry and topology of the boundary matrix ].
In the limit of GBC-(I) (t L(R) = 0, 1 , N = 0), the boundary matrix H B becomes Hermitian, where the winding number is given by the conventional Berry phase.In other GBCs, L 1 and L 2 can form two disjoint contours in the complex plane.At the topological phase transition of the winding number W non−Bloch , where h + B (z 1 , z 2 ) or h − B (z 1 , z 2 ) vanishes, the topological boundary states absorb into the continuum of the bulk bands in the complex momentum plane.The merging of the two generalized momenta manifests as the exceptional point, which is signatured by the coalescence of distinct eigenstates [See Fig. 2 Bottom left].We refer to the appearance of the exceptional point during the topological phase transition of the GBZ as the exceptional phase transition.We point out that the observed topological boundary mode defies the standard phenomenology of the conventional topological insulators as our model only contains a single band (See I E for explicit studies in the general boundary conditions).
In general non-Hermitian one-dimensional systems, the one-dimensional non-Hermitian Hamiltonian can be written as, where α, β takes values from [1, p] with p being total number of internal degrees of freedom.t m describes the Generalization to many-body systems-We consider the one-dimensional array of the Kuramoto oscillators as a generalization to the many-body systems [11].(7) Here, the dynamics of i-th oscillator is described by the phase θ i (t).K R,L is the coupling with the right and left nearest-neighbor oscillator.δθ ij = θ i − θ j is the phase difference.The model admits the stable limit cycles with the constant phase gradient δθ.Due to the PBC, the phase gradient is quantized as, ω = N i=1 δθ i+1,i ∈ Z. (See supplementary materials for the exact solutions and stability analysis).
The linearized excitation near the limit-cycle phase is described by the Jacobian of Eq. ( 7), The Jacobian of the limit cycle phases emulates the HN model with the PBC [Fig.3(a)].The eigenvalue spectra of the Jacobian form a center-shifted circle in a complex plane comprised of negative real eigenvalues [Fig.3(e)].
The negative eigenvalue spectrum of the Jacobian characterizes the stability of the limit cycle in non-linear systems [12,13].The phase transitions between the limit cycles with different phase gradients (say ω and ω ) accompany the phase slip of 2π(ω−ω ), where the phase slipping sites can be considered as the boundary.During π 2 < δθ N < 3π 2 , the Jacobian matrix undergoes the topological phase transition of GBZ and it is characterized by non-trivial winding number.The corresponding topological boundary state is energetically separated from the negativeeigenvalue bulk states and forms the positive real eigenvalue [See Fig. 3(d)], which physically represents the onset of the limit-cycle instability during the phase transitions.
In addition to the negative eigenvalue bulk states, we additionally identify the Goldstone mode associated with the global U(1) symmetry (θ i → θ i +φ).Regardless of the phases of oscillators, the Jacobian matrix forms a zeroline sum (ZLS) matrix [14], which has zero-sum of elements in each row and column i.e.
At the topological phase transition, δθ N = π 2 , 3π 2 , the Jacobian exhibits the OBC.The GBZ of the bulk and the topological boundary state meet each other [Fig.

3(c)]
. Correspondingly, we observe the coalescence of the Goldstone mode with the normal modes, observed by the vanishing phase rigidity, are the left (right) eigenmodes corresponding to the first eigenstate.The vanishing phase rigidity can be understood as the consequence of the NHSE where the left eigenstates and the right eigenstates are localized at the opposite boundaries, causing the vanishing wavefunction overlap (See supplementary materials).This coalescence of the Goldstone mode signifies the non-equilibrium phase transition, which has been referred to as ETs [8].The EP of the Goldstone mode proliferates the temporal fluctuations associated with U(1) symmetry.The physical consequence is the dynamical restoration of the spontaneously broken U(1) rotation symmetry.
We propose that the overall chiral velocity serves as a macroscopic order parameter that characterizes such topological non-reciprocal phase transitions, The finite chiral motion occurs for each limit-cycles due to the non-reciprocity of the couplings.In the previous work [8], a similar chiral motion has been studied in the presence of random noise.Here, we note the crucial difference that the stability of the chiral motion does not require random noise due to the peculiar property coming from periodic boundary conditions.Instead, the presence of random noise destabilizes the chiral motion by inducing phase slips.During the phase slip, the ET proliferates the Goldstone mode excitations inducing the abrupt change in the chiral velocity [See Fig. 4].In the long time limit, the system eventually stabilizes into the limit cycle with the zero winding number.
Discussions -The non-trivial topology of GBZ we reveal here can be readily realized in designable material platforms such as active matter [15,16], metamaterials in the fields of optics and photonics [17], acoustics [18], robotics [19], and electric circuit networks [20].In a more general sense, the excitations living on systems with the translational symmetry-breaking defect can be the subject of our study.We have generalized our result to manybody collective behaviors by proposing non-reciprocally coupled oscillators.Various dynamical systems described by coupled non-linear oscillators such as active matter, coupled Josephson junctions, and XY spin model belong to this class of dynamical systems [15,16,21,22].In condensed matter systems, recently proposed studies of spin dynamics with the spin current can be an interesting platform to realize the ET of collective phenomena [23].

SUPPLEMENTARY MATERIALS I. NON-BLOCH BAND THEORY IN GENERALIZED BOUNDARY CONDITIONS A. Classifications of generalized boundary conditions
The non-Bloch band theory [5,10] describes the band structures of non-Hermitian systems in generalized boundary conditions (GBC), where the conventional real Bloch wave vector is replaced by the generalized wave vector lies in the complex plane.Here, we define the generalized Brillouin zone (GBZ) as the set of the complex wave vector that constitutes the eigenstates.As a representative example, we consider the one-dimensional Hantano-Nelson model [9] in the GBC.The full Hamiltonian operator can be decomposed as the bulk term ( Ĥbulk ) and the boundary term ( Ĥboundary ) as follows, where |i ≡ ĉ † i |0 is the localized state at i-th site.The matrix elements of the boundary term demonstrate various GBCs.For example, (t R , t L , 1 , N ) = (t R , t L , 0 , 0 ) corresponds to the periodic boundary condition (PBC), while (t R , t L , 1 , N ) = (0, 0, 0 , 0 ) corresponds to the open boundary condition (OBC).In a general setting, different combinations of onsite potentials e 1 , e N and hopping terms t R , t L fully describe the possible sets of the GBCs.Table I summarizes the definitions of the different types of boundary conditions.

Open boundary condition (OBC)
t R = 0, t L = 0 and In this section, we describe the detailed methods of the non-Bloch band theory [5,10] to find the GBZ.The eigenstate of Eq. (S1), |ψ = i ψ i |i , with the eigenvalue E is required to simultaneously satisfy the bulk and the boundary equations.We analyze the bulk and the boundary equations separately.
Bulk equation: The bulk equations are given as, To solve for the bulk equation, we take the single wave vector ansatz of the wave function, ψ n ∝ z n j , where the bulk equation transforms as, There are the two solutions of the ansatz for a given eigenvalue E, which constitute the paired momenta as, In addition, regardless of the eigenvalue E, the two solutions are related as, which provide the relations between the two paired momenta.In general, the linear superposition of the two ansatz forms general solutions of the bulk equations, Boundary equation: The boundary equations are given as, With the help of bulk equations, we simplify the boundary equations which take the following form, By plugging in the bulk ansatz in Eq. (S6), the boundary equation [Eq.(S9)] can be compactly rewritten as the form of the eigenvalue equation, where Eigenstate Solutions : We derive the non-zero solution of the boundary equation [Eq.(S10)].The non-zero solutions are obtained by solving for the equation, det[H B ] = 0, which takes the explicit form as follows, By solving the above boundary equation with the complementary condition [Eq.(S5)], we derive the spectrum of the generalized Brillouin zones.The paired momenta satisfying Eq. (S5) can be expressed in the polar coordinates as follows, where φ ∈ C is a complex number.When the solutions of φ are real number, the GBZs of the two paired momenta coincides.Accordingly, Eq. (S11) can be rewritten as the following form, where In practice, we solve Eq. (S13) to find the solutions of the GBZ.
It is important to point out that, when φ = 0 or π, the two paired momenta are equal, z 1 = z 2 = r or −r.Accordingly, Eq. (S10) leads to the condition c 1 = −c 2 , which indicates the vanishing wave function, ψ n = 0 (i.e.non-physical solution).However, this conclusion is not correct.For the equal solutions z 1 = z 2 = z, the general solution of the bulk equation is given by the form of Consequently, the boundary equations will be modified which are expressed as follows, where ψ n = (c 1 + c 2 n)z n is the general solution of the eigenstates.

C. Representative examples of Generalized Brillouin zones
In this section, we present the representative solutions of the GBZs in GBCs.We determine the eigenvalues and corresponding eigenstates for a finite-size system in the generalized boundary condition.Using the paired momenta described by Eq. (S12), the eigenvalues E and corresponding wave functions can be expressed in terms of generalized momenta as follows, Fig. S2 compares the calculated GBZs with the result of numerical exact diagonalizations of the Hamiltonian.We find that the two methods agree well.We detail the results of some of the boundary conditions as follows: In the right panel of (d), the real eigenvalues which are marked with circles represent the localized boundary states due to GBC.These boundary states are characterized by the disjoint generalized momenta outside or inside GBZs which are also marked with circles in the left panel of (d).Parameters:

Periodic boundary condition
In this section, we describe the GBZs and eigenvalue spectra for periodic boundary condition (PBC) defined with which forms a closed curve (an ellipse) in the complex eigenvalue plane and agrees with that of numerical diagonalization (See Fig. S2(a)).The boundary equations together with GBZs determine the wave functions which take the following form, which form an open segment in the complex eigenvalue plane and collapse on the real eigenvalue axis for the present model (See Fig. S2 (c)).Also from boundary equations, we have c 1 = −c 2 and thus the wave functions become, It is clear from the above expressions that all the eigenstates localize at one of the boundaries of the system depending on whether r < 1 or r > 1.This anomalous localization phenomenon is known as the NHSE [1][2][3].The eigenvalue spectra and corresponding eigenstates differ from that of PBC, which shows the extreme sensitivity of non-Hermitian systems to the boundary conditions [1][2][3].

Generalized boundary conditions-(III)
In this section, we discuss GBZs and eigenvalue spectra for one of the generalized boundary conditions defined as , and 1 = N = 0 .In this case, Eq. ( S13) is written as which possesses N different physical solutions φ m ∈ C. We note that the solutions of Intermediate boundary condition-(i) have been discussed in Ref. [10].Similar to Ref. [10], Fig. S2 (b) clearly shows that the two GBZs start approaching each other for smaller values of t L /t L and t R /t R which corresponds to the flattening of the eigenvalue spectra (i.e.approaching to the OBC).

Generalized boundary conditions-(IV)
In this section, we consider the generalized boundary condition which is defined as

D. Topological classifications of GBZ
For each paired momenta (z 1 , z 2 ) in GBZ, the real space wave function is described as where the adiabatic change of the wave function can be described in terms of the effective spinor |χ = (c 1 , c 2 ) T .In this section, we show that the non-trivial topology of the GBZ manifests as the emergence of the topological boundary states.

Symmetry of boundary equation
Since the effective spinor satisfies the boundary equations H B (z i )|χ = 0, the symmetry group of the boundary matrix H B characterizes the topology of GBZ.The boundary equations can be rewritten as the eigenvalue problem, where In general, the non-Hermitian effective boundary matrix satisfies the following chiral symmetry by construction, {σ z , HB } = 0. (S27) In the presence of the chiral symmetry, for any right (left) eigenstates |χ R(L) with complex eigenvalue E (E * ), there always exists right (left) chiral eigenstates | χR(L) = σ z |χ R(L) with eigenvalues −E (−E * ) which satisfy the bi-orthogonality relations: In addition, for GBC-(I) (See Table I), HB becomes a Hermitian matrix since the meromorphic function h + B is unimodular (|h + B | = 1) together with det[H B ] = 0. Whereas, for other possible sets of generalized boundary conditions (See Table I), the meromorphic function h + B is not unimodular (|h + B | = 1), and hence HB remains a non-Hermitian matrix.

Winding number of GBZ
The winding number of the boundary matrix can be defined using Cauchy's argument principle, which relates the number of zeros and poles of a meromorphic function, h ± B (z 1 , z 2 ), to the contour integration of its logarithmic derivative as follows: where the integration contour L 1 × L 2 indicates the continuous contour of the GBZs satisfying goes along the GBZ.Z(P ) represents the total number of zeros(poles) inside GBZs counting the multiplicities.Fig. S3(a)-(c) shows exemplifies the GBZs with different winding numbers for GBC-(I) (See Sec.I E 2 for more details).In the topologically trivial GBZ [W = 0, Fig. S3(a)], the interior of the GBZ encloses an equal number of poles and zeros of h + B .As the system undergoes the topological phase transition [Fig.S3(b)], the zero of h + B touches the GBZ.Finally, the GBZ encircles a pair of poles, which manifest as the non-trivial value of the winding number [Fig.S3(c)].
In principle, the winding numbers W + and W − are unequal.Namely, we can separately define the winding numbers W and W as, As a result, W can only have non-zero numbers.In the next section, we discuss that there is an ambiguity in defining the topological invariant using Eq.(S29), and accordingly, we will re-define the topological invariant to characterize the topological phase transition of GBZ.

Ambiguity in defining winding numbers and non-Bloch topological invariant
Eq. (S28) possesses an ambiguity in defining the winding number.The boundary equations [Eq.(S10)] obtained with the ansatz ψ n = c 1 z n 1 + c 2 z n 2 are not unique.Since we can also write the ansatz of the wave function as . Accordingly, the modified boundary equations (Eq.(S10)) takes the following form, B ] = 0. Now the modified effective boundary equations are written as the following eigenvalue problem . Topological phase transition of GBZ in GBC-(I).Top panel: Vector plot of the meromorphic function h − B together with GBZ in the complex plane.We also show the schematics for the positions of poles and zeros of h + B together with the contour of GBZ.For N < tLr (See top panel (a)), there are an equal number of poles and zeros of h − B inside the GBZ.Accordingly, the non-topological invariant (defined in Eq. (S36)) vanishes.Hence the parameter regime ( N < tLr) defines the topologically trivial phase.In the parameter regime N > tLr (See top panel (c)), there are a pair of poles inside the GBZ.Accordingly, the non-Bloch topological invariant (defined in Eq. (S36)) takes a non-trivial value, and the parameter regime ( N < tLr) defines the topologically non-trivial phase.In this case, the non-zero poles and zeros of h − B coincide with the generalized momenta (disjoint from GBZ and marked by circles) characterizing the topological boundary state (See top panel (c)).The topological phase transition of GBZ occurs at N = tLr, where the disjoint generalized momenta, and the contour of the GBZ merge which characterizes the emergence of exceptional points in the eigenvalue spectra (See top panel (b)).Bottom panel: Non-Bloch topological invariant (See bottom panel (left)) takes non-trivial values at N = ±tLr characterizing the topological phase transitions of GBZs.The emergence of exceptional points at the topological phase transition points is also confirmed by the vanishing of the phase rigidity R1 and RN (See bottom panel (right)).Parameters: . Hence, the winding numbers (defined in Eq. (S28)) are also modified and take the following forms Accordingly, the topological invariant defined in Eq. (S29) will also be modified and is given by, It is evident from the above equations that the winding numbers and the topological invariant depend on m.We elaborate above arguments with the example of the open boundary condition in the following.In the OBC, the modified boundary equations (Eq.S30) are given by, Accordingly, the physical solutions (z 1 , z 2 ) are determined by the solving the following equation which is independent of m.The solutions of the above equations are of the form (z 1 , z 2 ) = (re iφ k , re −iφ k ) with φ k = kπ/(N + 1), k ∈ {1, 2, ..., N }.Therefore, the generalized Brillouin zones and eigenvalue spectra do not depend on m, however, c 1 , c 2 depend on it.Accordingly, the winding numbers associated with the meromorphic functions h + B , h − B depend on m, which are given by, W (m) ± = ∓2m.For instance, we can list the wave functions together with winding numbers for different values of m as follows: (i) for m = 0: + = −2l, and so on.
Despite of the ambiguity, the topological phase transition always accompanies a change in the total number of zeros and poles of h + B , h − B inside GBZs.Accordingly, we define our topological invariant as follows: We note here that the topological invariant W L can also be obtained by interchanging the parameters (t R , t L , 1 , N ) to (t L , t R , N , 1 ) in the eigenvalue equation or vice versa.Hence, we can think of the non-Bloch topological invariant W R as the topological invariant defined for the right bulk which is defined by the parameters (t R , t L , 1 , N ), and the non-Bloch topological invariant W L defined for left bulk by interchanging the parameters (t R , t L , 1 , N ) to (t L , t R , N , 1 ).We find that the non-Bloch topological invariant W non−Bloch (defined in Eq. (S36)) characterizes the topological phase transition of GBZ for different generalized boundary conditions which we discuss in some detail in the upcoming sections.

E. Exceptional and topological phase transitions of GBZ
In contrast to the degeneracy of eigenenergies, the degeneracy of two generalized momenta directly indicates the coalescence of the eigenstates, which manifests as an exceptional point in the spectra.The exceptional point is characterized by vanishing of phase rigidity, (See Fig. S3) here |Ψ R i is the right eigenstates of Ĥ (Eq.(S1)) with eigenvalue E i and |Ψ L j ) is the left eigenstates of Ĥ (Eq.(S1)) with eigenvalues E * i which satisfy bi-orthogonality relations i.e.Ψ R i |Ψ L j = δ ij .We calculate the phase rigidity of all the eigenstates by sorting them according to the real parts of the eigenvalues.We elaborate on these findings in upcoming sections when analyzing different types of GBC.In this section, we present an analytical study of the topological phase transition of GBZ.First, we consider the GBC-(I) (See Table I) with t R = t L = 0, 1 = 0 , N = 0 .In this case, the boundary equations become, The paired generalized momenta (z 1 , z 2 ) are obtained by solving the following equation, We describe the topologically trivial and non-trivial phases in the following manner: Topologically trivial phase: In the parameter regime N < t L r, two GBZs merge with each other similar to that of the OBC (See top panel (c) of Fig. S3).This parameter regime N < t L r describes topologically trivial GBZ since there is an equal number of zeros and poles of h − B inside GBZs, and hence the winding numbers, and non-Bloch topological invariant vanish (See Secs.I D 2 and I D 3).For N << t L r, the bulk states are described by the generalized momenta (z 1 , z 2 ) = (re iφ k , re −iφ k ) with φ k = kπ/(N + 1), k ∈ Z.Moreover, the effective spinor, and wave functions are given by, Topologically non-trivial phase: For N > t L r, again two GBZs merge with each other similarly to the topologically trivial phase, however, there exists one pair of the real solution (z 1 , z 2 ) = (t R / N , N /t L ) (See top panel (c) of Fig. S3) which describes the topological boundary state (ψ n ∝ z n 2 ) localized at the boundary of the chain.In this parameter regime (See top panel (c) of Fig. S3), the non-zero pole of h − B together with the pole at zero lies inside GBZs.Accordingly, the winding numbers and the non-Bloch topological invariant become (See Secs.I D 2 and I D 3) Therefore, the parameter regime N > t L r describes the topologically non-trivial phase which supports the emergence of the topological boundary state described by the disjoint generalized momenta (z 1 , z 2 ) = (t R / N , N /t L ).Hence the parameter regime N > t L r describes the topologically non-trivial phase.Moreover, for N >> t L r, the bulk equations are characterized by the generalized momenta (z In this case, the expanded spinor, and wave functions are given by, Similar to the previous section (Sec.(I E 2)), we also find that the non-Bloch topological invariant (the bottom left panel of Fig. S3) changes and takes non-trivial value at N = ±t L r characterizing the topological phase transition of GBZs.At this transition point the disjoint generalized momenta (z 1 , z 2 ) = (t R / N , N /t L ) describing the boundary state touches the GBZs (See top panel (b) of Fig. S3) which is the exceptional point in the eigenvalue spectra of the Ĥ.The appearance of the exceptional points is confirmed by the vanishing of the phase rigidity R 1 , R N in the large N −limit (see bottom right panel of Fig. S3).

Analytical examples (ii): Generalized boundary conditions-(I)
In this section, we present another analytical study of the topological phase transition of GBZ.Now, we consider the GBC-(I) with t R = t L = 0, 1 = 0 , N = 0 .In this case, the boundary equations become, The paired generalized momenta (z 1 , z 2 ) are obtained by solving the equation which is independent of m and may not be solved analytically.We can understand the topological phase transitions of GBZs by analyzing the boundary equations in the following.
Topologically trivial phase: For 1 > t L r, we find that the two GBZs merge with each other similar to that of OBC.However, in this parameter regime, there is an equal number of poles and zeros of the meromorphic function h + B inside GBZs.Accordingly, the winding numbers, and the non-Bloch topological invariant vanish (See Secs.I D 2 and I D 3).Therefore, the parameter regime 1 < t L r describes the topologically trivial phase.In the limiting case ( 1 << t L r), the approximate solutions of the boundary equations are (z 1 , z 2 ) = (re iφ k , re −iφ k ) with φ k = kπ/(N + 1), k ∈ Z.Moreover, the effective spinor and wave functions are given by, Topologically non-trivial phase: For 1 > t L r, we find that the two GBZs merge with each other similar to the previous case, and in addition to that, the boundary equations have one pair of the real solution (z 1 , z 2 ) = (t R / 1 , r 2  1 /t R ) in the large N limit.This pair of disjoint generalized momenta describes the topological boundary state (ψ n ∝ z n 2 ) localized at the boundary of the chain.In this case, there are a pair of poles of h + B inside the GBZs.Accordingly, the winding numbers and the non-Bloch topological invariant become (See Secs.I D 2 and I D 3) Therefore, the parameter regime 1 > t L r describes the topologically non-trivial phase which supports the emergence of the topological boundary state described by the disjoint generalized momenta (z 1 , z 2 ) = (t R / 1 , r 2 1 /t R ).For 1 t L r, the approximate solutions of the boundary equations are given by (z 1 , z 2 ) = (re iφ k , re −iφ k ) with φ k = kπ/(N ), k ∈ Z.In this case, the effective spinor, and wave functions are given by, The topological phase transitions of GBZs occur at the exceptional points 1 = ±t L r where the disjoint generalized momenta (z 1 , z 2 ) = (t R / 1 , r 2  1 /t R ) describing the topological boundary state touch the GBZs (See Fig. 2).The appearance of the exceptional points is also confirmed by the vanishing of the phase rigidity R 1 , R N in the large N −limit (See Fig. 2).Moreover, the non-zero zero and pole of h + B coincide with the disjoint generalized momenta in the topologically non-trivial phase and merge with the GBZs at the topological phase transition point.

Example (iii): GBC-(I)
In this section, we discuss the exceptional and topological phase transitions of GBZ for the GBC-(I) which is defined with t L = t R = 0, and 1 ( N ) is not equal to 0 at the same time but can have any value in general (See Table I).This type of boundary condition can be understood as an open boundary condition with different onsite potentials at the boundaries of the chain.In this case, the boundary equations are, Generalized Brillouin zones and boundary states: Solving the above boundary equations, we get N , N − 1, or N − 2 (depending on the strengths of 1 and N ) pairs of generalized momenta with |z 1 | = |z 2 | = r.Hence, two generalized Brillouin zones merge similarly to OBC (See Figs.S4(i)-(ii)).The wave functions corresponding to the bulk states can be written as , where |z 1 | = |z 2 | = r.These states describe localized states at one of the boundaries of the system (skin states).
The remaining one pair (or two pairs) of disjoint generalized momenta (See the evolution of GBZs and the emergence of disjoint generalized momenta in Figs.S4(i)-(ii)) lie on the real momenta axis, i.e. z 1 , z 2 ∈ R, (for r < 1).Accordingly, the wave functions corresponding to these disjoint momenta become ψ n ∝ z n 2 and describe the bound states localized at one of the system's boundaries.
Exceptional phase transitions of GBZ: We find exceptional points in the eigenvalue spectra of the Hamiltonian by varying 1 = t R , N = t L simultaneously which are dictated by the vanishing phase rigidity (defined in Eq. (S37)) R 1 , R 2 for positive values of , and R N , R N −1 for negative values of , respectively (See Fig. S4(iv)).The dips in the phase rigidity depend on the degree of non-reciprocity (r 2 = t R /t L ) and system size (N ) which become sharper with an increase in the system size.Importantly, we observe that at exceptional points, bound states described by disjoint generalized momenta in the complex plane emerge out of the bulk continuum (See Figs.S4(i)(c)-(e)).The numerical finding of exceptional points is confirmed by touching of disjoint generalized momenta to the contour of the generalized Brillouin zone (See Fig. S4(i)(c)-(e) and Fig. S4(ii)(e)).Therefore, exceptional points characterizing the emergence of bound states from the bulk continuum define the exceptional transitions of GBZ.
Topological phase transitions of GBZ: We find that boundary states in response to the GBC are characterized by the non-trivial values of the non-Bloch topological invariant W non−Bloch .The Fig. S4(iii) (left) dictates the change of the non-Bloch topological invariant at the exceptional points defined by the touching of disjoint generalized momenta to the contour of GBZ (See vanishing of phase rigidity R 1 , R 2 , R N −1 , and R N in Fig. S4(iv) (left)) which marks topological phase transitions of the GBZ.Therefore, topological phase transitions of GBZ accompany the exceptional points.
We also find that in GBC-(I), the effective boundary matrix HB (See Sec.I D 1) becomes Hermitian together with the chiral symmetry.Accordingly, in the topologically non-trivial phase, the effective spinor (c 1 , c 2 ) T projected on the Bloch sphere always lies on the equator for all generalized momenta of the bulk (See figures (right) of panels (iii) and (iv) of Fig. S4).Therefore, the associated Berry phase is quantized and serves the purpose of characterizing topological phase transitions.In this case, the non-Bloch topological invariant coincides with the winding number related to the conventional Berry phase associated with the effective spinor (c 1 , c 2 ) T .

Example (iv): GBC-(II)
In this section, we discuss the exceptional and topological phase transitions of GBZ in GBC-(II) which is defined as t L = t L , t R = t R , and 1 = N = 0 .This type of boundary condition can be understood as the periodic boundary conditions with different onsite potentials at two neighboring sites anywhere in the chain.In this case, the boundary equations are, Generalized Brillouin zones and boundary states: After solving the above boundary equations, we get N , N − 1 or N − 2 (depending on the strengths of 1 and N ) pairs of generalized momenta with |z 1 | = |z 2 | which form two GBZs.We find that similar to the case of PBC, in this case also, the two generalized Brillouin zones differ from each other.The remaining one pair (or two pairs) of disjoint generalized momenta lie on the real momenta axis, i.e. z 1 , z 2 ∈ R, (for r > 0).In this case, also, the wave functions corresponding to these disjoint generalized momenta become ψ n ∝ z n 2 which describe the bound states localized at the boundaries of the system.The evolution of GBZs and the emergence of disjoint generalized momenta by changing the onsite potentials are shown in Figs.S5(i)-(ii).We also find that at very large values of the onsite potentials ( 1 or N or both), the two GBZs merge, and the bulk eigenvalue spectra collapse on the real eigenvalue axis (See Fig. S6).
Exceptional phase transition of GBZ: We find two exceptional points in the eigenvalue spectra of the Hamiltonian by varying 1 , N simultaneously which are dictated by the vanishing phase rigidity R 1 , R 2 for positive values of 1 , and R N , R N −1 for negative values of 1 (See the left panel of Fig. S5(iv)).Importantly, we observe that after these exceptional transition points, there appear bound states localized at the boundaries of the system.The numerical finding of exceptional points due to GBC-(II) is confirmed by touching of disjoint generalized momenta to the generalized Brillouin zones (See Figs.S5(i)-(ii), and left panel of Fig. S5(iv)).Therefore, the exceptional points characterizing the emergence of bound states outside the bulk eigenvalue spectra define the exceptional transitions of the GBZ.Topological phase transition of GBZ: Similar to the GBC-(I), in this case also, the boundary states arising at the exceptional points in the eigenvalue spectra are characterized by non-trivial values of non-Bloch topological invariant (W non−Bloch ).Moreover, the topological phase transition of the GBZ occurs at these exceptional points which are marked by the quantized change in the non-Bloch topological invariant (See Fig. S5) (iii) (left)).The emergence of exceptional points at the topological phase transition points is also confirmed by the vanishing of phase rigidity R 1 , R 2 , R N −1 and R N (See left panels of Fig. S5(iii)-(iv)).
We also find that the conventional Berry phase associated with the effective spinor will not be quantized in this case because the effective boundary matrix HB is non-Hermitian (See Sec.I D 1).We determine the number of bound states (n b ) in the ( 1 , N )-parameter space ( 1 , N )-plane (See right panel of Fig. S5 (iii)).The non-Bloch topological invariant show a similar phase diagram in the ( 1 , N )-plane (See right panels of Figs.S5(iii)-(iv)) which characterizes the emergence of the boundary states.

Example (v): GBC-(III) and (IV)
In this section, we consider the generalized boundary condition by combining the ones described by (iii) and (iv) type of GBC mentioned in Table I: t L = t L , t R = t R , and 1 = N = 0 .In this case, the boundary equations are, where they merge with each other.The remaining one pair (or two pairs) of disjoint generalized momenta lie on the real momenta axis, i.e. z 1 , z 2 ∈ R, (for r > 0).In this case, also, the wave functions corresponding to these disjoint generalized momenta become ψ n ∝ z n 2 which describe the bound states localized at the boundaries of the system.The evolution of GBZs (L 1 , L 2 ) and the emergence of disjoint generalized momenta by changing the boundary hopping terms are shown in Fig. S7(i).Also, at very large values of the boundary hopping terms (both t R and t L ), the two GBZs merge, and the bulk eigenvalue spectra collapse on the real eigenvalue axis (See Fig. S8).
Exceptional transition of GBZ: We observe exceptional points by varying t R = t R , t L = t L simultaneously which are dictated by the vanishing phase rigidity R 1 and R N for positive and negative values of t R , respectively (See Fig. S7(iii)).Importantly, we find that after the appearance of exceptional points (for non-zero t R , t L ), there emerge bound states localized at the boundaries of the system.The boundary states do not arise for arbitrary small values of t L − t L or t R − t R , however, arise for |t R − rt R | > 0. These boundary states are characterized by the disjoint generalized momenta which touch the contours of GBZs at t R = rt R (See Fig. S7 (i)(c)-(d)).The touching of generalized momenta manifests as the exceptional points in the eigenvalue spectra (See Fig. S7(i)(c), (iii)).Thus, the numerical finding of exceptional points due to the GBC is confirmed by touching of disjoint generalized momenta to the GBZ (See Fig. S7 (i)(c), (iii)).Therefore, the emergence of the exceptional points in the eigenvalue spectra defines the exceptional transition of the GBZ.We note here that at t R = t L = 0, there are additional exceptional points corresponding to the coalescence of all the eigenvalues and eigenstates in the thermodynamic limit which are dictated by the vanishing phase rigidity for all the eigenstates and are confirmed by the merging of two GBZs (See Fig. S7 (i)(a)).At this point, the HN model realizes the OBC and we observe NHSE.
Topological phase transition of GBZ: Similar to the GBC-(I) and GBC-(II), in this case also, we find that the appearance of the boundary states after the exceptional points are marked by the non-trivial values of the non-Bloch topological invariant (W non−Bloch ) (See Fig. S7(ii)).Fig. S7(ii) dictates the quantized change in the non-Bloch topological invariant at the exceptional points which are marked by the vanishing of phase rigidity R 1 and R N (See Fig. S7(iii)).Therefore, the non-Bloch topological invariant characterizes the topological phase transition of GBZ.
In this case, the effective boundary matrix HB is non-Hermitian (See Sec.I D 1), and accordingly, the Berry phase associated with the effective spinor will not be quantized.We also determine the number of bound states (n b ) in the (t R , t L )-plane (See Fig. S5(iii) (right)).The non-Bloch topological invariant (W non−Bloch ) shows a similar phase diagram in the ( , )-plane (See Figs.S5(iii)-(iv) (right)) which characterizes the topological phase transition of GBZ where boundary states emerge out of bulk continuum.
It is interesting to compare these results with 1D Hermitian systems (t R = t L =⇒ r = 1) where the bound states appear for arbitrary small values of t L − t L .There is no exceptional point in this case, and the bound states arise at certain points in the parameter space when only eigenvalues (at least two) touch each other which is also known as the band-touching point or diabolic point.

II. GENERALIZATION TO MANY BODY DYNAMICAL SYSTEMS
To generalize our results to many body systems, we consider the Kuramoto model of coupled oscillators [11,24,25], which describes the synchronization phenomenon [26] in various physical systems including Josephson junction arrays [21], laser arrays [27][28][29].We consider the Kuramoto model of locally coupled oscillators (LCKM: locally coupled Kuramoto model) on a 1D lattice of a finite number of oscillators N with periodic boundary conditions.The Kuramoto model of locally coupled oscillators is given by [11,24,25], which describes the time evolution of the phase θ i (t) of oscillator i with natural frequency ω i for i ∈ {1, 2, ..., N } with δθ j,i (t) = θ j (t) − θ i (t) being the phase difference between the oscillators i, j ∈ {1, 2, ..., N }.K R and K L denote the strengths of forward and backward nearest neighbor couplings, and η i (t) denotes the random white noise of strength D with η i (t) = 0 and η i (t)η j (t ) = δ ij δ(t − t ).Here ... denotes the average over independent noise realizations.

A. LCKM for two oscillators
Before discussing the LCKM on a 1D lattice with a finite number of oscillators, we review the LCKM of two oscillators first.We consider the following dynamical equations for two oscillators θ 1 and θ 2 , Stationary synchronized solutions: For reciprocal coupling K R = K L = K, stationary synchronized solutions exist δθ 1,2 ≡ δθ = 0 for aligned phase (in phase synchronized solutions) and δθ = π for anti-aligned phase (out of phase synchronized solutions) [32].Similarly, for non-reciprocal couplings K R = K L , the oscillators align with each other for (K R + K L ) > 0 or anti-align with each other for (K R + K L ) < 0. Therefore, we have two possible synchronized phases which can be understood as follows: (i) for (K R + K L ) > 0, one of the oscillators is going to chase down the second one and they will eventually align (in phase synchronized solutions) (See Fig. S9).
(ii) for (K R + K L ) < 0, again one of the oscillators will win and they will eventually anti-align (out of phase synchronized solutions) (See Fig. S9).

FIG. S9. Synchronized solutions for LCKM for two oscillators (See text).
Time-dependent synchronized solutions: Interestingly, for exact anti-symmetric coupling i.e K R = −K L , the dynamical equations simplifies in the following manner, As a result, two oscillators move with finite and common drift velocity (K R − K L ) sin(δθ(0)) keeping constant phase difference δθ(0) which depends on the initial conditions.This type of synchronized stable phase has been called the chiral phase [8].Indeed, the system of two oscillators dynamically restores the U(1) (θ i → θ i + θ 0 ) global rotation symmetry.For K R = −K L , the system collapses to one of the two stable phases (in phase and out of phase).

B. LCKM for N oscillators with reciprocal couplings
In this section, we discuss the Kuramoto model of locally coupled oscillators.We start with the stationary solutions of LCKM with reciprocal couplings [12,13,32,33].

Stationary solutions
Consider the LCKM with reciprocal couplings K R = K L = K.The stationary solutions θ = ( θ1 , ..., θN ) T for the dynamical equations (S51) are given by, which immediately implies that LCKM with reciprocal couplings has two types of stationary solutions θ = ( θ1 , ..., θN ) T : (i) δ θi+1,i − δ θi,i−1 = 0 (mod 2π), ∀ i ∈ {1, 2, ..., N }, It is important to point out that the above stationary solutions described by Eqs.(S54) and (S55) will only hold with reference to the periodic boundary condition (i.e.θ N +i = θ i : the first oscillator is coupled to the last one in the same fashion as that of any θ i , i ∈ {2, 3, ..., N − 1}).In simple words, the stationary solutions of LCKM are those where the phase differences between all the pairs of oscillators are the same.This leads to two different synchronized states where the neighboring oscillators are aligned with |δ θi+1,i | < π/2 and anti-aligned with |δ θi+1,i | > π/2 depending on the coupling K.These stationary solutions also hold if we include the next nearest neighbor coupling (or any finite range coupling) in Eq. (S51).
Exact form of synchronized solutions: We can determine the exact form of the synchronized solutions analytically.For that, we write the synchronized solutions in the following general form [13,33]: here ω s is the common frequency and θ = ( θ1 , ..., θN ) T denote phases for the oscillators corresponding to the given synchronized solution.It has been shown in [13] that the synchronized solutions of LCKM with reciprocal couplings take the following form: θ with respect to constant nearest neighbor phase differences δ θ for twisted states with K > 0 and K < 0, respectively (See Sec.II B 2 for more details on the local stability of twisted states).Figures in the right of the top and bottom panels show the vector plots of stable twisted states (a < 0) with different winding numbers (Wp) for K > 0 and for K < 0, respectively.summation, the above equations immediately end up to Hence, the synchronized solutions are characterized by the integer 'p' which for δ θi+1,i << 1 is the same as the topological winding number W p ∈ Z.In other words, a unique synchronized solution θ(p) is defined via topological winding number: Therefore, the synchronized solutions have topological protection.Moreover, for small and equal phase differences, δ θi+1,i = δ θ = 2πW p /N, ∀ i, phases of individual oscillators are given by θ(p

Here θ(p)
0 is an arbitrary uniform phase shift.These synchronized solutions take the form of uniformly twisted waves with the number of complete phase twists equal to the value of the topological winding number.Accordingly, these synchronized solutions are also called W p -twisted states with topological charge or the number of complete phase twist W p ∈ Z [12,13] (See vector plots of stable twisted states in the top and bottom panels (right) of Fig. S10).We further divide different synchronized phases into the following two parts: (i) State of complete phase synchronization: In this case, starting from random initial conditions, after a sufficiently long time the phase differences among the oscillators vanishes i.e.
Therefore, this synchronized state of oscillators is called state of complete phase synchronization which is also known as the sync.In this case, the synchronized solutions are characterized by θ(p) with the topological winding number W p = 0 [12,13].
(ii) Phase-locked state: In this case, depending on different random initial conditions the phase differences among the oscillators approach to constant values i.e. |θ i (t → ∞) − θ j (t → ∞)| → const., ∀ i = j.Therefore, all the oscillators approach their different steady state value such that dθi dt = 0 ∀ i, and the oscillators are said to be in the phase-locked state.Moreover, for small and equal constant phase differences, the phase-locked states are known as W p -twisted states with the number of complete phase twist W p ∈ Z [12,13].This essentially means that stationary synchronized solutions θ(p) have topological protection identified by topological winding number W p ∈ Z − {0}.
2. Stability analysis of the synchronized solutions a. Jacobian of linear stability: One can analyze the linear stability of the above-synchronized solutions by perturbing the system by a small amount such that θ i (t) → θi + φ i (t), with |φ i (t)| 1.Then the evolution equation for φ i (t) is given by, We can also write the above equations in vector form as dφ dt ≈ Ĵ φ, with Ĵ being a N × N matrix (also known as the Jacobian or perturbation matrix) and φ = (φ 1 , ..., φ N ) T .
We note here the remarkable similarity between the Jacobian Ĵ and the one-dimensional tight-binding Hamiltonian with the generalized boundary conditions for 1D lattice systems with reciprocal hopping.We will come to this point later.We summarize the relationship between the Jacobian and local stability of synchronized solutions in the following manner [12,13]: We know that for the first-order ordinary differential equations, any initial conditions non-identical to a stable solution are perturbations for that stable solution.Moreover, the negative eigenvalues of the Jacobian or the perturbation matrix indicate the speeds of return to that stable solution for small perturbations.Therefore the initial conditions very close to the stable solution with negative eigenvalues of the Jacobian will approach the stable solution eventually.We also know that the basins of attraction are connected sets containing stable solutions.Based on these arguments, we assume the correspondence between the volume of the basin of attraction and the stability of the solution, which means that the solutions whose all eigenvalues are strongly negative are more stable and have a larger basin of attraction.Similarly, the solutions with maximal negative eigenvalues closer to zero are less stable and have smaller basins of attraction [12,13].We analyze the eigenvalue spectra of the Jacobian corresponding to different synchronized solutions in the following sections.

b. Local stability analyses:
The stability of the synchronized solutions of LCKM with reciprocal couplings has been studied analytically with the help of Gershgorin circle theorem [34] which states that every eigenvalue of a general N × N complex matrix A = [a ij ] i,j∈{1,2,...,N } with complex entries a ij lies within at least one of the closed Gershgorin discs D(a ii , R i ) ⊆ C, ∀ i ∈ {1, 2, ..., N } centered at a ii having radius R i = j =i |a ij | ∀ i ∈ {1, 2, ..., N }.For the Jacobian matrix (See Eq. (S61)), the center and radius of Gershgorin discs are as follows: the decrease in the maximal eigenvalue with the increase in the system size indicates that arbitrary small unbounded noise will destabilize the stable solutions with finite winding numbers in the thermodynamic limit.Through our numerical simulations, we observe that starting with the random initial conditions, the probability of finding stable synchronized solutions with higher finding numbers is low as compared to that of smaller winding numbers.We find that the stationary solution with the topological winding number W p = 0 is globally stable.Whereas other stable synchronized solutions with W p = 0 are only locally stable i.e. they have finite basins of attraction which again support the linear stability analyses.Hence the LCKM with reciprocal coupling hosts stable synchronized solutions for finite sizes of the system.These stable solutions can be found explicitly depending on the initial conditions and their number is proportional to the finite size of the system.

Stability of synchronized solutions in the presence of random white noise and disorder in frequency
In this section, we study the stability of the synchronized solutions in the presence of random white noise and also disorder in frequency.Firstly, we describe the numerical method which we have utilized in our numerical simulations.
a. Details about numerical simulations: In order to incorporate the white noise term correctly, we implement Euler-Maruyama's method and also Heun's method in our numerical simulation to integrate the stochastic differential equation [38].These algorithms are important to use because the standard methods like Euler, Runge-Kutta, etc. assume that the function governing the dynamics are well behaved and differentiable.However the white noise is not a differential function even for a single realization, hence its correct incorporation is necessary.The basic idea to implement these algorithms is to integrate the dynamical equation: although the values of the white noise and its derivatives are not well defined, the first integral of the white noise is, and it is the Wiener process which is a continuous function [38,39].In our case, using the Euler-Maruyama scheme, the iterative solutions for the dynamical equation are given by, here dt is the time interval, θ(t) = (θ 1 , θ 2 , ..., θ N ) T is a column vector describing the phases of all oscillators, f (θ(t), t) = ω + S T (t)×(K • S(t)) with S(t) = (cos(θ(t)), sin(θ(t))).K represents the coupling matrix describing the nearest neighbor couplings of oscillators on a lattice.ω = (ω 1 , ω 2 , ..., ω N ) T describes the natural frequencies for all oscillators and can take random values in the presence of frequency disorder.η is a vector of random numbers corresponding to the Gaussian white noise.It is important to note here that if the random number generated has a different variance than unity then the above formula should modify accordingly [39].The Euler-Maruyama method has a strong order of convergence of 0.5.We use a better method called Heun's method which has strong order of convergence is 1.0 and will produce better trajectories for our stochastic differential equation.Again, using Heun's scheme, the iterative solution for our stochastic equation becomes, here f 1 = f (θ(t), t), f 2 = f (θ(t) + dtf 1 + Dη √ dt, t + dt).The higher schemes for numerical integration become much more complicated and for comparison of numerical results for different schemes for Langevin-type stochastic differential equations has been provided in Ref. [38].Here we consider Heun's scheme for our numerical simulation to implement the Gaussian white noise.Also, we consider the fourth-order Runge-Kutta method in the absence of the random white noise together with Heun's scheme.For small enough time intervals, both schemes give similar results in the absence of noise.b.Noise induced non-equilibrium phase transitions: We have discussed in the previous section that the twisted states are characterized by distinct winding numbers W p ∈ Z.Therefore, the transition between two different winding number sectors (say W p and W p +1) is only possible when one of the phases (say θ N ) changes from 2πW p to 2π(W p +1) which implies that one of the phase differences (say δθ 1,N = θ 1 − θ N ) undergoes a change from δθ 1,N to δθ 1,N + 2π.Therefore, the transition between different winding number sectors representing the twisted states can occur through the abrupt change of one of the phase differences by 2π.This relatively rapid change of the phase difference by 2π is called phase slip (See Ref. [26] for more details).
In the presence of random white noise, the phase slips occur resulting in noise-induced phase transitions to the different twisted states (See left panel of Fig. S12).After a long enough time, the oscillators finally end up in the complete phase synchronized state with W p = 0 and the phase slips occur around this state W p = 0. We note here that for a given noise realization, both phase transitions from W p -twisted state to the (W p + 1)-twisted state and also W p -state to the (W p − 1)-twisted state can occur.These phase transitions are known as noise-induced non-equilibrium phase transitions.In addition, for twisted states with larger winding numbers, the phase slips occur rather quickly for small strengths of white noise due to the larger size of the basins of attraction.Whereas for twisted states with small winding numbers, the phase slips occur after longer times (See Fig. S12 (left)).We also observe similar nonequilibrium phase transitions as a consequence of phase slips due to random disorder in the frequency of Kuramoto oscillators (See Fig. S12 (right)).
c. Analytical modeling of noise-induced phase transitions: In order to understand the nature of the noise-induced phase transitions from one twisted state to the other, we model the adiabatic evolution of the phases as follows, here W p and W p are the winding numbers corresponding to different twisted states.At t = 0, the oscillators are in the W p -twisted state and evolve in an adiabatic manner to W p -twisted state at t = 1 as a consequence of the phase slip in the presence of noise.In the interval of the phase slip, one of the phase differences δθ 1,N undergoes a change from δθ 1,N to δθ 1,N + 2π, and it passes through the points δθ 1,N = π/2, 3π/2.At the non-equilibrium phase transition, the Jacobian corresponding to a given twisted state (W p ) starts with periodic boundary condition (PBC), then δθ 1,N = π/2, 3π/2 it realizes the open boundary condition (OBC) and finally, end up into the Jacobian with PBC corresponding to different twisted state (W p ).Interestingly, in the interval of the phase slip, the Jacobian is similar to the one-dimensional (1D) tight bonding Hamiltonian with generalized boundary conditions (See Sec.I A).During π/2 < δθ 1,N < 3π/2, the Jacobian has an eigenvalue with the positive real part which physically represents the appearance of the phase instability during the non-equilibrium phase transition from one stable twisted state to the other.From the perspective of the tight-binding model for the 1D chain, this eigenmode with positive eigenvalue is similar to the bound state localized at the edge of the system due to the generalized boundary conditions.The non-Bloch band theory (See Sec.I B) is also applicable here, however, we will not discuss this in detail (See for more detail in Secs.II C 3 and II C 4).

C. LCKM for N oscillators with non-reciprocal couplings
We consider the LCKM on a lattice with the periodic boundary condition consisting of N oscillators which are coupled through non-reciprocal couplings defined as K R = K L and denote ∆K = K R − K L .We first analyze the situation with ω i = 0, D = 0.In this case, the dynamical equations will modify to the following, looking at the above equations, it is straightforward to find out that the present problem is similar to the two oscillators problem.In this case, each oscillator moves with a finite common drift velocity having a constant phase difference with their neighbors.We also call these synchronized phases as chiral phases.These chiral phases are nothing but twisted states (See Sec.II B 1) with finite drift velocities which are characterized by distinct topological winding numbers.
Hence the total number of the chiral phases is proportional to the number of oscillators similar to the LCKM with reciprocal coupling.
It is important to note here that in contrast to the two oscillator problems, the stable chiral phases can appear for arbitrary values of non-reciprocal couplings.In this case of LCKM with a non-reciprocal coupling on a lattice with periodic boundary condition, the stable chiral phases appear as a consequence of chase and runaway motion and do not require random noise in contrast to globally coupled Kuramoto model with non-reciprocal couplings [8].We note here that for our case the periodic boundary condition is necessary to realize chiral phases which is also the case for realizing the twisted states of LCKM with reciprocal coupling (See Sec.II B 1). which reduces to for ω i = 0, and δ θi+1,i 1, ∀ i ∈ {1, 2, ..., N }.Here W p is the topological winding number.This is a remarkable result that indicates that synchronized solutions of this kind have quantized finite total velocity (in units of 2π∆K) equal to the topological winding numbers W p ∈ Z.
Chiral phases: Now we derive the expression of phases of the oscillators for the synchronized solutions which are characterized by the topological winding number.Considering ω i = 0, and f i = f 0 (const.)∀ i ∈ {1, 2, ..., N }, a. Jacobian of linear stability: We analyze the linear stability of the chiral phases via perturbing the system around them by a small amount such that θ i (t) → θi + φ i (t), with |φ i (t)| 1.Then the time evolution of φ i (t) is given by, In other words, we can write the above equations in vector form as dφ dt ≈ J φ, with J being a N × N matrix (also known as the Jacobian) and φ The above Jacobian matrix is similar to the HN model Hamiltonian for 1D lattice with non-reciprocal nearest neighbor hopping (See Secs.II C 3 and II C 4 for more discussion).b.Local stability analyses: We found through analytical and numerical analysis that the stability of the synchronized chiral phases depends not only on the phase differences between the neighboring oscillators (which was the case for the LCKM with reciprocal couplings) but also on the coupling strengths K R and K L .In order to analyze the stability of the chiral phase, we again utilize the Gershgorin theorem for the Jacobian Ĵ corresponding to Eq. (S51) with non-reciprocal couplings.In this case, the center a ii and radius R i of Gershgorin discs D(a ii , R i ) ⊆ C are as follows: The Jacobian for LCKM with non-reciprocal coupling always possesses a vanishing eigenvalue λ 1 = 0 (the Goldstone mode) due to the U(1) global rotation symmetry of Eq. (S51).Utilizing Gershgorin theorem together with numerical analyses for δ θi+1,i ∈ [−π, π], the results for the local stability analyses are summarized as follows: stability of the solutions may not be determined.
We have discussed that the chiral phases are the synchronized solutions with small and equal phase differences δ θi+1,i ≡ δ θ = 2πW p /N, ∀ i ∈ {1, 2, ..., N } which are characterized by the topological winding number.We show the behavior of the center of the Gershgorin disc a = −R i = −(K R +K L ) cos δ θ of the Jacobian matrix for chiral phases as a function of nearest neighbor phase difference δ θ in Fig. S14.It is easy to observe that for (K R + K L ) > 0, all centers of the Gershgorin discs a < 0 (a > 0) for chiral phases with |δ θ| < π/2 (> π/2).Accordingly, for (K R + K L ) > 0, all the non-zero eigenvalues of the Jacobian matrix for chiral phases with |δ θ| < π/2 (> π/2) have negative (positive) real parts (See the top and bottom panels (right) of Fig. S14).Therefore, for (K R + K L ) > 0, chiral phases with the eigenvalue spectra and hence the local stability also depend on the degree of non-reciprocity which is defined by ∆K = K R − K L .For example, in the case of fully non-reciprocal coupling (K R = −K L ), the eigenvalue spectra collapse on the Im[λ] axis consisting of degenerate eigenvalues which correspond to the marginal orbits.
Therefore, the stability analysis of the chiral phases is similar to that of the twisted states for LCKM with reciprocal coupling as both are described by the topological winding numbers.In nutshell, the LCKM with non-reciprocal coupling hosts stable time-dependent synchronized phases for finite sizes of the system.These stable solutions can be found explicitly depending on the initial conditions and their number is proportional to the finite size of the system.The synchronized solution with the topological winding number W p = 0 is globally stable.Whereas other stable synchronized solutions (chiral phases) with W p = 0 are only locally stable.We also observe the correlation between the winding number and stability of the chiral phases similar to the reciprocal case.

Stability of chiral phases: noise-induced non-equilibrium phase transitions
As discussed earlier that the different chiral phases are characterized by different winding numbers W p ∈ Z which are only locally stable having finite basins of attraction.Therefore, random and unbounded noise will induce phase transitions between different winding number sectors due to the phase slips (See Sec.II B 3 for the physical interpretation of phase slips) and eventually lead the system to the complete phase synchronized state.The main difference in comparison to the reciprocal case is that in the non-reciprocal case, the noise-induced phase transitions are characterized by the quantized jumps of the total common frequency v chiral = N i=1 dθi dt which is proportional to the topological eigenvalue form, Here J B is the boundary matrix.
Eigenstate solutions and generalized Brillouin zone: In order to have non-trivial solutions for c 1 and c 2 , Det[J B ] = 0 need to be satisfied together with z 1 z 2 = r 2 which takes the following form, The solutions of the above equation with z 1 z 2 = r 2 give rise to generalized Brillouin zones for the Jacobian of finite-size Kuramoto lattices with generalized boundary conditions.The paired generalized momenta (z 1 , z 2 ) can be expressed in the polar co-ordinates as follows In this section, we will prove that in the interval of the phase slip the Jacobian hosts exceptional points in its spectrum.
Results from numerical diagonalization: In the interval of the phase slip, δθ passes through the points π/2, 3π/2.At these points, the Jacobian realizes open boundary condition and we observe NHSE for the Jacobian.We also calculate the individual phase rigidity (defined in Eq. (S37)) for each of the eigenstates of the Jacobian matrix and find that the phase rigidity vanishes at δθ = π/2, 3π/2 in the large system size limit (See Fig. S19(iii) (left)).The vanishing phase rigidity can be understood as the consequence of the NHSE where the left eigenstates and the right eigenstates are localized at the opposite boundaries, causing the vanishing wavefunction overlap.At these points, the Goldstone mode (λ 1 = 0) coalesces with the normal mode which signifies the non-equilibrium phase transition and has been referred to as exceptional transition [8].
Results from non-Bloch band theory: In order to prove the existence of the exceptional transitions or exceptional points in the Jacobian spectra, we utilize the non-Bloch band theory and determine the GBZs in the full interval of the phase slip.In the interval of the phase slip, for t = 0, 1, the Jacobian realizes periodic boundary condition due to δθ = 2πW p /N, 2πW p /N and thus have two distinct GBZs (See Fig. S19 (i)(a)).Moreover, for δθ = π/2, 3π/2, we indeed observe the merging of two GBZs which defines the exceptional points in the spectrum of the Jacobian (See during the phase transitions between different chiral phases, the Jacobian realizes the exceptional transition of the generalized Brillouin zone as a consequence of the NHSE.
6. Noise-induced phase transitions: topological phase transition of GBZ for the Jacobian matrix: the emergence of positive eigenvalue bound state In the time interval of the phase slip for π/2 < δθ < 3π/2, we observe a positive eigenvalue of the Jacobian which shows the onset of the phase instability during the non-equilibrium phase transition.This eigenmode is described by the disjoint generalized momenta and hence localized at the boundaries of the system (See Fig. S19(i)(d)).Now we discuss the intrinsic topology of the GBZ for the Jacobian characterizing the positive eigenvalue eigenstate localized at the boundary of the system.It is clear from the discussion in the previous sections that the boundary value problem governing the Jacobian in the time-interval of the phase slip emulates the boundary value problem of the HN model in generalized boundary conditions with the following boundary equations, where

FIG. 2 .
FIG. 2. Top: (a), (b), and (c) dictating topological as well as exceptional phase transition of GBZs L1, L2 together with the complex vector plot of the meromorphic function related to non-Bloch topological invariant in the generalized boundary conditions -(I) (GBC-(I)) (See text).Middle left: non-Bloch topological invariant (W non−Bloch ) vs 1 which takes quantized value after 1 = tLr, and defines the topological phase transition of GBZs.Middle right: Non-Bloch topological invariant (W non−Bloch ) in ( 1, N )-plane.Here the value of W non−Bloch also represents the number of boundary states due to GBC-(I).Bottom left: Phase rigidity for eigenstates |ψ1 and |ψN vs onsite potential 1.The phase rigidity vanishes in large N -limit at 1 = tLr which defines the exceptional point and hence exceptional transition.Bottom right: Total phase rigidity in ( 1, N )-plane also captures the topological as well as exceptional phase transition of GBZs (See text).

FIG. 3 .
FIG.3.Bottom left: Eigenvalue spectrum of Jacobian for a chiral phase, bottom right: phase rigidity of the first eigenmode vs phase difference δθN during the interval of the phase slip (See text).Top panel: Generalized Brillouin zones for different δθN between [0, π] (the behavior of GBZs is similar in the other half).At δθN = π/2, the two GBZs merge with each other defining the exceptional transition.Parameters: KR = 2.5, KL = 1.9, N = 100.

FIG. 4 .
FIG. 4. Noise-induced phase transition: (a) time evolutions of phases in the presence of random white noise.The dotted lines mark the positions of phase slips which defines the phase transition in between different chiral phases.(b) Time evolution of the total chiral velocity (v chiral = v chiral /(KR − KL)) in the presence of random white noise.The discrete change of the total drift velocity serves as the physical observable of noise-induced non-equilibrium phase transition (See text).Parameters: KR = 2.5, KL = 1.9, D = 0.15, N = 100.

FIG
FIG. S1.Hatano-Nelson model for 1D chain in generalized boundary conditions.

) 2 .
Open boundary condition Now, we consider open boundary condition (OBC) with parameters t L = t R = 0, and 1 = N = 0 .In this case, Eq. (S13) can be simplified as sin[(N + 1)φ] = 0, (S21) which possesses N real physical solutions of φ m = mπ/(N + 1), where m ∈ {1, 2, ..., N }.Thus we have N pairs of generalized momenta (z 1 , z 2 ) = (re iφm , re −iφm ) with |z 1 | = |z 2 | = r, and as result the two GBZs coincide (See Fig. S2).The eigenvalue spectra are given by, and 1 = N = 0 .In this case, Eq. (S13) can be written as, sin[(N + 1)φ] − η 1 sin[(N − 1)φ] + η 2 sin[φ] = 0, (S25) which has N pairs of different physical solutions φ m ∈ C. In contrast to the IBC-(i), we find the localized boundary states in addition to the states of the bulk GBZs (See the left panel of Fig. S2 (d)).More specifically, among total N paired momenta, N − 2 pairs of generalized momenta with |z 1 | = |z 2 | form the continuum spectra in the complex plane.The two remaining pairs are isolatedly placed on the real momenta axis, i.e. z 1 , z 2 ∈ R, (for r < 1).As we explain later, these states correspond to the localized states characterized by the topological winding number.The eigenvalue spectra are again similar to that of PBC together with isolated real eigenvalues representing the localized boundary states due to GBC (See the right panel of Fig. S2 (d)).

1 .
Analytical examples (i): Generalized boundary conditions-(I) FIG. S4.Exceptional and topological phase transition of GBZ in the GBC-(I).Panels (i) and (ii) describe the evolution of GBZs and the emergence of the disjoint generalized momenta describing bound states as we tune the onsite potentials at systems boundaries.In the left panel of (iii), the emergence of boundary states is characterized by the quantized change of non-trivial values of the non-Bloch topological invariant (W non−Bloch ).The topological phase transition occurs when disjoint generalized momenta merge with GBZs (exceptional points) which is well captured by the change in W non−Bloch .The figure (right) of the panel-(iv) shows the vanishing of the phase rigidity (defined in Eq. (S37)) R1, R2, RN−1, RN at the topological phase transition points.The right panels (iii) and (iv) show the projection of the effective spinor on the GBZ and Bloch sphere in the topologically non-trivial phase (figure (ii)(e)) which clearly show the non-trivial winding and topology of GBZs (See Sec.I E 3 for more details).Parameters: N = 100, t L = 0, tL = 1.0, t R = 0, tR = 2.1, 1 = tR , N = tL , ∈ [−3, 3].

FIG. S5 .
FIG. S5.Exceptional and topological phase transition of GBZ in GBC-(II).Panels (i) and (ii) show GBZs and the disjoint generalized momenta for different values of the onsite potentials at the boundaries of the system.Left panels of (iii) and (iv) show the behavior of non-Bloch topological invariant and phase rigidity for the change in the onsite potentials(See Sec.I E 4 for more details).Parameters: N = 160, t L = tL = 2.1, t R = tR = 1.0, 1 = tR , N = tL , ∈ [−6, 6].The right panel of (iii) shows the number of bound states (n b ) in the ( , )-plane.In the right panel of (iv), the change in the non-Bloch topological invariant (W non−Bloch ) characterizes the topological phase transition points in the ( , )-plane, where topological boundary states emerge out of the bulk continuum (See Sec.I E 4 for more details).Parameters for right panels of (iii) and (iv): N = 80, t L = tL = 2.1, t R = tR = 1.0, 1 = tR , N = tL , with , ∈ [−6, 6].

)
Generalized Brillouin zones and boundary states: Similar to GBC-(II), solving the above boundary equations, we get N , N − 1, or N − 2 (depending on the strengths of t R and t L ) pairs of generalized momenta with |z 1 | = |z 2 | which form two GBZs.The two generalized Brillouin zones differ from each other except for the OBC (t R = t L = 0) FIG. S7.Exceptional and topological phase transition of GBZ in GBC-(III) and (IV).Panel (i) shows GBZs (L1, L2) together with disjoint generalized momenta for different values of boundary hopping terms.The left panels of (ii) and (iii) show the behavior of non-Bloch topological invariant and phase rigidity with respect to changes in the boundary hopping terms (t R = tR , t L = tL ) (See Sec.I E 5 for more details).Parameters: N = 160, tL = 1.0, tR = 2.1, 1 = 0, N = 0, ∈ [−6, 6].The right panels of (iii) and (iv) show the number of bound states (n b ) and non-Bloch topological invariant (W non−Bloch ) in the ( , )-plane with t R = tR , t L = tL .In the right panel of (iv), the change in the non-Bloch topological invariant characterizes the topological phase transition points in the ( , )-plane, where topological boundary states emerge out of the bulk continuum (See Sec.I E 5 for more details).Parameters for right panels of (iii) and (iv): N = 80, t L = tL , tL = 1.0, t R = tR , tR = 2.1, 1 = 0, N = 0, with , ∈ [−6, 6].
FIG. S8.GBC-(IV) with very large values of hopping terms at the boundaries.The top and bottom panels show the plots of eigenvalue spectra and GBZs for different values of boundary hopping terms t R = tR , t L = tL .For better visualization of eigenvalue spectra and GBZs, we have omitted the eigenvalue and generalized momenta corresponding to the bound states in (b), (c), (d), and (e).Parameters: N = 100, tL = e g , tR = e −g , 1 = N = 0, g = 0.1.
FIG. S10.Figures in the leftof the top and bottom panels show the variation of center of Gershgorin discs aii ≡ a = −2K cos δ θ with respect to constant nearest neighbor phase differences δ θ for twisted states with K > 0 and K < 0, respectively (See Sec.II B 2 for more details on the local stability of twisted states).Figures in the right of the top and bottom panels show the vector plots of stable twisted states (a < 0) with different winding numbers (Wp) for K > 0 and for K < 0, respectively.
FIG. S11.Left panel: behavior of maximal negative eigenvalue close to zero (Re[λ1 − λ2])of the Jacobian as a function of winding number (Wp) for different system sizes (N ).Right panel: Re[λ1 − λ2] vs N for a given winding number.Parameters: KR = 1.0,KL = 1.0.
FIG. S12.Stability of twisted phases: noise and disorder-induced phase transitions as a direct consequence of phase slips.Time evolution of the individual phases θi(t) (upper panel) and nearest neighbor phase differences δθi+1,i(t) (lower panel).The red dotted lines mark the noise (left panel) and disorder (right panel)-induced non-equilibrium phase transitions due to phase slips.It is clearly shown that in the time interval of the phase slips, one of the phase differences undergoes a change of 2π (See related discussion in the Sec.II B 3). Parameters: Left panel: KR = 1.0,KL = 1.0,D = 0.15, N = 100.Right panel: KR = 1.0,KL = 1.0, ∆ω = 0.25, N = 100..

FIG. S14 .
FIG. S14.Figures in the leftof the top and bottom panels show the variation of center of Gershgorin discs aii ≡ a = −(KR + KL) cos δ θ with respect to constant nearest neighbor phase differences δ θ for chiral phases with (KR + KL) > 0 and (KR + KL) < 0, respectively.Figures on the left of the top (bottom) panel show the eigenvalue spectra of the Jacobian matrix for stable and unstable chiral phases for (KR + KL) > 0 ((KR + KL) < 0) (See Sec.II C 2 for more details on the local stability of twisted states).
FIG. S16.Stability of chiral phases: noise-induced phase transitions as a direct consequence of phase slips and non-reciprocity.Left panel: time evolution of the individual phases θi(t) (upper panel) and total drift velocity v chiral = v chiral /(2π(KR − KL)), v chiral = N i=1 dθ i dt (lower panel).The red dotted lines mark the non-equilibrium phase transitions due to the phase slips (See Sec.II C 3 for more details).Right panel: time evolution of total drift velocity v chiral for different noise realizations (upper panel) and average total drift velocity v chiral (lower panel) (See Sec.II C 3 for more details).Parameters: Left panel: KR = 2.5, KL = 1.9, D = 0.15, N = 100.Right panel: KR = 1.9, KL = 1.0,D = 0.5, N = 100, number of noise realizations (N d = 200).
FIG. S17.Stability of the chiral phases: disorder-induced non-equilibrium phase transitions as a direct consequence of phase slip and non-reciprocity.Left panel: Time evolution of the individual phases θi(t) (upper panel) and total drift velocity v chiral = v chiral /(2π(KR − KL)) (lower panel).The red dotted lines mark the non-equilibrium phase transitions due to the phase slips (See Sec.II C 3 for more details).Right panel: Time evolution of nearest neighbor phase differences δθi+1,i(t) (upper panel) and total phase rigidity R = N i=1 Ri/N (lower panel) (See Sec.II C 3 for more details).Parameters: KR = 2.5, KL = 1.9, ∆ω = 0.5, N = 100.

TABLE I .
Definitions of different types of generalized boundary conditions.