††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# Stress control in non-ideal topological Maxwell lattices via geometry

Harry Liu  Department of Physics, University of Michigan, Ann Arbor, MI 48109,
USA  Siddhartha Sarkar  Department of Physics, University of Michigan, Ann
Arbor, MI 48109, USA  A. Nafis Arafat  Department of Physics, University of
Michigan, Ann Arbor, MI 48109, USA  City University of New York, New York, NY
10017, USA  Ethan Stanifer  Department of Physics, University of Michigan, Ann
Arbor, MI 48109, USA  Stefano Gonella  Department of Civil, Environmental, and
Geo- Engineering, University of Minnesota, Minneapolis, MN 55455, USA
Xiaoming Mao<EMAIL_ADDRESS>Department of Physics, University of Michigan,
Ann Arbor, MI 48109, USA

###### Abstract

Topological mechanical metamaterials have demonstrated exotic and robust
mechanical properties which led to promising engineering applications. One of
such properties is the focusing of stress at the interface connecting domains
of topological Maxwell lattices of opposite topological polarizations, which
protects the bulk of the material against fracturing. Here we generalize this
theory to non-ideal Maxwell lattices, incorporate real material features that
leads to interactions beyond previous ideal models. By quantitative analysis
of stress distributions of topological Maxwell lattices with self-stress
interfaces theoretically and computationally, we propose a design rule that
minimizes stress in the bulk of the material. This design rule can guide the
realization of stress focusing and fracturing protection in real materials.

##  I Introduction

How materials respond to applied forces has been a crucial question to
understand since ancient times when stable and robust machinery was desired.
In the past decades, the role geometry plays in mechanical responses has begun
to grasp attention, and a set of materials named mechanical metamaterials, in
which the mechanical properties are controlled by geometry, has been developed
[1, 2, 3, 4, 5]. These innovative designs are insensitive to the composition
and the scale because of the consistency in geometry, and they offer many
interesting engineering opportunities due to their unusual features such as
negative Poisson’s ratio [6] and negative stiffness [7]. Topological
mechanical metamaterials (TMMs) represent a special class of metamaterials in
which the mechanical responses are governed by the topology of the bulk phonon
bands, giving them robust properties immune to imperfections [8, 9, 10, 11].
TMMs can be realized in many different physical systems [12, 13, 14, 15, 16,
17, 18, 19], opening up opportunities to study fascinating physical phenomena
and to engineer intriguing functionalities [20, 21, 22, 23]. A feature of
particular interest is the focusing of stress at the domain wall when two
topologically distinct domains are connected together [20, 22]. The fact that
stresses in the structure are focused at the domain wall provides protection
against the development of fracturing sites in the structure [22].

In the linear approximation, it is understood that a lattice can have both
zero energy modes (ZMs), which allow the lattice to deform without energy
cost, and states of self-stress (SSSs), which allow the lattice to carry
stress without causing net force on the sites [8]. ZMs and SSSs are related by
the Maxwell-Calladine index, which is calculated from the number of degrees of
freedom and the number of constraints in the lattice [24, 21]. In topological
Maxwell lattices, ZMs and SSSs can be directed to become exponentially
localized at lattice edges and interfaces by tuning the unit cell geometry
[8]. The SSSs supply a backbone for the lattices in response to external
stresses [25] and offers an understanding of the aforementioned stress
focusing and fracture protection phenomena [22].

Besides the theoretical advancement in TMMs, efforts have been made toward the
engineering and manufacturing of TMMs [12, 13, 14, 17, 15, 16, 18, 19]. One
known difficulty in completing these tasks is the energy cost associated with
hinge bending [12, 26, 17, 15], which is not captured in the theory of TMMs,
where free joints are assumed. The non-zero bending stiffness of the hinges
has diluting effects on the mechanical properties predicted by theory. It
causes an up-shift of the ZM bands from zero frequency, causing hybridization
with bulk bands at long wavelength, which leads to lower contrast between the
predicted hard and soft boundaries of the lattice [14, 27, 11, 28, 16, 29]. It
also intensify the stresses at the joints away from where the SSSs are
localized, which masks the focusing effect of the stress responses [18, 19].
Despite the undesired engineering and manufacturing challenges, the bending
stiffness at the joints offers a new perspective for controlling the
transformations through which the lattice is able to switch between
topological domains [13, 15, 30].

In this paper, we report a systematic study of the effect of weak interactions
beyond the ideal Maxwell connectivity, such as bending stiffness at the hinges
which are ubiquitous in manufactured TMMs, on the topological stress focusing
effect in Maxwell lattice TMMs. Although these weak interactions cause
stresses in the bulk at the hinges, diluting the topological stress focusing,
we find a general design rule for minimizing stress in the bulk by tuning the
unit cell geometry, and we demonstrate this design rule via the agreement
between analytic theory and numerical simulations of mechanical loading on
these lattices. Our results will guide practical designs for Maxwell lattice
TMMs for stress focusing, fracture control, and impact mitigation.

##  II  Theory: topologically protected SSS and additional weak interactions

In this section, we theoretically analyze the effect of weak interactions in
non-ideal lattices on topological stress focusing. We first review fundamental
theories of topological polarization in ideal Maxwell lattices and how
topological SSSs control stress distribution in Maxwell lattices. We then show
how the addition of weak interactions, in terms of next-nearest neighbor (NNN)
bonds, leads to additional SSSs in the bulk that governs the stress
distribution the bulk.

###  II.1 Maxwell lattices, topological polarization, and states of self
stress

For a lattice of $N$ sites with equilibrium positions $\mathbf{r}_{i}$
($i=1,\dots,N$) and $N_{c}$ bonds (here and elsewhere we use bond and spring
interchangeably) in $d$ dimensions, the elastic energy can be generally
written as

$H=\sum_{\alpha\in\text{bonds}}\frac{k_{\alpha}}{2}\quantity(\mathbf{l}_{i_{\alpha}j_{\alpha}}-\mathbf{l}^{0}_{i_{\alpha}j_{\alpha}})^{2}$
(1)

where $i_{\alpha}$ and $j_{\alpha}$ are the indices of the sites that bond
$\alpha$ connects, $\mathbf{l}^{0}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$,
$\mathbf{l}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}+\mathbf{u}_{i}-\mathbf{u}_{j}$,
and $\mathbf{u}_{i}$ is the displacement of the site $i$. Up to quadratic
order in displacements $\mathbf{u}_{i}$, the elastic energy is

$H=\frac{1}{2}\mathbf{u}^{\intercal}\cdot\mathbf{D}\cdot\mathbf{u}$ (2)

where $\mathbf{u}$ is the column vector containing all displacement components
of all sites, ⊺ stands for transpose,
$\mathbf{D}=\mathbf{C}^{\intercal}\cdot\mathbf{k}\cdot\mathbf{C}=\mathbf{Q}\cdot\mathbf{k}\cdot\mathbf{Q}^{\intercal}$
is the dynamical matrix, in which $\mathbf{k}$ is the diagonal matrix with
springs constants as of all bonds as its entries, $\mathbf{C}$ is the
compatibility matrix that relates the $dN$-dimensional site displacement
column vector $\mathbf{u}$ to the $N_{c}$-dimensional column vector of bond
extensions $\mathbf{e}$, and $\mathbf{Q}=\mathbf{C}^{\intercal}$ is the
equilibrium matrix that relates the $N_{c}$-dimensional bond tensions
$\mathbf{t}$ to the $dN$-dimensional column vector of forces on sites [8]:

$\mathbf{C}\cdot\mathbf{u}=\mathbf{e},\ \ \ \ \
\mathbf{Q}\cdot\mathbf{t}=\mathbf{f}.$ (3)

The entries of the $\mathbf{C}$ matrix are determined by the equation for the
extension of the spring $\alpha$

$e_{\alpha}=\hat{l}_{\alpha}\cdot\quantity(\mathbf{u}_{i_{\alpha}}-\mathbf{u}_{j_{\alpha}}),$
(4)

where, once again, where $i_{\alpha}$ and $j_{\alpha}$ are the indices of the
sites that bond $\alpha$ connects, and $\hat{l}_{ij}$ is the unit vector
pointing from site $j$ to site $i$. The null space of $\mathbf{C}$ gives the
ZMs as it represents a set of site displacements that cause no bond
extensions, and the null space of the $\mathbf{Q}$ represents SSSs as it
represents a set of bond tensions that cause no net force on the sites.

For a homogeneous system with translation symmetry, it is easier to work with
the compatibility matrix in Fourier space $\mathbf{C}(\mathbf{q})$, which is
$n_{c}\times nd$ matrix, (where $n_{c}$ is the number of bonds per unit cell
and $n$ is the number of sites per unit cell). As shown in the seminal paper
by Kane and Lubensky [8], for a system with a domain wall between two
different Maxwell lattices, the number of SSSs localized at a the domain wall
can be calculated from the properties of the $\mathbf{C}(\mathbf{q})$ matrices
of each Maxwell lattice; below we describe this briefly and discuss the
consequences for our specific system shown in Fig. 1(a).

From the compatibility matrix $\mathbf{C}(\mathbf{q})$ of a homogeneous
lattice at the Maxwell point (i.e., $n_{c}=nd$) in Fourier space, the
dynamical matrix can be obtained as
$\mathbf{D}(\mathbf{q})=\mathbf{C}^{\dagger}(\mathbf{q})\mathbf{C}(\mathbf{q})$
(where we chose the spring constants and the masses of the particles to be 1
for simplicity). The normal mode frequencies of the lattice
$\omega^{2}(\mathbf{q})$ are the eigenvalues of the Dynamical matrix
$\mathbf{D}(\mathbf{q})$. Kane and Lubensky [8] defined a ‘square root’ of the
dynamical matrix, which in reciprocal space takes the following form:

$\mathcal{H}(\mathbf{q})=\begin{pmatrix}\mathbf{0}&\mathbf{C}^{\dagger}(\mathbf{q})\\\
\mathbf{C}(\mathbf{q})&\mathbf{0}\end{pmatrix}.$ (5)

For every nonzero eigenvalue $\omega^{2}(\mathbf{q})$ of
$\mathbf{D}(\mathbf{q})$, $\mathcal{H}(\mathbf{q})$ has two eigenvalues
$\pm\omega(\mathbf{q})$. The zero modes of $\mathcal{H}(\mathbf{q})$ include
nullspace of $\mathbf{C}(\mathbf{q})$ (ZMs) and nullspace of
$\mathbf{C}^{\dagger}(\mathbf{q})$ (SSSs). Maxwell-Calladine theorem [24, 9]
dictates that the number of ZMs ($n_{0}$) and number of SSSs ($n_{s}$) are
equal ($n_{0}=n_{s}$) for a Maxwell lattice. The matrix
$\mathcal{H}(\mathbf{q})$ has the property that
$S\mathcal{H}(\mathbf{q})S=-\mathcal{H}(\mathbf{q})$, where
$S=\text{Diag}\\{\mathbbm{1},-\mathbbm{1}\\}$. This property is known as
chiral (or sublattice) (anti)symmetry in the literature. Also, it is easy to
check that $\mathcal{H}(\mathbf{q})$ has time reversal symmetry:
$\mathcal{H}(\mathbf{q})=\mathcal{H}^{*}(-\mathbf{q})$, where ∗ is complex
conjugation. These two symmetries put the matrix $\mathcal{H}(\mathbf{q})$ in
BDI class of Altland Zirnbauer (AZ) classification [31]. According to the AZ
periodic table [32] BDI (or AIII if we exclude time-reversal symmetry) class
has a $\mathbbm{Z}$ topological index in 1D. Hence, a topological polarization
vector $\mathbf{R}_{T}$ of a 2D Maxwell lattice can be defined as

$\begin{split}\mathbf{R}_{T}&=\sum_{i}n_{i}\mathbf{a}_{i},\\\
n_{i}&=\frac{1}{2\pi
i}\oint_{\mathbf{q}\rightarrow\mathbf{q}+\mathbf{G}_{i}}d\mathbf{q}\cdot\mathbf{\nabla}_{\mathbf{q}}\log\det\mathbf{C}^{\dagger}(\mathbf{q}),\end{split}$
(6)

where $\mathbf{a}_{i}$ are the lattice vectors, $\mathbf{G}_{i}$ are the
reciprocal lattice vectors such that
$\mathbf{G}_{i}\cdot\mathbf{a}_{j}=2\pi\delta_{ij}$, and
$n_{i}\in\mathbbm{Z}$. Note that this definition of the polarization is only
applicable if the $\mathbf{C}(\mathbf{q})$ is a square matrix with $n_{c}=nd$
or in other words the lattice is Maxwell. Also, we assumed that there are no
bulk ZMs (except the uniform translations); hence, the integers $n_{i}$ do not
depend on $\mathbf{q}$. It was shown in [8] that at a domain wall between two
Maxwell lattices the total number of ZMs at the domain wall at every wave
number $q$ parallel to the domain wall is given by

$\begin{split}n_{dw}&=\frac{1}{2\pi}\mathbf{G}\cdot(\mathbf{R}_{T}^{L}-\mathbf{R}_{T}^{R}),\end{split}$
(7)

where the superscripts $L$ and $R$ denote systems on the left and right of the
domain wall respectively, $\mathbf{G}$ is smallest reciprocal lattice vector
pointing from left to right of the domain wall, and we assumed that the
polarizations for both systems were calculated in the same gauge. Note that if
this number is positive, it counts the number of ZMs, whereas if the number is
negative, its absolute value counts the number of SSSs [8].

Figure 1: Comparison between simulation and analytic results on an ideal (NN
only) Maxwell lattice under a vertical load. (a): An energy minimization
simulation by fixing the bottom boundary and giving a vertical displacement on
all sites at the top boundary with a magnitude $u=10^{-3}$. The left and right
boundaries are free, and sliding along $x$ is allowed for the top and bottom
boundaries. The triangle marks the fixed boundary and the arrows mark the
displaced boundary. The topological polarization vectors $\mathbf{R_{T}}$ for
the left and right domains are shown with arrows. The intensity of the red
color and the thickness of the bonds both increase with the magnitude of
stress on the bonds. The maximum stress in the response $t_{max}=2.6\times
10^{-4}$ in units of force. (b): Unit cells of the left and right domains; the
red circles are the sites inside the unit cell. The black lines show the bonds
in the unit cell. The blue arrows show the lattice vectors. (c-f): 4 states of
self-stress from unit cell analysis obtained by first solving for the four
values of $q_{x}$ from equation $\det\mathbf{Q}(q_{x},q_{y}=0)=0$
corresponding to the unit cell of the right domain, then plugging each of them
into the matrix $\mathbf{Q}(q_{x},q_{y}=0)$ and then obtaining the nullspace
in each case. (c), (d), (e), (f) correspond to $q_{x}=2\pi\mp 13.8174i,2\pi\mp
0.739487i,0,0$, respectively. The SSSs corresponding to the left domain are
mirror images of (c-f) w.r.t. a vertical mirror. (g): The lowest 4 eigenvalues
of $\tilde{\mathbf{D}}(q_{y})=\mathbf{C}(q_{y})\mathbf{Q}(q_{y})$ of the
super-cell shown in (h-k). The supercell is Bloch-periodic in $y$-direction
and open in $x$-direction. At large values of $|q_{y}|$ all 4 eigenvalues are
zero. At $q_{y}=0$ two of them are lifted to nonzero eigenvalue, the other two
remain at zero. (h-k): The lowest 4 eigenstates of
$\tilde{\mathbf{D}}(q_{y}=0)=\mathbf{C}(q_{y})\mathbf{Q}(q_{y})$ from
supercell analysis. The red and blue represent tension and compression for the
SSS respectively, and the thicknesses are proportional to the magnitudes on
the bonds. The modes in (h) and (i) have zero eigenvalue, thus they are SSSs,
but the modes in (j) and (k) are elevated to nonzero eigenvalue, hence they
are not SSSs.

We use the example shown in Fig. 1(a) throughout this section to illustrate
our points. Fig. 1(a) contains a domain wall between two deformed kagome
lattices which are mirror images of each other. A unit cell of the lattice on
the left (see Fig. 1(b)) consists of three sites with third and first sites
having positions $\mathbf{r}_{3}^{L}-\mathbf{r}_{2}^{L}=(-0.510684,-0.42569)$
and $\mathbf{r}_{1}^{L}-\mathbf{r}_{2}^{L}=(0.0105642,-0.297959)$ w.r.t. the
second site, whereas in the case of the lattice on the right, the third and
first sites have positions
$\mathbf{r}_{3}^{R}-\mathbf{r}_{2}^{R}=(-0.0105642,-0.297959)$ and
$\mathbf{r}_{1}^{R}-\mathbf{r}_{2}^{R}=(0.510684,-0.42569)$ w.r.t. the second
site. Both lattices have lattice vectors $\mathbf{a}_{1}=(1,0)$ and
$\mathbf{a}_{2}=(-0.5,0.7)$ (hence reciprocal lattice vectors
$\mathbf{G}_{1}=2\pi(1,1/1.4)$, $\mathbf{G}_{2}=2\pi(0,1/0.7)$). The three
sites are connected by bonds to each other inside the unit cell. Furthermore,
site 1 is connected to site 3 in the next unit cell in $\mathbf{a}_{1}$
direction, site 2 is connected to site 1 in the next unit cell in
$\mathbf{a}_{2}$ direction, site 3 is connected to site 2 in the next unit
cell in $-\mathbf{a}_{1}-\mathbf{a}_{2}$ direction (note that there are other
ways of choosing bonds, but this choice of bonds makes the unit cell
symmetric. This unit cell contains 3 sites, leading to 6 degrees of freedom
(DOFs). With the 6 nearest neighbor (NN) bonds, this defines an ideal Maxwell
lattice.

The polarization vectors of the two lattices are
$\mathbf{R}_{T}^{L}=-\mathbf{a}_{1}$, $\mathbf{R}_{T}^{R}=\mathbf{a}_{1}$, and
the smallest reciprocal lattice vector in the $x$-direction (perpendicular to
the domain wall) is $\mathbf{G}=2\mathbf{G}_{1}-\mathbf{G}_{2}$. Consequently
$n_{dw}=-4$, which implies that there are 4 SSSs at each wave number $q_{y}$
parallel to the domain wall ($y$-direction). The unit cell geometry is chosen
arbitrarily. As in the theory discussed above, as long as the topological
polarization is fixed, small changes in the geometry only results in change of
the localization length of topological ZMs and SSSs, without altering their
site of localization, as guaranteed by the topological protection.

Since under affine uniaxial strain in $y$-direction the strain energy only
depends on the projection of the strain onto the $q_{y}=0$ SSSs [9, 20, 22],
we only focus on those. One can obtain all the information on the SSSs at
$q_{y}=0$ from the equilibrium matrix $\mathbf{Q}(q_{y}=0)$. First, the
equation
$\det\mathbf{Q}(q_{x},q_{y}=0)=-0.352071-0.567495\cos(q_{x}/2)+0.919566\cos(q_{x})=0$
has four solutions $q_{x}=2\pi\mp 13.8174i,2\pi\mp 0.739487i,0,0$ for the left
and the right configurations, respectively. Plugging these solutions into
$\mathbf{Q}(q_{x},q_{y}=0)$ and solving for the nullspace of the matrix in
each case, we obtain the SSSs shown in Fig. 1(c-f). From this, the inverse
decay lengths of the four SSSs can be read off as
$\kappa=13.8174,0.739487,0,0$, which means that the first two are localized at
the domain wall, whereas the last two are extended throughout the system. It
is also worth noting that the real parts of the first two solutions are
$(\Re{q_{x}},q_{y})=(2\pi,0)=\mathbf{G}_{1}-\frac{1}{2}\mathbf{G}_{2}\equiv\frac{1}{2}\mathbf{G}_{2}$
(modulo a reciprocal lattice vector). The reason for this is that the
periodicity along of system along the domain wall is
$2\mathbf{a}_{2}+\mathbf{a}_{1}=(0,1.4)$. If we take a unit cell on either
side that is compatible with the domain wall, its lattice vectors would be
$\mathbf{a}_{1}^{\prime}=\mathbf{a}_{1}$ and
$\mathbf{a}_{2}^{\prime}=2\mathbf{a}_{2}+\mathbf{a}_{1}$ (this unit cell has
twice the area of the primitive unit cell, and has 6 sites inside it), and
hence the reciprocal lattice vectors are
$\mathbf{G}_{1}^{\prime}=\mathbf{G}_{1}-\frac{1}{2}\mathbf{G}_{2}$ and
$\mathbf{G}_{2}^{\prime}=\frac{1}{2}\mathbf{G}_{2}$, which means the point
$\mathbf{q}=\frac{1}{2}\mathbf{G}_{2}$ in reciprocal space “folds to”
$\mathbf{q}=\mathbf{0}$ for this choice of the domain wall compatible unit
cell.

###  II.2 Domain walls, finite system, stress projection, and localization

In the last subsection we mentioned that among the 4 SSSs at $q_{y}=0$, two
are extended in $x$-direction. Consequently, the last two are only really SSSs
for a system that is infinitely large in $x$-direction. Hence, in a finite
system with domain wall such as that in Fig. 1(a), these two modes are not
SSSs. To see this, we do a supercell analysis on the structure shown in Fig.
1(b-e) with Bloch-periodic boundary condition in $y$-direction, but open
boundary condition in $x$-direction. For this system, we plot the eigespectrum
of the matrix
$\tilde{\mathbf{D}}(q_{y})=\mathbf{C}(q_{y})\mathbf{C}^{\dagger}(q_{y})=\mathbf{C}(q_{y})\mathbf{Q}(q_{y})$
(the supersymmetric partner of the dynamical matrix $\mathbf{D}(q_{y})$ [8])
in Fig. 1(g). The eigenmodes of this matrix with zero eigenvalue are the SSSs.
As can be seen from the figure, at large $q_{y}$, there are 4 zero modes.
however at (and near) $q_{y}=0$, two of the modes get lifted to nonzero
eigenvalues, hence they are not SSSs. In Fig. 1(b-e), the lowest 4 eigenmodes
of $\tilde{\mathbf{D}}(q_{y}=0)$ are plotted. The ones in Fig. 1(h-i) are
localized at the domain wall with inverse decay lengths $13.8174\text{ and
}0.739487$ respectively, and they correspond to zero eigenvalues, hence they
remain SSSs. The two in Fig. 1(j-k) are extended in $x$-direction, and they
have nonzero eigenvalues, these are not SSSs for finite system size with open
boundary.

To understand the role of these two SSSs localized at the domain wall in
localizing stress distribution under uniaxial stress in $y$-direction, we
adopt the formulation in Refs. [9, 20, 22], by projecting the affine strain to
the linear space of these SSSs. We first note that under strain
$\epsilon_{yy}$, $\alpha^{\tiny\text{th}}$ bond gets stretched by amount
$e_{\text{aff},\alpha}=\epsilon_{yy}\sin^{2}\theta_{\alpha}$, where
$\theta_{\alpha}$ is the angle between $\alpha^{\tiny\text{th}}$ bond and the
$x$-direction. Defining $\mathbf{e}_{\text{aff}}$ as the vector of the affine
elongations $e_{\text{aff},\alpha}$, one can show that vector of the bond
tensions $\mathbf{t}$ in equilibrium to be

$\mathbf{t}=\sum_{i}\left(\left[\left(\mathbf{K}^{-1}\right)_{ss}\right]^{-1}\mathbf{e}_{\text{aff},s}\right)_{i}\mathbf{s}_{i}\xrightarrow{k=1}\mathbf{t}=\sum_{i}e_{\text{aff},s}^{i}\mathbf{s}_{i},$
(8)

where $\mathbf{s}_{i}$ are the SSSs at $q_{y}=0$, $\mathbf{e}_{\text{aff},s}$
is vector of projections
$e_{\text{aff},s}^{i}=\mathbf{e}_{\text{aff}}\cdot\mathbf{s}_{i}$,
$\left(\mathbf{K}^{-1}\right)_{ss}=[\mathbf{s}]^{T}\mathbf{K}^{-1}[\mathbf{s}]$
where $\mathbf{K}$ is the matrix of bond stiffness and $[\mathbf{s}]$ is the
matrix consisting of column vectors $\mathbf{s}_{i}$. The expression for bond
tensions simplifies to $\mathbf{t}=\sum_{i}e_{\text{aff},s}^{i}\mathbf{s}_{i}$
when all bond stiffnesses are equal and taken to be 1. Now, the localized SSS
in Fig. 1(i) has both tensions and compressions (red and blue colored bonds,
respectively); hence the projection $e_{\text{aff},s}^{i}$ of an affine
uniaxial strain on it is negligible compared to the projection of the affine
uniaxial strain in the SSS in Fig. 1(h). This is why bond tension distribution
(or stress distribution) of the finite system under affine strain
$\epsilon_{yy}$ should closely resemble the SSS in Fig. 1(h); since the SSS in
Fig. 1(b) is highly localized at the domain wall (recall that its inverse
decay length is $13.8174$), the stress distribution should be highly localized
at the domain wall.

We perform a numerical simulation on a finite system shown in Fig. 1(a) (see
Sec. III figure caption for details) to compare with the super-cell analysis.
The stress distribution in the simulation agrees well with the SSS shown in
Fig. 1(h), as we expect from the projection discussed above.

###  II.3 Addition of weak next nearest neighbor bonds

Figure 2: Comparison between simulation and the analytical results on a non-
ideal (with NNN bonds) Maxwell lattice under a vertical load. (a): An energy
minimization simulation as the one described in Fig. 1(a), but with the
inclusion of NNN spring constant $\kappa=10^{-4}$. The maximum tension on the
NN bonds is $t_{max}^{NN}=2.55\times 10^{-4}$ in the units of force, and
$t_{max}^{NNN}=1.56\times 10^{-7}$ on the NNN bonds. The NN bonds are labeled
in red, while the NNN bonds are represented in blue. The intensity of the
color increases with the stress magnitude. (b,c): The total stress on bonds
across columns in the lattice in linear (b) and log (c) scales, at different
$\kappa$. The peak in the middle correspond to the SSS domain wall. (d): Unit
cells of the left and right domains; the dark red circles are the sites inside
the unit cell, the light red circles are the sites at the adjacent unit cell.
The black solid and dashed lines show the NN and NNN bonds in the unit cell,
respectively. (e): Projection of affine strain $\epsilon_{yy}$ onto the $6$
new SSSs at $\mathbf{q}=\mathbf{0}$ due to the additional NNN bonds for a bulk
system with unit cell the same as on the right of the domain wall. Projection
of affine strain $\epsilon_{yy}$ onto the $6$ new SSSs for the bulk system on
te left of the domain wall is just a mirror image of this under a vertical
mirror. Here, we take $\kappa/k=1$. (f): Bond tensions in the supercell under
affine strain $\epsilon_{yy}$ obtained using Eq. (8) by projecting the strain
onto the $474$ SSSs at $q_{y}=0$ for the supercell (we take Bloch-periodicity
in the $y$-direction and open boundary condition in the $x$-direction). Away
from the the domain wall the bond tensions closely resemble that in (e). Here,
we take $\kappa/k=1$ to exaggerate the bond tensions in the NNN bonds. (g):
Stress in unit cells plotted across the columns from the supercell analysis at
different $\kappa$, which is analogous to (c) from the numerical simulation.

Real specimen of manufactured TMMs are never exactly at the Maxwell point,
which would require perfectly flexible hinges. To quantify the effect of hinge
bending stiffness on the stress localization in topologically polarized
Maxwell TMMs, we add weak NNN bonds to the Maxwell lattices, which to leading
order represents the angular bending stiffness at the hinges, and study the
stress distribution analytically and computationally. In particular, the
extensions/compressions of these NNN bonds are equivalent to bending at the
hinges/sites in real metmaterials, and stress concentration due to bending at
the hinges in the bulk of the lattice has been observed computationally [18].

With the addition of the $6$ NNN bonds per unit cell (Fig. 2(d)), the number
of constraints in each unit cell is 6 more than the number of degrees of
freedom. As a consequence, for the homogeneous system on either side of the
domains the matrix
$\mathbf{Q}(\mathbf{q})\equiv\mathbf{C}^{\dagger}(\mathbf{q})$ is $6\times
12$, which implies that there are $12-6=6$ extra SSSs at every $\mathbf{q}$
than before adding NNN bonds. Although these SSSs arise due to NNN bonds, the
tensions associated with these SSSs lie on both NN and NNN bonds.

Similar to the discussion at Eq. (8), the projection of the affine strain to
these $\mathbf{q}=\mathbf{0}$ SSSs determines the stress distribution in the
bulk of the lattice far from the domain wall. At $\mathbf{q}=\mathbf{0}$,
there are two more SSSs (in total $6+2=8$) due to the fact that there are two
ZMs that are uniform translations. However, there were already 2 SSSs for the
system with only NN bonds (the two extended ones shown in Fig. 1), which are
still SSSs of the system with both NN and NNN bonds. These two SSSs do not
remain SSSs for a finite system as discussed earlier. Hence, we are only
concerned with the 6 new SSSs at $\mathbf{q}=\mathbf{0}$ due to the NNN bonds
and the projection of the affine strain onto them for the calculation of
stress in the bulk. If the projections of the affine strain onto these SSSs
are nonzero, there will be finite extension/compression of bonds away from the
domain wall as we show in Fig. 2(e). It is worth noting that the $\mathbf{K}$
matrix is no longer an identity matrix here: the NN bonds stiffness $k$ and
NNN stiffness $k_{NNN}$ are typically different.  Furthermore, since these new
SSSs are due to the additional NNN bonds, the NNN bonds would get
extended/compressed under strain, which implies stress concentration at the
hinges in the bulk away from the domain wall in real metamaterials.

To show the fact that a linear combination of the 6 bulk SSSs at
$\mathbf{q}=\mathbf{0}$ are responsible for the NNN bond
stretching/compression under affine uniaxial strain $\epsilon_{yy}$, we first
perform a supercell analysis. To this end, we create a supercell with a domain
wall between two homogeneous systems (as shown in Fig. 2(f)) with $N_{0}=20$
domain wall compatible unit cells (note that a domain wall compatible unit
cell has twice as many DOFs as a primitive unit cell) on each side of the
domain wall. In total there are $6\times 2N_{0}+1=241$ sites (482 DOFs) and
$954$ bonds. We apply the open boundary condition in the $x$-direction and
Bloch-periodic boundary condition in the $y$-direction. Due to the excess of
bonds ($954-482=472$), at every wave number $q_{y}$ parallel to the domain
wall there are 472 bulk SSSs. At $q=0$, the number of SSSs is $474$ (the extra
two are due to translation symmetry). Among these $474$ SSSs, $2$ are
localized at the domain wall (these were already present in the nearest
neighbor model without the NNN bonds), and $472$ are bulk SSSs appearing due
to excess bonds at every unit cell. In principle, the affine uniaxial strain
can have projections onto all of these $474$ SSSs. However, if our argument in
the last paragraph is correct, among these $472$ extra SSSs, only $6$ of them
carry stress deep inside the bulk away from the domain wall and the edges. To
confirm this, we plot the tension obtained using Eq. (8) the projection of
affine strain onto the 474 SSSs from the supercell analysis in Fig. 2(f), the
same obtained from the projection of affine strain onto the 6 bulk SSSs
obtained from unit cell analysis at $\mathbf{q}=0$ discussed in the last
paragraph in Fig. 2(e) (in both cases, we take the bond stiffnesses for both
NN and NNN bonds to be 1); the bond tension patterns in these two figures
match quite well deep inside the bulk.

Furthermore, we perform a numerical simulation on a finite system with NNN
bonds, as shown in Fig. 2(a), where an axial displacement is applied (see Sec.
III for details). As can be seen from the color scale in Fig. 2(a), there is
stress localization near domain wall on the NN bonds (which is very close to
what we found in Fig. 1(a) for the NN only lattice), but NNN bonds are
stretched/compressed more or less homogeneously in the bulk away from the
domain wall, which is what we predicted above using the SSS argument. We also
plot the total tension on bonds in each column in the lattice in Fig. 2(b-c)
from the full simulation for different values of the ratio $\kappa/k$ of the
NNN bond stiffness and NN bond stiffness. Notice that for each value of
$\kappa/k$, the total bond tension is almost independent of the column number
in the bulk away from the domain wall and the edges. Moreover, this total bond
tension increases linearly with increasing $\kappa$ for small $\kappa/k$. This
can be understood using Eq. (8) in the following way. The bond tension in
equilibrium is given by
$\mathbf{t}=\sum_{i}\left(\left[\left(\mathbf{K}^{-1}\right)_{ss}\right]^{-1}\mathbf{e}_{\text{aff},s}\right)_{i}\mathbf{s}_{i}$.
The homogeneous bond tensions in the bulk away from the domain wall is due to
SSSs that appear due to the extra NNN bonds. For these SSSs,
$\left[\left(\mathbf{K}^{-1}\right)_{ss}\right]^{-1}\propto\kappa$ for small
$\kappa/k$; hence $|\mathbf{t}|\propto\kappa$. To confirm this, we plot
$\mathbf{t}=\sum_{i}\left(\left[\left(\mathbf{K}^{-1}\right)_{ss}\right]^{-1}\mathbf{e}_{\text{aff},s}\right)_{i}\mathbf{s}_{i}$
for different values of $\kappa/k$ under affine strain $\epsilon_{yy}$ in the
supercell analysis in Fig. 2(g). We find the same behavior of the total bond
tensions in each column as a function of $\kappa$ as we found in full
numerical simulation in Fig. 2(c).

##  III Simulation models

In this section we discuss two simulation models we use to study stress
distribution in topological Maxwell TMMs. The first model is the NNN bond
model which we introduced in the analytic theory, where NNN bonds are added to
the lattice capturing the bending stiffness at the hinges.

The second model is the angular spring model (AS). The model associates an
energy cost to the angular change from the rest angles with quadratic form

$H_{\theta}=\sum_{i,jk}\frac{\kappa_{\theta_{i,jk}}}{2}\quantity(\theta_{i,jk}-\theta^{0}_{i,jk})^{2}$
(9)

where $\kappa_{\theta_{i,jk}}$ and $\theta^{0}_{i,jk}$ are the angular spring
constant and the rest angle for angle $\theta_{i,jk}$ formed by the central
site $i$ and its neighboring sites $j,k$. From the AS energy function, the
forces on the particles $i,j,k$ from each angle

$F_{l\alpha}=-\sum_{\\{i,jk\\}}\kappa_{\theta_{i,jk}}\quantity(\theta_{i,jk}-\theta_{i,jk}^{0})\frac{\partial\theta_{i,jk}}{\partial
x_{l\alpha}},$ (10)

where $\alpha$ denotes $x,y$ directions, can be calculated by taking the first
derivative of Eq. (9) to reflect the bending stresses in the lattice. A
detailed full expression of the AS forces and the angle convention we use can
be found in the App. A. In our simulations, angular springs are added at all
sites for bond pairs that do not belong to the same triangle (the bending of
in-triangle bond pairs already cost energy in the NN-only lattice), and
$\kappa_{\theta_{i,jk}}$ is taken to be a constant for all angular springs.

The NNN and the AS model capture very similar effects where additional
constraints are added to the ideal Maxwell lattice to constrain DOFs that were
floppy. In deed, to leading order, the AS stiffness is related to the NNN
spring stiffness $\kappa$ by a geometric parameter which is determined by the
geometry of the unit cell, as it shows how much an NNN spring changes in
length as the corresponding angle changes:

$\kappa=-\frac{l_{ij}l_{ik}\sin\theta}{\sqrt{l^{2}_{ij}+l^{2}_{ik}-2l_{ij}l_{ik}\cos\theta}}\kappa_{\theta}.$
(11)

In the simulation, an axial displacement of $10^{-3}\ll 1$ in $y$-direction is
applied to all sites at the top boundary (which corresponds to strain
$\epsilon_{yy}\approx 7\times 10^{-5}$), and sites on the bottom boundary are
pinned in the $\hat{y}$-direction, while the left and right boundaries are
left open, and the sites on the top and bottom boundaries are allowed to slide
in $x$. We minimize the energy of the lattice under this boundary condition
(for details see App. C). Results from this protocol applied on the NNN model
and the AS model are shown in Fig. 2 and 3, respectively.

In both cases, we observe similar effects where at small $\kappa$, sharp
stress focusing occurs at the SSS domain wall. As $\kappa$ increases, this
topological focusing is diluted, with stress in the bulk growing linearly with
$\kappa$. This feature, as well as the detailed profile of stress in the unit
cell in the bulk, are well captured by our analytic theory based on unit cell
analysis and supercell analysis.

Figure 3: AS model under a vertical load (a): An energy minimization
simulation as the one described in Fig. 1(a), but with AS with spring constant
$\kappa_{\theta}=10^{-4}$ instead of NNN bonds. The maximum tension on the NN
bonds is $t_{max}^{NN}=2.59\times 10^{-4}$ in the units of force, and
$t_{max}^{AS}=2.48\times 10^{-7}$ on the AS. The NN bonds are labeled in red,
while the AS are represented in blue (as a short line connecting the centers
of the bond pairs). The intensity of the color increases with the stress
magnitude. (b,c): The total stress on bonds across columns in the lattice in
linear (b) and log (c) scales at different AS stiffness $\kappa$, where the
peak in the middle correspond to the stress focusing at the domain wall.

Although the AS model is closer to the actual bending stiffness in real TMMs,
its multi-body nature makes it more complicated to study analytically. Thus,
we chose to use the NNN model, which produces very similar results, in our
analytic theory study in Sec. II.

##  IV  Design rule for minimizing stress in the bulk

Figure 4: Generalized kagome lattice. (a): Unit cell of the kagome lattice,
the black solid lines represent the NN bonds and the gray dashed one the NNN
bonds. The red sites are the sites in this unit cell and the green and yellow
ones represent the same sites in the neighboring unit cells. (b): Lattice
vectors $\mathbf{a}_{1}$, $\mathbf{a}_{2}$ and
$\mathbf{a}_{3}=-\mathbf{a}_{1}-\mathbf{a}_{2}$. The vectors $\mathbf{p}_{1}$,
$\mathbf{p}_{2}$ and $\mathbf{p}_{3}$ are unit vectors rotated
counterclockwise by $\pi/2$ from $\mathbf{a}_{1}$, $\mathbf{a}_{2}$ and
$\mathbf{a}_{3}$ , respectively. (c): The deformation of the kagome lattice
when $x_{1}>0,x_{2}=x_{3}=0$. Site 1 and site 2 move by amount
$-\sqrt{3}\mathbf{p}_{1}$ and $x_{1}\mathbf{a}_{3}$ under this deformation.
The blue and black lines are the bonds before and after the deformation
respectively, whereas the solid and dashed lines show NN and NNN bonds
respectively. (d): Ternary phase diagram showing topological polarization
$\mathbf{R}_{T}$ for a particular value of $x_{1}+x_{2}+x_{3}>0$. The red
lines correspond to one of the $x_{i}$ being zero. The configurations inside
the red triangle are unpolarized. The polarization $\mathbf{R}_{T}$ in the six
sectors surrounding the red triangle are shown by blue arrows.

In Sec. II, we showed that addition of NNN bonds (equivalent to having bending
energy cost at the hinges) causes stress concentration on NNN bonds away from
the domain wall inside the bulk. Now, a natural question arises: how do we
minimize the stress in the NNN bonds in the bulk away from the domain wall,
and maintain strong stress focusing at SSS domain walls in these non-ideal
Maxwell lattices? This question is crucial as stress in the bulk at carried by
bending load on the thin hinges, making the lattice prone to damage.  The key
observation in this direction was made in Sec. IIC, that _the stress
concentration on the NNN bonds in the bulk is due to the projection of the
affine strain $\epsilon_{yy}$ onto the additional 6 bulk SSSs at
$\mathbf{q}=\mathbf{0}$._ This gives the design rule:

#### Design-rule:

To minimize the stress concentration on the NNN bonds, we need to design a
deformed kagome lattice such that projection of $\epsilon_{yy}$ on the $6$
extra SSSs at $\mathbf{q}=\mathbf{0}$ due to the additional NNN bonds is
minimized keeping the topological polarization intact.

In this section, we first review the parametrization of generalized kagome
lattice introduced in Ref. [8], and then search for the bulk-stress minimizing
configuration in this parameter space using the design rule mentioned above.

###  IV.1  Generalized Kagome Lattice

In the regular (undeformed) kagome lattice the positions of the three sites in
the unit cell are $\mathbf{r}_{1}=\mathbf{a}_{1}/2$,
$\mathbf{r}_{2}=-\mathbf{a}_{3}/2$, $\mathbf{r}_{3}=\mathbf{0}_{3}$ where
$\mathbf{a}_{n}=\left(\cos(2\pi(n-1)/3),\sin(2\pi(n-1)/3)\right)$ are the
lattice vectors (see Fig. 4(b)). The unit cell, the sites and the NN and NNN
bonds for regular kagome lattice are shown in Fig. 4(a). Now, this lattice can
be deformed in infinitely many ways. However, while deforming the kagome
lattice, if we keep the lattice vectors and one of the sites fixed, there are
only four degrees of freedom to deform the lattice which are the $x$ and $y$
coordinates of two of the sites, since we can always fix one of the sites.
Kane and Lubensky [8] used a different parametrization of these deformations
(of course, the number of of degrees of freedom of deformation is still 4), in
which the expression of the topological polarization as a function of the
parameters become simple. In this parametrization, the sites of the deformed
kagome lattice are given by

$\mathbf{R}_{i}=\mathbf{r}_{i}-\sqrt{3}x_{i}\mathbf{p}_{i}+x_{i-1}\mathbf{a}_{i+1}+\frac{z}{\sqrt{3}}\mathbf{p}_{i-1},$
(12)

where $\mathbf{p}_{n}=(-\sin(2\pi(n-1)/3),\cos(2\pi(n-1)/3))$ are unit vectors
rotated counterclockwise from $\mathbf{a}_{n}$ by angle $\pi/2$. The effect of
$x_{1}>0,x_{2}=x_{3}=0$ is shown in Fig. 4(c), the effect of the $x_{2}$ and
$x_{3}$ are similar. Parameter $z$ uniformly dilates the triangle formed by
the three sites. In this parametrization the, the expression for the
topological polarization is

$\mathbf{R}_{T}=\frac{1}{2}\sum^{3}_{i=1}\text{sgn}(x_{i})\mathbf{a}_{i},$
(13)

which is shown in the ternary phase diagram in Fig. 4(d). Since we want the
domain wall to be vertical (along $y$-direction) and the number of SSSs at
every $q_{y}$ to be $|n_{dw}|=4$, we want $\mathbf{R}_{T}$ on the left domain
to be $\mathbf{R}_{T}=-\mathbf{a}_{1}$ which implies
$(x_{1}<0,x_{2}>0,x_{3}>0)$, which is the middle left region in the ternary
phase diagram.

###  IV.2 Minimizing bulk stress

Figure 5: Minimization of stress on the NNN bonds in the bulk. (a): Colormap
of the 2-norm of NNN bond tensions in the $(x_{1},x_{2},x_{3})$ ternary phase
diagram for $x_{1}+x_{2}+x_{3}=0.15$ and $z=0$ obtained using the equation
$\mathbf{t}=\sum_{i=1}^{6}\left(\left[\left(\mathbf{K}^{-1}\right)_{ss}\right]^{-1}\mathbf{e}_{\text{aff},s}\right)_{i}\mathbf{s}_{i}$
where $\mathbf{s}_{i}$ are the 6 extra SSSs at $\mathbf{q}=\mathbf{0}$ due to
the additional NNN bonds. Here we used $\kappa/k=0.01$ and $\epsilon_{yy}=1$.
At each red line the coordinate perpendicular to the red line is 0 (for
example, $x_{1}=0$ along the vertical red line). The white arrows show the
topological polarization in different regions of the parameter space. (b): The
unit cells corresponding to the cyan and green dots in panel (a). The values
of the parameters are $x_{1}=-0.1,x_{2}=0.125,x_{3}=0.125,z=0$ for cyan and
$x_{1}=-0.03,x_{2}=0.09,x_{3}=0.09,z=0$ for green. (c): Stress in unit cells
(for $\kappa/k=0.01$) plotted across the columns from the supercell analysis,
for the configurations shown in panel (b).

Now that we have the design rule and a parametrization of the deformed kagome
lattice where the expression topological polarization is simple, we can do a
parameter sweep, and for each value of the set of parameters calculate the NNN
bond tensions using the equation
$\mathbf{t}=\sum_{i=1}^{6}\left(\left[\left(\mathbf{K}^{-1}\right)_{ss}\right]^{-1}\mathbf{e}_{\text{aff},s}\right)_{i}\mathbf{s}_{i}$
where $\mathbf{s}_{i}$ are the 6 extra SSSs at $\mathbf{q}=\mathbf{0}$ of the
unit cell due to the additional NNN bonds. Since the parameter space is
4-dimensional, visualization of the results is not simple. Hence, we restrict
ourselves to $x_{1}+x_{2}+x_{3}=0.15$ and $z=0$, these constraints make the
parameter space 2-dimensional. This is an arbitrary choice just to demonstrate
the optimization. More thorough parameter sweep can be done following the same
equations for practical designs.

We show the result in Fig. 5(a). For a left domain with
$\mathbf{R}_{T}=-\mathbf{a}_{1}$, we look at the central left region (the unit
cells of the right domain are mirror images of the unit cell of the left
domain). As can be seen from the colormap in Fig. 5(a), the stress in the NNN
bonds decreases as one goes toward left from the vertical red line (along
which $x_{1}=0$). We show the unit cells in Fig. 5(b) at two different points
in the parameter space marked by cyan and green dots in Fig. 5(a). We found
2-norm $|\mathbf{t}|$ of the NNN bond tensions
($|\mathbf{t}|=(\sum_{i=1}^{6}t_{i}^{2})^{1/2}$, where the sum is over the 6
NNN bonds) are $0.0357292$ and $0.0535986$ at the cyan and green points (note
that $\kappa/k=0.01$ and $\epsilon_{yy}=1$). Furthermore, the maximum tensions
in the NNN bonds are $0.0201862$ and $0.030927$ at the cyan and green points.

To verify these predictions, next we create supercells for these two
configurations analogous to those in Figs. 1 and 2 and evaluate the bond
tensions using the formula
$\mathbf{t}=\sum_{i=1}^{474}\left(\left[\left(\mathbf{K}^{-1}\right)_{ss}\right]^{-1}\mathbf{e}_{\text{aff},s}\right)_{i}\mathbf{s}_{i}$,
where $\mathbf{s}_{i}$ are all the 474 SSSs at $q_{y}=0$ of the supercell and
plot their norm per unit cell in Fig. 5(c). Clearly, there is stress
concentration at the domain wall, which is expected due to the two domain wall
localized topological SSSs. Also, due to the extra SSSs from the NNN bonds,
the bond tensions become constant away from the domain wall. However, the bond
tensions for the cyan configuration is lower than the green configuration.
This matches with our prediction from unit cell strain projection analysis in
Fig. 5(a). Also, the maximum NNN bond tensions for cyan and green
configurations are $0.0214706$ and $0.0316642$ respectively, which are in very
good agreement with our unit cell results reported earlier, verifying that the
unit cell SSS projection in non-ideal Maxwell TMMs serves as a efficient and
accurate prediction of stress in the bulk.

This brings us to the main findings of the article. We first argued that the
hinges in real metamaterials can be modeled as NNN bonds. Then using this
model we showed that the fact that there is stress concentration at the hinges
away from the domain wall inside the bulk in real metamaterials can be
attributed to the extra bulk SSSs coming from the extra NNN bonds. Using this
insight, we gave a design rule for Maxwell lattice to minimize the hinge
stresses away from the domain wall.

In principle, one can go even further to the left in the ternary phase diagram
in Fig. 5(a) for even smaller stress projection onto the NNN bonds; however,
at some point the system becomes self-intersecting (the NN bonds are physical,
they can not intersect each other), which is unphysical.

##  V Discussion

In this paper, we study the effect of weak interactions beyond the Maxwell
constraint-counting on the topological stress focusing effect in Maxwell TMMs,
and provide a design rule to minimize the diluting effect of these weak
interactions on topological stress focusing.

Ideal Maxwell lattices are known to focus stress when the system is under
external load, where stress exponentially localizes at SSS domain walls,
leaving the bulk virtually stress-free and protected against damage. However,
in real manufactured TMMs, weak interactions are always present, and are known
to lead to stress in the bulk, causing bending loads at hinges, making them
prone to damage.

We show that such stress in the bulk is well captured by a simple unit-cell
calculation based on the projection of the load to the SSSs caused by these
additional constraints in the bulk, and propose a design rule to optimize the
stress focusing of these TMMs, minimizing such stress in the bulk by tuning
the unit cell geometry. We test our theory using simulations of real-space
non-ideal Maxwell TMMs under load using two types of models based on NNN bonds
and AS springs, and show the agreement between the analytic theory and the
simulation.

This design rule provides an efficient and accurate way to find unit cell
geometries that optimize TMMs for stress focusing. Such stress focusing effect
can potentially be applied to many technologies from fracture protection,
stimuli response, and impact mitigation.

###### Acknowledgements.

This work is supported by the National Science Foundation NSF-CMMI-2026794
(H.L., S.G., and X.M.) and Office of Naval Research MURI N00014-20-1-2479
(H.L., S.S., N.A. E.S. and X.M.).

##  Appendix A Angular Spring Model Energy and Forces

Energy from the angular spring differs from usual springlike interactions
because it involves three sites, instead of being a pairwise interaction. We
represent these triples with $\\{i,jk\\}$ where $i$ represents the central
particle and $j$ and $k$ represent the side particles. A diagram of the angle
composed by central particle $0$ and particles $x$ and $y$ is shown in Fig. 6.

The energy of bending can be written as a springlike interaction on this
angle:

$U=\sum_{\\{i,jk\\}}\frac{\kappa_{i,jk}}{2}\left(\theta_{i,jk}-\theta_{i,jk}^{0}\right)^{2}$
(14)

where $\theta_{i,jk}$ and $\theta_{i,jk}^{0}$ are the current and rest angles
associated with the triple $\\{i,jk\\}$ and $\kappa_{i,jk}$ is the angular
stiffness of this bending bond.

This angle can be measured from the positions of the particles with
$\theta_{i,jk}=\arccos(\hat{r}_{ij}\cdot\hat{r}_{ik})$ where $\hat{r}_{ij}$
and $\hat{r}_{ik}$ are the unit vectors of the two bond vectors that form the
angle $\theta_{i,jk}$. As discussed, this suffers from a limited range of $0$
to $\pi$. We can investigate reflex angles via utilization of the cross
product as well as the dot product:

$\theta_{i,jk}=\pi+\text{sgn}\left(\overrightarrow{r}_{ij}\times\overrightarrow{r}_{ik}\right)\left(\arccos(\hat{r}_{ij}\cdot\hat{r}_{ik})-\pi\right).$
(15)

This measure has a range from $0$ to $2\pi$. Therefore,
$\left(\theta_{i,jk}-\theta_{i,jk}^{0}\right)$ ranges from
$-\theta_{i,jk}^{0}$ to $\left(2\pi-\theta_{i,jk}^{0}\right)$. It is also
important to note that $\theta_{i,jk}=2\pi-\theta_{i,kj}$.

The force on each particle from the AS via the first derivative of the energy
is:

$F_{l\alpha}=-\sum_{\\{i,jk\\}}\kappa_{i,jk}\left(\theta_{i,jk}-\theta_{i,jk}^{0}\right)\frac{\partial\theta_{i,jk}}{\partial
u_{l\alpha}}.$ (16)

This derivative of the angle with respect to particle motion is given by

$\frac{\partial\theta_{i,jk}}{\partial
u_{l\alpha}}=-\frac{\textrm{sgn}\left(\overrightarrow{r}_{ij}\times\overrightarrow{r}_{ik}\right)}{\sqrt{1-(\hat{r}_{ij}\cdot\hat{r}_{ik})^{2}}}\left(\frac{\delta_{jl}-\delta_{il}}{r_{ij}}(\hat{r}_{ik,\alpha}-(\hat{r}_{ij}\cdot\hat{r}_{ik})\hat{r}_{ij,\alpha})+\frac{\delta_{kl}-\delta_{il}}{r_{ik}}(\hat{r}_{ij,\alpha}-(\hat{r}_{ij}\cdot\hat{r}_{ik})\hat{r}_{ik,\alpha})\right).$
(17)

Figure 6: A bond pair in the AS model and its relation to NNN model. Angle
$\theta$ (shown in red) is the bond angle between bonds $\mathbf{r}_{ij}$ and
$\mathbf{r}_{ik}$. The NNN bond $\mathbf{r}_{jk}$ (shown in dashed blue line)
is a leading order approximation of the AS. When the site $j$ moves by a small
displacement $\mathbf{u}_{j}$, both the bond angle and the NNN bond length
change.

The force derivative $\derivative{u}$is taken with respect to the
displacements of the sites $\mathbf{u}$ away from their equilibrium positions,
and $\alpha$ labels the $x$ and $y$ components of the force. We do not need to
concern ourselves with the derivative of the sign function as this only gives
a delta function when the angle is $\pi$ due to the singularity.

Note that there is no force if the bond pair is at the rest angle. However, if
the bond is bent, there is a force not only on the external particles but also
on the central particle, which is only symmetric if the distances between
particles are the same.

We can simplify this equation using orthogonal vectors in 2D. For a given
vector $\overrightarrow{V}$, the orthogonal vector is provided by

$\overrightarrow{V}^{\perp}=\left(\begin{matrix}-V_{y}\\\
V_{x}\end{matrix}\right).$ (18)

The force on any particle due to the angular bonds is therefore given by

$F_{l\alpha}=\sum_{\\{i,jk\\}}\kappa_{i,jk}\left(\theta_{i,jk}-\theta_{i,jk}^{0}\right)\left(\frac{(\delta_{jl}-\delta_{il})\hat{r}_{ij,\alpha}^{\perp}}{r_{ij}}-\frac{(\delta_{kl}-\delta_{il})\hat{r}_{ik,\alpha}^{\perp}}{r_{ik}}\right).$
(19)

##  Appendix B $\kappa_{\theta}$ to $\kappa$ conversion

By examining an angle’s geometric relations, we can build a conversion
relation between the AS spring constant and the NNN spring constant. We do so
by examining the conversion between the small angle change $\delta\theta$ of
the angle $\theta$ between bonds $\mathbf{r}_{ij}$ and $\mathbf{r}_{ik}$ and
the change of length $\delta l$ of bond $\mathbf{r}_{jk}$ between sites $j$
and $k$, which is the NNN springs bond that the AS of angle $\theta$ is
approximated as. A diagram of the geometry is shown in Fig. 6.

As the Fig. 6 shows, the angle $\theta$ is changed by $\delta\theta$ when the
site $j$ goes through a small displacement $\mathbf{u}_{j}$, however, as the
displacement is small, we approximate the length
$r_{ij}=\absolutevalue{\mathbf{r}_{ij}}$ remains unchanged as the angle is
rotated. Then the length change $\delta l$ for the bond length
$r_{jk}=\absolutevalue{\mathbf{r}_{jk}}$ can be written as

$\delta
l=\sqrt{r^{2}_{ij}+r^{2}_{ik}-2r_{ij}r_{ik}\cos\quantity(\theta)}-\sqrt{r^{2}_{ij}+r^{2}_{ik}-2r_{ij}r_{ik}\cos\quantity(\theta-\delta\theta)}$
(20)

Taking the deformation angle $\delta\theta$ to be small and expand Eq. (20) to
linear order in $\delta\theta$, we have

$\delta
l\approx\frac{r_{ij}r_{ik}\sin\quantity(\theta)}{\sqrt{r^{2}_{ij}+r^{2}_{ik}-2r_{ij}r_{ik}\cos\quantity(\theta)}}\delta\theta,$
(21)

where $r_{ij}=\absolutevalue{\mathbf{r}_{ij}}$ and
$r_{ik}=\absolutevalue{\mathbf{r}_{ik}}$ are the lengths of the two bonds
which the angle $\theta$ is composed of. And by this geometric relation
between the length change of the NNN spring $\delta l$ and the angle change
$\delta\theta$, we conclude with the relation Eq. (11) between the AS constant
$\kappa_{\theta}$ and the NNN spring constant $\kappa$ in Sec. III.

##  Appendix C Numerical Methods

Two numerical methods were used in the simulations. For the simulations with
NN and NNN springs, a simple gradient descent method was used. This
implementation calculates the force of the nodes in the system and displaces
them following the force gradient until an error tolerance of $10^{-7}$ is
reached. As for the simulation including AS, a _Fast Inertial Relaxation
Engine (FIRE)_ algorithm is used to handle the non-linear multi-particle
interactions and to accelerate the computation. For the _FIRE_ algorithm, not
only the force gradient $\mathbf{F}$ is calculated, but also the velocity
$\mathbf{v}$ and its dot product with the force $P=\mathbf{F}\cdot\mathbf{v}$.
Then the direction and the speed of the descent are controlled by $P$.

## References

  * [1] Liu Z, Zhang X, Mao Y, Zhu Y, Yang Z, Chan CT, et al.  Locally resonant sonic materials.  Science. 2000;289(5485):1734-6.
  * [2] Huber SD.  Topological mechanics.  Nature Physics. 2016 6;12:621-3.  Available from: https://www.nature.com/articles/nphys3801.
  * [3] Bertoldi K, Vitelli V, Christensen J, Hecke MV.  Flexible mechanical metamaterials.  Nature Reviews Materials. 2017 10;2:1-11.  Available from: https://www.nature.com/articles/natrevmats201766.
  * [4] Yu X, Zhou J, Liang H, Jiang Z, Wu L.  Mechanical metamaterials associated with stiffness, rigidity and compressibility: A brief review.  Progress in Materials Science. 2018 5;94:114-73.
  * [5] Xin L, Siyuan Y, Harry L, Minghui L, Yanfeng C.  Topological mechanical metamaterials: A brief review.  Current Opinion in Solid State and Materials Science. 2020 10;24:100853.
  * [6] Lakes RS.  Negative-Poisson’s-Ratio Materials: Auxetic Solids.  https://doiorg/101146/annurev-matsci-070616-124118. 2017 7;47:63-81.  Available from: https://www.annualreviews.org/doi/abs/10.1146/annurev-matsci-070616-124118.
  * [7] Nicolaou ZG, Motter AE.  Mechanical metamaterials with negative compressibility transitions.  Nature Materials 2012 11:7. 2012 5;11:608-13.  Available from: https://www.nature.com/articles/nmat3331.
  * [8] Kane CL, Lubensky TC.  Topological boundary modes in isostatic lattices.  Nature Physics. 2013 12;10:39-45.  Available from: https://www.nature.com/articles/nphys2835.
  * [9] Lubensky TC, Kane CL, Mao X, Souslov A, Sun K.  Phonons and elasticity in critically coordinated lattices.  Reports on Progress in Physics. 2015 6;78:073901.
  * [10] Mao X, Lubensky TC.  Maxwell Lattices and Topological Mechanics.  Annual Review of Condensed Matter Physics. 2018 3;9:413-33.
  * [11] Sun K, Mao X.  Continuum Theory for Topological Edge Soft Modes.  Physical Review Letters. 2020 5;124:207601.  Available from: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.124.207601.
  * [12] Bilal OR, Süsstrunk R, Daraio C, Huber SD.  Intrinsically polar elastic metamaterials.  Advanced Materials. 2017;29(26):1700540.
  * [13] Rocklin DZ, Zhou S, Sun K, Mao X.  Transformable topological mechanical metamaterials.  Nature Communications. 2017 1;8:1-9.  Available from: https://www.nature.com/articles/ncomms14201.
  * [14] Ma J, Zhou D, Sun K, Mao X, Gonella S.  Edge Modes and Asymmetric Wave Transport in Topological Lattices: Experimental Characterization at Finite Frequencies.  Physical Review Letters. 2018 8;121:094301.  Available from: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.121.094301.
  * [15] Xiu H, Liu H, Poli A, Wan G, Sun K, Arruda EM, et al.  Topological transformability and reprogrammability of multistable mechanical metamaterials.  Proceedings of the National Academy of Sciences. 2022;119(52):e2211725119.
  * [16] Jolly JC, Jin B, Jin L, Lee Y, Xie T, Gonella S, et al.  Soft mechanical metamaterials with transformable topology protected by stress caching.  Advanced Science. 2023:2302475.
  * [17] Bergne A, Baardink G, Loukaides EG, Souslov A.  Scalable 3D printing for topological mechanical metamaterials.  Extreme Mechanics Letters. 2022 nov;57:101911.  Available from: https://doi.org/10.1016%2Fj.eml.2022.101911.
  * [18] Widstrand C, Hu C, Mao X, Labuz J, Gonella S.  Stress focusing and damage protection in topological Maxwell metamaterials.  International Journal of Solids and Structures. 2023;274:112268.
  * [19] Widstrand C, Mao X, Gonella S.  Robustness of stress focusing in soft lattices under topology-switching deformation.  arXiv preprint arXiv:231014972. 2023.
  * [20] Paulose J, Meeussen AS, Vitelli V.  Selective buckling via states of self-stress in topological metamaterials.  Proceedings of the National Academy of Sciences of the United States of America. 2015 6;112:7639-44.  Available from: https://www.pnas.org/doi/abs/10.1073/pnas.1502939112.
  * [21] Sun K, Souslov A, Mao X, Lubensky TC.  Surface phonons, elastic response, and conformal invariance in twisted kagome lattices.  Proceedings of the National Academy of Sciences of the United States of America. 2012 7;109:12369-74.  Available from: https://www.pnas.org/doi/abs/10.1073/pnas.1119941109.
  * [22] Zhang L, Mao X.  Fracturing of topological Maxwell lattices.  New Journal of Physics. 2018 6;20:063034.
  * [23] Zhou D, Ma J, Sun K, Gonella S, Mao X.  Switchable phonon diodes using nonlinear topological Maxwell lattices.  Physical Review B. 2020 3;101:104106.  Available from: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.101.104106.
  * [24] Calladine CR.  Buckminster Fuller’s ”Tensegrity” structures and Clerk Maxwell’s rules for the construction of stiff frames.  International Journal of Solids and Structures. 1978 1;14:161-72.
  * [25] Chen BGG, Upadhyaya N, Vitelli V.  Nonlinear conduction via solitons in a topological mechanical insulator.  Proceedings of the National Academy of Sciences of the United States of America. 2014 9;111:13004-9.  Available from: https://www.pnas.org/doi/abs/10.1073/pnas.1405969111.
  * [26] Guo XF, Ma L.  Periodic topological lattice with different indentation hardness on opposite surfaces.  Materials and Design. 2019 10;180:107953.
  * [27] Stenull O, Lubensky T.  Signatures of topological phonons in superisostatic lattices.  Physical review letters. 2019;122(24):248002.
  * [28] Saremi A, Rocklin Z.  Topological elasticity of flexible structures.  Physical Review X. 2020;10(1):011052.
  * [29] Charara M, McInerney J, Sun K, Mao X, Gonella S.  Omnimodal topological polarization of bilayer networks: Analysis in the Maxwell limit and experiments on a 3D-printed prototype.  Proceedings of the National Academy of Sciences. 2022;119(40):e2208051119.
  * [30] Xiu H, Frankel I, Liu H, Qian K, Sarkar S, MacNider B, et al.  Synthetically non-Hermitian nonlinear wave-like behavior in a topological mechanical metamaterial.  Proceedings of the National Academy of Sciences. 2023;120(18):e2217928120.  Available from: https://www.pnas.org/doi/abs/10.1073/pnas.2217928120.
  * [31] Altland A, Zirnbauer MR.  Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures.  Phys Rev B. 1997 Jan;55:1142-61.  Available from: https://link.aps.org/doi/10.1103/PhysRevB.55.1142.
  * [32] Chiu CK, Teo JCY, Schnyder AP, Ryu S.  Classification of topological quantum matter with symmetries.  Rev Mod Phys. 2016 Aug;88:035005.  Available from: https://link.aps.org/doi/10.1103/RevModPhys.88.035005.


††thanks<EMAIL_ADDRESS>
# Eigenvalues restricted by Lyapunov exponent of eigenstates

Tong Liu  Department of Applied Physics, School of Science, Nanjing University
of Posts and Telecommunications, Nanjing 210003, China  Xu Xia  Chern
Institute of Mathematics and LPMC, Nankai University, Tianjin 300071, China

(August 27, 2024)

###### Abstract

We point out that the Lyapunov exponent of the eigenstate places restrictions
on the eigenvalue. Consequently, with regard to non-Hermitian systems, even
without any symmetry, the non-conservative Hamiltonians can exhibit real
spectra as long as Lyapunov exponents of eigenstates inhibit imaginary parts
of eigenvalues. Our findings open up a new route to study non-Hermitian
physics.

###### pacs:

71.23.An, 71.23.Ft, 05.70.Jk

##  I Introduction.

In physics, the observability of physical quantities demands that the measured
eigenvalue must be a real number Shankar . Therefore, quantum mechanics
generally assumes that the observable quantity is the Hermitian operator.
However, with the breakthrough of quantum theory and the development of
experimental technology Lindblad ; Photonic1 ; ultracold1 ; Bender1 , it has
been proved that non-Hermitian systems can also host real eigenvalues and be
observed experimentally Xiao . Currently, the study of real energy spectrum
mainly focuses on two classes: (i) real energy spectrum of the symmetry
protection, such as parity-time ($\mathcal{PT}$) symmetry Bender2 ; (ii) real
energy spectrum induced by boundary conditions, such as open boundary
condition in non-Hermitian skin effect Yao1 ; Yao2 .

Thus, a question arises naturally: is there other physical mechanism for
generating the real energy spectrum of non-Hermitian systems? In this work, we
point out that the Lyapunov exponent (LE) of the eigenstate places
restrictions on the eigenvalue, hence a non-Hermitian system can exhibit the
real energy spectrum, independent of any symmetry and boundary conditions.
Previous studies Simon ; Avila  have shown that the absolutely continuous
spectrum of the Schrödinger operators is completely determined by the LE,
namely the absolutely continuous spectrum corresponds to the LE being zero.
Thus, LE of eigenstates and eigenvalues establish a certain relationship. More
explicitly, $\gamma$ of the eigenstate $\psi$ can always be expressed as a
function of its eigenvalue $E$, namely

$\gamma\sim f(E),$ (1)

while the Lyapunov exponent $\gamma$ has only two possibilities, $\gamma=0$
for extended/critical eigenstates in Fig. 1(a) and $\gamma>0$ for localized
eigenstates in Fig. 1(b). Hence $\gamma$ imposes restrictions on $E$, which
may lead to the imaginary part $\Im(E)=0$, non-Hermitian systems can exhibit
real energy spectrum.

To illustrate this point in detail, we analytically calculate a solvable
model, namely the following non-Hermitian difference equation,

$E\psi_{n}=\psi_{n+1}+\psi_{n-1}+Vi\cot(\pi\alpha n)\psi_{n},$ (2)

where $V$ is the complex potential strength, $E$ is the eigenvalue of systems,
and $\psi_{n}$ is the amplitude of wave function at the $n$th lattice. We
choose to unitize the nearest-neighbor hopping amplitude and a typical choice
for irrational parameter $\alpha$ is $\alpha=(\sqrt{5}-1)/2$. Firstly, the
complex potential $i\cot(\pi\alpha n)$ doesn’t satisfy $V(n)=V^{*}(-n)$ in the
discrete lattice, hence the model is independent of $\mathcal{PT}$ symmetry;
secondly, the left and right hopping amplitudes are the same, hence the model
is independent of non-Hermitian skin effect. Next we will calculate the
explicit expressions of $\gamma$ and $E$, and give the condition that $E$ are
real numbers.

##  II Lyapunov exponent and real energy.

Figure 1: (Color online) (a)Cartoon diagram of an extended/critical
eigenstate($\gamma=0$) versus $\Im(E)=0$. (b)Cartoon diagram of an localized
eigenstate($\gamma>0$) versus $\Im(E)\neq 0$.

---

Figure 2: (Color online) (a) $\operatorname{Re}(E)$ versus LE as a function of
the potential strength $V$. It clearly shows that when
$\operatorname{Re}(E)\in[V-2,2-V]$, $\gamma(E)=0$; otherwise $\gamma(E)>0$.
(b) $\Im(E)$ versus LE as a function of the potential strength $V$. It clearly
shows that when $\Im(E)=0$, $\gamma(E)=0$; otherwise $\gamma(E)>0$.

The definition of LE has many different forms, the following form is adopted
for a semi-infinite chain Thouless ,

$\gamma=-\lim_{n\rightarrow\infty}\frac{1}{n-n_{0}}\log\left|\frac{\psi_{n}}{\psi_{n_{0}}}\right|,$
(3)

where $\psi_{n_{0}}$ denotes the amplitude of wave function at the truncation
of left end , $\psi_{n}$ denotes the amplitude of wave function at the
infinity right end. Since Eq. (2) is a second-order difference equation,
$\left|\frac{\psi_{n}}{\psi_{n_{0}}}\right|$ is a matrix form, which is very
difficult to solve. Utilizing Avila’s global theory Avila , we can get the
explicit expression of LE,

$\displaystyle\gamma(E)=\max\\{$
$\displaystyle\operatorname{arcosh}\frac{\left|E+V+2\right|+\left|E+V-2\right|}{4},$
(4)
$\displaystyle\operatorname{arcosh}\frac{\left|E-V+2\right|+\left|E-V-2\right|}{4}\\}.$

Now, we get the explicit relation between $\gamma$ and $E$, a concrete example
of Eq. (1). To visually illustrate the above result, a plot of $\gamma$ versus
$E$ is shown in Fig. 2. Make $E=\operatorname{Re}(E)+i\Im(E)$, when $V\leq 2$,
$\gamma(E)=0$ corresponds to the region $E\in[V-2,2-V]$, while $\gamma(E)>0$
corresponds to the region $\operatorname{Re}(E)=0$ and
$\Im(E)\in\mathbb{R}^{*}$. When $V>2$, only $\gamma(E)>0$ exists, and $E$
belongs to the region $\operatorname{Re}(E)=0$ and $\Im(E)\in\mathbb{R}$. Here
we need to emphasize that Eq. (4) only gives the value range of $E$ when
$\gamma(E)=0$, and cannot prove that $E$ exactly corresponds to the eigenvalue
of the system.

Hence we need prove the eigenvalues of Eq. (2) are exactly the value range of
$E$ derived from $\gamma(E)$. We fist introduce the Fourier transformation,

$f_{\theta}=\frac{1}{\sqrt{L}}\sum_{n=1}^{L}e^{i\theta n}\psi_{n},$ (5)

thus the dual equation of Eq. (2) is given as,

$\displaystyle[\cos(\theta+2\pi\alpha)+V/2-E/2]f(\theta+2\pi\alpha)$ (6)
$\displaystyle=[-\cos(\theta-2\pi\alpha)+V/2+E/2]f(\theta-2\pi\alpha).$

By Sarnark’s method Sarnak , from Eq. (6), we can demonstrate that there are
two kinds of eigenvalues of the system: (i) pure real numbers, namely
$E\in[V-2,2-V]$; (ii) pure imaginary numbers, the value range of imaginary
part is the whole real axis. Consequently, we demonstrates that when
$\gamma(E)=0$, the non-Hermitian system hosts the real energy spectrum
$E\in[V-2,2-V]$, independent of any symmetry and boundary conditions, this is
the central innovation of our work.

##  III Numerical verification

To support the analytical result given above, we now present the numerical
verification, namely directly diagonalize Eq. (2) to find the eigenvalues and
eigenstates. With regard to the disordered system, the property of eigenstates
not only can be characterized by LE, but also can be measured by the inverse
participation ratio (IPR) conveniently Kohmoto . For any given normalized
eigenstate, the corresponding IPR is defined as
$\text{IPR}=\sum_{n=1}^{L}\left|\psi_{n}\right|^{4}$, which measures the
inverse of the number of sites being occupied by particles. The IPR of an
extended state scales like $L^{-1}$ which becomes zero in the large $L$ limit,
just as $\gamma=0$ in Fig. 1(a). While for a localized state, the IPR is
finite even in the large $L$ limit, just as $\gamma>0$ in Fig. 1(b). In Fig. 3
we show the diagram of the eigenvalues $\operatorname{Re}(E)$ versus the
complex potential strength $V$. The red eigenvalue curves denote the pure real
energy spectrum hosting critical eigenstates, and the black circle dots denote
the pure imaginary energy spectrum hosting localized eigenstates. These
numerical results are completely consistent with the theoretical results
derived by LE, which confirms the correctness of our theory.

Figure 3: (Color online) The real part of eigenvalues $\operatorname{Re}(E)$
of Eq. (2) as a function of $V$. The total number of sites is set to be
$L=610$. As shown in the figure, when $V\leq 2$, the system host the real
energy spectrum $E\in[V-2,2-V]$. In the inset, we plot the typical probability
density, as expected, they are the critical state (red) corresponding to the
pure real energy spectrum and the localized state (black) corresponding to the
pure imaginary energy spectrum, respectively.

##  IV Conclusion.

In this work, we uncover a new class of physical systems with pure real energy
spectrum, which is different from the known physical mechanism, namely
independent of any symmetry and boundary conditions. Our research shows that
this theory does not depend on any specific physical system, such as optics,
cold atoms or classical circuits. As long as the LE of the eigenstate is
determined, the eigenvalue of the system may have a pure real energy spectrum,
which means that it has observable effects. $\gamma=0$ leads to the real
eigenvalue in most cases, however, we need to emphasize that $\gamma=0$ is not
a necessary condition for the eigenvalue to be a real number, and $\gamma>0$
may also produce a real eigenvalue, as long as $\Im(E)=0$. Therefore, there
are still many academic gaps to be filled.

###### Acknowledgements.

This work was supported by the Natural Science Foundation of Jiangsu Province
(Grant No. BK20200737), NUPTSF (Grants No. NY220090 and No. NY220208), the
National Nature Science Foundation of China (Grant No. 12074064), and the
Innovation Research Project of Jiangsu Province (Grant No. JSSCBS20210521).
X.X. is supported by Nankai Zhide Foundation.

## References

  * (1) R. Shankar, Principles of Quantum Mechanics (Plenum Press, 1994).
  * (2) G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * (3) L. Feng, R. El-Ganainy, and L. Ge, Non-Hermitian photonics based on parity time symmetry, Nat. Photonics 11, 752 (2017)
  * (4) J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Observation of parity-time symmetry breaking transitions in a dissipative Floquet system of ultracold atoms, Nat. Commun. 10, 855 (2019).
  * (5) L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Non-Hermitian bulk-boundary correspondence in quantum dynamics, Nat. Phys. 16, 761 (2020)
  * (6) C. M. Bender and S. Boettcher, Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * (7) C. M. Bender, Making sense of non-Hermitian Hamiltonians, Rep. Prog. Phys. 70, 947 (2007).
  * (8) S. Yao and Z. Wang, Edge states and topological invariants of non-Hermitian systems, Phys. Rev. Lett. 121, 086803 (2018).
  * (9) S. Yao, F. Song, and Z. Wang, Non-Hermitian chern bands, Phys. Rev. Lett. 121, 136802 (2018).
  * (10) H. L. Cycon, R. G. Froese, W. Kirsch, and B. Simon, Schrödinger operators: With application to quantum mechanics and global geometry[M]. Springer, 2009.
  * (11) A. Avila, Global theory of one-frequency Schrödinger operators, Acta. Math. 1, 215, (2015).
  * (12) D. J. Thouless, Phys. Rep. 13, 93 (1974).
  * (13) P. Sarnak, Spectral behavior of quasi periodic potentials, Commun. Math. Phys. 84, 377 (1982).
  * (14) M. Kohmoto, Phys. Rev. Lett 51, 1198 (1983).