Gapping the Stacked Kagome Lattice

We examine the mechanics of metamaterial structures by using a lattice model for the displacements of nodes and strains of the links. We place a point mass at each vertex and connect the neighboring ones by linear springs. Such models capture the small-strain response of realistic structures that are either 3D-printed from soft polymers (35, 37) or assembled from construction sets (4) or laser-cut components (38). The mechanics of these ball–spring networks are captured via the linearized equation of motion 𝑿¨=−D𝑿, where 𝑿≡(𝒙1,…,𝒙N) is a d× N-dimensional vector containing the displacements of all N particles in d dimensions relative to their equilibrium positions. For a given lattice geometry, we calculate the dynamical matrix D, which relates the forces exerted by springs to displacements of the particles. For simplicity, we work in units, in which all particle masses and spring constants are one. In the linear regime, we can find the dynamical matrix by first relating the NB-dimensional vector of spring extensions 𝑺≡(s1,…,sNB) to the displacements via 𝑺=R𝑿. The rigidity matrix R contains dN×NB entries determined by the equilibrium positions of the particles and the connectivity of the lattice. In combination with Hooke’s law, the matrix R lets us calculate the dynamical matrix D via the relation D=R𝖳R (2, 11). From the relation 𝑺=R𝑿, we note that kerR contains soft modes [i.e., collective displacements that (to lowest order) do not stretch or compress any of the springs]. However, the matrix R𝖳 relates forces on particles 𝑿¨ to spring strains via 𝑿¨=−R𝖳𝑺. From this relation, one notes that kerR𝖳 contains combinations of spring tensions that do not give rise to particle forces. These configurations are dubbed states of self-stress, because they define load-bearing states, which put the lattice under a static tension (11).

To calculate soft modes (states of self-stress), we need to numerically solve the configuration-dependent equation R𝑿=𝟎 (R𝖳𝑺=𝟎). However, a mathematical result called the rank-nullity theorem lets us compute the difference between the number of zero modes N0≡nullR and the number of states of self-stress NB≡nullR𝖳 (2, 5, 11). This difference, the softness charge ν≡N0−NB, is given by the dimensionality of the R matrix ν=dN−NB and can only change if balls or springs are either added or removed. We focus on the special case of isostatic, or Maxwell, lattices defined by ν= 0 when the system is considered under periodic boundary conditions. These lattices are marginally rigid and exhibit a symmetry between zero modes and states of self-stress.

To design metamaterials based on simple, repeated patterns, we focus on periodic structures. Periodicity allows us to explicitly calculate zero modes and states of self-stress in a large sample. We begin with a highly symmetric, “undeformed” lattice and explore its configuration space by changing the positions (but not the connectivity) of the nodes. The specific geometry that we consider is illustrated in Fig. 1: this stacked kagome lattice has a coordination number z≡ 2NB/N= 6 and ν= 0; the lattice is isostatic. This lattice is based on the two-triangle unit cell shown in Fig. 1A, which corresponds to two unit cells of the kagome lattice stacked on top of each other. We deform the lattice via the orientations of the triangles, which are governed by the three angles (ϕ,θ,ψ) of rotation around the (x,y,z) axes, respectively. Within this structure, we explore a range for each of the three angles from −π/2 to π/2. In Fig. 1B, we show the deformed unit cell for a particular structure in this region, which corresponds to the choice (ϕ,θ,ψ)=(π/3,0,π/3). This unit cell builds the periodic geometry in Fig. 1C. We proceed to quantitatively show that this metamaterial has no Weyl lines.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Architecture of a 3D stacked kagome lattice. (A) We consider the doubled unit cell of a vertically stacked kagome lattice and rotate the top constituent triangle around three perpendicular axes by angles (ϕ,θ,ψ) [the bottom triangle is rotated by (−ϕ,−θ,−ψ)]. (B) For the choice (ϕ,θ,ψ)=(π/3, 0,π/3) of these parameters, we obtain a unit cell that corresponds to a (gapped) 3D lattice with a topological polarization P and without Weyl lines. (C) Architecture of the (gapped) metamaterial, with one unit cell emphasized in the lower left corner and a 2D kagome sublattice surface highlighted in green along the top.

We calculate the spectrum of normal modes by considering plane-wave solutions of the form xα,ℓ=eiℓ⋅𝒌xα, where the index α refers to a particle within one unit cell and the lattice index ℓ enumerates different unit cells within a lattice. The three-component wave vector k is periodic in each component with −π≤kx,ky,kz≤π. This collection of points forms the first Brillouin zone of the lattice. The Bloch representation of the rigidity matrix in this plane-wave basis is Rαβ(𝒌)=Rαβei(ℓα−ℓβ)⋅𝒌. A zero mode at wave number k is a vector of unit cell displacement xα that solves Rβα(𝒌)xα= 0. Thus, within this setting, zero modes correspond to the zeros of the complex function detR(𝒌). By examining this function, we find the zero modes for different configurations of the lattice.

We numerically evaluate zero modes for the stacked kagome lattice (details are in Generalized Stacked Kagome Lattice and Soft Directions) and for most values of the angles (ϕ,θ,ψ), we find collections of zero modes along compact loops in k space (Fig. 2A). These Weyl loops, shown in Fig. 2 B and C, are analogous to 1D nodal lines (43⇓–45) and experimentally observed 0D Weyl points (46, 47) in 3D electronic semimetals and photonic crystals. In isostatic lattices, the number of Weyl loops is always even, because the materials’ time-reversal symmetry maintains 𝒌→−𝒌 reflection symmetry in the Brillouin zone—each Weyl loop comes in a pair with its reflected partner (15). Furthermore, these loops attach to the origin of the Brillouin zone along soft directions (i.e., the lines tangent to the Weyl loops at the origin). We count the number of loops by looking at soft directions in the neighborhood of the origin, and in Fig. 2A, we plot this count in a slice of parameter space. In this phase diagram, we note regions with up to six different loops. Strikingly, the middle of the diagram displays a region in which no Weyl loops exist (Fig. 2 A and D).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Phase diagram of stacked kagome lattices. (A) The number of soft directions in ϕ–ψ parameter space for θ= 0 (Fig. 1A). The region in dark blue corresponds to a gapped topological lattice, whereas the rest corresponds to lattices containing two (cyan), four (green), or six (orange) Weyl loops. (B) Weyl loops along the line (θ,ψ)=(0,π/3) corresponding to magnified points in A. As the gapped region is approached, we observe Weyl loops shrinking toward the origin until they disappear. (C) An example of the soft mode structure in reciprocal space for a lattice with a single two-Weyl loop structure (in blue), which corresponds to the cyan point (ϕ,θ,ψ)=(5π/12, 0,π/4) in A. The integration of a phase around any contour (e.g., solid red circles) enclosing the Weyl loop gives an integer winding number that we use to check the topological nature of the Weyl loop. Note that this contour can be deformed into a pair of contours (Cz; dashed red lines) traversed in opposite directions: this illustrates how the winding number across the Brillouin zone changes as the Weyl loop is traversed. Near the origin (kx,ky,kz)=(0, 0, 0), the loop manifests itself via soft directions (teal). (D) The Brillouin zone of a lattice corresponding to the blue point (ϕ,θ,ψ)=(π/3, 0,π/3) in A. At the origin, the softest regions (blue region) align themselves with the soft directions of lattices with (θ,ψ)=(0,π/3) (red) or (ϕ,θ)=(π/3, 0) (green).

We conclude that, although Weyl loops are generic, the stacked kagome lattice also exhibits gapped configurations in the sense that they contain no Weyl loops (but do contain gapless Goldstone modes). In lattices with Weyl loops, the number of soft modes scales as the linear size of the system, whereas gapped lattices have only the three uniform Guest–Hutchinson modes (12, 13). We are not aware of any other realizations of a 3D isostatic lattice that is gapped in this sense. In the next section, we address the implications of the existence of a gap for the topological characterization of mechanical networks.

Topological Rigidity in Three Dimensions

For the gapped isostatic lattice, an integer topological invariant called the polarization can be computed from the bulk phonon spectrum (3). This winding number invariant is only well-defined in the gapped case: in Weyl lattices, the closure of the gap prevents a consistent definition. That is, when Weyl loops are present, the polarization changes depending on the choice of contour in the Brillouin zone. The mechanical consequences of polarization are found in the spectrum of the material’s soft surface waves. Polarization controls which surfaces have more soft modes: this is the mechanical version of the bulk-boundary correspondence principle.

For the polarized lattice, a polarization vector P can be computed via 𝑷=∑imi𝒂i−𝒅0, where 𝒂i are the three lattice vectors and in our case, the unit cell dipole 𝒅0=−3𝒂3/2. This dipole is computed via the expression 𝒅0≡−∑b𝒓b, where 𝒓b are the positions of the bond centers relative to the center of mass of the unit cell (3). The three coefficients mi are winding numbers computed by integrating the phase change of the complex function detR(𝒌) across a straight-line contour crossing the Brillouin zone:mi≡−12πi∮CidkiddkilndetR(𝒌),[1]where Ci is a closed contour in the 𝒌^𝒊 direction (3). This integral is well-defined for contours along which detR(𝒌)≠ 0. The integers mi depend smoothly on the choice of contour and are constant. In gapped lattices, the determinant is nonzero everywhere except for the origin 𝒌= 0. As a consequence, the winding numbers are independent of the chosen contour, and the polarization is a topological invariant. However, Weyl loops partition the space of straight-line contours into lines going through the inside (outside) of the loop. Because detR(𝒌)= 0 along the loops, contours on either side of the loop can have different winding numbers mi. Note that the combination of any two such contours taken in opposite directions can be smoothly deformed (without intersecting the Weyl line) into a small circle enclosing the Weyl line as shown in Fig. 2C. The winding number mW around this small circle is an invariant and equal to the difference between mi for Ci on two sides of the loop. This topological protection guarantees that a single Weyl loop cannot be destroyed, which explains why Weyl loops only vanish by shrinking to the origin within the phase diagram in Fig. 2A. In summary, Weyl lines are protected by a winding number, and gapped lattices can have a well-defined topological polarization.

Bulk-boundary correspondence states that the topological invariants computed in the bulk can have significant effects on the mechanics of a sample with boundaries. We show this correspondence by computing the topological invariants and the spectrum of soft edge states in the stacked kagome lattice. In the bulk of the Weyl lattice corresponding to (ϕ,θ,ψ)=(5π/12,0,π/4) (compare Fig. 2 A and C), the winding number around the loop is mW=−1. The m3 for contours inside and outside of the loop differs by one. To see the boundary counterpart of the correspondence, we compute the spectrum of soft edge states for this lattice with a stress-free surface parallel to the xy plane (Fig. 1C). In Fig. 3 A and B, we show the (signed) inverse penetration depth κ for soft surface modes in the 3D Brillouin zone for wave vector (kx,ky). These plane-wave solutions of detR= 0 correspond to soft modes of the form eiℓ⋅(𝒌+i𝜿)𝑿α. The middle mode of Fig. 3A (zoomed in view shown in Fig. 3B) changes sides: for κ> 0 (Fig. 3A, red region), the soft mode is attached to the bottom surface and decays upward into the bulk, whereas for κ< 0 (Fig. 3A, blue region), the mode is localized on the top surface. For the line corresponding to the projection of the Weyl loop onto the 2D Brillouin zone, the penetration depth is infinite, and κ= 0, because the Weyl lines are soft modes in the bulk. The difference in m3 between inside and outside of the Weyl loops is the bulk invariant mW=−1, which corresponds to the difference between the number of soft surface modes across the projected Weyl loop. This connection follows from Cauchy’s argument principle for Eq. 1, which states that the third winding numbers m3 count, up to a constant, the number of modes with zero energy (detR= 0) at a boundary. These observations confirm bulk-boundary correspondence.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Surface modes for a polarized lattice (C and D) compared with a lattice with Weyl lines (A and B). A surface mode in a plane perpendicular to the z axis (e.g., Fig. 1C) decays exponentially into the material bulk with amplitude e−κz/|𝒂3|, where κ is the inverse penetration depth. Lattice with Weyl lines (ϕ,θ,ψ)=(5π/12, 0,π/4). (A) Inverse penetration depth for three soft surface modes over the entire surface Brillouin zone (sBZ). For small wave vectors (kx,ky), these modes correspond to the three acoustic phonons. (B) Zoomed in view of the middle mode in A: the part of wave vector space for which the mode lives on the bottom side is colored in red, whereas the part corresponding to the top side is in blue. The projection of the Weyl line from the bulk onto the sBZ corresponds to a line with infinite penetration depth (κ= 0) and separates the red and blue regions. Polarized lattice (ϕ,θ,ψ)=(π/3, 0,π/3). (C) Inverse penetration depth for the same three modes as in A but for a different choice of parameters, corresponding to the blue region in Fig. 2A. (D) Zoomed in view of the middle mode in C, which (in contrast to the mode in B) is localized on the top of the lattice only: it is blue throughout. Significantly, the bulk (κ= 0) modes that correspond to the Weyl loop are absent for this lattice.

A similar correspondence exists in the polarized lattice. There, the invariant polarization pointing along the z axis is given by 𝑷=𝒂3/2, which we computed for parameters (ϕ,θ,ψ)=(π/3,0,π/3). To understand the effect of P on the boundary, note that, in analogy with electromagnetism, the softness charge νS in region S is related to the flux of polarization P through the region’s boundary ∂S. In a nearly uniform, polarized lattice (3),νS=∮∂SdAVcell𝒏^⋅𝑷,[2]where 𝒏^ is the boundary’s inward normal and Vcell=det(𝒂1,𝒂2,𝒂3) is the volume of a unit cell. The difference between the softness charges at the top vs. the bottom is determined by P: for the stacked kagome, we expect one more band of soft modes along the top surface relative to the bottom. The total number of soft surface modes is three as a result of three bonds being cut per unit cell (3). We plot the inverse penetration depth κ for this extra band in Fig. 3 C and D within the 2D Brillouin zone. Unlike the Weyl lattice, the whole middle band in Fig. 3C is blue: the top surface has an extra band of soft modes. These modes have amplitudes that decay exponentially away from the top surface and can be alternatively interpreted as modes attached to the bottom surface and exponentially amplified into the bulk, similar to the modes studied in disordered systems in ref. 48.

Eq. 2 shows that topological polarization acts analogously to an electric polarization. Inside a homogeneous polarized material, the charge is zero. However, if homogeneity is broken by the presence of boundary or defects, charge can accumulate at these spots.

Local Rigidity and Softness at Dislocations

In topological mechanics, regions with positive softness charge contain mechanisms for soft elastic deformation, whereas regions with negative charge are prone to buckling (37). Eq. 2 suggests that regions of nonzero polarization flux can occur in a structurally inhomogeneous material without modulating the composition (i.e., as an alternative to composites). We choose dislocations to provide this inhomogeneity.

The natural defects in 3D crystals are line dislocations: displacements of unit cells along straight lines. In many crystalline solids, such defects control mechanical deformations and plasticity. These defect lines are topologically protected and must either terminate at the material boundaries or form a closed dislocation loop (Fig. 4 A and B). When the dislocations terminate at boundaries, they provide inhomogeneities extended along a channel (or curve) in 3D space. The topological invariant corresponding to a dislocation is the Burgers vector given by 𝒃=∮Cd𝒖, where 𝒖(𝒙) is the displacement of a particle at point x in the periodic lattice because of the presence of the dislocation and C is any closed contour surrounding the dislocation line. The invariant b is contour-independent and therefore, constant along the dislocation. Together with dislocation line direction ℓ^, the Burgers vector expresses how the dislocation affects the surrounding lattice.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Dislocation loops can act as dipoles of topological charge. (A) Architectures of a dislocation loop in a periodic lattice: in a screw dislocation, the Burgers vector b is parallel to the dislocation line, whereas in an edge dislocation, the Burgers vector and the dislocation line are perpendicular. A dislocation loop can combine edge and screw dislocations along its contour. In a polarized lattice (with polarization P), a dislocation line separates in space an edge dislocation segment that carries zero modes (ZMs; blue) from an edge dislocation segment that carries states of self-stress [SSS; red; for example, via screw dislocations (green), which carry no such charges]. The net topological charge, defined as the difference between ZM and SSS, is zero when summed over a dislocation loop contour. Nevertheless, the dislocation loop carries a topological charge dipole, which is, in this example, parallel to the Burgers vector b. (B) Geometry of the dislocation loop. Each prism represents a unit cell with triangles oriented according to (ϕ,θ,ψ)=(π/3, 0,π/3) as shown in Inset. (C) Numerical results for the four lowest-frequency phonons (excluding the trivial translations) in a large (36×36×18 unit cells) polarized lattice that has a dislocation loop: warmer color signifies a larger displacement within that unit cell. Plotted are only those unit cells that have displacements above a cutoff of 20% of the maximum. Note the localization of the softest modes to the near side of the loop in accordance with A and the polarization P. (D) The common logarithm of lowest mode frequencies log⁡ω2 for the Nth lowest mode plotted vs. N, comparing two large samples (same size as C). The mode frequencies are significantly lower for the dislocated lattice (black) than in the nondislocated case (gray). Without a dislocation, the lowest eigenmodes are the extended plane-wave acoustic phonons, whereas with the dislocation, these modes include both the acoustic phonons and the modes localized along the dislocation loop (in blue, four lowest modes with eigenvectors that are plotted in C).

Line dislocations come in two primary types called edge and screw dislocations (Fig. 4A). For edge dislocations, b and ℓ^ are orthogonal: the displacement u pushes unit cells apart to insert a half-plane of unit cells that extends from the dislocation line in the direction ℓ^×𝒃. In this way, edge dislocations are 3D generalizations of 2D point dislocations. By contrast, for screw dislocations, the b and ℓ^ vectors are parallel: the displacement u pushes neighboring cells apart along the loop. In this way, screw dislocations give rise to an inherently 3D spiral structure. Along a dislocation loop, the Burgers vector b is constant, but the line direction ℓ^ changes. As Fig. 4A shows, a dislocation loop can contain both edges and screws.

The interplay between line properties (ℓ^,𝒃) and topological polarization P of the lattice leads to a softness charge. On evaluating Eq. 2 in a deformed region (Rigidity Charge Bound to Line Dislocations), we obtain the total topological charge to be νS=Vcell−1∫Sd3𝒙Pi[∂i,∂j]uj(𝒙), where Vcell is the unit cell volume of the uniform lattice. Significantly, the charge density is nonzero when the partial derivatives fail to commute, which happens only at the dislocation line itself and only when {𝑷,𝒃,ℓ^} are linearly independent. Then, the volume integral for the topological charge reduces to an integral along the dislocation line:νS=1Vcell∫0T(𝑷×∫Ctd𝒖)ℓ^dt,[3]where t= 0,…,T parameterizes the dislocation line and Ct is a family of contours enclosing this line at each t. The integrand of Eq. 3 gives the charge line density ρL at each t, which using the definition of the Burgers vector, evaluates toρL=1Vcell(𝑷×𝒃)⋅ℓ^[4](Rigidity Charge Bound to Line Dislocations has a detailed derivation). Equivalently, we could write ρL=Vcell−1det(𝑷,𝒃,ℓ^). In particular, screw dislocations do not carry charge—their Burgers vector points along the line direction. For edge dislocations, the inserted rows terminate at the dislocation line, creating a boundary within the material. If P points along (against) the half-plane direction ℓ^×𝒃, these rows carry extra polarization out (in). Note that, for any dislocation loop, the net charge ∮LdtρL is zero. However, dislocations loops do separate softness charges in space, resulting in a topological charge dipole, as shown in Fig. 4A. This dipole is oriented within the 𝒃–P plane and can be quantified by a dipole moment, which captures the amount of charge separated and the distance of separation. The dipole’s b (P) component is the area of the dislocation loop after it is projected onto the plane formed by the two vectors b (P) and 𝒃×𝑷 (Dislocation Loop Dipole Moment shows the derivation). For the loop in Fig. 4A composed of edge and screw dislocations, all of the charges are localized along the edge dislocations, and the dipole moment is along b.

A positive charge density corresponds to soft modes localized along part of the dislocation loop. We investigate this localized softness within a polarized lattice using the configuration shown schematically in Fig. 4A and plotted for a small sample in Fig. 4B. In Fig. 4C, we show that, for the softest (i.e., lowest-frequency nontranslational) modes of this lattice, the unit cells with the largest displacements are localized along the near side of the loop, in agreement with Eq. 4. This can be contrasted with the lowest-energy modes of a sample without a dislocation: in Fig. 4D, we show that the dislocated lattice has soft modes at lower frequencies. Whereas the lowest modes in the nondislocated sample are the largest-wavelength acoustic phonons that fit within the periodic box, the lowest modes of the dislocated lattice are a combination of these acoustic phonons and the many localized modes, such as the ones shown in Fig. 4C.

Our work suggests a design principle for materials that will localize either large strain or large stress along a boundary, line, or parts of a closed loop, leading to preprogrammed yielding and failure.