Skip to main content

Thank you for visiting You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Realization of active metamaterials with odd micropolar elasticity


Materials made from active, living, or robotic components can display emergent properties arising from local sensing and computation. Here, we realize a freestanding active metabeam with piezoelectric elements and electronic feed-forward control that gives rise to an odd micropolar elasticity absent in energy-conserving media. The non-reciprocal odd modulus enables bending and shearing cycles that convert electrical energy into mechanical work, and vice versa. The sign of this elastic modulus is linked to a non-Hermitian topological index that determines the localization of vibrational modes to sample boundaries. At finite frequency, we can also tune the phase angle of the active modulus to produce a direction-dependent bending modulus and control non-Hermitian vibrational properties. Our continuum approach, built on symmetries and conservation laws, could be exploited to design others systems such as synthetic biofilaments and membranes with feed-forward control loops.


Responsive materials in both biology and engineering distinguish themselves by their ability to respond to external stimuli in tailored ways1,2,3,4,5,6,7. For example, muscles contract in response to electrical signals8 and mechanocaloric solids undergo dramatic deformations in response to temperature changes9. Unlike computers or multicellular organisms with specialized functional components, the undifferentiated physical machinery available in a distributed material presents unique challenges for sensing, information processing, and response. Yet, the range of available functionalities is fundamentally extended when the materials possess distributed, local reservoirs of energy4,10,11,12,13,14,15,16,17,18,19. Such active materials exhibit responses not allowed by their passive, or energy conserving, counterparts.

Activity is intimately connected, but not equivalent, to a family of material symmetries collectively known as reciprocity20,21,22,23,24,25,26,27,28. Here we distinguish between two notions of reciprocity relevant for the design of mechanical metamaterials. The first notion is a generalization of Newton’s third law, which states that the forces between any two components of a mechanically isolated system must be equal and opposite. This notion of reciprocity can be formulated for other generalized momenta, such as angular momentum. Mechanical systems that violate this version of reciprocity must fundamentally be in contact with an external medium, such as a substrate or a background fluid, to act as a momentum sink. For systems described by a Lagrangian or Hamiltonian, translational (or rotational) symmetry gives rise to conservation of linear (or angular) momentum. Hence, systems with translational symmetry that violate this first notion of reciprocity must either be dissipative, driven, or active.

A second conceptually distinct notion of reciprocity is known as Maxwell-Betti reciprocity, which can be roughly defined as the symmetry between perturbation and response. When a mechanical system is deformed, the work done can schematically be written as \({{{\rm{d}}}}W=\mathop{\sum}\limits_{a}{\sigma }_{a}{{{{{\rm{d}}}}}}{u}_{a}\). Here, \({u}_{a}\) is a short hand notation for mechanical degrees of freedom such as displacements, rotations or strains, and \({\sigma }_{a}\) labels the conjugate forces, moments or stresses. For small perturbations about an undeformed state, we may write \({\sigma }_{a}={M}_{{ab}}{u}_{b}\). A generic medium is said to obey Maxwell-Betti reciprocity if and only if \({M}_{{ab}}\) is symmetric, i.e. \({M}_{{ab}}={M}_{{ba}}\). As long as the system’s linear response obeys Maxwell-Betti reciprocity, the forces can be derived from gradients of an energetic potential \({V=\frac{1}{2}M}_{{ab}}{u}_{a}{u}_{b}\). However, if the linear response of the medium violates Maxwell-Betti reciprocity, the internal energy is no longer a function of the coordinates \({u}_{a}\). In other words, the medium can perform nonzero work along a closed cycle of deformations27. Such a medium necessarily contains non-conservative forces that require an internal or external source of energy to be present.

While Maxwell-Betti reciprocity and momentum conservation are distinct, one can certainly design and build mechanical systems that harness linear or angular momentum sources to achieve violations of Maxwell-Betti reciprocity7,20,29. Such approaches, however, inherently involve mechanical coupling to an external medium. Here we take an alternative route. We explore the case where \({\sigma }_{a}\) represents shear and moment stresses, which conserve currents of linear and angular momentum. In this case \({u}_{a}\) represents geometric deformations and \({M}_{{ab}}\) represents the stiffness matrix containing all of the material’s elastic moduli. We refer to the antisymmetric components \(\left({M}_{{ab}}-{M}_{{ba}}\right)/2\) as odd elasticity27,30,31,32. A material displaying odd elasticity must violate Maxwell-Betti reciprocity, though it needs not rely on external sources of linear or angular momentum.

Recent advances in metamaterial design and prototyping have utilized active components to achieve functionalities such as sensing, lasing, and cloaking10,11,12,15, frequency dependent reflectivity33, unidirectional wave amplification7,34,35, energy harvesting16, and analog computation17. Nonetheless, all the active non-reciprocal metamaterials so-far realized exhibit either of the following fundamental limitations: the active non-reciprocal effects either vanish from the linear response in the quasistatic limit33 or they require the presence of background sources of linear or angular momentum7,29,36. As a result, their functionalities are largely restricted to finite-frequency control or fundamentally require the sample to be in contact with an additional medium that acts as a momentum sink or source.

Here, we report the design, construction, and experimental demonstration of a freestanding metamaterial whose elasticity is unattainable in passive media (Fig. 1a–d). The metamaterial is constructed with piezoelectric elements37,38,39,40,41,42,43,44,45,46 mounted on a beam and controlled by electrical circuits. Our approach enables an asymmetric coupling between bending and shearing in micropolar solids47,48. This results in an odd micropolar material that simultaneously breaks parity and Maxwell-Betti reciprocity. The spectrum of the metabeam exhibits a non-Hermitian topological index which results in the localization of vibrational modes at sample boundaries. We experimentally show the resulting unidirectional amplification/attenuation of waves propagating through the metamaterial. Our work sheds light on controlling non-reciprocal elasticity in artificial materials.

Fig. 1: Design and mechanics of an odd micropolar metabeam.

a A single unit cell featuring three piezoelectric patches mounted on a beam: one that acts as a sensor, and two that act as actuators. b A segment of the full metabeam. c Each unit cell has an electronic loop. The voltage Vs induced by the central piezoelectric is fed into a transfer function H(ω) = Va(ω)/Vs(ω) that sends opposing voltages Va and −Va to the piezoelectric actuators. d A photograph of the metabeam (horizontal) with the electronic circuits in the foreground. We note that the mechanical forces from the attached wires are negligible. The wires act only as sources of energy and computation, but not of linear or angular momentum. e The motion of the metabeam can be described by two independent fields, φ and h, which parameterize the angular and vertical displacements of the metabeam. Notice that under a reflection about the z-axis, we have φ → −φ and h → h. f When the beam bends, the center piezoelectric is stretched. g The antisymmetric electronic actuation then gives rise to a shearing stress proportional to the modulus P.


Design of the active metamaterial

The metamaterial we design takes the form of a thick beam, whose shape is characterized by two independent degrees of freedom (Fig. 1e): the height h(x) of the midplane and the angle φ(x) of the cross section with respect to the vertical. A single unit cell in the beam is equipped with three piezoelectric patches that enable shape-sensing and response. The central patch acts as a sensor that acquires a voltage proportional to the elongation or contraction of the top surface (Fig. 1f). The piezoelectric patches at the front and back of the unit cell serve as mechanical actuators that elongate and contract in response to an applied voltage49,50. A transfer function H(ω) processes the input voltage from the central patch and sends output signals to the actuating patches (Fig. 1c). For each unit cell, we implement H(ω) using a minimal electrical circuit (Fig. 1d). The electronics only couple piezoelectric patches within a single unit cell, thereby creating a control system that is both local and decentralized. The resulting system only needs to be connected to a voltage or, more generally, energy source. One practical advantage of this approach is that the voltage or power source can easily be housed inside the medium itself.

The active metamaterial we design is freestanding—it does not push or pull on an external medium. Indeed, it obeys Newton’s third law by preserving both angular and linear momentum as a traditional beam. However, the crucial difference between a traditional beam and the one we construct is the presence of internal energy sources. We design the feedback such that when the central patch experiences elongation or compression due to bending of the beam, the electronic loop produces output voltages that are antisymmetric (Fig. 1c), resulting in shear stresses (Fig. 1g). However, the ensuring shear strain does not stretch or compress the central piezoelectric patch. Therefore, the electromechanical control loop is entirely feed-forward: bending induces shear, while shear does not induce bending.

Odd micropolar elasticity

The feedback results in an elastic response that cannot be realized without an internal source of energy. This effect, apparent in the emergent continuum equations, may be deduced solely using symmetries and conservation laws based on classical beam theory (1D micropolar elasticity). Crucial to our design is the notion of a parity inversion P, defined here to be a mirror reflection of the beam about the z-axis that sends x to −x. Fig. 1e shows that under parity, the two independent degrees of freedom, h(x) and φ(x), transform as

$$h (x) \underrightarrow{{{\mathcal{P}}}} h (-x)$$
$${\varphi} \left(x\right)\underrightarrow{{{\mathcal{P}}}}-{\varphi} \left(-x\right)$$

Since φ(x) acquires a minus sign under parity, we say that φ(x) is a micropolar degree of freedom47,48. The equations of motion for a freestanding micropolar beam are then built out of conservation laws. Linear momentum conservation implies that

$$\rho \ddot{h}={\partial }_{x}{\sigma }_{{zx}}$$

where σzx is the shear stress and ρ is the mass density. Moreover, angular momentum conservation implies

$$I\ddot{\varphi }={\partial }_{x}M+{\sigma }_{{zx}}$$

where M is the bending moment and I is the cross-sectional moment of inertia. The moment M and stress σzx are themselves determined by the deformation of the beam via a set of constitutive relations. To leading order in gradients of h and φ, the internal geometry of the beam is approximated by two independent types of deformation: bending b(x) (Fig. 1f) and shearing s(x) (Fig. 1g), defined as

$$b\left(x\right)={\partial }_{x}\varphi$$
$$s\left(x\right)={\partial }_{x}h-\varphi$$

Under parity inversions, Eqs. (1, 2) imply that


Assuming a linear response, one may in general write the linear constitutive relations as:

$$\left[\begin{array}{c}{\sigma }_{{zx}}(t)\\ M(t)\end{array}\right]=\int_{-{{\infty }}}^{{{\infty }}}\left[\begin{array}{cc}{C}_{11}(t^{\prime} ) & {C}_{12}(t^{\prime} )\\ {C}_{21}(t^{\prime} ) & {C}_{22}(t^{\prime} )\end{array}\right]\left[\begin{array}{c}s\left(t-t^{\prime} \right)\\ b\left(t-t^{\prime} \right)\end{array}\right]{{{{{\rm{d}}}}}}t^{\prime}$$

or in terms of frequency:

$$\left[\begin{array}{c}{\sigma }_{{zx}}(\omega )\\ M(\omega )\end{array}\right]=\left[\begin{array}{cc}{C}_{11}(\omega ) & {C}_{12}(\omega )\\ {C}_{21}(\omega ) & {C}_{22}(\omega )\end{array}\right]\left[\begin{array}{c}s\left(\omega \right)\\ b\left(\omega \right)\end{array}\right]$$

We denote the matrix operator on the right-hand side of Eq. (9) as C(t). Causality implies that C(t) = 0 for t < 0 and reality of C(t) implies C(−ω) = C*(ω). Direct substitution of Eq. (10) into Eqs. (3, 4) reveals that the equations of motion are invariant under parity if and only if C21 and C12 are zero. Hence, we say that C21 and C12 are micropolar moduli.

It is useful to parameterize the C(ω) as

$$\left[\begin{array}{cc}{C}_{11}(\omega ) & {C}_{12}(\omega )\\ {C}_{21}(\omega ) & {C}_{22}(\omega )\end{array}\right]=\left[\begin{array}{cc}\mu (\omega ) & \alpha \left(\omega \right)+\beta (\omega )\\ \alpha \left(\omega \right)-\beta (\omega ) & B(\omega )\end{array}\right]$$

where μ is the shear modulus and B is the bending modulus, and α and β are the symmetric and antisymmetric components of the micropolar moduli. If the beam lacks an internal source of energy, then the total work done by the beam on any deformation process that begins and ends in the same configuration must be less than zero: ∆W ≤ 0. This energy condition places the following constraints on the finite frequency linear response coefficients (see “Methods”):

$$0\,\ge\, {{{{{\rm{Im}}}}}}\,[\mu \left(\omega \right)+B(\omega )]$$
$$0\,\le\, {{{{{\rm{Im}}}}}}\left[\mu \left(\omega \right)\right]{{{{{\rm{Im}}}}}}\left[B\left(\omega \right)\right]-{{{{{\rm{Im}}}}}}{\left[\alpha \left(\omega \right)\right]}^{2}-{{{{{\rm{Re}}}}}}{\left[\beta \left(\omega \right)\right]}^{2}$$

for all ω. Notice that the reality condition C(−ω) = C*(ω) implies Im[C] = 0 at ω = 0. As a consequence, we must have β(ω → 0) = 0 for any passive beams.

However, our active metamaterial has an internal source of energy and thus need not obey this constraint. The feed-forward coupling between bending and shearing suggests a linear response matrix of the form:

$$\left[\begin{array}{cc}{C}_{11}(\omega ) & {C}_{12}(\omega )\\ {C}_{21}(\omega ) & {C}_{22}(\omega )\end{array}\right]=\left[\begin{array}{cc}\mu (\omega ) & P(\omega )\\ 0 & B(\omega )\end{array}\right]$$

The electronic control loop introduces the coefficient P = 2α = 2β, which we refer to as the odd micropolar modulus. This modulus breaks two crucial symmetries. Since the electromechanical coupling violates parity, the active modulus P must occur in the off-diagonal entries. Moreover, P occurs only in the upper-right entry because the electro-mechanical coupling is feed-forward: bend causes shear, but shear does not cause bend. As a result, the matrix C is asymmetric, indicating that the beam violates Maxwell-Betti reciprocity, even at zero frequency.

For our metabeam, we measure the moduli via COMSOL simulations in which we apply controlled displacements at finite frequency to the front and back faces of a single unit cell. By measuring the reaction forces on these faces, we determine the resulting stresses and, consequently, the moduli. We empirically find that μ = 1.3 × 109 kg/ms2 and B = 0.112 × 106 kg/s2 are approximately independent of frequency, while P(ω) = ΠH(ω), where H(ω) is the transfer function and Π = 4.7 × 106 kgm/s2 is a material constant. From the metabeam geometry and materials, we compute the average volumetric density ρ = 5613 kg/m3and the average cross-sectional moment of inertia 5.9 × 10−3 kg/m. See Supplementary Notes 1 and 2 for additional characterization details.

Energy cycles

It is useful to consider the elastic limit of Eq. (11) in which C(ω) is real and approximately independent of frequency. For a passive beam in this limit, the matrix C is obtained by approximating the energetic cost of deformation by a quadratic function:

$$W=\frac{\mu }{2}{s}^{2}+\frac{B}{2}{b}^{2}+\alpha {sb}$$

The work (per unit volume) done in an infinitesimal deformation of the beam is given by:

$${{{{{\rm{d}}}}}}W={\sigma }_{{zx}}{{{{{\rm{d}}}}}}s+M{{{{{\rm{d}}}}}}b$$

From Eq. (16), we conclude \({\sigma }_{{zx}}=\frac{\partial W}{\partial s}\) and \(M=\frac{\partial W}{\partial b}\). Hence, we obtain the linear constitutive relations:

$$\left[\begin{array}{c}{\sigma }_{{zx}}\\ M\end{array}\right]=\left[\begin{array}{cc}\mu & \alpha \\ \alpha & B\end{array}\right]\left[\begin{array}{c}s\\ b\end{array}\right]$$

Since C in Eq. (17) is obtained via a second derivative of an energy function, it is symmetric. This is a manifestation of the Maxwell-Betti reciprocity theorem. In this case, the cumulative work done over a sequence of deformations depends only on the initial and final configurations: ∆W= ∫ dW = Wfinal − Winitial. In particular, if the procedure begins and ends at the same state, we have ∆W = 0.

Notice that the constitutive relation in Eq. (14) violates C = CT. Hence, the Maxwell-Betti reciprocity theorem implies that the form of C in Eq. (14) does not follow from a potential energy function. To see this explicitly, consider the differential of energy given by:

$${{{{{\rm{d}}}}}}W={\sigma }_{{zx}}{{{{{\rm{d}}}}}}s+M{{{{{\rm{d}}}}}}b$$
$$={{{{{\rm{d}}}}}}\left(\frac{\mu }{2}{s}^{2}+\frac{B}{2}{b}^{2}+\frac{P}{2}{sb}\right)+\frac{P}{2}\left(b{{{{{\rm{d}}}}}}s-s{{{{{\rm{d}}}}}}b\right)$$

Notice that the second term in Eq. (19) cannot be expressed as the differential of a potential. By Green’s theorem, integrating the work done dW over a closed loop in strain space yields:

$$\oint _{\partial V}{{{{{\rm{d}}}}}}W=\oint _{V}P{{{{{\rm{d}}}}}}s{{{{{\rm{d}}}}}}b=P\times {{{{{\rm{Ar}}}}}}{{{{{\rm{ea}}}}}}$$

where “Area” is the signed area of the region V enclosed by the path in the space of strains. Notice that the work done per cycle can be either positive or negative depending on the orientation of the path, and is independent of the rate of the process (in the quasistatic limit). Since the crucial ingredient for such cycles is an antisymmetry in C, we refer to this form of elasticity as “odd” (i.e., antisymmetric) elasticity 27.

To verify that our design displays this property at low frequencies, we perform COMSOL simulations of the beam with full piezoelectric coupling. As illustrated in Fig. 2a, we subject a single unit cell to a four-step protocol of shearing and bending by enforcing displacement-controlled boundary conditions at the two ends of the beam. We measure the reaction forces on the control surfaces to compute the work done by the beam, plotted in Fig. 2b. When the deformations are performed in a clockwise direction in strain space, as shown in the top panel of Fig. 2b, the cumulative work done is negative once the unit cell returns to its initial configuration. Hence, energy flows from the external agent into the internal power reserves of the medium. When the cycle is reversed, so is the flow of energy. The ability to inject or extract mechanical energy through quasistatic cycles is synonymous with odd elasticity, i.e. a quasistatic stress-strain relationship that does not follow from a potential energy. In the quasistatic limit, the total work done is given by P times the area enclosed in strain space. In the Methods, we also derive energy relations for cycles at finite frequency. For those cases, the total work done is related to |P|, arg(P), and the area enclosed in strain space.

Fig. 2: Quasistatic deformation cycles with odd micropolar elasticity.

a The state of the unit cell is tracked in the space of shear and bend. When a quasistatic closed path is traced out in this space, the unit cell performs work per unit volume that is proportional to the area times the modulus P. The z-displacement is provided in arbitrary units. b We numerically compute the work done for a clockwise (top) and a counterclockwise (bottom) path. The solid lines are predictions from the continuum theory, and the black dots result from finite element simulations of the unit cell. In the simulations, maximum amplitudes of bending and shearing are bmax ≈ 10−1 m−1 and smax ≈ 10−2, respectively. See Supplementary Note 1 for further details on the simulation.

Odd micropolar elastodynamics and the non-Hermitian skin effect

We now ask about the dynamic consequences of the active feed-forward control with parity violation. We first consider the elastic approximation in which C(ω) is real and frequency independent. The linearized continuum equations governing the motion of the beam are given by:

$$\rho \ddot{h}=\mu {\partial }_{x}^{2}h+P{\partial }_{x}^{2}\varphi -\mu {\partial }_{x}\varphi$$
$$I\ddot{\varphi }=\mu {\partial }_{x}h+B{\partial }_{x}^{2}\varphi +P{\partial }_{x}\varphi -\mu \varphi$$

Using Fourier transform, Eqs. (21, 22) may be cast in the form:

$$\omega \left[\begin{array}{c}{\tilde{p}}_{h}\\ {\tilde{p}}_{\varphi }\\ \tilde{s}\\ \tilde{b}\end{array}\right]=\underbrace{{\omega }_{1}\left[\begin{array}{cccc}0 & 0 &-k{l}_{1} & -\tilde{P}{l}_{1}k\\ 0 & 0 & i & i\tilde{P}-k{l}_{2}\\ -k{l}_{1} & -i & 0 & 0\\ 0 & -k{l}_{2} & 0 & 0\end{array}\right]} _{D(k)}\left[\begin{array}{c}{\tilde{p}}_{h}\\ {\tilde{p}}_{\varphi }\\ \tilde{s}\\ \tilde{b}\end{array}\right]$$

where k and ω are the wave number and frequency associated with the Bloch wave ei(kxωt). We have introduced the notation:

$$\widetilde{s}=\sqrt{\mu }\left({\partial }_{x}h-\varphi \right),\,\widetilde{b}=\sqrt{B}{\partial }_{x}\varphi$$
$${\widetilde{p}}_{\varphi }=\sqrt{I}\dot{\varphi },\,{\widetilde{p}}_{h}=\sqrt{\rho }\dot{h}$$

Here, \(\widetilde{s}\) represents the shear, \(\widetilde{b}\) represents the bending, \({\widetilde{p}}_{\varphi }\) is the angular momentum, and \({\widetilde{p}}_{h}\) is the z-component of the linear momentum. This parameterization is natural since the standard inner product

$$2e={\left|{\widetilde{p}}_{h}\right|}^{2}+{\left|{\widetilde{p}}_{\varphi }\right|}^{2}+{\left|\widetilde{s}\right|}^{2}+{\left|\widetilde{b}\right|}^{2}$$

is equal to twice the mechanical energy density e. The dynamical matrix D(k) depends on four parameters: \({l}_{1}\equiv \sqrt{I/\rho }\) is roughly the thickness of the metabeam; \({l}_{2}\equiv \sqrt{B/\mu }\) is the distance over which shearing and bending of equal transverse deflection cost equal amounts of energy; \({\omega }_{1}\equiv \sqrt{\mu /I}\) sets a frequency scale separating transverse flexural modes and high frequency shearing modes; finally, the parameter \(\widetilde{P}\equiv P/\sqrt{B\mu }\) is the normalized odd micropolar modulus. For our metabeam l1 ≈ 10−3 m, l2 ≈ 10−2 m, \({\omega }_{1}\,\)≈ 105 Hz, and \(\widetilde{P}\) 1.

Within the continuum theory, the vibrational dynamics can be captured by solving the secular equation detD(k)[− ω] = 0, which takes the form

$$0={\widetilde{\omega }}^{4}-\left[1-i\widetilde{P}k{l}_{2}+{k}^{2}\left({l}_{1}^{2}+{l}_{2}^{2}\right)\right]{\widetilde{\omega }}^{2}+{k}^{4}{l}_{1}^{2}{l}_{2}^{2}$$

where \(\widetilde{\omega }\) = ω/ω1. Notice that Eqs. (21, 22) are two coupled second order equations and hence permit a dispersion with four branches. For small wavenumber, the dispersion for the low frequency flexural bands is given by:

$${\omega }_{\pm }=\pm {\omega }_{1}\left[{l}_{1}{l}_{2}{k}^{2}\pm i\widetilde{P}{l}_{1}{l}_{2}^{2}{k}^{3}+{{{{{\mathscr{O}}}}}}\left({l}_{1}^{2}{l}_{2}^{2}{k}^{4}\right)\right]$$

As can be seen from Eq. (28), when P is nonzero and real, the periodic boundary spectrum acquires a nonzero imaginary component. The nonzero imaginary contribution arises since the active metabeam has the ability to physically introduce or remove mechanical energy. Moreover, the modulus P violates parity and hence breaks the symmetry between k → −k, in contrast to the design in ref. 7 where the parity violation is induced by a term in the dispersion that is linear in k.

The simultaneous breaking of parity and energy conservation allows our active micropolar metamaterial to selectively amplify and attenuate waves based on their direction of travel. In Fig. 3a, we show the spectrum of the flexural mode in the right half of the complex plane for P > 0. The solid line is the continuum theory, valid at small \(k\ll 1/\sqrt{{l}_{1}{l}_{2}}\), and the discrete points are the results of fully coupled COMSOL simulations. In the calculations, P = ±3Π, and there are no free fitting parameters. (The material constants Π, ρ, B, and μ in the continuum theory are determined from independent simulations, see Odd micropolar elasticity.) As illustrated in Fig. 3a, for P > 0, Im(ω) > 0 whenever Re(dω/dk) > 0 and Im(ω) < 0 whenever Re(dω/dk) < 0. Physically, this means that wave packets traveling to the right are amplified, while wave packets traveling to the left are attenuated. As an illustration, Fig. 3b shows the inverse penetration depth κ for ω taken to lie along the positive real axis. From the continuum theory, we compute the analytical formula for κ (see “Methods”)

$$\kappa =\sqrt{\frac{\rho }{B}}\frac{P}{4\mu }\omega$$
Fig. 3: Non-Hermitian skin effect via the odd micropolar elasticity.

a The vibrational spectrum for the flexural mode of a metabeam with periodic boundary conditions and odd micropolar modulus P = 3Π. The black line results from the continuum theory given by Eq. (28). The data points are obtained via fully piezoelectrically coupled simulations in COMSOL with the hue indicating the wavenumber kL, where L is the unit cell length. For the full spectrum plotted as a function of k in the continuum theory and in the numerics, see Fig. 7 and S2, respectively. The inset compares the continuum theory and simulations for small wavenumbers \(\ll 1/\sqrt{{l}_{1}{l}_{2}}\). b The inverse penetration depth κ for real ω in a medium with open boundary conditions. The points are the results of COMSOL simulations, the black lines are Eq. (29), and the dark lines are the result of the transfer matrix method, see Supplementary Note 3. c The localized states are connected to a topological index ν(ω). The periodic boundary spectrum for P > 0, P = 0, and P < 0 are represented schematically by the solid lines. The arrows indicate the direction of increasing k. For a given frequency ω, the winding number ν(ω) of the periodic boundary spectrum indicates the presence of a localized mode. d The localization of eigenmode at the value of ω denoted by the star in (c) is schematically illustrated.

In Fig. 3b, we compare Eq. (29) to the COMSOL simulations (see “Methods”) and a semi-analytical technique known as the transfer matrix method (see Supplementary Note 3).

This unidirectional amplification can be understood from the point of view of the non-Hermitian skin effect51,52,53,54,55,56,57,58,59,60,61,62. Given a complex frequency ω, we can define the following topological index:

$$\nu \left(\omega \right)=\frac{1}{2\pi i}\mathop{\sum}\limits_{\alpha }\oint _{-\pi /L}^{\pi /L}\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}k}{{{{{\rm{log }}}}}}\left[{\omega }_{\alpha }\left(k\right)-\omega \right]{{{{{\rm{d}}}}}}k$$

where L is the length of a single unit cell and ωα(k) is the frequency of the α band. Here, we take α to run over the flexural bands. The topological index ν(ω) indicates whether a system with semi-infinite boundary conditions will host a localized mode at the frequency ω. When ν(ω) > 0, a semi-infinite system with domain x [0,∞) will host a mode localized to its left boundary at frequency ω. Likewise, when ν(ω) < 0, a semi-infinite system with domain x (−∞,0], will host a mode localized to its right boundary.

As shown in Fig. 3c, d, the sign of κ can be rationalized by examining ν(ω) for an example frequency ω (denoted by the star) along the real axis. For P > 0, the periodic boundary spectrum (red) winds once counter-clockwise around the star, and hence ν(ω) < −1. However, for P < 0, the localization is reversed since the direction of the contour is reversed. In the Methods, we show how to compute ν(ω) directly from the continuum equations using a generalization of Eq. (30). Furthermore, in Supplementary Note 1 we discuss how the presence of additional vibrational bands affect the calculation and physical interpretation of ν(ω).

Pseudo-Hermitian dynamics and direction-dependent bending modulus

In the preceding section, we took P(ω) to be real and constant. However, by introducing a phase delay into the transfer function H(ω) we can control the complex argument of P. When arg(P) = ±π/2, the secular Eq. (27) has entirely real coefficients. Hence, the periodic boundary spectrum will consist of frequencies that come in real values or complex conjugate pairs. This additional symmetry is sometimes referred to as a generalized PT symmetry30,63,64, which arises if and only if there exists an antiunitary operator that commutes with D(k). We say that the PT symmetry is unbroken when the eigenvalues of D(k) are entirely real, and that it is broken otherwise. In the unbroken phase, D(k) is said to be pseudo-Hermitian. Pseudo-Hermiticity implies that each eigenvector of D(k) will individually conserve the mechanical energy density e in Eq. (26). However, when two or more eigenvectors are superimposed, e can oscillate in time, though remaining centered around a constant time-averaged value.

Since pseudo-Hermiticity constrains the periodic boundary spectrum of D(k) to lie along the real line, the unidirectional amplification vanishes when P is imaginary. Nonetheless, the effects of parity violation are still present. We note that Eq. (28) may be written in the form:

$${\omega }_{\pm }=\pm \sqrt{\frac{B}{\rho }}\left(1\pm i\frac{{Pk}}{\mu }\right){k}^{2}+{{{{{\mathscr{O}}}}}}\left({k}^{4}\frac{{I}^{2}{B}^{2}}{{\rho }^{2}{\mu }^{2}}\right)$$

When arg(P) = ±π/2, we can interpret the form of Eq. (28) as having a rescaled bending modulus:

$${B}_{{{{{{\rm{eff}}}}}}}=B{\left(1\mp \frac{\left|P\right|k}{2\mu }\right)}^{2}$$

Notice that the value of Beff depends on the sign of k, and hence the effective bending modulus is direction dependent. Fig. 4a shows the tilting of the dispersion for arg(P) = ±π/2. The tilt implies that the phase and group velocities for right and left traveling waves are unequal. Numerically solved modes are shown in Fig. 4b. We note that the pseudo-Hermiticity endowed by P(ω)  ± i must exist exclusively at finite frequency because P(ω) cannot be nonzero and imaginary at ω = 0 due to the requirement that P(−ω) = P*(ω).

Fig. 4: Pseudo-Hermitian dynamics.

a The spectrum is shown for the metamaterial with arg(P) = π/2. We note that the reality of the frequencies is maintained, while the modulus P breaks the k → −k symmetry. L is the unit cell length. b Transverse displacement wave fields for the waves traveling in different directions. The left and right traveling modes are excited at equal frequencies, but have differing wavenumbers due to the odd micropolarity. The red arrows indicate the direction of travel of the wave, and H is the transfer function such that P = ΠH.

A discrete model of the odd micropolar metabeam

To gain intuition into the mechanics of the metabeam, it is useful to consider a discrete model. As shown in Fig. 5a, b, the ith unit cell of the discrete model consists of a rod (gray) with moment of inertia J and total mass m, whose position and orientation are captured by a height hi and an angle φi−1. The rods are connected via mass-less, rigid frames and two springs. The top spring is a Hookean spring with tension

$${T}_{i}={k}_{\mu }({h}_{i-1}-{h}_{i}+L{\varphi }_{i-1})$$

and a bottom spring is a torsional spring that exerts an angular tension

$${\tau }_{i}={\kappa }_{B}({\varphi }_{i-1}-{\varphi }_{i})$$
Fig. 5: Discrete model for odd micropolar beam.

a A discrete model of a Timeoshenko beam consists of a central mass m with moment of inertial J, a Hookean spring of spring constant kμ and a torsional spring of spring constant κ, Band lattice spacing L. b The unit cell is described by the height of the mass hi and the angle φi of the black connecting rod. c, d The odd micropolar beam has an internal feedback P that senses the angle change of the torsional spring and actuates additional tension in the Hookean spring. The control loop is unidirectional: stretching or compressing the Hookean spring does not affect the torsional spring.

In Fig. 5c, d, we show the addition of an active element that senses the stretching of the bottom spring and actuates an additional tension in the top spring

$${T}_{i}^{a}=p({\varphi }_{i-1}-{\varphi }_{i})$$

Summing the forces in the vertical direction yields a dynamical equation for hi

$$ m{\ddot{h}}_{i}={T}_{i}-{T}_{i+1}+{T}_{i}^{a}-{T}_{i+1}^{a}$$
$$ = {k}_{\mu }\left({h}_{i+1}+{h}_{i-1}-2{h}_{i}\right)+L{k}_{\mu }\left({\varphi }_{i-1}-{\varphi }_{i}\right)\\ \quad +p({\varphi }_{i+1}+{\varphi }_{i-1}-2{\varphi }_{i})$$

Furthermore, summing the torques yields

$$ J{\ddot{\varphi }}_{i}={\tau }_{i}-{\tau }_{i+1}-L\left({T}_{i+1}+{T}_{i+1}^{a}\right)$$
$$ = {\kappa }_{B}\left({\varphi }_{i+1}+{\varphi }_{i-1}-2{\varphi }_{i}\right)+{pL}\left({\varphi }_{i+1}-{\varphi }_{i}\right)\\ \quad+L{k}_{\mu }({h}_{i+1}-{h}_{i}-L{\varphi }_{i})$$

Upon inspection, Eqs. (37) and (39) are precise discretizations of Eqs. (21) and (22) with ρ = m/L, I = J/L, μ = kμL, B = κBL, and P = Lp. See Fig. S4 for a comparison of the dispersion for the discrete model and continuum theory. Notice that the discrete model manifestly conserves linear and angular momentum since the linear and torsional springs exert equal and opposite forces and torques, respectively, on the units they connect. Even without externally applied torques, nontrivial internal angular momentum transfer occurs between the translation of its center of mass jLm\(\dot{{h}_{j}}\) and the rotation of the axis J\(\dot{{\varphi }_{j}}\), akin to a “spin-orbit” coupling. Nonetheless, Maxwell-Betti reciprocity is violated by the asymmetry in the relationship between the linear and torsional springs: bending of the torsional spring induces a tension Tai = p(φi−1φi) in the linear spring, while the deformation of the linear spring hi−1hi+i−1 has no response in the angular spring. This asymmetry implies that a cycle of alternating actuation and release of the linear and torsional spring is associated with a nonzero amount of work done. For additional discrete models illustrating the independence of Maxwell-Betti reciprocity and momentum conservation, see Supplementary Note 1.

Experimental demonstration

To probe the dynamic wave phenomena originating from the odd micropolar modulus P, we perform experiments in which we excite flexural waves in the metabeam using piezoelectric actuators, see Fig. 6a and “Methods”. In experiments, we implement the following electronic transfer function (cf. Fig. 1c):

$$H\left(\omega \right)=\frac{{H}_{0}}{{\left(i\omega /{\omega }_{0}\right)}^{2}+2i\zeta \omega /{\omega }_{0}+1}$$
Fig. 6: Experimental demonstration of skin modes and odd micropolar moduli.

a Experimental schematic. Flexural waves are generated in the active metabeam from either the right or left side using piezoelectric actuators (yellow), see “Methods”. A scanning laser Doppler vibrometer (SLDV) measures the transverse velocity of the surface of the active metabeam. b, c Unidirectional amplification of waves. A metamaterial consisting of 9 unit cells is actuated from either the right (blue) or left (red) with a 2 kHz tone burst signal (gray). The output velocity is normalized by the maximum velocity observed when the experiment is performed with no active feedback. d Observation of the non-Hermitian skin effect. Experiments are performed between 1.5 kHz and 4 kHz for right to left (blue) and left to right (red) traveling waves. A 2D FFT shows the intensity of the observed spectrum. The intensity is normalized by its maximum value. e The inverse decay length. In d, e the solid theoretical curves are based on the transfer matrix method. In d the gray dashed curves are theoretical predictions with no activity. f A plot of arg(P) as a function of frequency. At ω = ω0 (= 3 kHz), arg(P) = −π/2, indicating that the system is pseudo-Hermitian and accordingly we observe κ = 0 at ω = ω0.

Here, ζ = 0.48 and H0 = 3 are constants and ω0 = 3 kHz is the cutoff frequency at which arg(P) = −π/2. Other electrical circuit and geometric parameters of the metabeam can be found in Supplementary Note 2. We probe the vibrational dynamics of the beam by initiating waves from either the right or left side of the metamaterial via external piezoelectric elements. Fig. 6b, c shows the experiment at 2 kHz in which waves from the right are suppressed while waves from the left are amplified (see also Supplementary Movies 1 and 2).

To construct the full spectrum of the metamaterial, we perform the experiment with tone burst signals centered between 1.5 and 4 kHz. The transverse velocity wave fields are measured along the medium using a laser Doppler vibrometer. We apply fast-Fourier transforms in time and in space to extract the real k(ω) and imaginary part κ(ω) of the wavenumber as a function of frequency for right-going (red) and left-going (blue) waves (Fig. 6d, e). As described in Supplementary Note 3, the solid theoretical curves are produced using a semi-analytical technique known as the transfer-matrix method. The transfer-matrix method utilizes the beam geometry, known material parameters, and electronic feedback measured in simulations. No fitting parameters are used in the comparison between experiment and the transfer-matrix method curve.

As illustrated in Fig. 6f, the transfer function H(ω) is chosen such that −π/2 < arg(P) < 0 when ω < ω0. In this case, we find κ < 0 for both left- and right-propagating waves (panel e). A value of κ < 0 implies that waves propagating to the left are attenuated whereas waves propagating to the right are amplified, which is confirmed by the FFT intensity in panel d. Likewise −π < arg(P) < −π/2 for ω > ω0. In this case, we find κ > 0, indicating that waves propagating to the right are attenuated whereas waves propagating to the left are amplified. In addition, when ω = ω0, arg(P) = −π/2, and the waves propagating to the left and right display no attenuation or amplification. At this frequency, the effective dynamical matrix is pseudo-Hermitian, and the differences between left- and right-propagating waves reside only in the wavelengths and phase velocities. We note that, in experiments, the unidirectional amplification is limited by the maximum output voltage (±45 V) of our electric control system. In practice, to maximize the amplification ratio, one strategy is to use a modest value of the transfer function, say |H0| ≈ 3, and increase the number of unit cells over which the wave is amplified.


The metabeam presented here demonstrates active odd micropolar moduli and non-reciprocal responses absent in energy conserving media that are enabled by sensing, actuating and local computation. The minimal on-board electronics that power the active metabeam enable its multiple functions as an elastic engine, selective mode amplifier, and mechanical diode. We uncover an intrinsic relation between an odd micropolar modulus, the non-Hermitian skin effect, and a corresponding topological index. Numerical and experimental results show unidirectional amplification and attenuation of waves propagating through the metamaterial. Odd micropolarity extends the range of possible couplings between conventional strains/stresses and higher-order curvatures/moments by including antisymmetry in their relationship. The electronically assisted mechanical feedback provides an appealing solution to precisely modulate odd micropolar moduli without requiring changes to the metabeam’s structure, geometry, or passive moduli. Our design can be flexibly tuned through computer coding and scaled via microelectromechancial systems (MEMS)65,66. Our mechanical approach relies on a feed-forward control loop, a generic concept that can exist in both metamaterial and biological contexts. The continuum theory also makes our approach especially generalizable to the mechanics in other systems. Combining the principles illustrated here with disorder, nonlinearities, and strong dissipation suggests new approaches for the control of filaments and membranes arising in biological media4,18,19.


Sample fabrication

The metabeam is composed of three piezoelectric patches (STEMiNC PZT 5J: 6 mm × 4 mm × 0.55 mm) mounted via conductive epoxy onto a laser-cut stainless steel host beam. We achieve antisymmetric actuation without the use of an inverting voltage amplifier by mounting the two piezoelectric actuators such that their piezoelectric polarization directions are oppositely oriented.

Experimental procedures

In experiments, nine metamaterial unit cells are connected with control circuits, see Fig. 6a. Two piezoelectric transducers are attached on the left and right sides of the metamaterial to generate incident flexural waves. We employ ten-peak tone-burst signals with central frequencies ranging from 1.5 to 4.0 kHz in step sizes of 0.1 kHz. We generate and amplify incident wave signals via an arbitrary waveform generator (Tektronix AFG3022C) and a high voltage amplifier (Krohn-Hite), respectively. Transverse velocity wavefields are measured on the surface of the metamaterial by a scanning laser Doppler vibrometer (Polytec PSV-400). We note that the transfer-matrix method used to derive the theoretical curves in Fig. 6d, e rigorously assumes an infinite system. To experimentally approximate these conditions, we embed the active metamaterial within a larger host steel beam denoted by the gray region of Fig. 6a. When waves cross the boundaries from host beam to the metamaterial, the reflection at boundaries between the host beam and the metamaterial is be negligible, as evidenced by our numerical and experimental results. To suppress reflected waves at the free boundaries of the host beam, we bonded two layers of clay on the host beam with sufficient lengths. This way, waves can propagate through the metamaterial with approximated infinite boundary conditions. The decay length is then obtained by calculating the wave amplitudes at different points in the metamaterial.

Finite element simulations

We calibrate the transfer matrix method and continuum equations by conducting fully three-dimensional numerical simulations of the unit cell using the commercial finite element software COMSOL Multiphysics. In all the simulations, we model the piezoelectric patches via a three-dimensional linear piezoelectric constitutive law. The central piezoelectric patch acts as a sensor whose signal is obtained by integrating the free charge over the top surface of the piezoelectric sensor. The top and bottom surfaces of the piezoelectric sensor have zero electric potential. The bottom surfaces on the piezoelectric actuators are ground, and we apply electrical potentials on the top surfaces to act as actuating voltages. The actuating voltages are related to the sensing voltages via the electronic transfer function. For the wave dispersion computations in Fig. 3a, Floquet periodic boundary conditions are applied on the left and right boundaries of a metamaterial unit cell. We calculate eigenfrequencies of the unit cell with different real wavenumbers. To simulate the wave propagation with open boundaries (Fig. 3b), a metabeam composed of 15 unit cells is placed between two external beams. Two perfectly matched layers (PMLs) are attached to both ends of the external beams in order to suppress reflected waves from the boundaries. The incident flexural wave is generated by applying a harmonic force on the boundary of the host beam. The out-of-plane displacement is measured at the left- and right-hand sides of the metabeam. The penetration depth is calculated by comparing the amplitudes of the two extracted displacements.

Energy bounds on dynamic moduli

Here we discuss Eqs. (12, 13) in the main text. For simplicity, let us collect the stresses into a vector t = (σzx, M)T and the deformations into a vector u = (s, b)T. Suppose the beam is subject to a deformation procedure such that its initial and final configurations at times t = −∞ and t = ∞ are identical. Then the total work per unit volume done by the beam is given by

$$\Delta W=\int _{-{{{{{\rm{\infty }}}}}}}^{{{{{{\rm{\infty }}}}}}}\frac{{{{{{\rm{d}}}}}}{{{{{\bf{u}}}}}}}{{{{{{\rm{d}}}}}}t}\cdot {{{{{\bf{t}}}}}}{{{{{\rm{d}}}}}}t$$
$$=-i\int_{-{{{{{\rm{\infty }}}}}}}^{{{{{{\rm{\infty }}}}}}}\omega {{{{{{\bf{u}}}}}}}^{{{\dagger}} }\left(\omega \right)\cdot {{{{{\bf{C}}}}}}\left(\omega \right)\cdot {{{{{\bf{u}}}}}}\left(\omega \right){{{{{\rm{d}}}}}}\omega$$
$$=\int_{-{{{{{\rm{\infty }}}}}}}^{{{{{{\rm{\infty }}}}}}}\omega {{{{{{\bf{u}}}}}}}^{{{\dagger}} }\left(\omega \right)\cdot {{{{{\bf{M}}}}}}\left(\omega \right)\cdot {{{{{\bf{u}}}}}}\left(\omega \right){{{{{\rm{d}}}}}}\omega$$

In the final step, we have introduced the matrix M(ω) = i[C(ω) − C(ω)] and used the fact that u(−ω) = u*(ω) and C(−ω) = C*(ω). If the medium is passive, then we require that ∆W must be negative for all choices of u(ω). Therefore, we require that the matrix M(ω) be negative semidefinite. Using the parameterization in Eq. (11), this requirement implies Eqs. (12, 13). See refs. 67,68,69 for related discussions in two- and three-dimensional media.

Cycles at finite frequency

To gain intuition on the elastodynamics of the odd micropolar metabeam, it is useful to consider the notion of a cycle at finite frequency. At a finite frequency ω, the modulus P need not be real and we may write P = |P| eP. In this case, both the real and imaginary parts of P will contribute to the energy extracted. For example, consider a cyclic protocol that involves bending of amplitude |b| and shearing of amplitude |s| and relative phase delay φB. Applying Eq. (43), the total energy extracted per cycle is

$${{{{{\rm{Work}}}}}}=-\pi \left|P\right|\left|s\right|\left|b\right|{{{{{\rm{sin }}}}}}\left({\phi }_{B}+{\phi }_{P}\right)$$

Notice that Eq. (44) gives mechanistic insight into the amplification of the waves observed in the experiment. When P is nonzero, the eigenmodes of D(k) comprising a given plane wave will generically have a phase delay between bending and shearing. Hence, after one cycle, the nonconservative stresses will have converted stored electrical energy into mechanical energy. This conversion will cause the amplitude of the eigenvector to grow in proportion to its current amplitude, for which |s| |b| is a proxy. Hence, the mode will be exponentially amplified, as reflected by the imaginary component of the eigenfrequencies.

Topological index in the continuum

In this section, we discuss the index ν(ω) from the point of view of the continuum theory54. The spectrum as a function of k is plotted in Fig. 7a for P = 0. The spectrum is given by the roots of the secular Eq. (27) and contains four roots since the equations of motion are second order in time and involve two coupled fields. We explicitly compute the eigenvectors and eigenvalues for small k and P = 0 in Supplementary Note 1. The spectrum consists of a pair of Goldstone modes, which for small k and P = 0 represents a flexural motion of the beam. Additionally, the continuum equations imply two modes separated by a band gap ω1 ≈ 105 Hz. As shown in Supplementary Note 1, for small k and P = 0, these modes are dominated by a shearing motion. For the Goldstone mode, we expect a range of validity of the continuum theory at small k, since ω(k → 0) = 0. However, for the shear dominated mode, the continuum theory is not expected to self-consistently apply due to the finite gap. In Supplementary Note 1, we numerically compute the spectrum and eigenmodes for frequencies above the experimentally relevant range using COMSOL. There we also discuss why the winding number ν(ω) computed in the continuum is still physically relevant despite the presence of additional high frequency modes not captured by the continuum.

Fig. 7: Non-Hermitian band topology via odd micropolar elasticity.

a The spectrum for P = 0 features a pair of bending dominated bands (black) and shear dominated bands (gray) separated by a band gap ω1. (b-d) The spectrum is shown in the complex ω plane for P = 0, P > 0, and P < 0. The thick black lines represent the bending dominated band, while the thick gray lines represent the shear dominated bands, both with k [−R, R] for a finite R. The thin black lines represent the analytical continuation of the spectrum for k = Re for φ [0, π]. The arrows indicate the direction of increasing k. The numbers indicate the value of ν(ω) for ω in the corresponding colored regions of the complex plane. This number corresponds to the number of times that the spectrum winds around a given region. e For a semi-infinite system with a free boundary, the winding number of \(\widetilde{\nu }\) = 1 (\(\widetilde{\nu }\)  = 3) for our continuum theory indicates a mode localized to the right (left) boundary. The wave forms schematically depict the localization with A(x) representing amplitude. For a calculation of the precise eigenvectors, see Supplementary Note 1.

For a translationally invariant system, the difference between periodic and open boundary conditions is whether the differential operator that defines the equations of motion allows formal eigenvectors with complex wavenumber k. For a system with periodic boundary conditions, the spectrum consists of values of ω that solve det[D(k) − ω] = 0 for real k. However, for a system on a semi-infinite domain x [0,∞), we allow Im(k) > 0. These modes decay to the right as exp[(−Imk+iRek)x] and therefore maintain a finite L2 norm, which represents the mechanical energy density. To determine whether a given frequency ω is in the semi-infinite boundary spectrum, we first must count the number of solutions to det[D(k) − ω] = 0 with Imk > 0. We do so by applying Cauchy’s argument principle from complex analysis:

$$\widetilde{\nu }(\omega )=\mathop{\lim}\limits_{R\to {{{{{\rm{\infty }}}}}}}\frac{1}{2\pi i}{\oint \!\!}_{\Gamma \left(R\right)}\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}k}{{\log }}\,{{\det }}[{{{{{\boldsymbol{D}}}}}}\left(k\right)-\omega ]{{{{{\rm{d}}}}}}k$$

where Γ(R) is a counter clockwise curve in the complex plane given by [−R, R] together with Re for φ [0, π]. The winding number \(\widetilde{\nu }(\omega )\) itself does not directly determine the number of localized modes, since we must also consider the boundary conditions placed at x = 0. Suppose that γ independent homogeneous boundary conditions are placed at x = 0. Then the number of left localized modes will generically be given by:

$${{{{{\rm{left}}}}}}\,{{{{{\rm{localized}}}}}}\,{{{{{\rm{modes}}}}}}\,{{{{{\rm{at}}}}}}\,\omega =\tilde{v}(\omega )-\gamma$$

Likewise, for a system with boundary (0, ∞], the number of right localized modes is given by

$${{{{{\rm{right}}}}}}\,{{{{{\rm{localized}}}}}}\,{{{{{\rm{modes}}}}}}\,{{{{{\rm{at}}}}}}\,\omega =d-\tilde{v}(\omega )-\gamma$$

where d is the degree of det[D(k) − ω] as a polynomial in k. Whenever the right-hand side of Eqs. (46) or (47) is negative, the mode count is taken to be zero. In Supplementary Note 1, we provide a derivation of Eqs. (46) and (47) and explicitly define independent, homogeneous boundary conditions. Such boundary conditions include, for example, stress-free (M = σzx = 0) or motion-free (h = φ = 0) boundaries. Fig. 7b–d illustrates the computation of the winding number. The thick black and gray lines are the periodic boundary spectrum for k [−R, R], and the thin black lines are the analytical continuation for k = Re for φ [−π, π]. Colored regions are labeled by the value of \(\widetilde{\nu }(\omega )\). Suppose for example that the beam is given stress and moment free boundary conditions σzx = 0 and M = 0 at x = 0. In this case γ = 2, and therefore \(\widetilde{\nu }(\omega )\) = 1 indicates the presence of a right localized mode, \(\widetilde{\nu }(\omega )\) = 3 indicates the presence of a left localized modes, as shown in Fig. 7e. In the Supplementary Note 1, we explicitly compute examples of the eigenmodes in Fig. 7c.

Topological index from discrete models

We now discuss the topological index in Eq. (30), which is appropriate for discrete settings such as the discrete model and finite element simulations. Suppose the system is composed out of a unit cell of finite length L whose internal state is represented by a vector Ψ(x). The components of Ψ can represent, for example, the displacement and velocities of points in a finite element mesh. Here, x is a discrete label that takes values in integer multiples of L. The Fourier transform of the equations of motion now read:

$$\omega {{{{{\boldsymbol{\Psi }}}}}}\left(k\right)={{{{\boldsymbol{{{{{\mathscr{D}}}}}}}}}}\left(k\right)\cdot {{{{{\boldsymbol{\Psi }}}}}}(k)$$

where D(k) is the dynamical matrix for the discrete system, and k is the wave number assuming values in [−π/L, π/L]. For a system with periodic boundaries, the spectrum is given by solutions to det[D(k) − ω] = 0 for k [−π/L, π/L]. For a system with semi-infinite boundaries (e.g., x {0, L, 2 L,…}), we allow Imk > 0. As detailed in the Supplementary Note 1, we can invoke a similar application of the Cauchy argument principle to determine the number of eigenmodes at a given frequency ω. To do so, we write k = −ilogλ, where λ assumes values on the unit circle S in the complex plane. We can then apply Cauchy’s argument principle using S as a counterclockwise contour for λ. We have:

$$\nu \left(\omega \right)=\frac{1}{2\pi i}\int_{S}\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}\lambda }{{\log }}\,{{\det }}\left[{{{{\boldsymbol{{{{{\mathscr{D}}}}}}}}}}\left({{{{{\boldsymbol{-}}}}}}i{{\log }}\lambda \right){{{{{\boldsymbol{-}}}}}}\omega \right]{{{{{\rm{d}}}}}}\lambda$$
$$=\frac{1}{2\pi i}\int_{-\frac{\pi }{L}}^{\frac{\pi }{L}}\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}k}{{\log }}\,{{\det }}\left[{{{{\boldsymbol{{{{{\mathscr{D}}}}}}}}}}\left({{{{{\bf{k}}}}}}\right){{{{{\boldsymbol{-}}}}}}\omega \right]{{{{{\rm{d}}}}}}k$$
$$=\frac{1}{2\pi i}\mathop{\sum}\limits_{\alpha }\int_{-\pi /L}^{\pi /L}\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}k}{{\log }}[{\omega }_{\alpha }\left(k\right)-\omega ]{{{{{\rm{d}}}}}}k$$

where ωα(k) is the value of the periodic spectrum of the band α, in agreement with Eq. (30) from the main text. Notice that for the frequency denoted by the red star in Fig. 3c, we have ν(ω) = −1. However, for the same point in Fig. 7c (which lies in the red region to the right of the origin), we have \(\widetilde{\nu }(\omega )\) = 1. To see the relationship between these quantities, notice that \(\widetilde{\nu }(\omega )\) in Eq. (45) counts the number of zeros of f(k) = det[D(k) − ω] in the upper half plane. However, in computing ν from Eq. (49), the interior of S contains not only the zeros of F(λ) ≡ det[D(−ilogλ) − ω] but also a set of poles. Each of these poles physically represents a boundary condition that arises when transitioning between a Laurent operator to a Toeplitz operator (see Supplementary Note 1 for details). The number of poles depends on the precise discretization of the continuum equations. Hence the winding number ν(ω) as given by Eq. (30) represents the difference between the number of zeros (candidate modes) and the number of poles (boundary conditions). Therefore, one should compare ν(ω) to \(\widetilde{\nu }(\omega )\) − γ for an appropriate choice of γ. In practice, γ is determined by physically interpreting the boundary conditions implied by the discretization. In Supplementary Note 1, we show that the discretization used in the discrete model and the COMSOL simulations enforce displacement-free boundary conditions and hence corresponds to γ = 2. Hence \(\widetilde{\nu }(\omega )\) – γ = ν(ω) = −1 as required for physical consistency.

Calculation of penetration depth

In Fig. 3b, we show the penetration depth for modes associated with positive real ω. We can compute this expression by solving ω+(k) = ω in Eq. (28) for complex k. In particular, substituting k = q +  into Eq. (28) yields:

$$0=2k\kappa +\frac{P}{2\mu }({k}^{3}-3k{\kappa }^{2})$$
$$\omega =\sqrt{\frac{B}{\rho }}\left[\left(1-\frac{P\kappa }{\mu }\right)\left({k}^{2}-{\kappa }^{2}\right)-\frac{2P}{\mu }{k}^{2}\kappa \right]$$

Solving Eqs. (52, 53) to leading order in k, we obtain Eq. (29).

Data availability

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Information. Extended data and materials in the main text and the Supplementary Information are available upon request by contacting the corresponding authors.

Code availability

The computer code and algorithm that support the findings of this study are available from the corresponding author upon reasonable request.


  1. 1.

    Bertoldi, K., Vitelli, V., Christensen, J. & van Hecke, M. Flexible mechanical metamaterials. Nature Reviews. Materials 2, 17066 (2017).

    CAS  Google Scholar 

  2. 2.

    Huber, S. D. Topological mechanics. Nat. Phys. 12, 621–623 (2016).

    CAS  Article  Google Scholar 

  3. 3.

    Frenzel, T., Kadic, M. & Wegener, M. Three-dimensional mechanical metamaterials with a twist. Science 358, 1072–1074 (2017).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  4. 4.

    Needleman, D. & Dogic, Z. Active matter at the inter-face between materials science and cell biology. Nature Reviews. Materials 2, 17048 (2017).

    CAS  Google Scholar 

  5. 5.

    Coulais, C., Teomy, E., de Reus, K., Shokef, Y. & vanHecke, M. Combinatorial design of textured mechanical metamaterials. Nature 535, 529–532 (2016).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    Guseinov, R., McMahan, C., Perez, J., Daraio, C. & Bickel, B. Programming temporal morphing of self-actuated shells. Nat. Commun. 11, 237 (2020).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. 7.

    Brandenbourger, M., Locsin, X., Lerner, E. & Coulais, C. Non-reciprocal robotic metamaterials. Nat. Commun. 10, 4608 (2019).

    ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

  8. 8.

    Caruel, M. & Truskinovsky, L. Physics of muscle contraction. Rep. Prog. Phys. 81, 036602 (2018).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  9. 9.

    Wang, R. et al. Torsional refrigeration by twisted, coiled, and supercoiled fibers. Science 366, 216–221 (2019).

    ADS  CAS  PubMed  Article  Google Scholar 

  10. 10.

    Chen, Y. Y., Zhu, R., Barnhart, M. V. & Huang, G. L. Enhanced flexural wave sensing by adaptive gradient-index metamaterials. Sci. Rep. 6, 35048 (2016).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. 11.

    Fleury, R., Sounas, D. & Al’u, A. An invisible acoustic sensor based on parity-time symmetry. Nat. Commun. 6, 5905 (2015).

    ADS  PubMed  Article  Google Scholar 

  12. 12.

    Peng, B. et al. Loss-induced suppression and revival of lasing. Science 346, 328–332 (2014).

    ADS  CAS  PubMed  Article  Google Scholar 

  13. 13.

    Li, S. et al. Particle robotics based on statistical mechanics of loosely coupled components. Nature 567, 361–365 (2019).

    ADS  CAS  PubMed  Article  Google Scholar 

  14. 14.

    Miskin, M. Z. et al. Electronically integrated, mass-manufactured, microscopic robots. Nature 584, 557–561 (2020).

    ADS  CAS  PubMed  Article  Google Scholar 

  15. 15.

    Lin, Z. et al. Unidirectional invisibility induced by PT-symmetric periodic structures. Phys. Rev. Lett. 106, 213901 (2011).

    ADS  PubMed  Article  CAS  Google Scholar 

  16. 16.

    Yi, K., Collet, M., Chesne, S. & Monteil, M. Enhancement of elastic wave energy harvesting using adaptive piezolens. Mech. Syst. Signal Process. 93, 255–266 (2017).

    ADS  Article  Google Scholar 

  17. 17.

    Mohammadi Estakhri, N., Edwards, B. & Engheta, N. Inverse-designed metastructures that solve equations. Science 363, 1333–1338 (2019).

    ADS  MathSciNet  CAS  PubMed  PubMed Central  MATH  Article  Google Scholar 

  18. 18.

    Salbreux, G. & J ̈ulicher, F. Mechanics of active surfaces. Phys. Rev. E 96, 032404 (2017).

    ADS  PubMed  Article  Google Scholar 

  19. 19.

    Prost, J., Julicher, F. & Joanny, J.-F. Active gel physics. Nat. Phys. 11, 111–117 (2015).

    CAS  Article  Google Scholar 

  20. 20.

    Nassar, al. Nonreciprocity in acoustic and elastic materials. Nat. Rev. Mater. 5, 667–685 (2020).

  21. 21.

    Coulais, C., Sounas, D. & Al’u, A. Static non-reciprocity in mechanical metamaterials. Nature 542, 461–464 (2017).

    ADS  CAS  PubMed  Article  Google Scholar 

  22. 22.

    Fleury, R., Sounas, D. L., Sieck, C. F., Haberman, M. R. & Al’u, A. Sound isolation and giant linear nonreciprocity in a compact acoustic circulator. Science 343, 516–519 (2014).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. 23.

    Fruchart, M., Hanai, R., Littlewood, P. B. & Vitelli, V. Non-reciprocal phase transitions. Nature 592, 363–369 (2021).

    ADS  CAS  PubMed  Article  Google Scholar 

  24. 24.

    Saha, S., Agudo-Canalejo, J. & Golestanian, R. Scalar active mixtures: The nonreciprocal Cahn-Hilliard model. Phys. Rev. X, 10, 041009 (2020).

    CAS  Google Scholar 

  25. 25.

    You, Z., Baskaran, A. & Marchetti, M. C. Nonreciprocity as a generic route to traveling states. Proc. Natl Acad. Sci. USA 117, 19767–19772 (2020).

    ADS  MathSciNet  CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Gupta, R. K., Kant, R., Soni, H., Sood, A. K. & Ramaswamy, S. Active nonreciprocal attraction between motile particles in an elastic medium. arXiv:2007.04860. Preprint at (2020).

  27. 27.

    Scheibner, C. et al. Odd elasticity. Nat. Phys. 16, 475–480 (2020).

    CAS  Article  Google Scholar 

  28. 28.

    Faust, D. & Lakes, R. S. Reciprocity failure in piezo-electric polymer composite. Phys. Scr. 90, 085807 (2015).

    ADS  Article  CAS  Google Scholar 

  29. 29.

    Rosa, M. I. N. & Ruzzene, M. Dynamics and topology of non-Hermitian elastic lattices with non-local feedback control interactions. N. J. Phys. 22, 053004 (2020).

    MathSciNet  Article  Google Scholar 

  30. 30.

    Scheibner, C., Irvine, W. T. M. & Vitelli, V. Non-Hermitian band topology and skin modes in active elastic media. Phys. Rev. Lett. 125, 118001 (2020).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  31. 31.

    Zhou, D. & Zhang, J. Non-Hermitian topological meta-materials with odd elasticity. Phys. Rev. Res. 2, 023173 (2020).

    CAS  Article  Google Scholar 

  32. 32.

    Banerjee, D., Vitelli, V., Jülicher, F. & Surówka, P. Active viscoelasticity of odd materials. Phys. Rev. Lett. 126, 138001 (2021).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  33. 33.

    Trainiti, G. et al. Time-periodic stiffness modulation inelastic metamaterials for selective wave filtering: Theory and experiment. Phys. Rev. Lett. 122, 124301 (2019).

    ADS  CAS  PubMed  Article  Google Scholar 

  34. 34.

    Fink, M. et al. Time-reversed acoustics. Rep. Prog. Phys. 63, 1933–1995 (2000).

    ADS  Article  Google Scholar 

  35. 35.

    Sirota, L., Sabsovich, D., Lahini, Y., Ilan, R. & Shokef, Y. Real-time steering of curved sound beams in a feedback-based topological acoustic metamaterial. Mech. Syst. Signal Process. 153, 107479 (2021).

    Article  Google Scholar 

  36. 36.

    Sirota, L., Ilan, R., Shokef, Y. & Lahini, Y. Non-newtonian topological mechanical metamaterials using feedback control. Phys. Rev. Lett. 125, 256802 (2020).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  37. 37.

    Ouisse, M., Collet, M. & Scarpa, F. A piezo-shunted kirigami auxetic lattice for adaptive elastic wave filtering. Smart Mater. Struct. 25, 115016 (2016).

    ADS  Article  Google Scholar 

  38. 38.

    Yi, K. et al. Programmable metamaterials with digital synthetic impedance circuits for vibration control. Smart Mater. Struct. 29, 035005 (2020).

    ADS  CAS  Article  Google Scholar 

  39. 39.

    Chen, Y. Y., Huang, G. L. & Sun, C. T. Band gap control in an active elastic metamaterial with negative capacitance piezoelectric shunting. J. Vib. Acoust. 136, 061008 (2014).

    Article  Google Scholar 

  40. 40.

    Bergamini, A. et al. Phononic crystal with adaptive connectivity. Adv. Mater. 26, 1343–1347 (2014).

    CAS  PubMed  Article  Google Scholar 

  41. 41.

    Alan, S., Allam, A. & Erturk, A. Programmable mode conversion and bandgap formation for surface acoustic waves using piezoelectric metamaterials. Appl. Phys. Lett. 115, 093502 (2019).

    ADS  Article  CAS  Google Scholar 

  42. 42.

    Chen, Y., Hu, G. & Huang, G. A hybrid elastic metamaterial with negative mass density and tunable bending stiffness. J. Mech. Phys. Solids 105, 179–198 (2017).

    ADS  MathSciNet  Article  Google Scholar 

  43. 43.

    Wang, G., Cheng, J., Chen, J. & He, Y. Multi-resonant piezoelectric shunting induced by digital controllers for subwavelength elastic wave attenuation in smart metamaterial. Smart Mater. Struct. 26, 025031 (2017).

    ADS  Article  Google Scholar 

  44. 44.

    Casadei, F., Delpero, T., Bergamini, A., Ermanni, P. & Ruzzene, M. Piezoelectric resonator arrays for tunable acoustic waveguides and metamaterials. J. Appl. Phys. 112, 064902 (2012).

    ADS  Article  CAS  Google Scholar 

  45. 45.

    Chen, Y. Y., Hu, G. K. & Huang, G. L. An adaptive metamaterial beam with hybrid shunting circuits for extremely broadband control of flexural waves. Smart Mater. Struct. 25, 105036 (2016).

    ADS  Article  Google Scholar 

  46. 46.

    Sugino, C., Ruzzene, M. & Erturk, A. Merging mechanical and electromechanical bandgaps in locally resonant metamaterials and metastructures. J. Mech. Phys. Solids 116, 323–333 (2018).

    ADS  MathSciNet  Article  Google Scholar 

  47. 47.

    Eringen, A. C. Microcontinuum Field Theories (Springer, 1999).

  48. 48.

    Maugin, G. A. On the structure of the theory of polar elasticity. Philos. Trans.: Math., Phys. Eng. Sci. 356, 1367–1395 (1998).

    ADS  MathSciNet  CAS  MATH  Article  Google Scholar 

  49. 49.

    Chen, Y., Li, X., Nassar, H., Hu, G. & Huang, G. A programmable metasurface for real time control of broad-band elastic rays. Smart Mater. Struct. 27, 115011 (2018).

    ADS  Article  Google Scholar 

  50. 50.

    Chen, Y., Li, X., Hu, G., Haberman, M. R. & Huang, G. An active mechanical Willis meta-layer with asymmetric polarizabilities. Nat. Commun. 11, 3681 (2020).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. 51.

    Hatano, N. & Nelson, D. R. Localization transitions in non-Hermitian quantum mechanics. Phys. Rev. Lett. 77, 570–573 (1996).

    ADS  CAS  PubMed  Article  Google Scholar 

  52. 52.

    Bergholtz, E. J., Budich, J. C. & Kunst, F. K. Exceptional topology of non-Hermitian systems. Rev. Mod. Phys. 93, 015005 (2021).

    ADS  MathSciNet  Article  Google Scholar 

  53. 53.

    Ashida, Y., Gong, Z. & Ueda, M. Non-hermitian physics. Adv. Phys. 69, 249–435 (2020).

    ADS  Article  Google Scholar 

  54. 54.

    Trefethen, L. & Embree, M. Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators (Princeton University Press, 2020).

  55. 55.

    Kawabata, K., Shiozaki, K. & Ryu, S. Topological Field Theory of Non-Hermitian Systems. Phys. Rev. Lett. 126, 216405 (2021).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  56. 56.

    Ghatak, A., Brandenbourger, M., van Wezel, J. & Coulais, C. Observation of non-Hermitian topology and its bulk–edge correspondence in an active mechanical metamaterial. Proc. Natl Acad. Sci. 117, 29561–29568 (2020).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  57. 57.

    Lee, C. H. & Thomale, R. Anatomy of skin modes and topology in non-Hermitian systems. Phys. Rev. B 99, 201103 (2019).

    ADS  CAS  Article  Google Scholar 

  58. 58.

    Helbig, T. et al. Generalized bulk–boundary correspondence in non-Hermitian topolectrical circuits. Nat. Phys. 16, 747–750 (2020).

  59. 59.

    Yao, S. & Wang, Z. Edge states and topological invariants of non-Hermitian systems. Phys. Rev. Lett. 121, 086803 (2018).

    ADS  CAS  PubMed  Article  Google Scholar 

  60. 60.

    Okuma, N., Kawabata, K., Shiozaki, K. & Sato, M. Topological origin of non-Hermitian skin effects. Phys. Rev. Lett. 124, 086801 (2020).

    ADS  MathSciNet  CAS  PubMed  PubMed Central  Article  Google Scholar 

  61. 61.

    Hofmann, T. et al. Reciprocal skin effect and its realization in a topolectrical circuit. Phys. Rev. Res. 2, 023265 (2020).

    CAS  Article  Google Scholar 

  62. 62.

    Zhang, K., Yang, Z. & Fang, C. Correspondence between winding numbers and skin modes in non-Hermitian systems. Phys. Rev. Lett. 125, 126402 (2020).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  63. 63.

    Mostafazadeh, A. Physics of spectral singularities. (Kielanowski, P., et al. eds.) Geometric Methods in Physics, 145–165 (Springer International Publishing, Cham, 2015).

  64. 64.

    Bender, C. M. & Boettcher, S. Real spectra in non-Hermitian hamiltonians having PT symmetry. Phys. Rev. Lett. 80, 5243–5246 (1998).

    ADS  MathSciNet  CAS  MATH  Article  Google Scholar 

  65. 65.

    Whitaker, J. Microelectronics. Electronics Handbook Series (CRC Press, 2018).

  66. 66.

    Cui, H. et al. Three-dimensional printing of piezoelectric materials with designed anisotropy and directional response. Nat. Mater. 18, 234–241 (2019).

    ADS  CAS  PubMed  Article  Google Scholar 

  67. 67.

    Srivastava, A. Causality and passivity in elastodynamics. Proc. R. Soc. A: Math., Phys. Eng. Sci. 471, 20150256 (2015).

    ADS  MathSciNet  MATH  Article  Google Scholar 

  68. 68.

    Day, W. A. Time-reversal and the symmetry of the relaxation function of a linear viscoelastic material. Arch. Ration. Mech. Anal. 40, 155–159 (1971).

    MathSciNet  MATH  Article  Google Scholar 

  69. 69.

    Muhlestein, M. B., Sieck, C. F., Al’u, A. & Haberman, M. R. Reciprocity, passivity and causality in Willis materials. Proc. R. Soc. A: Math., Phys. Eng. Sci. 472, 20160604 (2016).

    ADS  MathSciNet  MATH  Article  Google Scholar 

Download references


The authors gratefully thank Prof. Hussein Nassar from University of Missouri for valuable discussions, and Michele Fossati for a critical reading of the manuscript. This work is supported by the Air Force Office of Scientific Research under Grant No. AF9550-18-1-0342 and AF 9550-20-0279 with Program Manager Dr. Byung-Lip (Les) Lee and the Army Research Office under Grant No. W911NF-18-1-0031 with Program Manager Dr. Daniel P Cole. V.V. was supported by the Complex Dynamics and Systems Program of the Army Research Office under grant W911NF-19-1-0268. C.S. was supported by the National Science Foundation Graduate Research Fellowship under Grant No. 1746045. Some of us benefited from participation in the KITP program on Symmetry, Thermodynamics, and Topology in Active Matter supported by Grant No. NSF PHY-1748958.

Author information




Y.C. and G.H. conceived the concept; X.L. conducted experiments; Y.C. and C.S. performed theoretical investigations; Y.C. performed numerical investigations; G.H. and V.V. supervised the research; all the authors discussed the results; all authors wrote the manuscript and interpreted the results.

Corresponding authors

Correspondence to Vincenzo Vitelli or Guoliang Huang.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Peer review information Nature Communications thanks Muamer Kadic and the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Chen, Y., Li, X., Scheibner, C. et al. Realization of active metamaterials with odd micropolar elasticity. Nat Commun 12, 5935 (2021).

Download citation


By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.


Quick links