Non-Hermitian Topological Metamaterials with Odd Elasticity

We establish non-Hermitian topological mechanics in one dimensional (1D) and two dimensional (2D) lattices consisting of mass points connected by meta-beams that lead to odd elasticity. Extended from the"non-Hermitian skin effect"in 1D systems, we demonstrate this effect in 2D lattices in which bulk elastic waves exponentially localize in both lattice directions. We clarify a proper definition of Berry phase in non-Hermitian systems, with which we characterize the lattice topology and show the emergence of topological modes on lattice boundaries. Interestingly, the eigenfrequencies of topological modes are complex due to the breaking of $\mathcal{PT}$-symmetry and the excitations could exponentially grow in time even in the damped regime. Besides the bulk modes, additional localized modes arise in the bulk band and they are easily affected by perturbations. We discuss realizations of these unique features in active metamaterials and potential applications in biological systems.

The study of topological band theory has been extended to open systems governed by non-Hermitian Hamiltonians , which are realized in classical systems subjected to intrinsic gain and/or loss such as optical , electric [86][87][88][89][90], and acoustic [91,92] structures. In contrast to the ordinary topological band theory, non-Hermitian systems exhibit unique features such as exceptional points, band structure which is sensitive to boundary conditions, and exponentially localized bulk modes ("non-Hermitian skin effect"). Among them the most fascinating subject is the interplay between non-Bloch bulk waves and the characterization of lattice topology. It is thus intriguing to ask: can exotic mechanical properties arise in energy non-conserving systems whose Hamiltonian is non-Hermitian?
In this letter, we study non-Hermitian topological lattices composed of mass points and meta-beams that lead to odd elasticity. Odd elasticity originates from active matters with microscopic interactions which do not conserve energy [93], and is ubiquitous in a broad range of natural (such as fiber networks with driving [94][95][96] and active microtubule networks [97]) and manmade materials (such as coupled gyroscopes [7,18,98,99] and active fluids [100,101]).
Up to date, major efforts of non-Hermitian mechanics (NHM) have been limited to 1D parity-time (PT ) symmetric systems [102,103] whose eigenvalues are real and the eigenmode amplitudes do not grow in time. Our research lifts this PT -symmetry by allowing the eigenvalues to be complex, meaning that in the dissipationless limit the lattice is unstable against infinitesimal stimulations. We show that a damping which counteracts the intrinsic energy production naturally stabilizes the lattice. We first realize the basic two-band models in 1D and 2D rotor lattices, and then in the third honeycomb lattice model we study a four-band non-Hermitian Hamiltonian. Interestingly, this four-band honeycomb lattice exhibits two features not observed in the Hermitian counterpart. First, the eigenvalues of topological modes are complex due to the breaking of PT -symmetry, and the excitation exponentially increases if the damping is not strong enough to counteract internal energy production. Second, besides the bulk modes, on lattice boundaries we observe additional localized modes whose eigenvalues are not separated from the bulk band. These non-topological localized modes may be essential for the in-band boundary softness in energy non-conserving systems such as biological structures [104][105][106][107][108]. We leave NHM of amorphous systems in future research.
The models-We model the non-conservative oddelastic interaction between particles as a pairwise force F ( u), where u is the relative displacement away from equilibrium between an interacting pair of particles. According to Ref. [93], to linear order, we use the force wheren is the unit vector along the connection orientation,φ is the unit vector rotated fromn by 90 • counterclockwise [ fig.1(a)], k is the spring constant, and k o represents the strength of the energy non-conserving force (dubbed odd elastic constant). The unit cells of all three lattices are subjected to fixed boundary conditions (FBCs), and are composed of two mass particles labeled A and B with mass m which are connected by oddelastic meta-beams. The cell is labeled by (n 1 , . . . , n d ), and the site displacements are denoted as u A,(n1,...,n d ) and u B,(n1,...,n d ) , where d is the spatial dimension, n i (1 ≤ n i ≤ N i ) is the cell labeling, and N i is the lattice length scale. We consider the Newtonian equation of motion and take the ansatz in every model, where β i is the decay rate of the non-Bloch bulk modes [27]. In all three models, |β i | stays a constant while Arg β i varies from 0 to 2π. We denote β = (β 1 , β 2 , . . . , β d ) for simplicity. The equation of motion in every lattice is simplified as an eigenvalue problem where D( β) is the dynamical matrix, u is the unit cell displacement field, and λ is the eigenvalue. Generalized Berry phase in non-Hermitian systems.-We first generalize the Berry phase for the N -band non-Hermitian dynamical matrix D( β) to characterize the lattice topology. We consider the simple case that all bands are separable bands [30] (i.e., eigenvalues λ n ( β) = λ n ( β ) for all β, β as long as 1 ≤ n = n ≤ N ). Thus the system is free of exceptional point. Based on the biorthogonality and completeness of the eigenbasis, we define the generalized Berry phase of band n as follows, where C i is the closed loop trajectory connecting Arg β i = 0 and 2π, A (n) i = i n L |∂ Arg βi |n R is the generalized Berry connection, and |n R ( β) and |n L ( β) are the right and left eigenvectors, respectively. We note that Eq.(4) is the proper definition for one of the Berry phases of Ref. [30]. For non-Hermitian systems higher than 1dimension, we can further define the generalized Berry curvature as Ω i . By summing over all energy bands, we can prove that In the later discussions it is useful to study the following form of non-Hermitian dynamical matrix: where h 1 , h 2 are 1 2 N × 1 2 N invertible matrices (N is even) which are free of exceptional point, σ ± = (σ x ± iσ y )/2, and σ x,y,z are the Pauli matrices. Since the eigenvalues appear in ±λ n pairs, we study the generalized Berry phase for all 1 2 N bands with −λ n eigenvalues. The generalized Berry phase is quantized by the symmetry operator Π = σ z ⊗ I 1 In SI(III) we show that the Berry phase can be expressed in the simple form, 1D and 2D rotor lattices.-We first realize the basic non-Hermitian topological mechanics in 1D rotor lattice. The unit cell is composed of mass particles that rotate freely about fixed pivots [3]. The nearest neighbor rotor particles are connected by meta-beams that possess odd elasticity [ fig.1(d)]. The chain consists of N 1 unit cells with the particles labeled A n1 and B n1 (1 ≤ n 1 ≤ N 1 ) and is subjected to FBCs [109,110] at B 0 and A N1+1 . The equilibrium configuration is that the A-rotors (B-rotors) make an angle θ A (θ B ) relative to the upward (downward) normals, and the meta-beams connecting A n1 and B n1 (B n1 and A n1+1 ) make an angle −θ 1 (θ 2 ) relative to the horizontal line. We study the eigenvalue problem Eq. (3), where d = 1 is the spatial dimension, h z = (a − a )/2, ω is the eigenfrequency, λ = mω 2 − (a + a )/2, u = (u A /|b| , a, b, c 1 and a , b , c 1 are constant parameters determined by k, k o , θ A,B and θ 1,2 (see SI(V)), and sgn(. . .) is the sign function. We require bb c i c i | i=1,...,d > 0 so that the bulk modes exist, and these bulk modes are exponentially localized near the lattice boundaries ("non-Hermitian skin effect") with the decay rate In general, the generalized Berry phase of the considered rotor lattice is not quantized. However, it can be quantized by imposing the condition h z = 0, which is practically achieved [109] by enabling the nearest neighbor springs perpendicular to each other. The topological phase is characterized by the generalized Berry phase, which is computed by Eq. (6). The topological phase transition occurs when |bb | = |c 1 c 1 |. The chain is topologically trivial when γ 1 = 0 and |bb | > |c 1 c 1 |, while the doubly degenerate topological modes [26,27] emerge on lattice boundaries when γ 1 = π and |bb | < |c 1 c 1 |. They both localize on the left boundary if |b| > |c 1 | and |b | < |c 1 | (right boundary if |b| < |c 1 | and |b | > |c 1 |), while they localize on different boundaries if |b| < |c 1 | and |b | < |c 1 |.
The unit cell of the 2D rotor lattice is shown in fig.2(a), where each rotor particle is connected to three odd-elastic meta-beams labeled i = 1, 2 and 3. The spring constants and the odd-elastic constants are denoted as k i and k o i , respectively. The lattice consists of N 1 × N 2 unit cells whose sites are labeled by A n1,n2 and B n1,n2 , and is subjected to FBCs. The equilibrium configuration is that the tangential directions of A-rotors and B-rotors make angles θ A and θ B to the x-axis, respectively, and the meta-beam labeled i makes an angle θ i to the x-axis. The dynamical matrix D is given by Eq. (7), where d = 2 is the spatial dimension (details see SI(V)). The bulk modes are exponentially localized on lattice boundaries with the decay rates β 1 and β 2 given by Eq. (8), which manifest non-Hermitian skin effect in 2D systems.
The generalized Berry phase of 2D rotor lattice is not quantized unless h z = 0, which is practically achieved by letting k i / sin(2θ i +1 − 2θ i +2 ) = Const. for i = 1, 2, 3. When the bands are separable, the topological phase is characterized by the generalized Berry phase. The lattice is free of topological edge modes when (γ 1 , γ 2 ) = (0, 0), while the topological modes emerge on the upper-left and bottom-right lattice boundaries when (γ 1 , γ 2 ) = (π, 0), or on the upper-right and bottom-left boundaries when (γ 1 , γ 2 ) = (0, π). In the non-separable region the two bands touch at points β (w) and β * (w) . It is notable that in the Hermitian case when b = b , c 1 = c 1 and c 2 = c 2 , these two band-touching points correspond to mechanical Weyl points [10]. The extension of Weyl points to non-Hermitian systems will be our future research.
Both 1D and 2D rotor lattices are stable in the branch bb > 0 and c i c i | i=1,...,d > 0. The topological modes can be excited by employing an external monochromatic shaking force. However, the lattices are unstable in the branch bb < 0 and c i c i | i=1,...,d < 0 in the dissipationless limit. Therefore, we ask all eigenmodes to decay in time by counteracting the intrinsic energy production with a damping f = −ηu, where However, due to the real eigenvalues of topological modes, the excitations are heavily damped, which is in sharp contrast to the complex eigenvalues of the topological modes in the following honeycomb lattice. Exceptional points emerge in both rotor lattices when bb < 0, c i c i | i=1,...,d < 0, and h z = 0. The eigenvalues coalesce at zero at the exceptional points β (e) and β * (e) , and only one eigenvector |n R ( β (e) ) (|n R ( β * (e) ) ) instead of two occurs at each point. In general this may not fulfill the reality condition of particle displacements. However, in SI(IV) we prove that, as long as the coalesce eigenvalue is real, the two coalesce eigenvalues λ( β (e) ) and λ( β * (e) ) must degenerate. The linear combination of |n R ( β (e) ) and |n R ( β * (e) ) is real for all particle displacements.  (1,15) with ωext = 1.13 and F = 10 −10 . We plot A/F versus t red (blue) curves for A (1,15) (B (1,15) ), where A is the x-direction oscillation amplitude and t varies from 0 to 2000 × (2π/ωext). Solid curves for η = ηc − δη and dashed curves for η = ηc + δη, where ηc = 0.0107 and δη = 0.001 ≈ ηc/10. Displacement-time relations are plotted from t0 = 1.10 × 10 4 to t0 + 2 × (2π/ωext).
2D honeycomb lattice.-In rotor lattices the generalized Berry phases are not quantized unless we fine-tune their parameters. Therefore, a lattice whose generalized Berry phase is quantized without fine-tuning is in need, and it is the honeycomb lattice as follows. Each site is connected to three odd-elastic meta-beams [ fig.3(a)]. The beam orientation, the spring constant, and the odd elastic constant are denoted as θ i , k i and k o i , i = 1, 2, 3, respectively. The lattice is subjected to FBCs. We study the lattice from Eq.(3), where λ = mω 2 , D is a four-band dynamical matrix and h β1,β2 is a 2 × 2 matrix specified in SI(V). β 1 and β 2 of bulk modes happen to lie on a unit circle, |β 1 | = |β 2 | = 1. Therefore, all bulk modes are extended in space instead of being localized on lattice boundaries. Hence, the bulk mode eigenvalues under FBCs are the same as those under PBCs.
To observe in-gap topological edge modes, we consider the dynamical matrix whose eigenvalues Re λ 1,2 (β 1 , β 2 ) < Re λ 3,4 (β 1 , β 2 ) for all (β 1 , β 2 ) and (β 1 , β 2 ), meaning that the first and second two bands are separated by the band gap ∆ [ fig.3(b)]. Different from the topological rotor lattices free of exceptional points, two exceptional rings arise in the (Arg β 1 , Arg β 2 ) parameter space [SI(VII)] when the eigenvalues of the second two bands coalesce. Interestingly, the first two bands have no exceptional point, providing a well-defined Berry phase to characterize the lattice topology. The dynamical matrix yields the symmetry property ID( β) is not quantized at each Arg β 2 , when averaged over Arg β 2 the Berry phase γ 1 is quantized [111,112]: γ 2 is also quantized to 0 or π when averaged over Arg β 1 . This model presents several fascinating features not observed in the Hermitian counterpart. First, the eigenvalues of topological modes are complex due to the breaking of PT -symmetry, meaning that in the dissipationless limit, the eigenmode excitations exponentially grow in time. An external damping can neutralize the internal energy production. However due to the large damping, both topological modes and bulk modes are excited when the lattice is shaken by an external force with frequency Re (ω topo ), where ω topo is the topological mode eigenfrequency. Hence, we ask the excitations of bulk modes to be small by imposing a lower bound of the bandgap ∆ η Re (ω topo ). The unique consequence of complex eigenvalues is that the excitations of topological modes exponentially increase in time if η η c (η c given by Eq. (9)), and they saturate if η η c [ fig.3(e)]. The second peculiarity is that besides the bulk modes which are extended in space, we observe additional modes whose eigenvalues are not separated from the bulk band, and they are localized on the lattice boundaries. Compared to topological modes, these localized modes are not topologically protected and are easily affected by perturbations [SI(VII)]. It is therefore interesting to seek exotic edge mode responses driven by external force with bulk mode frequencies.
Discussions.-In this letter we extend "non-Hermitian topological theory" to mechanical systems. Our study is based on but not limited to odd elasticity. Any active matter that injects energy could be the element of non-Hermitian metamaterials. Now we discuss potential applications in biological structures.
First, the dynamics of biological structures are challenged by the ubiquitous heavy dissipation [104][105][106]. It is worth applying NHM to active biological systems in which the heavy damping is counteracted by energy production, such as motor activities in fiber networks and motility of cell sheets. So far the study of non-Hermitian systems is limited to periodic lattice structures. Thus, the extension of NHM to amorphous structures is useful to explore energetic biological systems.
Second, in addition to topological modes, there also exist non-topological localized modes in non-Hermitian systems. It is thus intriguing to ask if these modes are associated with the in-band boundary softness or exotic mechanical responses of active biological systems [107,108].
Acknowledgement.-D. Z. would like to thank insightful discussions with Vincenzo Vitelli, Po-Yao Chang, Anton Souslov and Colin Scheibner. D. Z. is supported in part by the National Science Foundation (Grant No. NSF-EFRI-1741618).
Note added.-After this work was completed, we became aware of an independent work by C. Scheibner, W. Irvine, and V. Vitelli [117]. We thank them for bringing their research to our attention.

I. BIORTHOGONALITY AND COMPLETENESS OF GENERALIZED EIGENVECTORS IN NON-HERMITIAN MATRIX
We start by considering a non-Hermitian N ×N matrix D with eigenvalues denoted as λ n , 1 ≤ n ≤ N . The corresponding generalized right eigenvector and left eigenvector of rank l are denoted as |n R (l) and n L (l) |, respectively. They satisfy The generalized eigenvectors is a complete basis. In particular, ordinary eigenvectors (denoted as |n R and n L |) are rank 1 generalized eigenvectors. It is straightforward to prove that acting D − λ n I on the generalized right (left) eigenvector |n R (l) ( n L (l) |) gives the linear superposition of lower rank generalized right (left) eigenvectors: where R n,(k,l) (L n,(k,l) ) is the coefficient of the rank k generalized right (left) eigenvector of the eigenvalue λ n . It is obvious that R n,(l−1,l) = 0 (L n,(l−1,l) = 0). We denote the highest rank of left (right) generalized eigenvector of the eigenvalue λ n as l n (r n ). It is easy to prove that l n = r n . To this end, without loss of generality we assume r n > l n (i.e., r n ≥ l n + 1). We calculate 0 = n L (ln−l) |(H − λ n I) rn−1 |n R (rn) for all l = 0, 1, . . . , l n − 1 to find |n R = 0 which is not true. Therefore, in what follows we use l n to represent the highest rank of the left and right generalized eigenvectors. From now on we define R n,(l ,l) | l ≥l = L n,(l ,l) | l ≥l = 0 and then Eq.(S3) can be re-written as This can be proved by substituting Eq.(S3) into n L (l ) |(D − λ n I)|n R (l) repeatly in every step of the second principle of mathematical induction. Next, for the left and right generalized eigenvectors of the same eigenvalue λ n , the biorthogonality is presented as follows, This is proved by using the mathamatical induction and computing n L (l ) |(D − λ n I)|n R (l) in every step. Rescaling the right eigenvector |n R (l) → |n R (l) / n L (l ) |n R (l) , and summarizing Eqs. (S5, S6), we obtain the normalized biorthogonality of the generalized eigenvectors, n L (l ) |n R (l) = n R (l) |n L (l ) = δ n,n δ l+l ,ln+1 .
Based on Eq.(S7) it is easy to prove the completeness of the generalized eigenvectors, In particular, if all highest ranks of generalized eigenvectors are l n = 1, the generalized eigenbasis is the same as the ordinary eigenbasis, and there is no exceptional point in the non-Hermitian matrix D. Eqs.(S7, S8) are reduced to and Finally, based on Eq.(S7) we further show the coefficients R n,(k,l) and L n,(k,l) in Eq.(S4) are related. To this end, we calculate n L (l) |H − λ n I|n R (r) to find which is valid for all l n ≥ l, l ≥ 1.

II. GENERALIZED BERRY PHASE, BERRY CONNECTION AND BERRY CURVATURE
Newton's equation of motion is the second order derivative in time, while Schrodinger's equation is the first order derivative. In order to define the generalized Berry phase in Newtonian equation of motion, we consider the simple case that for the N × N non-Hermitian dynamical matrix D, all bands are separated from each other [30]. Hence the D-matrix is free of exceptional points, and the left and right eigenbasis are complete. We use the eigenvalue equation, where |n R is the right eigenvector corresponding to eigenvalue λ n . We then define an auxilliary wave function which evolves as the auxilliary parameter t advances.
It also satisfies the eigenvalue equation D|u R n (t ) = λ n |u R n (t ) . Since any N -dimensional wave function |ψ R (t ) can be expressed as the linear superposition of the complete basis {|u R n (t ) }, we find that |ψ R (t ) is subjected to the Schrodinger-like equation of motion Starting from Eq.(S14), we derive the generalized Berry phase of non-Hermitian mechanical systems. Based on Eq.(2) of main text, we use the ansartz ( u A,(n1,...,n d ) , u B,(n1,...,n d ) ) = β n1 1 . . . β n d d ( u A , u B ) to reduce the lattice equation of motion as D( β)u = λu, where the dynamical matrix D( β) is the function of β. We adiabatically evolve β(t ) as the auxilliary parameter t advances. In SI(V) we prove that, for all three models we study in this paper, the modulus |β i | keeps constant while the complex phase Arg β i varies from 0 to 2π. The starting eigenvector is denoted as |n R ( β(t = 0)) . According to the quantum adiabatic theorem [113][114][115], a system initially in one of its eigenstate |n R ( β(0)) will stay as an instantaneous eigenstate of the Hamiltonian D( β(t )) throughout the process. The only degree of freedom is the phase of the state. Write the state at t as and insert it into Eq.(S14), we obtain (dγ . (S16) Next, according to the normalized biorthogonality and completeness of the eigenbasis, we multiply the above equation on the left hand side with the left eigenvector n L ( β(t ))|, to obtain where C is the closed path such that the variable β(0) = β(T ) goes back to itself, and is the generalized Berry connection. Next, we define the generalized Berry curvature In what follows, we prove that in the non-Hermitian system which is free of exceptional points, the generalized Berry curvature satisfies To this end, we exploit the following relation to simplify the generalized Berry curvature as follows, By summing over all N bands we prove Eq.(S20).

III. GENERALIZED BERRY PHASE QUANTIZED BY SYMMETRIES
In this section we show how generalized Berry phase of the non-Hermitian dynamical matrix D( β) is quantized by three different symmetries [116].
In the first case, we consider a 1D lattice with the dynamical matrix D(β 1 ) subjected to the symmetry prop- where M x is the symmetry operator which satisfy M x = M −1 x , and β * 1 is the complex conjugate of β 1 . The eigenvalues of M x are constrained to ±1. The generalized Berry phase is quantized by this symmetry. To prove this, we notice which means where R (n) x (β 1 ) is the symmetry phase connecting M x |n R (β 1 ) and |n R (β * 1 ) . At high symmetry points when Arg β 1,hs = 0, π ("hs" is short for "high symmetry"), we find [M x , D(β 1,hs )] = 0, meaning that M x and D share the same eigenvectors |n R (β 1,hs ) , and it satisfies M x |n R (β 1,hs ) = ±|n R (β 1,hs ) . (S25) Thus at high symmetry points, the symmetry phase R (n) x (β 1,hs ) = 0 or π mod 2π. Similarly, the left eigenvector obeys the relation where L (n) x (β 1 ) is the symmetry phase factor connecting n L (β 1 )|M x and n L (β * 1 )|. It satisfies L (n) x (β 1 ) = R (n) x (β 1 ) mod 2π due to the normalized biorthogonality n L |n R = 1. Thus, the generalized Berry phase of band n is γ (n) 1 = R (n) x Arg β1=π Arg β1=0 = 0 or π mod 2π, (S27) It is notable that in the Hermitian case, β 1 = e iq1 and M x is the reflection symmetry operator.
In the second case we consider a 2-dimensional lattice with the Hamiltonian D( β = (β 1 , β 2 )) subjected to the symmetry property ID( β)I −1 = D( β * ), where the symmetry operator I = I −1 . The generalized Berry phase is quantized by this symmetry. To this end, we notice which means where R (n) I ( β) is the symmetry phase connecting I|n R ( β) and |n R ( β * ) . At high symmetry points β hs when (Arg β 1,hs , Arg β 2,hs ) = (0, 0), (0, π), (π, 0), (π, π), we find [I, D( β hs )] = 0, meaning that I and D share the same eigenvector |n R ( β hs ) which satisfies I|n R ( β hs ) = ±|n R ( β hs ) . (S30) Thus, at high symmetry points the symmetry phase R (n) I ( β hs ) = 0 or π mod 2π. The left eigenvectors follows the same transformation rule under the symmetry operator I, where it is straightforward to prove L   Since all bands are separable, we find γ  can be expressed in terms of symmetry operator eigenvalues: where Y = (0, π) and M = (π, π) stand for the high symmetry points of the parameter space (Arg β 1 , Arg β 2 ). Although γ (n) 1 (Arg β 2 ) is not quantized, it is quantized when averaged over the complex phase Arg β 2 , which is equivalent to For the bands which are separated from other bands, the summation of generalized Berry phase is also quantized, The above results holds true for γ 2 . It is notable that in the Hermitian case, β = (e iq1 , e iq2 ) and I is the inversion symmetry operator.
In the third case, we consider the Hamiltonian D( β) subjected to the symmetry property ΠDΠ −1 = −D where the symmetry operator Π = Π −1 . It follows that which means the eigenvalues of D come in ±λ n pairs (we let λ n > 0). Hence N/2 bands have +λ n eigenvalues and N/2 bands have −λ n eigenvalues, where N is an even number. In addition, by comparing to the eigenvalue equation we find the states |n R and | − n R are related by a phase factor e iR (n) Π , The same relation is true for the left eigenvectors, The generalized Berry phase γ 1 (Arg β 2 ) at each Arg β 2 is quantized under this symmetry, where the summation is to sum over all N/2 bands with −λ n eigenvalues. The same is also true for γ 2 (Arg β 1 ). We note that Π is the chiral symmetry in the Hermitian case.
A particularly interesting non-Hermitian dynamical matrix subjected to this Π symmetry is of the following form, where N is even, h 1 and h 2 are 1 2 N × 1 2 N matrices. The symmetry operator is Π = σ z ⊗ I 1 2 N × 1 2 N . In the Hermitian case we require h 1 = h † 2 , and this is lifted in non-Hermitian systems. We now attempt to express the generalized Berry phase in terms of h 1 and h 2 . To this end, we denote |n R = (|n R (S47) Together with the normalized biorthogonality n L |n R = 1, we have As an intermediate step, the Berry phase of the 1 2 N bands with −λ n eigenvalues is simplified as The first term of Eq.(S49) is zero modulo 2π since it's the summation of the Berry phase of all bands of the Hamiltonian h 1 h 2 . Exploiting the identity Finally, the Berry phase of all −λ n bands is given by
If the coalesce eigenvalue λ n ( β (e) ) is real, the eigenvalues of exceptional points β (e) and β * (e) must degenerate. Hence, the linear combination of the eigenmodes |n R ( β (e) ) +|n R ( β * (e) ) indicates that all particle displacements are real. Here we show calculation details of non-Hermitian skin effect in 1D rotor chain. Let us first denote The parameters a, b, c and a , b , c which appear in the D matrix of the 1D chain are displayed as follows, The displacement field of the chain is denoted as {u An , u Bn }, which streches the meta-beams connecting the rotor particles. The Newton's equation of motion is mü An = (bu Bn + c 1 u B,n−1 ) − au An subjected to fixed boundary conditions We use the ansartz (u An , u Bn ) = β n 1 (u A , u B ) to simplify the equation of motion as Du = λu, with u = (u A /|b| where h z = (a − a )/2 and λ = mω 2 − (a + a )/2. Given the eigenvalue λ, we find that det(D − λI) = 0 is the second order equation of β 1 . Solving this equation gives us two solutions β where A = bc 1 , B = bb + c 1 c 1 − λ 2 + h 2 z , and C = b c 1 . The eigenvector correspoding to β (i) 1 is denoted as u (i) . The general wave function is given by which is subjected to the fixed boundary conditions. Eliminating the eigenvectors gives us 1 ). (S60) We are concerned with the bulk modes for a long chain in the limit N → ∞, which requires |β (1) 1 | = |β (2) 1 | for the bulk modes [27]. If the bulk mode condition is not fullfilled, without loss of generality we let |β = 0, which gives a single β 1 solution instead of two, and this solution is independent of the chain length N . In order to have the bulk mode condition |β In general the decay rate |β 1 | of the bulk modes is not 1, a unique feature of non-Hermitian systems dubbed the "skin effect".

B. Non-Hermitian skin effect in 2D rotor lattice
To study the non-Hermitian skin effect in 2D honeycomb rotor lattice, we first denote a set of parameters which will be used later: subjected to fixed boundary conditions, for all n 1 = 1, 2, ..., N 1 and n 2 = 1, 2, ..., N 2 . By applying the ansartz (u An1,n2 , u Bn1,n2 ) = β n1 1 β n2 2 (u A , u B ), the Newton's equation of motion is simpified as Du = λu, where and u = (u A /|b| . Similar as the 1D rotor chain, given the eigenvalue λ and decay rate β 2 (β 1 ), we find det(D − λI) = 0 is the second order equation of β 1 (β 2 ). Solving this equation gives us two solutions β where A = c 1 (b + c 2 /β 2 ), B = bb + c 1 c 1 + c 2 c 2 + c 2 bβ 2 + c 2 b /β 2 − λ 2 + h 2 z , and C = c 1 (b + c 2 β 2 ). By denoting the eigenvector as u (i) which corresponds to β (i) 1 , we then express the general wave function as the linear superposition of the eigenvectors u (1) and u (2) : subjected to fixed boundary conditions. Eliminating u (i) gives us the following relationship, The bulk mode condition [27] requires |β Solving Eq.(S71) gives us the decay rate of β 2 : By substituting |β 2 | into Eq.(S68) we further obtain the decay rate β 1 , In general, these bulk mode decay rates are not 1, which are the manifestation of non-Hermitian skin effect in 2D systems.
(2) If λ is real, we can prove that a 1,2,3 , a 1,2,3 and a 4 are real, and they satisfy a 1 = a 1 , a 2 = a 2 , and a 3 = a 3 . (3) If det(D(β 1 , β 2 )−λI) = 0, then det(D(β −1 1 , β −1 2 )−λI) = 0 as well. Equipped with the above three properties, we are now ready to study the non-Hermitian skin effect. Given λ and β 2 , solving the equation det(D(β 1 , β 2 ) − λI) = 0 gives us two solutions β where A = a 1 + a 3 /β 2 , B = a 2 β 2 + a 2 /β 2 + a 4 and C = a 1 + a 3 β 2 . The eigenvector which corresponds to β (i) 1 is denoted as u (i) . Thus, the wave function of the general form is Employing fixed boundary conditions and eliminating u (i) gives us the simplified relation between β 1 and β Similar as the rotor lattices, the bulk mode condition requires |β 1 | = 1. In summary, we find Eq.(S83) is obtained by solving det(D(β 1 , β 2 ) − λI) = 0 in terms of β 1 , but solving this determinant in terms of β 2 offers us the same result. It means that instead of localizing on the lattice boundaries, all bulk modes are extended. The generalized Brillouin zone happens to be a unit circle. An important consequence of Eq.(S83) is that β * i = β −1 i , which is essential for the symmetry property ID( β)I −1 = D( β * ) which quantizes the generalized Berry phase. Despite the fact that all bulk modes are extended, the dynamical matrix of the honeycomb lattice is still non-Hermitian.

VI. HONEYCOMB LATTICE CONFIGURATION WITH COMPLEX EIGENFREQUENCIES OF TOPOLOGICAL MODES
In order to find a honeycomb lattice configuration which possesses complex eigenfrequencies of topological boundary modes, this four-band non-Hermitian dynamical matrix has to fullfill five criteria. (1) The first and the second two bands are separated by the band gap ∆.
(2) The generalized Berry phase of the first two bands (γ 1 , γ 2 ) = (π, 0) or (γ 1 , γ 2 ) = (0, π). (3) The eigenfrequencies of topological edge modes are complex. (4) The critical damping coefficient η c , the topological mode eigenfrequencies ω topo , and bulk mode eigenfrequencies ω bulk have to satisfy the following relationship Eq.(S84) implies that the excitation of topological boundary modes are marginal. The amplitude exponentially grows in time when η η c , while it yields a constant if η η c . The excitations of bulk modes are bounded in both cases. (5) Driven an external shaking with the frequency ω ext = Re (ω topo ), both of the topological modes and the bulk modes are excited due to the damping. We impose a lower bound to the band gap ∆ η Re (ω topo ) (here we let ∆ 5η Re (ω topo ) in the main text), and then the excitations of bulk modes are weak compared to the topological modes.
Most of the topological lattice configurations which fullfill criteria (1), (2) and (3) do not accomplish the rest two criteria (4) and (5). To seek the parameter region which satisfies all five criteria presented above, we start from a special configuration in which all eigenvalues of this non-Hermitian dynamical matrix are real and positive. Small parameter deviations from this configuration render eigenfrequencies with small imaginary parts, which means a small damping η is strong enough to counteract the energy injection, and then criterion (5) is assured. Finally, we go through the neighborhood of this special configuration to find the parameter region which satisfies criterion (4). One of such a set of parameters which satisfy all five criteria is depicted in fig.3(b) of the main text.
We now discuss how to find this special configuration, in which all eigenvalues λ = mω 2 are real and positive. This configuration can be achieved by allowing h(1, 1) of Eqs.(S74) to be an identity, which in turn demands Given the normal spring constants k 1 , k 2 and k 3 and bond orientations θ 1 , θ 2 and θ 3 , one can determine the odd elastic constants k o 1 , k o 2 and k o 3 accordingly. Through analytic calculations, we prove that all bulk mode eigenvalues are real and positive. We then numerically solve the topological boundary modes and find that their eigenvalues are real and positive, too. In summary, all eigenvalues are real and positive numbers in this special configuration.
Deviations of the parameters k 1 , k 2 , k 3 , θ 1 , θ 2 , θ 3 , k o 1 , k o 2 , and k o 3 lead to a wide range of lattice configurations with complex topological mode eigenfrequencies, which is depicted in fig.3(a).

VII. ADDITIONAL RESULTS
A. Finite region of parameters that leads to topological 2D rotor lattices Here we show that the topological configuration of 2D rotor lattice presented in fig.2(a) of the maintext is not unique. There is a finite region of parameters which leads to topological modes localized on lattice boundaries [ fig.S1]. To this end, we illustrate a 2-dimensional phase diagram by varying the B-site tangential direction θ B and the odd elastic constant (divided by k 3 ) k o /k 3 , while we keep the primitive vectors a 1 = (− √ 3 2 , 3 2 ) and a 2 = ( √ 3 2 , 3 2 ), the rotation radii of the sites l A = l B = 0.1, and the A-site tangential direction θ A = 172 • the same as the unit cell configuration shown in fig.2(a) of the maintext. The normal spring constants are given by k i / sin(2θ i +1 − 2θ i +2 ) = Const for i = 1, 2, 3 accordingly.

B. Finite region of parameters that leads to topological honeycomb lattices
Similar as the 2D rotor lattice, there is a finite parameter region which leads to topological honeycomb lattices.
In fig.S2(a) we show a phase diagram by displacing the A-site from (0, 0) to r a = (x a , y a ) while we keep the other parameters of the maintext fig.3(a), including site B position r b = (0, 1), primitive vectors a 1 = (−  fig.S3(a)]. Besides these in-gap topological boundary modes, we observe additional modes whose eigenvalues are not separated from the bulk bands and they also localize on lattice boundaries. These modes are not topologically protected and can be easily affected by perturbations. To prove this, under FBCs we plot the Im λ v.s. Re λ figure of the lattice whose particle mass of the n 1 = 1 layer are replaced by m A = m B = 0.7m [ fig.S3(c)]. We find that, while the eigenvalues of topological modes are insensitive to perturbations, the eigenvalues of these additional localized modes are easily changed by perturbations.