Phase Transitions in Hierarchical, Multi-stable Metamaterials

In this article, we consider the dynamics of transition waves in phase-transforming metamaterials with hierarchical architecture, i.e., 1D/2D periodic systems comprising a network of intersecting chains of elastically-coupled bi-stable elements. To this end, we develop continuum models of discrete 1D systems that, nevertheless, also elucidate the transition wave dynamics in 2D environments, which have received little attention in the literature. We find the potential driving and the wavelength relative to the hierarchical dimensions to play important roles in determining the wave mobility. The unique construction provokes some interesting results, including the growth of non-circular domains and the stabilization of domains of arbitrarily prescribed morphology; the latter representing an avenue toward reconfigurable performance via domain patterning. Altogether, in a break from the paradigm of homogeneity, the results not only elucidate the influence of hierarchy on the dynamics of phase-transforming metamaterials, but also its potential utility.


I. INTRODUCTION
Phase transformations 1,2 are a hallmark of materials whose microstructure possesses more than one stable equilibrium configuration (i.e., phase, state) as distinguished by a set of order parameters.In lowering the free energy of a sample, regions of homogeneous configuration (i.e., domains) emerge, separated by an interface (i.e., domain wall) that interpolates the adjoining phases and, in propagating, constitutes a transition wave.
The physical properties can vary drastically between domains of different phase and even within the domain walls, a quality exploited for a variety of applications 3,4 .Yet, transformation phenomena are not limited to natural materials systems and to the molecular scale; rather, phase-transforming mechanical metamaterials comprising arrays of coupled, multi-stable structural elements have been shown to not only mimic material-level processes at the structural level, but also serve as versatile platforms within which to engineer novel transformation characteristics 5 .
Prior investigation of phase-transforming metamaterials has revealed the peculiarities of domain wall motion within a variety of architectures, including those characterized by different substrates, spatio(-temporal) gradients, and hard defects [6][7][8][9][10][11][12][13] .These efforts elucidate the fundamental role of structure in shaping the observed transformation performance and provide a scholarly foundation upon which to develop design strategies toward the functionalization of transformation phenomena 9,[14][15][16][17][18][19][20][21] .However, prior examples of phase-transforming metamaterials exhibit only a single level of structure; yet, in the broader metamaterial literature, structural hierarchy 22 is already an intensely investigated and utilized characteristic of the architecture [23][24][25][26][27] .In the context of phase-transforming metamaterials, structural hierarchy not only represents a) Corresponding author email: mjfrazier@ucsd.edua means to expand the current design space within which to engineer the transformation performance, but also the opportunity to elicit new transformation behavior.
In the following, we consider the dynamics of transition waves in phase-transforming metamaterials with hierarchical architecture, i.e., 1D/2D periodic systems comprising a network of intersecting chains of elastically-coupled bi-stable elements.Within these systems, wave propagation becomes inherently inhomogeneous: the hierarchical structure functioning as soft material defects.The nature of the wave-structure interaction depends on the free energy wave driving and the wavelength relative to the hierarchical structure dimensions, affecting the propagation velocity.Unlike their homogeneous counterparts, in 2D hierarchical systems, a nucleus of the low-energy phase can grown into a domain with unusual, non-circular morphology or a prescribed domain pattern can be stabilized.Together, these results demonstrate the ability of hierarchical structure to influence the dynamics of phase-transforming metamaterials and to represent an additional design degree of freedom for applications.Applying both theoretical and numerical tools, we systematically study domain wall motion in and the characteristic domain morphology of representative 1D/2D hierarchical lattices.

A. Model
Figure 1a displays the schematic of a 1D hierarchical metamaterial which serves as our point of departure, comprising a main conduit and its periodically-distributed offshoots.
Each consists of bi-stable elements ascribed a single degree of freedom, u i,j (e.g., rotation, out-of-plane displacement), and coupled elastically to the nearest-neighbors.The elements are separated by a spacing, a, the fundamental spatial unit that also characterizes the offshoot length and periodicity, ℓ = pa and D = qa, p, q ∈ N, respectively.
We characterize the metamaterial as hierarchical due to the presence of multiple, potentially well-separated length scales: a, ℓ, and D.
The corresponding non-dimensional discrete governing equations with on-site viscosity, η, are given by [SI]   üi,0 + η ui,0 where δ Ii is the Kronecker delta and I = {i|i/D ∈ Z} collects the indices of junctions (i.e., elements in the main conduit to which the offshoots connect).Bi-stability stems from the on-site potential, ψ(u i,j ), with local minima at u i,j = u s1 and u i,j = u s2 , corresponding to the initial and final stable equilibrium phases, respectively; the local maxima at u i,j = u s0 is an unstable equilibrium.
The offshoots may be reguarded as soft defects having a purturbative affect on transisiton wave motion within an otherwise homogeneous system.Thus, in addition to numerical analyses of the discrete source model [Eq.( 1)], we also utilize the collective coordinate method in a semi-analytical approach to study the purturbed transition wave dynamics.To this end, we define an ansatz in the collective coordinates based on the transition wave solution of the continuum approximation of the homogeneous discrete system [SI]: where is a reduced variable with v 0 the steady-state velocity, For a given set of parameters, Fig. 1b.ii compares the numerical [Eq.( 1)] and analytical [Eq.( 3)] wave profiles where good agreement is observed, suggesting that results precipitating from the collective coordinate approach with Eq. ( 3) can be relevant to the discrete system.Here and in the following, we set η = 3/5 and k 0 = 1/10 while k 1 varies in order to manipulate ξ.In Fig. 1b.ii, k 1 = 4/50 yields ξ = 0.941.
In the following, we find the wavelength, λ to be a relevant parameter.We approximate λ as the span of the continuum solution centered at z = z 0 accounting for 95% of the total kinetic energy of the propagating wave.This criterion leads to the formulation [SI]:

B. Analysis
Before analyzing the 1D system with periodically-distributed offshoots as depicted in Fig. 1a, it is worth considering the effect of isolated offshoots (identically, D → ∞) on the domain wall motion.This choice simplifies the initial presentation of concepts and analytical procedures that extend to the more complicated systems to come.To this end, we excite a right-propagating transition wave in the discrete model with ξ = 0.941 (λ = 5.87) and, after attaining a steady-state velocity, observe its interaction with an isolated offshoot pair.Figure 2a offers a snapshot of the kinetic energy distribution in the vicinity of the junction when the propagation velocity is at a minimum.For ℓ < λ, the kinetic energy appears equally divided among the elements comprising the junction and offshoots as if their motion is that of a single unit; conversely, for ℓ > λ, the kinetic energy appears equally divided among the offshoots and conduit region beyond the junction as if the system is triplicated at the junction.From these observations, we formulate phenomenological descriptions of wave propagation within the main conduit, focusing on two limiting cases: (i) the short-offshoot condition (ℓ < λ) and (ii) the long-offshoot condition (ℓ > λ).
For ℓ < λ, we formulate a continuous governing equation that treats elements in the junction and offshoots collectively as a single, substitutional defect element with the summed local properties: where δ(x) denotes the Dirac-δ function and n is the number of offshoots per junction (presently, n = 2).For ℓ > λ, the two offshoots and the portion of the conduit beyond the junction respond identically, implying a force balance at the junction of the form u ,x | x=0 − = 3u ,x | x=0 + as if the junction was the interface between two different media.This inspires the governing equation: where H(x) is the Heaviside function.
To analyze transition wave motion governed by the Eqs.( 6) and ( 7), we follow a collective coordinate approach 28 using an ansatz base on Eq. (3): where X(t) and A(t) are the collective coordinates representing the transition wave position and shape, respectively.This leads to the governing equations in the collective coordinates [SI]: where q = [X A] T .The total Lagrangian, L, and generalized forces, Q X and Q A , are specific to the shortand long-offshoot models, Eqs. ( 6) and (7) [SI].Equation ( 9) is the reduced order model (ROM) of Eqs. ( 6) and (7), which are phenomenological descriptions of transition wave interaction with isolated offshoots in the discrete source model [Eq.( 1)].For validation, in Fig. 2, we compare the numerical results produced by the ROMs to that of the source model.We numerically integrate the ROMs from initial conditions, ] and q0 = [v 0 , 0], representing a transition wave approaching the junction from the left at steady-state velocity.We calculate the mean velocity, v m , in the region x ∈ [−2λ, 2λ] about the junction.As evidenced by the decrease in v m /v 0 with increasing ℓ displayed in Fig. 2b, the energy cost to the wave traversing the junction increases as ℓ/λ → 1.As a transition wave traverses a junction, it stimulates a duplicate wave in each of the offshoots; however, the duplicates cannot fully form until ℓ = ⌈λ⌉.As ℓ increases, more of the transition wave energy goes to creating the duplicates, decreasing v m /v 0 ; the energy and v m /v 0 stabilize once ℓ is sufficient for the duplicates to fully form.For ℓ > 4, the wave driven by ξ = 0.885 (λ = 5.98) can no longer overcome the energy barrier established by the offshoots and, therefore, is immobilized (i.e., pinned) at the junction.Conversely, for ξ = 0.941 (λ = 5.87), wave propagation persists for all ℓ; however, v m /v 0 rapidly converges to a constant.Apparently, v m depends on both ξ and ℓ/λ.Good agreement is observed between the numerical results from Eq. ( 1) and from the ROMs: the variation in v m /v 0 generated by the source model is well-approximated by the short-offshoot ROM (solid line) while ℓ < λ; by the long-offshoot ROM (dashed line) while ℓ > λ.   1), showing a clear bifurcation: if ξ is less than a critical value, ξ c , then the transition wave is pinned at the junction.The behavior can be anticipated from the time-independent solutions of Eqs. ( 6) and (7), respectively: While Eq. (10a) requires a root-finding algorithm to determine ξ c , the same can be determined analytically via Eq.(10b), giving (presently, n = 2 yields ξ c = 8/9).Apparently, these solutions plotted in Fig. 2c show excellent agreement with the simulation results.
For the case of the periodic hierarchical structure depicted in Fig. 1a, in addition to the offshoot length, ℓ, the offshoot period, D, affects the transition wave motion.In order to accommodate D, Eqs. ( 6) and (7), respectively, are modified as follows: where P (x; k) = (n + 1) k H(x − kD).Following the same collective coordinate approach, we determine the short-and long-offshoot ROMs for the periodic structure [SI].Because Eqs. ( 6) and ( 7) pertain to isolated offshoots, their periodic extensions [i.e., Eqs.(11)] and the corresponding ROMs assume D sufficiently large to approximate this condition, i.e., D > λ [SI].

III. 2D HIERARCHICAL METAMATERIAL
We also study the dynamics of transition waves propagating in 2D lattices formed from identical conduits of length, d = a(N − 1), where N is the number of constituent bi-stable elements.
While myriad (quasi-)crystalline lattices may be devised, the periodic square and hexagonal lattices (Fig. 3a,b) are amenable to analysis via only slight modification of Eqs.(11).Nevertheless, the dynamics observed in the simple square and hexagonal lattices extend qualitatively to more complex networks of 4-and 6-fold rotational symmetry [SI].The study proceeds along two tracks: (i) plane waves and (ii) domain growth from an initial nucleus.

A. Analysis: Propagation of Plane Waves
From the parent 2D lattice, we can isolate a 1D hierarchical structure whose dynamics are representative of plane wave propagation along the lattice symmetry directions.
Through these isolated structures, the previous 1D analysis can be extended to the 2D setting.
In particular, for horizontal plane wave propagation in the square lattice, we envision the vertical conduits bisected by horizontal axes (Fig. 3a.i).For even N , the two elements on either side of an axis displace synchronously and, thus, the interaction force between them is zero.Accordingly, we isolate the 1D structure situated between two adjacent bisecting axes, which resembles the hierarchical system in Fig. 1a; therefore, for d > λ which corresponds to ℓ ≥ λ and D = d, the ROM developed from Eq. (11b) is expected to describe transition waves propagating within the structure.As the 2D system emerges from appending copies of the 1D structure in the perpendicular direction, the results of the analysis are expected to be pertinent to horizontal plane wave propagation.Similarly, for diagonal plane wave propagation, we envision the junctions bisected by diagonal axes (Fig. 3a.ii).In isolating the 1D structure between adjacent axes, the bisected junction behaves as a defect element with local properties half that of the background.Consequently, the ROM developed from Eq. (11a), with nℓ replaced by −1/2, is expected to apply.
To isolate 1D hierarchical structures from the 2D hexagonal lattice, we follow a similar procedure as above.For horizontal plane wave propagation, we bisect the vertical conduits and isolated the 1D structure between two adjacent bisections, whose junctions are host to n = 1 offshoot (Fig. 3b.i); thus, we expect the relevant ROM to be that emerging from Eqs. (11b) with n = 1.For "diagonal" plane wave propagation, we bisect the conduits linking A-and B-type junctions (Fig. 3b.ii).The bisecting axes evenly splits diagonal conduits while leaving intact conduits which subtend consecutive bisections.Therefore, the local properties of elements in split conduits are effectively half those in subtending conduits.The ROM for this unique scenario emerges from: where To validate our approach, we simulate the propagation of a plane transition wave along the horizontal/diagonal axes of the square and hexagonal lattices, and then compare these results to those of the ROMs inspired by the isolated 1D hierarchical structures (Figs.3c,d).The simulated discrete lattices (source models) have the approximate dimensions, W × L = 35d × 35d, with constituent chains of length, d = 15.Periodic boundary conditions are applied along the edges perpendicular to the propagation direction (i.e., the L dimension).A right-propagating plane wave driven by a range of ξ is initiated by prescribing (u, u) = (u s2 , 0) to elements at the left boundary (i.e., along the W dimension).As the effective wavefront propagates in the either the horizontal or diagonal directions, the path taken by constituent transition waves within the lattice conduits is less direct.The mean velocity is plotted in Figs.3c,d is calculated from effective wavefront (equiv., by projecting the velocity of the constituent waves onto an axis parallel to the horizontal/diagonal direction.).
In general, excellent agreement is observed in Figs.3c,d, validating the ROMs and, by extension, the proposed method of extracting a 1D hierarchical structure from a 2D parent lattice for the purpose of analyzing plane transition wave propagation in the parent lattice.In Fig. 3c, the profile of v m (ξ) as determined by the ROM for horizontal propagation is supported by the simulation results; however, the critical ξ separating the propagating and non-propagating regimes differs between the ROM and source model.We attribute this to the enhanced stiffness of the continuum model underlying the ROM compared to the discrete system, exacting a higher energy penalty as the wave traverses a junction, i.e., requiring a higher ξ to overcome.
These results are the first demonstration of direction-dependent propagation for transition waves in a 2D phase-transforming metamaterial.Interestingly, as indicated by the intersection of horizontal and diagonal results in Figs.3c,d, adjusting ξ enables the exchange of the favored propagation direction (i.e., that in which the propagating wave is least impeded by the hierarchical structure).On one hand, compared to effective propagation horizontal direction, the hierarchical architecture establishes a diminished energy barrier for constituent waves propagating in the lattice conduits.On the other hand, compared to effective propagation horizontal direction, the propagation of constituent waves in the lattice conduits is less aligned with the diagonal direction, diminishing the effective propagation velocity.The exchange occurs as one of these competing influences dominates the other.

B. Morphology of Closed Domains
In the following, we numerically investigate the morphology of closed domains that evolve from a single, centrally-located nucleating element within 2D lattices.For the 2D system initially uniform in the meta-stable phase, u s1 , a small nucleus of the less energetic phase, u s2 , will trigger a phase transformation accommodated by the propagation of transition waves forming the expanding contour of a closed domain uniform in u s2 .Figure 4a illustrates three, ξ-dependent morphologies that the nucleated domain may acquire as it evolves within a square hierarchical lattice; indeed, any hierarchical lattice with 4-fold rotation symmetry: an octagon and square with one of two orientations with respect to the host lattice.
Similarly, as exemplified by the hexagonal hierarchical lattice in Fig. 4b, within a hierarchical lattice of 6-fold rotation symmetry, the domain morphology approximates either a dodecagon or a hexagon with one of two orientations with respect to the underlying lattice.
Conversely, in structurally homogeneous systems (or hierarchical systems with d/λ → 0), the wave velocity is (near-)independent of the propagation direction; therefore, the domain acquires a circular morphology (Fig. 4c).
Although the theoretical results and discussion of earlier sections are most pertinent to the propagation of transition waves along a single dimension, Figs.3c,d, nevertheless, aid in justifying the morphology of the closed domains in 2D hierarchical lattices.While the simulations utilize a single (point-like) nucleating element, the justifications are most straightforwardly provided by considering the evolution of a circular nucleus.For ξ c ≤ ξ < ξ x , Figs. 3c,d show that plane transition waves propagate along the diagonal axis at a greater velocity than along the horizontal axis.Considering the outward normals to the contour of the circular nucleus, as the system evolves, it may be expected that the portion of the contour with normals most aligned with the diagonal (horizontal) axis correlates to the fastest (slowest) expanding region of the nucleus.The difference in velocity deforms the nucleus; however, no portion of the contour can become concave since this would produce normals directed along the diagonal and the wavefront would quickly fill the concavity.
Therefore, within the square (hexagonal) lattice, the circular nucleus grows into a domain with square (hexagonal) morphology with symmetry axes aligned with those of the square (hexagonal) lattice.On the other hand, for ξ > ξ x , Figs. 3c,d now show the velocity to be greatest along the horizontal axis.Therefore, similar reasoning involving the outward normals to the contour of the circular nucleus suggests that the nucleus, once again, develops a square (hexagonal) morphology within a FIG. 5. Hierarchy and Domain Stability.For a sufficiently large potential difference, the low-energy domain will expand into and overtake the meta-stable domain, yielding a system of homogeneous phase.Below a critical value of ξ, except for small changes along the perimeter to balance forces at the junctions, an arbitrary initial phase distribution will not evolve in time.
square (hexagonal) lattice; although, the axes of the domain are rotated by π/4 (π/6) with respect to the underlying square (hexagonal) lattice.Finally, for ξ ≈ ξ x , Figs. 3c,d indicate propagation of near-equal velocity along the horizontal and diagonal axes.Applying our earlier reasoning, we expect the circular nucleus to grow into a domain with octagonal (dodecagonal) morphology within the square (hexagonal) lattice.Apparently, by changing ξ, it is possible to manipulate the morphology of nucleated domain.
As mentioned above in the context of conventional phase-transforming materials, the physical properties can vary drastically between domains of different phase, which may be exploited for applications.Although the literature contains examples of metamaterials with configuration-specific properties, only a few of these (restricted to 1D architectures) allow multiple domains to stably exist simultaneously 18,[29][30][31] ; otherwise, the excess free energy drives the system to a uniform low-energy phase.The ability to stabilize 2D domains of arbitrary morphology raises the prospect of a custom and re-configurable 2D domain patterning which, may open the door to functionalities and applications inherently out of reach of 1D systems.Here, for the first time, we demonstrate a stable 2D domain of complex morphology utilizing a hierarchical 2D lattice with ξ < ξ c (Fig. 5).
In the SI, an impromptu waveguide is constructed from a prescribed domain pattern.

IV. CONCLUDING REMARKS
In summary, we have introduced the concept of structural hierarchy to phase-transforming metamaterials through 1D/2D networks of intersecting chains of bi-stable elements and developed a framework for investigating the consequences thereof for transition wave propagation.For the purpose of analysis, the hierarchical structure is interpreted as a material heterogeneity in continuum approximations of discrete source models, which are shown to be predictive of the dynamic response.Physically, the hierarchical structure generates an energy barrier that either impedes or, for insufficient potential driving, halts the propagation.Interestingly, an understanding of transition wave propagation within 1D systems is shown to extend to 2D periodic systems along the symmetry axes.In particular, we are able to use the 1D results to justify the unusual domain morphologies that uniquely develop within 2D hierarchical systems from a smaller nucleus.Moreover, if the potential driving is insufficient for propagation, then any prescribed complex morphology is rendered stable, which with configuration-specific physical properties supports the notion of domain patterning for custom performance in potential, e.g., soft matter applications.Altogether, the results not only elucidate the influence of hierarchical structure -a new design degree of freedom -on the dynamics of phase-transforming metamaterials, but also its potential utility.

FIG. 1 .
FIG. 1. Hierarchical, Phase-transforming Metamaterial.(a) Schematic of a 1D periodic hierarchical system.(b) On-site non-convex potential characterized by two stable states, us1 and us2, and granting (ii) the formation and propagation of an anti-kink transition wave.The wave profile as determined numerically and analytically, respectively, from the homogeneous discrete model [Eq.(1)] and its continuum approximation [Eq.(3)].

FIG. 2 .
FIG. 2. Propagation in 1D System with Isolated Offshoots.(a) The kinetic energy distribution in the vicinity of the offshoot pair as the transition wave traverses the junction.(i) Elements comprising short offshoots move in tandem with the junction; (ii) new, out-going fronts are identically stimulated in long offshoots and the conduit region beyond the junction.(b) Comparison of the numerical results for the (ℓ, ξ)-dependent vm/v0 from the source model and the ROMs.(c) Mapping ξ and ℓ/λ to the vm.The critical potential difference, ξc, as predicted by the ROMs for the transition wave to successfully traverse the junction (i.e., not become immobilized).

Figure
Figure2cplots v m (ℓ, ξ) as determined from the simulation of Eq. (1), showing a clear bifurcation: if ξ is less than a critical value, ξ c , then the transition wave is pinned at the junction.The behavior can be anticipated from the time-independent solutions of Eqs.(6) and(7), respectively:

FIG. 3 .
FIG. 3. Plane Wave Propagation in 2D Hierarchical Lattices.Schematics of (a) square and (b) hexagonal lattices comprising intersecting chains of bi-stable elements.The lattice are oriented for a transition wave propagating (shaded region with arrow) along the designated horizontal and diagonal directions.1D hierarchical structures are isolated between two bisecting axes.(c,d) Comparison of numerical results for mean plane wave velocity as determined from the discrete 2D source model and the ROMs developed from the isolated structures.According to the ROMs, for ξ ≥ ξc, wave propagation is supported in both the horizontal and diagonal directions; for ξ = ξx, the propagation velocity is identical in both directions.

FIG. 4 .
FIG. 4. Hierarchy and Domain Morphology.Simulation snapshots of the closed domain that evolves from a nucleus of the low-energy phase.(a) Within a square hierarchical lattice, the nucleus grows into a domain with 4-or 8-fold rotational symmetry.(b) Within a hexagonal hierarchical lattice, the nucleus becomes a domain with 6-or 12-fold rotational symmetry.(a,b) The potential difference selects for the principal direction (horizontal or diagonal) along which transition wave propagation is favored, yielding a domain morphology with the symmetry of the underlying lattice.If both propagation directions are equally-favored, then the domain morphology develops a rotational symmetry twice that of the underlying lattice.(c) The nucleus initiated within (i) a homogeneous system or (ii) a hierarchical lattice with d/λ ≈ 1/3 will develop into a circular domain.