Topologically Protected Steady Cycles in an Ice-Like Mechanical Metamaterial

Competing ground states may lead to topologically constrained excitations such as domain walls or quasiparticles, which govern metastable states and their dynamics. Domain walls and more exotic topological excitations are well studied in magnetic systems such as artificial spin ice, in which nanoscale magnetic dipoles are placed on geometrically frustrated lattices, giving rise to highly degenerate ground states. We propose a mechanical spin-ice constructed from a lattice of floppy, bistable square unit cells. We compare the domain wall excitations that arise in this metamaterial to their magnetic counterparts, finding that new behaviors emerge in this overdamped mechanical system. By tuning the ratios of the internal elements of the unit cell, we control the curvature and propagation speed of internal domain walls. We change the domain wall morphology from a binary, strictly spin-like regime, to a more continuous, elastic regime. In the elastic regime, we inject, manipulate, and expel domain walls via textured forcing at the boundaries. The system exhibits dynamical hysteresis, and we find a first-order dynamical transition as a function of the driving frequency. We demonstrate a forcing protocol that produces multiple, topologically-distinct steady cycles, which are protected by the differences in their internal domain wall arrangements. These distinct steady cycles rapidly proliferate as the complexity of the applied forcing texture is increased, thus suggesting that such mechanical systems could serve as useful model systems to study multistability, glassiness, and memory in materials.


I. INTRODUCTION
Mechanical metamaterials are often created by compatible coupling of unit cells that exhibit a floppy mode, so as to produce a desired global deformation [1]. This design strategy has been employed to create materials exhibiting auxetic responses [2][3][4], collective pattern transformation [5][6][7][8][9][10], multistable structures [11,12], violations of mechanical reciprocity [13], and simple programmability [14,15]. Any incompatibility in cell arrangement implies geometric frustration [16][17][18]. Understanding how local frustration affects the global response is crucial, not only to avoid frustration which may hinder some functionality, but more interestingly for the controlled incorporation of frustration, which enables novel functionality [1]. For example, taking periodic or nonperiodic compatible mechanical metamaterials [15], topological defects may be introduced by rotating individual unit cells, and are useful for steering stresses along designed trajectories [19][20][21]. Frustration between competing antiferromagnetic ground state orientations has been exploited to produce programmable hysteresis in perforated elastic sheets [14], and frustration was found to generate multiple complex chiral patterns in triangular networks of beams [6]. In general, metamaterials incorporating frustration can often accommodate the conflicts in surprising ways, inducing complex and potentially useful mechanical responses.
The design of new properties emerging from frustration [22,23] has also been extensively studied in * carl.merrigan@fulbrightmail.org engineered magnetic nanomaterials based on the ice rule [24,25], and called artificial spin ice (ASI) [26][27][28]. ASI can be helpfully modeled as sets of classical Ising spins arranged along the edges of a lattice, where each spin lies in the plane of the lattice and points away from one vertex and towards another [29]. It is possible to map an ASI to a mechanical metamaterial in the following way (see Fig. 1): We map each vertex of the ASI to a mechanical unit cell which has a single floppy mode. Importantly, the mechanical floppy mode is such that the deformation of each unit-cell edge may be mapped to a binary spin within a ground-state ASI vertex. Neighboring cells that attempt to adopt incompatible floppy modes create local stress, raising the energy of the metamaterial. Thus, the rules for reducing or manipulating frustration in icelike systems relate to the rules of stress compatibility in metamaterials. However, we emphasize that the latter are still continuum, elastic systems. Therefore, we expect that mechanical spin-ice systems, unlike their magnetic counterparts, will accommodate the effects of frustration in ways that are not afforded to a truly binary system, leading to novel phenomenology.
In this paper, we put forward the first step in a broader program to relate concepts and models between spin systems and mechanical metamaterials. Here we begin with the simplest geometry-a two-dimensional mechanical analogue of square ice [30][31][32]. First we explore static configurations of the metamaterial, which we show may be tuned from a binary, pseudo-spin behavior to a continuous response. Then we go on to analyze the response to external driving. A novel and crucial property of our mechanical metamaterial is the possibility of manipulating it by acting on its boundary. Magnetic systems are typi- The cell is constructed from eight framework (k 1 ,black), four interaction (k 2 ,blue), and four compressive (k 3 ,red) springs with stiffness k 1 , k 2 , and k 3 and relaxed lengths l 1 , l 2 = √ 2l 1 , and l 3 < 2 √ 2l 1 , respectively. (b) Energy as a function of the right triangle rotation angle θ for l 3 /l 1 = 0.95 × 2 √ 2. The undeformed state (red) and the two ground states (green) are denoted by solid squares. (c) Ground state of a 3 × 3 network. A single unit cell is highlighted in green. (d) Mapping of the unit cell displacements to binary spins within square ASI.
cally driven by using a bulk magnetic field [33][34][35][36][37] which acts on all the spins at once. In contrast, for our mechanical system, it is natural to push and pull on the faces of the unit cells at the metamaterial boundaries. When driven in a compatible way at the boundaries, our mechanical metamaterial exhibits a dynamical phase transition from asymmetric to symmetric hysteresis loops, as a function of the driving frequency, as well as the existence of a wide variety of distinct, topologically-protected steady cycles under textured driving protocols. This nontrivial dynamics arises from the topology and motion of domain walls in the system and could be developed as a dynamical form of memory [38][39][40].

II. MECHANICAL SPIN ICE
We design elastic networks to mechanically mimic square ice. Four rigid right triangles, combined as a square, exhibit a floppy mode in which all the triangles may freely rotate together by some arbitrary angle θ without stretching any bonds [see Fig. 1(a)]. This gen-eral rigid mechanism is the basis for a variety of metamaterial designs [1], such as auxetic and chiral metamaterials [3,4,9], metamaterials with propagating internal domain walls [10], and also pattern-transforming perforated elastic sheets [5,7,14]. If the face nodes where the triangles meet are treated as pseudo-spins, then this global mechanism maps to the ground states in magnetic ASI [ Fig. 1(d)], for which each vertex has two opposing spins pointing inwards and two spins pointing outwards. This "ice-rule" constrains the ground states of spin ices to compatible combinations of vertices with zero topological charge, that is, the number of spins pointing towards the vertex equals the number of spins pointing away from the vertex [26].
Our square cell is composed of three types of springs [ Fig. 1(a)]. Two are standard in assuring the floppy modes: "Framework springs" of stiffness k 1 and relaxed length l 1 form the boundary of the square, while "interaction springs" with stiffness k 2 and relaxed length l 2 = √ 2l 1 connect the four faces, completing the four right triangles. In ASI, spins are strictly binary, with magnetization per spin that can point in two possible directions but which is fixed in magnitude [26,27]. However, in mechanical metamaterials the pseudo-spin magnitudes may vary continuously, and in typical designs the undeformed state (θ = 0 in Fig. 1(a)) with all spins equal to zero is the unique ground state. Furthermore, deformations may continuously vary over large length scales, so that crisp, topological-defect structures do not appear spontaneously. Consequently, usually frustration comes into play only when a mechanical system is driven [10,15,19,20]. To obtain bistability, we introduce diagonal "compressive springs" with stiffness k 3 and length l 3 . By requiring l 3 /l 1 < 2 √ 2, the square cells experience a local compressive stress that causes them to act as bistable hysterons [38], with two states ±θ * ( Fig. 1(b)). Each edge of the square cell prefers to point in or out by some fixed magnitude, causing the overall metamaterial to mimic a spin ice with discrete states. In order for the compressive springs and the framework springs to simultaneously achieve their relaxed lengths, the geometry determines the deflection angle to be θ * = cos −1 [l 3 /(2 √ 2l 1 )]. We consider square mechanical metamaterials composed of L x columns and L y rows of bistable square cells, and we take L x = L y . The low-energy, metastable configurations of this metamaterial are governed by tuning the stiffnesses of the three spring types used to construct the unit cell, of which one can choose two independent dimensionless ratios. For simplicity, we adopt units such that k 1 = 1, so that the values of k 2 and k 3 indicate the dimensionless relative stiffness of the interaction and compressive springs compared to the framework springs.
As in the ice-rule obeying ASI, we have two "antiferromagnetically" ordered ground states [ Fig. 1(c)] in which cells alternate orientations. We quantify the order in the Staggered magnetization M αβ for the corresponding configurations. As k 2 is increased, the domain wall curvature is reduced and the gradients in the magnetization along the boundary between competing orientations become more smooth. system using the following staggered magnetization: where ∆x αβ is the width of the unit cell, ∆y αβ is its height, and δ = l 1 sin(θ * ) is the displacement amplitude of the face nodes in the ground state [ Fig. 1(a)]. M has value 1 on all the cells of one ground state, −1 on the other, and 0 on undeformed cells, and might be greater than 1 under elastic deformation. Then, the globally averaged magnetization is (2)

III. METASTABLE CONFIGURATIONS
We numerically simulate our metamaterial within an overdamped dynamics, subsuming any physical viscosity in our definition of time, see Appendix A. In the numerical simulations shown below we use l 3 /l 1 = 0.95 × 2 √ 2, so that θ * = 18.2 • and δ/l 1 = 0.31.
Starting from a neutral state M = 0 where all cells are square, the system relaxes into metastable configurations corresponding to domains of different orientation of M αβ , separated by domain walls, see Fig. 2. Different lowenergy, metastable configurations of the metamaterial are obtained when uniformly distributed, random perturbations ∆x i , ∆y i ∈ [−χ, χ] are applied to each node of the lattice (here χ/l 1 = 0.4). Distinct minima are then found by repeating the simulation for different random initial conditions.
By gauging the dimensionless interaction stiffness k 2 , we can control how the metamaterial mimics its binary counterpart and therefore whether its excitations match the topologically charged vertex configurations characteristic of ASI [30][31][32]. Moreover, with our continuous, elastic system, we can also investigate regimes that go beyond strictly binary behavior, evincing new phenomena that would be absent in magnetic ASI.
First, consider k 3 1 k 2 . The large value of the compressive stiffness k 3 constrains the corners of each unit cell to move inwards, so that the sides of the square unit cell have length less than 2l 1 . Therefore, in order to relax the framework springs, the four nodes at each face need to buckle in or out with amplitude ±δ. For such configurations, only the interaction springs are strained: applying a spin approximation, 2 4 configurations of the unit cell are possible (see Appendix B, Fig. 7). The icerule is obeyed for two of these configurations and they correspond to the energy minimum, as in Fig. 1(b). The ice-rule is violated by excitations for which two neighboring faces both point inwards or outwards, causing the interaction springs to be strained. We note that with large k 3 compressing the interaction springs costs more energy than extending them, and this alters the vertex energy hierarchy compared to ASI (see Appendix B). When, k 2 ∼ 1, the pseudo-spin approximation breaks down and a wider variety of continuous excitations may appear in which stresses are shared between both the framework springs and the interaction springs.
The dependence of the excitation morphology on the dimensionless interaction stiffness k 2 is demonstrated in Fig. 2. The compressive stiffness is large, with value k 3 = 5 throughout, and k 2 varies from 0.01 to 2. An experimental realization would require component springs with a variable stiffness range of order 500 : 1, which may readily be achieved with standard materials. For a small dimensionless stiffness k 2 = 0.01 of the interaction springs, the domain walls have shapes very similar to ASI walls: one cell thick, and made of configurations in which two neighboring faces of the same cell point in the same direction [41]. We also see other topological defects, especially at the edges of the lattice, corresponding to configurations in which three faces deform in the same direction, and are analogous to ASI monopoles (see Appendix B, Fig. 8). This suggests that such topological charges may be expelled toward the boundaries of the system, similarly to what was shown for colloidal ice [18,[42][43][44]. As k 2 increases, such "monopoles" become more energetically costly and appear within the domain walls. The curvature of the latter grows smaller with increasing k 2 . Domain walls seen in Fig. 2(b), with k 2 ∼ 0.05, resemble domain walls seen in experiments on square, magnetic ASI [30,31]. As k 2 increases further, the domain walls become straight, and they become wider with strains distributed over neighboring rows or columns of unit cells.
During the energy quench, all domain walls tend to move and straighten for some time until stabilizing at some average curvature imposed by the value of k 2 . Further, any domain walls that connect two adjacent edges of the lattice and start near a corner tend to be expelled. Corner spanning domain walls are expelled at a decreasing speed as k 2 is decreased, until k 2 ∼ 0.18, where the corner spanning domain walls become stationary and are no longer expelled (see Appendix C, Fig. 9). At values of k 2 which approach 1, corner spanning domain walls tend to be eliminated rapidly, and the resulting configurations only contain straight domain walls, such as in Fig. 2(d).

IV. DYNAMICAL HYSTERESIS
We now investigate how this rich topological zoology affects dynamics. In square spin ice, monopoles and domain walls can be rearranged by external magnetic fields. While an analogue of a bulk magnetic field acting on all unit cells at once is not readily available for mechanical systems, we can focus on the experimentally realizable case of forcing the boundaries of the mechanical metamaterial. First, we apply forces to the nodes of the cells at the boundary in a compatible way, i.e by alternating the sign of the force among consecutive nodes (see inset of Fig. 3) and modulating the intensity harmonically in time F = F 0 sin(ωt), where F 0 is the force amplitude and ω is the dimensionless driving frequency. The characteristic relaxation time for the framework springs, T 1 = γ/k 1 = 1, where γ is the linear drag coefficient, is taken as the unit of time (see Appendix A), and we study the dependence on the dimensionless driving frequency approaching the quasistatic limit, ω 1. When driven in such a manner, the system undergoes a hystersis cycle, with the shape of the hysteresis cycle governed by the dimensionless driving frequency ω, see Fig. 3. When the intensity of the force F (t) overcomes a coercive value F c , the boundary cells all snap from M αβ = +1 to M αβ = −1, creating a domain wall loop sitting just inside the edge of the metamaterial. Once formed, this loop starts to contract to reduce its overall length and energy, see Fig. 3 inset. The subsequent internal dynamics depend on how the driving frequency ω compares to the speed of the internal domain wall contraction.
At high ω, the boundary force oscillates too rapidly for the system to mechanically respond, leading to a shallow penetration depth of stressed cells and leaving the bulk unchanged, with M (t) ∼ 1 at all times. For less rapid forcing, edge cells invert, creating a circumferential do-main wall. This domain wall begins to contract toward the center of the system, but it is subsequently annihilated by a second domain wall forming in the second half of the cycle, preventing it from penetrating into the bulk. Sample time series are shown in Appendix D, Fig. 10, and corresponding videos are provided in Appendix F, Videos 1-4. The fact that the second loop catches up with the first is evidence of an attractive interaction between the two loops, since otherwise each loop would follow the exact same trajectory, differing only by a time delay. Once the driving period τ = 2π/ω is just large enough for the first domain wall to propagate inwards without being caught by the subsequent domain wall, a first-order dynamical transition appears at which the hysteresis loops change abruptly from being asymmetric, to being symmetric about M (t) = 0. For the k 2 = 0.2 hysteresis loops shown in Fig. 3, the transition occurs as the driving frequency is reduced from ω = 4.5 × 10 −6 to ω = 3.5 × 10 −6 . This behavior resembles the dynamic phase transition studied in kinetic Ising models driven at finite frequency by an external magnetic field [33][34][35][36][37], except that the transition we observe appears to be discontinuous. We also note that it is natural to inject and manipulate domain walls via the boundaries of a mechanical metamaterial, whereas magnetic systems are typically driven via a bulk magnetic field. In the quasistatic limit, ω → 0, the forcing period τ is far larger than the time needed for the internal loop to contract. Thus, the hysteresis loop is square.
We quantify the hysteresis through its normalized area A(ω) = 1 4F0 M dF and the time-averaged magnetization over one period Q(ω) = 1 τ M dt [33][34][35][36][37]. Q = 0 when the cycle is symmetric, and Q = 0 indicates an asymmetric cycle. In Fig. 4(a,b), we plot A and Q vs. ω for a set of k 2 values ranging from the binary, spinice limit to the elastic regime. For k 2 ≥ 0.1, the dynamic transition from asymmetric to symmetric cycles is marked by a discontinuous jump in both A and Q, with Q remaining zero below the critical frequency ω c . ω c can be related to τ c , the time it takes a circular domain wall to collapse, as ω c = π/τ c . When ω < ω c full inversion is achieved by the collapse of the circular domain wall before a new half cycle begins. When ω > ω c a new domain wall is created before the previous can collapse, and M doesn't fully invert sign, leading to asymmetry.
Another interesting frequency is the resonant frequency ω r , which corresponds to the maximum in the area A. Here, just after reaching orientation M = −1, the system is driven to immediately return to orientation M = +1, once again. The maximum area, resonant frequency loop for k 2 = 0.2, ω r = 1.28 × 10 −6 is shown in pect model-A dynamics to describe the domain wall motion [45]. For a circular domain wall, it has been shown that the circle contracts at a rate proportional to its mean curvature: dR dt ∝ 1 R [46,47]. Consequently, the rate of change of the area inside the circle

V. TEXTURED FORCING PROTOCOL
We have demonstrated how the domain wall dynamics control the probe response of the system in the very simple case of a compatible drive. However, our system allows us to "inject" desired domain walls for different dynamical responses, by finely controlling the detail of the applied force pattern. For instance, wherever we skip the alternation of sign of the applied force along the boundary we inject a domain wall. This can only be done an even number of times, and each will correspond to a domain wall ending at the interface.
In the compatible forcing protocol, only one static metastable configuration matches the forcing texture at maximum (minimum) amplitude of the textured boundary forces, at times such that ωt = π/2 (ωt = 3π/2), namely, the global ground state M = +1 (M = −1). Many such matching, metastable states are made possible for an applied forcing texture composed of alternating M αβ = ±1 segments (see boundary arrows in Fig. 5). In Appendix E, we demonstrate that the number of these states is given by the Catalan numbers [48], which grow nearly exponentially as a function of the number of alternating force segments. Below, we demonstrate that different initial conditions with different internal domain wall configurations may be attracted to topologicallydistinct steady cycles. These cycles are always between a subset of the static metastable configurations (and their inverses) for which the boundary configurations and the arrangement of internal domain walls match compatibly with the forcing texture at its extremum values.
First, we study a relatively simple textured forcing protocol. Figure 5 illustrates three steady cycles which appear when the system is driven with 8 alternating forcing segments. Metastable states which match this forcing texture must contain 4 internal domain walls. There are 14 such static configurations: three of these distinct states appear in the second column of Fig. 5, panels (b), (f ) and (j), while their inverses appear in the last column: (h) is the inverse of (b), (d) is the inverse of (f), and (l) is the inverse of (j) rotated by π. There are in fact 6 states with nondegenerate domain wall configurations, shown in Appendix E, Fig. 11, and all 14 topologicallydistinct states may be generated by applying π/2 rotations to these 6 states.
We tested a large number of different initial conditions: for each initial condition we imposed 8 alternating M αβ = ±1 segments as fixed boundary conditions, thus giving one of the 14 topologically-distinct metastable configurations as an initial condition. At high frequencies, the domain wall penetration depth is shallow, and a variety of different steady states appear, reflecting the different internal domain wall configurations of the initial conditions. These can be seen by the different Q values at large ω in Fig. 5(n). Interestingly, in the quasistatic limit, the internal domain wall configurations are erased, yet the system tends to one of three steady-state M (t) times series, each giving a signature Q value. Figure 5(m) shows these three distinctive magnetization time series. The panels in each column show the magnetization configurations M αβ at times ωt = 0, π/2, π, and 3π/2 in the steady state, when the textured forcing amplitude is zero or extremal.
For Cycle 1, M (t) varies asymmetrically between 0.3 and 1, with time-averaged magnetization Q = 0.6. In Cycle 2, M (t) varies between −0.3 and −1, and Q = −0.6. This cycle mirrors Cycle 1, but with a phase shift of π. This phase shift can be understood by comparing the second and last columns of Fig. 5, when the force is at it maximum (minimum) value. There are two distinct topological configurations that match the boundary conditions, one with domain walls spanning the four corners [ Fig. 5(b) Fig. 5(j) give four distinct domain wall arrangements, so we consider Cycle 3 to have four variants, giving a total of 6 distinct steady state cycles arising for this simple textured forcing protocol.
For the compatible drive studied in Sec. IV, starting from an initial condition with complicated internal domain wall configurations will not lead to multiple cycles. When compatibly cycling initial states containing several domain walls, domain walls within the starting configurations are always annealed away after some transient number of cycles. Eventually, the system always returns to the same hysteresis loops of the previous section, and all memory of the starting state is erased. The duration of long-lived transients grows longer for larger ω. Similarly, if the 8 segment driving protocol of Fig. 5 were applied to a random configuration containing some complex domain wall pattern not necessarily compatible with the applied force, one of the 6 distinct steady states cycles will eventually be reached.
The ultimate cycle a given initial condition is attracted to depends critically on the initial domain wall topology: for example, an initial condition needs to have at least one domain wall spanning the width of the metamaterial to access Cycle 3. However, the initial topology does not guarantee the resulting cycle that the state will be attracted to, and finer details of the internal domain wall curvatures and small variations in the internal configuration determine the cycle which a particular initial arrangement is attracted to. For example, we found that different initial conditions each with two long domain walls [see Appendix E, Fig. 11(d)] sometimes reach Cycle 3, and sometimes reach Cycle 1 or Cycle 2. Finer details of the internal domain wall arrangements can also affect the stability of the distinct cycles. We note that Cycle 3 appears to be less stable: for some initial conditions, it appears to be the steady state over a wide range of frequencies, but then an abrupt switch may occur to Cycle 1 or Cycle 2. Except for the anomaly at ω ∼ 10 . L x = L y = 45 and k 2 = 0.5, and the cycles were found for fixed small frequency ω = 1.9 × 10 −6 (the M αβ colorbar is the same as in Fig. 2). (i) Magnetization time series corresponding to the snapshots shown above illustrating the rich variety of steady state hysteresis cycles arising from this textured forcing protocol.
at least down to ω ∼ 10 −8 . Sample videos of Cycle 1 and Cycle 3 are given in Appendix F, Videos 5 and 6. Finally, we demonstrate that the variety of the distinct hysteresis cycles can be easily increased by adding more alternating forcing segements along the boundary. Doubling the number of boundary defects, we tested a forcing protocol with 16 alternating M αβ = ±1 segments along the boundary. Compatible initial conditions for this forcing protocol must have 8 internal domain walls, and there are a total of 1, 430 distinct ways to place the domain walls (see Appendix E). Among these, each state has some trivial degeneracy of 1, 2, 4, or 8 from rotations or reflections (see Appendix E, Fig. 12). A simple lower bound may be attained by assuming the maximum degeneracy for each state, implying that there are at least 1, 430/8 ≈ 178 different domain wall configurations.
By fixing the boundary conditions and quenching the energy as in Sec. III, we manually identified 63 distinct states out of a sample set of 2, 000 trial cases, and we expect additional rare states would be found for large ensembles or could be directly imposed on the metamaterial. Snapshots of these energetically distinct domain wall configurations are given in Appendix E, Fig. 13 and Fig. 14. Taking these states as initial conditions, and driving them at low frequency ω = 1.9 × 10 −6 , we observe many distinct steady state cycles, 8 examples of which are shown in Figure 6. Two panels for each cycle show the topologically-distinct configurations which appear at the maximum and minimum textured forcing amplitude. Some of the cycles, like those shown in (a), (f) and (g), are relatively common and reached by many of the trial initial conditions. Others, like (b), (c), (d), and (e) are more rare. Sample videos of the cycles corresponding to Fig. 6(a-d) are provided in Appendix F, Videos 7-10. Additional cycles occurred which are rotational and phase-shifted counterparts to these example cycles, and we expect to find more cycles by testing a larger ensemble of initial configurations.

VI. DISCUSSION
We have demonstrated that our mechanical analogue of spin-ice exhibits real-space topological defects, similarly to magnetic and other realizations of ASI. These domain walls tend to escape to the system boundaries, but they can sometimes be immobilized depending on the ratios of spring stiffnesses used to construct the metamaterial or on imposed boundary conditions. As the dimensionless interaction stiffness k 2 is reduced, the domain walls become slower and more curved, until their motion ceases in the ice-like limit k 2 k 1 . Uniform inversion of the boundaries creates circular domain wall loops which contract upon themselves, setting up a dynamic hysteresis cycle between the competing ground state configurations. There appears to be a first-order dynamic phase transition associated with this uniform driving protocol. The system exhibits a novel dynamic complexity when textured forcing is applied at the boundary: different initial configurations of the network are attracted to multiple, topologically-distinct steady cycles. The number of cycles is closely related to the multiplicity of topologicallydistinct static states consistent with the texture of the boundary forcing. Further, the possibility of many distinct ways of connecting domain walls inside the material protects these cycles: once a steady cycle is achieved, the system cannot easily be pushed into a new cycle because of the large energy barriers which would need to be crossed in order to change the positions of domain walls from their current arrangement.
This potential for many steady cycles arising for a single driving texture relies on the ability to manipulate the boundary of the metamaterial, unlike ASI and other systems which are typically driven by bulk fields. Thus, our mechanical spin-ice system opens up further exploration of rich dynamical phenomenology which have no counterparts in magnetic systems. Our work opens a new direction of using periodic prestressed mechanical metamaterials in order to generate complex energy landscapes, which in turn lead to multiple steady states. It should be fruitful to explore additional driving protocols, as well as modified lattice geometries. Additional rich hysteresis behavior could be generated by driving opposite pairs of the boundaries out of phase from one another, or driving different sides with different amplitudes. Specific domain wall configurations could be imposed by programmed protocols that manipulate parts of the boundary one by one in a specific sequence. Moreover, the degeneracies of distinct metastable configurations inside the material can be engineered by changing the shape of the boundary or the spacing of forcing defects along the boundary, possibly allowing for more rich behaviors. Finally, it should be interesting to consider more elaborate spin-ice mechanical metamaterials, but based on lattices that are inherently frustrated and thus have extensively degenerate ground state configurations [49,50]. Studying such lattices could promote the understanding of the complex dynamics and memory of amorphous solids and other glassy systems. As such, it will be useful to study the defects and possible dynamic hysteresis in such generalized mechanical spin-ice systems. Our mechanical metamaterial is modeled as an elastic spring network simulated using overdamped dynamics or the method of steepest energy descent. For a network of nodes with equal mass m, the overdamped equation of motion for a given node with position r i is where the spring connecting node i to node j has spring constant k ij , relaxed length l ij , and current extension e ij = | r i − r j |, and γ is the linear drag damping force coefficient. We adopt units such that relaxation time scale T 1 = γ/k 1 = 1, and assume the system has an overdamped response such that the viscous damping time scale is much faster than the unit time m/γ γ/k 1 = 1. In the characteristic relaxation rate for interaction springs is T 2 −1 = k 2 and for the compressive springs is T 3 −1 = k 3 . Thus, the fastest relaxation rate in the system is given by k 3 −1 = 0.2. We integrated the equations of motion with variable time step ∆t = 0.01 − 0.1. For characterizing the various metastable conditions in Sec. III, we use free boundary conditions, with no constraints or additional forces on the nodes. In Sec. IV and Sec. V, we drive the lattice from the boundary by applying additional forces to the boundary nodes, perpendicular to the boundary, with magnitude |F i | = F 0 sin(ωt), and alternating sign as determined by the desired driving texture, taking F 0 = 1.5 throughout.

Appendix B: Spin Ice Vertex Classification
Within ASI systems, spin configurations at a vertex are classified by their energy and their topological charge q, which is defined as the number of spins pointing into the vertex minus the number of spins pointing away from the vertex. Similarly, in the limit that the interaction springs are much weaker than the framework springs, and both are much weaker than the compressive springs, k 2 1 < k 3 , we can obtain an approximate vertex hierarchy in energies for the square unit cells by treating the four faces of the unit cell as pseudo-spins. In this limit, we require that the two diagonal compressive springs and the eight framework springs must all be fully relaxed. This means that each of the four face nodes of the unit cell must bend inwards or outwards by an amount ±δ = l 1 sin(θ * ). As discussed in the main text, all of the springs are relaxed when two opposing faces of the unit cell bend inwards and the remaining two faces bend outwards, and the total elastic energy is zero. Excitations happen when two adjacent faces both point outwards or both point inwards, causing the interaction spring connecting the two faces to be stretched or compressed relative to its relaxed length l 2 = √ 2l 1 . Figure 7. shows the strain on all springs for all six possible vertex configurations. Unit cell configurations very similar to these idealized vertex configurations appear within the bulk of the metamaterial, as shown in Fig. 8. Due to the unit-cell geometry, the energetic cost of compressing an interaction spring is slightly larger than the cost of stretching the same interaction spring. For a given value of θ * = cos −1 (l 3 /2 √ 2), the energy of a single interaction spring is given by E/k 2 = 0.5∆l 2 = ( √ 1 ± sin 2θ * − 1) 2 , where the plus sign is for a stretched interaction spring and the minus sign for a compressed spring. Simulations in the main text were carried out using θ * = 18.2 • . Adding up the energy contribution from each stressed interaction spring, we find the following energy hierarchy for the vertices: E a /k 2 = 0, E b /k 2 = 0.1376, E c /k 2 = 0.2, E d /k 2 = 0.2626, E e /k 2 = 0.2752, and E f /k 2 = 0.5252. We may also assign the topological charges q a = 0, q b = −2, q c = 0, q d = +2, q e = −4, and q f = +4. Note that in magnetic ASI, the energy for the q = ±2 and q = ±4 monopoles is independent of the sign of the charge q, whereas this degeneracy is lifted due to the asymmetry between stretching and compressing the interaction springs in our square mechanical ice.

Appendix C: Domain Wall Speed
In order to quantify the domain wall speed, we compare the gradient descent dynamics for configurations starting with a long corner-spanning domain wall. We use a system of size L x = L y = 31, and prepare the initial configuration with a long domain wall spanning two adjacent sides and starting just above the diagonal of the metamaterial. During the dynamics, the domain wall decreases its length by escaping towards the corner of the metamaterial. In Fig. 9., we plot the time dependence of the contour length l(t)/l d , where the diagonal length is l d = 2 √ 2L x , of the M αβ = 0 contours for decreasing dimensionless interaction stiffness k 2 . The escape time grows larger for decreasing k 2 , until below a value k 2 ≈ 0.18, the initial corner-spanning domain wall remains in place forever without escaping to the edge. These results are consistent with the change in the nature of the hysteresis cycles seen in Fig. 4, where for k 2 = 0.01 − 0.05, the hysteresis cycle is always restricted to perturbations along the boundaries even in the quasistatic limit ω → 0.  Fig. 10. illustrates several time series of the magnetization M αβ for the hysteresis cycles shown in Fig. 3. The first four panels, (a)-(d), show the trivial high frequency behavior for which the domain wall never penetrates into the bulk of the metamaterial. For a frequency slightly larger than the critical frequency ω c , panels (e)-(h), the domain wall loop created on the second half of the cycle catches up and annihilates with the first domain wall loop, creating a quadrupolar pattern. Unit cells deep within the center of the lattice are never flipped, keeping their magnetization value M αβ = +1. In contrast, for a frequency just below the critical frequency, panels (i)-(l), the two circular loops move concentrically without colliding. At the resonant frequency, ω r , panels (m)-(p), the first circular loop contracts completely just as the second circular loop is nucleated at the boundary of the lattice. Finally, in the quasistatic limit ω → 0 (not shown), the loops contract rapidly compared to the driving period τ , and so the metamaterial spends most of the time waiting in either state M = +1 or M = −1.

Appendix E: Multiplicity of Topologically Distinct Configurations
It is possible to count the multiplicity of distinct metastable states as a function of the number N of alternating segments of orientation A and B along the boundary. Circling the boundary of the lattice, each defect marking the end of a domain wall is a transition from A to B or from B to A. Since each domain wall must both start and end on a boundary, each domain wall emanating from an AB defect must terminate at a BA defect, so that there are a total N/2 of each defect type. If there were an odd number of defects injected along the boundary, this would imply the existence of at least one domain wall in between that terminates in the bulk, which is impossible. Therefore, we know each AB transition is matched to a BA transition, and we get a simple upper limit on the number of states by counting all permutations between the two sets: Ω max (N ) = (N/2)!.
This upper limit overestimates the number of states because it includes states in which different domain walls intersect, which is forbidden. This no-crossing restriction naturally leads to a recursion relation for the number of allowed configurations. Assume a boundary condition with N + 2 total defects, and pick two of them to be connected. Since no additional domain wall can cross through this wall, the system is divided into two disjoint sections, one section containing m = 0, 2, 4, ..  (E4) Invoking the generating function definition then gives and solving this quadratic relation for g(z), we find that where the sign of the root is chosen so that lim z→0 g(z) = Ω(0) = 1. Finally, the coefficients of the Taylor series expansion of g(z) provide the exact multiplicity for any desired number of boundary segments N . The rapid growth of the number of possible configurations is shown in Fig. 11(g). We also checked that directly counting the topologically allowed Ω(N ) up N = 18 matches the calculated formula. Illustrations of the resulting internal domain wall configurations are given in Fig. 11(a)-(f) as well as in Fig. 12.
The generating function Eq. E6 is in fact known to be the generating function for a well known series known as the Catalan numbers [48], Ω(2n) = C n = 1 n+1 2n n . For n large, the series grows asmptotically as C n ∼ 4 n n −3/2 .
Though the Catalan numbers account for the number of topologically distinct configurations, these configurations may be grouped by energy. The number of energetically degenerate states depends on the shape of the boundary, as well as on the spacing of the defects, both of which will determine the lengths of the internal domain walls. For example, Figure 11(a)-(f) shows that for N = 8 there are a total of 6 energy levels. The remaining topologically distinct states are accounted for by π/2 rotations of the ones shown in Fig. 11. For a square boundary with N/4 defects on each side, states may be related by π/2 degree rotations, or by mirror reflection along the diagonal, so that all configurations for any N have degeneracy 1, 2, 4, or 8. Examples of degenerate states for N = 16 are given in Fig. 12.
Classifying the configurations by energy becomes exceedingly difficult for increasing N . Figure 11 shows that for N = 8, the Ω(8) = 14 topological distinct states can be grouped into 6 energetically degenerate states, 2 × 1 + 2 × 2 + 2 × 4 = 14. For N = 12, there are Ω(12) = 132, and we have checked these may be grouped into 48 energetically distinct states, 8 × 1 + 18 × 2 + 22 × 4 = 132. For the case N = 16, Ω(16) = 1, 430, and enumerating the energy levels becomes exceedingly difficult. Since each state has a maximum degeneracy 8, we obtain a lower bound of 1, 430/8 ≈ 178 distinct configurations. As a check, we have simulated an ensemble of 2, 000 trials, using N = 16, L x = L y = 45, and k 2 = 0.5. In each case we start from a random state, and minimize the energy by gradient descent. Sorting the results of these runs by energy, we find 63 energy levels, see Fig. 13 and Fig. 14, where the states are presented in order of increasing energy. Counting the possible rotations and reflections of these 63 groups, gives a total of 381 out of the expected 1, 430 configurations. The undiscovered states presumably have higher energies, making them difficult to generate via our energy minimization protocol.