Wave coarsening drives time crystallization in active solids

TL;DR

This paper uncovers a novel wave coarsening phenomenon in active elastic media, leading to the formation of time crystals through spontaneous growth of vibrations, highlighting the role of symmetries in non-equilibrium systems.

Contribution

It introduces wave coarsening as an inertial mechanism for time crystal formation in active solids, expanding understanding of non-equilibrium pattern dynamics.

Findings

Wave coarsening observed in 1D and 2D active solids

Dynamical scaling laws for wave coarsening established

Breaking momentum conservation reveals topological edge modes

Abstract

When metals are magnetized, emulsions phase separate, or galaxies cluster, domain walls and patterns form and irremediably coarsen over time. Such coarsening is universally driven by diffusive relaxation toward equilibrium. Here, we discover an inertial counterpart - wave coarsening - in active elastic media, where vibrations emerge and spontaneously grow in wavelength, period, and amplitude, before a globally synchronized state called a time crystal forms. We observe wave coarsening in one- and two-dimensional solids and capture its dynamical scaling. We further arrest the process by breaking momentum conservation and reveal a far-from-equilibrium nonlinear analogue to chiral topological edge modes. Our work unveils the crucial role of symmetries in the formation of time crystals and opens avenues for the control of nonlinear vibrations in active materials.

Equations

$$\ddot{\varphi }={\nabla }^2\left[(1-i\alpha )\varphi +\dot{\varphi }+|\varphi {|}^2\varphi \right].$$

$$\ddot{h}=-\partial^{2}_{x}\left[\partial^{2}_{x}h-\alpha\Bigg{(}\frac{\partial^{3}_{x}h}{\sqrt{1+(\partial^{3}_{x}h)^{2}}}\Bigg{)}+\partial^{2}_{x}\dot{h}\right].$$

$$\frac{\dot{|a_n|}}{|a_n|}=\frac{1}{2}\left(\frac{\alpha }{\sqrt{1+U(t)}}{q}_n^p-q_n^{2p}\right).$$

$${\ddot{a}}_n(t)$$

$${\dot{\sigma }}_n^2+{\ddot{\sigma }}_n=-q_n^2\left((1-i\alpha )+{\dot{\sigma }}_n+2\sum _m|a_m{|}^2\right).$$

$${\dot{\sigma }}_n={\dot{⟨\sigma ⟩}}_n+{\dot{\delta }}_n$$

$${\dot{⟨\sigma ⟩}}_n^2=$$

$${\ddot{⟨\sigma ⟩}}_n+\ddot{{\delta }_n}+2{\dot{⟨\sigma ⟩}}_n\dot{{\delta }_n}+{\dot{{\delta }_n}}^2=0.$$

$${\ddot{\delta }}_n+2{\dot{⟨\sigma ⟩}}_n{\dot{\delta }}_n$$

$${\dot{⟨\sigma ⟩}}_n$$

$${\dot{\delta }}_n$$

$${\dot{\sigma }}_n\approx -\frac{q_n^2}{2}+i{q}_n\sqrt{1+2\sum _m|a_m{|}^2}+\frac{\alpha q_n}{2\sqrt{1+2{\sum }_m|a_m{|}^2}}.$$

$${\omega }_n(t)=q_n\sqrt{1+2\sum _m|a_m(t){|}^2},$$

$$\dot{|a_n|}=\frac{1}{2}\left(\frac{\alpha }{\sqrt{1+2{\sum }_m|a_m{|}^2}}{q}_n-q_n^2\right)|a_n|,$$

$$\frac{d}{dt}(|a_n{|}^2)=\left(\frac{\alpha }{\sqrt{1+2{\sum }_m|a_m{|}^2}}{q}_n-q_n^2\right)|a_n{|}^2.$$

$$2\sum _m|a_m{|}^2$$

$$2\sum _m{q}_m|a_m{|}^2$$

$$2\sum _m{q}_m^2|a_m{|}^2$$

$$\dot{A}=\alpha Q(t)A^{1/2}(t)-Q(t{)}^{2}A(t).$$

$$t^{\chi -1}\sim t^{2\beta +\chi }+t^{\beta +\frac{1}{2}\chi },$$

$$⟨q⟩\sim t^{\beta }=t^{-1/2}$$

$$\sqrt{\sum |a_m{|}^2}\sim t^{\gamma /2}=t^{1/2}.$$

$$m{\partial }_t^{2}h+a^2\kappa {\partial }_x^{4}h-2a^3{\kappa }^a{\partial }_x^{5}h+a^2\gamma {\partial }_x^4{\partial }_{t}h=0,$$

$$\ddot{h}=-{\partial }_x^2\left[{\partial }_x^{2}h-\alpha {\partial }_x^{3}h+{\partial }_x^2\dot{h}\right],0\le x\le L,$$

$$\partial_{x}^{5}h\rightarrow\partial_{x}^{2}\Bigg{(}\frac{\partial_{x}^{3}h}{\sqrt{1+(\partial_{x}^{3}h)^{2}}}\Bigg{)},$$

$$\ddot{h}=-\partial^{2}_{x}\left[\partial^{2}_{x}h-\alpha\Bigg{(}\frac{\partial^{3}_{x}h}{\sqrt{1+(\partial^{3}_{x}h)^{2}}}\Bigg{)}+\partial^{2}_{x}\dot{h}\right].$$

$$\displaystyle\ddot{a}_{n}=-q_{n}^{4}\,a_{n}-q_{n}^{4}\dot{a}_{n}-\alpha\,q_{n}^{2}\widehat{\left(\frac{\partial_{x}^{3}h}{\sqrt{1+(\partial_{x}^{3}h)^{2}}}\right)}\Bigl{.}_{\!n}$$

$${\ddot{a}}_n=-q_n^4{a}_n-q_n^4{\dot{a}}_n+\alpha i{q}_n^5\frac{a_n}{\sqrt{1+2{\sum }_m{q}_m^6|a_m{|}^2}}.$$

$${\dot{\sigma }}_n^2+{\ddot{\sigma }}_n=-q_n^4\left(1+{\dot{\sigma }}_n-i{q}_n\frac{\alpha }{\sqrt{1+2{\sum }_m{q}_m^6|a_m{|}^2}}\right).$$

$${\dot{\sigma }}_n=-\frac{q_n^4}{2}+\frac{q_n^2}{2}\sqrt{q_n^4-4\left(1-i\alpha \frac{q_n}{\sqrt{1+2{\sum }_m{q}_m^6|a_m{|}^2}}\right)}.$$

$${\dot{\sigma }}_n=-\frac{q_n^4}{2}+i{q}_n^2+\frac{q_n^2}{2}\alpha \frac{q_n}{\sqrt{1+2{\sum }_m{q}_m^6|a_m{|}^2}}.$$

Full text

A striking manifestation of nonlinear physics is when microscopic periodic dynamics synchronize across space to produce coherent oscillations. Examples range from Faraday waves [1, 2] to parametrically pumped optomechanical lattices [3] and trapped ion chains [4]—which exhibit subharmonic oscillations arising from spontaneously broken discrete time-translation symmetry. More recently, analogous phenomena have been explored with broken continuous time-translation symmetry, such as fireflies lighting up in synchrony [5, 6], cold atomic gases oscillating under optical pumping [7, 8], as well as biological collectives [9, 10, 11, 12, 13] and metamaterials [14, 15, 16]—where self-sustained oscillations emerge without any periodic drive. These are called time crystals.

Time crystals are typically driven externally—for example, by accelerating a water tank to make Faraday waves [1], or forcing atoms, ions or resonators using light [4, 7, 8, 15, 17]. This external drive injects both energy and, crucially, momentum. An exception is active matter, which is energized from within. Active matter with elastic interactions has recently been shown to host time-crystalline phases [9, 10, 11, 13, 14, 16]. In this context, a key paradigm is non-reciprocal elasticity [18, 19], in which energy is injected through an elastic tensor which is not symmetric  [11, 20, 21, 22, 23, 16, 24]. Unlike traditional routes to time crystallization, these odd elastic solids inject energy whilst respecting momentum conservation. Here we ask: what is the role of conservation laws in the formation of time crystals?

We construct odd elastic metamaterials and discover a time-crystallization phenomenon that we dub wave coarsening: active elastic waves spontaneously emerge, and their wavelength, period, and amplitude all grow in a self-similar fashion. This process arises from a nonlinear coupling between modes, sustained by active internal stresses. These stresses are generated through the injection of energy at the material’s microscopic scale and cascade to ever larger scales. Wave coarsening fundamentally differs from the traditional synchronization towards time crystals via phase locking. It is also distinct from traditional coarsening phenomena that relax towards equilibrium via diffusion. Our results open new questions about the transient physics of systems with non-potential interactions [25, 15, 26, 27] and suggest an organizing principle for nonlinear waves in acoustics [28, 29], optomechanics [30, 31], living matter [11, 32], and soft robotics [16] in the presence of conservation laws [33, 34].

Our active solids are mechatronic honeycombs (Fig. LABEL:COARSENING:trihex). Each unit cell is made from a trio of elastic hexagons, whose strains are communicated locally to a custom printed circuit board which implements a nonreciprocal and non-pairwise relation between the deformations sensed and the forces actuated (SM §I). The active forces $F_{i}$ on each site $i$ depend on the difference in deformation $\delta L_{j}-\delta L_{k}$ between two neighboring sites $j$ and $k$ via

where $k^{a}$ sets the strength of the non-reciprocal interaction. As a consequence, the unit cell breaks Maxwell-Betti reciprocity and responds asymmetrically regardless of which hexagon is compressed (Fig. 1a), a property not allowed in absence of non-potential forces [18].

When $k^{a}$ is driven beyond a threshold value, the unit cell undergoes a Hopf bifurcation and exhibits a limit cycle (Fig. 1b). We now examine how this instability unfolds at the collective scale by arranging 37 unit cells into a honeycomb lattice (Fig. 1c). With the material placed on top of a glass plate, lubricated with polytetrafluoroethylene (PTFE) to minimize friction, we now drive the nonreciprocal interactions beyond the bifurcation. Upon nudging the material out of equilibrium, we observe amplifying self-oscillations of around 30Hz at the level of each unit cell, which we measure and represent as a strain field in Fig. 1de. As the material oscillates, a longer wavelength mode slowly amplifies and drowns the initial unit cell oscillations. This wave travels unidirectionally along the edge before the system eventually settles into a globally synchronized oscillation with stable frequency, characteristic of a time-crystalline state. We refer to this transient process which is marked by the simultaneous growth of wavelength, period and amplitude, as wave coarsening.

To quantify the evolution of strain patterns, we track the inverse correlation length $\xi^{-1}$ of strains over time. Averaged over space and across multiple runs, these data confirm that strain waves coarsen in wavelength over time (Fig. 2a). Moreover, the rate at which the process unfolds suggests a power-law behavior. However, in tandem with the growth of wavelength, the strain amplitudes also steadily grow (Fig. 2b), in contrast to gradient flow domain coarsening and synchronization phenomena.

To gain insight into the coarsening mechanism and its scaling behavior, we coarse-grain microscopic displacements at the unit cell level into a continuum field (Fig. LABEL:COARSENING:TriHexSchematic, SM §LABEL:Coarsening:coarsegrainingTrihex):

[TABLE]

Equation (1) describes the nonlinear vibrations of a dissipative medium that is powered by odd elasticity, where shear modes couple antisymmetrically. The complex field $\phi(x,t)$ is the elastic strain field, whose real and imaginary parts $S_{1}$ and $S_{2}$ correspond to the two linearly independent shear modes. The parameter $\alpha$ is the non-dimensional activity that relates to the microscopic non-reciprocity $k^{a}$ and hinge stiffness $\kappa$ of our plaquettes as $\alpha=\sqrt{3}/8k^{a}/\kappa$. Because energy is injected whilst conserving linear and angular momentum, the material does not undergo uniform translations or rotations but instead deforms.

These conservation laws are encoded in the fully Laplacian form of Eq. (1), whose dispersion features gapless modes of wavenumber $q$ with complex frequencies $\sigma\approx\frac{1}{2}q(\alpha-q)+iq$ (Figs. LABEL:fig:PhiDispersion LABEL:fig:SecondaryStability, SM §LABEL:COARSENING:PhiInstabilities). All modes oscillate at a frequency proportional to their wavenumber $q$, while the non-reciprocal term amplifies modes with $q<\alpha$, with maximal growth at $q^{*}=\alpha/2$. Amplification at $q^{*}$ aligns with our observation that in the initial stage of wave coarsening traveling patterns with a characteristic wavenumber emerge before nonlinear effects saturate growth.

To confirm that waves coarsen in our continuum model, we investigate its nonlinear dynamics by numerically integrating Eq. (1) with random initial conditions, on both one-dimensional and two-dimensional periodic domains. These simulations show that coherent wave patterns rapidly emerge and steadily coarsen in wavelength (Fig. 2c) with a characteristic power-law scaling while amplitudes steadily grow (Fig. 2d), regardless of the dimension (Fig. S5), mirroring our experimental results.

Our results suggest that wave coarsening is a robust route to time crystallization, powered by momentum conserving energy injection. To test this hypothesis, we now consider a distinct non-reciprocal platform with a different dispersion and non-linearity.

We construct a 1D periodic metamaterial in the form of an elastic ring embedded with 24 torque motors that respond antisymmetrically to angular deviations in neighboring sites (Fig. 3a, SM §VII). We place the ring upon an air table to let it float frictionlessly and drive the non-reciprocal coupling strength beyond its threshold of instability. Any small perturbation renders the ring unstable, causing waves with high mode number $n$ to emerge and grow in amplitude until eventually the motors saturate and deviate from their linearly non-reciprocal response. These waves then coalesce into longer wavelengths, deforming the ring into a series of unidirectionally morphing polygonal shapes with decreasing phase velocity and increasing amplitude (Fig. 3b).

Like the two-dimensional metamaterials described above, wavelength, period and amplitude all increase during the coarsening process. Yet, in our one-dimensional ring, the time evolution of the mode structure is clearer (Fig. 3c). While at early times, the $n=6$ mode quickly amplifies and briefly dominates the spectrum, soon the next highest mode $n=5$ catches up and overtakes. This victory for $n$th mode is again short-lived as it rapidly decays once its lower wavenumber neighbor $n-1$ takes over. As mode amplitudes transfer and mix, a travelling phase slip can be observed (Fig. 3b). Eventually, the ring finds its dynamical steady state on a slow cycle around the lowest accessible deformation mode in the shape of a pancake.

During the coarsening process, the wave velocity $\omega_{n}/q_{n}$, where $q_{n}=2\pi n/L$, consistently decreases proportionally to the mode number $n$ (Fig. 3d). This suggests the existence of a quadratic dispersion characteristic of flexural waves that is still valid in the nonlinear regime and coincides with a characteristic amplitude for each wavelength. Finally, the inverse correlation length $\xi^{-1}$ plotted over time reveals another power-law trend, which we fit at $\xi^{-1}\propto t^{-0.29}$ (Fig. 3e).

In order to consolidate these findings, we derive a coarse-grained description of a 1D non-reciprocal ring with a saturating nonlinearity, describing its real-valued deviations from its equilibrium state as a height displacement field $h(x,t)$ (SM §VII):

[TABLE]

Analogously to the 2D model Eq. (1), the right hand side of Eq. (2) is entirely inside the Laplacian and thus conserves momentum whilst breaking reciprocity. However, in this beam equation, the flexural deformations are inherently dispersive, and the nonlinearity takes a different form. To test whether coarsening persists under these conditions, we simulate Eq. (2) and observe a clear power-law growth of the inverse correlation length $\xi^{-1}$, with an exponent close to $-1/4$, which is consistent with our experimental findings (Fig. 3e). In summary, despite the differences in dimension, dispersion and non-linearity, waves in our 1D non-reciprocal ring grow in both wavelength and amplitude in a similar fashion to the 2D odd elastic metamaterial of Fig. 1, yet with a different dynamical exponent.

The data of Fig. 3 strongly suggests a picture of various modes interacting with one another but the precise nature of this coupling remains unclear. To gain analytical insight, we write our continuum theories of Eqs. (1), (2) in the Fourier domain and make a mean field approximation for the interaction between modes (Figs. S6, S7, S8). Further averaging over fast oscillations in favor of the slow coarsening dynamics allows us to derive an evolution equation for the mode amplitudes which remains valid at the nonlinear level:

[TABLE]

Here $|a_{n}|$ and $\dot{|a_{n}|}$ denote the mode amplitudes and their growth rate, $p=1$ for the odd elastic system of Eq. (1) and $p=2$ for the non-reciprocal beam of Eq. (2). Crucially, $U(t)=2\sum_{m}{q_{m}}^{6(p-1)}|a_{m}|^{2}$ represents an amplitude-dependent rescaling of the activity $\alpha$ that invites a compelling picture.

At the start of the coarsening process, modes grow according to their linear dispersion, and the most unstable wavelength is selected by $\alpha$. However, as modes inevitably amplify, eventually the nonlinear term $U(t)$ grows to become comparable to the linear terms, which shifts the peak of the dispersion to a lower wavenumber and growth rate simultaneously (Fig. S7). This new wavenumber now grows, albeit at a slower rate, and further increases $U(t)$, in turn pushing the dominant wavenumber and growth rate even lower. A feedback loop ensues until a finite system size prevents a further cascade of wavenumbers. The resulting progressive shift in the dispersion curve suggests a self-similar dynamics: as energy builds up in the system, non-reciprocal forces drive growth at ever-larger length scales and at ever-slower rates.

To understand the timescales at which wave coarsening occurs, we search for solutions to Eq. (3). Motivated by the experimental observation that even beyond the linear instability, waves remain quite coherent with a narrowly distributed mode structure, we use as as ansatz a wavepacket whose central wavenumber is larger than its width (SM §VI). This yields an evolution equation for the dominant length scale. Then, using a power-law ansatz $\xi^{-1}\propto t^{\beta}$ for the inverse correlation length, we find power laws $\beta=-1/(2p)=-1/4$ for Eq. (2), confirming the experimental and numerical results of Fig. 3e. In the 1D case of the nonlinearly odd elastic model of Eq. 1, we similarly find $\beta=-1/(2p)=-1/2$, which is in agreement with numerical data (Fig. 2c). It is remarkable that our mean-field picture successfully captures the coarsening dynamics despite neglecting the details of the phase slip events. However, our approximation breaks down in 2D, where faster coarsening occurs through a sequence of strain defect pair annihilation events. We interpret the difference in coarsening rates as follows. While in 1D, wave coarsening is driven by mode competition as described by Eq. (3), with phase slips occurring only sporadically, the 2D counterpart is dominated the dynamics of mobile and attracting defects, which accelerates coarsening.

Similar phase slips and defect annihilations occur in a plethora of coarsening and synchronizing systems [35, 36, 37, 38, 33, 39, 40]. Yet in these cases, the coarsening dynamics is effectively a gradient flow that minimizes a global potential. By contrast, here the patterns and defects observed are nonlinear excitations that can only be sustained by continuous energy injection. Instead of relaxing toward a fixed potential minimum, the system redistributes energy as the nonlinear dispersion relation progressively favors larger length scales.

Now that we have captured the essence of wave coarsening, we ask in the remainder of the paper how robust it is to perturbations. We consider two perturbations: (i) breaking of momentum conservation and (ii) edge effects.

First, we break momentum conservation. Our mean field analysis showed that wave coarsening continues until the wavelengths become system-sized, ultimately due to the fact that the linear growth rate $\real\sigma\propto q(\alpha-q)$ is positive for all wavenumbers $q<\alpha$ , a consequence of the momentum conserving nature of the energy injection. We now gap these low $q$ modes by breaking momentum conservation, which can be achieved by enforcing momentum exchange of the material with a substrate, i.e. through friction.

We now remove the PTFE lubrication between the 2D material of Fig. 1 and the substrate and observe how wave coarsening is arrested indefinitely while the system cycles around higher mode solutions without cascading further (Fig. 4a). To quantify this behavior, we numerically integrate Eq. (1) with an added friction term proportional to $\Gamma\dot{\phi}$, and find that the steady state correlation length $\xi^{-1}$ perfectly coincides with the marginally stable wavelength predicted by our stability analysis (Fig. 4b and Fig LABEL:fig:SecondaryStability).

Second, we consider the effect of boundary conditions. To this end, we employ large-scale discrete simulations of the 2D material. Starting with a disk consisting of $5813$ unit cells, these reveal that while coarsening initially proceeds homogeneously throughout the bulk, chiral edge waves soon appear and eventually dominate the dynamics. These are Rayleigh waves that are unidirectionally amplified by non-reciprocity and stabilised by nonlinearities (Fig. S9, SM §VIII). As their amplitude grows, so does their wavelength and penetration depth, causing them to progressively invade the bulk (Fig. 4c). At early times, coarsening is driven by the annihilation of free-roaming strain defects of opposite charge. Once edge waves emerge, however, they impose spatial order on the defects and confine them to be regularly spaced. Coarsening then proceeds through discrete events reminiscent of the phase slips seen in Fig. 2, as defects are gradually expelled at the system’s edge, each time increasing the edge wavelength.

To confirm the crucial role of edges in our system, we repeat our simulations on an annular domain and observe edge waves appear on both the inner and outer edge and propagate in opposite directions (Fig. 4d). As their amplitude and wavelength grows, these edge waves eventually meet and the outer edge overwhelms the inner one due its larger circumference. Their counterpropagating nature is reminiscent of chiral edge states in Chern insulators [41, 42, 43, 44, 45, 46, 47]. However, because these waves progressively invade the bulk we do not yet have an insulating state.

Inspired by this analogy we again include friction, drawing on the finding that it selects a steady-state length scale and arrests waves from coarsening. Since the penetration depth of Rayleigh waves is proportional to their wavelength (SM §VIII), we find that the system now maintains an insulating state with counter-propagating edges waves, representing a striking nonlinear, far-from-equilibrium counterpart to the quantum Hall effect (Fig. 4e). This demonstrates how wave coarsening irremediably occurs even in the presence of perturbations, but the final time-crystalline state is sensitive to both the breaking of momentum conservation and to boundary conditions.

In conclusion, we have uncovered a nonlinear feedback mechanism in active elastic media, where momentum-conserving energy injection drives waves to ever larger wavelength, amplitude and period. Looking ahead, we anticipate that wave coarsening represents a paradigmatic route towards time-crystallization in systems with additional symmetries and associated Goldstone modes, that can be leveraged as a control strategy for shape changes and directed motion in fields such as robotics, active matter and optomechanics with non-reciprocal phase transitions [25, 48].

Acknowledgments

We thank Tjeerd Weijers, Kasper van Nieuwland, Sven Koot, Melanie Pichon, for technical assistance and Colin Scheibner, Michel Fruchart and Gustavo Düring for valuable discussions. C.C. and J.V. acknowledge funding from the European Research Council under Grant Agreement No. 852587 and from the Netherlands Organisation for Scientific Research (NWO) under grant agreement VI.Vidi.213.131.3. J.B. acknowledges funding from the European Union’s Horizon research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 101106500. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP) and by the Dutch Institute for Emergent Phenomena (DIEP).

Supplementary Materials for

Jonas Veenstra1, Jack Binysh1, Vito Seinen1, Rutger Naber1, Damien Robledo Poisson1, Andres Hunt2, Wim van Saarloos3, Anton Souslov4, Corentin Coulais1 1Institute of Physics, Universiteit van Amsterdam, 1098 XH Amsterdam, The Netherlands 2Instituut-Lorentz, Universiteit Leiden, 2300 RA Leiden , The Netherlands 3TCM Group, Cavendish Laboratory, JJ Thomson Avenue, Cambridge, CB3 0US, United Kingdom

Contents

1. I Materials and methods

2. Construction of the two-dimensional active metamaterial

3. V Observation of coarsening dynamics in the $\phi$ model

4. VI Mean field approximation to the $\phi$ model

5. Reduction to the mean field model

6. Reduction from second to first order

7. Numerical verification of the second to first order reduction

8. VII The nonreciprocal beam and its coarsening exponents

9. Mean field and coarsening of the beam model

10. VIII Unidirectionally amplified Rayleigh waves in odd elastic media

I Materials and methods

Construction of the two-dimensional active metamaterial

The metamaterial shown in Fig. 1 consists of 3D-printed hexagonal structures (Formlabs Form 2 SLA printer with Tough1500 resin) with an embedded compliant flexure (E3D ToolChanger multi nozzle fused filament fabrication) that houses a custom-made voice coil actuator, as shown in Fig. LABEL:COARSENING:trihexa. The actuator consists of a weak iron cylindrical core attached on both ends to neodymium magnets aligned to repel (Fig. LABEL:COARSENING:trihexb). The resulting magnetic field lines meet in the middle and bend radially outwards perpendicularly to a copper coil that is wound around the weak iron core. The radius of the copper coil is larger than the radius of the iron core, which allows the coil to translate along the core (Fig. LABEL:COARSENING:trihexc).

When an electrical current is turned on, a Lorentz force is generated that pushes the iron core in the opposite direction of the copper coil along the vertical direction. The compliant flexure ensures that the coil remains aligned with the core and transmits forces to the outer passive structure. Attached to the right coil vertex is a Texas Instruments linear Hall effect sensor (DRV5057A3), which produces a voltage proportional to the magnetic field it senses. When the hexagon deforms along its actuator axis, the distance between the core and the Hall sensor changes, resulting in a voltage difference which we calibrate to the corresponding deformation of the hexagon.

We now connect sets of three hexagons together as shown in Fig. LABEL:COARSENING:trihexd, with the actuators pointing symmetrically along the three principle axes of the hexagons which guarantees isotropy at the large scale. Each hexagon contains a flexible circuit that connects its Hall sensor and the terminals of its copper coil to a custom made printed circuit board (pcb), that sits on top of the central vertex. Integrated in the pcb is a 32-bit micro controller (ATSAMD21E18A) that controls the three bidirectional motor drivers, which can manipulate the current through the voice coil. The central circuit board collects Hall voltages from each hexagon and converts these signals into the corresponding deformations $\delta L=L_{0}-L$ where $L_{0}=2l$ is the hexagon’s long diagonal at equilibrium in terms of the hexagon side length $l=2$cm. The pcbs are powered by a 30A current source. We now implement non-reciprocal interactions between the forces exerted by the actuators and the shear deformations of neighboring hexagons according to the following antisymmetric dynamical matrix:

=0, whose determinant gives the growth rates $\sigma$ of perturbations $k$ to a travelling wave at wavenumber $q$ (Fig. LABEL:fig:SecondaryStability). For finite friction $\Gamma$ we find a band of modes about the dominant wavenumber $q=\Gamma^{1/2}$ of width $\sim\Gamma^{1/2}$ which are stable. Waves outside this band are vulnerable to a Benjamin-Feir instability. As $\Gamma\rightarrow 0$ this stable band vanishes, and at zero friction we find that all modes are unstable towards phase slips.

V Observation of coarsening dynamics in the $\phi$ model

The above stability analyses suggest a two-step process. First, non-reciprocity will pick a travelling pattern of given wavenumber through a type-II linear instability. Next, nonlinear couplings will destabilize this mode. What is the eventual fate of this travelling pattern? In Fig. S5(a), we show simulation snapshots of the $\phi$-model in a 1D geometry under periodic boundaries. Above the threshold for linear instability, we initially see a dominant wavelength emerge. Over time, however, this mode undergoes a sequence of secondary instabilities, which cascade towards ever lower mode number, and ever larger amplitude. This observation is confirmed by the amplitude spectrum, whose peak shifts to lower modes and larger amplitudes over time [Fig. S5(b,c)]. In one dimension, the cascade of wavenumbers is mediated by phase slips, events where the dominant wavenumber suddenly decays and a lower mode takes over [Fig. S5(c)]. Tracking the inverse correlation length of such data for varying $\alpha$ in both 1D and 2D reveals that this coarsening process occurs via a power law, which we show in Fig. 2 of the Main Text.

VI Mean field approximation to the $\phi$ model

In the Main Text, we approximated the full mode equation corresponding to a nonlinear odd elastic solid by neglecting cross coupling between modes, and reducing from a second order dynamics to a first-order dynamics that captures the coarsening process but neglects fast oscillations. Here we provide numerical results that justify the neglect of cross couplings, effectively producing a mean field model. Next, we apply a multiple scales approach to reduce this second order mean field model to a first order model.

Reduction to the mean field model

We write Eq. (LABEL:eq:phi_modelSI) in mode space, $\phi=\sum_{n}a_{n}e^{iq_{n}x}$, and decompose the nonlinear term into an onsite nonlinearity, a bath interaction term and a phase dependent term:

[TABLE]

where we have parametrized the bath and phase dependent terms by $B,C\in\{0,1\}$. In order to assess the validity of the approximation in the main text that $C=0$, we integrate this equation using a spectral method for the cases that $B=C=0$, $B=1,C=0$, and $B=C=1$. In the case that $B=C=0$, the equation represents a collection of uncoupled nonlinear oscillators (Fig. S6a). Although the lowest mode eventually dominates, the absence of interaction terms means that there is a broad distribution of modes at all times (Fig. S6b).

When $B=1,C=0$, we recover the approximated model discussed in the main text which captures the features of the experimental data observed in Figs 1 and 4 (Fig. S6c). As one mode starts to dominate, the next one grows and eventually takes over until only the lowest mode is left. As demonstrated in the next section, the existence of a nonlinear growth rate implies that only a narrowly distributed band of modes is excited at any time. During coarsening, this band shifts at a rate proportional to $\propto t^{-1/2}$(Fig. S6d).

Finally, when $B=C=1$ we recover the full model where modes resonate during the mode evolution (Fig. S6e). Nonetheless, the dominant mode still coarsens over time and tends to the same $-1/2$ power-law as the mean-field model (Fig. S6f).

Reduction from second to first order

Having established that the wave coarsening can quantitatively be captured in the mean field setting, we analyze this simplified model. For ease of notation, in this section we use $q=q_{n}=2\pi n/L$ interchangeably where it is clear from context. Assuming that the on-diagonal bath term dominates couplings between modes in Eq. (S32) yields

[TABLE]

Eq. (S33) is second order in time. However, we now show that a multiple scales approach can reduce the dynamics to an effectively first order description. We make an exponential coordinate transformation, $a_{n}(t)=e^{\sigma_{n}(t)}$. The growth rate of amplitudes $\dot{|a_{n}|}$ and instantaneous frequencies $\omega_{n}=\dot{\theta}_{n}$ of $a_{n}=|a_{n}|e^{i\theta}$ are recovered via the relation $\sigma_{n}=\dot{|a_{n}|}/|a_{n}|+i\omega_{n}$. In terms of this transform, Eq. (S33) reads

[TABLE]

Here, $\dot{\sigma}_{n}$ is a time dependent function with instantaneous frequency $\omega_{n}:=\text{Im}(\dot{\sigma}_{n})$ and instantaneous exponential growth $g_{n}:=\text{Re}(\dot{\sigma}_{n})$. We note that there are two timescales in Eq. (S34). In the absence of the $\ddot{\sigma}_{n}$ term, Eq. (S34) is an algebraic relation for $\dot{\sigma}_{n}$. We hypothesize that this algebraic relation gives the dynamics of $\dot{\sigma}_{n}$ on a slow timescale. We therefore split $\dot{\sigma}_{n}$ as

[TABLE]

where we define $\dot{\langle\sigma\rangle}_{n}$ as a slow time variable which satisfies

[TABLE]

with a perturbative remainder $\dot{\delta}_{n}$ capturing the fast time dynamics. We look at the stability of small perturbations, such that we assume $\dot{\delta}_{k}$ is small. Substituting Eq. (S35) into Eq. (S34) and using Eq. (S36) yields

[TABLE]

We have that $\dot{\delta}^{2}$ is small, and by assumption of multiple timescales $\ddot{\langle\sigma\rangle}$ is also small. We obtain

[TABLE]

Equation (S38) is now solvable on fast timescales. In summary, we obtain

[TABLE]

Taylor expanding from time $t_{0}$ to obtain $\dot{\delta}_{n}(t_{0}+\Delta t)=\dot{\delta}_{n}(t_{0})e^{-2\dot{\langle\sigma\rangle}_{n}(t_{0})\Delta t}$. Since $g_{n}\approx\text{Re}(\dot{\langle\sigma\rangle}_{n})$ gives the instantaneous exponential growth and $\omega_{n}\approx\text{Im}(\dot{\langle\sigma\rangle}_{n})$, we have $\dot{\delta}_{n}(t_{0}+\Delta t)=\dot{\delta}_{n}(t_{0})e^{-2g_{n}(t_{0})\Delta t-2i\omega_{n}(t_{0})\Delta t}$. Therefore, for a growing mode, the perturbation decays exponentially fast in an oscillatory manner. The hypothesis encoded in Eq. (S35) is self-consistent: after a short time interval we find that $\dot{\delta}=0$, and in the coarsening regime $\dot{\sigma}_{n}\rightarrow\dot{\langle\sigma\rangle}_{n}$.

We now consider the limit of large nonlinearity $\sum_{m}|a_{m}|^{2}$. Our linear stability analysis Eq. (LABEL:eq:roots) shows that amplified modes are bounded by $\alpha$. We therefore Taylor expand the root in Eq. (S40) by taking $1+2\sum_{m}|a_{m}|^{2}$ to be large compared to $\alpha$ and $q_{n}$. We obtain

[TABLE]

We take the real part of $\langle\dot{\sigma}\rangle_{n}$ to obtain the instantaneous exponential growth rate and the imaginary part to obtain the instantaneous dispersion. We obtain for the frequency

[TABLE]

and for the amplitude $|a_{n}|$

[TABLE]

where Eq. (S43) corresponds to Eq. (3) of the Main Text for the case where $p=1$, and represents a first order dynamics for the nonlinear growth rate of mode amplitude $|a_{n}|$. The structure of Eq. (S43) is identical to the linear dispersion Eq. (LABEL:eq:Dispersion), except for a dynamic rescaling of $\alpha$ by a factor $\sqrt{1+2\Sigma_{m}|a_{m}|^{2}}$. Intuitively, as this bath term grows, the effective maximal growth rate of the linear dispersion is shifted towards $q\rightarrow 0$, and ever-lower modes are preferentially amplified at ever lower rates (Fig. S7). To solve Eq. (S43), we multiply by $|a_{n}|$ and obtain

[TABLE]

Appealingly Eq. (S44) suggests an argument for the coarsening exponents. Suppose coarsening occurs via the motion of a tightly peaked packet of wavenumbers, such that the mean wavenumber is much larger than the packet variance. We then convert Eq. (S44) into an equation involving the total amplitude and central wavenumber $Q$ only. We proceed by obtaining the first two moments of the normalized distribution $2|a_{n}|^{2}/A$ as

[TABLE]

Given the concentrated mode distribution, we have that its variance $s$ is much smaller than the square of the central wavenumber $Q$, such that $\langle q^{2}\rangle=Q^{2}+s^{2}\approx Q^{2}$. Using this, we sum both sides of Eq. (S44) over all wavenumbers $q$ and multiply by $2$ to obtain

[TABLE]

Finally we make the ansatz $A(t)\sim t^{\chi}$, $Q(t)\sim t^{\beta}$. We find

[TABLE]

which only has a solution if the powers are the same for each term. From the equation $\chi-1=2\beta+\chi$ we obtain $\beta=-1/2$. From the equation $\chi-1=\beta+\frac{1}{2}\chi$ we obtain $\chi=1$. We therefore predict a characteristic scale for the peak wavenumber, and amplitude, of

[TABLE]

Numerical verification of the second to first order reduction

We make several assumptions in deriving Eq. (S43) from Eq. (S33), which we now numerically verify. In Fig. S8 we compare the simulated dynamics of Eq. (S33) with $n=20$ modes to the simplified dynamics and coarsening exponents predicted in Eqs. (S42), (S43) and (S50). In Fig. S8a we plot instantaneous frequency $\omega_{n}(t):=\dot{\theta}_{n}(t)$ compared to the frequency slaving prediction $\omega_{n}\sim\sqrt{\Sigma|a_{m}|^{2}}$ given by Eq. (S42). We find excellent agreement after the initial linear instability. Eq. (S50) further predicts power law decay of $\langle q\rangle\sim t^{-1/2}$ and mean amplitude growth of $\sqrt{\Sigma|a_{m}|^{2}}\sim t^{1/2}$. In Fig. S8b we plot these quantities for the second order dynamics Eq. (S33) and the first order dynamics Eq. (S43). After the initial linear instability we find good agreement with our predicted exponents.

VII The nonreciprocal beam and its coarsening exponents

In this section we apply the analysis performed above on the $\phi$ model to a non-reciprocal beam theory of the ring pictured in Fig. 3, with the continuum model given by Eq. (2). We begin with the linearized dynamics of a non-reciprocal beam derived in Ref. [16]:

[TABLE]

where $h$ describes the transverse deflection of the vertices, which have mass $m$, are spaced apart by distance $a$ and are connected with torsional bonds with stiffness $\kappa$, viscous damping $\gamma$ and non-reciprocal interactions $\kappa^{a}$. The sign of of $\kappa^{a}$ in Eq. (S51) is chosen to amplify right-moving waves.

Non-dimensionalizing this equation yields

[TABLE]

with $\alpha=(\kappa^{a}/\kappa)(\kappa m)^{1/4}\gamma^{-1/2}a^{1/2}$ the dimensionless non-reciprocity and $L$ the ring size. To mimic the saturation of the torque motors embedded in the metamaterial of Fig. 3, we replace the linear active term $\propto\partial_{\lambda}^{5}h$ by a saturating nonlinear function that breaks spatial inversion symmetry whilst retaining the Laplacian form of Eq. (S52). We opt for the following substitution:

[TABLE]

and find Eq. (2):

[TABLE]

Mean field and coarsening of the beam model

We now perform a series of simplifications to the beam model Eq. (S54) which are identical to those performed on the $\phi$ model Eq. (LABEL:eq:phi_modelSI) in §VI. Fourier transforming as $h=\sum_{n}a_{n}e^{i{q_{n}}x}$ we obtain

[TABLE]

where the Fourier transform \widehat{\left(\frac{\partial_{x}^{3}h}{\sqrt{1+(\partial_{x}^{3}h)^{2}}}\right)}\Bigl{.}_{\!n} equals $\int^{L/{2}}_{-L/2}\frac{\partial_{x}^{3}h}{\sqrt{1+(\partial_{x}^{3}h)^{2}}}e^{i{q_{n}}x}dx$. We want to obtain only the phase-independent terms, which we can do before the integration over $x$. For a real function $h$ we have $(\partial_{x}^{3}h)^{2}=(\sum_{l}(i{q_{l}})^{3}a_{l}e^{i{q_{l}}x}+\textrm{c.c.})(\sum_{m}(i{q_{m}})^{3}a_{m}e^{i{q_{m}}x}+\textrm{c.c.})=2\sum_{m}q_{m}^{6}a_{m}a_{-m}+\text{phase dependent terms}$. We insert this sum over modes into the Fourier transform of the nonlinear term and obtain

[TABLE]

We now choose to write $a_{n}=e^{\sigma_{n}}$, as in §VI. We insert this transformation into Eq. (S56) and multiply with $e^{-\sigma_{n}}$ to obtain

[TABLE]

We now repeat the logic described in §VI for the $\phi$ model by writing $\dot{\sigma}_{n}=\dot{\langle\sigma\rangle}_{n}+\dot{\delta}_{n}$. We obtain an expression that is analogous to Eq. (S40),

[TABLE]

We again are interested in the coarsening, highly nonlinear regime. Taking $\sum_{m}q_{m}^{6}|a_{m}|^{2}$ to be much larger than $\alpha$ we Taylor expand to obtain

[TABLE]

Here the imaginary component gives the instantaneous frequency $\omega_{n}$, such that

[TABLE]

The real part gives the instantaneous growth rate such that

[TABLE]

which corresponds to Eq. (3) of the main text for the case where $p=2$. We proceed by obtaining the exponents for the total amplitude $A$ and central wavenumber $Q$. We define the total amplitude as $A:=\sum_{m}2|a_{m}|^{2}q_{m}^{6}$, while we define the $j$-th moment as $\langle q^{j}\rangle:=\sum_{m}q_{m}^{j}(2|a_{m}|^{2}q_{m}^{6})/A$. We again assume a concentrated distribution of amplitude, such that the standard deviation $s\ll Q$, such that $\langle q^{j}\rangle\approx Q^{j}$. We again perform the same analysis as for the $\phi$ model, except that we now also multiply both sides by $q_{l}^{6}$, and we obtain

[TABLE]

We again insert the ansatz $Q\propto t^{\beta}$ and $A\propto t^{\chi}$. We likewise see that by equating the powers in $t$, Eq. (S62) can only hold when $\beta=-1/4$ and $\chi=1$. We therefore obtain

[TABLE]

The prediction $\langle q\rangle\sim\xi^{-1}\sim t^{-1/4}$ from Eq. (S64) is shown in Fig. 3e of the Main Text.

VIII Unidirectionally amplified Rayleigh waves in odd elastic media

In this section we derive the dispersion relation for Rayleigh waves in an odd elastic medium, and show that they are unidirectionally amplified. Rayleigh waves have been previously been shown to edge localize under open boundaries, with some numerical evidence of their amplification [16, 59].

Our starting point is linear odd elastodynamics Eq. (LABEL:eq:LinearDisplacement). We consider the situation of zero internal torque ($A=0$), zero viscous dissipation $(\eta=0)$ and zero substrate friction ($\Gamma=0)$. Our analysis nevertheless provides insight into a viscously dissipative medium, but only as $q\rightarrow 0$. Finite wavenumbers of $O(q^{2})$ will undergo viscous damping not considered here, which will stabilize the dispersion at high $q$. We consider a bulk medium occupying the $y>0$ half-plane, with a no stress boundary at $y=0$:

[TABLE]

We take a wave ansatz

[TABLE]

In general, $q_{x},q_{y},\omega\in\mathbb{C}$. Here we focus on the case $q_{x}\in\mathbb{R}$, but $\omega,q_{y}\in\mathbb{C}$, corresponding to a wave propagating along $x$ with complex growth rate $\omega$ and penetration depth $q_{y}=i\kappa$. We nondimensionalize by $q_{x}$, defining a dimensionless penetration depth $r=q_{y}/q_{x}$ and Rayleigh wave speed $v=\omega/q_{x}$. Our bulk dispersion relation then reads

[TABLE]

We now require the determinant of Eq. (S67) to vanish, giving us a secular equation

[TABLE]

with roots

[TABLE]

where $a,b=\pm 1$ index the four possible roots. Each root comes with a bulk wave eigenvector

[TABLE]

We now combine these eigenvectors to construct Rayleigh waves that satisfy the no-stress boundary conditions. In Fourier space, these boundary conditions read

[TABLE]

Given a displacement vector $\tilde{u}$, this matching condition can be written as a linear map

[TABLE]

We now consider this displacement $\tilde{u}$ to be a linear mix of two candidate Rayleigh wave eigenvectors Eq. (S70) which come from solving the bulk dispersion:

[TABLE]

with the two modes indexed by the two tuples $(a_{1},b_{1})$ and $(a_{2},b_{2})$. Our stress-matching condition Eq. (S73) can then be written in matrix form

[TABLE]

where each of the $S(n)$ are now known functions. The $\binom{4}{2}=6$ different tuple pairs $(a_{1},b_{1}),(a_{2},b_{2})$ encode all the possible stress matches coming from mixing eigenvectors. We require the determinant of $S$ to vanish,

[TABLE]

The solvability condition Eq. (S76) gives an algebraic equation in $v$. We simplify this expression by considering the incompressible limit $B\rightarrow\infty$. Inspection of the resulting polynomial reveals the Hermitian Rayleigh wave to be captured by $(a,b)=(+1,+1),(1,-1)$. With this choice Eq. (S76) reads

[TABLE]

where $s=v/\sqrt{\frac{\mu}{\rho}}$ gives the complex Rayleigh wave speed in terms of the shear wave speed. To lowest order in $K^{o}/\mu$ we find that the solutions to this polynomial which have positive penetration depths are

[TABLE]

where $s_{r},s_{i}$ are numerical approximations to exact radical expressions. Here $s_{r}=0.955\ldots$ is the standard speed of an incompressible Rayleigh wave. By contrast, $s_{i}=0.269\ldots$ gives the imaginary part of $s$, and represents wave amplification. Note that only the phase velocity $s_{+}$ with the same sign as $K^{o}$ is amplified, with $s_{-}$ attenuated (Fig. S9).

Using Eq. (S69), we find the decay lengths of these roots to order $K^{o}/\mu$ are

[TABLE]

in terms of the dimensionless penetration depth $\kappa=-ir$. The real parts of $\kappa$, $\mathrm{Re}(\kappa)=1,\ \kappa_{r}=0.296\ldots$ give the penetration depth of an incompressible Rayleigh wave in a passive medium. In active media we find that these penetration depths also acquire a spatially oscillatory correction $\kappa_{i}=0.871\ldots$. We conclude that Rayleigh wave states in elastic media continue to exist in the odd elastic context, with the first order correction of non-reciprocity being amplification of these waves in time. Moreover this amplification is thresholdless: even for values of $K^{o}$ at which bulk modes are stable, surface modes are still amplified. We observe these surface-bound states in the coarsening dynamics shown in Fig. 4 of the Main Text.