Topological mechanics [1] is a rapidly growing research field that studies classical analogs of topological effects in quantum many-body physics [2]. A prime example are spectrally gapped mechanical systems that can host topologically protected zero modes at their boundaries [3–5], similar to localized electronic excitations in the quantum spin Hall effect [6]. Another important class of examples are solid- or fluid-mechanical systems with broken time-reversal symmetry, which can exhibit chiral edge modes at finite frequency [7–10], analogous to the (anomalous) quantum Hall effect [11, 12]. Because these edge modes are topologically protected and robust against the introduction of defects, they may provide a powerful tool for the resilient localized transmission of sound signals in elastic materials [13].

Over the last 5 years, substantial progress has been made in the understanding of topological phenomena in a wide variety of classical systems, ranging from mechanical systems with lattice symmetry inspired by quantum analogs [14, 15] and amorphous networks [9, 16] to active systems [7, 17, 18], electrical circuits [19–21], and even ocean waves [22]. Many of the recently discovered mechanical topological insulators rely on a known underlying lattice structure [7, 14] or curvature [22] to induce the required gaps in their excitation spectra. From a practical perspective, it would be interesting to design and build more general structures with desired topological properties.

Complementing recent work aimed at engineering continuum topological insulators [23], we consider here the design of topological excitations in bandgap-optimized [24] mass-spring networks (MSNs). Specifically, we will demonstrate that MSNs with an inversely designed bandgap can host topologically protected finite-frequency edge modes, and convert non-robust non-topological edge modes into robust topological edge modes when time-reversal symmetry is broken. While many traditional topological materials, including those based on a hexagonal lattice like the Haldane model [25], do not possess the mode conversion property, this desirable feature is frequently encountered in our designed MSNs.

In the remainder, we focus on the dynamics of periodic crystals of 2D mechanical balls-and-springs networks. In all cases, the spring stiffnesses of these MSNs were numerically tuned such that the excitation spectrum exhibits a band gap (using the algorithm introduced in [24]). In formal analogy with quantum Hall systems [11, 12], we will then break time reversal symmetry by placing the MSN into a rotating frame, with the Coriolis forces acting equivalently to an external magnetic field. To study and characterize the topological phase transition and the emerging protected chiral edge modes in detail, we will (i) numerically calculate the non-zero Chern invariant associated to the topological phase, (ii) demonstrate the dynamics of the localized chiral edge excitations in numerical simulations, and (iii) explicitly identify those dynamical edge modes that are related to topological protection. The underlying inverse-design framework [24] uses a generic linear response optimization and is, therefore, broadly applicable. Promising candidates for experimental implementations are mechanically coupled phase oscillators under Coriolis acceleration, such as hydrodynamic spin lattices of walking droplets [26] or gyroscopic mechanical metamaterials [8, 10].

MSNs provide a generic modeling framework for many physical systems. The potential energy of an MSN with E springs is given by

where k_e is the stiffness of spring e, ℓ_e is its current length and ℓe(0) is its preferred rest length. Here, we are interested in the dynamics near the equilibrium configuration where all springs are at their rest lengths, ℓe=ℓe(0), corresponding to the masses being at positions xi(0). Expanding in small deviations ui=xi-xi(0) and neglecting frictional effects, we obtain the linearized equations of motion,

where K is the stiffness matrix of the network and m is the mass of the balls (we assume identical masses throughout). The vector u generally has dN components, where d is the dimension of space and N is the number of masses. From now on, we specialize to the case d = 2. The stiffness matrix can be further decomposed as K=Qk^Q⊤, where Q is the equilibrium matrix encoding the network geometry and k^=diag(k1,k2,…,kE) is the diagonal matrix of spring stiffnesses [27]. Neglecting thermal fluctuations throughout, the subsequent discussion focuses on macroscopic topological metamaterials, similar to those realized experimentally in Chen et al. [28]. In principle, it is possible to incorporate thermal or non-thermal noise [18] and/or more general nonlinear potentials, such as in the elastic Lennard-Jones model [29], provided these admit linearizations in the form of Equation (1).

The MSN dynamics, specifically its harmonic response and its phononic modes, are encoded in the eigenmodes

where ω_j are the eigenfrequencies. If the network is a crystal consisting of N_c periodically repeated unit cells with lattice vectors R_ℓ, the dynamical problem can be simplified by performing a lattice Fourier transform [27],

where we decompose the rest positions x_i = R_ℓ + v_n. Here, ℓ indexes the unit cell and n indexes the degree of freedom within the unit cell. The wavevector k=∑i=1dbiNiKi with b_i ∈ ℤ lies in the first Brillouin zone, N_i is the number of unit cells in the ith dimension, and the reciprocal lattice vectors satisfy K_i · R_j = 2πδ_ij for the primitive lattice vectors R_j. The Fourier transform decouples the eigen-problem Equation (2) for different wavevectors k and leads to a phononic band structure ω_n(k).

The band structure ω_n(k) is of great interest both scientifically and for engineering applications because it efficiently encodes the elastic response of the infinite network. Specifically, band structure engineering allows for explicit tuning of wave propagation in acoustic materials and can be used to design, for instance, waveguides, acoustic cloaks, or selective sound suppression. Whereas in a generic band structure, the shape and frequency of the acoustic modes depends strongly on the details of the dynamics, topological modes are protected by an integer invariant, which cannot change through continuous changes of the interaction parameters.

Although the physical realizations of topological insulators are vast already in the quantum case [2, 12], the possible invariants and topological phases have been completely classified [30]. In linear topological mechanics, a similar scheme exists as long as the dynamical matrix is positive definite [1]. For our MSNs, this condition is always satisfied. In the following, we shall focus on one particular class of classical topological band structures for two-dimensional systems.

Topological band structure is intimately related to the theory of Berry phases, or geometrical phases [31]. While a full account of the underlying theory is beyond the scope of this paper, the fundamental result can be stated for a linear dynamical system iψ˙=H(r)ψ with Hermitian matrix H(r) which depends on some parameter r. The band structure topology is then encoded in the eigenstates of this effective ‘Hamiltonian’ H, and can be characterized by calculating an integer topological invariant, the Chern number. We now give a brief sketch of this calculation.

If the system is prepared in an instantaneous eigenstate H(r)u(r) = λ(r)u(r) and the parameter r is varied adiabatically along a closed curve C in parameter space, then the solution will always remain in the instantaneous eigenstate. After traversing the curve, the solution will pick up a phase factor eiγC with

This is the celebrated Berry phase with A(r) the Berry connection (superscript H denotes the Hermitian transpose). While the Berry connection changes under reparametrizations of the curve (gauge transformations), the phase is invariant up to 2π, and therefore in principle a physical observable. One particular parameter space of interest is the Brillouin zone of a crystal. In two dimensions, the BZ has the topology of a torus, such that any curve connecting k and k + K is closed (because wavevectors k and k + K are equivalent if K is a reciprocal lattice vector). By Stokes' theorem, Equation (3) can then be expressed as a surface integral independent of the curve,

where

is called the Berry curvature. Equation (4) defines the Chern number χ, which is an integer modulo 2π, and characterizes the eigenstates {u(k)}_k∈BZ. Thus, because in a crystal each eigenstate parametrized by the wavevector k corresponds to a band, it is possible to assign a topological Chern number χ_n to each band n. This Chern number does not change under perturbations of the matrix H(k), unless bands cross. Then, the eigenstates are no longer non-degenerate and the above analysis fails.

The Chern number defined by Equation (4) is nonzero only if the dynamics are not time-reversal invariant. If the system has time reversal invariance, Ω(k) is an odd function of k, and the integral over the Brillouin zone vanishes.

For systems with many bands and a gap between bands n′ and n′ + 1, the key insight [32] is then that one can associate an invariant to the gap itself, namely

which can only change if the gap closes due to a perturbation of H(k). The gap-Chern number C characterizes the bulk of a gapped crystal. Near a boundary to another gapped crystal with a different C or to the vacuum, the topology of the system must therefore change locally by closing the gap. This argument implies the existence of modes that are localized to the boundary between different topological phases and located in the bulk gap. Because these modes are tied to the bulk topological invariants, they are robust and must always exist, regardless of the specific shape of the boundary. We note that while for historical reasons the notion of adiabatic changes of parameters was invoked to define the Chern number, no actual adiabatic processes are necessary for it to exist, and it makes sense for any Fourier-transformed Hamiltonian.

For numerical purposes, the above integrals can be discretized while retaining their gauge-invariant characteristics [33]. This way, Chern numbers can be computed robustly and quickly with reasonably coarse discretizations of the Brillouin zone. In addition, any Chern number numerically computed in this way will automatically be an integer.

In the remainder, we demonstrate that such topologically protected edge modes can indeed exist in mechanical networks which have been tuned to exhibit bandgaps at specified frequencies, opening up an inverse-design pathway toward explicitly programmable topology.

There are many mechanical systems that possess topological gaps by virtue of their lattice structure. Here, we consider a different approach by tuning a desired gap into the spectrum of a mechanical network through numerical Linear Response Optimization (LRO) [24]. Starting from a basic lattice topology such as a triangular grid or a randomized unit cell topology defining mass points and springs (Figure 1), the spring stiffnesses k_e are numerically optimized to produce a gapped material between two desired bands. Applying the numerical LRO approach introduced and described in detail in Ronellenfitsch et al. [24], we minimize the average response of the network at frequency ω,

where Gω(k^)=(mω21-Qk^Q⊤)-1 is the linear response matrix to harmonic forcing with frequency ω and Tr(·) is the matrix trace. Numerically minimizing Equation (5) over the individual spring stiffnesses k^ while fixing a certain ω_n < ω < ω_n+1 for eigenmodes ω_n is then equivalent to maximizing a spectral gap between the nth and (n + 1)th eigenvalue. Generalizing from spectral gaps to bandgaps, since the Fourier transform is a linear map that block-diagonalizes Gω(k^), the trace in Equation (5) is replaced by a sum over the traces over the responses at each individual wavevector k, Gω(k^,k). For practical purposes, this sum is truncated, and only traces over a small number of wavevectors are actually used in the numerical optimization. To avoid the spring stiffnesses converging to either zero or infinity, we additionally impose bound constraints 0.1 ≤ k_e ≤ 1.0, and employ the Limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm [34] to perform the numerical optimization. Particle masses are set to unity (m = 1).

The above LRO approach generalizes to arbitrary network topologies and dimensions [24]. Throughout this paper, we will illustrate general ideas by focusing on three specific examples of bandgap-tuned networks: One with a regular triangular grid unit cell topology, and two different randomized unit cell topologies (Figures 1A,C,E). All three networks were optimized to exhibit a bandgap at some predetermined frequency. Despite some notable differences between them, their band structures all show features reminiscent of band inversion (Figures 1B,D,F), a characteristic that is often (but not always) present in topological band structures [35–37].

Adopting band inversion as an indicator for the potential existence of a topological transition, all that remains to do is to break time-reversal invariance of the system dynamics by introducing a suitable interaction. In the case of electronic systems, an externally applied magnetic field can provide such a symmetry-breaking interaction [12]. A classical formal counterpart considered in the remainder is the Coriolis force [38] which breaks the time reversal symmetry of the MSN dynamics when the mechanical network is placed in a rotating frame [14].

To sketch the general procedure for formulating the MSN dynamics in a rotating frame, we first consider a point mass in a harmonic potential with stiffness K confined to the x–y plane, and under the influence of a constant rotation perpendicular to the plane, Ω = (0, 0, Ω). Let x be the position of the point mass as measured from the rotational axis. Then, Newton's equations of motion in the rotating frame take the form

The Coriolis force is

where

encodes the cross product and we introduced the 2D vector x′. Similarly, the centrifugal force is -Ω∧(Ω∧x′) = Ω∧(ΩΓx′,0) = -Ω²(Γ²x′,0) = Ω²(x′,0). Clearly, the fictitious forces lie in the plane of rotation, so that from now on we can analyze the system in 2D. Dropping the primes, Equation (6) then yields the in-plane equations of motion

We can now generalize from a single particle to the full MSN dynamics by collecting all the x coordinates of the point masses in the network into the first N components of the 2N-component vector x, and all the y coordinates into the second N components. Then the matrix Γ takes the form

where 1 is the N × N identity matrix, and K now denotes the stiffness matrix such that Equation (7) remains formally unchanged.

We would like to express these equations in terms of small displacements around an equilibrium configuration. In doing so, we need to take into account that the equilibrium configuration is changed by the rotation. To find the new equilibrium positions x^* in the rotating frame, we set x¨=x˙=0 and solve for x^*,

Thus, a steady state exists unless the rotation frequency Ω² resonantly matches one of the eigenfrequencies of the stiffness matrix K. In the absence of resonance, we introduce displacements u = x − x^*, and find their equations of motion,

Here, the stiffness matrix was shifted due to the centrifugal force, and a new Coriolis term has appeared.

In the following, we further assume slow rotations compared to the smallest eigenmode of interest, typically the frequency of the gap, and neglect the term proportional to Ω² ≪ 1. This leads to the final equations of motion,

The Coriolis term proportional to u˙ is responsible for breaking time-reversal symmetry in this classical system (the transformation t ↦ −t maps u˙ to -u˙ but leaves all other terms invariant), analogous to the Lorentz force in a quantum electron gas [12]. Because the eigenmodes of Equation (8) cannot be computed directly by inserting a harmonic ansatz, we must resort to the equivalent first order system

where y = (u, w)^⊤, w=u˙, and D is the dynamical matrix, the eigenmodes of which can be readily computed.

Equation (9) can be brought into a form manifestly equivalent to the Schrödinger equation by introducing the change of variables [1],

where the matrix square root K is well-defined because K is positive-semidefinite. Under this change of variables, the dynamics becomes

where the “Hamiltonian” H is manifestly Hermitian. This form makes explicit the connection between classical mechanical and quantum systems, as now the machinery of quantum mechanics is applicable to Equation (10).

Below, we illustrate and analyze the generic consequences of time-reversal symmetry breaking via rotation for three distinct mechanical networks based on the inversely designed unit cells in Figure 1. We will see that the corresponding MSNs undergo a topological phase transition when the rotation frequency exceeds a critical value, resulting in topologically protected gapless modes that are exponentially localized at the boundary of samples.

The three mechanical networks from Figure 1 exemplify typical phenomena encountered with mechanical Chern networks. For each of them, a topological phase transition occurs at some finite 0 < |Ω_c| < 0.1, independent of the sign of Ω. This is due to the fact reversing the sign of the rotation frequency Ω is equivalent to reversing time t ↦ −t, and therefore mirrors the band structure, ω(k, Ω) = ω(−k, −Ω). In particular, this means that one can use the sign of Ω to control the unidirectional propagation of excitations: A wave packet will reverse direction when the sign of Ω is flipped. In the topological phase |Ω|> Ω_c, all of the considered networks have a gap-Chern invariant C = ±1, which we calculated using the numerical procedure outlined in Fukui et al. [33].

Generally, edge bands can be visualized by taking an infinite periodic crystal in 2D and restricting to a ribbon-like slice that is finite in one direction with open boundary conditions (Figure 2G). The resulting 1D crystal now possesses a one-dimensional band structure in which localized edge modes are directly visible. Mode localization can be measured by the participation ratio λ=(∑n|un|4)/(∑n|un|2)2 of the eigenvector u(k). The ratio λ is large if the mode is localized to few elements of the vector, and small if it is spread over many elements of the vector.

For all three example networks from Figure 1, the corresponding 1D crystals exhibit two bands of localized modes in the bulk gap in the topological phase |Ω| > Ω_c (Figures 2B,D,F). The two bands host wave packets with opposite group velocity v_g = dω/dk, and are localized at opposite edges of the semi-infinite ribbon system. They thus correspond to one single chiral edge excitation. The match between the bulk gap-Chern number C = ±1 and the number of edge excitations (more precisely, the difference between clockwise and counter-clockwise edge modes) is a direct manifestation of the celebrated bulk-boundary correspondence [32, 39].

We further note that although the existence of |C| protected edge bands is guaranteed in the topological regime |Ω| > Ω_c, this does not preclude unprotected edge states in the trivial phase |Ω| < Ω_c. To illustrate this fact explicitly, consider the example in Figures 2E,F. The band structure for Ω = 0 in Figure 2E is topologically trivial (C = 0) but exhibits features two localized edge bands, which are converted into the topologically protected bands in Figure 2F as one crosses the phase transition at finite |Ω| = Ω_c > 0.

All three networks analyzed in Figures 1, 2 have in common that they support only a single chiral edge mode, the direction of which can be reversed by changing the sign of Ω. Additional simulation scans suggest that this is typical of mechanical networks designed with the LRO scheme: Among all bandgap-designed networks that exhibited a topological transition, we never observed a case with |C| > 1. This empirical finding is consistent with results from previous studies which reported that larger Chern numbers are typically associated with materials that possess long-range interactions or with systems that are periodically quenched or driven [40]. Mechanical networks with long-range interactions could, in principle, be designed by introducing additional bonds that connect beyond the nearest neighbor unit cells. While certainly intriguing, such “non-local” networks are beyond the scope of the present study.

The wave packets hosted by the topological edge bands of our short-range MSNs can be excited dynamically by forcing a semi-infinite or a finite network near the boundary at a frequency inside the bulk gap. As a specific example showcasing this generic effect, we consider the mechanical network from Figure 1C and construct a finite realization consisting of 12 × 12 unit cells. To demonstrate the robustness of the topological modes, we remove three unit cells from the left side boundary to introduce a boundary perturbation (Figure 3A). We then numerically simulate the forced dynamics

where the forcing vector f = (1, 0, …, 0, 1, 0, …, 0)^⊤ is zero except for the x and y components of one single node near the bottom left corner. We pick Ω = 0.15 such that the network is in the topological phase, and ω = 0.71 inside the bulk gap. The window function h(t) = sin(πt/150)Θ(150 − t), where Θ(t) is the Heaviside Theta function, slowly turns on the forcing at t = 0, and turns it off entirely at t = 150. The forcing injects energy into the network at the frequency ω, which preferentially excites edge modes and creates a wave packet that travels unidirectionally along the edge of the network (Figures 3B–E, Supplemental Video 1). In particular, due to the topological protection of the edge modes, the precise shape of the boundary does not matter for the existence of these wave packets. Back-scattering modes are suppressed, and the wave packet is able to travel around the perturbation in the boundary (Figures 3B,C). As anticipated at the beginning of this section, the chirality of these wave packets is controlled by the sign of the rotation rate Ω (Supplemental Video 2). If the network is put in the topologically trivial regime, no edge modes exist and the energy injected by forcing does not create a chiral traveling wave packet (Supplemental Video 3).

The dynamical behavior described above is encoded in a set of eigenmodes u with Du = iωu that are exponentially localized to the boundaries of the system, and where ω lies in the bulk gap. For the three networks shown in Figure 1, we again constructed finite realizations consisting of many unit cells in a square array, and computed the eigenmodes of the finite dynamical matrix D from Equation (9). For all three networks, we identified modes inside the bulk gap which were then found to be localized at the boundary (Figures 4A–C). To demonstrate exponential localization in each case, we analyzed a slice of the eigenmodes in y direction. Plotting the logarithm of the average node displacement ui,x2+ui,y2 as a function of the x position of the node confirms an exponential decay of the node displacement with distance from the boundary (Figure 4D).

We have demonstrated the existence of topologically protected chiral edge modes in the gaps of inversely designed mechanical networks, and have characterized their dynamical properties. For the network realizations considered here, we found that band inversion near the gap was a robust predictor for a topological phase transition induced by sample rotation. The direction of rotation enables control over the chirality of the edge excitation, and topological protection of the edge excitations was confirmed in direct numerical simulations and through calculations of an appropriate Chern invariant.

We hope that the present work can serve as a stepping stone toward the precise inverse programming of topological features into discrete disordered metamaterials. Instead of constructing gapped materials on the basis of known lattices by using certain features of the band structure (e.g., band inversion) as indicators of potential topological transitions, we envision that Linear Response Optimization [24] may eventually allow the direct tuning of such properties by implementing the desired topological characteristics into the optimization objectives.

The datasets generated for this study are available on request to the corresponding author.

HR and JD conceived the study and wrote the paper. HR performed calculations, numerical simulations, and visualized results.