Soft topological modes protected by symmetry in rigid mechanical metamaterials

Topologically protected modes possess novel properties and extraordinary robustness stemming from their dual nature: these modes appear at the boundaries, yet are generated by bulk properties huber2016topological ; bertoldi2017flexible ; mao2018maxwell ; lu_topological_2014 ; ozawa2019topological ; ashida2020 ; shankar2020 ; hasan2010colloquium . First realized in electronic states von_klitzing_quantized_1986 ; thouless_quantized_1982 ; haldane_model_1988 ; kosterlitz_ordering_1973 ; hasan2010colloquium , this topological bulk-boundary correspondence has since been extended to the mechanics, acoustics, and photonics of structured matter rechtsman_photonic_2013 ; lu_topological_2014 ; ozawa2019topological ; nash_topological_2015 ; wang_topological_2015 ; peano_topological_2015 ; susstrunk_observation_2015 ; khanikaev_topologically_2015 ; paulose2015selective ; he_acoustic_2016 ; chen_tunable_2016 ; chen2016topological ; souslov_topological_2017 ; murugan_topologically_2017 ; dasbiswas_topological_2018 ; vila_observation_2017 ; trainiti_optical_2017 ; mitchell_amorphous_2018 ; souslov2019topological ; susstrunk_classification_2016 ; huber2016topological ; bertoldi2017flexible ; mao2018maxwell ; Tauber2020 . All of these systems are characterized by topological invariants, quantized numbers associated with a physical state. The type of invariants that a particular system can exhibit, or whether such topological character can exist at all, depends on its symmetries and has been classified for conventional, Hermitian Hamiltonians via the tenfold way ryu_topological_2010 . More recently, topological concepts have been generalized to open quantum and classical systems in the presence of external drive and dissipation using fundamental ideas from non-Hermitian physics, which nevertheless focus on square Hamiltonians yao_edge_2018 ; esaki_edge_2011 ; liang_topological_2013 ; lee_anomalous_2016 ; leykam_edge_2017 ; menke_topological_2017 ; xu_weyl_2017 ; gonzalez_topological_2017 ; xiong_why_2018 ; shen_topological_2018 ; hu_exceptional_2017 ; kawabata2019symmetry ; yoshida_exceptional_2019 ; brandenbourger2019non ; Baardink2020 ; scheibner2020 ; ashida2020 .

For the mechanics of ball-and-spring networks, Ref. kane_topological_2014 predicts a number of localized floppy (zero-energy) modes proportional to a local flux of a bulk topological polarization. These topologically protected modes have been realized in mechanical metamaterials along interfaces in one, two, and three dimensions, as well as at dislocation defects kane_topological_2014 ; chen2014nonlinear ; rocklin2016mechanical ; paulose2015topological ; bilal2017intrinsically ; baardink2018localizing . These realizations require a fine balance (called isostaticity) between the numbers of degrees of freedom and constraints (e.g., springs) to define the underlying topological polarization. Isostaticity allows for this topological invariant by enforcing a one-to-one mapping between degrees of freedom at sites and the bonds between them. These distinct quantities can then be related via a non-Hermitian rigidity matrix (analogous to a Hamiltonian), which is square only for isostatic systems. However, isostatic materials are inherently unstable guest2003determinacy making them susceptible to deformations due to thermal fluctuations. This makes realizing isostatic topological lattices in atomic, molecular, or colloidal crystals especially challenging.

For the non-isostatic case, topologically protected soft modes are a consequence of topological invariants distinct from topological polarization. For example, overconstrained (i.e., hyperstatic) systems can possess low-dimensional topological boundary modes at the corners of two-dimensional systems saremi_controlling_2018 . For modes at interfaces (instead of corners), Ref. roychowdhury_classification_2018 includes an exhaustive classification scheme evocative of the tenfold way but for non-square non-Hermitian Hamiltonians, which could be applied to design non-isostatic systems. However, the lack of a bulk-boundary correspondence principle for the topological classification in Ref. roychowdhury_classification_2018 leaves open the problem of realizing topologically protected interface modes in overconstrained systems.

In this letter, we focus on mechanically stable lattices, which are overconstrained, as are nearly all naturally occurring crystals. They are robust to thermal fluctuations and can be realized on atomic, molecular or colloidal scales. Generically, these crystals are not expected to have any soft modes. Building on the classification scheme in Ref. roychowdhury_classification_2018 , we design materials in which a certain symmetry class can guarantee the presence of soft modes at any interface between topologically distinct states. We dub this the generalized inversion symmetry and design a one-dimensional hyperstatic lattice which respects this symmetry. We show for the first time, using exact solutions and numerical calculations, that topological modes localized at interfaces between topologically distinct lattice configurations arise in an overconstrained mechanical system. Furthermore, we show that these topologically protected interface gap modes have low energies rather than the zero energy modes of Ref. kane_topological_2014 . We show that at interfaces for which the bulk band-gaps on the two sides of the interface are sufficiently different, these topological modes are absent. Our work contributes to the understanding of non-Hermitian topology by extending design principles for topological modes beyond square Hamiltonians to rectangular matrices.

Generalized inversion symmetry.—The linear deformations of a mechanical system may be described via a rigidity matrix $\mathcal{R}$, a linear map $\mathbf{e}=\mathcal{R}\cdot\mathbf{u}$ determined from the system geometry that maps the displacements of sites $\mathbf{u}$ (or more general degrees of freedom) onto the extensions of springs $e$ (or more general violations of constraints) and hence may be used to generate a potential energy. In real space, this rigidity matrix is real but not necessarily square. For periodic systems, the matrix can be written in Fourier space in terms of the wavenumber $k$, resulting in a block diagonal matrix with blocks $\mathcal{R}(k)$ as Laurent polynomials in powers of the phase factor $\exp(ik)$ with real coefficients $\mathbf{R}_{n}$: $\mathcal{R}(k)=\sum_{n}\mathbf{R}_{n}\,\exp(ikn)$. Zero edge modes appear at complex wavenumbers $k$, for which the phase factor $\exp(ik)$ becomes a general complex number $z$ lubensky2015phonons ; mao2018maxwell .

Generalized inversion symmetry is defined by the existence of a basis in which $\mathcal{R}(k)$ is real for real wavenumbers $k$, i.e., that there exist unitary matrices $U,W$ such that $U^{\dagger}\cdot\mathcal{R}(k)\cdot W$ is a real matrix for all real $k$. A consequence of this symmetry is that zero modes come in pairs: for every zero mode at complex number $z$, where $z=\mathrm{e}^{ik}$, there is also a zero mode at $z^{-1}$, as shown in the Appendix. In other words, a lattice with generalized inversion symmetry has equal numbers of zero modes localized on the left and the right interfaces. Lattices with generalized inversion symmetry can be classified into topologically distinct phases roychowdhury_classification_2018 , even when the number of constraints differs from the number of degrees of freedom.

The canonical example of a mechanical lattice with topologically protected modes: the one-dimensional chain of rotors and springs studied in Ref. kane_topological_2014 , does not obey generalized inversion symmetry. Its zero mode is localized on either the right or the left interface indicating the presence of a topological polarization.

The topological invariant that distinguishes lattices in the classification that we use is calculated from the Singular Value Decomposition (SVD, a generalization of the eigenvalue decomposition) of the rigidity matrix. In SVD, the Fourier-transformed rigidity matrix $\mathcal{R}(k)$ is written as $\mathcal{R}=\mathcal{U}\Lambda_{R}\mathcal{V}^{\dagger}$, with $\mathcal{U},\mathcal{V}$ unitary and $\Lambda_{R}$ a rectangular matrix with only non-negative so-called singular values along the diagonal. $\mathcal{R}(k)$ can be transformed into its SVD-flattened version $\mathcal{Q}(k)$ by replacing every nonzero element of $\Lambda_{R}$ by $1$. For an isostatic lattice, in the basis in which $\mathcal{R}(k)$ is real, its SVD-flattened version $\mathcal{Q}(k)$ is a real orthogonal matrix, which can be classified into topologically distinct classes according to the homotopy groups of such matrices. Even a hyperstatic lattice in which the number of constraints per unit cell exceeds the number of degrees of freedom per unit cell by one can be similarly classified. We do this by adding to the SVD-flattened rigidity matrix $\mathcal{Q}(k)$ a column that is orthogonal to all its column vectors and thus transforming $\mathcal{Q}(k)$ to be orthogonal.

Maxwell lattices.— So far, we have only considered the classification of rigidity matrices, in line with the matrix classification scheme from Ref. roychowdhury_classification_2018 . Now, we proceed beyond classification to realising rigidity matrices for topological materials with generalized inversion symmetry. For an isostatic lattice with generalized inversion symmetry and $2$ sites per unit cell, the SVD-flattened rigidity matrix $\mathcal{Q}(k)$ is equivalent to a two-dimensional rotation matrix. The topological invariant classifying such a lattice is the integer winding number of its rotation angle around the unit circle, as the wavenumber $k$ goes from $0$ to $2\pi$. At an interface where this topological invariant changes, there appear topologically protected gap modes, as shown in the appendix.

Similarly, for an isostatic lattice with generalized inversion symmetry and $N=3$ sites per unit cell, the SVD-flattened rigidity matrix $\mathcal{Q}(k)$ is equivalent to a three-dimensional rotation matrix. To calculate its topological invariant, we represent $\mathcal{Q}(k)$ by a point in a solid sphere of radius $\pi$ whose antipodal points are identified, where the radius vector of the point encodes the rotation angle in its magnitude and the rotation axis in its direction. Then the topological invariant is $0$ or $1$ depending on the contractibility of the loop traced out by $\mathcal{Q}(k)$, as the wavenumber $k$ goes from $0$ to $2\pi$. At an interface where this topological invariant changes, there appears a topologically protected gap mode, as described in the Appendix. Maxwell lattices obeying generalized inversion symmetry with more than $3$ sites per unit cell are similarly characterized by a $\mathbb{Z}_{2}$ topological invariant.

Topological modes in a hyperstatic lattice.—The study of topological modes has hitherto been almost exclusively in the realm of Maxwell lattices. As we proceed to show, the advantage of generalized inversion symmetry is that it allows us to construct hyperstatic lattices that have topologically protected modes. We restrict our analysis to hyperstatic lattices with one degree of freedom and two constraints per unit cell. For these lattices, the topological invariant has a simple interpretation as the winding number around the origin of the real two-component vector $\vec{R}(k)$ as $k$ crosses the Brillouin zone, as shown in Fig. 1(b,d). For instance, consider a hyperstatic lattice with the following Fourier transformed rigidity matrix: $\mathcal{R}(k)=\left(c-\cos{k}\,,\,\sin{k}\right)^{T}$, where $c>0$ is a dimensionless parameter determined via the structure’s geometry. Its winding number is $1$ or $0$ for $c<1$ and $c>1$ respectively. When $c=1$, the loop traced out by $\vec{R}(k)$ passes through the origin of the 2D plane, indicating that the lattice must have a bulk zero mode as it crosses over from one topological phase to the other.

The above rigidity matrix $\mathcal{R}(k)$ can be transformed to an equivalent rigidity matrix $\tilde{\mathcal{R}}(k)$ which is real in real space, as follows:

The rigidity matrix $\tilde{\mathcal{R}}(k)$ is realized by the hyperstatic rotor chain in Fig. 1(a,c) where $c=(a+2r\sin\theta)/(a-2r\sin\theta)$, $a$ is the lattice spacing, $r$ is the distance between the fixed point and the rotor head, and $\theta$ is the rotor angle measured from the vertical. The hyperstatic rotor chain with $c>1$ ($\Leftrightarrow\theta>0$) and $c<1$ ($\Leftrightarrow\theta<0$) belong to topologically distinct phases with winding numbers $0$ and $1$, respectively, as shown in Fig. 1.

The equation of motion for the angular displacement $u_{n}$ of a rotor at lattice site $n$ becomes:

where the dimensionless parameter $c$ is taken to be a function of the lattice site $n$ ($c\rightarrow c_{n}$). For a detailed derivation, see Appendix.

We now analytically and numerically study the modes localized at an interface between topologically distinct hyperstatic rotor chains. These interface modes lie in the gap, i.e., have energy lower than the minimum bulk mode energy $(c_{n}-1)^{2}$, but are not necessarily soft. The interface modes are guaranteed to be soft when the lowest bulk mode energy, $(c_{n}-1)^{2}$, is small. These interface modes are reminiscent of the zero mode that the hyperstatic rotor chain must pass through as it is continuously deformed from the topological phase on one side of the interface ($c>1$) to the topological phase on the other side of the interface ($c<1$).

Sharp interface.—At a sharp interface between topologically distinct phases of the hyperstatic rotor chain, the parameter $c$ jumps from below $1$ to above $1$ or vice-versa. Using normal modes of the generalized Bloch form $u_{n}(t)=u_{0}\,z^{n}\,\mathrm{e}^{i\omega t}$, we exactly solve the equation of motion, Eq. (20) (see Appendix for details). Using our exact solution, we calculate the values of $c$ on the left $(c_{L})$ and right $(c_{R})$ of the interface, for which the exact solution admits a localized interface mode. As shown in Fig. 2(d), we find that the existence of a topological interface mode in the gap is guaranteed only when the gap frequencies on both sides of the interface overlap. This overlap occurs when the energies for the lowest bulk modes on either side of the interface coincide: $(c_{L}-1)^{2}=(c_{R}-1)^{2}$. For the interface to obey this condition and be topologically non-trivial (i.e. $c_{L}\neq c_{R}$), it must satisfy: $c_{L}+c_{R}=2$. Taking $c_{L}=1+m_{0}\,,\,c_{R}=1-m_{0}$ and substituting in the exact solution, the energy of the localized mode is $\omega^{2}=m_{0}^{2}-m_{0}^{4}/(4-m_{0}^{2})$, which is lower than the lowest bulk mode energy: $(c-1)^{2}=m_{0}^{2}$. Furthermore, the mode decay rates are: $z_{L}=(2+m_{0})/(2+m_{0}-{m_{0}}^{2})\,,\,z_{R}=(2-m_{0}-{m_{0}}^{2})/(2-m_{0})$. That is, for $|m_{0}|<1$, the mode amplitude is right-growing on the left of the interface ($|z_{L}|>1$), and right-decaying on the right side of the interface ($|z_{R}|<1$), and hence localized at the interface.

To confirm our exact analysis, we numerically calculate the modes of a periodic hyperstatic chain of $N=4000$ rotors with two sharp interfaces. One of the interfaces with $c_{L}=1.3\,,\,c_{R}=0.7$ (i.e. $m_{0}=0.3$), is shown in Fig. 2(a). The topological mode localized at that interface is shown in Fig. 2(b), in red, with the mode decay rates in agreement with the calculated values: $z_{L}=1.04$ and $z_{R}=0.95$. The density of states shown in Fig. 2(c) shows the two topological modes localized at the two interfaces as having the lowest energy (in red), in agreement with the calculated value: $\omega^{2}=0.0879$. They also have the lowest participation ratio, where the participation ratio wegner_inverse_1980 of a normalized mode $u_{n}$ is: $PR=1/(N\sum_{n}|u_{n}|^{4})$, indicating their localizated nature. As the gap frequency increases, the topological interface mode persists even when the gap frequencies on the two sides of the interface are not equal, i.e. when $c_{L}+c_{R}\neq 2$, as shown in Fig. 2(d).

Smooth interface.—To study the modes localized at a smooth interface between the two topologically distinct phases of the hyperstatic rotor chain, we take the continuum limit of Eq. (20) with $c_{n}\to c(x)=1+m(x)$ and $u_{n}(t)\to u(x,t)$, $u(x,t)=u(x)\,\mathrm{e}^{i\omega t}$ to get

with $u^{\prime}$ etc. denoting spatial derivatives.

In the limit $m(x)\ll 1$, i.e., in the region where $c(x)$ is close to $1$, the above equation of motion becomes: $u^{\prime\prime}(x)+(\omega^{2}-m^{2}){u}=0$, bearing a close resemblance to the time-independent Schrödinger equation with the energy $E=\omega^{2}$, the potential $U(x)=m^{2}(x)=(c(x)-1)^{2}$, and $2M/\hbar^{2}=1$. In this analogy, a smooth interface corresponds to a potential well with a minimum value of $U(x)=0$ when $m(x)=0,\,c(x)=1$, and depth dictated by the asymptotic values of $U(x)=m^{2}(x)$ on either side of the well. The potential $U(x)$ is symmetric about its minimum when ${m_{L}}^{2}={m_{R}}^{2}$: the asymptotic values of $m^{2}(x)$ on the left and right side of the interface are equal. This occurs when the lowest bulk mode energies on the two sides of the interface are equal. For a topologically non-trivial interface, this requires that $c_{L}+c_{R}=2$, where $c_{L},c_{R}$ are the asymptotic values of $c(x)$ on the left and right sides of the interface. When this condition is satisfied, at least one bound state solution exists irrespective of the depth of the potential well kocher_criteria_1977 , ensuring a localized mode at the interface.

The analogy with the Schrödinger equation enables an exact solution of Eq. (31) for $m(x)=m_{0}\tanh(x/W)$ landau_quantum_1981 . For sufficiently large width $W$ and depth $m_{0}^{2}$ of the potential well, there are multiple localized modes at each interface, their number given by $\lfloor s(m_{0},W)\rfloor+1$, where $s(m_{0},W)=(-1+\sqrt{1+(4m_{0}^{2}\,W^{2})})/2$. The energy of the localized interface modes is $\omega_{n}^{2}=m_{0}^{2}-(s(m_{0},W)-n)^{2}/W^{2}$, for $n=0,1,2,\ldots,\lfloor s\rfloor$, which is less than the lowest bulk mode energy $m_{0}^{2}$.

To confirm our exact analysis, we numerically calculate the modes of a periodic hyperstatic chain of $N=4000$ rotors with $m_{0}=0.3$, $W=6$, and two smooth interfaces for which $m(x)=m_{0}\tanh(x/W)$, one of which is shown in Fig. 3(a). The variation of $c(x)=1+m(x)$ (yellow) and the effective potential $U(x)=m^{2}(x)$ (blue) across the interface is shown in Fig. 3(b). The two localized mode profiles at the interface plotted in Fig. 3(c) (red) are as predicted by the exact solution described in the Appendix. The density of states shown in Fig. 3(d) shows two soft modes (red) per interface having the lowest energies, at the values predicted by the exact solution. These modes also have the lowest participation ratios indicating their localized nature.

The case when the lowest bulk mode energies on the two sides of the interface are unequal, i.e. $c_{L}+c_{R}\neq 2$, corresponds to an asymmetric potential well in the analogy with the Schrödinger equation. An approximate criterion from Ref. kocher_criteria_1977 says that a localized interface mode exists if $W\,(m_{L}^{2}+m_{R}^{2})\gtrsim 2\sqrt{2|m_{L}^{2}-m_{R}^{2}|}$, where $m_{L}^{2}$, $m_{R}^{2}$ are the asymptotic values of $m^{2}(x)$ on the left and right sides of the interface, and $W$ is the interface width. The region of parameter space $\{c_{L},c_{R}\}=\{1+m_{L},1+m_{R}\}$ for which the above criterion is fulfilled is plotted for different values of $W$ in Fig. 3(e), showing that topological interface modes persist even when $c_{L}+c_{R}\neq 2$, with wider and deeper potential wells allowing for greater deviations from the symmetric case.

Conclusions—We have presented the first theoretical study of topological mechanical interface modes derived from an overconstrained structure, with exact solutions and numerical calculations. These models rely on generalized inversion symmetry to define topological invariants and create robust soft modes at interfaces between topologically distinct lattices. For isostatic topological structures, experimental proposals chen_triblock_2011 ; socolar_mechanical_2017 have yet to be realized on submicron scales. By contrast, the structures that we propose are overconstrained and otherwise rigid, which potentially makes them more accessible to fabrication via existing techniques at scales down to the submicron cha_experimental_2018 . Such an architecture will be robust to thermal fluctuations and hence amenable to miniaturization to the micro- and nano-scale. Designing topologically protected soft modes in overconstrained materials may lead to future applications from cushioning using soft regions bilal2017intrinsically to controlled failure at topological interfaces paulose2015selective .

Acknowledgments—A.S. acknowledges the support of the Engineering and Physical Sciences Research Council (EPSRC) through New Investigator Award No. EP/T000961/1. H.K. thanks Jeremy England for his support.