Universal Continuum Theory for Topological Edge Soft Modes

TL;DR

This paper develops a continuum topological elasticity theory that generalizes the Maxwell-Calladine index theorem, enabling prediction of edge zero modes in continuous media and extending topological mechanics beyond discrete lattices.

Contribution

It introduces a gauge-invariant bulk topological index for continuous elastic media, bridging topological lattice theories and continuum elasticity.

Findings

Defines a gauge-invariant topological index for continuous media.

Predicts the number of edge zero modes from bulk properties.

Extends topological zero modes to media away from the Maxwell point.

Abstract

Topological edge zero modes and states of self stress have been intensively studied in discrete lattices at the Maxwell point, offering robust properties concerning surface and interface stiffness and stress focusing. In this paper we present a continuum topological elasticity theory and generalize the Maxwell-Calladine index theorem to continuous elastic media. This theory not only serves as a macroscopic description for topological Maxwell lattices, but also generalizes this physics to continuous media. We define a \emph{gauge-invariant} bulk topological index, independent of microscopic details such as choice of the unit cell and lattice connectivity. This index directly predicts the number of zero modes on edges, and it naturally extends to media that deviate from the Maxwell point, depicting how topological zero modes turn into topological soft modes.

Equations

$$E=\int \mathbf{dr}\sum _{i=1}^{d(d+1)/2}\frac{{\lambda }_i}{2}$$

$$+O(\partial \partial u{)}^2,$$

$${\mathcal{C}}_{ij}{u}_j=0\text{for}{\lambda }_i>0,$$

$$E=\frac{1}{2}\int \mathbf{dr}{\mathbf{u}}^T{\overleftarrow{\mathcal{C}}}^T\Lambda \overrightarrow{\mathcal{C}}\mathbf{u},$$

$$N_{T/B}({\bar{q}}_x)=\frac{1}{2\pi i}{\oint }_{T/B}d{q}_y\mathrm{tr}\left({\mathcal{C}}^{-1}{\partial }_{q_y}\mathcal{C}\right),$$

$$E=\int \mathbf{dr}$$

$$\displaystyle+\mathcal{M}^{(1)}_{ijk}\partial_{j}\phi_{k}+\ldots\big{)}^{2}+O[(\partial\partial u)^{2},(\partial\phi)^{2}],$$

$$\displaystyle E=\frac{1}{2}\int\mathbf{dr}\sum_{i=1}^{d(d+1)/2}\bigg{[}\lambda_{i}\left(\mathcal{M}^{(0)}_{ijk}\partial_{j}u_{k}+\mathcal{M}^{(1)}_{ijk\ell}\partial_{j}\partial_{k}u_{\ell}+\mathcal{M}^{(2)}_{ijk\ell m}\partial_{j}\partial_{k}\partial_{\ell}u_{m}+\ldots\right)^{2}$$

$$+\sum _{i=1}^{d^2(d+1)/2}{\tilde{\lambda }}_i{\left({\tilde{\mathcal{M}}}_{ijk\ell }^{(1)}{\partial }_j{\partial }_k{u}_{\ell }+{\tilde{\mathcal{M}}}_{ijk\ell m}^{(2)}{\partial }_j{\partial }_k{\partial }_{\ell }{u}_m+\dots \right)}^2$$

$$+\sum _{i=1}^{(d+2)(d+1)d^2/6}{\tilde{\tilde{\lambda }}}_i{\left({\widetilde{\tilde{\mathcal{M}}}}_{ijk\ell m}^{(2)}{\partial }_j{\partial }_k{\partial }_{\ell }{u}_m+\dots \right)}^2$$

$$\displaystyle+\ldots\bigg{]},$$

$$f(x,y)=a(x+dy{)}^2+e{y}^2$$

$$E=\frac{1}{2}\int \mathbf{dr}\left(\begin{array}{cccc}{\mathbf{u}}^{(1)}& {\mathbf{u}}^{(2)}& {\mathbf{u}}^{(3)}& \dots \end{array}\right)\left(\begin{array}{cccc}{\mathcal{K}}^{(1,1)}& {\mathcal{K}}^{(1,2)}& {\mathcal{K}}^{(1,3)}& \dots \\ {\mathcal{K}}^{(2,1)}& {\mathcal{K}}^{(2,2)}& {\mathcal{K}}^{(2,3)}& \dots \\ {\mathcal{K}}^{(3,1)}& {\mathcal{K}}^{(3,2)}& {\mathcal{K}}^{(3,3)}& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right)\left(\begin{array}{c}{\mathbf{u}}^{(1)}\\ {\mathbf{u}}^{(2)}\\ {\mathbf{u}}^{(3)}\\ ⋮\end{array}\right)$$

$$\left(\begin{array}{cccc}{\mathcal{K}}^{(1,1)}& {\mathcal{K}}^{(1,2)}& {\mathcal{K}}^{(1,3)}& \dots \\ {\mathcal{K}}^{(2,1)}& {\mathcal{K}}^{(2,2)}& {\mathcal{K}}^{(2,3)}& \dots \\ {\mathcal{K}}^{(3,1)}& {\mathcal{K}}^{(3,2)}& {\mathcal{K}}^{(3,3)}& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right)=\left(\begin{array}{cccc}{\mathcal{L}}^{(1,1)})& 0& 0& \dots \\ {\mathcal{L}}^{(2,1)})& {\mathcal{L}}^{(2,2)}& 0& \dots \\ {\mathcal{L}}^{(3,1)})& {\mathcal{L}}^{(3,2)}& {\mathcal{L}}^{(3,3)}& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right)\left(\begin{array}{cccc}{\mathcal{D}}^{(1,1)}& 0& 0& \dots \\ 0& {\mathcal{D}}^{(2,2)}& 0& \dots \\ 0& 0& {\mathcal{D}}^{(3,3)}& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right){\left(\begin{array}{cccc}{\mathcal{L}}^{(1,1)}& 0& 0& \dots \\ {\mathcal{L}}^{(2,1)}& {\mathcal{L}}^{(2,2)}& 0& \dots \\ {\mathcal{L}}^{(3,1)}& {\mathcal{L}}^{(3,2)}& {\mathcal{L}}^{(3,3)}& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right)}^T$$

$$u_i(x,y)=a_{ix}x+a_{iy}y+\frac{1}{2}(a_{ixx}{x}^2+2a_{ixy}xy+a_{iyy}{y}^2)+\dots$$

$$E=\frac{1}{2}\int \mathbf{dr}{\mathbf{u}}^{(1)}{\mathcal{K}}^{(1,1)}{\mathbf{u}}^{(1)}=\frac{1}{2}\int \mathbf{dr}\sum _{i=1}^{d_c}{\lambda }_i{\left({\mathcal{O}}_{ij}{u}_j^{(1)}\right)}^2$$

$${\mathcal{O}}_{ij}{u}_j^{(1)}={\mathcal{O}}_{ij}{\mathcal{P}}_{jk\ell }{\partial }_k{u}_{\ell }=0\text{for}{\lambda }_i>0$$

$$E=\frac{1}{2}\int \mathbf{dr}({\mathcal{K}}_{ij}^{(1,1)}{u}_i^{(1)}{u}_j^{(1)}+2{\mathcal{K}}_{ij}^{(1,2)}{u}_i^{(1)}{u}_j^{(2)}+{\mathcal{K}}_{ij}^{(2,2)}{u}_i^{(2)}{u}_j^{(2)})$$

$$E=\frac{1}{2}\int \mathbf{dr}\left[\sum _{i=1}^{d_c}{\lambda }_i({\mathcal{O}}_{ij}{u}_j^{(1)}{)}^2+2({\mathcal{K}}_{kl}^{(1,2)}{\mathcal{O}}_{ki}^T{u}_l^{(2)})({\mathcal{O}}_{ij}{u}_j^{(1)})+{\mathcal{K}}_{ij}^{(2)}{u}_i^{(2)}{u}_j^{(2)}\right]$$

$${\lambda }_i({\mathcal{O}}_{ij}{u}_j^{(1)}{)}^2+2({\mathcal{K}}_{kl}^{(1,2)}{\mathcal{O}}_{ki}^T{u}_l^{(2)})({\mathcal{O}}_{ij}{u}_j^{(1)})={\lambda }_i{\left({\mathcal{O}}_{ij}{u}_j^{(1)}+\frac{{\mathcal{K}}_{kl}^{(1,2)}{\mathcal{O}}_{ki}^T}{{\lambda }_i}{u}_l^{(2)}\right)}^2-\frac{{\mathcal{K}}_{mn}^{(1,2)}{\mathcal{O}}_{mi}^T{\mathcal{K}}_{kl}^{(1,2)}{\mathcal{O}}_{ki}^T}{{\lambda }_i}{u}_l^{(2)}{u}_n^{(2)}.$$

$$E=\frac{1}{2}\int \mathbf{dr}\left[\sum _{{\lambda }_i>0}{\lambda }_i({\mathcal{O}}_{ij}{u}_j^{(1)}+\frac{{\mathcal{K}}_{kl}^{(1,2)}{\mathcal{O}}_{ki}^T}{{\lambda }_i}{u}_l^{(2)}{)}^2+({\mathcal{K}}_{ln}^{(2)}-\sum _{{\lambda }_i>0}\frac{{\mathcal{K}}_{mn}^{(1,2)}{\mathcal{O}}_{mi}^T{\mathcal{K}}_{kl}^{(1,2)}{\mathcal{O}}_{ki}^T}{{\lambda }_i})u_l^{(2)}{u}_n^{(2)}\right].$$

$$E=\frac{1}{2}\int \mathbf{dr}{\mathbf{u}}^T{\overleftarrow{\mathcal{C}}}^T\Lambda \overrightarrow{\mathcal{C}}\mathbf{u},$$

$${\mathcal{C}}^T\mathbf{w}=0.$$

$$E=\frac{1}{2}\int \mathbf{dr}{\left({\Lambda }^{-1}\mathbf{w}+\mathcal{C}\mathbf{u}\right)}^T\Lambda \left({\Lambda }^{-1}\mathbf{w}+\mathcal{C}\mathbf{u}\right).$$

$${\sigma }^{(m)}=\frac{\delta E}{\delta {\mathbf{u}}^{(m)}}\ne 0$$

$$\displaystyle\mathbf{f}_{\mathbf{u}=0}=-\frac{\delta E}{\delta\mathbf{u}}\bigg{|}_{\mathbf{u}=0}=-\mathcal{C}^{T}\mathbf{w}=0$$

$$N_0=$$

$$N_s=$$

$$N_0-N_s=N_d-N_c,$$

$${\mathbf{u}}^{(\alpha )}(\mathbf{R})=\left(u_1^{(\alpha )}(\mathbf{R}),u_2^{(\alpha )}(\mathbf{R}),\dots ,u_d^{(\alpha )}(\mathbf{R})\right).$$

$$E=\sum _{\mathbf{R}}{E}_{\mathbf{R}},$$

$$E_{\mathbf{R}}=\frac{1}{2}\sum _m{\kappa }_m{\left[{\xi }_m(\mathbf{R})\right]}^2,$$

$${\xi }_m(\mathbf{R})=\sum _{\mathbf{t},i}{\mathcal{M}}_{mi}^{(\alpha )}(\mathbf{t})u_i^{(\alpha )}(\mathbf{R}+\mathbf{t}).$$

Full text

Introduction—Topologically-protected edge zero modes (ZMs) and states of self stress (SSSs) have been intensively studied recently in lattice models at the verge of mechanical instability Kane and Lubensky (2014); Lubensky et al. (2015); Mao and Lubensky (2018), i.e., Maxwell lattices, leading to a rich collection of interesting phenomena including topological solitons Chen et al. (2014), mechanical Weyl modes Rocklin et al. (2016); Stenull et al. (2016), switchable surface stiffness Rocklin et al. (2017), programmable buckling Paulose et al. (2015a) and fracturing processes Zhang and Mao (2018), and beyond Zhou et al. (2018, 2019); Chen et al. (2016); Paulose et al. (2015b); Süsstrunk and Huber (2016); Socolar et al. (2017); Rocklin (2017); Saremi and Rocklin (2018).

These topological phenomena are all based on discrete lattice and network models, and the physical principles that govern them are dependent of microscopic details, e.g., the connectivity of the networks. However, topological edge ZMs and SSSs are found to exist at long wave length (e.g., as opposed to isolated points in the Brillouin zone), where microscopic details become unimportant, which indicates that these topological phenomena may not be exclusive for lattice models. Instead, they may also arise in continuum elasticity. Preliminary continuum elasticity analysis using linear Lubensky et al. (2015) and nonlinear strain tensor Rocklin et al. (2017) identified two different regimes based on the intrinsic soft strain of the lattice: a dilation dominant regime where FMs localize on all edges with inverse decay length $\sim q_{x}$ (where $q_{x}$ is the wavevector parallel to the edge), and a shear dominant regime where low frequency ($\omega\sim q^{2}$ where $q$ is the wavevector) soft modes exist in the bulk. To obtain topological edge modes, which decay as $\kappa\sim O(q_{x}^{2})$, and their topological index, lattice theory with optical bands was necessary. It remains an open question whether a universal description independent of microscopic mechanisms for such topological phenomena exists.

Remarkably, it is shown recently that these topological edge modes still survive when the system deviates from the idealized discrete lattice limit, in specimen obtained from 3D printing Bilal et al. (2017), cutting Ma et al. (2018), and finite-element Ma et al. (2018, 2019) and discrete lattice Stenull and Lubensky (2019) analysis. The original definition of the topological polarization $R_{T}$, based on the Maxwell condition, is no longer valid in these systems. However edge soft modes similar to topological Maxwell lattices are observed, e.g., asymmetric edge stiffness and edge waves. Are these systems truly topological?

In this paper we present a universal continuum theory that can be directly applied to both continuous and discrete elastic media and predict topological soft modes and SSSs. We extend the “Maxwell condition” of the verge of mechanical instability from lattice models, where it takes the form of $\langle z\rangle=2d$ with $\langle z\rangle$ being the mean coordination number and $d$ being the spatial dimension) Jacobs and Thorpe (1995); Liu et al. (2010); Mao et al. (2010); Ellenbroek and Mao (2011); Broedersz et al. (2011); Sun et al. (2012); Mao et al. (2015); Zhang et al. (2015); Zhang and Mao (2016), to the continuum in terms of macroscopic elastic moduli. We name elastic systems that satisfy this macroscopic Maxwell condition “Maxwell media”, which include many Maxwell lattices and beyond. We further define a universal way to measure how far a system is away from this macroscopic Maxwell condition. For systems at or not far away from the Maxwell condition, a topological index $N$ can be defined, which predicts the number of ZMs for each edge. Remarkably, this topological index is insensitive to any microscopic details such as choice of the unit cell and is determined purely by macroscopic elastic constants. We also develop methods to obtain the elastic constants in the continuum theory from first principles, and verify our theory through good agreement with results on lattice models at and close to the Maxwell point.

*Elastic energy with higher-order gradients—*Within linear elasticity, the elastic energy is a quadratic form of displacement-field gradients Landau and Lifshitz (1986). Typically, only the first-order gradients are included, while higher-order derivative terms only contribute to short-distance physics Mindlin (1964) and thus can be ignored at long distance. However, as mentioned above, for Maxwell systems, the leading-order theory is insufficient, and thus a more generic setup with higher-order gradients becomes essential even for the long-distance limit. As proved in the SI (Sec. S-1), the elastic energy of any stress-free elastic medium must take the following form

[TABLE]

where $u$ is the displacement field. This sum contains $d(d+1)/2$ terms where $d$ is the spatial dimension. In standard elasticity theory (using the Voigt notation) Kittel et al. (1976), elastic moduli form a $[d(d+1)/2]\times[d(d+1)/2]$ symmetric matrix. Once this matrix is diagonalized, it gives the leading order terms of Eq. (1), with $\lambda_{i}$ and $\mathcal{M}^{(0)}$ being the eigenvalues and eigenvectors of the Voigt elastic moduli matrix. For topological indices, it is often sufficient to just include one higher order term into the theory, $\mathcal{M}^{(1)}$.

In addition to these $d(d+1)/2$ squares, the elastic energy also contains some higher-order terms [$O(\partial\partial\mathbf{u})^{2}$ and above], which will be ignored for now. These ignored terms are non-negative, and thus for ZMs (e.g., in Maxwell lattices), they must vanish. For a general elastic medium, they lead to subleading corrections to the elastic energy [$\propto O(q^{4})$], which will be calculated below via a variational method.

*Maxwell’s counting in continuum—*In continuum elasticity, the number of degrees of freedom is $d$ for each point in real space, because $u$ contains $d$ components. For each positive $\lambda_{i}>0$, the corresponding square term in Eq. (1) enforces one constraint at this point

[TABLE]

where $\mathcal{C}_{ij}=\mathcal{M}^{(0)}_{ikj}\partial_{k}+\mathcal{M}^{(1)}_{ik\ell j}\partial_{k}\partial_{\ell}+\ldots$ If a deformation field satisfies all the constraints, it is a ZM.

For a rigid solid, all $\lambda_{i}$’s are positive and thus the constraint number is $d(d+1)/2$, which always exceeds $d$ for $d>1$, i.e., the medium is over constrained. This conclusion agrees with the generalized counting argument for continuous media Sun and Mao (2019), and the Janet-Cartan embedding theorem Janet (1926); Cartan (1927). For a system to be at the verge of mechanical instability, some $\lambda_{i}$ may vanish. Most importantly, if the matrix has $d$ positive eigenvalues, while the rest are [math], the number of constraints coincides with that of the degrees of freedom, which will be referred to as the macroscopic Maxwell condition, and this type of systems will be called Maxwell media. In contrast to Maxwell lattices/networks, Maxwell media only require the Maxwell condition at the macroscopic level. In a Maxwell medium the $\mathcal{C}$ matrix defined above becomes a square matrix and this matrix serves the role of the compatibility matrix in Maxwell lattices.

*Zero modes and states of self stress—*Let us consider generic elastic media without assuming the Maxwell condition. The elastic energy [Eq. (1)] can be written as

[TABLE]

where $\Lambda$ is the diagonal matrix composed of positive $\lambda_{i}$’s. Because $\Lambda$ is positive definite by definition, Eq. (3) implies that any ZMs must satisfy $\mathcal{C}\mathbf{u}=0$. It should be emphasized that in Eq. (3), derivatives in $\mathcal{C}$ act to the right on $\mathbf{u}$, while $\mathcal{C}^{T}$ to the left on $\mathbf{u}^{T}$, as marked by the arrows above. This differs from the dynamical matrix $\mathcal{D}=-\overrightarrow{\mathcal{C}}_{\partial\to-\partial}^{T}\Lambda\overrightarrow{\mathcal{C}}$, where all derivatives act to the right. In $\mathcal{D}$, derivatives in $\mathcal{C}^{T}$ are flipped from left to right via integral by parts, which introduces boundary terms. These boundary terms are irrelevant for bulk properties but crucial for physics at the edge. In particular, one needs to distinguish two conditions (a) $\mathcal{C}\mathbf{u}=0$ and (b) $\mathcal{D}\mathbf{u}=0$. From Eq. (3), it is easy to check that (a) is the sufficient and necessary condition of ZMs (i.e., $E=0$). In contrast, (b) is the sufficient and necessary condition for the system to be force free, because force in an elastic medium is $f_{i}=-\delta E/\delta u_{i}=D_{ij}u_{j}$. The condition (b) is weaker than (a), because force-free deformations can be classified into two categories: (1) ZMs ($\mathcal{C}\mathbf{u}=0$) and (2) modes that leave the bulk in force balance but there is force on the boundary ($\mathcal{D}\mathbf{u}=0$ but $\mathcal{C}\mathbf{u}\neq 0$). Furthermore, similar to the lattice theory, $\mathcal{C}^{T}$ plays the role of the equilibrium matrix and its null space contains SSSs of the medium, which yields all possible ways the medium can carry residual stress. Interestingly, the Maxwell-Calladine index theorem that relates the number of ZMs and SSSs still holds [see the SI (Sec. S-1)].

*The topological index—*Similar to lattices, edge ZMs in a Maxwell medium also come from nontrivial topology in the bulk and in this section we define the topological index, which dictates the number of ZMs on each edge. We demonstrate the physics in a 2D system, which can be easily generalized to other spatial dimensions.

Consider an infinitely long strip parallel to the $x$ axis with two open edges (top/bottom). As shown in the SI (Sec. S-3), for such a Maxwell system, zero-energy edge states may arise, with deformation fields $\mathbf{u}=\mathbf{A}e^{iq_{x}x+iq_{y}y}$, where $q_{x}$ is real and $q_{y}$ is complex, whose imaginary part describes how fast the amplitude decays into the bulk. For each fixed $q_{x}=\bar{q}_{x}$, these ZMs are characterized by a bulk topological index.

[TABLE]

where $\mathcal{C}_{ij}=i\mathcal{M}^{(0)}_{ikj}q_{k}-\mathcal{M}^{(1)}_{ik\ell j}q_{k}q_{\ell}+\ldots$ is the Fourier transform of the real-space $\mathcal{C}$ matrix in Eq. (2). For any closed contour, this index is quantized to integer values, because it is the winding number of the phase of $\det\mathcal{C}$. Here, the contour is a semi circle on the complex $q_{y}$ plane as shown in Fig. 1(a). If the semicircle is in the upper (lower) half of the plane, the integral defines $N_{T}$ ($N_{B}$), which gives the number of edge modes for the top (bottom) edge at $q_{x}=\bar{q}_{x}$. Similar to Maxwell lattices, these two indices characterize the topological structure of the bulk bands and their values remain invariant, unless the system undergoes a topological phase transition, which closes and the reopens the bulk phonon gap at $q_{x}=\bar{q}_{x}$.

The radius of the semicircle is the ultraviolet momentum cutoff $\Lambda_{q}$, which arises in (almost) any continuum theory, separating universal properties at long distance ($q<\Lambda_{q}$) from microscopic details at short distance ($q>\Lambda_{q}$). Although a cutoff is needed for its definition, this topological index is cutoff independent, as long as $1/\Lambda_{q}$ is much shorter (longer) than any macroscopic (microscopic) length scales. In our theory, the microscopic length scale can be estimated as $\mathcal{M}^{(1)}/\mathcal{M}^{(0)}$, where $\mathcal{M}^{(0)}$ and $\mathcal{M}^{(1)}$ are the characteristic elastic constants.

*Gauge invariance—*It is worthwhile to emphasize that the topological index defined above differs slightly from the topological polarization defined for Maxwell lattices, due to differences in the integral contours. In Maxwell lattices, it is known that the topological polarization is gauge sensitive, e.g., by choosing a different unit cell, its value may change by an integer. Such a constant shift doesn’t change the topological classification. However, as one tries to connect topological polarization with physical observables, e.g., the number of ZMs at an edge, the gauge dependent part needs to be carefully deducted. In Maxwell lattices, this is usually achieved utilizing a local polarization, a microscopic quantity that depends on local lattice connectivity at the edge Kane and Lubensky (2014); Lubensky et al. (2015); Mao and Lubensky (2018). For a continuum theory which doesn’t rely on microscopic details, a local gauge-fixing is unavailable. Thus a universal gauge-fixing scheme independent of microscopic details becomes necessary. Remarkably, the integral contour defined above offers a gauge independent topological index, which directly gives the number of edge states for each edge.

This can be seen by comparing the contours of the continuum [Eq. (4)] and the lattice theory. In Maxwell lattices, because $q_{y}$ has a periodicity associated with the reciprocal lattice vector (the first Brillouin zone), the integral contour is usually chosen to be along the real $q_{y}$ axis over one period (e.g., from $q_{y}=-\pi$ to $q_{y}=+\pi$). A change of variable $z=e^{iq_{y}}$ is then defined, under which the contour becomes the unit circle in the complex $z$ plane, and components of the topological polarization is given by $n_{i}=\oint_{i}\frac{dz}{2\pi i}\ \operatorname{tr}(\mathcal{C}^{-1}\partial_{z}\mathcal{C})$. As shown in the SI (Sec. S-3), $\operatorname{tr}(\mathcal{C}^{-1}\partial_{z}\mathcal{C})$ is analytic for the entire complex $z$ plane, except for a finite number of poles. Thus the integral above only comes from residue at each pole, and there are two types of poles. Each ZM contributes one pole with residue $1$, and their contributions are gauge invariant. This is the reason why $R_{T}$ can probe ZMs. However, there is a second type of poles at $z=0$ and $z=\infty$. Residue at these two poles are gauge dependent. These two poles are not related with physical ZMs, and they are natural consequences of analytic continuation. Thus, their contribution needs to be removed. In Maxwell lattices this is done via the local polarization $R_{L}$.

In this continuum theory, this gauge fixing is automatically achieved, because these two points are beyond the ultraviolet cutoff [Fig. 1(b)]. This scheme can also be applied to Maxwell lattices, utilizing integral contours in Fig. 1(c).

It needs to be emphasized that this gauge-invariant topological index does not detect edge modes with zero penetration depth, such as a dangling bond. This is because zero penetration depth corresponds to $|\operatorname{Im}q_{y}|=\infty$ (i.e., $z=0$ or $z=\infty$), which has been removed in our definition. Because such a zero-penetration-depth edge mode is purely due to local and microscopic features at the edge, a bulk topological index is not expected to detect them. Instead, they can only be detected via microscopic local probes, such as a local polarization Kane and Lubensky (2014); Lubensky et al. (2015); Mao and Lubensky (2018).

*Optical phonon modes—*In certain soft materials, such as the regular kagome lattice Sun et al. (2012), low-energy optical modes may exist, whose energy can become comparable with acoustic modes. These soft optical modes can be included in the continuum theory to more accurately capture the low-energy physics, especially at intermediate $q$. With soft optical modes, the elastic energy becomes

[TABLE]

where $\phi_{i}$ represents the soft optical modes, and the $m$th order derivatives of $\mathbf{u}$ should be treated as the same order as the $(m-1)$th order derivatives of $\phi_{i}$. The number of degrees of freedom is now increased to $d+n_{O}$ at each point in real space, with $n_{O}$ representing the number of soft optical modes included. The number of eigenvalues ($\lambda_{i}$) is now $d(d+1)/2+n_{O}$. Same as before, if all $\lambda_{i}$’s are positive, the system is over constrained (for $d>1$), and when $d(d-1)/2$ of the eigenvalues vanishes, the Maxwell condition is achieved, at which the topological indices and all the physics can be easily generalized. The same construction [Eq. (5)] also applies to pinned systems, e.g., the pinned 1D chain model discussed in Ref. Kane and Lubensky (2014), where all phonon modes are optical-like, and $\mathbf{u}$ should be removed.

*Comparison with Maxwell lattices—*In the SI (Sec. S-2), we present a first-principle method to determine the coefficients in our continuum theory, based on lattice models or finite-element method (FEM). The key idea in this method is to integrate out high-energy optical modes, which are irrelevant in the long-distance low-energy continuum theory. For the dynamic matrix, this is a standard practice Mao and Lubensky (2011). However, because the topological information is stored in the $\mathcal{C}$ matrix instead, here we need to develop a new approach based on the $\mathcal{C}$ matrix. Remarkably, it can be proved that as far as ZMs are concerned, our procedure of integrating out high-energy modes is an exact method, although the number of degrees of freedom is reduced. Next, we then take the long wavelenghth limit, which leads to the continuum theory.

With this first-principle technique, we can now directly compare the continuum theory with lattice models or FEM. In Fig. 2, we show the decay rates of the edge modes in a topological kagome lattice and compare it with the continuum theory [SI (Sec. S-5)]. Fig. 2(a) shows the continuum theory with only acoustic modes [Eq. (1)], while (b) contains one extra soft optical mode [Eq. (5)]. Both of them provide the correct topological index, fully agreeing with the lattice model, and the decay rate from the continuum theory is asymptotically exact at long wavelength. Remarkably, the continuum theory even captures fine details at intermediate $q_{x}$. Near $\pi/4$, the two edge modes happen to become degenerate with the same $q_{y}$, which is also seen in the continuum theory. The two continuum models in (a) and (b) share identical asymptotic behavior at small $q_{x}$, while including the soft optical branch improves the accuracy at intermediate $q_{x}$, where energy of the soft optical mode becomes comparable with the acoustic ones.

*Maxwell medium beyond Maxwell lattices—*For the deformed kagome lattice, the continuum theory predicts two edge ZMs, agreeing with the microscopic Maxwell counting  Sun et al. (2012); Kane and Lubensky (2014); Lubensky et al. (2015); Mao and Lubensky (2018). For this lattice, the microscopic counting guarantees that the energy of these modes remain zero, even when higher-order gradients beyond Eq. (1) are included. However, for a generic Maxwell medium, terms ignored in Eq. (1) [i.e., $O(\partial\partial u)^{2}$ terms] may push the energy of the edge modes up (i.e., they become soft modes instead of ZMs), even into the bulk bands. As a result, for these systems, in general, after using Eq. (1) to obtain the topological index and edge modes, we need to put these edge-mode candidiates back into the full elastic energy to check whether their energy is indeed below bulk phonon bands. Only if the energy is below the bulk bands, the candidiate becomes as a real edge soft mode.

As shown in the SI (Sec. S-4), this calculation can be easily implemented via a variational method. We utilize each edge-mode candidiate obtained from the topological theory as a trial wavefunction and minimize the full elastic energy. As a variational method, this minimization is exact for ZMs but tends to overestimate the energy of finite-frequency modes  Griffiths (2004). Thus if it reveals an edge soft mode with energy lower than the bulk band, the accurate energy of this mode must be even lower, and thus it must be indeed an edge mode.

Beyond Maxwell medium— Maxwell media require some of the $\lambda_{i}$’s to be exactly [math]. In a real elastic solid all $\lambda_{i}$’s are finite and positive. Here, we can order the eigenvalues from high to low and compute the ratio between the $(d+n_{O}+1)$th largest eigenvalue and the $(d+n_{O})$th, $r_{M}=\lambda_{d+n_{O}+1}/\lambda_{d+n_{O}}$, where $n_{O}$ is the number of soft optical modes included in the continuum theory. $r_{M}$ offers a universal measurement to determine how far away a real system is from the ideal Maxwell limit. Its value is between [math] and $1$, while [math] implies a Maxwell medium. If this quantity is small enough, one can set all $\lambda_{i}$’s with $i>d+n_{O}$ to zero in the leading-order treatment, and perform the topological analysis as shown above. Then, similar to the previous section, we compute the energy of edge-mode candidates using the variational method.

In Fig. 3, we study a deformed kagome lattice with next-nearest-neighbor springs, whose spring constant is 1/1000 times the nearest-neighbor ones. This system is over constrained, but is not far away from the Maxwell point. Thus, it exhibits similar edge modes. We compute the phonon band structure on an infinitely-long strip with two open edges ($200$ unit cells wide), and compare it with predictions from the continuum theory and the variational method, and find good agreement. The continuum theory used in Fig. 3(a) includes only acoustic modes, while (b) contains one extra optical mode. It must be emphasized that in contrast to the lattice numerical solution, which is performed on a strip with open edges, curves from the continuum theory are computed in the wave vector space without any physical edge, but this wave vector space treatment correctly captures the long-distance physics at an open edge.

Acknowledgements— We thank S. Gonella for enlightening discussions which inspired this work. This work was supported by the National Science Foundation under Grant No. NSF-EFRI- 1741618.

Supplementary Information

I Generic elastic theory

In this section, we discuss the elastic theory of a generic elastic medium. We will stay in the linear response regime. However, higher-order derivative terms of the deformation field are allowed and will be included, which play a crucial role in the topological description. The key result of this section is that for any stress-free elastic medium, the elastic energy can always be organized in the following form

[TABLE]

where $u$ is the displacement field (i.e., each point moves $r\to r+u(r)$), and $i,j,k,\ldots$ are the Cartesian indices. Here, we rewrite the elastic energy as sets of square terms with non-negative coefficients ($\lambda_{i}\geq 0$, $\widetilde{\lambda}_{i}\geq 0$, $\widetilde{\widetilde{\lambda}}_{i}\geq 0$, $\ldots$). The first set contains $d(d+1)/2$ terms, each of which is the square of some linear combination of deformation-field gradients. The second sets contains $d^{2}(d+1)/2$ squares. Here, each square term is also composed of a linear combination of deformation-field gradients, but the derivatives starts from the second order (i.e., no first-order derivative terms). The third set contains $(d+2)(d+1)d^{2}/6$ terms, and the linear combination starts from third-order gradient terms, etc.

I.1 A simple example

We first use a simple example to demonstrate the key ideas. Consider a real quadratic form $f(x,y)=ax^{2}+2bxy+cy^{2}$, where $a$, $b$ and $c$ are three real constants, and $x$ and $y$ are two real variables. Here, we show here that if the quadratic form is semi-positive definite, i.e. $f(x,y)\geq 0$ for any real $x$ and $y$, then it can always be written as the sum of two squares

[TABLE]

with $a\geq 0$ and $e\geq 0$.

For a semi-positive definite quadratic form, $a$ must be non-negative, i.e., $a\geq 0$. If $a>0$, Eq. (7) can be achieved by absorbing the cross term $bxy$ into $x^{2}$, which gives $d=b/a$ and $e=c-b^{2}/a$. If $a=0$, the $x^{2}$ term vanishes and thus, we couldn’t use it to absorb the $xy$ term. However, for a semi-positive-definite quadratic form, it is easy to verify that if $a=0$, $b$ must also vanish. Thus, when $a=b=0$, $f(x,y)=cy^{2}$ and Eq. (7) automatically holds as long as we set $e=c$.

It must be emphasized that Eq.(7) only applies to semi-positive-definite quadratic forms. This is very important, because it sets an important constraint to our theory. As will be shown below, our theory only applies to stress-free elastic media, because stress-free is the necessary and sufficient condition to ensure that the quadratic form used in our theory is semi-positive definite.

I.2 Elastic energy with higher order derivative terms

In this section, we use $\mathbf{u}$ to represent the displacement field. All the $m$th-order gradients of the displacement field will be made into a column vector, which will be dubbed as $\mathbf{u}^{(m)}$. Here, the first-order terms are jus the linearized strain tensor (presented in the Voigt notation), e.g., in 2D $\mathbf{u}^{(1)}=(u_{xx},u_{yy},u_{xy})^{T}$ where $T$ stands for transpose. For higher order terms with $m>1$, $\mathbf{u}^{(m)}$ contains $d(d+m-1)!/[(d-1)!m!]$ components, where $d$ is the spatial dimension and $!$ represents the factorial, e.g., in 2D, $\mathbf{u}^{(2)}=(\partial_{x}^{2}u_{x},\partial_{x}\partial_{y}u_{x},\partial_{y}^{2}u_{x},\partial_{x}^{2}u_{y},\partial_{x}\partial_{y}u_{y},\partial_{y}^{2}u_{y})$.

As will be proved in the next section, for any stress-free system, the elastic energy (within linear elasticity) must be a quadratic form of these $\mathbf{u}^{m}$’s,

[TABLE]

The matrix in Eq. (8), which will be dubbed the $\mathcal{K}$ matrix is a semi-positive-definite symmetric real matrix, composed of elastic constants.

Below we prove that a semi-positive definite $\mathcal{K}$ can be decomposed into the following form $\mathcal{K}=\mathcal{L}\mathcal{D}\mathcal{L}^{T}$

[TABLE]

This decomposition is a block version of the so-called LDLT decomposition, where the matrix $\mathcal{D}$ in the middle is a diagonal matrix. $\mathcal{L}$ is a block lower-triangular matrix and $\mathcal{L}^{T}$ is its transpose. After this decomposition, it is easy to verify that the elastic energy takes the form of Eq. (6).

I.3 Stress-free systems and the $\mathcal{K}$ matrix

In this section, we prove that the elastic energy of a stress-free system must take the form of Eq. (8) and the $\mathcal{K}$ matrix must be semi-positive definite. For the leading order theory (i.e., $u^{(1)}$ only), this is obvious. However, because higher-order derivatives of the displacement field are included here, this conclusion needs to be proved with this more generic setup. In particular, because higher-order gradient terms can be obtained from lower-order ones (by taking derivatives), globally, $u^{(m)}$ and $u^{(n)}$ (with $m\neq n$) are not independent, and thus a careful proof becomes necessary.

Here, we define stress-free based on stability against cutting. We cut a system into smaller pieces. If every piece can maintain its original shape after cutting, no matter how we cut it, the system is called stress-free. With this definition, it is easy to check that for a stress-free system, its elastic energy needs to be locally semi-positive definite, i.e., for any deformation, the elastic energy density at every single point needs to increase or stay the same. This is in contrast to a stressed system, whose elastic energy only needs to be globally semi-positive definite, i.e., for any deformation, the total elastic energy must increase or stay the same, but the elastic energy in some local region could decrease. Obviously, if we cut this local region out from the system, it cannot maintain its shape. Instead, it will deform to further reduce the elastic energy.

Early on, we mentioned that globally, $u^{(m)}$ and $u^{(n)}$ (with $m\neq n$) are related, because we can get higher-order derivatives terms from lower-order ones by taking derivates. Here, we prove that $u^{(m)}$’s are locally independent, i.e., in a local region, we can design the displacement field to make $u^{(m)}$’s into any values that we want. This can be achieved via the Taylor series. Without loss of generically, here we consider a local region around the origin. Here, we can define a displacement field in the form of a Taylor series

[TABLE]

Near the origin, coefficients in this Taylor series are just the values of $u^{(m)}$’s, and thus by selecting the values of these coefficients, we can achieve arbitrary $u^{(m)}$’s. It must be emphasized that this is a local construction and is only valid in a domain near the origin. As we move far away from the origin, this deformation may hit singularities, at which the series stops to converge.

The local construction is already enough for us, because as shown above for a stress-free system, the elastic energy density needs to be semi-positive definite for any local regions in the system. This immediately implies that for any points in the system, the elastic energy density cannot contain linear terms of $u^{(m)}$’s. If some linear $u^{(m)}$ terms exists, we can always find some $u^{(m)}$’s that make the energy density at this point negative, and as shown above, such a set of $u^{(m)}$’s can always be realized using the construction above. This means that if we cut this local region out, it will deform to lower its energy, and thus the system is stressed. This proves that for a stress-free system, the elastic energy density must start from quadratic order of $u^{(m)}$’s. If we ignore higher order terms (i.e., in the linear response regime), this is a quadratic form as shown in Eq. (8).

Following the same line of thinking, we can prove that the quadratic form must be semi-positive definite in a stress-free system. If the $\mathcal{K}$ matrix contains negative eigenvalues, then we can design the stress field as shown above to make $u^{(m)}$’s the eigenvector of this negative eigenvalue in this local region. This will make the local energy density negative, which is prohibited for stress free systems.

In summary, for a stress-free system, $\mathcal{K}$ is semi-positive definite. Thus, we can use the idea demonstrated in Sec. I.1 to absorb cross terms into squares.

I.4 Leading-order theory

At the leading order (i.e., $\mathbf{u}^{(1)}$ only), Eq. (8) recovers the standard elastic theory presented in the Voigt notation. In the Voigt notation, the stress tensor is written as a $d_{c}=d(d+1)$ vector, $\mathbf{u}^{(1)}$, where $d_{c}=d(d+1)/2$ is known as the Cartan dimension Cartan (1927). The elastic energy is a quadratic form of the stress tensor

[TABLE]

Here, $\mathcal{K}^{(1,1)}$ is the $d_{c}\times d_{c}$ real symmetric matrix, which contains $d_{c}(d_{c}+1)/2$ elastic moduli, among which $d(d-1)/2$ of them can be made zero by choosing proper coordinate orientation Landau and Lifshitz (1986). As shown in Eq. (11), via diagonalizing the $\mathcal{K}^{(1,1)}$ matrix, the elastic energy can be written as the sum of $d_{c}=d(d+1)$ terms, each of which is the square of certain linear combination of $\mathbf{u}^{(1)}$. Here, $\mathcal{O}$ is an orthogonal matrix, which diagonalizes $\mathcal{K}^{(1,1)}$ [$\mathcal{O}\mathcal{K}^{(1,1)}\mathcal{O}^{T}=\textrm{diag}(\lambda_{1},\lambda_{2},\ldots\lambda_{d_{c}})$ with $\lambda_{i}$ representing the eigenvalues of $\mathcal{K}^{(1,1)}$]. For stability reasons, $\mathcal{K}^{(1,1)}$ must be semi-positive definite and thus the eigenvalues $\lambda_{i}$ must be non-negative. Following Maxwell’s counting argument, for each positive eigenvalue $\lambda_{i}>0$, the corresponding square term enforces one constraint

[TABLE]

Here, we rewrite the strain tensor $u^{(1)}_{j}$ as the gradient of the displacement field with proper coefficients represented by $\mathcal{P}_{jk\ell}$. If a displacement field satisfies all the constraints [Eq (12)], it is a zero mode. For a rigid solid with positive-definite $\mathcal{K}^{(1,1)}$, all the eigenvalues are positive and thus we have $d_{c}=d(d+1)/2$ constraints. On the other hand, the number of degrees of freedom is $d$, because the displacement field has $d$ components. Thus, above 1D, a rigid body is always over-constraint ($d_{c}>d$). This conclusion is in agreement with the generalized counting argument for continuous media Sun and Mao (2019) as well as the embedding theorem of Janet and Cartan Janet (1926); Cartan (1927).

As shown in the main text, for systems at the verge of mechanical stability, some of the $\lambda_{i}$’s may vanish. If $d(d-1)/2$ of them vanish and the remaining $d$ eigenvalues are positive, the system is at the Maxwell point, with the same number of constraints and degrees of freedom, within this leading-order theory.

I.5 Second order terms

The second-order derivatives of the displacement field form a rank-3 tensor $\mathbf{u}^{(2)}=\{\partial_{i}\partial_{j}u_{k}\}$, and the first two indices in this tensor are symmetric because $\partial_{i}\partial_{j}u_{k}=\partial_{j}\partial_{i}u_{k}$. As a result, $\mathbf{u}^{(2)}$ has $d^{2}(d+1)/2=dd_{c}$ components.

By considering both first- and second-order derivative terms, the elastic energy now takes the following form

[TABLE]

where we used Einstein’s notation to sum over repeated indices and the middle term is the cross term between $\mathbf{u}^{(1)}$ and $\mathbf{u}^{(2)}$. As shown in the previous section, $\mathcal{K}^{(1,1)}$ is a $d_{c}\times d_{c}$ matrix. $\mathcal{K}^{(2,2)}$ and $\mathcal{K}^{(1,2)}$ are $(dd_{c})\times(dd_{c})$ and $d_{c}\times(dd_{c})$ matrices respectively.

Same as shown in the previous section, the $\mathcal{K}^{(1,1)}$ matrix can be diagonalized by an orthogonal matrix $\mathcal{O}$. With the help of this orthogonal matrix, the elastic energy can be written as

[TABLE]

In the second term, we used the fact that $\mathcal{O}$ is an orthogonal matrix and thus $\mathcal{O}^{T}_{ki}\mathcal{O}_{ij}=\delta_{kj}$.

As demonstrated in Sec. I.1, for $\lambda_{i}>0$, we can absorb $2(\mathcal{K}_{kl}^{(1,2)}\mathcal{O}^{T}_{ki}u_{l}^{(2)})(\mathcal{O}_{ij}u^{(1)}_{j})$ into $\lambda_{i}(\mathcal{O}_{ij}u^{(1)}_{j})^{2}$, i.e.,

[TABLE]

Here, we don’t sum over $\lambda_{i}$, but all other repeated indices are summed over. For $\lambda_{i}=0$, the corresponding cross term $2(\mathcal{K}_{kl}^{(1,2)}\mathcal{O}^{T}_{ki}u_{l}^{(2)})(\mathcal{O}_{ij}u^{(1)}_{j})$ must vanish, because this quadratic form is semi-positive definite. In summary, we can always rewrite the elastic energy as

[TABLE]

Up to the second-order derivative terms, this is precisely the structure in Eq. (6). We can repeat the same procedure as higher order derivative terms are included into the elastic energy, which leads to Eq. (6) and  (8).

I.6 States of self stress and the Maxwell-Calladine index theorem

In this section, we show that similar to discrete systems, the transpose of the $\mathcal{C}$ matrix in the continuum theory dictates states of self stress.

As shown in the main text, the elastic energy in the continuum theory is

[TABLE]

and the null space of $\mathcal{C}$ (i.e., every $\mathbf{u}$ that satisfies $\mathcal{C}\mathbf{u}=0$) gives zero modes. In the main text, only the first set of squares in Eq. (6) (with coefficients $\lambda_{i}$) was considered. But Eq. (17) can be easily generalized to the full elastic energy (i.e. the entire Eq. (6) including all terms with coefficients $\lambda_{i}$, $\widetilde{\lambda}_{i}$, $\widetilde{\widetilde{\lambda}}_{i}$, $\ldots$). This can be achieved by adding more rows to the $\mathcal{C}$ matrix, and each row of the $\mathcal{C}$ matrix corresponds to one square (with a nonzero coefficient) in Eq. (6), which enforces one constraint to the system. If all the constraints are satisfied, $\mathcal{C}\mathbf{u}=0$, the displacement field is a zero mode.

Here, we study the null space of $\mathcal{C}^{T}$, which is also known as the left null space of $\mathcal{C}$, and show that it gives the states of self stress. Assume that there exists a $\mathbf{w}$ such that

[TABLE]

With this $\mathbf{w}$, the system has a state of self stress. Under this self stress, the elastic energy is

[TABLE]

With $\mathcal{C}^{T}\mathbf{w}=0$, it is easy to check that the cross term vanishes (i.e., $\mathbf{u}^{T}\mathcal{C}^{T}\mathbf{w}=0$). Thus, up to a constant $\int\mathbf{dr}\;\mathbf{w}^{T}\Lambda^{-1}\mathbf{w}/2$, this elastic energy is identical to Eq. (17).

It is easy to check that this elastic energy describes a system under stress, i.e., the stress

[TABLE]

even at zero deformation. Here, because the elastic energy contains higher-order gradients of the displacement field, their corresponding high-order stress fields are introduced here $\sigma^{(m)}$. For $m=1$, $\sigma^{(m)}$ recovers the standard stress tensor. It is easy to verify that based on cutting stability (as defined above), a stress-free system must have $\sigma^{(m)}=0$ for all $m$. And thus, $\mathbf{w}$ indeed induces stress in an elastic system. More importantly, this system is under perfect force balance, i.e., the force

[TABLE]

Thus this is a state of self stress. Because every $\mathbf{w}$ defines a state of self stress. the number of states of self stress is determined by the dimension of the null space of $\mathcal{C}^{T}$.

Because the null spaces of $\mathcal{C}$ and $\mathcal{C}^{T}$ give zero modes and states of self-stress, the Maxwell-Calladine theorem still holds in the continuum. This theorem has been proved in discrete lattices/networks. Our continuum theorem now generalizes it to continuous media. Here we briefly outline the proof in the continuum. As mentioned in the main text and above, the number of rows of $\mathcal{C}$ is the number of constraints $N_{c}$, while the number of columns is the number of degrees of freedom $N_{d}$. Thus, $\mathcal{C}$ is a $N_{c}\times N_{d}$ matrix, while its transpose $\mathcal{C}^{T}$ is $N_{d}\times N_{c}$. For an arbitrary matrix, the dimension of the null space is the the number of columns minus the rank of the matrix. Thus, the null space of $\mathcal{C}$ has dimension $N_{d}-\operatorname{rank}{(\mathcal{C})}$, while the null space of $\mathcal{C}^{T}$ has $N_{c}-\operatorname{rank}{(\mathcal{C}^{T})}$. Because the dimension of the null space of $\mathcal{C}$ ($\mathcal{C}^{T}$) is the number of zero modes (states of self stress), this means that

[TABLE]

where $N_{0}$ and $N_{s}$ represent the numbers of zero modes and states of self stress respectively. Because the rank of a matrix always equals to the rank of its transpose, $\operatorname{rank}{(\mathcal{C})}=\operatorname{rank}{(\mathcal{C}^{T})}$, the two equations above imply that

[TABLE]

which is known as the Maxwell-Calladine index theorem. In sharp contrast to lattice systems, here $N_{d}$, $N_{s}$, $N_{0}$ and $N_{s}$ are numbers in the continuum theory and are independent of microscopic details. For example, in the continuum theory constructed above (based on displacement field gradients), $N_{d}$ is the number of spatial dimensions, regardless of lattice structures and/or how many nodes a lattice unit cell have. While in lattice systems, its number of degrees of freedom per unit cell is $d\times N_{n}$, where $N_{n}$ is the number of nodes in a cell. The lattice degrees of freedom is very sensitive to lattice structure, and its number can become very big, if each unit cell contains a high number of nodes (and the number of nodes per cell diverges in a continuum medium), while in contrast, the continuum theory $N_{d}$ is always finite ($d$) and its value is independent of microscopic details.

I.7 Low-energy optical modes

In this section, we consider an elastic system that contains soft optical phonon modes, i.e., optical phonon modes with a small gap. If the gap of such a optical phonon mode (at zero wavevector $q=0$) is zero, (e.g., in a regular kagome lattice with nearest springs Sun et al. (2012)), the energy of this optical mode is comparable to (and even lower than) the acoustic ones, and thus it must be included in the continuum theory to correctly describe the low-energy physics. Generically, such a soft optical phonon mode has a small but finite gap (e.g., twisted or deformed kagome lattices Sun et al. (2012); Kane and Lubensky (2014)). For these gapped soft optical modes, their contributions to low-energy physics becomes negligible in the long-wave-length (small $q$) and low-energy (small $\omega$) limit. And thus if we focus on small wavevector ($q$ close to [math]), the continuum theory do not need to include these optical modes. However, at intermediate wavevector, at which the energy of the soft optical mode becomes comparable with that of the acoustic ones, including these soft modes provides a more accurate description.

If we include $n_{O}$ such soft optical modes in the continuum theory, $\phi_{i}$ with $i=1$, $2$, $\ldots$, $n_{O}$, we just need to introduce some more components for each $\mathbf{u}^{(m)}$. For example, $\mathbf{u}^{(1)}$ is now a $[d(d+1)/2+n_{O}]$-component vector. The first $d(d+1)/2$ components are the strain tensor (same as before), and in addition it now contains $n_{O}$ more components, which are these soft optical modes $\phi_{i}$. We do the same for all $\mathbf{u}^{(m)}$, i.e., introducting the $(m-1)$th order derivatives of $\phi_{i}$ as new components of $\mathbf{u}^{(m)}$. Then, we can follow the same procedure shown above to reorganize the elastic energy and to define the $\mathcal{C}$ matrix.

II Low-energy theory

In this section, we present a generic framework to compute coefficients in Eq. (6) from first principles. This objective is achieved through integrating out high energy degrees of freedom (i.e., optical modes). In a continuum theory of elasticity, this is a standard procedure. However, the conventional approach usually focuses on the dynamic matrix Mao and Lubensky (2011), which doesn’t carry topological information about Maxwell systems. Here, we develop a more general approach, which can be applied to both the $\mathcal{C}$ matrix and the dynamic matrix. For the dynamic matrix, it gives the same results as the conventional approach. Remarkably, as far as zero modes are concerned, although our method dramatically reduces the number of degrees of freedom, it contains the same zero modes identical to the full theory.

II.1 Elastic theory in the real space

In this section we consider an elastic medium with a periodic structure. It must be emphasized that here we don’t assume the system is an ideal lattice model composed of ideal springs. Instead, the scope of the study also includes more realistic setups used in experiments, which requires finite element analysis to describe.

In an ideal discrete lattice model, the displacement field is defined on a set of discrete nodes. For a more-realistic continuous medium, we can still utilize a discrete set of nodes to describe the system through the finite-element method (FEM). For elastic properties, the FEM is a reliable approximation, which becomes asymptotically exact as the mesh becomes denser and denser. Below, we will use this discrete setup to define and compute the topological index. Because a topological index measures the topological property of a system, which is a qualitative property insensitive to quantitative details, as long as the mesh is not too sparse, the FEM is expected to produce the exact value of the topological index.

Here, we consider an infinite-largely periodic structure. In each unit cell, we introduce the same set of discrete nodes, which will be labeled by an integer $\alpha=1$, $2$, $\ldots$, $N_{n}$ with $N_{n}$ being the number of nodes per cell. We will label each unit cell by its real-space coordinate $\mathbf{R}$, which is a lattice vector $\mathbf{R}=\sum_{i=1}^{d}n_{i}\mathbf{a}_{i}$ with $n_{i}$ being integers and $\{\mathbf{a}_{i}\}$ representing a set of primitive lattice vectors. With this setup, the deformation of a node is represented by

[TABLE]

Here, $\mathbf{u}$ is the displacement vector of the $\alpha$th node in the unit cell located at $\mathbf{R}$. Because $\mathbf{u}$ is a $d$-dimensional vector, below we will use a subindex $i=1$, $\ldots$, $d$ to label its components.

Without loss of generality, here we assume the mass matrix is an identity matrix, which can always be achieved by changing the basis. As a result, the frequency of a mode is just the square root of the energy. With a periodic structure, the elastic energy of the system can be written as

[TABLE]

where $E_{\mathbf{R}}$ is the elastic energy associated with the unit cell located at $\mathbf{R}$. Due to the lattice translational symmetry, the elastic energy associated with two different unit cells (e.g. $E_{\mathbf{R}_{a}}$ and $E_{\mathbf{R}_{b}}$) are related by a lattice transition with a translation vector $\mathbf{t}=\mathbf{R}_{b}-\mathbf{R}_{a}$. It is also worthwhile to mention that $E_{\mathbf{R}}$ depends not only on nodes inside the unit cell $\mathbf{R}$, but also nodes from neighboring cells, because the elastic energy in general contains cross terms between nodes in different unit cells.

Within linear elasticity, each $E_{\mathbf{R}}$ is a quadratic form of the displacement field, and for a stress-free system, it is semi-positive definite. For a semi-positive-definite quadratic form, it can always be written as the sum of a set of squares.

[TABLE]

where $\kappa_{m}>0$ are positive coefficients and $\xi_{m}$ are linear combinations of the displacement field

[TABLE]

Here, $\mathcal{M}_{mi}^{(\alpha)}(\mathbf{t})$ are some real coefficients and $\mathbf{t}$ is a lattice vector. For nodes in the unit cell $\mathbf{R}$, $\mathbf{t}=0$, but because nodes from neighboring cells will also arise here, it is necessary to include nonzero $\mathbf{t}$ in $\xi_{m}(\mathbf{R})$.

Mathematically, the reason that we can always written $E_{\mathbf{R}}$ in this square sum form [Eq. (27)] is because each quadratic form can be represented by a real symmetric matrix, and once we diagonalize the matrix and use its eigenvectors as the new basis to rewrite the quadratic form, it turns into the sum of a set of squares. Physically, from the view point of Maxwell’s counting argument, each square term (with a positive $\kappa_{m}$) enforces one constraint. In an ideal model composed of balls connected by ideal springs, each square term corresponds to the elastic energy of one spring and $\kappa_{m}$ is the spring constant.

II.2 Wave vector space formula

Now, we go to the wave vector space via a Fourier transformation,

[TABLE]

Here $\mathbf{u}^{(\alpha)}(\mathbf{q})$ is the displacement of the node $\alpha$ in the wave vector space, while $\mathbf{u}^{(\alpha)}(\mathbf{R})$ is the real space displacement of the node $\alpha$ in the unit cell $\mathbf{R}$. In Eq. (29), we sum over all lattice vectors, and the integral in Eq. (30) is over one Brillouin zone with $\Omega_{B.Z.}$ beging the volume of a Brillouin zone.

Notice that in this Fourier transformation, the phase factor is chosen to be $\exp(-i\mathbf{q}\cdot\mathbf{R})$ with $\mathbf{R}$ being a lattice vector. This convention has two consequences. First of all, it ensures that all functions defined in the wave vector space are periodic with respect to the Brillouin zone, e.g., $\mathbf{u}^{(\alpha)}(\mathbf{q}+\mathbf{G})=\mathbf{u}^{(\alpha)}(\mathbf{q})$ where $\mathbf{G}$ is a reciprocal lattice vector. This periodicity in the wave vector space is crucial for the definition of the topological index, same as in lattice systems. Secondly, it is easy to notice that this phase factor of the Fourier transformation depends on how we choose the unit cell. For a periodic structure, there exist infinitly many equivalent ways to define a unit cell Kittel et al. (1976). Each option is known as a gauge choice. Such a gauge choice will not change any physical observables, i.e., physics observable should be gauge invariant. However, if we choose a different unit cell, the same node may now be assigned to a different cell (with a different $\mathbf{R}$) and thus the corresponding phase factor $\exp(-i\mathbf{q}\cdot\mathbf{R})$ will be shifted. This phase factor turns out to play a very important role in the definition of the topological index, and most importantly, for Maxwell lattices, the conventional way of computing the topological index is sensitive to the gauge choice, and different gauge choices will shift the index by an integer constant. As a result, if we want to connect the topological index with a physics observable, e.g., number of edge modes, it is necessary to deduct the gauge sensitive part.

In the wave vector space, the elastic energy becomes

[TABLE]

where

[TABLE]

The coefficients $\kappa_{m}$ are defined in Eq. (27), and $\xi_{m}(\mathbf{q})$ is the wave vector space version of Eq. (28)

[TABLE]

This elastic energy can be written in a matrix form. Here, we write all displacement field (in the wave vector space) as a column vector $\mathbf{u}(\mathbf{q})=\{u^{(\alpha)}_{i}(\mathbf{q})\}$, and define a matrix $\mathcal{C}(\mathbf{q})$ as

[TABLE]

where the row index $m$ runs over all the square terms [defined in Eq. (27)], and the column index $(\alpha)i$ goes over all the components of the displacement vectors for all nodes. Notice that here we absorbed $\kappa_{m}$ into the definition of $\mathcal{C}(\mathbf{q})$, which differs slightly from the convention used in Maxwell lattices. There, a separate matrix is utilized for $\kappa_{m}$, which corresponds to the spring constant and is usually set to unity.

With the $\mathcal{C}$ matrix defined above, the elastic energy can be written as

[TABLE]

From this formula, it is easy to realize that the dynamic matrix $\mathcal{D}({\mathbf{q}})$ is

[TABLE]

For an idea lattice, this $\mathcal{C}(\mathbf{q})$ recovers the $\mathcal{C}$ matrix defined for discrete lattices, and the matrix becomes a square matrix for a Maxwell lattice. For a continuous medium, in general $\mathcal{C}(\mathbf{q})$ is not a square matrix. Instead, we expect the matrix to have more rows than columns, which implies that the system is over-constrained at the microscopic level.

It is worthwhile to emphasize here that same as in ideal lattice models, the $\mathcal{C}(\mathbf{q})$ matrix contains more information than the dynamic matrix. As shown in Eq. (36), $\mathcal{C}(\mathbf{q})$ uniquely determines the dynamic matrix. However, from the dynamic matrix, although it can be mathematically decomposed into the form of Eq. (36), the decomposition is not unique and thus $\mathcal{C}(\mathbf{q})$ cannot be uniquely identified from $\mathcal{D}({\mathbf{q}})$. Below, we will keep working with the $\mathcal{C}$ matrix, without condense its information into $\mathcal{D}$, which is crucial for the topological index.

II.3 Separating high- and low- energy modes.

For simplicity, we now drop the wave vector $(\mathbf{q})$ used in Eqs. (35) and (36). Here, we split displacement fields $\mathbf{u}({\mathbf{q}})$ into the low-energy and high-energy parts. This is achieved via changing of basis. We diagonalize the dynamical matrix at $\mathbf{q}=0$ and using its eigenvectors as the new basis. We emphasize here that the eigenvectors utilized here are obtained at $\mathbf{q}=0$, and thus this changing of basis is $q$-independent. Here, we choose to order the eigenvectors according to their eigenvalues from low to high. The first $d$ eigenvectors here are the acoustic modes with eigenvalues [math]. The rest are optical modes, which are usually gapped. Because the continuum theory focuses on small wave vectors ($\mathbf{q}$ close to [math]), this new basis naturally separates low- and high-energy modes

[TABLE]

Here, the first a few components will be our low-energy modes. Typically, we keep the lowest $d$ modes, i.e., acoustic modes. However, in some systems, especially systems near the verge of mechanical stability, the energy of certain optical modes may be close to zero. These soft optical modes can also be included in $\mathbf{u}_{L}$ as well. As shown in the main text, including soft optical modes can help improve the accuracy of the continuum theory, especially near the wave vector at which the energy of the soft optical mode becomes comparable with acoustic modes.

Once we split the displacement fields into the low- and high-energy sectors, the $\mathcal{C}$ matrix define in Eq. (34) is also splitted into two blocks

[TABLE]

The matrices $\mathcal{A}$ and $\mathcal{B}$ have the same number of rows as $\mathcal{C}$ (i.e. same number of constraints), while the number of columns of $\mathcal{A}$ ($\mathcal{B}$) is the dimensions of $\mathbf{u}_{L}$ ($\mathbf{u}_{H}$). With this block form, the elastic energy [Eq. (35)] becomes

[TABLE]

In Eq. (39), the low- and high- energy modes are mixed together through the cross terms in the elastic energy. These cross terms can be absorbed into the high energy mode through defining a new basis, and the elastic energy [Eq. (39)] becomes

[TABLE]

where

[TABLE]

and

[TABLE]

Here, we used the fact $(\mathcal{B}^{\dagger}\mathcal{B})$ is invertible. As shown in Eq (39), in the dynamical matrix, $(\mathcal{B}^{\dagger}\mathcal{B})$ is part of the dynamic matrix associated with the high-energy modes. Because these modes are gapped optical phonons, $(\mathcal{B}^{\dagger}\mathcal{B})$ must be positive definite and thus it is invertible.

This new basis decouples the low- and high-energy modes. Thus, for low-energy physics, we can now focus on the low-energy sector and set

[TABLE]

In comparison with Eq. (35), here the $\mathcal{C}$ matrix is reduced into its low-energy sector $\widetilde{\mathcal{A}}$, while the degrees of freedom is reduced from all displacement fields $\mathbf{u}$ into low-energy deformations only $\mathbf{u}_{L}$. The physical meaning of this reduction in degrees of freedom can be viewed as integrating out high-energy modes, which is a standard procedure in elastic theory, especially in the continuum theory. However, the conventional approach usually reduces the dynamic matrix into its low-energy sector. Here, we generalized the approach to the $\mathcal{C}$ matrix, making it possible to study topological properties via low-energy degrees of freedom.

Another way to understand this reduction is by trying to find low-energy deformations via minimizing the elastic energy. From Eq (40), it is easy to notice that this minimization can be achieved by setting $\widetilde{\mathbf{u}}_{H}$ to zero, and then search for low-energy excitations in the reduced theory [Eq. (43)].

To conclude this section, we emphasize that this $\mathcal{\widetilde{A}}$ is not a square matrix, in direct contrast to the $\mathcal{C}$ matrix utilized in discrete Maxwell lattices. This is because our reduction of degrees of freedom only removes high-energy degrees of freedom, while all the constraints are maintained in the reduced theory. As a result, $\mathcal{\widetilde{A}}$ and $\mathcal{C}$ matrices have the same number of rows (i.e., the same number of constraints), but $\mathcal{\widetilde{A}}$ contains fewer number of columns (i.e., fewer degrees of freedom).

II.4 Exactly solution for zero modes

In this section, we show that there exists a one-to-one correspondence between zero modes in the reduced theory, i.e., $\mathcal{\widetilde{A}}\mathbf{u}_{L}=0$ and zero modes in the full theory $\mathcal{C}\mathbf{u}=0$. Therefore, as far as zero modes are concerned, our reduced theory offers the exact solutions with no error or approximation.

If $\mathbf{u}=(\mathbf{u}_{L},\mathbf{u}_{H})$ is a zero mode satisfying $\mathcal{C}\mathbf{u}=0$, we have

[TABLE]

By multiplying to both sides $(\mathcal{B}^{\dagger}\mathcal{B})^{-1}\mathcal{B}^{\dagger}$, we get $\mathbf{u}_{H}=-(\mathcal{B}^{\dagger}\mathcal{B})^{-1}\mathcal{B}^{\dagger}\mathcal{A}\mathbf{u}_{L}$. If we substituted this $\mathbf{u}_{H}$ into Eq. (44), we get

[TABLE]

Thus, $\mathcal{C}\mathbf{u}=0$ must implies $\mathcal{\widetilde{A}}\mathbf{u}_{L}=0$.

We now start from a zero mode in the reduced theory with $\mathcal{\widetilde{A}}\mathbf{u}_{L}=0$. By setting $\mathbf{u}_{H}=-(\mathcal{B}^{\dagger}\mathcal{B})^{-1}\mathcal{B}^{\dagger}\mathcal{A}\mathbf{u}_{L}$, it is easy to check that

[TABLE]

II.5 The continuum theory

As mentioned in the previous section, the $\mathcal{\widetilde{A}}$ matrix is not yet the $\mathcal{C}$ matrix of the continuum theory. One more procedure is still needed to get the coefficients in Eq. (6).

Here, we take a Taylor expansion and write $\mathcal{\widetilde{A}}$ as a power series of $\mathbf{q}$. For the continuum theory, although keeping more terms improves its accuracy at larger $\mathbf{q}$, two leading-order terms are usually sufficient as far as the topology at long-wavelength is concerned. For acoustic modes, the leading order theory is linear in $\mathbf{q}$, and we also need to keep the next order terms ($\propto\mathbf{q}^{2}$) as well. For soft-optical modes (if they are included in the low-energy sector), their leading order terms are constants independent of $\mathbf{q}$, and the next order terms, which are proportional to $\mathbf{q}$, are also needed.

To obtain the continuum theory, now we go back to the real space (replacing $\mathbf{q}$ with $-i\nabla$). Here, $\mathcal{\widetilde{A}}\mathbf{u}_{L}$ becomes a linear combination of displacement fields and their derivatives. And in real space, the elastic energy is

[TABLE]

where the arrow above $\mathcal{\widetilde{A}}$ implies that the derivatives in $\mathcal{\widetilde{A}}$ are for $\mathbf{u}_{L}$ on the right side of it, while the derivatives in the transpose matrix $\mathcal{\widetilde{A}}^{T}$ are for $\mathbf{u}_{L}^{T}$ on the left side.

Now, Eq. (47) recovers the generic form of the elastic energy [Eq. (8)], and we can follow the procedure discussed in that section to reorganize the elastic energy into the form shown in Eq. (6). As shown in the main text, from there, the $\mathcal{C}$ matrix and topological indices in the continuum theory can be defined.

III Topological index, zero modes and gauge choice

In this section we discuss the connection between the topological index and zero modes, as well as how to make the index gauge independent.

Discussions in Sec. III.1 applies to both Maxwell lattices and Maxwell media, while Secs. III.2, III.3 and III.4 focus on discrete systems. And Sec. III.5 is devoted to the continuum theory.

III.1 Real-space zero modes and zeros in $\det\mathcal{C}(\mathbf{q})$

In this section, we use a discrete setup to demonstrate the physics, but the same conclusions can be easily generalized to the continuum theory using the continuum elastic energy and its $\mathcal{C}$ matrix defined in the main text. Here, we show that zero modes in the real space is related with zeros in the determinant of the wave vector space $\mathcal{C}$ matrix, $\det\mathcal{C}(\mathbf{q})=0$.

As shown above, for an elastic medium with a periodic structure, the total elastic energy can be written as the sum of $E_{\mathbf{R}}$, which is the elastic energy associated with each unit cell as shown in Eq. (27).

Here, we consider a deformation pattern,

[TABLE]

where $q_{x}$ is real and $q_{y}$ is complex. This setup describes bulk waves (if $q_{y}$ is real) and edge modes along edges parallel to the $x$ axis. From Eq. (27), the elastic energy for this deformation can be easily obtained

[TABLE]

where $q_{y}^{\prime}$ and $q_{y}^{\prime\prime}$ are the real and imaginary part of $q_{y}$, $q_{y}=q_{y}^{\prime}+iq^{\prime\prime}_{y}$. Although this formula may look like Eq. (35), its physical meaning is different. Here, we are computing the elastic energy of a deformation in the real space, while Eq. (35) is a wave vector space formula. From this $\langle E_{R}\rangle$, it is easy to realize that a zero energy mode must satisfy

[TABLE]

For a Maxwell system, $\mathcal{C}(\bar{q}_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$ is a square matrix, and thus to have a nontrivial solution, we have to require

[TABLE]

III.2 $\det\mathcal{C}$ in discrete lattices

In this section, we focus on discrete Maxwell lattices and show that at a fixed $q_{x}$, as long as interactions are short-ranged, $\det\mathcal{C}(\bar{q}_{x},q_{y})$ must take the following functional form

[TABLE]

where $z=\exp(-i2\pi q_{y}/b_{2})$ is a complex variable with $b_{2}$ being the length of a reciprocal lattice vector (see below for details). $c$ and $\{z_{i}\neq 0\}$ are complex constants, while $N$ and $M$ are integers. Here, $z_{i}$’s are nonzero and $M\geq 0$, while $N$ can be positive, zero or negative. As will shown below, this implies that as a function of $z$, the function utilized to define the topological index, $\operatorname{tr}\mathcal{C}^{-1}\partial_{z}\mathcal{C}$, is an meromorphic function for the entire complex $z$ plane (i.e., analytic except for some discrete poles). In addition, the number of poles is finite, which equals to $M+2$. Among these $M+2$ poles, $M$ of them ($z=z_{i}$) correspond to zero modes and their contributions to the topological index are gauge invariant. The other $2$ poles ($z=0$ and $z=\infty$) don’t corresponds to physical zero modes, and their contributions make the topological index gauge dependent.

Here, we use $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ to represent the two primitive lattice vectors, and $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ as their corresponding reciprocal lattice vectors. For a Maxwell lattice, typically, the edge is chosen to be along a lattice vector, and here we set the edge to be along $\mathbf{a}_{1}$, which is our $x$ axis. Because the $y$ axis is perpendicular to the edge ($\mathbf{a}_{1}$), by definition, $\mathbf{b}_{2}$ must be along the $y$ axis.

The matrix elements of the $\mathcal{C}$ matrix are shown in Eq. (34), and its $q_{y}$ dependence comes from the exponential $\exp(-i\mathbf{q}\cdot\mathbf{t})$, where $\mathbf{t}$ is a lattice vector $\mathbf{R}=m\mathbf{a}_{1}+n\mathbf{a}_{2}$. Thus

[TABLE]

Here we use the fact that $\mathbf{q}=q_{x}\mathbf{\hat{x}}+q_{y}\mathbf{\hat{y}}$ with $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ being unit vectors along $x$ and $y$. We now focus on the second phase factor in Eq. (53), which is associated with $q_{y}$. Because $\mathbf{\hat{y}}\parallel\mathbf{b}_{2}$, we have $\mathbf{\hat{y}}\cdot a_{1}=0$, and thus

[TABLE]

Here, we define $z=\exp(-i2\pi q_{y}/b_{2})$ with $b_{2}$ being the length of $\mathbf{b}_{2}$. And by definition, $n$ is an integer.

This proves that as a function of $z$, each element of the $\mathcal{C}$ matrix is a Laurent polynomial (i.e. $\sum_{n}c_{n}z^{n}$ with the power $n$ being an integer, which is not necessarily positive). More importantly, this Laurent polynomial must have a finite length, i.e. the sum contains only finite number of terms, instead of an infinitely long power series.

[TABLE]

with $c_{n}$ be some complex coefficients and the integer power $n$ must stay in a finite window $M_{\textrm{min}}\leq n\leq M_{\textrm{max}}$. Here $M_{\textrm{min}}$ and $M_{\textrm{min}}$ are finite integers, and they could be positive, negative or zero. This is because in the sum of Eq. (34), each $\mathbf{t}$ implies a cross term in the elastic energy between two nodes whose unit cells are separated by distance $\mathbf{t}$. For a local theory with short-distance interactions, this $\mathbf{t}$ cannot be large, and thus it sets a cap to the absolute value of the integer power $n$.

Because the determinant is the sum of products of matrix elements, if each matrix element is a finite Laurent polynomial [Eq. (55)], the determinant must also be a finite Laurent polynomial

[TABLE]

The $c_{n}$ here is a set of coefficients different from the ones in Eq. (55). And the integer power $n$ is between $N_{\textrm{min}}$ and $N_{\textrm{max}}$. Here, we used the fact that the $\mathcal{C}$ has a finite dimension, which is true for discrete Maxwell lattices where the dimension of the matrix is $N_{d}\times N_{d}$ with $N_{d}$ being the number of the degrees of freedom per cell.

This $\det{\mathcal{C}}(z)$ can be rewritten as

[TABLE]

Here $\{z_{i}\}$ is a set of nonzero complex constants, at which $\det{\mathcal{C}}(z_{i})=0$. $c=c_{N_{\textrm{max}}}$ is a complex constant, and the two integers $N=N_{\textrm{min}}$ and $M=N_{\textrm{max}}-N_{\textrm{min}}$ ($M\geq 0$ but $N$ cant be positive, negative or zero).

As shown in Sec. III.1, each $z_{i}$ corresponds to one zero modes with $\det{\mathcal{C}}(z_{i})=0$. It is worthwhile to mention that in addition to these nonzero solutions ($z_{i}$), $\det{\mathcal{C}}(z)=0$ may have one additional solution at $z=0$ (if $N>0$) or $z=\infty$ (if $M+N<0$). As will be shown below, they don’t correspond to physical zero modes, and their contributions makes the topological polarization gauge dependent..

III.3 Topological index and gauge choice

As mentioned in Sec. II.2, by choosing a different gauge, i.e., a different unit cell convention, the elastic energy associated with a unit cell $E_{R}$ [Eqs. (26) and (27)] will change accordingly, which changes the $\mathcal{C}$ matrix [Eq. (34)]. This change of $\mathcal{C}$ can be represented by multiplying a diagonal matrix $\mathcal{D}_{1}$ from the left side of $\mathcal{C}$, i.e., $\mathcal{C}\to\mathcal{D}_{1}\mathcal{C}$, where the diagonal elements of this diagonal matrix $\mathcal{D}_{1}$ are $z$ to some integer powers. In addition, because some nodes will be assigned to a different cell as we change the unit cell convention, the effect also changes the associated phase factor in $\mathcal{C}$, which corresponds to multiplying an diagonal matrix $\mathcal{D}_{2}$ from the right side, i.e. $\mathcal{C}\to\mathcal{C}\mathcal{D}_{2}$. Same as $\mathcal{D}_{1}$, the diagonal elements of $\mathcal{D}_{2}$ are also $z$ to some integer powers.

As a result, changing unit cell convention changes $\det\mathcal{C}$ into $\det\mathcal{D}_{1}\det\mathcal{C}\det\mathcal{D}_{2}$, i.e

[TABLE]

where $\Delta N$ is some integer. This immediately implies that in Eq. (57), the integer power $N$ is gauge dependent, where the product part $\Pi_{i=1}^{M}(z-z_{i})$ is gauge invariant. Depending on how we choose a unit cell, the value of $N$ can be shifted by an arbitrary integer, but the product part always remain. As will be shown below, the product part is why the topological index can detect edge modes, while the $z^{N}$ part is the reason why the topological polarization is gauge dependent.

In a Maxwell lattice, the topological index, i.e., each component of topological polarization, is defined as

[TABLE]

where the contour is the unit circle of the complex $z$ plane. Using Eq. (57), it is easy to verity that

[TABLE]

and thus

[TABLE]

Here, we emphasize again that $z_{i}\neq 0$. As discussed above, the first integral is gauge invariant, but the second integral is gauge dependent, because the value of the integer $N$ changes if a different unit-cell definition is adopted.

In Maxwell lattices, the integral contour is usually set as the unit circle on the complex $z$ plane, and this integral can be easily computed utilizing the Cauchy integral theorem.

[TABLE]

where $N_{B}$ is the number of $z_{i}$’s inside the unit circle. Because each $z_{i}$ corresponds to a zero mode, as shown above, $N_{B}$ counts the number of zero modes whose $z_{i}$’s are located inside the unit circle (aside from at the origin). By definition, $z_{i}$ inside the unit circle means that the imaginary part of $q_{y}$ is positive, and thus according to Eq. (48), such a zero mode is localized at the bottom edge of a strip. In other words, $N_{B}$ gives the number of zero modes for the bottom edge. This is the key reason why $n_{i}$ is sensitive to the number of edge modes for each edge.

However, in addition to $N_{B}$, $n_{i}$ also contains a gauge-dependent contribution $N$, whose value is sensitive to the unit-cell convention. This gauge-dependent $N$ is the reason why the value of $n_{i}$ changes if a different unit cell is adopted. As shown in the loop integral above, this gauge-dependent part comes from the pole of $\partial_{z}\ln\det\mathcal{C}$ at $z=0$. Thus, to obtain an gauge independent topological index, this pole should be removed from the integral, which can be achieved utilizing the contour shown in Fig. 3(c) in the main text (i.e., the boundary of the region B).

If we reverse the direction of the unit-circle contour and then applies the Cauchy integral theorem, then $-n_{i}$ should equals to the total residue summed over all poles outside the unit circle. For each $z_{i}$ with $|z_{i}|>1$, it gives a pole with residue $1$, and thus contributes $1$ to the integral. Because $|z_{i}|>1$ implies $q^{\prime\prime}_{y}<0$, these are edge modes on the top edge, i.e., each edge mode on the top edge gives $1$ to the integral. In addition, $\partial_{z}\ln\det\mathcal{C}(z)$ also have a pole at $z=\infty$ and its residue is $-(N+M)$. This is because at $z=\infty$, $1/z$ gives a residue $-1$, and in addition each $1/(z-z_{i})$ also has a residue $-1$ at $z=\infty$. Thus, from Eq. (61), the total residue at $z=\infty$ is $N-M$. After adding all the poles ($z=z_{i}$ with $|z_{i}|>1$ and $z=\infty$), we get

[TABLE]

where the first term comes from zero-energy modes at the top edge, and the second term is from the pole at $z=\infty$. Again, the first term here is gauge invariant, but the second term is gauge dependent. It is worthwhile to mention that this result agrees with $n_{i}=N_{B}+N$ that we obtained above. This can be easily checked by noticing that $M=N_{B}+N_{T}$. Here, we can avoid the gauge-dependent contribution from the pole at $z=\infty$ by using the contour shown in Fig. 3(c) in the main text (i.e. the boundary of the region T). Once $z=\infty$ is removed, the gauge-dependent second term ($N-M$) disappear and the topological index directly gives $N_{T}$.

In kagome lattices, one frequently used gauge has $N=-1$. In the non-polarized phase, where each edge has one zero mode per unit cell ($N_{B}=N_{T}=1$), this leads to $n_{i}=N_{B}+N=0$. For the polarized phase, both two zero modes are located on the bottom edge ($N_{B}=2$ and $N_{T}=0$), and thus $n_{i}=N_{B}+N=1$. By keeping $z=0$ and $z=\infty$ out of the integral contour, the topological index directly gives $N_{B}$ and $N_{T}$ without the gauge dependent part $N$.

In summary, in this section we showed that $\operatorname{tr}\mathcal{C}^{-1}\partial_{z}\mathcal{C}$ is analytic for the entire complex plane, except for $M+2$ poles. Among these poles, $M$ of them comes from zero modes, each of which has residue $1$ and is gauge invariant. In addition to physical zero modes, there are two additional poles located at $z=0$ and $z=\infty$. As will be shown in the next section, these two poles are not due to zero modes. Instead, they are required by analytic continuation (as we extend the definition of $\det\mathcal{C}$ from real $q_{y}$ to the entire complex $z$ plane), and their contribution to the topological index is gauge dependent. By choosing a new set of contours, their contributions can be removed, which leads to a gauge independent topological index.

Finally, it is worthwhile to mention that in the calculation above, we assumed that $\det\mathcal{C}(z)\neq 0$ on the unit circle ($|z|=1$). Because the unit circle corresponds to real $q_{y}$, this means that the bulk phonon modes (at this fixed nonzero $q_{x}$) are fully gapped without bulk zero mode. If this is not the case, the topological index becomes ill-defined and the system (at this $q_{x}$) is at the phase boundary between two different topological phases (i.e., a topological phase transition Rocklin et al. (2017) or mechanical Weyl points Rocklin et al. (2016)).

III.4 Dangling bonds and analytic continuation

In this section, we discuss two subjects (a) edge modes from dangling bonds and (b) why analytic continuation necessarily leads to poles that are not related with edge modes.

Here, we use the terminology “dangling bonds” broadly to refer to all zero-energy edge modes with zero penetration depth, i.e., the imaginary part of $q_{y}$ is infinity. Due to its local nature (i.e. zero penetration depth), these zero modes are associated local degrees of freedom at the edge, instead of a bulk origin. Thus, bulk quantities like bulk topological indices are not expected to capture this physics. In terms of $z$, these “dangling bonds” corresponds to $z=0$ (for $q^{\prime\prime}_{y}=+\infty$) and $z=\infty$ (for $q^{\prime\prime}_{y}=-\infty$). Thus, because bulk physics is not expected to capture zero modes at $z=0$ and $z=\infty$, contributions from poles at these to points are not expected to carry correct information about the edge, and this is the reason why they are gauge dependent. Furthermore, when a Maxwell lattice is cut at different depths at an edge (e.g., shifting the cut up or down for a horizontal edge), the number of edge zero modes that have a topological origin (i.e., due to bulk polarization) is invariant, and only the number of zero modes due to dangling bonds changes.

Here, we demonstrate that due to analytic continuation, the poles of $\operatorname{tr}\left(\mathcal{C}^{-1}\partial_{z}\mathcal{C}\right)$ cannot all correspond to physical edge modes. As mentioned above, each zero modes corresponds to a pole with residue $1$ in $\operatorname{tr}\left(\mathcal{C}^{-1}\partial_{z}\mathcal{C}\right)$, which is the reason why the topological index is sensitive to zero mode. However, analytic continuation requires this function to have additional poles, beyond those who come from zero modes. This is because according to the Cauchy integral theorem, the total residue from all poles in the entire complex plane must be zero (we showed above that $\mathcal{C}$ has no other singularity beyond poles for the entire complex plane). Thus, there must be some poles with negative residue to counter balance the poles from edge modes, which is the only way to keep total residue zero. Where can these additional poles be located? As shown above, $z=0$ and $z=\infty$ corresponds to local degrees of freedom at the boundary, and thus for a bulk topological index, they are the perfect place for these extra poles.

III.5 Continuum theory

In the continuum theory, same as shown in Sec. III.1, zero modes are also associated with zeros in $\det\mathcal{C}$. Thus, in terms of the complex variable $z$ [i.e., Fig 2.(b) in the main text], again we can search for edge modes by finding zeros of $\det\mathcal{C}(z)$.

Different from lattice models, here, the continuum theory is only expected to hold at long distance (small $q$), thus in the complex $z$ plane, the theory only applies to a domain near $z=1$ (i.e., $q_{y}=0$). In the vicinity of $z=1$, $\det\mathcal{C}(z)$ should be an analytic function and thus can be written as a power law expansion,

[TABLE]

In contrast to lattices, this sum contains no negative power. In addition, this power-law expansion is only valid in a domain around $z=1$. Outside this domain, this function in general is not meromorphic for the entire complex plane, and essential singularities that cannot be represented as poles or branch cuts may arise as we move far away from $z=1$. This is the mathematical reason why a cut off needs to be introduced, which keeps us inside this domain.

Because we always stay inside the convergence domain of the power law expansion, similar to the lattice models discussed above, one can easily prove that here, if $\det\mathcal{C}(z_{i})=0$ and $z_{i}$ is inside the convergence domain, $\operatorname{tr}\left(\mathcal{C}^{-1}\partial_{z}\mathcal{C}\right)$ shall have a pole at $z=z_{i}$ with residue $1$. Thus, the same integral shown above counts the number of zero modes inside the integral contour.

Because our power expansion here has no negative power and because our integral contour is alway inside the convergence domain, there is no unphysical pole, i.e., poles inside our contour has a one-to-one correspondence with physical zero modes and their contribution to the integral is gauge invariant. This is why this topological index is gauge independent and directly measures the number of zero modes.

Finally, we emphasize that the integral contour utilized in the continuum theory only covers part of the unit circle. Thus, for a Maxwell lattice, it is possible that there might exist zero modes whose $z_{i}$’s are outside the integral contour of the continuum theory. For those zero modes, because their $z_{i}$’s is far away from $1$, i.e. $q_{y}$ far away from [math], they involve short-distance physics, which the continuum theory is not expected to capture.

IV Beyond Maxwell lattices

In this section, we present the variational method discussed in the main text. For systems not far from the Maxwell point, we can first treat the system as a Maxwell medium by ignoring terms associated with small eigenvalues, as shown in the main text. Within this zeroth order approximation, the continuum theory is at the Maxwell point, and we can compute the topological index and obtain the decay rate of the edge modes. Afterwards, we can use this zeroth order results as input and compute the energy expectation value of this mode.

Similar to the ground state of a quantum system, at a fixed wave vector, the deformation configuration that minimizes the elastic energy is the normal mode with the lowest eigen-energy. Instead of performing a full minimization, here we use the idea of the variational principle. Instead of searching all possible real space configurations, we limit our deformation to the following functional form

[TABLE]

Here we fix one component of the wave vector along the edge direction ($q_{x}=\bar{q}_{x}$), and the amplitude $\mathbf{A}$ as well as $q^{\prime}_{y}$ and $q^{\prime\prime}_{y}$ are variational parameters, which we will adjust in search for the lowest-energy mode. Such a minimization will produce an estimation of the normal mode with lowest eigen energy at $q_{x}=\bar{q}_{x}$. If $q^{\prime\prime}_{y}=0$, it is a bulk mode. If $q^{\prime\prime}_{y}\neq 0$, it is an edge mode. In a rigid solid, the true lowest-energy mode (e.g. Rayleigh waves) don’t take this form. Instead, they are often composed of linear combination of multiple exponentials, each of which contains a different set of $q^{\prime}_{y}$ and $q^{\prime\prime}_{y}$. However, in our leading order theory at the Maxwell point, the mode with lowest energy indeed take this form, and thus as long as a system is not far from the ideal Maxwell limit, this waveform provides a good approximation.

It should be emphasized that such a constrained minimization is not expected to produce the exact answer. Instead, it always over estimate the energy for the lowest energy mode. However, if the initial guess we made are reasonably close to the real minimum, this approach gives a good upper bound for energy. In our calculation here, we will further simplify the problem by setting $q^{\prime}_{y}$ and $q^{\prime\prime}_{y}$ to be what we obtained from the zeroth order theory. As mentioned early on, at the zeroth order, we can use the topological theory to obtain edge modes, as well as their $q^{\prime}_{y}$ and $q^{\prime\prime}_{y}$, which are the values that we will use here. The idea here is that we are using the profile revealed by the zeroth-order topological theory to guide our variational calculation. As long as the system is not far from the ideal Maxwell point, this zeroth order theory should produce reliable estimation about the edge modes.

By plugin this zeroth order profile into the elastic energy Eq. (27), it is easy to realize that the energy expectation value is

[TABLE]

where the $\mathcal{C}$ matrix is defined in Eq. (34), and the amplitude for each node $\mathbf{A}$ is our variational parameter. Because the matrix $[\mathcal{C}(\bar{q}_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})]^{\dagger}\mathcal{C}(\bar{q}_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$ is unitary, minimizing $\langle E_{R}\rangle$ here corresponds to finding the lowest eigenmode of this matrix. Thus, within this variational method, the energy of the lowest-energy mode in our system is the lowest eigenvalue of $[\mathcal{C}(\bar{q}_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})]^{\dagger}\mathcal{C}(\bar{q}_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$, and its square root gives a good estimation about the frequency of the mode, which is what we plotted in Fig. 3. Because the zeroth order theory gives two edge modes for each $q_{x}$ (with different $q^{\prime}_{y}$ and $q^{\prime\prime}_{y}$ in general), the procedure discussed above are performed for each of them, which gives the two edge-mode branches in Fig. 3. As can be seen in Fig. 3, in comparison with full scale numerical simulations, these two curves slightly overestimate the mode frequency, which is expected due to the variational nature of our method.

It must be emphasized that here we are searching for the lowest eigenvalue of $[\mathcal{C}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})]^{\dagger}\mathcal{C}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$. This should not be confused with the lowest eigenvalue of the dynamic matrix $\mathbf{D}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$. For bulk modes (i.e., real $q_{x}$ and $q_{y}$), $\mathbf{D}(q_{x},q_{y})=[\mathbf{C}(q_{x},q_{y})]^{\dagger}\mathbf{C}(q_{x},q_{y})$ as have been shown above, but it is not the case for edge modes with complex $q_{y}$. This is because in reality, the dynamic matrix is defined as

[TABLE]

instead of $\mathcal{C}^{\dagger}\mathcal{C}$. For real $q_{x}$ and $q_{y}$, $[\mathbf{C}(-q_{x},-q_{y})]^{T}$ is identical to $[\mathbf{C}(q_{x},q_{y})]^{\dagger}$. However, because Eq. (67) is what is really utilized in the definition of the dynamic matrix, as we generalize $q_{y}$ into complex values (for the edge modes), the dynamic matrix becomes

[TABLE]

which differs sharply from the matrix $[\mathcal{C}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})]^{\dagger}\mathcal{C}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$ that we used above, Most importantly, in contrast to our hermitian matrix $[\mathcal{C}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})]^{\dagger}\mathcal{C}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$, the matrix $\mathbf{D}(q_{x},q^{\prime}_{y}+iq^{\prime\prime}_{y})$ is not hermitian with a complex $q_{y}$, and thus its eigenvalues have no direct physical meanings. This is the key reason why edge-mode information is not obtained from the dynamic matrix.

V Lattice models

In this section, we present details about the lattice modes that we used for Figs 2 and 3 in the main text and the coefficients of their continuum theory.

The lattice we use is shown in Fig. 4, which is a topological kagome lattice with point-like sites connected by central-force bonds (balls connected by ideal springs). Sides of red triangles are of lengths $(1,1,1)$ and sides of the blue triangles of at lengths $(1,1/\sqrt{3},1/\sqrt{3})$. The angle $\theta$ between the two triangles is $\theta=5\pi/12$. Each unit cell contains 3 sites (e.g., the 3 vertices of the red triangles), and the lattice primitive vectors are $a_{1}=(\sqrt{2+\sqrt{2+\sqrt{3}}},0)$, $a_{2}=(-\sqrt{(2+\sqrt{2})/3},\sqrt{(2+\sqrt{2+\sqrt{3}})/3})$, as shown in the figure. With nearest neighbor bonds only, this lattice has $z=4$ and is Maxwell. When next-nearest neighbor bonds are included, each site has an extra coordination of 4 above the Maxwell point.

For Fig. 2 in the main text, we used this lattice with only nearest neighbor bonds, and we took all spring constants to be 1. The topological polarization is $R_{T}=a_{2}$, implying that the top edge has two zero modes per edge unit cell, while the bottom edge is rigid with no zero mode.

For Fig. 3 in the main text, we used a same lattice but next-nearest-neighbor springs are introduced, which makes the system over constrained. The spring constant of the next-nearest-neighbor springs are set to $1/1000$, which is much smaller than the spring constant of the nearest-neighbor springs ($1$). Because the next-nearest-neighbor springs are weak, this system is not far from the ideal Maxwell point, and thus the top edge contains two soft modes, whose frequencies are no longer exactly zero, while the bottom edge remains rigid.

Following the procedure shown above, we computed the elastic constants for the continuum theory of these two models. For the Maxwell lattice (Fig. 2), we used two different continuum model, one with only two acoustic bands and the other contains also a soft optical band. Both of them yield topological index agreeing with the discrete lattice theory.

For the first model where we keep two acoustic bands [Fig. 2(a)], the elastic energy is

[TABLE]

where $\lambda_{1}=0.731617$ and $\lambda_{2}=0.130026$, and

[TABLE]

At small $q_{x}$, the $q_{y}$ of the two edge modes have the following asymptotical form

[TABLE]

and the imaginary part of $q_{y}$ is plotted as the dotted lines in Fig. 2(a).

For the second model where we also include the lowest optical band [Fig. 2(b)], the elastic energy is

[TABLE]

where $\lambda_{1}=0.922319$, $\lambda_{2}=0.476694$ and $\lambda_{3}=0.126092$ and

[TABLE]

Here $\phi$ is a softest optical phonon mode. In contrast to the elastic theory with only acoustic modes, here, the $\partial_{y}u_{x}$ term and the $\partial_{x}u_{y}$ term are not symmetric and they have different coefficients. This is expected because we included an optical phonon mode. Without optical modes, $\partial_{y}u_{x}$ and $\partial_{x}u_{y}$ must have the same coefficient, because the antisymmetric combination $\partial_{y}u_{x}-\partial_{x}u_{y}$ corresponds to rigid body rotation and thus shall not arise in the elastic energy. In our elastic energy above, this is still true and the antisymmetric part $\partial_{y}u_{x}-\partial_{x}u_{y}$ can be absorbed into a redefinition of $\phi$, if we replace $\phi$ with $\widetilde{\phi}=\phi+0.274063(\partial_{y}u_{x}-\partial_{x}u_{y})$. After this redefinition, the $\partial_{y}u_{x}$ term and the $\partial_{x}u_{y}$ term will become symmetric.

At small $q_{x}$, the $q_{y}$ of the two edge modes have the same asymptotical form as shown in the two bands model

[TABLE]

and the imaginary part of $q_{y}$ is plotted as the dotted lines in Fig. 2(b). The agreement of the two-band and three-band model at the small $q_{x}$ limit serves as a consistency check for our calculations.

After the next-nearest-neighbor springs are introduced, the two-band continuum theory becomes

[TABLE]

with $\lambda_{1}=0.734928$, $\lambda_{2}=0.131478$, and $\lambda_{3}=0.00869398$. Because the $\lambda_{3}$ is very small, we will drop it as the zero order approximation, i.e., $r_{M}=\lambda_{3}/\lambda_{2}=0.066125$. This zeroth order theory is at the Maxwell point and thus we can perform topological analysis as shown above. And afterwards, we use the variational method to compute the energy of these topological modes. The square root of the energy, which is the frequency, is shown in Fig. 3(a).

In Fig. 3(b), we included an optical mode and the continuum elastic energy becomes

[TABLE]

where $\lambda_{1}=0.925252$, $\lambda_{2}=0.478491$, $\lambda_{3}=0.1272$, and $\lambda_{4}=0.00470203$ and $r_{M}=\lambda_{4}/\lambda_{3}=0.0369656$. Because $r_{M}$ is small, we can ignore $\lambda_{4}$ as the leading order theory, which is at the Maxwell point. The topological indices and zero modes are obtained from this leading order theory, whose energy is then computed using the variational method. Its square root, i.e., the frequency, is plotted in Fig. 3(b).