Localizing softness and stress along loops in three-dimensional topological metamaterials
Guido Baardink Anton Souslov Jayson Paulose Vincenzo Vitelli
cond-mat.soft cond-mat.mtrl-sci





Instituut-Lorentz, Universiteit Leiden, 2300 RA Leiden, The Netherlands

Instituut-Lorentz, Universiteit Leiden, 2300 RA Leiden, The Netherlands

Instituut-Lorentz, Universiteit Leiden, 2300 RA Leiden, The Netherlands
Departments of Physics and Integrative Biology, University of California, Berkeley, CA 94720, USA 

Instituut-Lorentz, Universiteit Leiden, 2300 RA Leiden, The Netherlands




Topological states can be used to control the mechanical properties of a material
along an edge or around a localized defect. The surface rigidity of elastic
networks is characterized by a bulk topological invariant called the
polarization; materials with a well-defined uniform polarization display a
dramatic range of edge softnesses depending on the orientation of the polarization relative to the
terminating surface. 
However, in all three-dimensional mechanical metamaterials proposed to date, the topological edge modes are mixed with bulk soft modes and so-called Weyl loops. Here, we report the design of a gapped
3D topological metamaterial with a uniform polarization that displays a corresponding asymmetry
between the number of soft modes on opposing surfaces and, in addition, no bulk
soft modes. We then use this construction to localize topological soft modes in interior
regions of the material by including defect structures—dislocation
loops—that are unique to three dimensions. We derive a general formula that
relates the difference in the number of soft modes and states of self-stress localized along the dislocation loop to the handedness of the 
vector triad formed by the lattice polarization, Burgers vector, and dislocation-line direction. Our
findings suggest a novel strategy for pre-programming failure and softness localized
along lines in 3D, while avoiding extended periodic failure modes
associated with Weyl loops.




Localizing softness and stress along loops in three-dimensional topological metamaterials
    Vincenzo Vitelli
    December 30, 2023
=========================================================================================






§ INTRODUCTION

Mechanical metamaterials can control softness via 
a balance between the number of degrees of freedom of their components or nodes
and the number of constraints due to connections or links <cit.>.
This balance, first noted by Maxwell <cit.> and later explored by Calladine <cit.>, is termed isostaticity.
In isostatic materials, softness can manifest itself via large-scale deformations, for example
as uniform Guest-Hutchinson modes <cit.> or via periodic soft deformations, corresponding to so-called Weyl modes <cit.>.
Uniform softness can be exploited to create extraordinary mechanical response <cit.>, such as materials with a negative Poisson's ratio <cit.>.
Alternatively, localized softness has been programmed into isostatic materials in one and two dimensions
via a topological invariant called the polarization <cit.> that controls mechanical response and stress localization <cit.> at an edge [including the edge of a disordered sample <cit.>], an interface, or bound to a moving soliton <cit.>.
These mechanical <cit.> examples of topological metamaterials <cit.> exhibit a general feature of topological matter <cit.>: a correspondence between integer invariants in the bulk and response at a boundary.
Large-scale and localized deformations are deeply intertwined, as can be seen in demonstrations in which topological edge softness is created or destroyed by applying large uniform strains <cit.>.
A combination of topological polarization and localized defects can be used to program in softness or failure at a specified region in the material <cit.>.




Although there are a number of examples of isostatic periodic structures in one, two, and three dimensions, the three-dimensional case is unique because all prior realizations of three-dimensional isostatic lattices include large-scale periodic deformations along continuous lines in momentum space <cit.>.
These Weyl lines define families of periodic soft modes in the material bulk
and contain a number of modes, which scales with the linear size of the structure.
As Ref. <cit.> explores, Weyl lines can be useful to create a metamaterial surface with anisotropic elasticity, but in order to create a material whose top surface is much softer than the bottom, it proves necessary to collapse two Weyl lines on top of each other. 
An alternative would be to find a metamaterial without Weyl lines.
However, these Weyl lines are generic and have a topological character which ensures that they cannot be annihilated locally—a single Weyl loop can only be destroyed by shrinking it to a point.
This presents a challenge in three-dimensional isostatic metamaterial design: 
to achieve a gapped topological material analogous to those in two dimensions <cit.>, in which softness can be controlled and localized.


In this work, we design gapped topological materials by exploring the parameter space of the generalized stacked kagome lattice. 
The mechanical response of this gapped structure is characterized by a nonzero topological polarization oriented along the z-axis. 
This topological polarization P can be exploited 
to localize soft modes in the material bulk by introducing topological defects within the lattice structure called dislocation loops.
These dislocations are characterized by a topological invariant called the Burgers vector b.
Along the dislocation, we show that the topological charge characterizing the softness or rigidity of the lattice (with unit cell volume V_cell) depends on the orientation ℓ̂ of the dislocation line and is given by P·( b×ℓ̂)/V_cell per unit length.









§ 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. Within this set-up, softness corresponds to displacements with small strains, and rigidity to strains that induce only small displacements.
We place a particle (or ball) at each vertex and connect the neighboring particles by linear springs. Such models capture the small-strain response of realistic structures that are either 3D-printed from soft polymers <cit.> or assembled from construction sets <cit.> or laser-cut components <cit.>. 

The mechanics of this ball-spring model are captured via the linearized equation of motion 
Ẍ = - D X, where X≡(x_1,…,x_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 N_B-dimensional vector of spring extensions S≡(s_1,…,s_N_B) to the displacements via S=RX.
The rigidity matrix R contains dN× N_B 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 <cit.>. Although the rigidity theory in terms of the matrix R offers an equivalent description to the dynamical matrix, it contains additional physical insight.

The null spaces (kernels) of the rigidity matrix and its transpose give us information about lattice mechanics. From the relation S=RX, we note that R contains soft modes, i.e., collective displacements that (to lowest order) do not stretch or compress any of the springs. On the other hand, the matrix R^𝖳 relates forces on particles Ẍ to spring strains via Ẍ = - R^𝖳S. From this relation, one notes that R^𝖳 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 <cit.>. Localized states of self-stress can concentrate uniformly-applied external loads, which can lead to selective buckling in pre-programmed regions <cit.>.

To calculate soft modes (states of self-stress), we need to numerically solve the configuration-dependent equation RX = 0 (R^𝖳S = 0). 
However, a mathematical result called the rank-nullity theorem lets us compute the difference between the number of zero modes N_0 ≡null R and the number of states of self-stress N_B ≡null R^𝖳 <cit.>. This difference, the rigidity charge ν≡ N_0 - N_B, is given by the dimensionality of the R matrix: ν = dN - N_B,
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≡ 2N_B/N = 2d= 6 and ν = 0: the lattice is isostatic.
This lattice is based on the two-triangle unit cell shown in Fig. <ref>A, 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 fix the connectivity and explore a range for each of the three angles from -π/2 to π/2. In Fig. <ref>B, 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. <ref>C. We proceed to quantitatively demonstrate that this metamaterial is gapped.


We calculate the spectrum of normal modes by considering plane-wave solutions of the form x_α,ℓ=e^iℓ·kx_α, 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 wavevector k is periodic in each component with -π≤ k_x, k_y,k_z≤π. 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_αβ(k)=R_αβe^i(ℓ_α-ℓ_β)·k. A zero mode at wave-number k is a vector of unit-cell displacements x_α that solves R_βα(k)x_α=0. Thus, within this setting, zero modes correspond to the zeros of the complex function R(k). 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 (see Supporting Information [SI] for details) and for most values of the angles (ϕ,θ,ψ) find collections of zero modes along compact loops in k-space [Fig. <ref>A]. These Weyl loops, shown in Fig. <ref>B–C, 
are analogous to one-dimensional nodal lines <cit.> and experimentally observed zero-dimensional Weyl points <cit.> in three-dimensional electronic semimetals and photonic crystals.
In isostatic lattices, the number of Weyl loops is always even, because the materials' time-reversal symmetry maintains k→-k reflection symmetry in the Brillouin zone—each Weyl loop comes in a pair with its reflected partner <cit.>. 
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. <ref>A, 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. <ref>A,D].

We conclude that although Weyl loops are generic, the stacked kagome lattice also exhibits gapped configurations which contain no Weyl loops. 
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 <cit.>. 
We are not aware of any other realizations of a three-dimensional isostatic lattice which is gapped.
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 <cit.>. 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 apparent by looking at the spectrum of the material's soft surface waves. The polarization controls which surfaces have more soft modes: this is the mechanical version of the bulk-boundary correspondence principle.





For the gapped polarized lattice, a polarization vector P can be computed via P=∑_i m_i a_i - d_0, where a_i are the three lattice vectors 
and, in our case, the unit-cell dipole d_0 = -3  a_3/2. This dipole is computed via the expression d_0 ≡ - ∑_br_b, where r_b
are the positions of the bond centers relative to the center of mass of the unit cell <cit.>.
The three coefficients m_i are winding numbers computed by integrating the phase change of the complex function R( k) across a straight-line contour crossing the Brillouin zone:

    m_i ≡ - 1/2π i∮_C_id k_i d/d k_iln R( k),

where C_i is a closed contour in the k̂_i direction <cit.>. This integral is well defined for contours along which R( k) ≠ 0. The integers m_i depend smoothly on the choice of contour and are constant. In gapped lattices, the determinant is non-zero everywhere except for the origin k=0. As a consequence, the winding numbers are independent of the chosen contour and the polarization is a topological invariant. On the other hand, Weyl loops partition the space of straight-line contours into lines going through the inside (outside) of the loop. Because R( k)=0 along the loops, contours on either side of the loop can have different winding numbers m_i. 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. <ref>C. The winding number m_W around this small circle is an invariant and equal to the difference between m_i for C_i 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. <ref>A. 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 demonstrate 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) (c.f., Fig. <ref>A,C), the winding number around the loop is m_W=- 1. The m_3 for contours inside and outside of the loop differ by one. To see the boundary counter-part of the correspondence, we compute the spectrum of soft edge states for this lattice with a stress-free surface parallel to the xy-plane (see Fig. <ref>C). In Fig. <ref>A–B,
we show the (signed) inverse penetration depth κ for soft surface modes in the two-dimensional Brillouin zone for wavevector (k_x,k_y). These plane-wave solutions of R = 0 correspond to soft modes of the form e^iℓ·(k+iκ)x_α. The middle mode of Fig. <ref>A (see zoom in Fig. <ref>B) changes sides: for κ > 0 (red region), the soft mode is attached to the bottom surface and decays upwards into the bulk, whereas for κ < 0 (blue region), the mode is localized on the top surface. 
For the line corresponding to the projection of the Weyl loop onto the two-dimensional Brillouin zone, the penetration depth is infinite and κ = 0, because the Weyl lines are soft modes in the bulk.
The difference in m_3 between inside and outside of the Weyl loops is  the bulk invariant m_W = - 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. (<ref>), which states that the third winding numbers m_3 count, up to a constant, the number of modes with zero energy (R = 0) at a boundary. These observations confirm bulk-boundary correspondence.


A similar correspondence exists in the gapped lattice. There, the invariant polarization pointing along the z-axis is given by P =  a_3/2, which we computed for parameters (ϕ,θ,ψ)=(π/3,0,π/3). 
To understand the effect of P on the boundary, note in analogy with electromagnetism that the rigidity charge ν^S in region S is related to the flux of polarization P through the region's boundary ∂ S. In a nearly-uniform, gapped lattice <cit.>

    ν^S=∮_∂ SdA/V_celln̂· P,

where n̂ is the boundary's inward normal and V_cell=( a_1, a_2, a_3) is the volume of a unit cell. 
The difference between the rigidity charges at the top versus 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 <cit.>.
We plot the inverse penetration depth κ for this extra band in Fig. <ref>C–D within the two-dimensional Brillouin zone. 
Unlike the Weyl lattice, the whole middle band in Fig. <ref>C is blue: the top surface has more soft modes and is, therefore, softer.

Equation (<ref>) 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. 



figure-1





§ LOCAL RIGIDITY AND SOFTNESS AT DISLOCATIONS

For many applications such as cushioning <cit.>, programmed assembly <cit.>, or controlled failure <cit.>, materials need to be designed with build-in softness or rigidity. The previous section illustrates how a topological invariant can be programmed into a material to design softness at boundaries. We now proceed to show how a flux of this topological polarization can be harnessed to create soft regions in the bulk of the material. Equation (<ref>) suggests that a nonzero polarization flux can be achieved by considering an inhomogeneous material.
We choose dislocations to provide this inhomogeneity.

The natural defects in three-dimensional crystals are line dislocations: displacements of unit cells along straight lines which terminate at material boundaries. In many crystalline solids, such defects control mechanical deformations and plasticity. As these lines are topologically protected, the only way to confine a line dislocation to a finite region is to form a closed dislocation loop (Fig. <ref>A–B).  The topological invariant characterizing a dislocation is the Burgers vector b, which expresses the effect of a dislocation on the surrounding lattice. If a particle at point x is displaced to x+ u( x) according to displacement u( x), then the Burgers vector is b=∮_C d u, where C is any closed contour surrounding the dislocation line.

Line dislocations come in two primary types called edge and screw dislocations, which are distinguished by the orientation of the Burgers vector b relative to the dislocation line direction ℓ̂ (see Fig. <ref>A). For edge dislocations, these vectors are orthogonal: the displacement u pushes unit cells apart in order to insert a half-plane of unit cells that extends from the dislocation line in the direction ℓ̂× b. In this way, edge dislocations are three-dimensional generalizations of two-dimensional 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 three-dimensional spiral structure. Along a dislocation loop, the Burgers vector b is constant, but the line direction ℓ̂ changes. As Fig. <ref>A shows, a dislocation loop can contain both edges and screws.

The interplay between line properties (ℓ̂,b) and the topological polarization P of the lattice leads to a rigidity charge that can be expressed in a simple formula (derived in the SI) for the charge line density ρ^L localized at the dislocation line: 

    ρ^L = 1/V_cell P·( b×ℓ̂),

where  V_cell is the unit cell volume of the uniform lattice.
Equivalently, we could write ρ^L = V_cell^-1 ( P, b,ℓ̂). Consequently, for any configuration in which these three vectors are not linearly independent, the charge density is zero. In particular, screw dislocations do not carry charge—their Burgers vector points along the line direction. For edge dislocations, the charge is negative (positive) if P points along (against) the half-plane direction ℓ̂× b: intuitively, these rows carry extra polarization out or in.
Note that for any dislocation loop, the net charge ∮_L ρ^L is zero. 
However, dislocations loops do separate rigidity charges in space, resulting in a topological charge dipole, as shown in Fig. <ref>A. 
This dipole is oriented within the b-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 b×P (see SI for derivation). For the loop composed of edge and screw dislocations
shown in Fig. <ref>A, all of the charges are localized along the edge dislocations and the dipole moment lies 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. <ref>A and plotted for a small sample in Fig. <ref>B. In Fig. <ref>C, we demonstrate that for the softest (i.e., lowest-frequency non-translational) modes of this lattice, the unit cells with the largest displacements are localized along the near side of the loop, in agreement with Eq. (<ref>).
This can be contrasted with the lowest-energy modes of a sample without a dislocation: in Fig. <ref>D, we show that the dislocated lattice has soft modes at lower frequencies.
Whereas the lowest modes in the non-dislocated 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. <ref>C.




§ CONCLUSIONS AND OUTLOOK

We demonstrated that dislocation loops in polarized three-dimensional lattices
can localize soft modes at one part of the loop. Similarly, the opposite part of the loop localizes states of self-stress. Although the total topological charge along the loop is zero, the loop acts as a topological charge dipole. The loop can be used to localize topologically protected soft modes in one region of the material, inside the material bulk. Furthermore, because the lattice we have designed is gapped, these localized soft modes do not compete with the many extended, periodic soft modes characteristic of Weyl loops.

Let us conclude by discussing applications for both three-dimensional gapped lattices and for localized softness. Three-dimensional topological materials may be useful in the design of devices for cushioning, including safety devices such as helmets, in which a hard side necessary to resist stress exists in combination with a soft side necessary to cushion the body <cit.>.
Localized softness can be used to pre-program large displacements and to isolate the rest of the material from strain. In active mechanical metamaterials <cit.>, for example in self-folding origami <cit.>, these mechanisms could be actuated using motors to strain the material in a well-controlled way <cit.>. Furthermore, localized states of self-stress can be used to control material failure via either buckling or fracture, which then isolates the rest of the material from failure even under a large external load. 




§ ACKNOWLEDGMENTS

We thank Paul Baireuther and Bryan G. Chen for fruitful discussions and T. C. Lubensky for a critical reading of the manuscript. 
We gratefully acknowledge funding from FOM, NWO, and Delta Institute for theoretical physics.



seccntformat#1
  c@#1@̧sectionthe#1  









§ SUPPORTING INFORMATION FOR “LOCALIZING SOFTNESS AND STRESS ALONG LOOPS IN THREE-DIMENSIONAL TOPOLOGICAL METAMATERIALS”











§ GENERALIZED STACKED KAGOME LATTICE


The generalized stacked kagome lattice is built by repeating a given unit cell along the following lattice vectors

    [        a_1=a(1,0,0) a_2=a(1/2,√(3)/2,0)        a_3=a(0,0,1) ],

where a defines the lattice spacing. In the undeformed kagome lattice the equilibrium positions of the particles in the unit cell is given by the triangular prism

    [        r(1)=-14 a_3  r(2)=12 a_1-14 a_3  r(3)=12 a_2-14 a_3;        r(4)=+14 a_3  r(5)=12 a_1+14 a_3 r(6)=12 a_2+14 a_3, ]

where particles s=1,2,3 (s=4,5,6) constitute the bottom (top) triangle of the prism.

Within the unit cell we connect the particle pairs (2,3), (3,1), (1,2), (5,6), (6,4), (4,5), (1,4), (2,5), (3,6) by springs. Horizontally between unit cells we connect pairs (3,2), (1,3), (2,1), (6,5), (4,6), (5,4) according to the nearest neighbor principle. Finally we connect (4,1), (5,2) and (6,3) between the unit cells in the positive vertical direction.

We let a deformation act on the unit cell by rotating the triangles around three orthogonal axes attached to the triangles. The first two axes are coplanar to the triangle, with the first one bisecting one of its angles. The third axis is perpendicular to the triangle. The three axes meet in the center of the triangle. We then define our deformation by rotating the top triangle around the three axis by angles ϕ, θ, and ψ respectively, while rotating the bottom triangle by -ϕ, -θ, and -ψ. Note that by taking our rotation axes relative to the triangles, our rotation group is commutative and we do not need to worry about the order of rotations. Explicitly, the equilibrium positions of the particles in the generalized stacked kagome can be given as


    r(1)   =a/2√(3)[    s_ψ c_θ-√(3)c_ψ c_ϕ-c_ψ s_ϕ s_θ; -c_ψ c_θ-√(3) s_ψ c_ϕ- s_ψ s_ϕ s_θ;           - c_ϕ s_θ+ √(3)s_ϕ+ √(3) ]+ar/2[ c_α; s_α;   0 ]
     r(2)   =a/2√(3)[    s_ψ c_θ+√(3)c_ψ c_ϕ-c_ψ s_ϕ s_θ; -c_ψ c_θ+ √(3)s_ψ c_ϕ- s_ϕ s_ϕ s_θ;           - c_ϕ s_θ- √(3)s_ϕ+ √(3) ]+ar/2[ c_α; s_α;   0 ]
     r(3)   =a/√(3)[ - s_ψ c_θ+c_ψ s_ϕ s_θ;  c_ψ c_θ+ s_ψ s_ϕ s_θ;      c_ϕ s_θ + √(3)/2 ]+ar/2[ c_α; s_α;   0 ]
     r(4)   =a/2√(3)[  - s_ψ c_θ-√(3)c_ψ c_ϕ-c_ψ s_ϕ s_θ; -c_ψ c_θ+√(3) s_ψ c_ϕ+ s_ψ s_ϕ s_θ;             c_ϕ s_θ- √(3)s_ϕ- √(3) ]-ar/2[ c_α; s_α;   0 ]
     r(5)   =a/2√(3)[  - s_ψ c_θ+√(3)c_ψ c_ϕ-c_ψ s_ϕ s_θ; -c_ψ c_θ- √(3)s_ψ c_ϕ+ s_ϕ s_ϕ s_θ;             c_ϕ s_θ+ √(3)s_ϕ- √(3) ]-ar/2[ c_α; s_α;   0 ]
     r(6)   =a/√(3)[  s_ψ c_θ+c_ψ s_ϕ s_θ; c_ψ c_θ- s_ψ s_ϕ s_θ;   - c_ϕ s_θ - √(3)/2 ]-ar/2[ c_α; s_α;   0 ],


where s_·=sin(·) and c_·=cos(·).


Following the above order of bond enumeration, the compatibility matrix takes the form:

    R^†( k)=[                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   [  b̂^𝖳_1;  b̂^𝖳_2;  b̂^𝖳_3;  b̂^𝖳_4;  b̂^𝖳_5;  b̂^𝖳_6;  b̂^𝖳_7;  b̂^𝖳_8;  b̂^𝖳_9; b̂^𝖳_10; b̂^𝖳_11; b̂^𝖳_12; b̂^𝖳_13; b̂^𝖳_14; b̂^𝖳_15; b̂^𝖳_16; b̂^𝖳_17; b̂^𝖳_18 ] [              0             -1              1         b̂^𝖳_1;              1              0             -1         b̂^𝖳_1              O;             -1              1              0         b̂^𝖳_1;         b̂^𝖳_1              0             -1              1;         b̂^𝖳_1              O              1              0             -1;         b̂^𝖳_1             -1              1              0;         b̂^𝖳_1;         b̂^𝖳_1             -I              I;         b̂^𝖳_1;              0       e^-i k_1             -1         b̂^𝖳_1;             -1              0       e^-i k_2         b̂^𝖳_1              O; e^-i (k_1-k_2)             -1              0         b̂^𝖳_1;         b̂^𝖳_1              0       e^-i k_1             -1;         b̂^𝖳_1              O             -1              0       e^-i k_2;         b̂^𝖳_1  e^-i(k_1-k_2)             -1              0;         b̂^𝖳_1;         b̂^𝖳_1      e^-i k_3I             -I;         b̂^𝖳_1;                ] ]


where the initial spring direction b_i is defined as the normalization of r(s')- r(s), where (s,s') is the i-th spring. In this notation, each number in the matrix grid represents a multiple of the triple on the far left in the same row. This turns the matrix into an 18×18 square matrix. As the initial bond directions depend smoothly on the initial particle positions, each element of Q depends smoothly on the deformation parameters ϕ, θ, and ψ.




§ SOFT DIRECTIONS


Weyl lines are formed by collections of points satisfying R( k)=0. This equation is often quite hard to solve. However, with observational evidence suggesting that the Weyl loops generally connect to the origin, we can formally expand the determinant around k≡|| k||=0 as

    R( k)
    =∑_n∈ℤ^3c_ ne^-i n· k
    =i k^3𝒫_3(k̂)+k^4𝒫_4(k̂)+𝒪(k^5),

for finitely many nonzero coefficients c_ n∈ℝ, where 𝒫_m is an m-th order homogeneous polynomial with real coefficients in three variables over the unit sphere ||k̂||=1. The vanishing constant, linear, and quadratic polynomials reflect the triple multiplicity of the zero at k=0.

A Weyl line can only attach to the origin at a direction k̂ where both 𝒫_3(k̂) and 𝒫_4(k̂) vanish. 
Suppose, for instance, that 𝒫_4(k̂)  0 and note that the error term is dominated by M k^5 for all k<δ, where M and δ are fixed positive numbers.
As a consequence, the relation

    | R( k)|> k^4(|𝒫_4(k̂)|-M k)

holds for all k<δ. In particular, this means | R( k)|>0 for all k<min(δ,1/M|𝒫_4(k̂)|). We call these common zeroes of 𝒫_3(k̂) and 𝒫_4(k̂) over the unit sphere the soft directions of R. 

As any Weyl loop connects to the origin along two directions, we find that a configuration with 2n soft directions can accommodate no more than n pairs of Weyl loops. 
In the main text we show the θ=0 slice of soft-mode phase space, containing a region where the lattice does not exhibit any soft directions. This absence of soft directions is rare. Sampling of the 0<ϕ,θ,ψ<π/2 phase space with step size π/60 we obtain the following data for the prevalence of different numbers of soft directions:






§ RIGIDITY CHARGE BOUND TO LINE DISLOCATIONS


The dislocation acts on the polarization and the unit cell volume through the lattice vectors a_i. Letting the displacement function u act on both the initial and terminal point of the vector a_i, we obtain the lattice vector field in the dislocated lattice:

    ã_i( x)
    =[ x+ a_i+ u( x+ a_i)]-[ x+ u( x)]
    = a_i+( a_i·∇) u( x),

where we assume that higher order differentials of u are negligible.
The polarization changes as

    P̃( x)
    =∑_i m_iã_i
    = P+( P·∇) u( x).

The unit cell volume changes as

    Ṽ_cell( x)
    =( ã_1,ã_2,ã_3)
    =V_cell(1+∇ u( x)).

We start from Eq. (1) in the main text. In a canonical basis with coordinates (x^1,x^2,x^3), a surface element is given by dx^ℓ∧ dx^k=ϵ^iℓ kdS_i, where d S is positively oriented. In this notation, the change in rigidity charge in the region S from the ν^S=0 periodic case is given by

    ν̃^S-ν^S
       =∫_∂ S(Ṽ_cell^-1 P̃-V_cell^-1  P)·(-d S)
       
    =V_cell^-1∫_∂ S P_i∂_ju_j-P_j∂_ju_i dS_i
       
    =V_cell^-1∫_∂ SP_i ϵ_ijk∂_ℓ u_j ϵ^mℓ kdS_m
       
    =V_cell^-1∫_∂ S( P× d u)_k∧ dx^k.


Consider a straight dislocation line segment L in direction ℓ̂ and let S be the cylinder with axis L and radius ℛ. Without loss of generality, we choose canonical coordinates such that L={tℓ̂: 0<t<T} with ℓ̂ pointing along the third basis vector. Furthermore, we take Tℓ̂ to be an integer combination of lattice vectors a_i. If the displacement function respects the periodicity of the lattice, the top (t=T) and bottom (t=0) surfaces of S will cancel by orientation. It remains to calculate Eq. (<ref>) along the vertical boundary of S. First consider the k=1 term:

    ∫_∂ S( P× d u)_1∧ dx^1=∫_0^2π( P×∫_L_θd u)_1 dℛcosθ=0,

where ∂ S=⋃_θ L_θ is a decomposition of the cylindrical boundary in vertical lines at different angle θ. The last equality holds since ∫_L_θd u is independent of θ as can be shown by contour integration

    ∫_L_θd u-∫_L_θ'd u=∫_θ^θ'd u|^t=T-d u|^t=0=0.

The k=2 term vanishes similarly. Hence we find that the rigidity charge is given entirely by the term corresponding to the line direction. Let ∂ S=⋃_t C_t be the decomposition of the cylindrical boundary in horizontal circles at different height t. Recalling the definition of the Burgers vector, we then obtain

    ν̃^S
    =V_cell^-1∫_0^T( P×∫_C_td u)_3  dt=1/V_cell( P× b)· Tℓ̂.

The rigidity charge density along the line can then be given by ρ_L≡ν_L/T=V_cell^-1( P, b,ℓ̂), where the three vectors P, b and ℓ̂ form the columns of a 3×3 matrix.



§ DISLOCATION LOOP DIPOLE MOMENT


Let r(s) be an arc length parametrization of the loop L. The total charge along the loop is given by

    ν^L=∮_Lρ^L( r)ds∝∮_L r'(s)ds=0.

As the loop separates charges of different sign, it is sensible to consider the associated dipole moment

    d^L≡∮_L  rρ^L( r)ds=∮_L  r(s)( P× b)· r'(s)/V_cellds.

Without loss of generality we can choose our basis (ê_1,ê_2,ê_3) such that P× b points in the direction of ê_3. We then obtain:

    d^L=|| P× b||/V_cell∮_L r  dr_3
    =|| P× b||/V_cell(𝒜_13,𝒜_23, 0),

where 𝒜_ij is the signed area of the projection of the dislocation loop on the plane generated by the basis (ê_i,ê_j).




10

Bertoldi2010
Bertoldi K, Reis PM, Willshaw S, Mullin T (2010) Negative poisson's ratio
  behavior induced by an elastic instability.
Advanced Materials 22(3):361–366.

Sun2012
Sun K, Souslov A, Mao X, Lubensky TC (2012) Surface phonons, elastic response,
  and conformal invariance in twisted kagome lattices.
Proceedings of the National Academy of Sciences
  109(31):12369–12374.

Kane2014
Kane CL, Lubensky TC (2013) Topological boundary modes in isostatic lattices.
Nature Physics 10(1):39–45.

Chen2014
Chen BGg, Upadhyaya N, Vitelli V (2014) Nonlinear conduction via solitons in a
  topological mechanical insulator.
Proceedings of the National Academy of Sciences
  111(36):13004–13009.

Lubensky2015
Lubensky TC, Kane CL, Mao X, Souslov A, Sun K (2015) Phonons and elasticity in
  critically coordinated lattices.
Reports on Progress in Physics 78(7):073901.

Meeussen2017
Meeussen AS, Paulose J, Vitelli V (2016) Geared topological metamaterials with
  tunable mechanical stability.
Phys. Rev. X 6(4):041029.

Rocks2017
Rocks JW, et al. (2017) Designing allostery-inspired response in mechanical
  networks.
Proceedings of the National Academy of Sciences
  114(10):2520–2525.

Yan2017
Yan L, Ravasio R, Brito C, Wyart M (2017) Architecture and coevolution of
  allosteric materials.
Proceedings of the National Academy of Sciences
  114(10):2526–2531.

Coulais2017
Coulais C, Sounas D, Alù A (2017) Static non-reciprocity in mechanical
  metamaterials.
Nature 542(7642):461–464.

Maxwell1864
Maxwell JC (1864) On the calculation of the equilibrium and stiffness of
  frames.
Philosophical Magazine 27(182):294–299.

Calladine1978
Calladine C (1978) Buckminster fuller's “tensegrity” structures and clerk
  maxwell's rules for the construction of stiff frames.
International Journal of Solids and Structures 14(2):161 –
  172.

Guest2003
Guest S, Hutchinson J (2003) On the determinacy of repetitive structures.
Journal of the Mechanics and Physics of Solids 51(3):383 –
  391.

Kapko2009
Kapko V, Treacy MMJ, Thorpe MF, Guest SD (2009) On the collapse of locally
  isostatic networks.
Proceedings of the Royal Society of London A: Mathematical,
  Physical and Engineering Sciences 465(2111):3517–3530.

Rocklin2016
Rocklin DZ, Chen BGg, Falk M, Vitelli V, Lubensky TC (2016) Mechanical Weyl
  Modes in Topological Maxwell Lattices.
Physical Review Letters 116(13):135503.

Stenull2016
Stenull O, Kane C, Lubensky T (2016) Topological phonons and weyl lines in
  three dimensions.
Physical Review Letters 117(6):068001.

Florijn2014
Florijn B, Coulais C, van Hecke M (2014) Programmable mechanical metamaterials.
Phys. Rev. Lett. 113(17):175503.

Lakes1987
Lakes R (1987) Foam structures with a negative poisson’s ratio.
Science 235(4792):1038–1040.

Rocklin2016a
Rocklin D (2016) Directional mechanical response in the bulk of topological
  metamaterials.
arXiv preprint arXiv:1612.00084.

Sussman2016
Sussman DM, Stenull O, Lubensky TC (2016) Topological boundary modes in jammed
  matter.
Soft Matter 12(28):6079–6087.

Nash2015
Nash LM, et al. (2015) Topological mechanics of gyroscopic metamaterials.
Proc. Natl. Acad. Sci. USA 112(47):14495–500.

Susstrunk2015
Susstrunk R, Huber SD (2015) Observation of phononic helical edge states in a
  mechanical topological insulator.
Science 349(6243):47–50.

Kariyado2015
Kariyado T, Hatsugai Y (2015) Manipulation of dirac cones in mechanical
  graphene.
Scientific reports 5:18107.

Mousavi2015
Mousavi SH, Khanikaev AB, Wang Z (2015) Topologically protected elastic waves
  in phononic metamaterials.
Nature communications 6:8682.

Khanikaev2015
Khanikaev AB, Fleury R, Mousavi SH, Alù A (2015) Topologically robust
  sound propagation in an angular-momentum-biased graphene-like resonator
  lattice.
Nature Communications 6(May):8260.

Yang2015
Yang Z, et al. (2015) Topological Acoustics.
Physical Review Letters 114(11):114301.

Huber2016
Huber SD (2016) Topological mechanics.
Nature Physics 12(7):621–623.

Rechtsman2013
Rechtsman MC, et al. (2013) Photonic Floquet topological insulators.
Nature 496(7444):196–200.

Lu2014
Lu L, Joannopoulos JD, Soljačić M (2014) Topological photonics.
Nature Photonics 8(11):821–829.

Wang2015
Wang P, Lu L, Bertoldi K (2015) Topological Phononic Crystals with One-Way
  Elastic Edge Waves.
Physical review letters 115(10):104302.

Hasan2010
Hasan MZ, Kane CL (2010) Colloquium: Topological insulators.
Reviews of Modern Physics 82(4):3045–3067.

Bernevig2013
Bernevig BA, Hughes TL (2013) Topological Insulators and Topological
  Superconductors.

Rocklin2017
Rocklin DZ, Zhou S, Sun K, Mao X (2017) Transformable topological mechanical
  metamaterials.
Nature communications 8:14201.

Paulose2015
Paulose J, Chen BGg, Vitelli V (2015) Topological modes bound to dislocations
  in mechanical metamaterials.
Nature Physics 11:153–156.

Paulose2015a
Paulose J, Meeussen AS, Vitelli V (2015) Selective buckling via states of
  self-stress in topological metamaterials.
Proceedings of the National Academy of Sciences
  112(25):7639–7644.

Bilal2017
Bilal OR, Susstrunk R, Daraio C, Huber SD (2017) Intrinsically polar elastic
  metamaterials.
Advanced Materials p. 1700540.

Volovik2007
Volovik GE (2007) Quantum Phase Transitions from Topology in Momentum
  Space, eds. Unruh WG, Schützhold R.
(Springer Berlin Heidelberg, Berlin, Heidelberg), pp. 31–73.

Burkov2011
Burkov AA, Hook MD, Balents L (2011) Topological nodal semimetals.
Phys. Rev. B 84(23):235126.

Lu2013
Lu L, Fu L, Joannopoulos JD, Soljačić M (2013) Weyl points and line
  nodes in gyroid photonic crystals.
Nature photonics 7(4):294–299.

Xu2015
Xu SY, et al. (2015) Discovery of a weyl fermion semimetal and topological
  fermi arcs.
Science 349(6248):613–617.

Lu2015
Lu L, et al. (2015) Experimental observation of weyl points.
Science 349(6248):622–624.

Hawkes2010
Hawkes E, et al. (2010) Programmable matter by folding.
Proceedings of the National Academy of Sciences
  107(28):12441–12445.

Souslov2017
Souslov A, van Zuiden BC, Bartolo D, Vitelli V (2017) Topological sound in
  active-liquid metamaterials.
Nature Physics.

Silverberg2014
Silverberg JL, et al. (2014) Using origami design principles to fold
  reprogrammable mechanical metamaterials.
Science 345(6197):647–650.

Chen2016
Chen BGG, et al. (2016) Topological Mechanics of Origami and Kirigami.
Physical Review Letters 116(13):135501.




