Signatures of Topological Phonons in Superisostatic Lattices

TL;DR

This paper investigates how topological surface phonons in idealized lattices manifest as finite-frequency modes in realistic superisostatic lattices, with implications for designing mechanical metamaterials.

Contribution

It introduces a study of finite-frequency topological surface phonons in superisostatic lattices with added interactions, bridging ideal models and real-world materials.

Findings

Finite-frequency surface phonons are present in superisostatic lattices.

Signatures of topological modes persist despite added interactions.

Results are relevant for designing mechanical metamaterials.

Abstract

Soft topological surface phonons in idealized ball-and-spring lattices with coordination number $z=2d$ in $d$ dimensions become finite-frequency surface phonons in physically realizable superisostatic lattices with $z>2d$. We study these finite-frequency modes in model lattices with added next-nearest-neighbor springs or bending forces at nodes with an eye to signatures of the topological surface modes that are retained in the physical lattices. Our results apply to metamaterial lattices, prepared with modern printing techniques, that closely approach isostaticity.

Figures

Equations

$$\delta {\mathrm{r}}_1(X)$$

$$\delta {\mathrm{r}}_2(X)$$

$$\delta {\mathrm{r}}_3(X)$$

$${\omega }^2{\mathrm{u}}^L={\mathrm{D}}_{10}{\mathrm{u}}^{L-1}+{\mathrm{D}}_{00}{\mathrm{u}}^L+{\mathrm{D}}_{01}{\mathrm{u}}^{L+1},$$

$$\text{Det}\big{[}\bm{{\rm D}}_{10}Z^{-1}+\bm{{\rm D}}_{00}+\bm{{\rm D}}_{01}Z-\omega^{2}\,\bm{{\rm\delta}}\big{]}=0\,,$$

$${\mathrm{u}}^L=\sum _{n=1}^4{A}_n{\mathrm{w}}_n{Z}_n^L{e}^{iqx-i\omega t},$$

$$B_{mn}=\sum _{k=1}^6{\left[{\mathrm{D}}_{00}+{\mathrm{D}}_{01}{Z}_n-{\omega }^2\delta \right]}_{mk}{w}_{n,k},$$

$$E=\frac{1}{2}\sum _{NN}{\left(s_b^{NN}\right)}^2+\frac{v}{2}\sum _{NNN}{\left(s_b^{NNN}\right)}^2,$$

$$s_b^{NN}={\mathrm{u}}_b\cdot {\mathrm{a}}_b,$$

$$\displaystyle\bm{{\rm C}}(\bm{{\rm q}})=\left(\begin{array}[]{cccccc}a_{1,x}&a_{1,y}&-a_{1,x}&-a_{1,x}&0&0\\ 0&0&a_{2,x}&a_{2,y}&-a_{2,x}&-a_{2,y}\\ -a_{3,x}&-a_{3,y}&0&0&a_{3,x}&a_{3,y}\\ -a_{4,x}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{1}}&-a_{4,y}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{1}}&a_{4,x}&a_{4,y}&0&0\\ 0&0&-a_{5,x}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{3}}&-a_{5,x}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{3}}&a_{5,x}&a_{5,y}\\ -a_{6,x}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{2}}&-a_{6,y}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{2}}&0&0&a_{6,x}&a_{6,y}\\ b_{1,x}&b_{1,y}&-b_{1,x}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{3}}&-b_{1,y}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{3}}&0&0\\ 0&0&b_{2,x}&b_{2,y}&-b_{2,x}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{1}}&-b_{2,y}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{1}}\\ b_{3,x}&b_{3,y}&0&0&-b_{3,x}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{1}}&-b_{3,y}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{1}}\\ -b_{4,x}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{2}}&-b_{4,y}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{2}}&b_{4,x}&b_{4,y}&0&0\\ 0&0&-b_{5,x}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{2}}&-b_{5,y}e^{-i\bm{{\rm q}}\cdot\bm{{\rm T}}_{2}}&b_{5,x}&b_{5,y}\\ -b_{6,x}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{3}}&-b_{6,y}e^{i\bm{{\rm q}}\cdot\bm{{\rm T}}_{3}}&0&0&b_{6,x}&b_{6,y}\end{array}\right)$$

$${\mathrm{a}}_1$$

$${\mathrm{a}}_3$$

$${\mathrm{a}}_5$$

$${\mathrm{b}}_1$$

$${\mathrm{b}}_3$$

$${\mathrm{b}}_5$$

$$\displaystyle\bm{{\rm C}}_{00}(q)=\left(\begin{array}[]{cccccc}a_{1,x}&a_{1,y}&-a_{1,x}&-a_{1,x}&0&0\\ 0&0&a_{2,x}&a_{2,y}&-a_{2,x}&-a_{2,y}\\ -a_{3,x}&-a_{3,y}&0&0&a_{3,x}&a_{3,y}\\ -a_{4,x}e^{-iq}&-a_{4,y}e^{-iq}&a_{4,x}&a_{4,y}&0&0\\ 0&0&0&0&a_{5,x}&a_{5,y}\\ 0&0&0&0&a_{6,x}&a_{6,y}\\ b_{1,x}&b_{1,y}&0&0&0&0\\ 0&0&b_{2,x}&b_{2,y}&-b_{2,x}e^{-iq}&-b_{2,y}e^{-iq}\\ b_{3,x}&b_{3,y}&0&0&-b_{3,x}e^{iq}&-b_{3,y}e^{iq}\\ 0&0&b_{4,x}&b_{4,y}&0&0\\ 0&0&0&0&b_{5,x}&b_{5,y}\\ 0&0&0&0&b_{6,x}&b_{6,y}\end{array}\right)$$

$$\displaystyle\bm{{\rm C}}_{01}(\bm{{\rm q}})=\left(\begin{array}[]{cccccc}0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&-a_{5,x}e^{iq/2}&-a_{5,x}e^{iq/2}&0&0\\ -a_{6,x}e^{-iq/2}&-a_{6,y}e^{-iq/2}&0&0&0&0\\ 0&0&-b_{1,x}e^{iq/2}&-b_{1,y}e^{iq/2}&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ -b_{4,x}e^{-iq/2}&-b_{4,y}e^{-iq/2}&0&0&0&0\\ 0&0&-b_{5,x}e^{-iq/2}&-b_{5,y}e^{-iq/2}&0&0\\ -b_{6,x}e^{iq/2}&-b_{6,y}e^{iq/2}&0&0&0&0\end{array}\right)$$

$$f=\frac{\lambda }{2}{u}_{ii}^2+\mu u_{ij}{u}_{ij}.$$

$$\lambda$$

$$\mu$$

$$f$$

$$+2K_{13}{u}_{xx}{u}_{xy}+2K_{23}{u}_{yy}{u}_{xy}.$$

$$B{K}_{11}$$

$$\displaystyle+3381681827843928199167638v^{3}+1357752557780034687526437v^{4}+12534744296675727990045v^{5}\Big{)}\,,$$

$$B{K}_{22}$$

$$\displaystyle+1115965371795836672225418v^{3}+493549308640875173434839v^{4}+5898796907214956230395v^{5}\Big{)}\,,$$

$$B{K}_{33}$$

$$\displaystyle+5156348205989863308398v^{3}+7211220231643636602303v^{4}+1105991231576795293275v^{5}\Big{)}\,,$$

$$B{K}_{12}$$

$$\displaystyle+1200964451099927268925858v^{3}+485931485533103278158399v^{4}+4410736097612180096943v^{5}\big{)}\,,$$

$$B{K}_{13}$$

$$\displaystyle-2468306508061388354586v^{3}+2760356898057381688611v^{4}+1222235230936930676739v^{5}\Big{)}\,,$$

$$B{K}_{23}$$

$$\displaystyle+1213564070972292895811718v^{3}+488017109696386616047611v^{4}+1164113231256862985007v^{5}\Big{)}\,,$$

$$B$$

$$\displaystyle+80387187958650600842997v^{3}+796574499957021428370v^{4}\Big{)}\,.$$

$$B{K}_{11}$$

$$\displaystyle+2973129913162249070121170v^{3}+1310973960535841431703559v^{4}+82311191396248539177999v^{5}\Big{)}\,,$$

$$B{K}_{22}$$

$$\displaystyle+151989412941048260843482v^{3}+103607540769847961350959v^{4}+20808603404027818221843v^{5}\Big{)}\,,$$

$$B{K}_{33}$$

$$\displaystyle+23376116240867084249846v^{3}+9184980163653051188739v^{4}+1988538118416508451247v^{5}\Big{)}\,,$$

Full text

Recent work Kane and Lubensky (2014); Lubensky et al. (2015); Mao and Lubensky (2018) laid the foundation for a theory, akin to the topological band theory of electronic materials such as quantum Hall systems Halperin (1982); Haldane (1983) and topological insulators Kane and Mele (2005a, b); Bernevig et al. (2006); Moore and Balents (2007); Fu et al. (2007); Hasan and Kane (2010); Qi and Zhang (2011), of topological mechanics of periodic ball-and-spring isostatic lattices with average coordination number $z$, under periodic boundary conditions, equal to twice the spatial dimension, $2d$. This theory predicts the existence of zero-energy surface-modes at every surface wavenumber with the number of these modes on different surfaces depending on the topological properties of the bulk phonon spectrum. It has been applied to a variety of systems and phenomena Stenull and Lubensky (2014); Paulose et al. (2015a, b); Sussman et al. (2016); Rocklin et al. (2016); Stenull et al. (2016); Chen et al. (2016); Meeussen et al. (2016); Baardink et al. (2018); Zhou et al. (2018a, b) from random and jammed systems to stress concentration at topological domain walls. Our focus here is on periodic fully gapped systems in which the only bulk zero modes are those imposed by translational invariance at wavenumber $\bm{{\rm q}}=\bm{{\rm 0}}$. Naturally occurring crystals always have an effective coordination number greater than $2d$ (because forces between sites have a range greater than the inter-site separation) or stabilizing bending forces favoring particular angles between bonds incident on a given site, and they are not candidates to exhibit topological mechanics. On the other hand, with the aid of modern printing and cutting techniques, metamaterials with $z=2d$ consisting of vertices connected by thin nearest-neighbor (NN) elastic beams can be designed Stenull et al. (2016); Baardink et al. (2018) and constructed Bilal et al. (2017); MaZhou2018 to minimize bending forces and, thereby, closely approach the isostatic limit to which the topological theory of Refs. Kane and Lubensky (2014); Lubensky et al. (2015); Mao and Lubensky (2018) applies.

Here we study generalized kagome lattices (GKLs) to which weak next-nearest-neighbor (NNN) springs or bending forces Mao and Lubensky (2011) are added [Fig. 1], and we focus on how their surface modes evolve as the magnitudes $v$ of these forces are increased from zero. In the presence of either such force, the originally isostatic lattices become stable elastic materials whose long-wavelength excitations are described by continuum elasticity, which predicts identical Rayleigh waves Landau and Lifshitz (1986) on opposite surfaces of a strip [see Supplemental Material (SM)]. It would be natural to expect that these Rayleigh waves evolve from zero-energy surface modes of the isostatic lattice, and this is indeed the case for non-topological lattices, which have the same number of zero modes on all pairs of opposite parallel surfaces Kane and Lubensky (2014). But this cannot be the case for topological lattices, which have opposite parallel surfaces with different numbers of zero modes –at the extreme no zero modes on one and an associated excess of zero modes on the opposite surface. In what follows, we discuss Rayleigh waves in weakly superisostatic lattices in the context of topological phonons, and we detail how the dilemma posed by the topological lattices is resolved.

For the sake of generality, we consider generic non-topological ($X_{\text{nt}}$) and topological ($X_{\text{t}}$) GKLs that have the lowest possible plane crystallographic symmetry, p1. To study surface modes, we assume that a free surface parallel to the $x$-axis exists as indicated in Fig. 1 so that the network as a whole is semi-infinite with a parallel opposite surface located at infinite distance. For the standard GKL with $v=0$, liberating these two surfaces from the constraints of periodic boundary conditions amounts to removing 2 bonds or 4 bonds and one site per surface unit cell. Both choices lead to two zero-surface-modes per surface wavenumber $q$ distributed on the combined lower and upper surfaces, but the latter, which we consider, has smoother upper and lower surfaces as shown in Fig. 1. The topological polarization $\bm{{\rm R}}_{T}$ Kane and Lubensky (2014) calculated from the bulk phonon spectrum is zero for $X_{\text{nt}}$, and it is non-zero and pointing towards the bottom surface, $\bm{{\rm R}}_{T}=-\frac{1}{2}(1,\sqrt{3})$, for $X_{\text{t}}$. As a consequence, there is one zero-surface-mode per $q$ on either surface for $X_{\text{nt}}$ and two (zero) zero-surface-modes per $q$ on the bottom (top) surface for $X_{\text{t}}$.

Turning to $v>0$, we will first review our results and then present some details about how we obtained them. Because of space constraints and for concreteness, we center our discussion on the case with NNN forces. Further details, model elastic energies etc., and results for the case with bending forces are provided in the SM. Figure 2 summarizes our major results about changes in the phonon band structure as the strength of the NNN coupling increases from zero and, in particular, how long-wavelength Rayleigh waves with the same speed develop on opposite surfaces and how the zero-energy surface states at $v=0$ evolve with increasing $v$. At $v=0.1$, both $X_{\text{nt}}$ and $X_{\text{t}}$ have one acoustic surface mode on each surface at each wavenumber $q$ in the surface Brillouin zone (SBZ). At small $q$, the modes reduce to the elastic Rayleigh waves with dispersion $\omega_{R}(q)=c_{R}q$ on opposite surfaces with the same surface $c_{R}$, predicted by elastic theory. The situation at $v=0.001$ is very similar to that at $v=0.1$ for $X_{\text{nt}}$ except that the acoustic surface-mode frequency $\omega_{s}(q)$ is smaller at every $q$, indicating an approach to a single zero mode at each $q$ on each surface as $v\rightarrow 0$. Figure 2 (g) shows that $c_{R}\propto\sqrt{v}$ as follows from the observation that $\omega_{s}^{2}(q)$ must be linearly proportional to some combination of spring constants and be equal to zero at $v=0$. The situation for $X_{\text{t}}$ is more complex. The bottom surface has an acoustic mode that stretches across the SBZ and reduces to the expected Rayleigh wave at small $q$ and, in addition, a low-frequency optical mode whose frequency, $\omega_{\text{opt}}$ is proportional to $\sqrt{v}$ across the SBZ and that vanishes into the continuum at a critical wavenumber $q_{0}$. The top surface, on the other hand, below the lowest bulk band only has a Rayleigh wave with the same velocity as that of the bottom surface, that disappears into the bulk continuum at a wavenumber that vanishes as $v\rightarrow 0$. The two zero-frequency modes of the topological $v=0$ lattice on the bottom surface are then the limits of the acoustic mode and the low-frequency optical mode. At $v=0$, the bottom of the band of bulk states, $\omega_{\text{band}}$, is proportional to $q^{2}$ rather than $q$ as can be calculated from the envelope of the bulk dispersion, which has the form $\omega^{2}_{\text{bulk}}=(q_{y}^{2}-\beta q^{2})+O(q^{4})$ Kane and Lubensky (2014); Lubensky et al. (2015) at small $q$. Thus equating $\omega_{\text{opt}}$ to $\omega_{\text{band}}$ yields $q_{0}\propto v^{1/4}$ in agreement with our numerical calculations.

The finite $v$ frequency dispersions of surface states in both the $X_{\text{nt}}$ and $X_{\text{t}}$ lattices (both with p1 symmetry) depicted in Fig. 2 are in general different on the top and bottom surfaces, as one would expect because opposite surfaces in lattices with such low symmetry are not equivalent. However, consistent with elastic theory, the small $q$ Rayleigh waves on both surfaces are the same and do not reflect p1 symmetry. All of the higher frequency modes do however. The high-$q$ frequencies of the acoustic modes of both lattices are different on the two surfaces as are all the higher-frequency optical surface modes [see SM].

The approach of the finite-frequency phonons to the topological phonons is also reflected in their inverse penetration depths $\kappa$ shown in Figs. 2 (e) and (f). For both $X_{\text{nt}}$ and $X_{\text{t}}$, $\kappa$ is the same for $v>0$ at sufficiently small $q$ on the bottom and top surfaces as predicted by elastic theory, but differences between the bottom and top surfaces arise as $q$ becomes larger. As observed in the dispersion curves, the $v\to 0$ limit unfolds differently in the two lattices. For $X_{\text{nt}}$, the $\text{Re}(\kappa(q))$ curves of the two surfaces approach one another as $v$ vanishes and eventually become identical across the entire SBZ. For $X_{\text{t}}$, the $\text{Re}(\kappa(q))$ curves of the top surface terminate at values of $q$ that decrease with $v$ whereas the $\text{Re}(\kappa(q))$ curves of the acoustic and the lowest optical mode on the bottom surface approach each other to produce a two-fold degenerate zero-frequency mode at $v=0$. Note that the penetration depth of the most dominant contribution to this mode diverges for $q\to 0$. The inset to Fig. 2 emphasizes the extremely small (but which we have verified is nonetheless positive) value of $\text{Re}(\kappa(q))$ throughout the region that the surface acoustic mode exists on the top surface indicating a very large penetration depth.

Our results for the GKL with bending forces are very similar [see SM]. The only notable difference is that the interaction strength $v$ is effectively larger than in the NNN model due to factors mandated by the rotational invariance of the bending energies. Apart from that, the approach of the finite-frequency phonons to the topological phonons is qualitatively the same.

We now outline how these results were obtained. The GKLs are derived from the standard kagome lattice (KL) by displacing Kane and Lubensky (2014) the 3 KL unit cell sites $\bm{{\rm r}}_{1}=(0,0)$, $\bm{{\rm r}}_{2}=(1/2,0)$, and $\bm{{\rm r}}_{3}=(1/4,\sqrt{3}/4)$ by

[TABLE]

where $X\equiv(\chi_{1},\chi_{2},\chi_{3})$ 111Our convention is equivalent to that of Ref. [1,2] with a change in the signs of the $\chi$’s, so that the topological polarization is $\bm{{\rm R}}_{T}=-\frac{1}{2}\sum_{\mu}\bm{{\rm T}}_{\mu}\,\text{sign}\chi_{\mu}$, where $\bm{{\rm T}}_{\mu}$ are the primitive translation vectors.. $\bm{{\rm a}}_{b}$ are the normalized NN bond vectors of the KL: $\bm{{\rm a}}_{1}=(1,0)$, $\bm{{\rm a}}_{2}=1/2\,(-1,\sqrt{3})$, $\bm{{\rm a}}_{3}=1/2\,(-1,-\sqrt{3})$. $\bm{{\rm e}}_{b}$ are unit vectors perpendicular to the $\bm{{\rm a}}_{b}$: $\bm{{\rm e}}_{1}=(0,-1)$, $\bm{{\rm e}}_{2}=1/2\,(\sqrt{3},1)$, $\bm{{\rm e}}_{3}=1/2\,(-\sqrt{3},1)$. The displacements are designed Kane and Lubensky (2014); Lubensky et al. (2015) so that making one of the $\chi_{b}$’s nonzero causes filaments (i.e., sample traversing straight lines of bonds) parallel to $\bm{{\rm a}}_{b}$ to zigzag while keeping the remaining filaments straight. The crystallographic symmetry of the resulting GKL depends on $X$. For example, the twisted KL with $X=(\chi,\chi,\chi)$ (where $\chi$ is some reasonable positive or negative number) has p31m symmetry [see SM]. For $X=(0,\chi,\chi)$ and $X=(-\chi,\chi,\chi)$, the symmetry is reduced to cm and pm, respectively. Our generic GKLs have deformation parameters $X=X_{\text{nt}}=(0.1,0.15,0.2)$ and $X=X_{\text{t}}=(0.1,0.15,-0.2)$. We have chosen these parameters so that the resulting GKLs have the lowest possible (p1) symmetry, and moderate distortions relative to the KL. Otherwise, these choices are arbitrary, and manifolds of alternative choices lead to qualitatively the same results.

The vibrational modes of an elastic network are governed by its dynamical matrix $\bm{{\rm D}}$. In the bulk GKL, the equation of motion is simply $\omega^{2}\,\bm{{\rm u}}(\bm{{\rm q}})=\bm{{\rm D}}(\bm{{\rm q}})\,\bm{{\rm u}}(\bm{{\rm q}})$, where $\bm{{\rm u}}=(u_{1x},u_{1y},u_{2x},u_{2y},u_{3x},u_{3y})$ is the displacement vector of the basis sites, $\bm{{\rm q}}$ is the wave vector, and $\omega$ is the angular frequency. $\bm{{\rm D}}=\bm{{\rm Q}}\bm{{\rm S}}\bm{{\rm C}}$ is the $6\times 6$ lattice dynamical matrix (for unit mass at sites), with $\bm{{\rm Q}}$ the equilibrium matrix, $\bm{{\rm C}}=\bm{{\rm Q}}^{\dagger}$ the compatibility matrix, and $\bm{{\rm S}}=\text{Diag}(1,1,1,v,v,v)$ the spring constant matrix (see Ref. Lubensky et al. (2015) for background information). In the elastic (continuum) limit, $\bm{{\rm u}}$ turns into a 2-component displacement field and $\bm{{\rm D}}$ turns into a $2\times 2$ effective dynamical matrix. Details about the dynamical matrices in the two theories are given in the SM.

To get a comprehensive picture, we use both lattice and elastic theory. In our elastic theory, we adapt the standard textbook calculation Landau and Lifshitz (1986) of the decay lengths and sound velocities of acoustic surface phonons in isotropic continua to our anisotropic GLKs [see SM]. This approach applies only to the longest wavelength acoustic phonons. Our lattice-based calculations are a generalization to discrete lattices of the standard Rayleigh-wave continuum calculations Landau and Lifshitz (1986). Like the latter calculations, they are done on semi-infinite systems that clearly separate top and bottom surfaces, yet they allow access to wave vectors ranging across the entire SBZ. To carry out our calculations, we break the lattice into one-cell-thick layers $L$, with $L=0$ the surface layer, $L=1$ the next layer into the bulk, and so on, stacked in the $y$-direction and with periodic boundary conditions along $x$. The equilibrium matrix has non-vanishing components $\bm{{\rm Q}}_{L,L}\equiv\bm{{\rm Q}}_{00}$ and $\bm{{\rm Q}}_{L,L-1}\equiv\bm{{\rm Q}}_{10}$ connecting sites in layer $L$ to bonds in layers $L$ and $L-1$, respectively; and the compatibility matrix has non-vanishing components $\bm{{\rm C}}_{LL}\equiv\bm{{\rm C}}_{00}$ and $\bm{{\rm C}}_{L,L+1}\equiv\bm{{\rm C}}_{01}$ connecting bonds in layer $L$ to sites in layers $L$ and $L+1$, respectively. The dynamical matrix then has components $\bm{{\rm D}}_{L,L-1}=\bm{{\rm D}}_{10}=\bm{{\rm Q}}_{10}\bm{{\rm S}}\bm{{\rm C}}_{00}$, $\bm{{\rm D}}_{LL}=D_{00}=\bm{{\rm Q}}_{00}\bm{{\rm S}}\bm{{\rm C}}_{00}+\bm{{\rm Q}}_{10}\bm{{\rm S}}\bm{{\rm C}}_{01}$, and $\bm{{\rm D}}_{L,L+1}=\bm{{\rm Q}}_{00}\bm{{\rm S}}\bm{{\rm C}}_{01}$. The equation of motion for any layer $L>0$ then reads

[TABLE]

which is solved by $\bm{{\rm u}}^{L+1}=Z\,\bm{{\rm u}}^{L}$ provided that

[TABLE]

where $\bm{{\rm\delta}}$ is the unit matrix, and $Z$ determines the the inverse decay length $\kappa$ in the $y$-direction via $Z=\exp(-\kappa)$ (with $\kappa$ in general complex). Solutions $Z(v,\omega,q)$ of Eq. (3) come in pairs with reciprocal magnitude. Solutions with $|Z(v,\omega,q)|=1$ correspond to bulk modes, whereas solutions with $|Z(v,\omega,q)|<1(>1)$ decay away from the bottom (top) surface and correspond to surface modes. The points in $\omega$-$q$-space where bulk modes exist, i.e., points for which there is at least one pair of solutions with magnitude 1, form bands akin to the projected band structures in electronic systems (see Fig. 2). Surface modes can exist only within the bulk band gaps as the solutions of the equations of motion must obey the conditions imposed by the surface. Sites 1 and 2 of the surface unit cell lie directly in the free surface (note that site 3 does not). The force on these surface sites comes only from NN and NNN bonds $1$ to $4$ [Fig. 1] in the zeroth layer, and as a result for the free boundary condition we impose, the first four components of the force vector satisfy $\bm{{\rm f}}^{0}=\bm{{\rm D}}_{00}\bm{{\rm u}}^{0}+\bm{{\rm D}}_{01}\bm{{\rm u}}^{1}=\omega^{2}\bm{{\rm u}}^{0}$. For $v>0$, there are a total of eight zeros at any point in $\omega$-$q$-space. This implies that, at any point in a band gap, there are four modes with $|Z(v,\omega,q)|<1(>1)$ that decay away from the bottom (top) surface. The boundary conditions can be satisfied by superimposing these decaying modes,

[TABLE]

where the $A_{n}$ are mode amplitudes and $\bm{{\rm w}}_{n}=\bm{{\rm w}}_{n}(v,\omega,q,Z_{n})$ are polarization vectors. The band gap points for which the determinant of the $4\times 4$ boundary matrix $\bm{{\rm B}}$, defined by

[TABLE]

vanishes determine the dispersion relation of the surface modes. To find these points, we use the standard secant-method for computing zeros with an array of starting points that sweeps the band gaps.

Modern $3d$ printing and cutting techniques now produce bespoke materials, including regular lattices, with almost arbitrary designs. In particular, these techniques can produce mechanical lattices, whose geometry is almost identical to isostatic mechanical NN topological lattices. To fully understand and control these lattices, it is important to know how their properties - elastic energy, bulk- and surface- mode structure, etc. - differ from those of the ideal NN isostatic lattice. Our formalism treats semi-infinite systems exactly and can easily be used to calculate linearized response, for example to a localized force at a surface.

The main result of the present work is the unravelling of the apparent dilemma of topological lattices where topological-phonon theory predicts for the example we are studying two soft surface modes on one surface (soft, bottom) and zero on the other (hard, top) whereas elasticity theory mandates that there be one Rayleigh wave per surface wavenumber on each surface and that the two waves have equal speeds. Our work shows that the resolution of the dilemma is as follows: as $v\to 0$, the domain of existence of the Rayleigh wave on the top surface shrinks to zero. On the bottom surface, there is a low-energy optical surface mode, whose domain grows to the full SBZ and which approaches the bottom surface Rayleigh wave as $v\to 0$. These two together produce the two surface-zero modes predicted by topological-phonon theory.

Our results provide guidance for interpreting results of experiments on metamaterials targeting topological phonons. Reference MaZhou2018 reports experiments and finite element analysis on kagome-like lattices that show an asymmetric bulk phonon spectrum in a topological lattice but a symmetric one in a non-topological lattice. They verify the existence, in the same geometry we study, of the two low-energy surface modes on the soft surface that emerge from the two zero modes of the ideal topological lattice which they interpret as “an interesting departure from the conventional case of Rayleigh waves”. Curiously neither the finite element analysis nor the measurements show any evidence of the acoustic Rayleigh wave on the hard surface mandated by elasticity theory. It would be interesting to see additional experiments that specifically target the evolution of the hard surface Rayleigh wave with increasing bending rigidity.

I Supplemental Material

I.1 Symmetries of the GKL

The topological properties of Maxwell lattices, and our GKLs in particular, are not determined by their geometric symmetry, even though there are symmetry changes for lattices with $X=(\chi_{1},\chi,\chi)$ as $\chi_{1}$ changes sign as can be seen from (a) to (c) in Fig. 3. These lattices, the gapped non-topological [$X=(-0.1,-0.1,-0.1)$] and topological [$X=(0.1,-0.1,-0.1)$] lattices and the critical [$X=(0,-0.1,=0.1$)] lattice in which the gaps along $q_{y}$ vanish, all have different symmetries. All three of these lattices can, however, be continuously distorted into “generic” lattices with the lowest polar p1 symmetry, as shown in Figs. 3 (d) to (f), without changing their gap structure simply by allowing the magnitudes of $\chi_{1}$, $\chi_{2}$, and $\chi_{3}$ to be different. It should be noted that all topological lattices with a non-vanishing topological polarization have a geometric polar symmetry [wallpaper groups p1 or pm] but both non-topological and critical lattices can also have this symmetry. In the main text, we focused on the surface band structure of generic lattices.

I.2 GKL with NNN stretching forces

I.2.1 Model energy

To adapt the GKLs to the superisostatic situation typically found in the lab, we augment them here with NNN springs. This leads to the ball-and-spring model elastic energy

[TABLE]

where the first sum runs over the 6 NN bonds and the second sum over the 6 NNN bonds of the unit cell shown in Fig. 1 (a) of the main paper.

[TABLE]

is the stretch of NN bond $b$, where $\bm{{\rm u}}_{b}=\bm{{\rm u}}_{i}-\bm{{\rm u}}_{j}$ is the difference in the elastic displacements of lattice sites $i$ and $j$ connected by that bond which has a normalized bond vector $\bm{{\rm a}}_{b}$. The NNN-bond stretch $s_{b}^{NNN}$ is defined in a similar, obvious manner. For simplicity, we have set the spring constant of the NN bonds and the masses of the sites equal to 1.

I.2.2 Lattice theory – equilibrium, compatibility and dynamical matrixes

The equilibrium, compatibility and dynamical matrixes are elementary to the lattice description of elastic networks. For any $d$-dimensional central-force elastic network with $n$ sites and $n_{B}$ bonds, the $n_{B}\times dn$ compatibility matrix $\bm{{\rm C}}(\bm{{\rm q}})$ relates bond displacements $\bm{{\rm u}}(\bm{{\rm q}})$ to bond extensions $\bm{{\rm e}}(\bm{{\rm q}})$ via $\bm{{\rm C}}(\bm{{\rm q}})\bm{{\rm u}}(\bm{{\rm q}})=\bm{{\rm e}}(\bm{{\rm q}})$. The null space of $\bm{{\rm C}}(\bm{{\rm q}})$ constitutes the zero modes of the network. The $dn\times n_{B}$ equilibrium matrix $\bm{{\rm Q}}(\bm{{\rm q}})=\bm{{\rm C}}^{\dagger}(\bm{{\rm q}})$ relates bond tensions $\bm{{\rm t}}(\bm{{\rm q}})$ to site forces $\bm{{\rm f}}(\bm{{\rm q}})$ via $\bm{{\rm Q}}(\bm{{\rm q}})\bm{{\rm t}}(\bm{{\rm q}})=\bm{{\rm f}}(\bm{{\rm q}})$. Its null space constitutes the states of self-stress of the network. The $dn\times dn$ dynamical matrix governing the phonon spectrum is related to the equilibrium and compatibility matrixes by $\bm{{\rm D}}(\bm{{\rm q}})=\bm{{\rm Q}}(\bm{{\rm q}})\bm{{\rm S}}\bm{{\rm C}}(\bm{{\rm q}})$, where $\bm{{\rm S}}$ is the spring constant matrix.

The bulk compatibility matrix of our model lattice with NNN bonds reads

[TABLE]

in $\bm{{\rm q}}$-space, where $\bm{{\rm T}}_{1}=(1,0)$ and $\bm{{\rm T}}_{2}=\frac{1}{2}(1,\sqrt{3})$ are the primitive translation vectors we are using, and $\bm{{\rm T}}_{3}=\bm{{\rm T}}_{1}-\bm{{\rm T}}_{2}=\frac{1}{2}(1,-\sqrt{3})$. $\bm{{\rm T}}_{1}$, $\bm{{\rm T}}_{2}$, and $\bm{{\rm T}}_{3}$ are chosen so that their sum is zero. Note that the primitive translation vectors are independent of $X$. $\bm{{\rm a}}_{b}$ and $\bm{{\rm b}}_{b}$ are the normalized NN and NNN bond vectors of the GKL, respectively. These depend on the deformation parameters $X$. For $X=X_{\text{nt}}$, for example,

[TABLE]

and

[TABLE]

Calculating the nullspaces of the equilibrium and compatibility matrixes for $X_{\text{nt}}$ and $X_{\text{t}}$, we find that there are 8 states of self-stress and the 2 inevitable trivial zero modes for $\bm{{\rm q}}=\bm{{\rm 0}}$ which is consistent with the Maxwell counting. The dynamical matrix of the lattice theory is readily obtained by taking the product of the equilibrium and compatibility matrixes.

In the presence of a planar surface, it is useful to decompose equilibrium, compatibility and dynamical matrixes into layer matrixes describing springs respectively connecting sites in the same and ones in different surface-parallel layers. For our choice of having a free surface parallel to the $x$-direction, we have

[TABLE]

for the intra-layer compatibility matrix and

[TABLE]

for the extra-layer compatibility matrix, where $q=q_{x}$.

I.2.3 Elastic theory – Lagrange elastic energy

Under imposed external strain, basis sites $\alpha$ undergo displacements $u_{\alpha,i}=\eta_{ij}x_{\alpha,i}+\delta u_{\alpha,i}$ for $i=x,y$, where $\eta_{ij}$ is the imposed macroscpic deformation, and $\delta u_{\alpha,i}$ is the nonaffine part of the displacement. Minimizing our model elastic energy over $\delta u_{\alpha,i}$, we obtain an effective elastic energy density that can be expressed in terms as the usual Lagrange strain tensor $u_{ij}$.

For the conformations of the GKL with higher symmetry, the effective elastic energy can be very simple. The GKL with $X=(\chi,\chi,\chi)$ , for example, corresponds to the twisted kagome lattice which is macroscopically isotropic. Hence is Lagrange energy density is of the form

[TABLE]

The Lame coefficients of this lattice with $\chi=0.1$ are, e.g., given by

[TABLE]

Note that the bulk modulus $B=\lambda+\mu$ vanishes for $v\to 0$ as it should for the twisted kagome lattice without NNN bonds.

For our generic lattices $X_{\text{nt}}$ and $X_{\text{t}}$, the Lagrange energy density is considerably more complicated because there are six independent elastic constants:

[TABLE]

After Fourier transformation of $f$, we can straightforwardly extract the dynamical matrix for $X_{\text{nt}}$ and $X_{\text{t}}$ in elastic theory by taking second derivatives with respect to the components of the elastic displacement.

For $X_{\text{nt}}$, the six elastic constants are given by

[TABLE]

For $X_{\text{t}}$, the elastic constants read

[TABLE]

I.2.4 Calculation of elastic Rayleigh waves in general anisotropic crystals

We seek Rayleigh waves on edges parallel to the $x$-axis and decaying exponentially into the bulk for $y>0$. The elastic dynamical matrix, $\overleftrightarrow{D}$, is homogeneous in $q_{x}$ and $q_{y}$, and we can assume that $q_{y}=pq_{x}\equiv pq$, where $p$ must have a positive imaginary part. In this case, we can scale $\overleftrightarrow{D}$ via $\overleftrightarrow{D}=K_{11}q^{2}\overleftrightarrow{R}$, where the components of $\overleftrightarrow{R}$ are

[TABLE]

where $k_{ab}=K_{ab}/K_{11}$. Then,

[TABLE]

where $\vec{u}=(u_{x},u_{y})$, $s=\omega^{2}/(K_{11}q^{2})$, and where

[TABLE]

determines the relation between $p$ and $\omega$. This is a quartic equation in $p$ whose solutions are either real or part of a complex-conjugate pair. Two decaying solutions, i.e., solutions with positive imaginary parts for $p$, are needed to meet the decay constraint and the surface boundary conditions, so in parameter regions where Rayleigh waves exist, there are two complex conjugate pairs. This means that two solutions have positive imaginary parts and two have identical negative imaginary parts implying that the Rayleigh waves on opposite surfaces will have exactly the same energy and penetration depths in spite of the fact that opposite surfaces are not equivalent in systems with polar p1 symmetry.

To determine $s$, and thus $p$, we impose the boundary condition of zero stress at the edge $y=0$:

[TABLE]

where $S_{ab}=\sigma_{ab}/(iqK_{11})$ is the reduced stress tensor. The solutions to Eq. (55) are

[TABLE]

where $a=x,y$, which when inserted into Eq. (57b) yield

[TABLE]

where

[TABLE]

The Rayleigh wave sound speed is determined by $|\det\overleftrightarrow{M}(\omega,p_{1}(\omega),p_{2}(\omega),q)|=0$. This program is easily implemented numerically.

I.2.5 Top surface and deformation parameters with flipped signs

The top surface of our model network can be conveniently studied by flipping the signs of the deformation parameters $X$ while keeping the surface at the bottom. For the top surface of $X_{\text{nt}}$, we can instead do our actual calculation for the bottom surface of $\bar{X}_{\text{nt}}=(-0.1,-0.15,-0.2)$. In the topological case, we can likewise use $\bar{X}_{\text{t}}=(-0.1,-0.15,0.2)$ instead of $X_{\text{t}}$. Comparison of Fig. 4 (a) [(b)] with Fig. 1 (c) [(d)] of our main paper demonstrates that the bottom surface of $\bar{X}_{\text{nt}}$ [$\bar{X}_{\text{t}}$] is equivalent to the top surface of $X_{\text{nt}}$ [$X_{\text{nt}}$] up to an inconsequential rotation of the entire system by $180^{\circ}$.

I.2.6 Full band structure

In the main text, our focus lies on the lowest frequency surface modes that approach that become topological zero modes in the isostatic limit. Our lattice calculation approach, however, allows us to go beyond to low-frequency limit and calculate the full surface mode structure. Figure 5 presents as an example the full surface mode structure of the NNN GKL with $X_{\text{nt}}$ and $X_{\text{t}}$ at $v=0.5$. Note that the dispersions of the optical surface modes are different for the bottom and top surfaces.

I.2.7 Inverse penetration depth

Each of our surface modes consist of a superimposition of four normal modes that decay away from the respective surface and hence each of our surface modes is associated with four $\text{Re}(\kappa(q))$ curves. In Fig. 2 (e) and (f) of the main text, we focus on the two longest ranging contributions to each surface mode, i.e., we display only the two lowest $\text{Re}(\kappa(q))$ curves for each surface mode to make the plots less busy. For the sake of completeness, we show here in Fig. 6 the full set of $\text{Re}(\kappa(q))$ curves pertaining to the surface modes in Fig. 2 (a) to (d) of the main text.

I.3 GKL with bending forces

I.3.1 Model energy

In the usual GKL, the lattice sites act as free hinges, i.e., there is no preferred angle between any pair of bonds that meet at a given site. Here, we extend the GKL to include bending energies that penalize deviations of bond-pair angles from their equilibrium values, see Fig. 1(b) of our main text. Our model elastic energy reads

[TABLE]

where the NN stretching contribution with the bond stretch $s_{b}=\bm{{\rm u}}_{b}\cdot\hat{\bm{{\rm a}}}_{b}$, where $\bm{{\rm u}}_{b}=\bm{{\rm u}}_{i}-\bm{{\rm u}}_{j}$ is the difference in the elastic displacements of lattice sites $i$ and $j$ connected by bond $b$, is identical to that of our model with added NNN forces. In the second term, the bending contribution, the sum runs over the 12 bond pairs associated with the angles defined in Fig. 1 (b) of our main text. $\theta_{B}$ measures the deviation of the angle of bond pair $B$ from its equilibrium value $\theta_{B}^{0}$, and $v$ is the bending stiffness. The specific form of $\theta_{B}$ depends on the value of $\theta_{B}^{0}$. For $\theta_{B}^{0}=0$, i.e., for a pair of bonds $b$ and $b^{\prime}$ that is straight in equilibrium,

[TABLE]

with $\bm{{\rm w}}_{b,\perp}=P_{b}\bm{{\rm w}}_{b}$. Here, $\bm{{\rm w}}_{b}=\bm{{\rm u}}_{b}/|\bm{{\rm x}}_{b}|$, where $|\bm{{\rm x}}_{b}|$ is the equilibrium length of the bond and $P_{b}=\bm{{\rm\delta}}-a_{b}a_{b}$, with $\bm{{\rm\delta}}$ the unit matrix, the projector on the direction perpendicular to it. For the generic GKLs that we focus on in our present work, all bond pairs are bent to some degree at equilibrium, $\theta_{B}^{0}>0$, so that

[TABLE]

for all bond pairs. Note that, by construction, $\theta_{B}$ is invariant under global rotations, and that the four $\theta_{B}$’s about any given node sum up to zero.

For our actual calculations, it is more convenient to rewrite the model elastic energy as

[TABLE]

where

[TABLE]

and

[TABLE]

Note that this effective bending stiffness is larger than the bare $v$. For our $X_{\text{nt}}$ and $X_{\text{t}}$ lattices, the average of $v_{B}$ over all 12 angles per unit cell is roughly 100 times larger than $v$. This must be taken into account when comparing results for the GKL with NNN and bending energies, respectively, see below.

I.3.2 Lattice theory – compatibility matrix

The bulk compatibility matrix of our model lattice in $\bm{{\rm q}}$-space reads

[TABLE]

where $\bm{{\rm a}}_{b}^{\perp b^{\prime}}=P_{b^{\prime}}\bm{{\rm a}}_{b}$ is the projection of the bond vector of bond $b$ onto the direction perpendicular to bond $b^{\prime}$. Note that this compatibility matrix is based on the model elastic energy as written in Eq. (64), i.e., the corresponding spring constant matrix has $1$ and the $v_{B}$ on its diagonal rather than $1$ and the bare $v$. From here on, the lattice and elastic theory calculations proceed exactly as for the NNN GKL.

I.3.3 Results

Having the compatibility matrix for the bending GKL, we can proceed exactly as for the NNN GKL. Inter alia, we can readily contract from it the bulk dynamical matrix of the lattice theory and then calculate the bulk spectrum. We can decompose it into the layer matrixes that provide the foundation for our lattice theory approach for calculating the surface modes. And, we can extract from it the 2 by 2 effective dynamical matrix of the elastic theory.

Figure 7 compiles our main results for the GKL with bending forces. It shows the low-frequency mode structure for $X_{\text{nt}}$ and $X_{\text{t}}$, the accompanying results for the inverse penetration depths and our our results for the sound velocities $c_{R}$ for both $X_{\text{nt}}$ and $X_{\text{t}}$, as well as for $X_{\text{t}}$ the vertical gap $\Delta\omega$ at $q=\pi$ between the acoustical surface mode and the lowest optical mode on the bottom surface and the onset $q_{0}$ of the latter.

As explained above, one should expect that a favorable comparison between the model lattices with NNN and bending forces requires a rescaling of $v$ because it gets, in the model with bending, effectively renormalized to larger values through factors stemming from the rotational invariance of the bending interaction. Comparing Fig. 7 to Fig. 2 of the main paper, we see that the bulk and surface mode frequencies for $X_{\text{nt}}$ are almost identical in both model lattices when $v$ is rescaled in the bending model by a factor of $10^{-1}$. The inverse penetration depths are also very similar in this case. For $X_{\text{t}}$, the results become very similar when $v$ by a factor that is closer to $10^{-2}$. The upshot is that apart from this trivial rescaling, the results for the GKL with NNN and bending forces are very similar, and the signatures of the topological phonons in both are qualitatively the same.