Programmable and robust static topological solitons in mechanical metamaterials

Solitary, persistent wave packets called solitons hold potential to transfer information and energy across a wide range of spatial and temporal scales in physical, chemical, and biological systems. Mechanical solitons characteristically emerge either as a single wave packet or uncorrelated propagating topological entities through space and/or time, but these are notoriously difficult to control. Here, we report a theoretical framework for programming static periodic topological solitons into a metamaterial, and demonstrate its implementation in real metamaterials computationally and experimentally. The solitons are excited by deformation localizations under quasi-static compression, and arise from buckling-induced kink-antikink bands that provide domain separation barriers. The soliton number and wavelength demonstrate a previously unreported size-dependence, due to intrinsic length scales. We identify that these unanticipated solitons stem from displacive phase transitions with periodic topological excitations captured by the well-known \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varphi }^{4}$$\end{document}φ4 theory. Results reveal pathways for robust regularizations of stochastic responses of metamaterials.


Supplementary Note 1: Details of finite element (FE) simulations Material constitutive model
The solid material constituting the metamaterial is an incompressible hyperelastic elastomer. In the FE modeling, its constitutive behavior is assumed to follow the Mooney-Rivlin hyperelastic model, with the strain energy density given by 1

FE simulation and spring stiffness
The experimental and simulated responses of the metamaterial to uniaxial compression are represented in Supplementary Fig. 2a, showing excellent agreement with each other.
Note that, to stabilize the simulations, an artificial inter-dissipation factor of 1 × 10 −4 is introduced in the FE modeling. This specific choice is made by trial and error: Much more or less dissipation would result in non-realistic deformed patterns or unstable simulations, respectively. The nominal stress versus engineering strain curves (σ ∼ ε) indicate instability occur during loading, which is confirmed by the sharp drop in σ beyond the threshold strain ε c (see Supplementary Fig. 2b). To capture the symmetry breaking and the formation of periodic solitons in the metamaterial, the characteristic rotation angle of the n-th unit cell in the middle row is defined as θ n ≜ (θ n,l − θ n,r )/2, with θ n,l = arctan ) and ) being the rotations of the lift and right midpoint pairs (i.e., l u ∼ l d and r u ∼ r d ) marked on the vertical necks ( Supplementary Fig. 2c), respectively.
In the experiments, we use a digital camera (SONY-FDR-AX40) to record the deformation process. The rotation angles of the midpoint pairs can be calculated by employing the image processing technology via MATLAB(R2017b, The MathWorks Inc.). Similarly, in the simulation, the coordinates of the midpoint pairs can be extracted, and the corresponding rotation angles can be obtained. We can further introduce the notion of polarization to the metamaterial consisting of bi-stable unit cells, by harnessing the outward and inward polarizations to the deformed unit cell with θ n > 0 and θ n < 0, respectively. Accordingly, the characteristic rotation angle θ depicts the states of the unit cells and, in fact, serves as the order parameter during the structural phase transition.
Based on the experiments and simulations, as shown in Supplementary Fig. 2, the deformations of the metamaterial can be characterized by the stretching of the plates, and the bending and shearing of the necks in the unit cells. As the ligament thickness t is significantly small, the deformation of the metamaterial is mainly concentrated in the neck region 2, 3 . We  a-c The FE model to calculate the stiffness of tension K, torsional C i (i = 1, 2) and shear springs C s , respectively, with the equivalent rod-spring elements and boundary conditions graphed on them. d-e The reaction force/moment versus applied loads curves. The stiffness can be obtained by a linear fitting of the FE results within our considered range of strains, and in present cases, K = 1.478 N · mm −1 , C 1 = 1.150 N · mm, C 2 = 1.069 N · mm and C s = 0.3230 N · mm −1 , respectively. develop rod-spring elements to represent the deformation features ( Supplementary Fig. 3ac). The corresponding stretching, torsional and shearing stiffness are defined as K = F /δ, C 1,2 = M 1,2 /θ and C s = F /u, respectively. Since analytical solutions for these parameters are unavailable due to the irregular neck geometries, they can be calculated by FE simulations. In the simulation, plane strain elements (i.e., 4-nodes with hybrid formulation and

Supplementary Note 3: General static soliton framework On site potential of the unit cell
We calculate the on-site potentials of the unit cell via both numerical and theoretical procedures. In our simulations, two static loading steps are conducted. First, a uniform compressive displacement u y = 2δ is applied to the unit cell ( Supplementary Fig. 7a), which results in the deformed configuration shown in Supplementary Fig. 7b. We denote the corresponding rotation angle and strain energy by θ 0 and U FE cell (θ 0 ), respectively. Secondly, we keep the compressive load and twist the four plates to change the rotation angle (consider the first compressive state as the initial configuration ( Supplementary Fig. 7b)). During the torsion, the current configuration ( Supplementary Fig. 7c), the corresponding rotation angle θ and strain energy U FE cell (θ) can be recorded and extracted via Abaqus/standard software. Setting the configuration with θ 0 as the zero potential state, the on-site potentials of the unit cell can be, therefore, approximately calculated by Based on the thin ligament assumption mentioned in the main text, we propose an effective rod-spring model to capture the deformation feature of the unit cell ( Fig.2 in main text and Supplementary Fig. 7a). It consists of elastic rods connected by neck springs (torsional and shearing springs) whose stiffness can be calculated by FE simulations (Supplementary Fig. 3). As t/a 1 is small enough, the bending of the rods can be ignored. Under the compression u y = 2δ, the strain energy of the unit cell is resulted from the stretching of the rods, torsion and/or shear of the necks springs. Ignoring the asymmetric modes of the unit cell, the energy function can be written as ). In fact, there are numerous geometry-allowed configurations (various θ as shown in Supplementary Fig. 7d) of the rod-spring model under the compression u y = 2δ. However, the actual states (i.e., θ = θ 0 ) should be determined by minimizing U cell via ∂U cell /∂θ | θ=θ 0 = 0. Therefore, the deformations of the unit cell can be regarded as a 'quasiparticle' sited in a potential filed, seeking and occupying the stable state. Analogous to the definition of numerical on-site potential via the simulation procedure, this theoretical on-site potential of the unit cell can be defined by Deformed configurations of the unit cell are related to the features of P cell . As Θ cell ⩽ 0, P cell (δ; θ) only has one stable state θ = 0, which indicates a compression deformation. On the contrary, Θ cell > 0 or equivalently has two stable states at θ = ∓ √ Θ cell , which corresponds to the inward and outward polarizations, respectively. The numerical and theoretical on-site potentials are presented in Supplementary Fig. 7e, showing good agreement with each other to some extent. It should be noted that ε ≜ δ/H 0 and in this Note, θ = 0 is selected as the zero potential state to facilitate the comparisons. Setting Θ cell = 0, we obtain the critical strain ε cell = δ cell /H 0 , above which the on-site potential will be nonconvex.

Structural phase transition and static periodic-soliton
To uncover the physical mechanism of the excitations of static periodic solitons in the metamaterials, a quasi-1D rod-spring model system is constructed (Fig. 3a in the main text and Supplementary Fig. 8a). We consider a finite size model consists of N x × 1 unit cells, as illustrated in Supplementary Fig. 8a. The unit cells are connected by bending and shearing springs (magenta dots in Fig. 3a and Supplementary Fig. 8a). The vertical compressive load u y = 2δ is uniformly applied at the top boundary of the rod-spring model system.

Deformation.
Here we retain the assumption t/a 1 ≪ 1 to obtain analytic solutions.
Therefore, stretching of the rods and, bending and shear of the necks dominate the deformation of this quasi-1D rod-spring model system. We further notice that only bending occurs at the yellow neck due to the symmetry in the vertical direction, while bending and shear are coupled at the magenta necks owing to the incompatible deformations of adjacent unit cells.
These three kinds of deformations at the necks are governed by bending stiffness C i (i = 1,2) and shear stiffness C s , respectively.

Klein-Gordon equation and λφ 4 theory.
For simplicity, we tacitly ignore the asymmetric modes of the unit cell and, therefore, a characteristic rotation angle θ n (n = 1, 2, ..., N x ) can be specified to the n-th unit cell, as defined in Supplementary Fig. 2. Based on this assumption, three typical deformed configurations in the rod-spring model system under the compression u y = 2δ are listed in Supplementary Fig. 8b-d. The total strain energy of the system is related to these deformations by where the first term works as the interactions between the adjacent (n − 1)-th and n-th unit cells, and the last two terms result from the bending and stretching of the n-th unit cell itself. It is noted that this strain energy E sys is dependent on the numerous geometry-allowed configuration of the rod-spring system. However, the actual stable states will be found at stationary points of E sys . Analogous to the analysis in Supplementary Fig. 7, we can regard the unit cells as quasiparticles and introduce a particle-chain model placed in the discrete potential fields to equilibrate the loading process (Fig. 3b in main text). The Lagrangian of the chain system can be written as where the interaction energy E inter ≜ C s s 2 n−1,n = C s d 2 (sin θ n−1 − sin θ n ) 2 and the strain energy of the unit cell and ). Note that the constant D 2 is θ-independent and can be omitted in the following derivations. Therefore, the effective on-site potential of the particle-chain model can be defined as Considering the Taylor expansions sin θ = θ − θ 3 /6 + o(θ 4 ) as θ ≪ 1, the equilibrium equation can be found by ∂L/∂θ = 0, that is, The deformations of the chain model are governed by the features of P eff (θ n ): as Θ eff < 0, P eff (θ n ) only has one stable state and the system undergoes uniform deformations under compression; On the contrary, Θ eff > 0, or equivalently )/ 2 (based on Eq. (S12)). Hence, P eff (θ n ) is a double well potential which indicates the structural phase transitions are allowable in the metamaterials. The critical load δ s = δ | Θ eff =0 can be obtained by the losing of convex in P eff (θ n ). We develop a phase diagram to identify whether the unit cells of the metamaterials can possess monostable or bi-stable on-site potentials. Here, we are interested in the bi-stable phase (Green region in Supplementary Fig. 9a). Note that the effective on-site potential P eff of unit cells in this region is evidently dependent on the loads. Gradually increasing the load beyond the threshold δ s , P eff varies form monostable to bi-stable ( Supplementary Fig. 9b).
As δ > δ s , we can expect two limiting physical regimes based on Eq. (S13): the order-disorder and displacive phase transitions, which essentially depend on the relative magnitudes of the coupling energy and the on-site potential-barrier height.
one has a collection of weakly coupled nonlinear oscillators, randomly displaced to θ ≈ ± √ Θ eff . In this situation, the order-disorder phase transition takes place and induces random localizations in the metamaterials.
• Displacive phase transitions. In contrast, as C s d 2 ≫ (K(H 0 − δ)δ − 2(C 1 + C 2 ))/4 rotation angle θ changes slowly and the standard continuum approximations can be used to reduce Eq. (S13) to where sn is the Jacobi elliptic function, a = √ 2m/(1 + m)Θ eff and b = √ 2/(1 + m)Θ eff . The parameter m can be determined by the boundary conditions, and in present case m = 0.47. As to a finite system, the periodically localized configuration is usually termed as a soliton-lattice with kink width w = √ being the elliptic integral of the first kind.
• Other trivial or more complex excitations. In addition to the two limiting regimes, the trivial excitations such as θ n = ± √ Θ eff are also the solution of Eq. (S13). Evidently, it is the lowest-energy state of the system because all unit cells (quasiparticles) silently reside at the bottoms of a potential well. The present displacive regime reported in the main text is, in fact, one kind of low-energy yet intrinsically nonlinear excitations above the lowest-energy state. Other more complex nonlinear excitations in the multistable metamaterial systems may also be expected, which is still an open question.

General soliton framework to program periodic localizations.
Remarkably, the structure phase transition mechanism revealed in this study, in fact, can serve as a general framework to excite static periodic-solitons and thereby program ordered localizations in mechanical metamaterials. Designing unit cells with multi-stabilities, for example, satisfying condition of KH 2 0 > 8(C 1 + C 2 ) and C s d 2 ≫ (K(H 0 − δ)δ − 2(C 1 + C 2 ))/4, the compressive load can trigger the displacive phase transitions in the metamaterials (see Fig. 5a,b in main text and Supplementary Fig. 9). Based on Eq. (S15), gradually increasing δ beyond the threshold δ s , the effective on-site potential P eff switches from monostable to bi-stable, and correspondingly, the deformation of the metamaterial undergoes a transition from uniform compression to ordered localization. This displacive phase transition accompanied by the spontaneous symmetry-breaking of the characteristic rotation (order parameter θ ) further induces the excitations of static kink, as schematically illustrated in Supplementary Fig. 9b,c. The topological constraint of the double well potential indicates each kink should necessarily be followed by an antikink in our metamaterial systems. Our experiments and simulations clearly demonstrate this evolution process (see Fig.   1 in main text, Supplementary Fig. 2 and Movie 1). To

Supplementary Movies
The