Martin Brandenbourger

Colin Scheibner

Jonas Veenstra

Vincenzo Vitelli
email: vitelli@uchicago.edu

Corentin Coulais
email: coulais@uva.nl

Limit

Limit cycles turn active matter into robots

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

low it to roll autonomously uphill by virtue of its own deformation, as demonstrated by our prototypes. Similarly, familiar scattering phenomena, like a ball bouncing off a wall, turn into basic robotic manipulations when either the ball or the wall is odd. Using continuum mechanics, we reveal collective robotic mechanisms that steer the outcome of collisions or influence the absorption of impacts in experiments. Beyond robotics [1-3, 5, 6, 21-27],

work-generating limit cycles can also control the non-linear dynamics of active soft materials [START_REF] Palacci | Living crystals of lightactivated colloidal surfers[END_REF][START_REF] Aubret | Metamachines of pluripotent colloids[END_REF][START_REF] Sanchez | Spontaneous motion in hierarchically assembled active matter[END_REF][START_REF] Bricard | Emergence of macroscopic directed motion in populations of motile colloids[END_REF][START_REF] Bililign | Motile dislocations knead odd crystals into whorls[END_REF][START_REF] Keber | Topology and dynamics of active nematic vesicles[END_REF][START_REF] Giomi | Spontaneous division and motility in active nematic droplets[END_REF][START_REF] Baconnier | Selective and collective actuation in active solids[END_REF], biological systems [START_REF] Tan | Development drives dynamics of living chiral crystals[END_REF][START_REF] Peyret | Sustained oscillations of epithelial cell sheets[END_REF][START_REF] Gilpin | The multiscale physics of cilia and flagella[END_REF][START_REF] Marder | Central pattern generators and the control of rhythmic movements[END_REF][START_REF] Lavi | Deterministic patterns in cell motility[END_REF] and driven nanomechanical devices [START_REF] Miskin | Electronically integrated, mass-manufactured, microscopic robots[END_REF][START_REF] Miskin | Graphene-based bimorphs for micron-sized, autonomous origami machines[END_REF][START_REF] Mcdonald | Exponentially-enhanced quantum sensing with non-hermitian lattice dynamics[END_REF][START_REF] Mathew | Synthetic gauge fields for phonon transport in a nanooptomechanical system[END_REF].

We start by asking a simple question: what is a robot? Here, we take the following perspective: an autonomous robot-one that follows its own internal rules without receiving an external set of instructions-can be viewed as an autonomous dynamical system-one in which the time variable t does not appear in the governing equation [START_REF] Aguilar | A review on locomotion robophysics: the study of movement at the intersection of robotics, soft matter and dynamical systems[END_REF][START_REF] Strogatz | Nonlinear Dynamics and Chaos[END_REF]. In this picture, cyclic robotic functionalities, such as climbing or locomotion, are enabled by robust limit cycles in the underlying dynamical system [START_REF] Asano | Limit cycle gaits[END_REF]. Limit cycles arising in living systems [START_REF] Gilpin | The multiscale physics of cilia and flagella[END_REF][START_REF] Marder | Central pattern generators and the control of rhythmic movements[END_REF][START_REF] Lavi | Deterministic patterns in cell motility[END_REF][START_REF] Brückner | Stochastic nonlinear dynamics of confined cell migration in two-state systems[END_REF], active matter [START_REF] Marchetti | Hydrodynamics of soft active matter[END_REF], and robotic walkers [START_REF] Asano | Limit cycle gaits[END_REF][START_REF] Conradt | Distributed central pattern generator control for a serpentine robot[END_REF] can be driven by external stimuli, constant energy fluxes, or inertia. Here, we investigate limit cycles that emerge due to non-conservative, or odd, mechanical forces that naturally provide the work needed to perform the desired robotic functionalities. Nonconservative (odd) matter is made of repeated components that obey simple local interactions, very much like ordinary matter composed of atoms and molecules, but with a crucial difference. In odd matter, the interaction forces between basic constituents are not constrained to be the gradient of a potential-hence by taking the system along a cycle work is generated [START_REF]is distinct from other possible usages, such as odd as in parity violating, odd stress as in continuous media with local torques, and odd response coefficients that violate Onsager reciprocity[END_REF].

As an example, consider the familiar concept of bond-bending stiffness. Ordinary angular or bond-bending interactions, such as those in pendula or molecules, are approximated by a potential energy E = κ 2 i δθ i 2 where the subscript i labels each bond, δθ i is the corresponding deviation from the equilibrium angle and κ is the bending stiffness. The corresponding torque, obtained from gradients of the potential, is then given by the (generalized) force-displacement relation τ i = -κδθ i . For a simple pendulum κ is set by gravity. We now revisit this familiar concept by presenting a minimal non-conservative generalization of bond-bending stiffness illustrated in Fig. 1a-c and Supplemental Video 1. Consider a system of rigid bars described by two bond angles θ 1 and θ 2 (schematics in Fig. 1d-f) whose force-displacement relation reads:

  τ 1 τ 2   =   -κ -κ a κ a -κ     δθ 1 δθ 2   (1)
or equivalently τ i = M ij δθ j . The coefficient κ appears in the diagonal entries of M ij and represents the familiar bond-bending stiffness. In addition to κ the matrix M ij also has an off diagonal component κ a , which represent a coupling between the two angles. In a passive system, e.g., two pendula coupled by a spring, the off diagonal couplings would be equal, so M ij would be symmetric. [42]. Here, we remove this restriction and allow an antisymmetric contribution in the matrix, i.e. M ij = M ji . We call this contribution represented by the coefficient κ a odd, since it is antisymmetric under exchange of the indexes i and j in M ij . Figure 1 illustrates the physical meaning of this oddness: for κ a > 0 when θ 1 is contracted, θ 2 contracts too (Fig. 1b,e); when θ 2 is contracted, θ 1 expands instead (Fig. 1c,f). The crucial point here is that this asymmetric relationship cannot exist in passive matter because the force law in Eq. ( 1) is nonconservative. When the vertices are deformed the work done by the vertices W = κ a (δθ 2 δ θ1 -δθ 1 δ θ2 ) dt depends on the path taken. Hence, if the path is cyclic, the work need not be zero; it can either be positive or negative depending on how one cycles around the loop. This implies that the force law in Eq. ( 1) requires a source of energy to be sustained [43].

When odd force laws such as Eq. ( 1) are combined with inertia, intuition alone suggests their dynamical fate: any initial perturbation will grow, powered by the non-conservative forces, until non-ideal facets of the system (e.g. friction, self-collisions, or breaking) stabilize the dynamics.

In some cases, the dynamical system defined by the generalized coordinates and their conjugate momenta will proceed to a limit cycle in which the energy injected by the non-conservative force balances with dissipation. As an illustration, consider the equations of motion for the non-conservative bonds in Fig. 1:

Iδ θi = -κδθ i + f i (δθ 1 , δθ 2 ) -Γδ θi ( 2
)
where I is the moment of inertia, Γ is a dissipation coefficient, and f i is a nonlinear and nonconservative function. A generic form of nonlinearity, relevant also for our prototypes, is a saturation of the motors at a maximum torque τ max (Fig. 1h)-other geometric or dissipative sources of nonlinearities may play a similar role in non-conservative active and biological systems [START_REF] Tan | Development drives dynamics of living chiral crystals[END_REF][START_REF] Fruchart | Non-reciprocal phase transitions[END_REF]. Using the complex variable

z = δθ 1 + iδθ 2 and f = f 1 + if 2 we can model the nonlinear function as f (z) = -iκ a z if |z| < τ max /κ a , and f (z) = -iτ max z/|z| if |z| > τ max /κ a ,
see Fig. 1h. Notice that f is a nonlinear extension of Eq. [START_REF] Aguilar | A review on locomotion robophysics: the study of movement at the intersection of robotics, soft matter and dynamical systems[END_REF]. Using the ansatz z = z lim e -iωt , we obtain the solution [START_REF]Solving for the amplitude |z lim | yields: |z lim | = iτ max /(-Iω 2 + κ -iΓω). The constraint that |z lim | must be positive and real restricts the allowed frequencies to ω[END_REF]

δθ 1 + iδθ 2 = τ max Γ I κ cos κ I t + i sin κ I t (3)
up to an arbitrary phase. Crucially, Eq. ( 3) is only a valid solution for sufficiently large drive.

The dimensionless parameter ξ = κ a /Γ I/κ captures the ratio of drive (κ a ) and inertia (I)

to dissipation (Γ) and restoring forces (κ). As shown in Fig. 1i, when |ξ| < 1, dissipation and restoring forces win and the process terminates at its rest configuration z = 0. As shown in Fig. 1j, when |ξ| > 1, a Hopf bifurcation occurs and there exists a cycle of finite amplitude z lim = τ max /Γ I/κ. Mathematically, the trajectory in Eq. ( 3) is a limit cycle of the dynamical system in Eq. [START_REF] Aguilar | Collective clog control: Optimizing traffic flow in confined biological and robophysical excavation[END_REF]. We refer to this type of limit cycle as a nonlinear work cycle because it arises mechanically as the balance between dissipation and a non-conservative, configurationdependent force that performs the work necessary to sustain it [46].

We first illustrate that a robust tendency to cycle can be harnessed to perform robotic functionalities when the system is in contact with an environment (see Supplemental Video 2). In Fig. 2a, six such vertices are connected on a closed hexagon and are levitated on an air-table

(see Methods). In addition to solid-body translations and rotations, the hexagon has three independent deformation modes: two shear modes S 1 and S 2 and one breathing mode B (Fig. 2b).

Due to the odd stiffness κ a , when the hexagon is initially poked, it undergoes a spontaneous cycle of large shape change oscillating between the two shear modes (Fig. 2c). When the table is tilted, the hexagon comes in contact with an inclined railing. The repetitive contact provides a means of locomotion: the hexagon acts as an odd wheel. It overcomes gravity and locomotes uphill on an inclined ramp (Fig. 2d) or on a rough terrain (Extended Data Fig. 1). Unlike an ordinary wheel that is driven by external torques, the odd wheel is perched on the brink of an instability that allows it to roll autonomously uphill by virtue of its own deformation. The dominant mechanism is an interplay between inertial rolling and internal shape change: the bottom surface of the hexagon tends to thrust left against the railing, propelling the center of mass to the right (see Fig. 2e and Supplemental Video 2). This mechanism emerges from the minimal electronic feedback, rather than a directly engineered [START_REF] Chiu | deformable wheel"-a self-recovering modular rolling track[END_REF][START_REF] Shen | Multimode locomotion via superbot robots[END_REF][START_REF] Shen | Multimode locomotion via superbot reconfigurable robots[END_REF][START_REF] Sastra | Dynamic rolling for a modular loop robot[END_REF][START_REF] Chiu | Multifunctional superbot with rolling track configuration[END_REF][START_REF] Shen | Hormone-inspired adaptive communication and distributed control for conro self-reconfigurable robots[END_REF] or on-the-fly [START_REF] Li | Particle robotics based on statistical mechanics of loosely coupled components[END_REF][START_REF] Oliveri | Continuous learning of emergent behavior in robotic matter[END_REF] locomotion algorithm. In contrast to limit-cycle walking based purely on inertia [START_REF] Asano | Limit cycle gaits[END_REF], the limit cycles on display in Fig. 2 results from a balance of energy injection with dissipation, and hence can power uphill motion.

Beyond cyclic functionalities, such as locomotion, we show that odd matter can perform basic robotic manipulations such as steering motion and forces. As an example, we turn to a most basic probe of matter equally familiar to scientists or kids alike: smashing things into each other. We consider what happens when a ball collides with a wall when either the ball or the wall are odd. In Fig. 3ab a twelve-sided projectile with a positive odd stiffness κ a > 0 is launched vertically towards an inelastic passive surface. Unlike a carefully swatted tennis ball, the incoming odd projectile has no initial rotation. Yet it deforms asymmetrically upon impact and thereby steers its outgoing trajectory to an angle determined by the sign of the odd stiffness κ a (Extended Data Fig. 1a and Supplemental Video 3). To make sense of this robotic functionality, we examine the transient dynamics inside the odd projectile. In Fig. 3bc, we reconstruct the full internal degrees of freedom of the projectile and plot the angular deviation of each vertex as a function of time in the associated heat map. The qualitative similarity between the active and passive bending stiffness in Eq. ( 1) suggests an approach to understanding the collision process: collective modes. In Fig. 3d, we project the vertex dynamics onto the lowest two Fourier modes of the ring. Physically, these Fourier modes represent a shearing deformation analogous to the S 1 and S 2 deformation of the hexagon. The out-spiraling motion in Fig. 3d captures the lowest mode contributions to the unidirectional wave propagating counterclockwise in Fig. 3c. This out-spiraling behavior is powered by activity and is eventually stabilized by self intersection, dissipation, and motor nonlinearities.

We now show that treating the odd projectile as a non-conservative continuous material provides basic insights into the outcome of the impact. We model the chain as a continuum fiber with periodic boundary conditions (Fig. 3e-top inset). By linearizing the equations of motion and taking the continuum limit of a 1D chain of active vertices (see S.I.), we arrive at the following distinctive wave equation:

(∂ 2 τ + ∂ 4 λ + 2κ a ∂ 5 λ + Γ∂ 4 λ ∂ τ )h(λ, τ ) = 0, (4)
where h(λ, τ ) is the amplitude of the fiber's deflection as a function of time τ and space λ, while Γ is the dissipation, κa = κ a /N κ, and N is the number of unit cells, see Methods. When κ a = 0, Eq. ( 4) is the standard wave equation for a thin beam with bending stiffness. Crucially, κ a = 0, a new term is introduced that violates λ → -λ. The broken parity and Hermiticity imply that the waves will be amplified in one direction relative to the other, an effect known as the non-Hermitian skin effect [START_REF] Shankar | Topological active matter[END_REF][START_REF] Brandenbourger | Non-reciprocal robotic metamaterials[END_REF][START_REF] Martinez Alvarez | Non-hermitian robust edge states in one dimension: Anomalous localization and eigenspace condensation at exceptional points[END_REF][START_REF] Yao | Edge states and topological invariants of non-hermitian systems[END_REF][START_REF] Ghatak | Observation of non-hermitian topology and its bulk-edge correspondence in an active mechanical metamaterial[END_REF][START_REF] Helbig | Generalized bulk-boundary correspondence in nonhermitian topolectrical circuits[END_REF][START_REF] Xiao | Non-hermitian bulk-boundary correspondence in quantum dynamics[END_REF][START_REF] Chen | Realization of active metamaterials with odd micropolar elasticity[END_REF][START_REF] Coulais | Topology and broken hermiticity[END_REF][START_REF] Bergholtz | Exceptional topology of non-hermitian systems[END_REF][START_REF] Scheibner | Non-hermitian band topology and skin modes in active elastic media[END_REF][START_REF] Zhou | Non-hermitian topological metamaterials with odd elasticity[END_REF].

The unidirectional propagation is associated with shape-changes that control the bounce.

Initially, when the projectile compresses under its own inertia, the deformation induced in the impact corresponds to the ring's second harmonic, of wavelength L/2. Starting from this initial condition, the analytical solution for the curvature ∂ 2 λ h over time (Fig. 3e-bottom inset) captures the angular deflection at each vertex appearing in the experimental color-map of Fig. 3c. This unidirectional wave is responsible for the asymmetry and the enhancement of the rebound.

While in contact with the surface, the bumps of negative curvature (blue in Fig. 3f inset) on the sides propagate in one direction: the left bump travels down, while the right bump travels up.

The downward push on the left propels the center of mass to the right and gives it an additional vertical thrust. As implied by symmetry considerations, in Fig. 3f, the vertical momentum transfer p y,out /p y,in is even in κ a while the transverse momentum transfer p x,out /p y,in is odd.

Using a minimal model of a solid block only capable of shearing, we are able to derive analytical expressions (solid lines in Fig. 3f, see Methods) that capture the general experimental trends (red and blue markers for κ a > 0 and κ a < 0 respectively).

Finally, we examine an additional robotic functionality: an odd wall can steer post-impact vibrations and the outgoing trajectory of the bullet. We consider the scenario of a passive projectile striking against a two-dimensional odd wall with positive odd stiffness κ a > 0 (Fig. 4a) and negative odd stiffness κ a < 0 (Extended data Fig. 2). The odd wall consists of 17 odd hexagons (120 powered vertices in total) patterned into a honeycomb lattice (see S.I.). Upon impact, we observe a rotation and deflection of the passive projectile and a high degree of asymmetry of the vertical displacement field u y (x, y) in the odd solid, reconstructed from experiments in Fig. 4b.

Fig. 4c shows u y (x, t) averaged in the y direction, revealing that the wall undergoes stronger vibrations on its left edge during and after the impact (See Supplemental Video 4). In Fig. 4d, we see that the spatially averaged components (S 1 and S 2 ) of the shear strain in the wall exhibits a cyclic coupling reminiscent of the locomotion in Fig. 2f and impact in Fig. 3d .

Since the impact absorption relies on the coordinated motion of many interacting components, we ask whether the emergent features can be captured in a way agnostic to microscopic details, i.e. via continuum elasticity. Passive, isotropic 2D elasticity is summarized by just two parameters: the Young's Modulus E and the Poisson's ratio ν. For media with nonconservative interactions, an additional parameter exists: the odd ratio ν o [START_REF] Scheibner | Odd elasticity[END_REF]. The odd ratio captures the nonconservative relationship between shear stress and shear strain, see Methods. In the present realization, the odd ratio, a macroscopic quantity, relates to the microscopic odd stiffness as

ν o ≈ 0.58(κ a /κ) (see S.I.
). We can utilize this non-conservative continuum theory, known as odd elasticity [START_REF] Scheibner | Odd elasticity[END_REF], to rationalize basic features of the internal dynamics inside the wall under impact.

We quantify the steering of the post-impact vibrations by measuring the center of vibration defined in Eq. (M64) as a function of ν o in Fig. 4e. To rationalize the dependence of on ν o , we consider the ideal case of a half-infinite continuum medium with ν o = 0. In the Methods, we show that such a medium can host unidirectionally amplified Rayleigh waves (Fig. 4e-top inset). Crucially when ν o = 0, these waves acquires a nonzero exponential amplification along the boundary, quantified by an inverse growth length q x , such that u y (x) ∝ e qxx (Fig. 4e-bottom inset). This unidirectional amplification, reminiscent of the non-Hermitian skin effect [START_REF] Shankar | Topological active matter[END_REF][START_REF] Brandenbourger | Non-reciprocal robotic metamaterials[END_REF][START_REF] Martinez Alvarez | Non-hermitian robust edge states in one dimension: Anomalous localization and eigenspace condensation at exceptional points[END_REF][START_REF] Yao | Edge states and topological invariants of non-hermitian systems[END_REF][START_REF] Ghatak | Observation of non-hermitian topology and its bulk-edge correspondence in an active mechanical metamaterial[END_REF][START_REF] Helbig | Generalized bulk-boundary correspondence in nonhermitian topolectrical circuits[END_REF][START_REF] Xiao | Non-hermitian bulk-boundary correspondence in quantum dynamics[END_REF][START_REF] Chen | Realization of active metamaterials with odd micropolar elasticity[END_REF][START_REF] Coulais | Topology and broken hermiticity[END_REF][START_REF] Bergholtz | Exceptional topology of non-hermitian systems[END_REF][START_REF] Scheibner | Non-hermitian band topology and skin modes in active elastic media[END_REF][START_REF] Zhou | Non-hermitian topological metamaterials with odd elasticity[END_REF], is what determines the steering of post impact vibrations, i.e., the sign of in Fig. 4e main panel.

To intuit how the odd wall manipulates the outward trajectory of the passive bullet, we examine the asymmetry in the displacement field of the wall at the deepest point of impact. We consider an idealized situation in which a point force is exerted on the boundary of the semiinfinite odd medium (Fig. 4f). In this case, the continuum solution for the vertical component of the static displacement field u y takes the simple form (see Methods):

u y = 2F πE [2 log r + 1 -cos 2φ + ν o (sin 2φ + 2)] (5)
where r and φ are polar coordinates about the point of impact, F is the applied force, which is taken to be the mass of the projectile times its deceleration. The odd ratio, ν o in Eq. ( 5)

introduces a term that violates the φ → -φ symmetry, a continuum manifestation of the experimentally observed asymmetry. In Fig. 4g, we plot the outgoing angle α of the projectile as a function of ν o . The sign of α can be inferred qualitatively via the S 1 -S 2 shear coupling at the point of contact, as sketched in Fig. 4g (inset) for a single motorized hexagon and highlighted in Fig. 4a (inset) for the two powered hexagons touching the projectile. Beyond the first layer, the pattern continues smoothly further away from the point of impact Fig. 4b. This suggests that the discrete behavior should be evident in the continuum solution. Indeed, the sign of the slope (∂ x u y ∝ ν 0 , see Fig. 4f dashed line) at the point of contact in the continuum agrees with the sign of the tilt imparted to the projectile inferred from the discrete mechanism and observed in experiment.

We have shown how macroscopic robotic functionalities can emerge from simple nonconservative interactions. Identifying and engineering these cycles offers an approach to tailoring functionalitiies across platforms: the key ingredients to support these work-generating limit cycles are present in biological systems such as spinning embryos or bacteria [START_REF] Tan | Development drives dynamics of living chiral crystals[END_REF][START_REF] Petroff | Fast-moving bacteria self-organize into active twodimensional crystals of rotating cells[END_REF], biomembranes [START_REF] Needleman | Active matter at the interface between materials science and cell biology[END_REF][START_REF] Jülicher | Hydrodynamic theory of active matter[END_REF][START_REF] Banerjee | Active viscoelasticity of odd materials[END_REF][START_REF] Kole | Layered chiral active matter: Beyond odd elasticity[END_REF], or driven colloidal systems [START_REF] Bililign | Motile dislocations knead odd crystals into whorls[END_REF][START_REF] Yan | Rotating crystals of magnetic Janus colloids[END_REF][START_REF] Grzybowski | Dynamic self-assembly of magnetized, millimetre-sized objects rotating at a liquid-air interface[END_REF][START_REF] Aragones | Elasticity-induced force reversal between active spinning particles in dense passive media[END_REF], as well as micro-and nanomechanical devices [START_REF] Miskin | Electronically integrated, mass-manufactured, microscopic robots[END_REF][START_REF] Miskin | Graphene-based bimorphs for micron-sized, autonomous origami machines[END_REF][START_REF] Mathew | Synthetic gauge fields for phonon transport in a nanooptomechanical system[END_REF][START_REF] Liu | Micrometer-sized electrically programmable shape-memory actuators for low-power microrobotics[END_REF]. In analogy with how a sequence of amino acids folds into a functional protein, we envision an organic design process in which minimal active units evolve towards a functional dynamic state represented by a limit cycle [START_REF] Thubagere | A cargo-sorting dna robot[END_REF].

It is an open question how such principles generalize on scales in which thermal or environmental fluctuations are relevant

II. MATERIALS AND METHODS

A. Construction of the robotic materials

The linkages shown in Fig. 1a-c are composed of motorized vertices connected by plastic arms. Each vertex consists of a DC coreless motor (Motraxx CL1628) embedded in a cylindrical heatsink, an angular encoder (CUI AMT113S), and a microcontroller (ESP32) connected to a custom electronic board. The electronic board enables power conversion, interfacing between the sensor and motor, and communication between vertices. Each vertex has a diameter 50 mm, height 90 mm, and mass 0.2 kg. The power necessary to drive the motor is provided either by two Vapcell 16340 T8 batteries (Fig. 2d) or by an external 48 V DC power source (Figs. 3,4). Rigid 3D-Printed arms connect each motor's drive shaft to the heat sink of the adjacent unit. The angle formed between the two arms at vertex i is denoted by θ i . The on-board sensor measures θ i at a sample rate of 100 Hz and communicates the measurement to nearest neighbors.

In response to the incoming signal, vertex i exerts an active torsional force τ a i

τ a i = κ a (δθ i-1 -δθ i+1 ) . ( M1
)
where δθ i = θ iθ 0 and θ 0 is the rest angle. The active stiffness κ a was programmed by the In Fig. 2a-c, a device consisting of 6 vertices (shown in Fig. 2a) was subjected to a manually applied initial compression. Regardless of the exact initial perturbation, a stable cyclic motion emerged. We parameterize the internal state via the three degrees of freedom

S 1 = 1 √ 2 (δθ 1 -δθ 3 ) (M2) S 2 = 1 √ 5 (δθ 1 + δθ 3 -2δθ 5 ) (M3) B = 1 √ 3 (δθ 1 + δθ 3 + δθ 5 ) (M4)
The projection of experimental data onto the S 1 , S 2 axis is shown in Fig. 2c.

In Fig. 2d and Extended Data Fig. 1, the air table was tilted by 6 • with respect to the horizontal plane. This creates an effective gravitational constant of g eff ≈ 1 m s -2 that points towards the aluminum profile and tends to push the robotic hexagon against the profile. In turn, the aluminum profile (initially parallel to the floor) is itself rotated 5 • about the normal of the table. The total increase in altitude of the robotic hexagon when it climbs across the profile is thereby 1.5 m × sin 6 • × sin 5 • = 1.4 cm, see Supplemental Fig. S4. In Fig. 2d, strips of vinylpolysiloxane (Elite Double 32) are placed along the aluminum profile to provide additional friction between the robotic material and the boundary.

The experiments shown in Fig. 3 and Fig. 4 are performed on a horizontal air table. In Fig. 3, a projectile consisting of 12 active vertices is launched towards an aluminum profile. Two initial velocities, 0.18 m/s and 0.26 m/s, are used. In Fig. 4, a projectile of mass 8.2 kg was propelled at 0.5 m/s towards an active hexagonal lattice of vertices. The active feedback rule in Eq. (M1)

was implemented on vertices connected by the double line in Fig. 4b (see also Supplemental Fig. S3). The bottom row of the wall is anchored to the aluminum profile to prevent slipping.

For all experiments, the incoming projectile velocity was set by an aluminum profile connected to a conveyor belt and a stepper motor (Beckhoff AM8112 and controller EL7211-9014). In

C. Active projectile: wave mechanics and center of mass motion

Here we detail the analysis of the active projectile under impact. When the projectile makes contact with the wall, it is deformed under its own inertia into an elliptical shape, which represents a shearing of the enclosed 2D area. We model the collision by treating the elliptical deformation as as a second harmonic distortion of a periodic ring.

Non-Hermitian wave mechanics

In Section S1 of the S.I., we provide the nonlinear equations of motion for an idealized 1D chain of vertices of arbitrary shape. The simplest reduction of this nonlinear equation is to assume that the chain is initially straight, has periodic boundaries, and total length L. The linearized continuum equation of motion to order a 3 for the chain is given by

m∂ 2 t h + a 2 κ∂ 4 x h + 2a 3 κ a ∂ 5 x h + a 2 Γ∂ 4 x ∂ t h = 0 (M5)
where h is the transverse deflection of the vertices. It is useful to write the equations in dimensionless form:

∂ 2 τ h + ∂ 4 λ h + 2κ a ∂ 5 λ h + Γ∂ 4 λ ∂ τ h = 0 (M6)
where we have introduced the quantities

κa ≡ κ a N κ (M7) Γ ≡ Γ N 2 a √ κm (M8) ω 0 ≡ 1 N 2 κ ma 2 (M9) τ ≡tω 0 (M10) λ ≡ x L (M11)
where N = L/a is the number of unit cells. From Eq. (M6), the dispersion is given by:

ω± = - i Γ 2 q4 ± - Γ 2 2 q8 + q4 + 2iκ a q5 (M12)
where ω ≡ ω/ω 0 is the dimensionless frequency and q = qL is the dimensionless wave number.

In the experiment, N = 12, a = 7.5 cm, m = 0.2 kg, and κ a is the control parameter. We use κ = 20 mN m/rad and Γ = 4 mN m s/rad, determined as fitting parameters between theory and experiment.

In Fig. 3a and Extended Data Fig. 2, we show runs of the experiment for κ a = ±9.9 mN m/rad. The heat maps show the angular deviation δθ i for each vertex as a function of time. In the continuum model, the curvature c(x) = a 2 ∂ 2 x h is the coarse-grained variable corresponding to the local angular deflection δθ i . In inset of Fig. 3e, we show the curvature distribution corresponding to the solution:

h(x, t) = A e i( 4πx L -ω + t) + e i( 4πx L -ω -t) + c.c. ( M13
)
where A is an overall amplitude and c.c. stands for complex conjugate, and ω is evaluated at q = 4π/L. In Fig. 3e, we provide vertical, gray dashed lines where the imaginary part of the spectrum

ω Im = ω 0 • Im(ω) q=sgn(κ a )•4π (M14)
passes from positive to negative. When ω Im > 0 (ω Im < 0), the q = 4π mode is amplified (attenuated).

To quantitatively compare between model and experiment, we compute the real and imaginary parts of the lowest harmonic modes from the experimental data:

S 1 (t) = - 1 6 12 j=1 cos(jπ/3)δθ j (t) (M15) S 2 (t) = 1 6 12 j=1
sin(jπ/3)δθ j (t) (M16)

The deformations S 1 and S 2 over time are plotted in Fig. 3d. Geometrically, S 1 and S 2 correspond to independent shearing deformations of the area enclosed by the ring. We then define

A(t) = S 2 1 (t) + S 2 2 (t)
and ϕ(t) = arctan(S 2 /S 1 ) to be the polar coordinates (i.e. amplitude and angle) in Fig. 3d. The flux is defined to be

Ψ(t c ) = φe Ȧtc/A (M17)
where φ and Ȧ/A are averaged over all data points between the release from the wall and the end of the recorded data (which corresponds to the projectile exiting the field of view, or deformations large enough that self intersection between the vertices occurs). Here t c = 1 s is a characteristic time scale for the experiment. Notice that the flux quantity captures the unidirectionality of the wave via the prefactor φ as well as the amplification or attenuation via the exponent Ȧ/A. In Fig. 3e, data points denote weighted averages of all experiments and the error bars indicate standard deviation uncertainties arising from the particle tracking resolution and linear fits of the growth rate and phase velocity.

Center of mass motion

To capture the center of mass motion, we model the projectile as a square of side length R and we assume its only modes of deformation are the two shearing modes S 1 = u xxu yy and S 2 = 2u xy , where u ij is the strain tensor (assumed to be uniform throughout the square). The force on the center of mass of the block is proportional to the shear stress along the bottom surface. Hence we have:

  F x F y   = α 1   σ xy σ yy   =   σ 2 -σ 1   (M18)
where α 1 is a proportionality constant and σ 1 = σ yyσ xx and σ 2 = 2σ xy are the shear stresses corresponding to the stress tensor σ ij . The coupling κ a induces an asymmetric shear coupling parameterized by an overall magnitude α 2 and angle θ

  σ 1 σ 2   = α 2   cos θ sin θ -sin θ cos θ     S 1 S 2   (M19)
If we denote the center of mass of the block by (x, y), then we have S 1 = -2y/R and S 2 = 2x/R. Thus we obtain the equation of motion:

  ẍ ÿ  = 2α 1 α 2 LM   cos θ -sin θ sin θ cos θ     x y   (M20)
where M is the mass of the block.

In Fig. 3f, we plot the ratios p y,out /p y,in and p x,out /p y,in for the simple box model. These quantities are defined as follows. We integrate Eq. (M20) using as initial conditions (x, y) = 0 and ẏ = p y,in /(M Σ). Here Σ is a prefactor that sets the inelastic loss in the experiment that occurs upon impact. We then terminate the integration when y = 0 (meaning that that block is departing the floor) and define the outgoing momentum (p x,out , p y,out ). Since Eq. (M20) is linear, the quantities (p x,out /p y,in , p y,out /p y,in ) are independent of the numerical prefactors and only depend on the value of θ and Σ. The value of (p x,out /p y,in , p y,out /p y,in ) is easily determined numerically for all θ ∈ [-π/2, π/2]. For θ 1, one can also write the leading order expression:

p y,out p y,in =Σ + O θ 2 (M21) p x,out p y,in =Σθ + O θ 2 (M22)
To compare to experiments, the value Σ can be determined by setting the ratio of p y,out /p y,in at κ a /κ = 0. Finally, one needs to determine the relationship between θ and κ a /κ. Notice that the eigen-frequencies of Eq. (M20) are:

ω 2 ∝ 1 + i tan θ (M23)
Similarly, ignoring dissipation, the eigenvalues of the continuum model take the form:

ω2 ∝ 1 + i 8πκ a N κ (M24)
Since one cycle in shear space should correspond to one cycle in the vibration of the chain, we equate Re(ω)/ Im(ω) = Re(ω)/ Im(ω) and obtain

θ = arctan 8π N κ a κ (M25)
The box model only introduces one additional parameter Σ, which sets the overall inelastic loss upon impact but before rebound. By examining the rebound at κ a /κ = 0, we find Σ = 0.42. Here we derive the point response in the continuum limit of the wall. In the continuum, the internal forces are captured by the stress tensor σ ij and the geometric deformation is captured by the (linearized) strain tensor

The theoretical predictions are compared against experimental measurements in

u ij = 1 2 (∂ i u j + ∂ j u i )
, where u i is the displacement field. For small deformations, the stress-strain relationship can be approximated as linear σ ij = C ijmn u mn . The form of the elastic modulus tensor C ijmn can be deduced on the basis of symmetries and conservation laws. Since angular momentum is conserved (i.e. there are no external sources of angular momentum in the bulk of the system), we have C ijkl = C jikl . Moreover, since the internal forces are not triggered by solid body rotation, we have C ijkl = C ijlk . Since the microscopic structure has C 3 symmetry, the macroscopic elastic modulus tensor must be isotropic. These three conditions imply that the elastic modulus tensor takes the following form [START_REF] Scheibner | Odd elasticity[END_REF]:

C ijkl =Bδ ij δ kl + µ(δ il δ jk + δ ik δ jl -δ ij δ kl ) + K o E ijkl (M26)
where E ijkl = 1 2 [ il δ jk + ik δ jl + jl δ ik + jk δ il ] and where δ ij and ij are the Kronecker and Levi-Civita symbols, respectively. Here B and µ are the well known bulk and shear moduli, respectively, while the parameter K o is an example of an odd elastic modulus. In a pictoral notation, we may write the stress and strain in the following basis:

=σ xx + σ yy =u xx + u yy (M27) =σ xx -σ yy =u xx -u yy (M28) =σ xy + σ yx =u xy + u yx (M29)
Here, positive is equal to negative S 2 and is equal to positive S 1 in the notation of the main text. Using this notation, Eq. (M26) may be written as:

          = 2      B 0 0 0 µ K o 0 -K o µ                (M30)
It is useful to cast the equations into an overall stiffness E, known as the Young's modulus, as well as two dimensionless parameters: the Poisson's ratio ν and the so-called odd ratio ν o .

These parameters are given by:

E = 4B[µ 2 + (K o ) 2 ] µ(µ + B) + (K o ) 2 → 4µ 1 + K o µ 2 (M31) ν = µ(B -µ) -(K o ) 2 µ(B + µ) + (K o ) 2 → 1 (M32) ν o = BK o µ(B + µ) + (K o ) 2 → K o µ (M33)
The honeycomb lattice is approximately incompressible (to linear order). Hence in Eq. (M31-M33), we show the limiting expressions as B → ∞.

To compute the response to a point force, we must solve the bulk equation for static equilibrium ∂ i σ ij = 0. To do so, it is useful to introduce the airy stress function χ defined via

σ ij = il jk ∂ l ∂ k χ. Then ∂ i σ ij = 0 becomes ∆ 2 χ = 0.
We consider a point force of magnitude F on the boundary of a semi-infinite medium with stress free boundaries. Then the airy stress function is given by [START_REF] Landau | Theory of Elasticity[END_REF]:

χ = - F π rφ cos φ (M34)
where φ is the angle with respect to the direction of the applied force. The corresponding stress is: σ rr = -2F πr cos φ and σ rφ = σ φφ = 0. The strain is obtained by inverting Eq. (M26):

u ij = 1 E (1 -ν)δ ij δ mn + (1 + ν)(δ im δ jn + δ in δ jm -δ ij δ mn ) -2ν o E ijmn σ mn (M35)
Eq. (M35) yields the differential equations:

u rr =∂ r u r = - 4F πEr cos φ (M36) u φφ = 1 r ∂ φ u φ + u r r = ν 4F πEr cos φ (M37) 2u rφ =∂ r u φ - u φ r + 1 r ∂ φ u r = -ν o 8F πEr cos φ (M38)
The solution to these equations is given by:

u r = 4F πE -log r cos φ + ν -1 2 sin φ -ν o (sin φ + φ cos φ) (M39) u φ = 4F πE 1 + ν + 2 log r 2 sin φ + ν -1 2 φ cos φ + ν o φ sin φ (M40)
Upon setting ν = 1, as it is appropriate for the incompressible limit, we obtain the expressions

u r = 4F πE [-log r cos φ -ν o (sin φ + φ cos φ)] (M41)
u φ = 4F πE [(1 + log r) sin φ + ν o φ sin φ] (M42)
The the y-component of displacement is given by u y = -u r cos φ + u φ sin φ, yielding Eq. ( 5)

of the main text.

Rayleigh wave dispersion

Here we derive the Rayleigh wave dispersion shown in Fig. 4e. The bulk equations of motion reads

ρ∂ 2 t u j = ∂ i σ ij (M43)
We assume that the medium occupies the y < 0 half plane, which implies the boundary conditions

σ yi (y = 0) = 0 and u i (y → -∞) = 0 (M44)
where u i is the displacement field. The stress tensor σ ij is given by σ ij = C ijmn ∂ m u n . It is useful to parameterize the elastic modulus tensor in terms of the Young's modulus E, Poisson's ratio ν, and odd ratio ν o . We can then non-dimensionalize the equations of motion by defining

r ≡ k y /k x and kx = k x E ρω 2 ,
where k is the wave-vector. In this notation, Eq. (M43) takes the form:

  u x u y   = k2 x 2(1 -ν)   (r 2 +1)(1-ν 2 ) (ν+1) 2 +4(ν o ) 2 + 1 2(r 2 +1)(1-ν)ν o (ν+1) 2 +4(ν o ) 2 + r -2(r 2 +1)(1-ν)ν o (ν+1) 2 +4(ν o ) 2 + r (r 2 +1)(1-ν 2 ) (ν+1) 2 +4(ν o ) 2 + r 2   M (r, kx)   u x u y   (M45)
Hence, we obtain a secular equation

0 = det M (r, kx ) -1 (M46)
which is equivalent to

0 = (r 2 + 1) 2 k4 x -[2(1 -ν 2 ) + (ν + 1) 2 + 4(ν o ) 2 ](r 2 + 1) k2 x + 2(1 -ν)[(ν + 1) 2 + 4(ν o ) 2 ] (M47)
The solutions to Eq. (M46) take the form:

r ± 1 = ± 2(1 -ν 2 ) + (ν + 1) 2 + 4(ν o ) 2 -(1 + ν) 4 -8(1 -6ν + ν 2 )(ν o ) 2 + 16(ν o ) 4 2 k2 x -1 (M48) r ± 2 = ± 2(1 -ν 2 ) + (ν + 1) 2 + 4(ν o ) 2 + (1 + ν) 4 -8(1 -6ν + ν 2 )(ν o ) 2 + 16(ν o ) 4 2 k2 x -1 (M49)
Each solution comes with a normalized vector n i (r, kx ) ∈ ker[M (r, kx )-1]. Each n i represents a candidate mode for a Rayleigh wave. We are looking for linear combinations of candidate modes that satisfy Eq. (M44). Moreover, we require Im(r kx ) < 0 for all modes in the superposition in order to ensure that each mode decays into the bulk. The condition Im(r kx ) < 0 implies that linear combinations are of the form:

n ++ i ( kx ) ≡a 1 n i (r + 1 , kx ) + a 2 n i (r + 2 , kx ) (M50) n +- i ( kx ) ≡a 1 n i (r + 1 , kx ) + a 2 n i (r - 2 , kx ) (M51) n -+ i ( kx ) ≡a 1 n i (r - 1 , kx ) + a 2 n i (r + 2 , kx ) (M52) n -- i ( kx ) ≡a 1 n i (r - 1 , kx ) + a 2 n i (r - 2 , kx ) (M53)
We can express the stress at the boundary by:

  σ yx σ yy   =   S ±± x1 ( kx ) S ±± x2 ( kx ) S ±± y1 ( kx ) S ±± y2 ( kx )     a 1 a 2   (M54)
where we have introduced the matrix S ±± ( kx ) that is determined by computing

σ yj = C ijmn k i k m n ±± n .
Finally, the stress free boundary condition Eq. (M44) can be expressed algebraically as:

0 = g ±± ( kx ) ≡ det S ±± ( kx ) (M55)
For the special case that ν = 1, which corresponds to an incompressible solid, the algebraic expressions may be written explicitly. One finds:

r ± 1 = ± i (M56) r ± 2 = ± 4 + 4(ν o ) 2 -k2 x k2 x (M57)
Moreover, g ±± up to an unimportant denominator reads

g ++ =0 =⇒ 0 = 4[1 + (ν o ) 2 ] + k4 x -2(iν o + 2) k2 x -(2iν o kx + k3 x ) k2 x -4[1 + (ν o ) 2 ] (M58) g +-=0 =⇒ 0 = 4[1 + (ν o ) 2 ] + k4 x -2(iν o + 2) k2 x + (2iν o kx + k3 x ) k2 x -4[1 + (ν o ) 2 ] (M59) g --=0 =⇒ 0 = 4[1 + (ν o ) 2 ] + k4 x + 2(iν o -2) k2 x + (2iν o kx -k3 x ) k2 x -4[1 + (ν o ) 2 ] (M60) g -+ =0 =⇒ 0 = 4[1 + (ν o ) 2 ] + k4 x + 2(iν o -2) k2 x -(2iν o kx -k3 x ) k2 x -4[1 + (ν o ) 2 ] (M61)
Notice that g ++ = 0 and g +-= 0 are both equivalent to the condition

k6 x -(6 -4iν o ) k4 x + 4(1 -iν o )(2 -iν o ) k2 x -4(1 -iν o ) 2 = 0 (M62)
while g --= 0 and g -+ are both equivalent to the condition

k6 x -(6 + 4iν o ) k4 x + 4(1 + iν o )(2 + iν o ) k2 x -4(1 + iν o ) 2 = 0 (M63)
Inspection of the roots of Eq. (M62) and Eq. (M63) along with the condition Im(r kx ) < 0 reveals that the roots of Eq. (M62) should be taken for right traveling waves (Re( kx ) > 0) and the roots of Eq. (M63) should be taken for left traveling waves (Re( ky ) < 0). In the main text, we define q x = -Im(k x ).

Experimental strains and vibrations

In Fig. 4, we measure the local displacement field u i of each node. The colormap overlaying the reconstruction in Panels a and c shows u y averaged over each hexagon. The space-time heat maps show the average value of u y along averaged along vertical slices. For Panel d, we compute the strain tensor u ij = 1 2 (∂ i u j + ∂ j u i ) from the displacement field and average over all space. Shear 1 is defined as S 1 = (u xy + u yx )/2 and Shear 2 is defined as S 2 = (u yyu xx )/2.

In Panel e, the center of vibration is defined as

= x Vt [u y (x, t)] x Vt [u y (x, t)] x (M64)
where the variance V is over time and the average is over space. In these experiments, the value of the passive bond stiffness is taken at the calibrated value κ = 48 mN m/rad, see Section S2.

Equations (S2, S3, S4) constitute a closed system of equations.

Next, we take the continuum limit of Eq. (S2). We let x become a function of a continuous parameter x(λ) where λ ∈ [0, aN ] and approximate the differences by a Taylor expansion:

x α+1 -x α = ax (α) + a 2 2 x (α) + . . . (S5)
Expanding the discrete equations of motion to order a 3 yields

ẍ = - • x |x | 2 a 2 κ x × x (4) |x | 2 + 3 x × x |x | 2 + 2 (x × x ) 3 |x | 6 + 6 (x • x )[(x • x )(x × x ) -x × x ] |x | 8 + 2a 3 κ a x × x (5) |x | 2 -4 (x • x )(x × x (4) ) + (x × x )(x • x (4) ) |x | 4 -6 (x × x )(x • x ) |x | 4 + 12 [(x • x ) 2 -(x × x ) 2 ](x × x ) + 2(x • x )(x • x )(x × x ) |x | 6 -24 (x × x )(x • x )[(x • x ) 2 -(x × x ) 2 ] |x | 8 - • x |x | 2 -2 • x (x • x ) |x | 4 a 2 κ x × x |x | 2 -2 (x × x )(x • x ) |x | 4 + 2a 3 κ a x × x (4) |x | 2 + 3 x × x |x | 2 + 2 (x × x ) 3 |x | 6 + 6 (x • x )[(x • x )(x × x ) -x × x ] |x | 8 + ka 2 x (x • x ) |x | 2 - x |x | 3 (1 -|x |) + O a 4 (S6)
To simplify the equations of motion, we take k → ∞, which imposes the constrain equation |x | = 1. In this case, we obtain:

ẍ = -• x a 2 κ x × x (4) + 3x × x + 2(x × x ) 3 + 2a 3 κ a x × x (5) -4(x × x )(x • x (4) ) -12(x × x ) 2 (x × x ) -• x a 2 κ(x × x ) + 2a 3 κ a x × x (4) + 3x × x + 2(x × x ) 3 + d dλ (x T ) (S7)
where T (λ) is a tension field that enforces the constraint |x | = 1. (The field T is analogous to the pressure field p in an incompressible fluid). In principle, one must solve for T (λ) to produce a nonlocal equation entirely in terms of x and its derivatives. For example, let us denote the coefficient of • x in Eq. (S7) as g(λ) and likewise the coefficient of x as f (λ). Then the

condition |x | 2 = 1 implies 0 =| ẋ | 2 + ẍ • x (S8)
When Eq. (S8) is combined with Eq. (S7), we obtain a differential equation for T :

| ẋ | 2 -(x × x )(g + f ) -(x × x )f + (x × x ) 2 T + T = 0 (S9)
In principle, we can represent the solutions of Eq. (S9) as a non-local kernel T ( ẋ , x , x , . . . )

and place the nonlocal kernel into Eq. (S7) to obtain a closed equation of motion. We note that in deriving Eq. (S6) we used

θ i = √ 2 1 - (x α+1 -x α ) • (x α-1 -x α ) |x α+1 -x α ||x α-1 -x α | (S10)
instead Eq. (S4) in order to simplify the expansion. Moreover, for simplicity, we considered the dissipation free case in the derivation above, but the role of hinge dissipation Γ may be added in a similar fashion.

Geometry 1: straight line

Since we have the full nonlinear equations, we may linearize them around various background geometries. We will do this in two cases: an infinite straight line and a circle. For the straight line, we write:

x(λ) =   λ h(λ)   (S11)
We note that |x | = 1 + O(h 2 ), so the ansatz is consistent with the inextensibility condition.

The linearization of Eq. (S7) yields: ḧ = a 2 κh (4) + 2a 3 κ a h (5) (S12) Fig. S1. An illustration of the q = 0, 1, 2, 3 modes on a circular chain. Notice that q = 0 is a solid body rotation and q = 1 is a solid body translation, and hence are exact zero modes of the system. The q = 2 mode is the lowest order deformation.

as written in Eq. ( 4). The spectrum is

ω 2 = a 2 κq 4 + 2ia 3 κ a q 5 (S13)
which has positive and negative imaginary branches. We note that Eq. (S9) is satisfied by T = 0, so the tension does not enter the straight-line problem to leading order h.

Geometry 2: circle

Next, we perform a linearization about a circular geometry. For simplicity of notation, we will measure all distances in terms of the radius of the circle R, e.g. x → x/R and λ → λ/R.

Without loss of generality, we may write x as:

x =   cos λ sin λ   [1 + h(λ)] +   -sin λ cos λ   g(λ) (S14)
To leading order in g and h, requiring |x | = 1 yields the requirement h = -g . Thus we write

x = r(λ)[1 -g (λ)] + φ(λ)g(λ) (S15)
where r(λ) = (cos λ, sin λ) and φ = r . The equations of motion become:

g =φ • ẍ = A φ + T (S16) -g =r • ẍ = A r -T (S17)
where A r =12κ a gκg + 10κ a g (4) κg (5) -2κ a g (6) (S18)

A φ = -κg -2κ a g -κg (4) -2κ a g (5) (S19)
We can use Eqs. (S16-S17) to eliminate the tension, and we obtain: 5) a 2 κg (6) -2a 3 κ a g (7) (S20)

(1 -∂ 2 λ )g = A r + A φ = -a 2 κg + 10a 3 κ a g -2a 2 κg (4) + 8a 3 κ a g (
Hence, the spectrum is given by: ω 2 = a 2 κ(q 2 -2q 4 + q 6 )aiκ a (10q 3 -8q 5 -2q 7 ) 1 + q 2 (S21)

where q takes integer values since g(λ) must be periodic under λ → λ + 2π. As shown in Fig. S1, q = 0 corresponds to a solid body rotation and q = 1 corresponds to a solid body translation. Hence we obtain ω 2 = 0 for q = 0 and q = 1. The q = 2 is the lowest geometric deformation of the object, and we obtain a finite-frequency response. Notice that the expression for q = 2 matches exactly that of the straight chain. Finally, for q 1, we recover the result for the 1D chain ω 2 = a 2 κq 4 + 2ia 3 κ a q 5 , as expected.

B. Non-conservative forces and Hodge decomposition

Here we review how the notions of a nonconservative work cycle generalize to more arbitrary degrees of freedom. First we consider a system with N discrete degrees of freedom, which we can represent as a manifold M of dimension N . Recall that the tangent space of T p M at a point p ∈ M is a vector space of all derivatives evaluated at point p. The cotangent space T * p M is a vector space consisting of all linear maps from T p M to the real line. Given a local set of coordinates (θ may be identified through its coordinates λ α is the basis expansion λ = λ α dθ α . Likewise for vectors λ ∈ T p M may be written as λ = λ α ∂ ∂θ α . For mechanical systems, we are particularly interested in a distinguished vector field τ = τ α dθ α that represents the configuration-dependent forces. This force field is defined by the property that for any trajectory γ : [t i , t f ] → M, the quantity P = τ [ γ] represents the physical power exerted by the given forces. (Recall that γ t = dθ α dt t ∂ ∂θ α ∈ T γ(t) M.) The total work done over the trajectory γ is given by:

W [γ] = t f t i P (t)dt = γ τ (S22)
Vector fields such as τ can be classified via their behavior under integration. This classification scheme, known as Hodge decomposition, relies on two distinct notions of a derivative. Recall that an arbitrary differential k-form ω may be written as

ω = α 1 <•••<α k ω α 1 ...α k dθ α 1 ∧ • • • ∧ dθ α k . ( S23
)
where ∧ is the wedge product. The coefficients ω α 1 ...α k are function of the point p ∈ M. The exterior derivative d is a map from k forms to k + 1 forms given by dω =

α 1 <•••<α k γ ∂ω α 1 ...α k ∂θ γ dθ γ ∧ dθ α 1 ∧ • • • ∧ dθ α k . (S24)
Next, the Hodge dual, denoted by , is a vector space isomorphism between differential k forms and Nk forms. Its definition is most succinctly stated on a single basis element:

dθ α 1 ∧ • • • ∧ dθ α k = √ g (N -k)! α 1 •••α k β 1 ...β N -k dθ β 1 ∧ • • • ∧ dθ β N -k (S25)
forces. As alluded to before, the Hodge decomposition theorem allows us to classify vector fields τ based on their integral properties. Suppose that the trajectory γ(t) forms a closed loop that is also contractible to a point (i.e. it is a topologically trivial curve in parameter space). In this case, we may write γ = ∂V. Then we have

W = γ τ = V dτ = V dδA (S29)
where we used the fact that ddV = dh = 0. However, if γ enclosed a non-contractible loop (i.e.

it encloses a whole in parameter space), then h can contribute to the work done. In this sense h can be interpreted as the force due to a multi-valued potential. In certain mechanical contexts, the vertices evolve according to a first order equation of motion. To apply the uniaxial compression, the y (vertical) coordinate of the top and bottom row of particles is controlled, while the x coordinate is allowed to evolve according to the internal forces, thereby creating sliding boundary conditions. Upon compression, the strain data was sampled over all hexagons. When the tilt u xy is plotted as a function of u yy , the odd ratio is given by ν o = -duxy duyy evaluated at u yy = 0 (Fig. S3B). In Fig. S3C, the odd and Poisson ratio are plotted as a function of κ a /κ.

To find the odd modulus K o = dσxy duyy , the transverse stress σ xy was calculated from force data of a compression simulation with top and bottom row particles constrained in both directions, and its slope evaluated at u yy = 0 (Fig. -S3D). The code for the simulations is uploaded on the public repository (see Main Text for link).

FIG. 1 .FIG. 2 .FIG. 3 .FIG. 4 .

 1234 FIG. 1. Odd forces induce nonlinear work generating limit cycles (a) Three rigid linkages are connected by motorized vertices. (b) A hand pushes in on the left, and the right vertex contracts. (c) A hand pushes in on the right, and the left vertex expands.(d-f) Schematics summarizing the asymmetric (or odd) stiffness for κ a > 0. (g) The nonconservative force is proportional to the angular deflection δθ for small amplitudes but saturates at a value τ max at large amplitudes. (h-i) The dynamics are parameterized by the dimensionless quantity ξ = κ a /Γ I/κ capturing the strength of the odd forces. For |ξ| < 1 the system relaxes to its rest configuration. For |ξ| > 1, the odd stiffness, non-linearity, and dissipation conspire to produce a limit cycle at finite amplitude. Color indicates phase angle. See also Supplemental Video 1. The distance between two vertices is 7.5 cm.

  I. EXTENDED DATA FIGURESExtended Data Fig.1. Locomotion over complex terrain. The odd wheel from Fig.2climbs over a complex terrain. The line tracks the center of mass and the color bar is the same as in Fig. 2. See Methods for experimental details. See also Supplemental Video 2. Extended Data Fig. 2. Impact of an odd projectile against passive wall (a) An odd projectile with κ a < 0 is vertically incident against a wall. Scale bar: 10 cm. (b-c) The angular deformation of each vertex is plotted as a function of time, revealing a unidirectional wave propagating in the opposite direction of Fig. 3bc. See also Supplemental Video 3. Extended Data Fig. 3. Impact of a passive projectile on an odd wall (a) A passive projectile is launched at a wall of 120 motorized vertices with κ a < 0. Scale bar: 10 cm. (b) Corresponding reconstructions of the wall, where the vertical displacement u y of each hexagon is plotted. The total length is L = 128 cm. (c) The displacement is averaged over each column u y and is shown as a function of time and of the horizontal coordinate x. The deformation over time is opposite to the one observed in Fig. 4. See also Supplemental Video 4.

  microcontroller and calibrated by measuring torque-displacement slopes for different values of the electronic feedback. Each arm has a length of a = 7.5 cm. The coreless motor saturates at a maximum torque of τ max = 12 mN m. Adjacent vertices are also connected by rubber bands (blue in Fig.1a) which provide a passive elastic torsional stiffness denoted by κ = 48 mN m/rad. See Section S2 of the S.I. for calibration details.

B

  . Experimental protocolExperiments take place on top of a custom-made low friction air table. The table consists of two 1.5 m × 1.5 m plexiglas plates that sandwich 5 mm-wide air channels. The top plate is pierced by an array of holes (pitch 10 mm, diameter 1 mm) conveying air pressurized at 10 bars.Each motorized vertex floats on a thin layer of pressurized air without making contact with the table. The experiments are filmed via a Nikon D5600 camera equipped with a 50 mm lens, recording 2 Mpx images at 50 frames per second. A white marker is placed on the top of each vertex, and we use the python package OpenCV to detect and track the position of each vertex with a spatial resolution of 0.8 mm.

Fig. 3 and

 3 Fig. 3 and Fig. 4, the value of κ a assumes values in the range [-22, 22] mN m/rad. As detailed in Section S2 of the S.I., simulations are used to determine the odd ratio ν 0 as a function of κ a /κ.

  Fig. 3k. Here, the markers indicate averages over all runs of the experiment and the error-bars indicate the full range. D. Impact of an odd wall 1. Odd elasticity and response to a point force

∑

  Fig. S2. (A) The restoring torque τ measured over a range of twist angles δθ. The data shown is for rest angle θ 0 = 5π/6 (blue points) and θ 0 = 2π/3 (green points). The slope of the linear regression yields a value κ = 48 mN m/rad. (B) For different values of the active feedback, the motorized torque is approximately linear up to saturation of the DC motor at τ max ≈ 12 mN m. The parameter κ a is taken to be the slope of the line connecting the saturation points. (C) The decaying oscillations of the building block were fitted to the solution of a damped harmonic oscillator indicating a quality factor of Q = 2.6.

  1 , . . . , θ N ) on M in a neighborhood containing p, we may use the coordinate derivatives { ∂ ∂θ 1 , . . . , ∂ ∂θ N } as a basis for T p M and define a basis {dθ 1 , . . . , dθ N } for T * p M via the relationship dθ α [ ∂ ∂θ β ] = δ α β . In this basis, a given covector λ ∈ T * p M

Acknowledgments.

We thank Daan Giesen, Tjeerd Weijers, Ronald Kortekaas and Kasper van Nieuwland for technical assistance and Michel Fruchart for valuable discussions. M.B. acknowledges funding from the Fonds de la Recherche Scientifique-FNRS. C.C. and J.V. acknowledge funding from the European Research Council under Grant Agreement No. 852587. C.S. and V.V. acknowledge support from the Simons Foundation, the Complex Dynamics and Systems Program of the Army Research Office under grant W911NF-19-1-0268, and the University of Chicago Materials Research Science and Engineering Center, which is funded by the National Science Foundation under Award No. DMR-2011854. C.S. acknowledges funding from the National Science Foundation Graduate Research Fellowship under Grant No. 1746045.

Author contributions. Code and data availability.

The data and analysis codes are available at https://doi.org/10.5281/zenodo.5502579 Supplementary Information.

Supplementary Information is available for this paper.

Supplementary Information for:

Limit cycles turn active matter into robots S1. MATHEMATICAL DERIVATIONS

A. Wave equation for 1D Chain

In this section, we derive the equations of motion for a 1D chain of motorized vertices. We consider a 1D chain of point particles of mass m = 1 located at positions x α , α = 0, 1, . . . , N -1. We assume each bond has a rest length a, and we denote the angular tension in each vertex α by τ α . Additionally, we will at first assume that the bonds are extensible and have a Hookean spring constant k. (We will subsequently take k → ∞ to obtain inextensibility constraints.)

The equations of motion are constructed via Newton's laws:

where is the Levi-Civita symbol. The first term in Eq. (S2) are the forces due to the angular tensions, while the second term is the Hookean resistance to stretching. The functions τ α (x)

determine the vertex bending properties. For the experiments, we model the vertices via

Here, θ 0 is the rest angle (assumed to be identical for all vertices) and θ α is the angle of the vertex at point α:

In Eq. (S25), the quantity g is the determinant of the metric g αβ on M and γ 1 •••γ N is the totally antisymmetric Levi-civita symbol (whose indices may be raised using the metric tensor). The metric g αβ does not have an inherent physical meaning in the current context-nonetheless, the decomposition theorem stated below applied unambiguously regardless of the (arbitrary) choice of metric. Notice that is a duality in the sense that

We may use to induce an inner product on k forms on a closed manifold. Given ω, λ ∈ Λ k and taking M to be closed, we define:

Finally, we use the Hodge start to introduce a second type of derivative, δ, known as the codifferential. Unlike d, which is a linear map from k → k + 1 forms, the codifferential is a map from k → k -1 forms. The codifferential is explicitly given by δ = -1 d , and it readily verified that ω, dλ = δω, λ for arbitrary k form ω and k -1 form λ. Hence δ is the adjoint of d with respect to the inner product •, • . We note that d and δ both have the property

Finally, we are prepared to state the Hodge decomposition theorem. The simplest formulation of the Hodge decomposition theorem states that an arbitrary k-form τ on a closed manifold M may be uniquely decomposed into three pieces: The motorized building blocks described in Section II A are modeled via an angular stiffness κ, active feedback κ a , and angular dissipation Γ. As shown in the inset of Fig. S2A, we determine κ a and κ by placing a single robotic vertex in a rotational testing device (Instron E3000 linear-torsional) equipped with a torsional load cell (angular resolution: 2.5 arcsecs; torque resolution: 1 mN m). To determine the passive stiffness κ provided by the rubber band, the torque was measured as a function of angular deviation. The value of κ is taken to be the slope of the linear fit. We find κ = 48 mN m/rad for both the hexagonal ring with rest angle θ 0 = 2π/3 and the twelve-sided projectile with rest angle θ 0 = 5π/6.

To calibrate the torque provided by the motors, we remove the rubber band and measure the torque as a function of applied voltage, see Fig. S2B. When implemented in the robotic devices, the applied voltage is proportional to the angular deviation δθ of the neighboring vertices. The programmed proportionality constant gives rise to different torque-angle curves. The value of κ a is taken to be the slope of the straight line connecting the saturation points of the motors.

To determine Γ, we place the rubber band on the linkage system, and we apply an angular perturbation to the free edge. We measure the angular displacement and extract Γ from the decay rate, see Fig. S2C.

B. Calibration of odd ratio

To connect the ratio κ a /κ to the macroscopic odd ratio ν o and odd modulus K o in Fig. 4, we perform numerical simulations of the 2D robotic medium. As depicted in Fig. S3A, a simulated wall is subjected to a uniaxial compression. To reach a force balanced configuration,


Didier Francisco J Valdés-Parada

A compact

teaching and research institutions in France or abroad, or from public or private research centers.


G Janodet

Ph Rivière

J.-M Lamet

V Rialland

L Tessé

A Soufiani

MODELLING RADIATIVE PROPERTIES OF GAS MIXTURES IN NONEQUILIBRIUM HIGH-ALTITUDE ROCKET PLUMES

This paper presents the development of a Statistical Narrow Band model (SNB) in a nonequilibrium vibrational state for HCl and CO molecules. The population densities of the energy levels are obtained by a multi-temperature approach to compute nonequilibrium Line By Line (LBL) spectra. The SNB parameters are obtained by fitting the curves of growth from the LBL approach by a least squares error minimization using a Newton method for pure Lorentz and Doppler broadening regimes. The model is tested in Voigt broadening regime using a mixing rule and agrees well with the LBL approach. Finally, spectral correlation problems between η σ /κ σ and κ σ , where η σ and κ σ are the emission and absorption coefficients respectively, have been exhibited for CO 2 .

INTRODUCTION

The Local Thermodynamic Equilibrium (LTE) assumption is no longer valid when internal relaxation time scales are higher or comparable to hydrodynamic ones in gas flows. It can occur, for example, in high-altitude plumes when the flow undergoes a strong expansion at the nozzle exit. The gas is then under a Non-Local Thermodynamic Equilibrium (NLTE) state, which affects the molecule's radiative properties and thus needs to be considered for plume radiation [START_REF] Binauld | CO 2 IR radiation modelling with a multi-temperature approach in flows under vibrational nonequilibrium[END_REF]. However, using a high-resolution spectral approach to calculate the radiative properties makes solving the radiative transfer equation (RTE) too costly in the 3D coupled radiative-fluid computational framework. This study investigates the radiative properties of gaseous species of high-altitude rocket engine plumes. In this case, the gaseous molecules that mainly contribute to radiation are CO 2 and H 2 O. However, the spectral contributions of CO and HCl are not negligible, as shown in Fig. 1. In such expanding flows, CO 2 , CO and HCl are in a nonequilibrium vibrational state [START_REF] Binauld | CO 2 IR radiation modelling with a multi-temperature approach in flows under vibrational nonequilibrium[END_REF]. For H 2 O, the vibrational energy levels thermalize quickly with the translational levels and can therefore be modeled in the LTE assumption [START_REF] Taylor | Survey of vibrational relaxation data for processes important in the CO 2 -N 2 laser system[END_REF]. In this paper, nonequilibrium modeling of the radiative properties of CO, HCl and CO 2 molecules is considered, but only the SNB model for diatomic molecules is presented here.

SPECTROSCOPIC DATA AND MULTI-TEMPERATURE MODEL

The emission and absorption coefficients for wavenumber σ can be expressed under the assumption of independent lines as :

η σ = u→l A ul 4π hcσ ul n u f ul (σ -σ ul ) , (1)
κ σ = u→l (n l B lu -n u B ul ) hσ ul f ul (σ -σ ul ) , (2)
where n u and n l are the populations of energy levels u and l of the radiating isotopologue, A ul , B ul and B lu are the Einstein coefficients for spontaneous emission, induced emission, and absorption, σ ul is the central transition wavenumber, and f ul (σ -σ ul ) is the spectral line profile (assumed to be the same for the three radiative mechanisms) of the transition u → l. Spectroscopic data were taken from HITEMP-2010 for CO 2 [START_REF] Rothman | the high-temperature molecular spectroscopic database[END_REF], HITEMP-2019 for CO [START_REF] Rothman | the high-temperature molecular spectroscopic database[END_REF] and HITRAN-2020 [START_REF] Gordon | The HITRAN2020 molecular spectroscopic database[END_REF] for HCl. The population densities of the energy levels k are expressed using the multi-temperature approach and are defined as follows for the diatomic molecules [START_REF] Nagnibeda | Non-equilibrium reacting gas flows[END_REF] :

n k,HCl/CO = n HCl/CO Ab I g k Z I (T, T v ) exp - E rot {k} k B T - E vib {k} k B T v , (3)
where n HCl/CO is the total population of HCl or CO respectively, Ab I is the abundance of the isotopologue I considered, g k is the degeneracy of level k. The partitioning scheme of the energy is composed of a rotational part E rot {k} and a vibrational part E vib {k} such as :

E vib {k} = E vib {v k , J k } = min J k E{v k , J k }, (4)
E rot {k} = E{v k , J k } -E vib {v k , J k }, (5)
where v k and J k are the vibrational and rotational quantum numbers used for identifying the energy level k of energy E{v k , J k } in HITRAN/HITEMP spectroscopic databases. Z I (T, T v ) is the two-temperature partition function of the isotopologue I and its expression is

Z I (T, T v ) = k g k exp - E rot {k} k B T - E vib {k} k B T v . ( 6
)
For CO 2 , the population of level k, using the three-temperature approach of [START_REF] Binauld | CO 2 IR radiation modelling with a multi-temperature approach in flows under vibrational nonequilibrium[END_REF], is given by

n k,CO 2 = n CO 2 Ab I g k Z I (T, T 12 , T 3 ) exp - E rot {k} k B T - E v 12 {k} k B T 12 - E v 3 {k} k B T 3 , (7)
where T , T 12 , T 3 are the temperature associated to, respectively, rotation-translation modes, coupled vibrational modes ν 1 , ν 2 and antisymmetric vibrational mode ν 3 . Details concerning the energy partitioning scheme (E rot {k}, E v12 {k}, E v3 {k}) and the modelling of the three-temperature partition function Z I may be found in [START_REF] Binauld | CO 2 IR radiation modelling with a multi-temperature approach in flows under vibrational nonequilibrium[END_REF].

The models and data presented in this section are used to perform LBL computations for HCl, CO, and CO 2 molecules for the construction of SNB models. For diatomic molecules, a line wing cutoff of 500 cm -1 and a resolution of 10 -2 cm -1 were chosen for Lorentz profiles and 10 cm -1 and 10 -3 cm -1 respectively for Doppler profiles. For CO 2 , a line wing cutoff of 50 cm -1 and a resolution of 10 -2 cm -1 were selected for Lorentz profiles and 10 cm -1 and 10 -3 cm -1 respectively for Doppler profiles.

SNB MODEL FORMULATION FOR DIATOMIC MOLECULES

This section aims to develop a RTE formulation under nonequilibrium conditions suitable for applying statistical narrow-band models. The intensity escaping from a column (z 1 , z 2 ), neglecting any incoming intensity in z 1 , may be expressed as

I σ (z 1 → z 2 ) = z 2 z 1 η σ (z)τ σ (z → z 2 )dz, (8)
where

τ σ (z → z 2 ) = exp(- z 2 z κ σ (z)dz) = 1 κ σ (z) ∂τ σ ∂z (z → z 2 ) is the transmissivity of the column (z, z 2 ).
In order to use SNB models for solving the RTE, we need to express the intensity averaged over the spectral band ∆σ. However, the emission coefficient η σ and the transmissivity τ σ are strongly correlated. To correctly account for self-correlation, the intensity averaged over the spectral band ∆σ is expressed as

I σ ∆σ (z 1 → z 2 ) = z 2 z 1 η σ κ σ (z) ∂τ σ ∂z (z → z 2 ) ∆σ dz, (9)
where η σ κ σ is the analog of Planck's function in LTE conditions.

To compute Eq.( 9) with SNB models, a further simplification is necessary, which is

I σ ∆σ (z 1 → z 2 ) ≈ z 2 z 1 η σ κ σ (z) ∆σ ∂τ σ ∆σ ∂z (z → z 2 )dz, (10)
assuming that η σ κ σ and ∂τ σ ∂z are not correlated. Under LTE conditions, this assumption is always true since Planck's function can be considered constant inside a narrow spectral band ∆σ. The most critical test of this decorrelation assumption is to compare η σ ∆σ with η σ /κ σ ∆σ × κ σ ∆σ . This comparison is performed for CO (Fig. 2) and for HCl (Fig. 3) and shows that η σ κ σ and κ σ are perfectly decorrelated for significant nonequilibrium cases. The random SNB model of Mayer and Goody [START_REF] Goody | Atmospheric radiation[END_REF] has been selected in order to build an accurate model of mean transmissivities τ ∆σ . This model is based on statistical assumptions concerning line positions and intensities within a narrow band of width ∆σ and allows to express the transmissivity of a uniform column of length L, averaged over ∆σ [START_REF] Goody | Atmospheric radiation[END_REF] and the Malkmus tailed inverse-exponential distribution [START_REF] Malkmus | Random Lorentz Band Model with Exponential-Tailed S -1 Line-Intensity Distribution Function[END_REF]. They are special cases (respectively for α = -1 et α = 0) of the generalized Malkmus distribution function of [START_REF] Ph | Generalized Malkmus line intensity distribution for CO 2 infrared radiation in Doppler broadening regime[END_REF].

τ ∆σ (L) = 1 ∆σ ∆σ exp(-κ σ L)dσ = exp(- W δ ), (11)
P α (S) = C α S S M S α exp - S S M -exp - rS S M , (12)
where C α is a normalization constant.

In this study, we use the Malkmus distribution in the Lorentz broadening regime leading to

W L δ = 2γ L δ L   1 + k L xpLδ L γ L -1   . (13)
In the Doppler broadening regime, P α distribution leads to the mean equivalent line width expression

W D δ = β D H α ( k D xpL β D ), (14)
with

H α (y) = 1 α √ π +∞ -∞ [(1 + y exp(-ξ 2 )) α -1]dξ for α ̸ = 0 (15) H 0 (y) = 1 √ π +∞ -∞ ln(1 + y exp(-ξ 2 ))dξ for α = 0. (16)
The parameters k L and k D are the reduced mean absorption coefficients and can be considered equal

(k D = k L = k = κ/(xp)
) with good accuracy in the significant spectral bands of the spectrum for radiative transfer, γ L represents the mean Lorentz half width at half maximum of the lines, δ L is a modified mean line spacing, β D characterizes the degree of line overlapping and x the molar fraction of the absorbing species.

For each broadening regime, parameters k have been obtained directly by averaging over each narrow band (∆σ = 25 cm -1 ) high-resolution absorption coefficients provided by LBL calculations. The parameters δ L and β D have been generated from adjustments enabling the corresponding analytical expression of W L/D /δ, as a function of the column length L, fits the logarithm of the mean transmissivity obtained with LBL calculations (called the curve of growth). This adjustment is made by minimizing the least squares error on transmissivities thanks to the Newton method. The Lorentz half width at half maximum γ L for CO has been taken from [START_REF] Ph | Updated band model parameters for H 2 O, CO 2 , CH 4 and CO radiation at high temperature[END_REF] and we used the estimated values from HITRAN-2020 for the coefficients of γ L,HCl = (p/p s ) 0.075x HCl (T s /T ) 0.92 + 0.011(1 -x HCl )(T s /T ) 0.5 . The parameters δ L and k are assumed independent of the pressure and are generated at p = 1bar for air-CO/HCl mixtures with x = 0.1. For the curve of growth adjustments in the Lorentz broadening regime, 20 transmissivities are prescribed from 0.02 to 0.95, and the associated column lengths are computed from these transmissivities. For the adjustments in Doppler broadening regime, 21 column lengths are provided from xpL = 5 Pa.m to xpL = 5000 Pa.m with a regular logarithmic step. To quantify the errors introduced by the SNB model, a mean quadratic error is defined as follows for each spectral band

s(σ) = 1 N T .N T v .N L   N T i=1 N Tv j=1 N L k=1 τ ∆σ LBL (σ, T i , T v,j , L k ) -τ ∆σ SN B (σ, T i , T v,j , L k ) 2   , ( 17
)
where N T , N T v , N L are the number of points of the variable T, T v , L respectively. This error estimator in both broadening regimes is presented in Fig. 4. and Fig. 5 for the diatomic molecules. For each molecule in the Doppler broadening regime, the α parameter was optimized by testing different values. It shows that the best α value on average across all temperature conditions is -1 for HCl and 0 for CO, which correspond to the exponential and the Malkmus distributions, respectively. The fit errors are maximum when performed in the molecule's absorption bands, which is expected since the curve of growth is linear for small optical thicknesses. The errors for Lorentz broadening regime appear higher than for the Doppler regime. However, transmissivities are prescribed in Lorentz broadening, which means that the adjustment is made for a wide range of optical thicknesses. Therefore, the errors displayed for Lorentz broadening involve an adjustment on all transmissivity ranges, unlike for Doppler broadening, where the adjustments are made on the optically thin or medium part of the curve of growth.

In the context of high-altitude rocket plumes, we have to consider the Voigt broadening regime for which the mean black equivalent line width is calculated by combining W L /δ and W D /δ expressions, according to

W V δ = kxpL 1 -Ω -1/2 with Ω = 1 - 1 kxpL W D δ 2 -2 + 1 - 1 kxpL W L δ 2 -2 -1 (18)
as proposed in [START_REF] Ludwig | Handbook of infrared radiation from combustion gases[END_REF].

To validate the SNB model in this regime, we compare the radiative intensity at the outlet of a uniform column of length L = 1000 m from the model and the LBL approach in a wide range of pressure. This allows us to check the model in all broadening regimes and ensure that the mixing formula does not downgrade the model's accuracy. The comparison in Fig. 6 and Fig. 7 shows a good agreement between the SNB model and the LBL approach. Interestingly, the SNB model produces larger errors for the emitted radiative intensity in some band wings because the pressure-independent assumption of k is no longer valid in these spectral regions. Indeed, the reduced absorption coefficient k becomes proportional to the pressure in the band wings where the absorption coefficient is primarily due to far line wings. However, the contribution to the radiative transfer of these spectral bands is negligible, so these errors do not affect the relevance of the model. The multi-temperature model thus takes into account three temperatures instead of two for CO and HCl, which considerably increases the complexity of the spectral correlations between η σ and κ σ . The same decorrelation test in Voigt broadening regime, as shown in Fig. 2 and3, was performed in Fig. 8 and showed an important correlation between η σ /κ σ and κ σ in the 2.7µm and 4.3µm bands for conditions where T 12 > T 3 and P < 1000 Pa. To solve this issue, a vibrational splitting of the spectroscopic database is in progress and will be presented in future communication.

CONCLUSION

A nonequilibrium vibrational SNB model for CO and HCl has been developed and presented in this paper. The calculation of population densities for these diatomic molecules is based on a twotemperature approach (T, T v ). Since Kirchoff's law is no longer valid under these conditions, it is necessary to check the decorrelation between η σ /κ σ and κ σ beforehand. SNB model parameters were then obtained by adjusting the curves of growth obtained from the LBL approach for Lorentz and Doppler broadening regimes. The quality of the adjustments was then discussed quantitatively using error statistics. Finally, the SNB model was validated for the Voigt broadening regime by comparing the radiative intensity escaping from uniform columns for a wide range of pressures. For CO 2 , the calculation of population densities is based on a three-temperature approach (T, T 12 , T 3 ). Nevertheless, unlike the diatomic molecules, decorrelation between η σ /κ σ and κ σ is no longer valid for crucial spectral bands. To overcome this issue, a vibrational splitting of the spectroscopic database is under development and will be presented in an incoming publication.

Figure 1 .

 1 Figure 1. Radiation intensity emitted by a 10 m long column of exhaust gases at equilibrium with T = 1000 K, p = 100 Pa, and typical molar fractions x H2O = 0.0539, x CO2 = 0.0136, x CO = 0.4397, x HCl = 0.0191. The high resolution intensities are averaged over 5 cm -1 for legibility reasons.

Figure 2 .Figure 3 .

 23 Figure 2. Comparison of η σ ∆σ and η σ /κ σ ∆σ × κ σ ∆σ for x = 0.1, p = 100 Pa, T = 1000 K, T v = 2000 K (left) and p = 100 Pa, T = 2000 K, T v = 1000 K (right) for CO

Figure 4 .Figure 5 .

 45 Figure 4. Quadratic mean error of the transmissivities in Doppler broadening regime along wavenumbers for CO (left) and HCl (right)

Figure 6 .

 6 Figure 6. Radiative intensity for a uniform column at L = 1000 m for x = 0.1, T = 1000 K, T v = 2000 K (left) and T = 2000 K, T v = 1000 K (right) for CO. The pressures are 1 Pa (black), 10 Pa (red), 100 Pa (green), 1000 Pa (blue), 10000 Pa (yellow).

Figure 7 .

 7 Figure 7. Radiative intensity for a uniform column at L = 1000 m for x = 0.1, T = 1000 K, T v = 2000 K (left) and T = 2000 K, T v = 1000 K (right) for HCl. The pressures are 1 Pa (black), 10 Pa (red), 100 Pa (green), 1000 Pa (blue), 10000 Pa (yellow).

Figure 8 .

 8 Figure 8. Spectral correlation tests for x = 0.5, p = 100 Pa, T = 300 K, T 12 = 1000 K, T 3 = 2000 K (left) and x = 0.5, P = 100 Pa, T = 2000 K, T 12 = 1000 K, T 3 = 300 K (right) for CO 2

ACKNOWLEDGMENTS

The authors acknowledge the financial support of this study from the French Agence de l'Innovation de Défense (AID). This work was granted access to the HPC resources of IDRIS under the allocation 2022-A0102B00209 attributed by GENCI (Grand Equipement National de Calcul Intensif).


A Souissi A-B

M Arghir

D Lasseux

L Amami

E Chaduiron

Analysis of carbon brush seals with long bristles

Keywords: brush seal, carbon fiber bristles, theoretical model, experimental validation

Introduction

The present work introduces the analysis of a dynamic seal employed in rotating machinery to limit gas leakage between high and low pressure zones. The seal is comprised of a carbon brush consisting of closely serrated, long bristles of small diameter. The free length of the bristles exceeds the radial clearance between the rotor and the stator in which they are installed. The brush flexes in the axial direction along the rotor and fills the radial clearance between the rotor and the stator. Thus, a porous barrier is formed, impeding the flow of leakage. The theoretical estimation of the leakage flow is particularly challenging due to the complexity of the porosity estimation. The behavior of brush depends not only on the density defined as number of bristles per circumferential length, but also by its deformation. This is accounted for by considering the bristles as beams that undergo large deformations. Contact interactions occur between the bristles and the rotor, as well as among the bristles themselves. Consequently, the bristles experience deformations not only due to contacts but also due to fluid forces acting upon them. By achieving a state of equilibrium between the flow forces and the deformed brush structure, a constant permeability is then established.

Objective and methods:

The objective of this study is the prediction of the leakage flow through the brush seal. A simplified model is established by considering only the contacts of a first row of bristles with the rotor and by neglecting the contacts between bristles. The distance between bristles varies between zero at their clamped end and an imposed spacing at their free end. A Darcy law describes the flow through the heterogeneous bristles. The permeability coefficients are estimated from the homogenized flow within representative cells orthogonal to and along the deformed bristles. The theoretical model predictions are compared with experimental results obtained on dedicated test rig.

Macroscopic Model and bristle deformation:

The bristles have very large aspect ratio between length (over 1.9mm) and diameter (5 µm) and are submitted to large non-linear deformations. A Kirchhoff beam theory is therefore employed for describing the deformed beam together with a finite element approach using third order Hermitian functions. A penalty formulation is used for taking into account the contact with the rotor. The first row of bristles is in contact with the rotor. The free end of the next rows is progressively shifted to larger radial values. The spacing between the free ends of the bristles is a parameter of approximation (Dp). The resulting form of the deformed characteristic brush section is shown in Figure1.

: A boundary fitted curvilinear coordinate system is then associated to the deformed bristles. The coordinate system follows the first upstream, the last downstream bristle and the local normal direction. The porosity then varies along the coordinate system from its lowest value at the clamped end to its largest value at its free end. Two local problems are then solved in each point of the brush for the homogenized flow along and normal to the bristle. Local permeabilities are then calculated from the curvilinear coordinate system aligned with the bristles toward a cartesian coordinate system. The Darcy equation is then solved in a domain encompassing the upstream zone, the porous brush and the downstream zone as shown in Figure 1. A test rig was also developed and is depicted in Figure 2. It consists of a rotor deprived of rotation and enables testing carbon brush seals with five different interferences between the bristles and the rotor. Measured leakage flow rates are depicted in Figure3. Comparisons of the Position 3 between the theoretical model and experimental results are depicted in Figure 4. The numerical values of the Dp parameter, ranging from 0.1 times bristle diameter to 0.15 times bristle diameter, show consistency with the experimental results. However, the variation of the experimental leakage flow rate cannot be closely matched by the approach for two primary reasons. First, the deformation of the bristles considers only the rotor interference and disregards the influence of fluid forces. Second, the present approach assumes constant distance between deformed bristles, which is proportional to Dp parameter. These are strong assumptions that lust be reviewed in future developments of the model.

Figure 1 :

 1 Figure 1: Macroscopic model domain

Figure 2 :Figure 3 :

 23 Figure 2: Test rig


D Lasseux

F J Valdés-Parada

Macroscopic dynamic capillary pressure for two-phase flow in porous media

Keywords: Two-phase flow, dynamic capillary pressure, upscaling

Introduction

A closed predictive model for the macroscopic pressure difference, often referred to as the macroscopic capillary pressure, for two-phase flow in homogeneous porous media is of major importance to form a complete set of governing macroscale equations. This has been a long-lasting subject of research during the past 40 years [START_REF] Hassanizadeh | Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries[END_REF][START_REF] Marle | From the pore scale to the macroscopic scale: equations governing multiphase fluid flow through porous media. Flow and transport in porous media[END_REF], and the majority of the existing models remain empirical at some level of their derivation. A new closed macroscopic dynamic capillary pressure equation is reported in this work on the basis of the existence of a representative (periodic) unit cell able to locally describe microscale momentum transport in the system. This is carried out with an adjoint method together with a Green's formulation, requiring no other simplifying assumption.

Upscaled dynamic capillary pressure

The general framework of the analysis is the incompressible, creeping flow of two Newtonian fluids α and β saturating a homogeneous porous medium whose rigid and spatially homogeneous solid skeleton is the σ-phase.

Pore-scale and adjoint problems

The pore-scale flow problem can be written in a periodic unit cell representative of the domain as follows (α = β, γ)

∇ • v α = 0, in V α , (1a)
0 = ∇ • T pα -∇P α , in V α , (1b)
v β = v γ , at A βγ , (1c)
n βγ • T pβ = n βγ • T pγ + n βγ p β β -p γ γ r βγ + 2Hγn βγ + ∇ s γ, at A βγ , (1d)
v α = 0, at A ασ , (1e)
v α | S - αi = v α | S + αi , i = 1, . . . , N, (1f)
(n α • T pα ) S - αi = -(n α • T pα ) S + αi , i = 1, . . . , N, (1g)
pα α = 0. (1h)
Here, v α is the velocity in the α-phase, in which the pressure, p α , is decomposed as

p α = p α α + pα , T pα = -p α I + µ α ∇v α + ∇v T α , ∇P α = ∇ p α α -ρ α b α , H
is the double mean curvature of the fluidfluid interface, A βγ , which unit normal directed from the β phase to the γ phase is n βγ , µ α and γ are the dynamic viscosity and interfacial tension and ∇ s γ stands for possible Marangoni effects; S ± αi are the entrances/exits of the α-phase at the ith edge of the periodic unit cell.

For the derivation of the upscaled model, an adjoint problem is now considered that writes

∇ • f β = 1 V β , in V β , ∇ • f γ = - 1 V γ , in V γ , (2a)
0 = ∇ • T fα , in V α , (2b)
FICPM -Rueil-Malmaison, France -October 17-19, 2023.

f β = f γ , n βγ • T f β = n βγ • T fγ , at A βγ , (2c)
f α = 0, at A ασ , (2d)
f α | S - αi = f α | S + αi , (n α • T fα ) S - αi = -(n α • T fα ) S + αi , i = 1, . . . , N, (2e)
f α = 0, at r 0 α . (2f)

Macrosocpic capillary pressure

Following the procedure detailed in [START_REF] Lasseux | Upscaled dynamic capillary pressure for two-phase flow in porous media[END_REF], a Green's formula can be used to combine the physical flow problem and the adjoint (closure) problem that leads to the following macroscopic equation for the dynamic capillary pressure

∆P = p γ γ | xγ -p β β x β = -ψ β • ∇ p β β + ρ β b β • ϕ β -ψ γ • ∇ p γ γ + ρ γ b γ • ϕ γ + s βγ , (3)
where

ϕ α = Vα f α dV , ψ α = ϕ α - A βγ n ακ • f α z α dA, and s βγ = A βγ (2Hγn βγ + ∇ s γ) • f β dA.
Here, z α denotes the position of a point located at A βγ relative to the α-phase barycenter in the unit cell.

Results

This predictive equation is tested in the simple case of two-phase flow within a square pattern of parallel cylinders of circular cross section (see figure 1 in [START_REF] Lasseux | Upscaled dynamic capillary pressure for two-phase flow in porous media[END_REF]) with flow along a principal axis of the structure orthogonal to the cylinders axes. The wetting fluid flows under the form of sheets aligned with the mean flow and attached to the cylinders, whereas the non-wetting phase flows in the core, between the wetting fluid sheets sheets. The test consists in comparing the results on the dimensionless macroscopic dynamic capillary pressure ∆P * obtained from the solution of the pore-scale flow problem (DNS), on the one hand, and from the upscaled equation (3) (UE), on the other hand. Results are reported in figure 1 versus the wetting fluid saturation S β for a porosity ε = 0.8, two viscosity ratios µ * = µ β /µ γ and two values of the capillary number Ca = h 2 /γ, h being the pressure gradient magnitude in both phases and the periodic unit cell size. No volume force and no Marangoni effect are considered. The classical Laplace pressure P * cL = 2 H * βγ /Ca, where H * βγ = 1 A * βγ A βγ H * dA * and A * βγ is the dimensionless measure of A βγ , is also reported in these figures. These results show the excellent agreement between DNS and the uspcaled equation, and highlight the fact that the classical Laplace relationship cannot represent a satisfactory estimate of ∆P * . The influence of the pressure gradient in each phase (dynamic effects) is also pinpointed.

Figure 1 :

 1 Figure 1: * obtained from direct numerical simulation (DNS) and predicted from the upscaled equation (3) (UE), along with the interfacial term contribution, s * βγ , and classical Laplace capillary pressure, P * cL (see text), versus the wetting phase saturation, S β . a) µ * = 0.1, Ca = 10; b) µ * = 10, Ca = 0.1. ε = 0.8.