Fracturing of topological Maxwell lattices

We present fracturing analysis of topological Maxwell lattices when they are stretched by applied stress. Maxwell lattices are mechanical structures containing equal numbers of degrees of freedom and constraints in the bulk and are thus on the verge of mechanical instability. Recent progress in topological mechanics led to the discovery of topologically protected floppy modes and states of self stress at edges and domain walls of Maxwell lattices. When normal brittle materials are being stretched, stress focuses on crack tips, leading to catastrophic failure. In contrast, we find that when topological Maxwell lattices are being stretched, stress focuses on states of self stress domain walls instead, and bond-breaking events start at these domain walls, even in presence of cracks. Remarkably, we find that the stress-focusing feature of the self-stress domain walls persists deep into the the failure process, when a lot of damages already occurred at these domain walls. We explain the results using topological mechanics theory and discuss the potential use of these topological Maxwell lattice structures as mechanical metamaterials that exhibit high strength against fracturing and well controlled fracturing process.


Introduction: topological mechanics and Maxwell lattices
In recent years, there has been substantial advances in applying the conceptual framework of topological states of matter to classical mechanical systems which are governed by Newton's laws . These mechanical systems can acquire exotic mechanical behaviors, such as one-way wave transport [5,6,[10][11][12][13][14][15][16][17], nonlinear soliton [4], switchable stiffness [21], and selective buckling [22], that originate in the topological states of their phonon band structures. Many of these systems belong to the class of Maxwell lattices, lattices that contain equal number of degrees of freedom and constraints in the bulk and hence are on the verge of mechanical instability. For any ideal mechanical system that consists mass points connected by harmonic springs on a periodic lattice in d-dimensional space, the condition for a Maxwell lattice is z = 2d, where z is the mean coordination number [3,[23][24][25][26][27][28][29][30][31]. This condition comes from balancing the degrees of freedom per site, d, with the average number of constraints per site, z /2.
For general mechanical structures with point-like particles and central-force bonds one can apply the Maxwell-Calladine counting rule [1,3,23,32], where N s is the number of sites in the lattice, N b is the number of springs, N 0 is the counting of floppy modes (modes of zero elastic energy) and N SS is the counting of states of self stress (eigenstates of tension and compression on bonds with no net force on any site). For a periodic Maxwell lattice which has no boundaries, we always have N s d − N b = 0 and therefore N 0 = N SS . A finite Maxwell lattice has fewer springs on the lattice boundaries and therefore N s d − N b > 0 and N 0 > N SS so there must be floppy modes. While the Maxwell-Calladine counting rule give us the relation between N 0 and N SS , it doesn't tell us where the floppy modes and states of self stress are in space, and whether they are localized or extended in the Maxwell lattice, which depend on the actual lattice architecture. A topological invariant, R T , called "topological polarization", was introduced by Kane and Lubensky in Ref. [1], to characterize localization of floppy modes and states of self stress. We will discuss more details of the topological polarization in the Theory section below. Similar to the topological states of matter in other systems, mechanical behaviors that are rooted in this topological invariant remain robust even when the system is locally perturbed. This makes topological mechanics a promising approach for designing mechanical metamaterials that are insensitive to impurities and defects.
In this paper, we study how Maxwell lattices fracture when they are macroscopically deformed by a uniaxial strain applied to the boundaries. We are especially interested in the fracturing process of lattices that are topologically polarized, i.e. R T = 0, and contain topologically protected localized states of self stress. It was shown in Ref. [22] that when Maxwell lattices with self-stress domain walls buckle under pressure, the buckling regions localize to the self-stress domain walls because stress focuses on these domain walls. Similarly, we find that fracturing of these Maxwell lattices under external tension starts at these self-stress domain walls due to this stress-focusing effect. More interestingly, even after a significant part of the self-stress domain walls is damaged during the fracturing process, stress still robustly focuses to these domain walls. After the first bond-breaking event, the lattice mean coordination number decreases below the Maxwell point, z < 2d, and thus simple topological mechanics theory that predicts localized states of self stress no longer directly apply. The fact that stress still focuses on the self-stress domain walls when the lattice is damaged originates from the robustness of these topologically protected edge states of these Maxwell lattices.
The stress and damage focusing effect of self-stress domain walls in Maxwell lattices provides remarkable protection on the rest of the structure. As we show below, these self-stress domain walls guide stress away from preexisting cracks in the structure. Thus, unlike normal materials where stress focuses on crack tips leading to catastrophic failure [33], Maxwell lattices with self-stress domain walls exhibit a much more gradual fracturing process where only a small number of bonds break at each step of strain increase, without any large avalanches. These observations reveal the great potential of using Maxwell lattices to design high-strength metamaterials in which failure occurs gradually at predicted locations.
To study the fracturing of Maxwell lattices, we choose two types of two-dimensional (2D) Maxwell lattices, the deformed square lattice and deformed kagome lattice, which are commonly used in the research of topological mechanics [1,20,23]. Both types of these Maxwell lattices can acquire nonzero R T by varying site positions in the unit cell. We first generate "phase diagrams" of R T in the parameter space of site positions and then determine the site positions that we need to polarize the lattice with desired R T . In fig. 1, we show these lattices with their topological polarization R T . The configurations of their unit cells will be used throughout this paper. To acquire states of self stress that are protected by topology, we introduce multiple domains of these lattices with different R T , which are separated by domain walls of floppy modes and states of self stress. Such structures are shown in fig. 2 for both the deformed square and the deformed kagome lattices. The connection between R T , domain walls and protected states of self stress will be explained in Sec. 2.

Theory: topological polarization and exponential decay of topological surface modes
In Maxwell lattices, the coordination number z = 2d ensures that the counting of constraints and degrees of freedom always equal locally, except near open boundaries of the lattices, where there are fewer constraints than degrees of freedom, giving rise to floppy modes. Whether an edge of a Maxwell lattice hosts exponentially localized floppy modes is characterized by the topological polarization R T , which is defined through calculating winding numbers of the determinate of the compatibility matrix C(q) over closed paths C ( fig. 1) across the first Brillouin zone, i.e. [1], The compatibility matrix C(q), as we discuss below in this Section, is a linear operator that maps lattice site displacements to bond extensions. We follow two paths along the reciprocal lattice vectors q 1 and q 2 and calculate the topological winding numbers n 1 and n 2 for any 2D periodic lattice. The topological polarization for a periodic Maxwell lattice is then defined as R T = n 1 r 1 + n 2 r 2 , where {r i } are the real space primitive vectors that correspond to the reciprocal vectors {q i } in momentum space.
To determine which boundaries have exponentially localized floppy modes or states of self stress, we need to calculate the surface integral of R T over a surface (contour for 2D lattices) S that encloses the lattice boundary of interest,.i.e. [1], where V 0 is the volume of the unit cell, G is the reciprocal vector orthogonal to the surface S. With symmetric choice of unit cells where the number of cut bonds are the In each of the four lattices, unit cells are marked by the black dotted boundary, and bonds that belong to the marked unit cell are marked as thick blue solid lines (other bonds in the lattice are thin black solid lines). Each unit cell in the deformed square lattice contains four sites, and each unit cell in the deformed kagome lattice contains three sites. The integral paths for determining R T are also marked in the first Brillouin zone (insets on the right) of the lattices. same on opposite edges (this is characterized by the local count ν S L as defined in Ref. [1] and symmetric unit cell means ν S L = 0), ν S T directly gives the difference between the numbers of floppy modes and states of self stress inside the surface S. In this paper we construct the deformed square and deformed kagome lattices with such symmetric unit cells. As a result, T is the counting of how many more localized floppy modes are inside S than localized states of self stress. Since both N 0 and N SS are non-negative, when ν S > 0, there are guaranteed floppy modes inside S that are protected by the topology of the bulk phonon structure of the lattice; When ν S < 0, there are guaranteed states of self stress inside S.
The formulation described above applies equally to lattice surfaces (outer boundaries, with vacuum considered R T = 0) and domain walls (inner boundaries between lattice domains with different R T ). Examples of lattices with domain walls are shown in fig. 2. When ν S > 0(< 0) for a surface that encloses a domain wall, there will be topologically protected floppy modes (states of self stress) localized on the domain wall. In these lattices, because the two domains separated by the domain wall are of homogeneous (but different) R T 's, the value of ν S is proportional to the height of the domain wall. Floppy modes and states of self stress that are localized on boundaries exponentially decay into the lattice bulk because their wave vectors q have nonzero imaginary parts, i.e. q = k +iκ. For example, a floppy modes on a lattice boundary u = A exp (iq · r) = A exp (ik · r − κ · r) decays into the lattice bulk with a decay rate |κ|. The sign of κ is dictated by the topological polarization R T (which is why R T determines which edge or domain wall floppy modes and states of self stress localizes to), whereas the magnitude of κ depends on the detailed geometry of the unit cell. To ensure localization of floppy modes or states of self stress on edges or domain walls, the lattice domain depth has to be greater than 1/|κ|. To determine the value of κ, we need to solve q from the condition for all floppy modes using the compatibility matrix, Taking the component of q along the boundary to be any real wave number, we can solve for κ that controls the decay rate perpendicular to the boundary. The decay rate of states of self stress can be similarly solved by requiring where Q = C T is the equilibrium matrix of the lattice. Since det Q = det C, the decay rate of states of self stress is the same as the decay rate of the floppy modes in a given domain. Finally, we need to connect the states of self stress to the linear response of the lattice to macroscopic deformations. We start with the equilibrium matrix Q and compatibility matrix C of a lattice [32,34] f = Qt, where t is the force on each springs, u is the displacement of each lattices sites, f is the net force on each lattice sites and e is the extension on each springs. It is straightforward to show that Q = C T . In addition, t is related to e as t = Ke where K is the diagonal matrix of the spring stiffness. For the study of lattice fracturing under external stress, we make explicit distinction between lattice boundaries with controlled displacement (denoted as ∂V in subscript) and the lattice bulk (V ) so that we can impose deformations on the boundaries and equilibrate the lattice in the bulk. The linear equations then become where boundary displacements u ∂V are given. We can solve t by setting f V = 0 (force balance on all internal sites) and simplify the above equations as Therefore t is a null vector of Q V and thus must be a superposition of states of self stress in the lattice bulk, i.e. t = i a i s i , where {s i } is the complete orthonormal basis for the null space of Q V . We can further express t as which shows that when bond extension caused by the boundary deformation u ∂V has nonzero overlap with states of self stress in the bulk, bonds in the lattice will have nonzero tension originated from these bulk states of self stress. We show the derivation of eq. (11) in detail in the Appendix.
Up to now we haven't consider the possibility that the lattices may have topologically protected Weyl points that can complicate our definition for R T . It has been shown that a deformed square lattice with four lattice sites in a unit cell can have zero or two Weyl points in the first Brillouin zone, depending on geometry of the lattice [20]. For a kagome lattice with three sites in the unit cell, there are no Weyl points. This paper is only concerned with fracturing of lattices with no Weyl points.

Simulation: lattice fracturing due to of external load
In simulation, we build up our system by connecting nearest neighbor (NN) point masses with harmonic springs that are free to rotate around both of their joints. The hamiltonian of our system is therefore where k is the spring stiffness which is set to one in the simulation, r i is the position of site i and l i,j is the rest length of the spring between site i and site j. The geometry of sites and springs are shown in fig.1 and fig.2 for both the deformed square and the deformed kagome lattices.
To simulate the fracturing process of the lattice, we impose uniaxial strain in the vertical direction to the lattice boundary by holding the lattice sites on the bottom boundary stationary vertically while displacing the lattices sites on the top boundary with a distance δh = hγ up, while h is the height of the lattice and γ is the strain that we impose. The lattice sites on top and bottom boundaries are allowed to slide along the horizontal directions of the boundaries. The lattices have a width w and open boundary conditions on the left and the right boundaries. With the boundary condition defined above, we relax the lattice to an energy minimum such that the force of springs are balanced on all internal lattice sites [35]. Thus, the vertical degrees of freedom of the top-and bottom-boundary sites correspond to ∂V , and their horizontal degrees of freedom, along with all internal sites correspond to V in our discussion of the theory.
At every strain step (given γ), we examine force on springs and compare the amplitude of the force to a compressive strength f c and a tensile strength f t . Both f t and f c define the limit force that a spring can bear, beyond which the spring breaks and is removed from the lattice. When we need to break and remove springs after relaxation, we do so and retake the relaxation method for the new lattice and further remove springs that are beyond compressive/tensile limits. We repeat this process until force balance is reached with all springs within limits. We then proceed with the next strain step. We set f t = f c ≡ f 0 = 10 −2 1 in our simulation for all springs in the lattices.

Results
We begin by testing the lattices' linear response to uniaxial load, as shown in fig. 3. The simulation protocol is discussed in Sec. 3, and to obtain linear response, we use very small strain (γ = 10 −3 1) so no bonds are beyond compression/tensile force limit. We not only include perfect lattices as shown in fig.2 but also include lattices that have cracks in the bulk. In fig. 3, we show the stress distribution of linear response for the lattices. Consistent with theory, high stress appears at the domain walls on the right of the lattice, where the topological polarization ensures a topological charge of ν S ∝ −h < 0, an indication that the domain wall acquires localized states of self stress that are protected by the lattice topology. There is no stress localization on the left domain wall, where the topological polarization gives a topological charge of ν S ∝ h > 0, indicating protected floppy modes instead of localized states of self stress. We also find that unlike normal materials, in which stress is amplified at crack tips [33,36], in these topological Maxwell lattices there is no obvious stress at the tips of the cracks. Thus, these cracks are protected by the self-stress domain walls. In addition, we observe no significant difference of stress distribution between the deformed square and deformed kagome lattices, and between lattices with and without cracks. The overall stress distribution in the linear regime is mainly determined by the topology of the lattices bulk, regardless of the microscopic details and small defects that may exist in the lattices.
We then test the entire process of fracturing of the lattices in the quasi-static limit, where kinetic energy of the lattices is quickly dissipated, much faster than the relaxation process of the lattices. We increase the strain by small steps such that at each new strain at most one bond breaks initially (there may be subsequent avalanches when we recalculate the stress distribution after the initial bond breaking but in our observation all avalanches are small events with number of breaking bonds of O(1)). For all lattices shown in fig. 3, we find that the fracture begins in one of the bonds at the self-stress domain wall where stress localizes, as one expect from the linear theory.
Interestingly, the next bond breaking events continue to concentrate near the self-stress domain walls ( fig. 4), even when the lattice is rather deep into the failure process and the coordination number becomes z < 2d. This is a manifestation of the robustness of topological protection -with small damages at the domain wall, the phonon band topology of the bulk of the lattice is unchanged and still dictates the exponential localization of states of self stress and floppy modes.  It is worth noting that localized stress distributions that appear to be similar to those shown in fig.3 can also emerge without a topological origin. An example of this is shown in fig. 5(a), where we show a lattice with R T = 0 in all domains but still shows localized states of self stress. The structure contains twisted kagome lattices [3] of opposite direction of twisting in its neighboring domains. In particular, the pointing up triangles rotate clockwise in the two domains on the left and right and rotate counterclockwise in the middle domain, leading to two domains walls. There is stress localization at the left domain wall. This, however, does not conflict with the theory of topological mechanics. In this structure, the left domain wall has topological charge ν S = 0, which indicates equal numbers of floppy modes N 0 and states of self stress N SS . As long as N 0 = N SS , they can both localize at the domain wall. Indeed, our calculation through the null spaces of matrix C and Q show pairs of states of self stress and floppy modes at this domain wall. We marked one of the floppy modes in fig.5(a) with blue arrows in the kagome lattice. To explicitly show that the nonzero N SS and N 0 here are not a consequence of the lattice topology, but rather the local geometry of the lattice in the region, we cut a thin strip of the lattice that contains the localized states of self stress but doesn't include the bulk of the lattice. As shown in fig.5(b), the linear response of the thin strip still contains the localized states of self stress, despite the fact that it is cut from the lattice bulk. For comparison, we take the deformed kagome lattice shown in fig.3(c) that has nonzero R T and a domain wall of ν S < 0, and cut a thin strip containing this domain wall out of the lattice bulk. We then test the linear response of the thin strip alone, which is shown in fig. 3(d). Interestingly, the states of self stress in the thin strip disappear when it's cut off from the lattice bulk. This is consistent with the theory of topological mechanics that the topologically protected states of self stress comes from the bulk phonon band topology, and when the domain wall is isolated out the exponentially localized states of self stress no longer exist. A similar demonstration for finite-frequency edge states has been shown in Ref. [6] in a gyroscopic system.
Next, we test the effect of embedding multiple domain walls in the lattices. Topologically protected states of self stress are characterized by their decay rate κ as we discussed in Sec. 2. Only when the bulk domain depth is greater than 1/κ, the state of self stress is well defined to be exponentially localized. Therefore, two domain walls that are separated by a distance ∆w 2/κ can be treated as independent. This allows us to design N DW domain walls that are independent to each other in a single lattice provided that the lattice depth w > N DW ∆w.
The advantages of designing mechanical metamaterials with multiple self-stress domain walls are the following. First, the self-stress domain walls "attract" stress and protect regions in between. This is shown in fig. 6(a-c) where we compare a lattice with multiple self-stress domain walls and a deformed kagome lattice with the same type of unit-cell geometry but only one domain. Under the same vertical tension, the lattice with domain walls only have high stress at the domain walls and the regions in between bear very low stress, whereas the one-domain lattice is homogeneously stressed. In fig. 6(c) we show quantitatively the bond force comparison between the two lattices. Second, by controlling the density of these domain walls we can also control the elastic modulus, E ∝ N DW , under the condition that domain walls are independent, as shown in eq. (11) and in fig. 6(d). Third, because the exponentially localized states of self stress are topologically protected, as we discussed above, they continue to attract stress deep into the fracturing process, leading to a gradual failure, in contrast to catastrophic failure in conventional brittle materials ( fig.7). Videos of our simulations of the fracturing process of Maxwell lattices with one or multiple domain walls, as well as brittle failure of the regular kagome lattice are included in the Supplementary Materials.

Conclusion
To summarize, we investigated how Maxwell lattices with domain walls hosting topologically protected states of self stress and floppy modes fracture under applied stress. We find that bond breaking events concentrates near self-stress domain walls, providing protection to the lattice bulk, even deep into the failure process.
Our results open the door to the design of high strength mechanical metamaterials based on topological mechanics. By controlling the line density of the domain walls, we can control both the elastic moduli and the fracturing process of the structure. We show that in the presence of self-stress domain walls, the bulk of the lattice is protected from fracturing, even when small cracks exist in the bulk. This is a useful property that can be used in structures where perfect periodicity in the bulk needs to be protected for functions, such as wave-manipulating acoustic metamaterials.
In this Appendix, we show the derivation of eq. (11) for the spring tensions t in the bulk of the lattice. It has been shown in the main text that t must be a superposition of states of self stress in the bulk, i.e. where {r i } are the set of eigenvectors for nonzero eigenvalues of matrix Q V . We then plug in eq. (9) in the main text to rewrite e using u V and u ∂V , i.e.
This can be simplified by decomposing the spring constant matrix K into the null space and the orthogonal complement space of Q V into K = K rr K rs K sr K ss , (A. 6) where (K rr ) ij ≡ r i · Kr j , (K rs ) ij = (K sr ) ji ≡ r i · Ks j and (K ss ) ij ≡ s i · Ks j . Note that K rr is invertible because its eigenvalues are nonzero. We then simplify our expression by denoting a i with a vector a such that (a) i ≡ a i , and have the equations for the linear coefficients as a = K sr P r C V u V + (K sr P r + K ss P s )C ∂V u ∂V , (A.7) 0 = K rr P r C V u V + (K rr P r + K rs P s )C ∂V u ∂V , (A. 8) where P s is the projection operator from the original bond label space into the states-ofself-stress space (null space of Q V ), and P r is the projection operator from the original bond label space into the orthogonal complement space. We have also used the property that C V s i = 0 for all states of self stress i. From eq. (A.8), we obtain P r C V u V = −K −1 rr (K rr P r + K rs P s )C ∂V u ∂V . (A.9) Finally, by plugging this relation into eq. (A.7), we obtain our linear coefficients a = −K sr K −1 rr (K rr P r + K rs P s ) + K sr P r + K ss P s C ∂V u ∂V , (A. 10) which can be further simplified as a = (K −1 ) ss P s C ∂V u ∂V , (A.11) using a matrix identity that (K −1 ) ss = K ss − K sr (K −1 rr )K rs . Putting the coefficient back into eq. (A.1), we arrive at eq. (11) in the main text.