††thanks: These three authors contributed equally††thanks: These three authors
contributed equally††thanks: These three authors contributed equally

# Backscattering-free edge states below all bands in two-dimensional auxetic
media

Wenting Cheng<EMAIL_ADDRESS>Department of Physics, University of
Michigan, Ann Arbor, MI 48109, USA  Kai Qian<EMAIL_ADDRESS>Department of
Mechanical and Aerospace Engineering, University of California San Diego, La
Jolla, CA 92093, USA  Nan Cheng<EMAIL_ADDRESS>Department of Physics,
University of Michigan, Ann Arbor, MI 48109, USA  Nicholas Boechler
<EMAIL_ADDRESS>Department of Mechanical and Aerospace Engineering,
University of California San Diego, La Jolla, CA 92093, USA  Program in
Materials Science and Engineering, University of California San Diego, La
Jolla, CA 92093, USA  Xiaoming Mao<EMAIL_ADDRESS>Department of Physics,
University of Michigan, Ann Arbor, MI 48109, USA  Kai Sun<EMAIL_ADDRESS>Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

###### Abstract

Unidirectional and backscattering-free propagation of sound waves is of
fundamental interest in physics, and highly sought-after in engineering.
Current strategies utilize topologically protected chiral edge modes in
bandgaps, or complex mechanisms involving active constituents or nonlinearity.
Here we propose a new class of passive, linear, one-way edge states based on
spin-momentum locking of Rayleigh waves in two-dimensional media in the limit
of vanishing bulk modulus, which provides $100\%$ unidirectional and
backscattering-free edge propagation at a broad range of frequencies instead
of residing in gaps between bulk bands. We further show that such modes are
characterized by a new topological winding number that is analogous to
discrete angular momentum eigenvalues in quantum mechanics. These passive and
backscattering-free edge waves have the potential to enable a new class of
phononic devices in the form of lattices or continua that work in previously
inaccessible frequency ranges.

##  I Introduction

Topologically-protected chiral edge modes have been the focus of extensive
studies over the past few decades, where nontrivial topological structures of
bulk waves give rise to edge and interface waves exhibiting features such as
one-way transport, quantized transport coefficients, and robustness against
disorder [1, 2, 3, 4, 5, 6, 7, 8]. These intriguing phenomena have been
observed in a range of systems from quantum electronic topological materials
to classical phononic or photonic metamaterials [9, 10, 11, 12, 13, 14, 15].
However, such exotic physics is subject to one constraint: these edge modes
(including quantum Hall and spin, valley, or mirror Hall types) can only arise
in frequency gaps between bulk bands.

Is it possible to break this constraint? For quantum Hall edge states, this is
impossible, because they are associated with a gauge anomaly, which cannot be
achieved below the lowest band (or above the highest) [1, 16]. In contrast,
one-way transport of topological edge modes of spin, valley, or mirror Hall
types come from (pseudo) spin-momentum locking [1, 2, 3, 4, 5, 6, 7, 8]. Thus,
if the same type of spin-momentum locking can be realized through a different
mechanism, it is in principle possible to achieve similar physical properties
without the in-gap constraint.

For sound waves, partial spin-momentum locking has long been known [17, 18,
19, 20, 21, 22]. One example is Rayleigh waves [23] on the surface of a solid,
where the surface exhibits elliptical motions (also known as “ground roll” in
seismology [24]). Interestingly, waves propagating in opposite directions
exhibit opposite rolling polarizations. The physics behind this is the
breaking of inversion symmetry. Rotations and linear motions are characterized
by two different types of physical quantities: pseudo-vectors (e.g., angular
momentum) and vectors (e.g., linear momentum), which can only couple when
inversion symmetry is broken, due to their opposite parity. The presence of a
physical boundary serves this role of breaking symmetry, and thus allows these
two different types of motions to be coupled. Although such locking naturally
arises for surface modes, this is generally a partial locking, in contrast to
topological edge modes of quantum spin [25, 26, 27] or valley [28, 29, 30, 31,
32, 33, 34, 35] Hall systems, where the locking is 100%. As with all other
elliptically polarized waves, Rayleigh waves contain both left- and right-
circularly polarized components, just with different amplitudes and phases.

Figure 1: One-way transport and spin-momentum locking of Rayleigh waves.
Simulated spatial energy distribution for 2D (a) and 3D (b) continua with
different $B_{3D}/\mu_{3D}$ values, using a prescribed particle displacement
at boundary point $r$ following $\psi_{k}(r)$. Green arrows denote the spin
direction (parallel to the surface) and blue arrows denote the strong
propagation direction. The relative displacement magnitude of the red points
marked in (a) and (b) is evaluated to measure $|C_{-k}/C_{k}|$, shown as
yellow points in (c) and (d). The shear modulus ($\mu_{3D}$) is kept constant
at $0.01$ GPa throughout all the simulations and the bulk modulus ($B_{3D}$)
varied. (c,d) Dependence of the one-way propagation metric on the elastic
moduli ratio in 2D (c) and 3D (d). The blue curves are analytic results, and
the red lines are bounds. (e-f) The dependence of polarization of the Rayleigh
wave at various depths when $B_{3D}/\mu_{3D}=0$ (normalized using the decay
length $\kappa_{t}$ of the transverse polarized component of the edge mode),
depicted as the ratio of the short to long axes ($b/a$) of the ellipse of
polarization (analytic). The inset illustrates the polarization (blue
ellipses) and $\vec{\psi}_{-k}(r)^{*}\cdot\vec{\psi}_{k}(r)$ (color) of the
Rayleigh waves at different depth for 2D (e) and 3D (f). Sizes of the ellipses
are proportional to the amplitude of the Rayleigh waves at the ellipses’
center. (g-h) Dependence of the longitudinal ($C_{l}$), transverse ($C_{t}$),
and Rayleigh ($C_{R}$) wavespeeds on the Poisson’s ratio ($\nu$) in 2D (g) and
3D (h) (analytic). (i) A summary of the one-way transport metric in 2D and 3D
media. (j) Dispersion of longitudinal, transverse, and Rayleigh waves in 2D
continua at $B_{3D}/\mu_{3D}=0.1$ (analytic). The shaded area represents bulk
modes.

In this article, we show that 100% spin-momentum locking can be achieved using
edge waves of two-dimensional (2D) isotropic continua/lattices in the limit of
vanishing bulk modulus $B\to 0$ (which implies auxeticity, where Poisson’s
ratio $\nu=-1$). These modes are below the bulk wave bands (i.e. lower
wavespeeds), and the system doesn’t need gaps between bulk bands. We show that
the spin-momentum locking shown herein is beyond that of classical Rayleigh
waves—which applies to semi-infinite bodies with flat surfaces —and supports
backscattering-free propagation (Fig. 1) under _arbitrary geometries_ ,
regardless of microscopic details (Fig. 2). We verify this theory in
experiments using auxetic Maxwell lattices. Although they have no connection
to valley, spin, or mirror Chern numbers, these waves are protected by the
Kane-Lubensky index defined for Maxwell lattices [36, 37, 38] and continua
[39] when the system is close to the $B\to 0$ limit. Furthermore, we show that
a new topological winding number can be defined for these waves, which plays
the role of discrete angular momentum eigenvalues, generalized to arbitrary
geometries, and protects the backscattering-free propagation of these waves.

##  II Results

###  II.1 One-way transport and spin-momentum locking of Rayleigh waves

How to excite a wave that only travels in $+k$ and not $-k$? In principle this
can be achieved by applying forces on the media to exactly follow the $+k$
mode in space-time. A more interesting question, however, is how to create
one-way propagation when one can only apply a periodic force $F$ on one point
$r$ on the edge of the medium. For simplicity, we first consider the simple
case where the material is a semi-infinite plane, and we will discuss general
geometries next. This point source will generate a response in all modes at
(angular) frequency $\omega$, and the amplitude of each mode
$\vec{\psi}_{\alpha}$ will be proportional to
$\vec{\psi}_{\alpha}(r)^{*}\cdot\vec{F}$, where ∗ denotes complex conjugate
(see Supplementary Materials (SM) for details). Given that $r$ is on the edge,
this overlap is much greater for edge localized modes than bulk waves.
Furthermore, these edge modes appear in pairs of $\vec{\psi}_{k}$ and
$\vec{\psi}_{-k}$, where $k$ is the momentum (wavevector) along the edge, and
a generic point source $\vec{F}$ excites both $\vec{\psi}_{k}$ and
$\vec{\psi}_{-k}$.

To approach “one-way transport” in the $+k$ direction, one needs to maximize
the ratio of the resulting amplitudes of $\vec{\psi}_{k}$ and
$\vec{\psi}_{-k}$. This can be achieved by choosing
$\vec{F}=\vec{\psi}_{k}(r)$. A metric for asymmetric edge transport can then
be defined as

$\Big{|}\frac{C_{-k}}{C_{k}}\Big{|}=\Big{|}\frac{\vec{\psi}_{-k}(r)^{*}\cdot\vec{F}}{\vec{\psi}_{k}(r)^{*}\cdot\vec{F}}\Big{|}=\Big{|}\frac{\vec{\psi}_{-k}(r)^{*}\cdot\vec{\psi}_{k}(r)}{\vec{\psi}_{k}(r)^{*}\cdot\vec{\psi}_{k}(r)}\Big{|},$
(1)

which vanishes in the limit of one-way transport. This ratio not only
describes how external forces at the boundary point $r$ excite left- vs right-
going waves, but also describes impurity scattering as the wave propagates in
the bulk close to the edge, where impurities creates new sources $F$. One
example where $|C_{-k}/C_{k}|=0$ is the quantum spin Hall effect, where spin-
polarized topological edge states lead to orthogonal states for $k$ and $-k$.

Here, we consider instead Rayleigh waves,

$\vec{\psi}_{\pm k}(x,y)=e^{i(\pm kx-\omega
t)}\left(\frac{1}{2}\xi\sqrt{1+\frac{\mu}{B}},\pm\frac{1}{4}i\xi(2-\xi^{2})\sqrt{\frac{1+(\mu/B)}{1-\xi^{2}}}\right)^{T}$
(2)

at the edge ($y=0$) of a semi-infinite plane, and $\xi=\xi(B/\mu)$ is the
ratio of Rayleigh to transverse wave speeds and a function of the ratio of
bulk modulus $B$ and shear modulus $\mu$ (see SM for details). This pair of
waves are linear combinations of longitudinal and transverse waves, satisfying
stress-free boundary conditions at $y=0$, and their longitudinal and
transverse components decay into the bulk with different decay lengths [23].
Here $B$ and $\mu$ are moduli for 2D isotropic continuum elasticity. To make
this discussion relevant for experiments, we relate them to 3D elastic moduli
(subscripts 3D) of the material that consists this 2D plate (under plane
stress condition), which reads
$B=3B_{3D}\mu_{3D}/(B_{3D}+\frac{4}{3}\mu_{3D})$ and $\mu=\mu_{3D}$, and use
3D moduli in our computational results.

Interestingly, at the limit of $B/\mu\to 0$, (equivalently $\xi\to 0$), this
pair of waves become

$\vec{\psi}_{\pm k}(x,y)=\frac{1}{\sqrt{2}}e^{i(\pm kx-\omega t)}\left(1,\pm
i\right)^{T},$ (3)

which are circularly polarized. In this case, $|C_{-k}/C_{k}|=0$, leading to
$100\%$ one-way transport. The general case of $B/\mu>0$ leads to ellipsoidal
polarization and $|C_{-k}/C_{k}|>0$. In the limit of $B/\mu\to\infty$,
$|C_{-k}/C_{k}|\to 0.4203<1$, meaning that the propagation along $\pm k$ is
still biased (Fig. 1(a,c) with simulation described in Methods).

The circular polarization of the Rayleigh waves in the limit of $B/\mu\to 0$
suggests that a spin-momentum locking mechanism can be defined for this one-
way transport. Indeed, similar to the spin angular momentum of photons, where
photons of $\pm 1$ spins correspond to quanta of right (R) and left (L)
polarized electromagnetic waves, a spin angular momentum can be defined for
phonons [40]. In the SM, we show that the rotational symmetry of the media
leads to conserved angular momentum density, following Neother’s theorem, and
this angular momentum $\vec{j}=\vec{l}+\vec{s}$ can be expressed as a sum of
an orbital ($l$) and spin ($s$) angular momentum density of the wave, where

$\displaystyle\vec{l}$
$\displaystyle=\vec{r}\times\vec{p},\quad\vec{s}=\rho\vec{u}\times(\partial_{t}\vec{u}),$
(4)

$\rho$ is the mass density, and
$\vec{p}=-\rho\partial_{t}u_{i}\vec{\gradient}u_{i}$ is the density of linear
momentum. Following this, the Rayleigh waves carry physical spin angular
momentum pointing along $\hat{z}$ (normal to the 2D plane), and under standard
canonical quantization, the spins of these phonons take $\pm 1$.

What about 3D Rayleigh waves? A similar calculation can be done, giving a pair
of Rayleigh waves for any propagation direction on the 2D surface of a semi-
infinite 3D solid (see SM for details) [41, 17]. However, the asymmetric edge
transport ratio never reaches zero ($|C_{-k}/C_{k}|>0.0509$, Fig. 1(b,d)).
Thus, $100\%$ one-way transport, in the manner defined herein using Rayleigh
waves, is only achieved for 2D solids. More importantly, 3D Rayleigh waves are
far from “perfect” in terms of spin-momentum locking, even in the $B\to 0$
limit. This is because in a 3D solid, although the Rayleigh-wave mode shape is
close to circular polarization at the surface layer, as one moves into the
bulk, the Rayleigh-wave mode shape quickly becomes elliptical [17] (Fig.
1(f)), i.e., $|C_{-k}/C_{k}|>0$ in the bulk, and thus any impurity would
result in backscattering. In contrast, 2D Rayleigh waves have
$|C_{-k}/C_{k}|=0$ at all depths (Fig. 1(e)), making it immune to
backscattering.

Figure 2: Backscattering-free propagation of one-way edge waves in arbitrary
geometries. (a) Simulated energy distribution for 2D continua with different
values of $B_{3D}/\mu_{3D}$: $0.01,0.02,0.05,$ and $0.1$ (at different times).
The clockwise actuation (Gaussian wave packet) excites right-going edge waves.
To characterize the backscattering, a sharp defect is created at the boundary
(black arrow). (b) Energy of points $A$ (actuation point), $B$, $C$, and $D$
(each marked in (a)) as functions of time. (c) Energy ratio between points $C$
and $A$ (cumulative over time) at different $B_{3D}/\mu_{3D}$ demonstrates
vanishing backscattering in the limit of $B/\mu\to 0$.

###  II.2 Backscattering-free one-way transport in arbitrary geometries

The Rayleigh wave formulation discussed above offers a concise picture to
characterize how spin-momentum locking arises in these edge waves in the
simple geometry of semi-infinite media. However, the observed one-way
transport appears to be robust in materials with very irregular geometries
with no backscattering in the limit of $B/\mu\to 0$ (Fig. 2). To reveal the
physical origin of this phenomenon, we propose a new formulation based on
mapping of zero modes in the $B/\mu\to 0$ limit in 2D to complex analytic
functions and perturbatively extending them to low-frequency edge waves at
small $B/\mu$. By doing this, we show that such spin-momentum locking is a
general phenomenon in 2D, small $B/\mu$ materials of any shape, and these edge
waves are characterized by a new topological winding number that plays the
role of discrete angular momentum eigenvalues despite the arbitrary geometry.

We first use the complex plane to rewrite the 2D elasticity problem, by
defining coordinate $z=x+iy$ and displacement field $\psi=u_{x}+iu_{y}$. In
the limit of $B/\mu\to 0$, the wave equation and boundary conditions can then
be written as (see SM for details)

$\displaystyle\rho\frac{\partial^{2}\psi(z,z^{*},t)}{\partial t^{2}}$
$\displaystyle=4\mu\frac{\partial^{2}\psi(z,z^{*},t)}{\partial z\partial
z^{*}},\quad z\in\Omega,$ (5)
$\displaystyle\frac{\partial\psi(z,z^{*},t)}{\partial z^{*}}$
$\displaystyle=0,\quad z\in\partial\Omega,$

where $\Omega$ and $\partial\Omega$ represent the bulk and edge of the system,
respectively. Solutions of Eq. (5) can be obtained via separation of variables
$\psi_{\pm}(z,z^{*},t)=\phi(z,z^{*})\,e^{\pm i\omega t}$. Because $e^{+i\omega
t}$ rotates the displacement vector counterclockwise and $e^{-i\omega t}$
rotates it clockwise over time, these two solutions correspond to spin up and
down states. This equation has $\omega=0$ solutions $\phi_{0}$ given by
$\partial\phi_{0}/\partial z^{*}=0$, i.e., arbitrary bounded analytic
functions $\phi_{0}(z)$. These solutions as static zero modes have been
characterized in Ref. [42], and they can be derived from the fact that zero
modes in 2D materials at $B/\mu=0$ are conformal transformations, which map to
analytic functions [42, 43].

Here we are interested in the dynamics of waves at small but nonzero $B/\mu$,
which can be solved perturbatively from the zero modes. Same as the
perturbation theory for the Schrödinger equation [44], to the leading order,
small perturbations only modify the frequencies (eigenvalues), while
wavefunctions remain analytic functions. Thus, at small $B$, these zero modes
acquire finite frequencies $\omega\propto\sqrt{B/\mu}$,

$\psi_{\pm}(z,z^{*},t)=\phi_{0}(z)\,e^{\pm i\omega t}.$ (6)

Figure 3: Eigenfrequencies and phase-winding of one-way edge waves in
arbitrary geometries. (a) Computed (green dots) and analytic (red line,
following Eq. (10)) eigenfrequencies of a 2D media of irregular shape plotted
against the mode index at $B_{3D}/\mu_{3D}=0.01$. Modes above the blue dashed
line are bulk modes. (b) Profiles of modes $n=4,6,8,10$ (corresponding to
winding numbers $N_{0}=1,2,3,4$) as examples of waves described by Eq. (6).
Pink dots mark $\phi=0$ points in the bulk, the number of which equals the
winding number $N_{0}$ (Eq. (S97)). The black arrows represent the
displacement vectors, which rotate $N_{0}$ times around the boundary, and
vanish at $N_{0}$ points in the bulk.

How do these waves propagate in space and time? The irregular geometry we
consider here exhibits no translational or rotational symmetries, so linear
momentum $k$ and orbital angular momentum $l$ are no longer conserved
quantities. Instead, a winding number can be defined by applying the argument
principle on the analytic functions $\phi_{0}$, which we suggest plays the
role of “quantized” orbital angular momentum eigenvalues,

$\frac{1}{2\pi i}\int_{\partial\Omega}d\ln(\phi_{0}(z))=\frac{1}{2\pi
i}\int_{\partial\Omega}\frac{\phi^{\prime}_{0}(z)}{\phi_{0}(z)}dz=N_{0}-N_{\infty},$
(7)

where $N_{0}$ and $N_{\infty}$ are the numbers of zeros and poles of $\phi(z)$
inside the bounded domain $\Omega$. The function $\phi_{0}(z)$ being bounded
in $\Omega$ implies $N_{\infty}=0$. On the other hand, edge localized modes
$\phi_{0}(z)$ must have zeros in $\Omega$ (Rouché’s theorem, as detailed in
SM), so $N_{0}\geq 1$. This indicates that all zero modes $\phi_{0}$ wind its
phase by at least one cycle of $2\pi$ as one goes along the boundary
counterclockwise. For a finite 2D material at $B/\mu=0$, there are infinitely
many such zero modes, and they each exhibit a winding number $N_{0}$, which is
also the number of points in the bulk where the amplitude of the zero mode
reaches zero (counting multiplicity), as shown in Fig. 3(b). In addition, zero
modes with $N_{0}$ zeros can be expanded as analytic functions as

$\phi_{0,n}(z)=a_{N_{0}}z^{N_{0}}+...+a_{1}z+a_{0}.$ (8)

It is thus clear that zero modes with larger $N_{0}$ are more localized on the
edge, which is consistent with predictions from conformal symmetry [42].

This immediately leads to spin-momentum locking at small $B/\mu$: constant
phase points of $\psi_{+}(z,t)=\phi_{0}(z)\,e^{+i\omega t}$ move clockwise and
constant phase points of $\psi_{-}(z,t)=\phi_{0}(z)\,e^{-i\omega t}$ move
counterclockwise in time. Therefore, this winding number is the extension of
discrete angular momentum eigenvalues from perfect circular geometry to
arbitrary geometries in the limit of $B/\mu\to 0$. It is worth noting that
such an extension is not generally defined for any 2D material, where wave
eigenstates are usually affected by irregular shapes of the boundary and do
not exhibit a clear propagation direction or momentum (see examples in Fig.
2). It is the $B/\mu\to 0$ limit that permits such a clean definition of an
orbital angular momentum.

The two waves $\psi_{\pm}(z,t)$ are orthogonal, satisfying the one-way
transport criterion $|C_{-k}/C_{k}|=0$ defined in Eq. (1). This can be seen by
returning from the complex plane notation to the conventional wave notation
with $x,y$ components, where

$\phi_{0}(z)\,e^{\pm i\omega t}\to\Re{\phi_{0}(x+iy)e^{i\omega t}(1,\pm
i)^{T}},$ (9)

which are circularly polarized, with left- (L) and right- (R) waves orthogonal
to one another. This orthogonality is immune to shape irregularities of the
boundary, as it stems from exact zero modes in the $B/\mu\to 0$ limit, and
does not rely on any notion of plane waves or Rayleigh waves of semi-infinite
geometry. Thus, these circularly polarized waves approach backscattering-free
as $B/\mu\to 0$ for any geometry (Fig. 2).

The spectrum of these edge modes also resemble waves characterized by discrete
angular momentum. For 2D materials with shapes of no particular symmetry, to
leading order in $B/\mu$, these modes have frequencies

$\omega_{n}\simeq\frac{2\pi N_{0}(n)v_{phase}}{L}\simeq\frac{\pi
nv_{phase}}{L}$ (10)

where $n$ labels the modes, $L$ is the perimeter of the bounded domain
$\Omega$, $v_{phase}=\xi(B/\mu)\sqrt{\mu/\rho}$ is the phase velocity of
Rayleigh waves of the given material. The winding number is related to the
mode index by $N_{0}(n)=\lceil(n-3)/2\rceil$ where $\lceil\cdot\rceil$ is the
ceiling function (note that each zero mode at $B=0$ generates a pair of edge
modes at $B>0$, and there are three rigid body modes that remain at
$\omega=0$), and this reduces to $N_{0}(n)\simeq n/2$ at large $n$, as shown
in Fig. 3(a). Frequencies of these edge modes extend to the point of the first
bulk mode, which is at the scale of $\sqrt{\mu/\rho}/L$, noting that the
longitudinal and transverse wavespeeds converge as $B/\mu\to 0$ (see Fig.
1(g)).

Figure 4: Experimental verification of spin-momentum locking and one-way edge
waves in an auxetic Maxwell lattice. (a) The self-dual twisted kagome lattice
[45]. The red star marks the actuation point. Red and green solid circles
denote the triangles’ centroids, which are used for image tracking. The 2D
lattice has primitive basis vectors $(\vec{a}_{1},\vec{a}_{2})$ with a lattice
constant of $a=5.71$ cm. Insets: (left) circular actuator (with a radius of
$5$ mm), and linear actuator; and (right) unit-cell geometry (gray: triangles
with ideal pin-joint; blue: triangles with beam-like elastic connectors) and
assembly diagram. (b) Super-cell band structures simulated in FEM for the
elastic lattice in (a) with edge bands marked in red. (c) $|{C_{-k}}/{C_{k}}|$
of the lowest edge band. Shaded area indicates $|{C_{-k}}/{C_{k}}|<0.1$.
Dashed line marks $k=0.28{\pi}/{a}$. (d) Two simulated edge mode profiles at
$k=0.28{\pi}/{a}$, also showing the geometry simulated in (b). (e-g) Time
domain FEM simulations and (h-j) experimental results. Black circles and
arrows indicate the location and the direction of the excitation,
respectively. Solid circle sizes are proportional to velocity magnitudes and
black lines indicate velocity vector directions. Orange (light blue) color
indicates the sign of the rotation of the triangle centroids’ velocity vector,
corresponding to counterclockwise (clockwise) spins. ${E_{L}}/{E_{R}}$
indicates the energy ratio between left (L) and right (R) propagating modes.

###  II.3 One-way transport in twisted kagome lattices

To demonstrate this spin-dependent edge transport experimentally, we choose
the twisted kagome lattice, which is known to exhibit $B/\mu=0$ in the ideal
pin-joint limit (no bending stiffness between triangles) [42]. The lattice we
consider herein (Fig. 4(a)) deviates from this limit, as we use beam-like
elastic connectors instead of pin joints. We measure the Poisson’s ratio of
this lattice to be $\nu\simeq-0.488$ (see SM), corresponding to $B/\mu=0.344$.

Pairs of branches of exponentially localized edge zero modes are known to
arise in twisted kagome lattices under mixed (periodic horizontal and open
vertical) boundary conditions, at opposite open edges [42]. They show up as
the red branch in (Fig. 4(b)) via FEM super-cell analysis of the elastic
kagome lattice, which are actually two nearly overlapping branches At
wavelengths approaching the width of the simulated domain (Fig. 4(d)) or low
$k$, the two localized modes generally hybridize into Lamb-like modes [46]
(see SM). The effect of this hybridyzation can be seen in Fig. 4(c), where
$|{C_{-k}}/{C_{k}}|$ exhibit a peak near $k\to 0$. An important exception to
this is the $B/\mu\to 0$ case, wherein the two edge modes are exactly
orthogonal and thus, do not hybridize and retain their asymmetric behavior.
When $k=0.28\pi/a$, the decay length is short (about $1.4$ unit cells) and the
hybridyzation is weak, so we choose the frequency of this pair of modes ($5.5$
Hz, Fig. 4(d)) for FEM and experiment of edge waves.

We verify spin-momentum locking and one-way edge wave propagation by actuating
at the centroid of one triangle at the edge. With circularly polarized
excitation, we observe spin-dependent one-way edge wave propagation in both
simulation (Fig. 4(e,f)) and experiment (Fig. 4(h,i)). With linear actuation,
we observe nearly symmetric edge-wave propagation, wherein the wave’s spin is
matched to the direction of propagation (Fig. 4(g,j)). For each excitation
case, we also calculate the energy ratio ${E_{L}}/{E_{R}}$ (Fig. 4(e-j))
between L and R propagating edge modes. The energy for each mode is calculated
by integrating the average velocity squared (over $0.01-0.15$ s) of triangles
on the edge to the left and right of the excitation point (see SM for specific
triangles). Both the spatial response and the ${E_{L}}/{E_{R}}$ ratio shows a
close agreement between simulation and experiment (Fig. 4).

##  III Discussion

We have reported a new type of one-way and backscattering-free edge wave based
on spin-momentum locking in 2D auxetic materials with moduli $B/\mu\to 0$.
These modes are protected by a new topological winding number that generalizes
discrete angular momentum eigenvalues to arbitrary geometries. Unlike other
topologically protected edge modes, these protected edge modes are not
required to reside in a gaps between bulk bands and instead exist for a broad
range of frequencies. This is important, because, compared to the bandgap-tied
mechanisms (quantum Hall, quantum spin, valley, or mirror Hall types): i) it
opens previously inaccessible ranges of low frequencies for topologically
protected transport, and ii) it does not necessitate any precise periodic or
locally-resonant structures, and therefore expands the range of usable
frequencies in a given material, simplifies fabrication, and minimizes
dissipation.

As per its broadband functionality, the one-way edge mode transport mechanism
shown herein is applicable to essentially any auxetic lattice or continuum in
the $B/\mu\to 0$ limit, regardless of microscopic details. As auxetic
materials have been studied for decades [47], there exists a broad potential
set of candidate materials, including top-down manufactured structural
materials, such as 3D printed re-entrant honeycomb lattices [48], and
processed bottom-up manufactured materials, such as thermomechanically
compressed polyurethane foams [49]. Some naturally occurring auxetic materials
have also been found, exhibiting anisotropic negative Poisson’s ratios (albeit
still far from $\nu\to-1$), such as $\alpha$-cristobalite and some metals
[50]. Our theory predicts that asymmetric edge transport exists in these
materials, although perfect backscattering-free one-way transport calls for
proximity to $\nu\to-1$ and low dissipation.

Finally, we emphasize that the “spin” of these edge waves are real, in the
sense that they are directly connected to angular momentum, instead of the
“pseudo spins” typically unitized to introduce the spin Hall effect in
phononic systems [25, 26, 27]. Because of this physical spin, it has been
previously suggested that such circularly polarized waves may be able to
couple to spin degrees of freedom of other signals [18] such as photons [51,
52, 53] and magnons [54]. Similarly, we speculate that such coupling and
topological protection may find use in new types of optomechanical device
strategies [55].

##  IV Methods

###  IV.1 Maxwell lattice manufacturing

The self-dual kagome lattice with $90^{\circ}$ rotated solid equilateral
triangles connected via long and slender beams shown in Fig. 4(a) was designed
to be assembled together from repeating sub-lattice elements. The sub-lattice
elements (design shown in SM) were machined from $6.2$ mm thick polycarbonate
sheets via desktop CNC mill then assembled together via press fits.

###  IV.2 Experimental edge wave characterization

A misaligned shaft coupling ($5$ mm central offset) is installed on a stepper
motor to induce a circular displacement path to the actuated triangle. The
motor is controlled by a RAMBo board via MATLAB R2022b. A ball bearing is used
at the end of the shaft to minimize rotation of the actuated triangle. The
excitation frequency is $5.5$ Hz. An electrodynamic exciter (The Modal Shop,
Smartshaker, K2007E01) is driven sinusoidally at $5.5$ Hz by a signal
generator (Tektronix, AFG3022C) to perform the linear actuation. A digital
camera records the displacements of the lattice under actuation from above at
a sampling rate of $240$ Hz. The locations of the red and green markers in
Fig. 4(a) are tracked from the recorded image sequence by using the
imfindcircles command in MATLAB R2022b. The spatial measurement resolution is
calculated to be $0.37$ mm for the circular actuation cases, and $0.48$ mm for
the linear actuation cases. The location coordinates are smoothed in the time
domain with a Gaussian-weighted moving average filter using the smoothdata
command with a window size of $20$ in MATLAB R2022b.

###  IV.3 FEM simulation

The FEM simulations are performed using the COMSOL Multiphysics Structural
Mechanics Module. The 2D and 3D continua depicted in Figs. 1, 2, and 3 use
user-defined material properties with a mass density of $960$ kg/m3 and a
shear modulus of $0.01$ GPa, and differing bulk moduli, with extra fine mesh
(Fig. 1: Mapped; Figs. 2,3: Free Triangular, minimum element size $0.0075$ cm
and maximum element size $2$ cm). The plane stress approximation has been
employed in the 2D simulations. The dimensions of these systems are
approximately $1$ m $\times~{}1$ m for 2D ($1$ m $\times~{}1$ m $\times~{}1$ m
for 3D). Figures 1 and 2 are simulated in the time domain, with a time step of
$0.0001$ s, and the simulations are carried out from $0$ s to $0.05$ s. Figure
3 is simulated in the frequency domain. The governing equation is given by
$\rho\frac{d^{2}u}{dt^{2}}=\nabla\sigma+\vec{F(r,t)}$, where $\vec{F(r,t)}$
represents the time-dependent excitation at a boundary point, expressed in the
form of $\vec{\psi}_{k}(r,t)$ (noting that a prescribed displacement is
applied). In the 2D continua, the excitation spin direction is out-of-plane,
while in the 3D continua, it is parallel to the surface. In Fig. 1, a free
boundary condition is applied to the boundary with excitation, and low-
reflecting boundary conditions are employed on the remaining boundaries. This
setup allows us to evaluate the ratio $|C_{-k}/C_{k}|$. This ratio is
calculated by taking the ratio of the square root of the energy at the marked
points, which are located on opposite sides of the excitation, over an
extended period of time. In Figs. 2 and 3, free boundary conditions are
applied. The excitation frequency is adjusted to $1000$ Hz for different
$B_{3D}/\mu_{3D}$ in 3D. The excitation frequency varies depending on the
value of $B_{3D}/\mu_{3D}$ in 2D, with frequencies of $500$ Hz, $1500$ Hz, and
$2000$ Hz used for $B_{3D}/\mu_{3D}=0.01$, $B_{3D}/\mu_{3D}=0.1$, and
$B_{3D}/\mu_{3D}=1$, respectively. The simulated results are smoothed in the
time domain following the same process as the experimental data.

The Maxwell lattices shown in Fig. 4 are filled with polycarbonate material,
which has a mass density of $1190$ kg/m3, a Young’s modulus of $0.751$ GPa,
and a Poisson’s ratio of $0.3182$. In the super-cell analysis, periodic
boundary conditions are applied along $y$ direction to conduct parametric
sweeps of wave vectors. The ratio $|{C_{-k}}/{C_{k}}|$ in Fig. 4 (c) is
computed using equation 1, utilizing the eigenvectors (displacement)
$\vec{\psi}{k}(r)$ and $\vec{\psi}{-k}(r)$ obtained for the first band. Here,
$r$ represents the centroid of the first left triangle, as the energy of the
first band localizes at the left edge when $k$ is not in the small limit. In
Fig.4 (e-g), the excitation is applied at a single triangle centroid on the
boundary, with a frequency $f=\omega/(2\pi)=5.5$ Hz. This frequency
corresponds to the lowest edge mode frequency in the super-cell analysis when
$k=\pm 0.28\pi/a$, where $a$ is the lattice constant.

## Acknowledgement

W. C., N. C., K. S., and X. M. acknowledge support from the Office of Naval
Research (MURI N00014-20-1-2479). K. Q., N. B., and X. M. acknowledge support
from the US Army Research Office (Grant No. W911NF-20-2-0182).

## Data and code availability

The data and code that support the findings of this study are available from
the corresponding authors on reasonable request.

## Author contributions

W. C. performed the full-wave simulations, K. Q. performed the fabrication,
and experimental characterization. N. C., X. M. and K. S. developed the
mathematical framework. W. C. and K. Q. analyzed and interpreted the data. All
authors wrote the manuscript. N. B., X. M., and K. S. supervised the project.

## Competing interests

The authors declare no competing interests.

## Supplemental Material

The supplemental materials contain detailed mathematical discussion on 2D and
3D Rayleigh wave mode shape, criterion for perfectness of Rayleigh wave one-
way transport in continuum systems and discrete lattices, the spin and
momentum of elastic waves as well as full spin and momentum locking of
Rayleigh wave in 2D material with $B/\mu\ll 1$.

## References

  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized hall conductance in a two-dimensional periodic potential, Physical review letters 49, 405 (1982).
  * Sheng  _et al._ [2005] L. Sheng, D. Sheng, C. Ting, and F. Haldane, Nondissipative spin hall effect via quantized edge transport, Physical review letters 95, 136602 (2005).
  * Prodan [2009] E. Prodan, Robustness of the spin-chern number, Physical Review B 80, 125327 (2009).
  * Sheng  _et al._ [2006] D. Sheng, Z. Weng, L. Sheng, and F. Haldane, Quantum spin-hall effect and topologically invariant chern numbers, Physical review letters 97, 036808 (2006).
  * Teo  _et al._ [2008] J. C. Teo, L. Fu, and C. Kane, Surface states and topological invariants in three-dimensional topological insulators: Application to bi 1- x sb x, Physical Review B 78, 045426 (2008).
  * Takahashi and Murakami [2011] R. Takahashi and S. Murakami, Gapless interface states between topological insulators with opposite dirac velocities, Physical review letters 107, 166805 (2011).
  * Zhang  _et al._ [2011] F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. MacDonald, Spontaneous quantum hall states in chirally stacked few-layer graphene systems, Physical review letters 106, 156801 (2011).
  * Ezawa [2013] M. Ezawa, Spin valleytronics in silicene: Quantum spin hall–quantum anomalous hall insulators and single-valley semimetals, Physical Review B 87, 155415 (2013).
  * Haldane and Raghu [2008] F. D. M. Haldane and S. Raghu, Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry, Physical review letters 100, 013904 (2008).
  * Raghu and Haldane [2008] S. Raghu and F. D. M. Haldane, Analogs of quantum-hall-effect edge states in photonic crystals, Physical Review A 78, 033834 (2008).
  * Wang  _et al._ [2008] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Reflection-free one-way edge modes in a gyromagnetic photonic crystal, Physical review letters 100, 013905 (2008).
  * Wang  _et al._ [2009] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Observation of unidirectional backscattering-immune topological electromagnetic states, Nature 461, 772 (2009).
  * Prodan and Prodan [2009] E. Prodan and C. Prodan, Topological phonon modes and their role in dynamic instability of microtubules, Physical review letters 103, 248101 (2009).
  * Nash  _et al._ [2015] L. M. Nash, D. Kleckner, A. Read, V. Vitelli, A. M. Turner, and W. T. Irvine, Topological mechanics of gyroscopic metamaterials, Proceedings of the National Academy of Sciences 112, 14495 (2015).
  * Wang  _et al._ [2015] P. Wang, L. Lu, and K. Bertoldi, Topological phononic crystals with one-way elastic edge waves, Phys. Rev. Lett. 115, 104302 (2015).
  * Haldane [1988] F. D. M. Haldane, Model for a quantum hall effect without landau levels: Condensed-matter realization of the ”parity anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * Lipsett and Beltzer [1988] A. W. Lipsett and A. I. Beltzer, Reexamination of dynamic problems of elasticity for negative poisson’s ratio, The Journal of the Acoustical Society of America 84, 2179 (1988).
  * Long  _et al._ [2018] Y. Long, J. Ren, and H. Chen, Intrinsic spin of elastic waves, Proceedings of the National Academy of Sciences 115, 9951 (2018).
  * Shi  _et al._ [2019] C. Shi, R. Zhao, Y. Long, S. Yang, Y. Wang, H. Chen, J. Ren, and X. Zhang, Observation of acoustic spin, National Science Review (2019).
  * Bliokh and Nori [2019] K. Y. Bliokh and F. Nori, Transverse spin and surface waves in acoustic metamaterials, Physical Review B 99, 020301 (2019).
  * Long  _et al._ [2020a] Y. Long, H. Ge, D. Zhang, X. Xu, J. Ren, M.-H. Lu, M. Bao, H. Chen, and Y.-F. Chen, Symmetry selective directionality in near-field acoustics, National Science Review 7, 1024 (2020a).
  * Long  _et al._ [2020b] Y. Long, D. Zhang, C. Yang, J. Ge, H. Chen, and J. Ren, Realization of acoustic spin transport in metasurface waveguides, Nature communications 11, 4716 (2020b).
  * Rayleigh [1885] L. Rayleigh, On waves propagated along the plane surface of an elastic solid, Proceedings of the London mathematical Society 1, 4 (1885).
  * Telford  _et al._ [1990] W. M. Telford, W. Telford, L. Geldart, and R. E. Sheriff,  _Applied geophysics_ (Cambridge university press, 1990).
  * Süsstrunk and Huber [2015] R. Süsstrunk and S. D. Huber, Observation of phononic helical edge states in a mechanical topological insulator, Science 349, 47 (2015).
  * Mousavi  _et al._ [2015] S. H. Mousavi, A. B. Khanikaev, and Z. Wang, Topologically protected elastic waves in phononic metamaterials, Nature communications 6, 8682 (2015).
  * Miniaci  _et al._ [2018] M. Miniaci, R. Pal, B. Morvan, and M. Ruzzene, Experimental observation of topologically protected helical edge modes in patterned elastic plates, Physical Review X 8, 031074 (2018).
  * Lu  _et al._ [2016] J. Lu, C. Qiu, M. Ke, and Z. Liu, Valley vortex states in sonic crystals, Physical review letters 116, 093901 (2016).
  * Lu  _et al._ [2017] J. Lu, C. Qiu, L. Ye, X. Fan, M. Ke, F. Zhang, and Z. Liu, Observation of topological valley transport of sound in sonic crystals, Nature Physics 13, 369 (2017).
  * Pal and Ruzzene [2017] R. K. Pal and M. Ruzzene, Edge waves in plates with resonators: an elastic analogue of the quantum valley hall effect, New Journal of Physics 19, 025001 (2017).
  * Vila  _et al._ [2017] J. Vila, R. K. Pal, and M. Ruzzene, Observation of topological valley modes in an elastic hexagonal lattice, Physical Review B 96, 134307 (2017).
  * Zhu  _et al._ [2018] H. Zhu, T.-W. Liu, and F. Semperlotti, Design and experimental observation of valley-hall edge states in diatomic-graphene-like elastic waveguides, Physical Review B 97, 174301 (2018).
  * Liu and Semperlotti [2018] T.-W. Liu and F. Semperlotti, Tunable acoustic valley–hall edge states in reconfigurable phononic elastic waveguides, Physical Review Applied 9, 014001 (2018).
  * Wu  _et al._ [2018] Y. Wu, R. Chaunsali, H. Yasuda, K. Yu, and J. Yang, Dial-in topological metamaterials based on bistable stewart platform, Scientific reports 8, 112 (2018).
  * Chaunsali  _et al._ [2018] R. Chaunsali, C.-W. Chen, and J. Yang, Subwavelength and directional control of flexural waves in zone-folding induced topological plates, Physical Review B 97, 054307 (2018).
  * Kane and Lubensky [2014] C. L. Kane and T. C. Lubensky, Topological boundary modes in isostatic lattices, Nature Physics 10, 39 (2014).
  * Lubensky  _et al._ [2015] T. Lubensky, C. Kane, X. Mao, A. Souslov, and K. Sun, Phonons and elasticity in critically coordinated lattices, Reports on Progress in Physics 78, 073901 (2015).
  * Mao and Lubensky [2018] X. Mao and T. C. Lubensky, Maxwell lattices and topological mechanics, Annual Review of Condensed Matter Physics 9, 413 (2018).
  * Sun and Mao [2020] K. Sun and X. Mao, Continuum theory for topological edge soft modes, Phys. Rev. Lett. 124, 207601 (2020).
  * Francesco  _et al._ [2012] P. Francesco, P. Mathieu, and D. Sénéchal,  _Conformal field theory_ (Springer Science & Business Media, 2012).
  * Lim and Lim [2015] T.-C. Lim and T.-C. Lim, Wave propagation in auxetic solids, Auxetic Materials and Structures , 367 (2015).
  * Sun  _et al._ [2012] K. Sun, A. Souslov, and X. Mao, Surface phonons, elastic response, and conformal invariance in twisted kagome lattices, Proceedings of the National Academy of Sciences 109, 12369 (2012).
  * Czajkowski  _et al._ [2022] M. Czajkowski, C. Coulais, M. van Hecke, and D. Z. Rocklin, Conformal elasticity of mechanism-based metamaterials, Nature Communications 13, 211 (2022).
  * Griffiths [2017] D. Griffiths,  _Introduction to Quantum Mechanics_ (Cambridge University Press, 2017).
  * Fruchart  _et al._ [2020] M. Fruchart, Y. Zhou, and V. Vitelli, Dualities and non-abelian mechanics, Nature 577, 636 (2020).
  * Lamb [1917] H. Lamb, On waves in an elastic plate, Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character 93, 114 (1917).
  * Saxena  _et al._ [2016] K. K. Saxena, R. Das, and E. P. Calius, Three decades of auxetics research- materials with negative poisson’s ratio: a review, Advanced Engineering Materials 18, 1847 (2016).
  * Yang  _et al._ [2015] L. Yang, O. Harrysson, H. West, and D. Cormier, Mechanical properties of 3d re-entrant honeycomb auxetic structures realized via additive manufacturing, International Journal of Solids and Structures 69, 475 (2015).
  * Li and Zeng [2016] Y. Li and C. Zeng, On the successful fabrication of auxetic polyurethane foams: Materials requirement, processing strategy and conversion mechanism, Polymer 87, 98 (2016).
  * Greaves  _et al._ [2011] G. N. Greaves, A. L. Greer, R. S. Lakes, and T. Rouxel, Poisson’s ratio and modern materials, Nature materials 10, 823 (2011).
  * Temnov [2012] V. V. Temnov, Ultrafast acousto-magneto-plasmonics, Nature Photonics 6, 728 (2012).
  * Guddala  _et al._ [2021] S. Guddala, F. Komissarenko, S. Kiriushechkina, A. Vakulenko, M. Li, V. Menon, A. Alù, and A. Khanikaev, Topological phonon-polariton funneling in midinfrared metasurfaces, Science 374, 225 (2021).
  * Choi  _et al._ [2022] W. J. Choi, K. Yano, M. Cha, F. M. Colombari, J.-Y. Kim, Y. Wang, S. H. Lee, K. Sun, J. M. Kruger, A. F. de Moura, _et al._ , Chiral phonons in microcrystals and nanofibrils of biomolecules, Nature Photonics 16, 366 (2022).
  * Holanda  _et al._ [2018] J. Holanda, D. Maior, A. Azevedo, and S. Rezende, Detecting the phonon spin in magnon–phonon conversion experiments, Nature Physics 14, 500 (2018).
  * He  _et al._ [2016] L. He, H. Li, and M. Li, Optomechanical measurement of photon spin angular momentum and optical torque in integrated photonic devices, Science advances 2, e1600485 (2016).
  * Landau  _et al._ [1986] L. Landau, E. Lifshitz, A. Kosevich, and L. Pitaevskii,  _Theory of Elasticity: Volume 7_, Course of theoretical physics (Elsevier Science, 1986).
  * Strauss [2007] W. Strauss,  _Partial Differential Equations: An Introduction_ (Wiley, 2007).
  * Gamelin [2003] T. Gamelin,  _Complex Analysis_, Undergraduate Texts in Mathematics (Springer New York, 2003).

Backscattering-free edge states below all bands in two-dimensional auxetic
media: Supplementary Materials

##  V 3D Rayleigh Waves

Here we review the linear elastic wave equation and Rayleigh wave formula in
3D. The main reference of this formulation is Ref. [56]. Consider a
homogeneous elastic material in $x^{3}\leq 0$ region in the
$x^{1}$-$x^{2}$-$x^{3}$ space with bulk modulus $B_{3D}$, Lamé constant
$\mu_{3D}$ and $\lambda_{3D}=B_{3D}-\frac{2}{3}\mu_{3D}$. Let $\vec{v}$ be the
displacement vector field. The strain tensor is given by

$u_{jk}=\frac{1}{2}(\delta_{ik}\frac{\partial v^{i}}{\partial
x^{j}}+\delta_{ij}\frac{\partial v^{i}}{\partial x^{k}}).$ (S1)

The elastic energy density for homogeneous material is

$\epsilon=\frac{1}{2}\lambda_{3D}\delta^{ij}\delta^{kl}u_{ij}u_{kl}+\mu_{3D}\delta^{ij}\delta^{kl}u_{ik}u_{jl}$
(S2)

The stress tensor $\sigma^{ij}$ is defined as

$d\epsilon=\sigma^{ij}du_{ij},$ (S3)

From this we obtain the stress tensor

$\sigma^{ij}=\lambda_{3D}\delta^{ij}\frac{\partial v^{a}}{\partial
x^{a}}+\mu_{3D}\delta^{ib}\frac{\partial v^{j}}{\partial
x^{b}}+\mu_{3D}\delta^{jc}\frac{\partial v^{i}}{\partial x^{c}}.$ (S4)

Let $\rho$ be the mass density of the elastic material. The equation of motion
is followed by Newton’s second law

$\rho\frac{\partial^{2}v^{i}}{\partial
t^{2}}=\frac{\partial\sigma^{ij}}{\partial x^{j}},$ (S5)

which is equivalent to

$\rho\frac{\partial^{2}\vec{v}}{\partial
t^{2}}=(\lambda_{3D}+\mu_{3D})\vec{\nabla}\vec{\nabla}\cdot\vec{v}+\mu_{3D}\nabla^{2}\vec{v}.$
(S6)

From Helmholtz decomposition, we can write

$\vec{v}=\vec{v_{l}}+\vec{v_{t}},$ (S7)

where $\vec{v_{l}}$ is the curl free and $\vec{v_{t}}$ is divergence free. By
taking curl on both sides of (S6), we get

$\vec{\nabla}\times(\rho\frac{\partial^{2}\vec{v_{t}}}{\partial
t^{2}}-\mu_{3D}\nabla^{2}\vec{v_{t}})=0.$ (S8)

Since $\vec{v_{t}}$ is divergence free, the vector field
$\rho\frac{\partial^{2}\vec{v_{t}}}{\partial
t^{2}}-\mu_{3D}\nabla^{2}\vec{v_{t}}$ is both gradient free and curl free,
hence

$\rho\frac{\partial^{2}\vec{v_{t}}}{\partial
t^{2}}-\mu_{3D}\nabla^{2}\vec{v_{t}}=0.$ (S9)

Similarly, by taking gradient on (S6), we have

$\rho\frac{\partial^{2}\vec{v_{l}}}{\partial
t^{2}}-(2\mu_{3D}+\lambda_{3D})\nabla^{2}\vec{v_{l}}=0.$ (S10)

Without loss of generality, we may assume the surface wave propagates in the
$x^{1}$ direction. Let the longitudinal and transverse components of the
Rayleigh wave be

$\displaystyle\vec{v_{l}}=e^{i(kx^{1}-\omega
t)+\kappa_{l}x^{3}}(\alpha^{1},\alpha^{2},\alpha^{3})^{T},$ (S11)
$\displaystyle\vec{v_{t}}=e^{i(kx^{1}-\omega
t)+\kappa_{t}x^{3}}(\beta^{1},\beta^{2},\beta^{3})^{T}.$

Plugging (S11) into (S9), (S10), we have

$\kappa_{l}=\sqrt{k^{2}-\frac{\omega^{2}}{c_{l}^{2}}},\;\kappa_{t}=\sqrt{k^{2}-\frac{\omega^{2}}{c_{t}^{2}}},$
(S12)

where $c_{l}=\sqrt{\frac{\lambda_{3D}+2\mu_{3D}}{\rho}}$ and
$c_{t}=\sqrt{\frac{\mu_{3D}}{\rho}}$ are the speed of longitudinal wave and
transverse wave. $x^{3}=0$ surface being an open boundary gives the boundary
conditions

$\displaystyle\sigma^{13}=\mu_{3D}(\frac{\partial v^{3}}{\partial
x^{1}}+\frac{\partial v^{1}}{\partial x^{3}})=0,$ (S13)
$\displaystyle\sigma^{23}=\mu_{3D}(\frac{\partial v^{3}}{\partial
x^{2}}+\frac{\partial v^{2}}{\partial x^{3}})=0,$
$\displaystyle\sigma^{33}=\lambda_{3D}\frac{\partial v^{1}}{\partial
x^{1}}+\lambda_{3D}\frac{\partial v^{2}}{\partial
x^{2}}+(\lambda_{3D}+2\mu_{3D})\frac{\partial v^{3}}{\partial x^{3}}=0.$

The conditions $\vec{\nabla}\times\vec{v_{l}}=0$,
$\vec{\nabla}\cdot\vec{v_{t}}=0$ and $\sigma^{23}=0$ give

$\displaystyle-$ $\displaystyle\kappa_{l}\alpha^{1}+ik\alpha^{3}=0,$ (S14)
$\displaystyle ik\beta^{1}+\kappa_{t}\beta^{3}=0$
$\displaystyle\alpha^{2}=\beta^{2}=0$

which says there exists $a$ and $b$ such that

$\displaystyle\alpha^{1}=ka,\alpha^{3}=-i\kappa_{l}a,$ (S15)
$\displaystyle\beta^{1}=\kappa_{t}b,\beta^{3}=-ikb.$

Combining (S11), (S13) and (S15) we get

$\displaystyle b(k^{2}+\kappa_{t}^{2})+2ak\kappa_{l}=0,$ (S16) $\displaystyle
2b\kappa_{t}k+a(k^{2}+\kappa_{t}^{2})=0.$

In order that (S16) has a non-zero solution, the coefficient matrix need to be
of determinant zero, which gives

$(k^{2}+\kappa_{t}^{2})^{2}=4k^{2}\kappa_{t}\kappa_{l}.$ (S17)

Let $\omega=c_{t}k\xi$, we obtain

$\xi^{6}-8\xi^{4}+8\xi^{2}(3-\frac{6}{3h_{3D}+4})-16(1-\frac{3}{3h_{3D}+4})=0,$
(S18)

where $h_{3D}=\frac{B_{3D}}{\mu_{3D}}$. The corresponding modes at $x^{3}=0$
surface are

$\displaystyle\psi_{k}(x^{1},x^{2},0)=e^{i(kx^{1}-\omega
t)}(-\frac{\xi^{2}\sqrt{1-\xi^{2}}}{2-\xi^{2}},0,i(\frac{2\sqrt{1-\xi^{2}}}{2-\xi^{2}}\sqrt{1-\frac{3\xi^{2}}{3h_{3D}+4}}-1))^{T},$
(S19) $\displaystyle\psi_{-k}(x^{1},x^{2},0)=e^{i(-kx^{1}-\omega
t)}(-\frac{\xi^{2}\sqrt{1-\xi^{2}}}{2-\xi^{2}},0,-i(\frac{2\sqrt{1-\xi^{2}}}{2-\xi^{2}}\sqrt{1-\frac{3\xi^{2}}{3h_{3D}+4}}-1))^{T}.$

##  VI 2D Rayleigh Waves

Here we review the elastic wave equation and Rayleigh wave formula for a flat
thin homogeneous elastic sheet within a plane stress condition, with bulk
modulus $B_{3D}$ and Lamé constant $\mu_{3D}$ and
$\lambda_{3D}=B_{3D}-\frac{2}{3}\mu_{3D}$. We choose a coordinate system such
that the $x^{1}$ and $x^{2}$ directions lie in the plane of the sheet and
$x^{3}$ direction is normal to the plane. On the top and bottom surfaces of
the sheet, we have

$\sigma^{13}=\sigma^{23}=\sigma^{33}=0.$ (S20)

By the continuity of stress tensor and assuming the thickness is small, we can
assume $\sigma^{13}=\sigma^{23}=\sigma^{33}=0$ holds everywhere. Making use of
(S20) and (S2), the energy density in this situation is

$\epsilon=\frac{1}{2}\lambda\delta^{ij}\delta^{kl}u_{ij}u_{kl}+\mu\delta^{ij}\delta^{kl}u_{ik}u_{jl},$
(S21)

where $i,j$ ranges from 1 to 2 and the effective 2D Lamé constants are

$\lambda=\frac{2(B_{3D}-\frac{2}{3}\mu_{3D})\mu_{3D}}{B_{3D}+\frac{4}{3}\mu_{3D}},~{}\mu=\mu_{3D}.$
(S22)

Thus, the effective 2D bulk modulus is

$B=\lambda+\mu=\frac{3B_{3D}\mu_{3D}}{B_{3D}+\frac{4}{3}\mu_{3D}}.$ (S23)

Let $\vec{v}=v^{1}\vec{x}^{1}+v^{2}\vec{x}^{2}$ be the in-plane displacement,
noting (S21) has the same form as (S2), the equation of motion for $\vec{v}$
is then

$\rho\frac{\partial^{2}\vec{v}}{\partial
t^{2}}=(\lambda+\mu)\vec{\nabla}\vec{\nabla}\cdot\vec{v}+\mu\nabla^{2}\vec{v}.$
(S24)

We simplify the equation again using the Helmholtz decomposition, giving

$\displaystyle\rho\frac{\partial^{2}\vec{v_{t}}}{\partial
t^{2}}-\mu\nabla^{2}\vec{v_{t}}=0,$ (S25)
$\displaystyle\rho\frac{\partial^{2}\vec{v_{l}}}{\partial
t^{2}}-(2\mu+\lambda)\nabla^{2}\vec{v_{l}}=0.$

Now we calculate Rayleigh waves in the aforementioned thin sheet, but modified
into a half-sheet laying in $x^{2}\leq 0$. Let the longitudinal and transverse
components of the Rayleigh wave be

$\displaystyle\vec{v_{l}}=e^{i(kx^{1}-\omega
t)+\kappa_{l}x^{2}}(\alpha^{1},\alpha^{2})^{T},$ (S26)
$\displaystyle\vec{v_{t}}=e^{i(kx^{1}-\omega
t)+\kappa_{t}x^{2}}(\beta^{1},\beta^{2})^{T}.$

Plugging (S26) into (S25), we have

$\kappa_{l}=\sqrt{k^{2}-\frac{\omega^{2}}{c_{l}^{2}}},\;\kappa_{t}=\sqrt{k^{2}-\frac{\omega^{2}}{c_{t}^{2}}},$
(S27)

where $c_{l}=\sqrt{\frac{\lambda+2\mu}{\rho}}$ and
$c_{t}=\sqrt{\frac{\mu}{\rho}}$ are the speeds of longitudinal wave and
transverse wave, respectively. The $x^{2}=0$ line being an open boundary gives
us the boundary conditions

$\displaystyle\sigma^{12}=\mu(\frac{\partial v^{1}}{\partial
x^{2}}+\frac{\partial v^{2}}{\partial x^{1}})=0,$ (S28)
$\displaystyle\sigma^{22}=\lambda\frac{\partial v^{1}}{\partial
x^{1}}+(\lambda+2\mu)\frac{\partial v^{2}}{\partial x^{2}}=0.$

Plugging (S26) into (S28), we get

$\displaystyle\kappa_{l}\alpha^{1}+ik\alpha^{2}+\kappa_{t}\beta^{1}+ik\beta^{2}=0,$
(S29)
$\displaystyle\frac{ik\lambda}{\lambda+2\mu}\alpha^{1}+\kappa_{l}\alpha^{2}+\frac{ik\lambda}{\lambda+2\mu}\beta^{1}+\kappa_{t}\beta^{2}=0.$

The curl free condition of $\vec{v_{l}}$ and divergence free condition of
$\vec{v_{t}}$ give us

$\displaystyle-$ $\displaystyle\kappa_{l}\alpha^{1}+ik\alpha^{2}=0,$ (S30)
$\displaystyle ik\beta^{1}+\kappa_{t}\beta^{2}=0.$

We have 4 linear equations for 4 variables
$\alpha^{1},\alpha^{2},\beta^{1},\beta^{2}$. In order that the equations have
non-zero solutions, the corresponding coefficient matrix needs to be of
determinant zero, which gives us

$\xi^{6}-8\xi^{4}+(24-\frac{16c_{t}^{2}}{c_{l}^{2}})\xi^{2}+\frac{16c_{t}^{2}}{c_{l}^{2}}-16=0,\;\xi=\frac{\omega}{kc_{t}}.$
(S31)

In terms of $h={B}/{\mu}$, $\xi$ satisfies

$\xi^{6}-8\xi^{4}+(8+\frac{16h}{1+h})\xi^{2}-\frac{16h}{1+h}=0.$ (S32)

The corresponding modes at the $x^{2}=0$ line are

$\displaystyle\psi_{k}(x^{1},0)=e^{i(kx^{1}-\omega
t)}(\frac{1}{2}\xi\sqrt{1+\frac{1}{h}},\frac{1}{4}i\xi(2-\xi^{2})\sqrt{\frac{1+h}{h-h\xi^{2}}})^{T},$
(S33) $\displaystyle\psi_{-k}(x^{1},0)=e^{i(-kx^{1}-\omega
t)}(\frac{1}{2}\xi\sqrt{1+\frac{1}{h}},-\frac{1}{4}i\xi(2-\xi^{2})\sqrt{\frac{1+h}{h-h\xi^{2}}})^{T}.$

##  VII Linear response to point source on the edge

Here we derive the linear response theory for both 2D continuum and discrete
systems (the same derivation applies to 3D systems) for systems subject to a
point source of vibration on the edge.

###  VII.1 Continuum version

Consider the same homogeneous elastic thin sheet as described in Sec. VI, we
apply a local driving force $\vec{F}(r,t)=Re(\vec{F}e^{-i\omega t})\delta(r)$
to the boundary (the $(0,0)$ point) of the system. The complex equation of
motion is then

$\displaystyle\rho\frac{\partial^{2}\vec{v}}{\partial
t^{2}}=(\lambda+\mu)\vec{\nabla}\vec{\nabla}\cdot\vec{v}+\mu\nabla^{2}\vec{v}+\vec{F}\delta(r)e^{-i\omega
t},$ (S34) $\displaystyle\vec{v}(r,0)=\frac{\partial\vec{v}}{\partial
t}(r,0)=0.$

Let $\vec{\psi}_{k}(r)e^{-i\omega_{k}t}$ be Rayleigh waves with wave vector
$k$ and $\vec{\phi}_{s}(r,t)$ be bulk waves with index $s$ (since the bulk
waves are not important here we use an arbitrary index to identify them). The
solution of (S34) can then be decomposed into

$\vec{v}(r,t)=\int_{k}dkC_{k}(t)\vec{\psi}_{k}(r)e^{-i\omega_{k}t}+\int_{s}dsG_{s}(t)\vec{\phi}_{s}(r,t).$
(S35)

Plugging (S35) into (S34), applying $\vec{\psi}_{k}(r)e^{-i\omega_{k}t}$ to
both sides, and doing an integration over time and position, we have

$\displaystyle\frac{d^{2}C_{k}(t)}{dt^{2}}-$ $\displaystyle
2i\omega_{j}(k)\frac{dC_{k}(t)}{dt}=\vec{\psi}_{k}(0)^{*}.\vec{F}\;e^{-i(\omega-\omega_{j}(k))t},$
(S36) $\displaystyle\left.\frac{dC_{k}(t)}{dt}\right|_{t=0}=C_{k}(0)=0.$

For a given $k$, if $\omega\neq\omega_{k}$, $C_{k}(t)$ is bounded. However, if
$\omega=\omega_{k}$, we have

$C_{j}(k,t)=\frac{i\vec{\psi}_{k}^{*}(0).\vec{F}}{\omega}t+q_{k}(t),$ (S37)

where $q_{k}(t)$ is bounded. Therefore, after long times, most energy is in
the two modes $\vec{\psi}_{k}(r)e^{-i\omega_{k}t}$ and
$\vec{\psi}_{-k}(r)e^{-i\omega_{k}t}$. To maximally excite the $+k$ Rayleigh
wave, the force applied to system needs to follow the trajectory of the motion
corresponding to $\vec{\psi}_{k}(r)e^{-i\omega_{k}t}$, which means
$\vec{\psi}_{k}(0)=\vec{F}$. In this case, the ratio $|{C_{-k}}/{C_{k}}|$,
whose square quantifies the ratio between energy going in $+k$ direction and
$-k$ direction, is then

$\Big{|}\frac{C_{-k}}{C_{k}}\Big{|}=\Big{|}\frac{\vec{\psi}_{-k}(0)^{*}.\vec{\psi}_{k}(0)}{\vec{\psi}_{k}(0)^{*}.\vec{\psi}_{k}(0)}\Big{|}.$
(S38)

Fig. 1(c) and Fig. 1(d) in the main text show $|{C_{-k}}/{C_{k}}|$ as a
function of ${B_{3D}}/{\mu_{3D}}$.

###  VII.2 Discrete Version

Consider a spring mass network that is periodic, and infinite in the $x^{2}$
direction and finite in the $x^{1}$ direction. For simplicity, we assume the
lattice spacing in the $x^{2}$ direction to be $1$. Such system has two edge
bands, and a finite number of bulk bands. We index the bands from low to high
frequency by $1,2,3...2N$ (assuming N nodes per unit cell) and also index the
nodes by $1,2,3...$ (The entire system has infinitely many nodes). Let
$\vec{u}^{m}$ be the displacement vector of the $m^{th}$ node and
$\vec{u}=(\vec{u}^{1},...,\vec{u}^{n}...)$ be the displacements of the entire
system. We also use $\vec{F}^{i}$ to denote the force on the $i^{th}$ node and
$\vec{F}=(\vec{F}^{1},...,\vec{F}^{n}...)$ be the forces on the entire system.
We assume an external force $\vec{F}^{i}=(F^{i}_{1}\cos\omega
t,F^{i}_{2}\cos(\omega t+\phi))$, where the frequency $\omega$ lies in the
range of the edge bands, is applied to the $i^{th}$ node at the time $t=0$.
Then, the equation of motion of the system together with initial conditions is

$m\vec{u}^{\prime\prime}=-D\vec{u}+\vec{F},\;\vec{u}(0)=\vec{0},\;\vec{u}^{\prime}(0)=\vec{0}.$
(S39)

To solve this equation, we rewrite $\vec{F}=Re(\vec{f}e^{-i\omega t}$) and
solve the complex differential equation

$m\vec{u}^{\prime\prime}=-D\vec{u}+\vec{f}e^{-i\omega
t},\;\vec{u}(0)=\vec{0},\;\vec{u}^{\prime}(0)=\vec{0}.$ (S40)

The real part of the solution of (S40) is the solution of (S39). Let
$\vec{\psi}_{j}(k)$ be the normalized Bloch state (where their inner product
is the $\delta$ function) in the $j^{th}$ band with wave vector $k$ and
corresponding eigenfrequency $\omega_{j}(k)$. The solution of the equation of
motion can be expanded as

$\vec{u}(t)=\int_{-\pi}^{\pi}C_{j}(k,t)e^{-i\omega_{j}(k)t}\vec{\psi}_{j}(k).$
(S41)

Plugging this expansion into (S39), we have

$\displaystyle\frac{d^{2}C_{j}(k,t)}{dt^{2}}-$ $\displaystyle
2i\omega_{j}(k)\frac{dC_{j}(k,t)}{dt}=\langle\vec{\psi}_{j}(k)|\vec{f}\rangle
e^{-i(\omega-\omega_{j}(k))t}$ (S42)
$\displaystyle\left.\frac{dC_{j}(k,t)}{dt}\right|_{t=0}=C_{j}(k,0)=0.$

For a given $j,k$, if $\omega\neq\omega_{j}(k)$, $C_{j}(k,t)$ is bounded.
However, if $\omega=\omega_{j}(k)$, we have

$C_{j}(k,t)=\frac{i\langle\vec{\psi}_{j}(k)|\vec{f}\rangle}{\omega}t+q_{j}(k,t),$
(S43)

where $q_{j}(k,t)$ is bounded. Therefore, at long times, most energy is in the
four modes $\vec{\psi}_{1}(k_{1})$, $\vec{\psi}_{1}(-k_{1})$,
$\vec{\psi}_{2}(k_{2})$, and $\vec{\psi}_{2}(-k_{2})$. Without loss of
generality, we can assume the $i^{th}$ node is on the left edge of the system
and we want most energy localized at the edge propagating in $+x^{2}$
direction. The force applied to the $i^{th}$ node needs to follow the
trajectory of the motion corresponding to some linear combination of modes
$\vec{\psi}_{1}(k_{1})$ and $\vec{\psi}_{2}(k_{2})$ at node $i$.
Mathematically, that means $\vec{f}$ maximizing

$\frac{|\langle\vec{\psi}_{1}(k_{1})|\vec{f}\rangle|^{2}+|\langle\vec{\psi}_{2}(k_{2})|\vec{f}\rangle|^{2}}{|\langle\vec{f}|\vec{f}\rangle|^{2}}=\frac{|\langle\vec{\psi}_{1}(k_{1})^{i}|\vec{f}^{i}\rangle|^{2}+|\langle\vec{\psi}_{2}(k_{2})^{i}|\vec{f}^{i}\rangle|^{2}}{|\langle\vec{f}^{i}|\vec{f}^{i}\rangle|^{2}}.$
(S44)

##  VIII Orbital and spin angular momentum of phonons

In this section, we discuss the general formulation of angular momentum of
phonons as vector fields and their orbital and spin components. The main
reference of this formulation is Ref. [40].

###  VIII.1 Noether’s theorem

Let’s start with a brief review of Noether’s theorem in field theory, where
each continuous symmetry leads to a conserved current.

Consider a generic action $S$ of a set of fields $\phi_{a}$, where
$a=1,\ldots,n$, such that

$S=\int dx\,\mathcal{L}[\phi_{a},\partial_{\mu}\phi_{a}],$ (S45)

where we have kept in the Lagrangian to leading order derivatives of $\phi$.
To make the notation compact, we use the 4D space-time coordinate
$x^{\mu}\equiv(t,x^{1},x^{2},x^{3})$, where the first component is time and
the latter three are coordinates in 3D space. A generic infinitesimal
transformation can be written as

$\displaystyle x^{\mu}$ $\displaystyle\to x^{\prime\mu}=x^{\mu}+\Delta
x^{\mu},$ $\displaystyle\phi_{a}(x)$
$\displaystyle\to\phi^{\prime}_{a}(x^{\prime})=\phi_{a}(x)+\Delta\phi_{a}(x),$
(S46)

where $\Delta x$ and $\Delta\phi$ are both caused by the same transformation
$w$, which we can generally write as

$\displaystyle\Delta x^{\mu}=\frac{\partial x^{\mu}}{\partial w_{b}}w_{b},$
$\displaystyle\Delta\phi_{a}=\frac{\partial\phi_{a}}{\partial w_{b}}w_{b}.$
(S47)

Examples of such transformations can be found in the following subsections.

The transformation of the spatial derivatives of the field follows

$\partial_{\mu}\phi_{a}\to\partial_{\mu}^{\prime}\phi^{\prime}_{a}=\partial_{\mu}\phi_{a}-\partial_{\nu}\phi_{a}\partial_{\mu}\Delta
x^{\nu}+\partial_{\mu}\Delta\phi_{a},$ (S48)

where the second term on the right hand side is from the variation of $x$, and
the last term is from the variation of $\phi$, as results of the
transformation.

At the same time, the integral over $dx^{\prime}$ can be written as the
integral over $dr$ with a Jacobian

$dx\to
dx^{\prime}=dx\det(\partial_{\mu}x^{\prime\nu})=dx(1+\partial_{\mu}\Delta
x^{\mu}).$ (S49)

Putting them all together, we have the transformed action

$S^{\prime}=\int dx(1+\partial_{\mu}\Delta
x^{\mu})\mathcal{L}\left[\phi_{a}(x)+\Delta\phi_{a}(x),\partial_{\mu}\phi_{a}-\partial_{\nu}\phi_{a}\partial_{\mu}\Delta
x^{\nu}+\partial_{\mu}\Delta\phi_{a}\right].$ (S50)

This action can be expanded to leading order in the transformation $\Delta
x,\Delta\phi$, so that

$S^{\prime}=\int dx\left[\mathcal{L}+\partial_{\mu}\Delta
x^{\mu}\mathcal{L}+\frac{\partial\mathcal{L}}{\partial\phi_{a}}\Delta\phi_{a}(x)-\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\partial_{\nu}\phi_{a}\partial_{\mu}\Delta
x^{\nu}+\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\partial_{\mu}\Delta\phi_{a}(x)\right].$
(S51)

We can also use the Euler-Lagragian equation (the equation of motion of this
field theory)

$\frac{\partial\mathcal{L}}{\partial\phi_{a}}=\partial_{\mu}\left(\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\right)$
(S52)

to simplify $S^{\prime}$ so that

$S^{\prime}-S=\int
dx\left\\{\left[\delta^{\mu}_{\nu}\mathcal{L}-\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\partial_{\nu}\phi_{a}\right]\partial_{\mu}\Delta
x^{\nu}+\partial_{\mu}\left[\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\Delta\phi_{a}(x)\right]\right\\},$
(S53)

where the first bracket denotes the effect of (external) transformation of
coordinate, and the second bracket denotes the effect of (internal)
transformation of the field.

Because both of these transformations are caused by the same source $w(x)$, we
can use Eq. VIII.1 to rewrite everything in terms of $w$. In particular, we
have terms carrying $w$ and terms carrying $\partial_{\mu}w$.

Since Eq. VIII.1 can be viewed as simply a coordinate transformation
(homogeneous or inhomogeneous), we must have $S^{\prime}=S$ by definition. On
the other hand, when the system has a continuous symmetry such as a
homogeneous spatial translation or rotation, the Lagragian (not just the total
action) must remain the same before and after the transformation. Thus the $w$
term must identically vanish everywhere. Therefore, to make $S^{\prime}=S$, we
must also have the integral of the $\partial_{\mu}w$ terms vanish, which are

$0=\int
dx\left\\{\left[\delta^{\mu}_{\nu}\mathcal{L}-\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\partial_{\nu}\phi_{a}\right]\frac{\partial\Delta
x^{\nu}}{\partial
w_{b}}+\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\frac{\partial\Delta\phi_{a}}{\partial
w_{b}}\right\\}\partial_{\mu}w_{b}.$ (S54)

Because $w$ can be an arbitrary transformation, we move the partial
differential via integral by parts and find conserved currents

$j^{\mu}_{b}=\left[-\delta^{\mu}_{\nu}\mathcal{L}+\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\partial_{\nu}\phi_{a}\right]\frac{\partial\Delta
x^{\nu}}{\partial
w_{b}}-\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_{a})}\frac{\partial\Delta\phi_{a}}{\partial
w_{b}},$ (S55)

which satisfies $\partial_{\mu}j^{\mu}_{b}=0$.

Now we have proved Noether’s theorem, where each continuous symmetry $w_{b}$
leads to a conserved current $j^{\mu}_{b}$ satisfying
$\partial_{\mu}j^{\mu}_{b}=0$. Next we discuss its consequence on spatial
translations and rotations, the cases of interest in this paper.

###  VIII.2 Translational symmetry

The case of spatial translations can be simply written as translations in the
three spatial directions $w=(w_{1},w_{2},w_{3})$. Thus $\frac{\partial\Delta
x^{\nu}}{\partial w_{b}}=0$ for $\nu=0$ (time), and $\delta_{b}^{\nu}$ for
$\nu=1,2,3$ (space). The $\frac{\partial\phi_{a}}{\partial w_{b}}$ terms all
vanish as the field itself doesn’t transform under spatial translation.

Plugging these into the formula of conserved current [Eq. (S55)], we find

$\displaystyle j_{b}^{0}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial(\partial_{0}\phi_{a})}\partial_{b}\phi_{a},$
(S56) $\displaystyle j_{b}^{i}$
$\displaystyle=-\delta_{b}^{i}\mathcal{L}+\frac{\partial\mathcal{L}}{\partial(\partial_{i}\phi_{a})}\partial_{b}\phi_{a},$
(S57)

where index $i=1,2,3$ runs through spatial indices.

Now we restore the notation where time and space are separated (instead of the
4D notation), and since all indices are in Euclidean space, we will have them
all as sub-indices. We have the resulting continuity equation

$\partial_{t}g_{m}+\partial_{n}\tau_{mn}=0,$ (S58)

which describes _momentum conservation_ of this field theory. Here the
momentum density $g$ and momentum flow $\tau$ are given by

$\displaystyle g_{m}\equiv-j_{m}^{0},$ (S59) $\displaystyle\tau_{mn}\equiv-
j_{m}^{n}.$ (S60)

This continuity equation [Eq. (S58)] tells us that as a result of spatial
translational symmetry, the total momentum of the system

$\displaystyle P_{m}=\int d^{3}xg_{m}(x)$ (S61)

is a conserved quantity.

###  VIII.3 Rotational symmetry

The case of spatial rotations, which is of central interest to our discussion
here, is slightly more complicated. An infinitesimal spatial rotation can be
written as $x\to x^{\prime}=(E+I\cdot w)x$, where $E$ is identity matrix, $I$
denotes the 3 generators of $SO(3)$, and $w=(w_{1},w_{2},w_{3})$ are the
angles of rotation in 3D along $x,y,z$, respectively. This can be written in a
more explicit form

$x^{\prime}_{i}=x_{i}+\varepsilon_{ibj}\,r_{j}\,w_{b},$ (S62)

where $\varepsilon$ is the Levi-Civita symbol. Thus, we have
$\frac{\partial\Delta x_{i}}{\partial w_{b}}=\varepsilon_{ibj}\,x_{j}$ for the
spatial degrees of freedom and $\Delta x^{0}=0$.

How does the field $\phi$ transform under spatial rotation? This depends on
the type of field we consider. Scalar fields won’t change. Vector fields, the
case of interest here as elasticity is concerned with displacements of
components in a solid, should rotate as a vector. This tells us that
$\frac{\partial\Delta\phi_{a}}{\partial w_{b}}=\varepsilon_{abc}\,\phi_{c}$.

Plugging these into the formula of conserved current [Eq. (S55)], we find that
for the case of rotation

$\displaystyle j_{b}^{0}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial(\partial_{0}\phi_{a})}\partial_{\nu}\phi_{a}\varepsilon_{\nu
bm}\,x_{m}-\frac{\partial\mathcal{L}}{\partial(\partial_{0}\phi_{a})}\varepsilon_{abm}\,\phi_{m},$
(S63) $\displaystyle j_{b}^{i}$
$\displaystyle=\left[-\delta_{\nu}^{i}\mathcal{L}+\frac{\partial\mathcal{L}}{\partial(\partial_{i}\phi_{a})}\partial_{\nu}\phi_{a}\right]\varepsilon_{\nu
bm}\,x_{m}-\frac{\partial\mathcal{L}}{\partial(\partial_{i}\phi_{a})}\varepsilon_{abm}\,\phi_{m},$
(S64)

where index $i=1,2,3$ runs through spatial indices (and we again wrote them
all as sub-indices for simplicity). From the continuity equation,
$\partial_{\mu}j^{\mu}=0$, we know that $j_{b}^{0}$ offers a new conserved
quantity. In particular, this can be written into two parts: the first term,
which is from the rotation of the coordinate, and the second term which is
from the rotation of the field,

$\displaystyle-j_{b}^{0}$ $\displaystyle=l_{b}+s_{b},$ (S65)

where

$\displaystyle l_{b}$ $\displaystyle\equiv\varepsilon_{bm\nu}x_{m}g_{\nu},$
(S66) $\displaystyle s_{b}$
$\displaystyle\equiv\varepsilon_{bma}\phi_{m}\frac{\partial\mathcal{L}}{\partial(\partial_{t}\phi_{a})}.$
(S67)

These two quantities are the orbital and spin angular momenta of the vector
field $\phi$. This can also be written as

$\displaystyle\vec{l}$ $\displaystyle\equiv\vec{x}\times\vec{g},$ (S68)
$\displaystyle\vec{s}$
$\displaystyle\equiv\vec{\phi}\times\frac{\partial\mathcal{L}}{\partial(\partial_{t}\vec{\phi})}.$
(S69)

###  VIII.4 Angular momentum of elastic waves

Now we apply this formalism to elasticity, to derive the expressions for the
orbital and spin angular momenta.

The Lagrangian of a general isotropic solid can be written as

$\mathcal{L}=\frac{\rho}{2}\partial_{t}u_{i}\partial_{t}u_{i}-\frac{\lambda}{2}(\partial_{i}u_{i})^{2}-\frac{\mu}{2}(\partial_{i}u_{j}\partial_{i}u_{j}+\partial_{i}u_{j}\partial_{j}u_{i})^{2},$
(S70)

where $\rho$ is the mass density, $\mu,\lambda$ are the Lame coefficients.
Identifying $\vec{u}$ as the vector field $\vec{\phi}$ discussed above, this
leads to densities of momentum $g$, orbital angular momentum $l$ and spin
angular momentum $s$, where

$\displaystyle g_{i}$ $\displaystyle=-\rho\partial_{t}u_{j}\partial_{i}u_{j},$
(S71) $\displaystyle l_{i}$ $\displaystyle=\varepsilon_{ijk}x_{j}g_{k},$ (S72)
$\displaystyle s_{i}$
$\displaystyle=\rho\varepsilon_{ijk}u_{j}\partial_{t}u_{k}.$ (S73)

In the vectorial form, the spin angular momentum can also be written as

$\displaystyle\vec{s}$ $\displaystyle=\rho\vec{u}\times(\partial_{t}\vec{u}).$
(S74)

From this formula it is clear that the spin angular momentum of these elastic
waves characterize the angular momentum of mass points rotating around their
mean positions. In contrast the orbital angular momentum comes from the wave
momentum $\vec{k}$.

We can apply them to a generic plane-wave field

$\vec{u}(\vec{x},t)=\real\left(\vec{A}\,e^{i(\vec{k}\cdot\vec{x}-\omega
t)}\right),$ (S75)

where the amplitude $A$ can be a complex vector denoting possible phase
difference of different wave components ($u_{x},u_{y},u_{z}$). This leads to

$\displaystyle\vec{g}(\vec{x},t)$
$\displaystyle=\rho\omega\left|\imaginary\left(\vec{A}e^{i(\vec{k}\cdot\vec{x}-\omega
t)}\right)\right|^{2}\vec{k}$ (S76) $\displaystyle\vec{l}(\vec{x},t)$
$\displaystyle=\rho\omega\left|\imaginary\left(\vec{A}e^{i(\vec{k}\cdot\vec{x}-\omega
t)}\right)\right|^{2}\vec{x}\times\vec{k},$ (S77)
$\displaystyle\vec{s}(\vec{x},t)$
$\displaystyle=\frac{1}{2}\rho\omega\imaginary(\vec{\overline{A}}\times\vec{A}),$
(S78)

where $\real$ and $\imaginary$ denote the real and imaginary parts,
respectively, and $\overline{A}$ is the complex conjugate of $A$. From this
form, it is clear that $\vec{\overline{A}}\times\vec{A}$ has to be complex in
order for this wave to carry a spin angular momentum.

###  VIII.5 Spin angular momentum of elastic waves in 2D and 3D

We first use these formula to examine the spin of elastic waves in the bulk of
isotropic materials. In 3D, it is obvious that longitudinal waves don’t carry
spin: their $\vec{A}$ is along the wave vector $\vec{k}$ so there is only one
direction of vibration, leading to $\vec{\overline{A}}\times\vec{A}=0$.

Transverse waves in 3D isotropic solids can be generally written as
$\vec{A}=\\{A_{1},A_{2},0\\}$, taking the direction of $k$ to be the third
dimension. Without losing generality we can normalize the amplitude to be
unity and factorize out a common phase, so that
$\vec{A}=\\{1,e^{i\theta},0\\}$. Therefore
$\vec{\overline{A}}\times\vec{A}=\\{0,0,\sin\theta\\}$, which is along
$\vec{k}$ with an amplitude between $[-1,1]$. Thus
$\vec{s}=\frac{1}{2}\rho\omega\sin\theta e_{\vec{k}}$ where $e_{\vec{k}}$ is
the unit vector along $\vec{k}$. It is straightforward to see that circularly
polarized transverse waves maximize $|\vec{s}|$.

In 2D, there are only 2 directions for vibrations, and the wave vector
$\vec{k}$ is in the same plane. The wave equation can be generally written as
(S25), where it is clear that, for materials with finite bulk modulus $B>0$
and shear modulus $\mu>0$, the longitudinal and transverse waves have
different dispersion relations and can’t combine to carry any angular momenta
in the bulk.

In the special case of $B=0$ in 2D, the two waves have the same wavespeeds
(and linear dispersion relations). They can be written in a similar form as
the 3D case, leading to

$\vec{s}=\frac{1}{2}\rho\omega\sin\theta e_{\vec{z}},$ (S79)

describing spin angular momentum along the normal direction of the 2D plane.

When quantizing this problem, the spin angular momentum of phonons as quanta
of these waves take values of $(-\hbar,0,+\hbar)$. In 3D, this is similar to
spins of photons, but the definition of $S_{z}$ in 2D is unique to $B=0$
materials.

Interestingly, near the edge of a 2D sheet, Rayleigh waves combine both
longitudinal and transverse components, and carry spin angular momentum, as we
discuss in the main text. It also becomes clear here that the same wavespeeds
of longitudinal and transverse waves in the $B=0$ case allows the spin of the
Rayleigh wave extend into the bulk (Fig. 1 in the main text).

##  IX Backscattering-free edge waves and winding number in 2D $B/\mu\ll 1$
continuum

In this section, we discuss in more detail the complex analytic function
formulation we used to describe backscattering-free edge waves in 2D $B/\mu\ll
1$ continuum, and the definition of the new winding number.

Consider a 2D homogeneous elastic material with effective bulk modulus $B=0$
and effective Lamé constants $\lambda=-\mu$ in a bounded simply connected
region $\Omega\subset\mathbb{C}$. The equation of motion together with stress
free boundary conditions under $B=0$ is

$\displaystyle\rho\frac{\partial^{2}\vec{v}}{\partial
t^{2}}=\mu\nabla^{2}\vec{v},$ (S80) $\displaystyle
n^{1}\sigma^{11}+n^{2}\sigma^{12}=0,$ $\displaystyle
n^{1}\sigma^{12}+n^{2}\sigma^{22}=0,$

where $\vec{n}=(n^{1},n^{2})$ is the vector normal to the boundary (pointing
outward). The in plane displacement vector field
$\vec{v}=v^{1}\vec{x}_{1}+v^{2}\vec{x}_{2}$ can be identified with a complex
scalar field $\psi$ via $\psi=v^{1}+iv^{2}$. We can rewrite the boundary
condition as

$n^{1}\sigma^{11}+n^{2}\sigma^{12}+i(n^{1}\sigma^{12}+n^{2}\sigma^{22})=0.$
(S81)

Under the identification above, the equation of motion together with stress
free boundary conditions under $B=0$ is

$\displaystyle\rho\frac{\partial^{2}\psi(z,t)}{\partial t^{2}}$
$\displaystyle=4\mu\frac{\partial^{2}\psi(z,t)}{\partial z\partial
z^{*}},z\in\Omega,$ (S82) $\displaystyle\frac{\partial\psi(z,t)}{\partial
z^{*}}$ $\displaystyle=0,z\in\partial\Omega,$

where $\frac{\partial}{\partial z}=\frac{1}{2}(\frac{\partial}{\partial
x}-i\frac{\partial}{\partial y})$, $\frac{\partial}{\partial
z^{*}}=\frac{1}{2}(\frac{\partial}{\partial x}+i\frac{\partial}{\partial y})$,
and $\partial\Omega$ is the boundary of $\Omega$. Note that (S82) admits
separation of variables

$\psi(z,t)=\phi(z)e^{i\omega t}$ (S83)

and

$\psi(z,t)=\phi(z)e^{-i\omega t},$ (S84)

where (S83) corresponds to spin-up states (rotate counterclockwise) and (S84)
corresponds to spin-down states (rotate clockwise). The profile $\phi(z)$
satisfies

$\displaystyle-\omega^{2}\rho\phi(z)$
$\displaystyle=4\mu\frac{\partial^{2}\phi}{\partial z\partial
z^{*}}(z),z\in\Omega,$ (S85) $\displaystyle\frac{\partial\phi}{\partial
z^{*}}(z)$ $\displaystyle=0,z\in\partial\Omega.$

To solve (S85), we let

$h(z)=\frac{\partial\phi}{\partial z^{*}}.$ (S86)

By taking $\frac{\partial}{\partial z^{*}}$ on both sides of (S85), we see
$h(z)$ satisfies Helmholtz equation with Dirichlet boundary conditions [57]

$\displaystyle-\frac{\omega^{2}\rho}{4\mu}h(z)$
$\displaystyle=\frac{\partial^{2}h}{\partial z\partial z^{*}}(z),z\in\Omega,$
(S87) $\displaystyle h(z)$ $\displaystyle=0,z\in\partial\Omega.$

Since $\Omega$ is a bounded region, it is follows that the eigenvalues
$-\frac{\omega^{2}\rho}{4\mu}$ are discrete. When $\omega\neq 0$, $h(z)$ has
large variation in the bulk and
$\phi(z)=-\frac{4\mu}{\omega^{2}\rho}\frac{\partial h}{\partial z}$ has large
value throughout the bulk, hence it is not an edge mode. When $\omega=0$,
$h(z)$ has unique solution $h(z)\equiv 0$, hence $\frac{\partial\phi}{\partial
z^{*}}(z)=0$, which is equivalent to $\phi(z)$ being a bounded analytic
function. Here, bounded means

$\exists R>0,|\phi(z)|<R,\forall z\in\Omega$ (S88)

and analytic means the limit

$\lim_{u\rightarrow 0}\frac{\phi(z+u)-\phi(z)}{u}$ (S89)

exists for all $z\in\Omega$. For $B>0$, the equation of motion together with
the boundary condition becomes

$\displaystyle\frac{\rho}{\mu}\frac{\partial^{2}\psi(z,t)}{\partial t^{2}}$
$\displaystyle=\frac{B}{\mu}(2\frac{\partial^{2}\psi(z,t)}{\partial z\partial
z^{*}}+\frac{\partial^{2}\psi^{*}(z,t)}{\partial
z^{2}})+4\frac{\partial^{2}\psi(z,t)}{\partial z\partial z^{*}},z\in\Omega,$
(S90) $\displaystyle\frac{\partial\psi(z,t)}{\partial z^{*}}$
$\displaystyle=\frac{B}{2\mu}\frac{n^{1}-in^{2}}{n^{1}+in^{2}}(\frac{\partial\psi(z,t)}{\partial
z}+\frac{\partial\psi^{*}(z,t)}{\partial z^{*}}),z\in\partial\Omega.$

In the $B\ll\mu$ limit, the $B/\mu$ term can be viewed as a perturbation. Let
$\psi_{\omega(B/\mu)}(z,t)$ be an eigenmode with frequency $\omega(B/\mu)$,
where

$\lim_{B/\mu\rightarrow 0}\omega(B/\mu)=0.$ (S91)

We have

$\lim_{B/\mu\rightarrow 0}\psi_{\omega(B/\mu)}(z,t)=\phi(z)$ (S92)

for some bounded analytic function $\phi(z)$. If we start from $B/\mu=0$, then
allow some small $B/\mu$, since $\phi(z)e^{i\omega(B/\mu)t}$ and
$\phi(z)e^{-i\omega(B/\mu)t}$ serve as zeroth order (in $B/\mu$) solution to
(S90), we see

$\psi_{\omega(B/\mu)}(z,t)\approx\phi(z)e^{i\omega(B/\mu)t}$ (S93)

or we have

$\psi_{\omega(B/\mu)}(z,t)\approx\phi(z)e^{-i\omega(B/\mu)t}.$ (S94)

Now we prove the following claim: for modes $\phi(z)$ giving rise to waves
$\psi_{\omega(B/\mu)}(z,t)$ localized at the boundary, we can find
$z_{0}\in\Omega$ such that $\phi(z_{0})=0$ and $\phi(z)\neq 0$ when
$z\in\partial\Omega$. The claim is a consequence of Rouché’s theorem [58],
which is:

######  Theorem IX.1.

Let $\Omega$ be a bounded domain with piecewise smooth boundary
$\partial\Omega$. Let $\phi(z)$ and $g(z)$ be two analytic functions that can
be analytically extended to an open set containing $\Omega\cup\partial\Omega$.
If $|g(z)|<|\phi(z)|,\forall z\in\partial\Omega$, then $\phi(z)$ and
$\phi(z)+g(z)$ have the same number of zeros in $\Omega$, counting
multiplicities.

Mode localized at the boundary has small displacement in the bulk and large
displacement at all boundaries, and hence we can find $z_{1}\in\Omega$ such
that

$0<|\phi(z_{1})|<|\phi(z)|,\forall z\in\partial\Omega.$ (S95)

We then define $g(z)\equiv-\phi(z_{1})$ and consider function

$f(z)=\phi(z)+g(z).$ (S96)

Since $\phi(z)$ and $g(z)$ satisfy (S95), we can apply Rouché’s theorem and
see if $f(z)$ and $\phi(z)$ have the same number of zeros (counting
multiplicity) in $\Omega$. Since $z=z_{1}$ is a zero of $f(z)$ in $\Omega$, we
know $\phi(z)$ has at least one zero in $\Omega$.

Now we look at the behavior of profile $\phi(z)$ at the boundary. By argument
principle [58],

$\frac{1}{2\pi i}\int_{\partial\Omega}d\ln(\phi(z))=\frac{1}{2\pi
i}\int_{\partial\Omega}\frac{\phi^{\prime}(z)}{\phi(z)}dz=N_{0}-N_{\infty},$
(S97)

where $\frac{1}{2\pi i}\int_{\partial\Omega}d\ln(\phi(z))$ is the number of
counterclockwise turns $\phi(z)$ rotates when walking counterclockwise along
the boundary (this means $\Omega$ is always on the left hand side), $N_{0}$ is
the number of zeros of $\phi(z)$ (counting multiplicity) inside $\Omega$, and
$N_{\infty}$ is the number of poles of $\phi(z)$ (counting multiplicity)
inside $\Omega$. We note that $\phi(z)$ being bounded in $\Omega$ implies
$N_{\infty}=0$ and $\phi(z_{0})=0$ implies $N_{0}\geq 1$. Hence when walking
counterclockwise along the boundary, $\phi(z)$ rotates counterclockwise on
average. From this, we can see the phase zero point on the boundary in the
spin-up solution $\phi(z)e^{i\omega(B/\mu)t}$ moves clockwise and the phase
zero point on the boundary in the spin-down solution
$\phi(z)e^{-i\omega(B/\mu)t}$ moves counterclockwise, spin and momentum are
coupled in $0<B\ll\mu$ material. Note that this argument doesn’t depend on the
shape of the material, implying the wave will not be backscattered by defects
on the boundary.

##  X Mode analysis in 2D $B/\mu\ll 1$ continuum materials

Here we derive qualitative behaviors of the eigenfrequencies and eigenmode
shapes for a homogeneous 2D $B/\mu\ll 1$ continuum material, with Lamé
constants $\lambda=-\mu$, within an arbitrary shaped, bounded simply connected
region $\Omega\subset\mathbb{C}$.

First, we examine its high frequency behavior. When $B/\mu=0$, from (S87) we
see $\rho\omega^{2}/\mu$ is the eigenvalue of the laplacian operator
$-\nabla^{2}$ with Dirichlet boundary conditions. Weyl’s law states that

$\lim_{n\rightarrow\infty}\frac{N(E)}{E}=\frac{|\Omega|}{2\pi},$ (S98)

where $N(E)$ is the number of eigenvalues below $E$ and $|\Omega|$ is the area
of $\Omega$. Let $\omega_{n}$ be the $n^{th}$ eigenfrequency, it follows from
Weyl’s law [57] that

$\omega_{n}\approx\sqrt{\frac{2\pi n\mu}{|\Omega|\rho}}\propto\sqrt{n},$ (S99)

when $n$ is sufficiently large. When $0<B/\mu\ll 1$,
$\omega_{n}\propto\sqrt{n}$ should still hold as $0<B/\mu\ll 1$ is only a
small perturbation. When the nonzero eigenvalue $\rho\omega^{2}/\mu$ is
present, the eigenfunction $h(z)$ exhibits a large value in the bulk. From
(S86), we can deduce that the mode shape $\phi(z)$ also exhibits a large value
in the bulk, hence confirming it as a bulk mode. Let $E_{1}$ be the smallest
eigenvalue of (S87), in nonzero but $B/\mu\ll 1$ limit, the bulk mode first
appears at frequency $\sqrt{\frac{E_{1}\mu}{\rho}}$.

Second, we exam frequency behavior below $\sqrt{\frac{E_{1}\mu}{\rho}}$ for
$0<B/\mu\ll 1$. Let $k_{n}$ be the wave vector of the $n^{th}$ mode along the
boundary, we then have

$k_{n}L=2\pi N_{0}(n),$ (S100)

where $L$ is the perimeter of $\Omega$ and $N_{0}(n)$ is the number of turns
$\phi_{n}(z)$ (the $n^{th}$ mode) winds. From

$\frac{\omega_{n}}{k_{n}}=v_{phase},$ (S101)

where

$\displaystyle v_{phase}=\xi\left(\frac{B}{\mu}\right)\sqrt{\mu/\rho},$ (S102)

is the phase velocity satisfying

$\displaystyle\xi\left(\frac{B}{\mu}\right)^{6}-8\xi\left(\frac{B}{\mu}\right)^{4}+(8+\frac{16\frac{B}{\mu}}{1+\frac{B}{\mu}})\xi\left(\frac{B}{\mu}\right)^{2}-\frac{16\frac{B}{\mu}}{1+\frac{B}{\mu}}=0,$
(S103)

for Rayleigh waves, we know

$\omega_{n}=\frac{2\pi N_{0}(n)v_{phase}}{L}.$ (S104)

Since the shape doesn’t have any spatial symmetry, $\omega_{n}$, in general,
can only have double degeneracy (corresponding to the degeneracy of spin up
and spin down states) and the $(2n-1)^{th}$ edge mode and the $2n^{th}$ edge
mode share a common mode shape $\phi_{n}(z)$. Assuming there is only one
$\phi_{n}(z)$ for each $n$ (we have not been able to establish a rigorous
relation for this), we then have $N_{0}(n)=\frac{n-2}{2}$ for even $n>3$,
$N_{0}(n)=\frac{n-3}{2}$ for odd $n>3$ and $N_{0}(n)\approx\frac{n}{2}$.
Therefore, at frequency below $\sqrt{\frac{E_{1}\mu}{\rho}}$, we have

$\omega_{n}=\frac{\pi nv_{phase}}{L}\propto n.$ (S105)

The simulation in Fig. 2(a) of the main text uses $L=1.48m$,
$B_{3D}/\mu_{3D}=0.01$. In that case, $v_{phase}=21.33$ m/s and

$\omega_{n}=\frac{\pi nv_{phase}}{L}=45.36n,$ (S106)

which agrees with the simulated fitting slope ($43.33$) of the green dots in
Fig. 3(a). Since we have $N_{0}(n)\approx\frac{n}{2}$, the mode shape of the
$(2n-1)^{th}$ and $2n^{th}$ edge mode $\phi_{n}(z)$ has $n$ zeros according to
(S97). Hence,

$\phi_{n}(z)=a_{n}z^{n}+...+a_{1}z+a_{0}$ (S107)

with some large $a_{n}$, because the $n$ zeros lie inside $\Omega$. Since the
function $z^{n}$ is more localized at $\partial\Omega$ than $z^{m}$ for $m<n$,
we know that high frequency edge modes are more localized, which agrees with
Fig. 3(b) in the main text.

##  XI Long wavelength behavior in Maxwell lattices

Figure S1: Dispersion relations of a 2D material strip in the low-$k$ regime
for (a) a continua strip ($B_{3D}/\mu_{3D}=0.1$) and (c) a discrete lattice
($B_{3D}/\mu_{3D}=0.173$), where the strip is free on left and right
boundaries and periodic on the top and bottom. (b,d) Mode profiles at
$k=0.05{\pi}/{a}$.

In Fig. 4(b) in the main text, we have identified two seemingly overlapping
branches, which we denoted as Rayleigh modes. We note that, particularly at
long wavelengths (low $k$), they can also be considered as Lamb modes due to
the presence of the two free surfaces (right and left sides of the simulated
domain shown in main text Fig. 4(d)). In Fig. S1(a,c), we show the dispersion
for two 2D strips of moduli $B_{3D}/\mu_{3D}>0$ with periodic boundary
conditions on the top and bottom and free boundary conditions on the left and
right surfaces (simulated geometries shown in Fig. S1(b,d)). Figures S1(a,b)
are of a continuum and Figs. S1(c,d) are of the same Maxwell lattice strip
considered in the main text Fig. 4(b-d), but focusing in on low frequencies
and wavenumbers. At low $k$, we can now see that the two Rayleigh modes split
into two Lamb modes: one characterized by antisymmetric behavior and parabolic
dispersion, and another characterized by symmetric behavior with linear
dispersion. At low $k$, the symmetric mode approaches the longitudinal
wavespeed. The higher branches, which we referred to as bulk modes, are also
commonly known as “higher order Lamb modes.”

There are a few noteworthy properties concerning these Lamb waves at small
$k$. First, a significant distinction exists between 2D and 3D scenarios. In
3D, Lamb waves always exhibit quadratic dispersion at small $k$, regardless of
the value of $B$. However, in 2D, Lamb waves only arise when $B$ is finite. At
$B=0$ in 2D, there is no hybridization between the Rayleigh waves at the two
edges, and thus Lamb waves are not observed. Second, the issue of
hybridization between modes localized to two independent edges or domain walls
is not unique to our study or system. Such hybridization arises in various
systems, including quantum Hall and quantum valley Hall edge states. When two
1D edges or domain walls are in close proximity, their hybridization leads to
the formation of symmetric and anti-symmetric modes, subsequently resulting in
a loss of one-way transport. Often such effects are neglected, parameters and
geometries are chosen where this hybridization is weak. In our study, we
follow the same approach, ensuring the parameter range studied avoids
significant hybridization, namely by avoiding small $k$ and approaching $B=0$.
However, it is important to note that except for the continuum case idealized
to $B=0$ or infinitely small wavelengths, such hybridization will contribute
to reduction of the edge mode backscattering immunity.

##  XII Additional experimental details

###  XII.1 Mortise and tenon assembly

An example of CNC machined polycarbonate parts for lattice assembly via press
fit is shown in Fig.S2.

Figure S2:  Example of CNC machined sublattice elements with mortise and tenon
structures to form solid triangles.

###  XII.2 Energy ratio calculation

The energy ratio calculated in experiments between the edge waves propagating
in opposite directions is calculated by comparing the sums of the squared
velocities of all edge triangles’ centroids on the left and right sides of the
actuator (shown in Fig. S3) over the time period of $0.01-0.15$ s, before the
two edge waves encounter each other. It should be noted that this ratio does
not correspond to the ratio $|C_{-k}/C_{k}|$, which predicts the relative
amplitude of the opposite propagating modes in the $t\to\infty$.

Figure S3: Triangles on the edge to the left and right of the excitation point
(black circle) used for calculating energy ratio $E_{L}/E_{R}$ shown in Fig.
4(e-j) and plotting average kinetic energy ($v^{2}$) in Fig. S5. Numbers in
the red and green solid circles are indices used for data processing in
MATLAB. $S1$ and $S2$ mark the triangles on the edges equidistant from the
actuated point used for tracking the trajectory phases shown in Fig. S5.

###  XII.3 Poisson’s ratio measurement

The Poisson’s ratio of our manufactured Maxwell lattice was measured based on
manual compression of the lattice and image tracking of axial compression
(perpendicular to the edge) and transverse (lateral) expansion. The fitted
slope gives $\nu\simeq-0.488$ (Fig. S4), corresponding to $B/\mu=0.344$
($B_{3D}/\mu_{3D}=0.173$).

Figure S4: Measured Poisson’s ratio of the auxetic Maxwell lattice.
$\epsilon_{axial}$ and $\epsilon_{trans}$ represent the strains along and
perpendicular to the compression, respectively.

###  XII.4 Additional experiment and simulation comparison

In Fig. S5, we present the raw data for the simulations and measurements of
spin-momentum locking and one-way edge waves in Maxwell lattices, which were
reported in the main text Fig. 4 (e-j). As explained in the main text, it was
observed that the opposite propagating waves exhibit opposite spins. To
further analyze this phenomenon, the velocity phases ($arg(v)$) of $S1$ and
$S2$ were examined. It was consistently found that S1 rotates in a clockwise
direction, while S2 rotates counterclockwise before the two edge waves
encounter each other. This observation confirms the presence of spin-momentum
locking in the edge waves. To quantify the energy distribution between the
left- (L) and right- (R) propagating edge modes, the energy ratio
${E_{L}}/{E_{R}}$ was calculated. This ratio was obtained by integrating the
average squared velocity (over the time range $0.01-0.15$ s) of the triangles
located on the edge to the left and right of the excitation point, as depicted
in Fig. S3. Considering that the group velocity at $5.5$ Hz is approximately
$0.03$ m/s and the distance from the actuation source to the intersection of
the two waves is around $0.45$ m, the two edge waves encounter each other
after $0.15$ s.

Figure S5: Experiments and simulation raw data for one-way edge waves in
Maxwell lattices. (a) Time domain simulations and (b) experimental
measurements of one-way edge waves in auxetic Maxwell lattices. In each panel,
the top figures display the velocity vector phases ($arg(v)$, which are
shifted to start from $0$) of the centroids of triangles $S1$ and $S2$ (as
indicated in Fig.S3) in the time domain under different excitation types. The
bottom figures illustrate the average kinetic energy ($v^{2}$) in the same
time range, calculated based on the triangles on the edge, as shown in Fig.
S3. The term ’ratio’ represents the energy ratio $E_{L}/E_{R}$.

##  XIII Additional Maxwell lattice scattering simulations

Additional FEM simulations are performed, where we deliberately introduce
disorder in our simulated Maxwell lattice by removing or adding unit cells at
the boundary or creating cavities within the bulk. As depicted in Fig. S6, the
results demonstrate that this does not significantly disrupt the one-way edge
transport along the edge.

Figure S6: (a-c) Robustness of edge wave transport in the auxetic Maxwell
lattice by simulating various types of defects, where three triangles were
added to the left edge ((a)), three triangles were removed from the left edge
((b)), and four triangles were removed from the bulk ((c)). The excitation
position is denoted by a red asterisk, and actuation is carried out in a
counterclockwise direction.