# Nonreciprocity and mode conversion in a spatiotemporally modulated elastic
wave circulator

Benjamin M. Goldsberry  Applied Research Laboratories, The University of Texas
at Austin, Austin, Texas 78758, USA  Samuel P. Wallen  Michael R. Haberman
Applied Research Laboratories, The University of Texas at Austin, Austin,
Texas 78758, USA  Walker Department of Mechanical Engineering, The University
of Texas at Austin, Austin, Texas 78712-1591, USA

###### Abstract

Acoustic and elastic metamaterials with time- and space-dependent material
properties have received great attention recently as a means to break
reciprocity for propagating mechanical waves, achieving greater directional
control. One nonreciprocal device that has been demonstrated in the fields of
acoustics and electromagnetism is the circulator, which achieves unirotational
transmission through a network of ports. This work investigates an elastic
wave circulator composed of a thin elastic ring with three semi-infinite
elastic waveguides attached, creating a three-port network. Nonreciprocity is
achieved for both longitudinal and transverse waves by modulating the elastic
modulus of the ring in a rotating fashion. Two numerical models are derived
and implemented to compute the elastic wave field in the circulator. The first
is an approximate model based on coupled mode theory, which makes use of the
known mode shapes of the unmodulated system. The second model is a finite
element approach that includes a Fourier expansion in the harmonics of the
modulation frequency and radiation boundary conditions at the ports. The
coupled mode model shows excellent agreement with the finite element model and
the conditions on the modulation parameters that enable a high degree of
nonreciprocity are presented. This work demonstrates that it is possible to
create an elastic wave circulator that enables nonreciprocal mode-splitting of
an incident longitudinal wave into outgoing longitudinal and transverse waves.

##  I Introduction

Acoustic and elastic wave reciprocity is a fundamental principle that requires
the symmetry of propagating waves between two points in space when propagation
takes place in linear time-invariant (LTI) media [1, 2]. The simplest
expression of reciprocity is that the propagation of waves from source to
receiver will be identical in both phase and magnitude if the source and
receiver are interchanged. Current research in nonreciprocal wave propagation
therefore aims to break this symmetry between source and receiver, enabling
direction-dependent wave phenomena such as one-way wave propagation and
vibration isolation [3]. Various techniques that have been investigated in
enabling nonreciprocity have often relied on violating the assumption of a LTI
system. For example, nonlinear wave propagation in geometrically-asymmetric
media have been shown to exhibit nonreciprocity [4, 5, 6]. In addition,
spatiotemporal modulation of the material properties [7, 8, 9, 10, 11, 12, 13,
14] via the utilization of active elements [15, 16, 17, 10], including
modulation of paired loss-gain media in non-Hermitian systems [18, 19],
enables nonreciprocal effects for linear waves.

Many nonreciprocal mechanical systems have been inspired by analogous multi-
port electronic devices [20]. An isolator is a two-port system that allows for
the transmission of waves in one direction, and has been realized for acoustic
waves in systems that utilize fluid flow [21], spatiotemporal modulation of
material properties [7], and nonlinearity [22, 4]. The gyrator is another two-
port system that manipulates the phase of a wave in one direction, and has
recently been achieved in a waveguide with fluid flow [23]. Finally,
circulators are three-port devices that transmit waves in a unirotational
fashion, and have been realized using both fluid flow [24] and spatiotemporal
modulation of the volumes of coupled cavities [25]. Nonreciprocal circulators,
which are typically composed of three ports with three-fold rotational
symmetry connected to coupled acoustic or elastic resonators, are the focus of
this work. These devices are primarily utilized to isolate source and receiver
channels for full-duplex communication devices [26]. Previous works have
limited the analysis of mechanical circulators to the acoustic (scalar) case
[24, 25], lumped-element mechanical systems that only admit one unique
resonance [25, 27], and elastic devices utilizing an angular momentum bias
imposed by steady-state rotation [28]. However, little attention has been
given to more general elastic wave circulators with spatiotemporally modulated
material properties, whose dynamics may be greatly enriched by coupling
between modes, including multiple wave polarizations in the circulator and all
ports.

We investigate an elastic wave circulator composed of a thin elastic ring
connecting three semi-infinite elastic waveguides, creating a geometrically-
symmetric three-port network. Nonreciprocity is achieved for both longitudinal
and transverse waves by modulating the elastic modulus of the ring in a
rotating fashion. Details of the circulator geometry and modulation of the
elastic modulus is given in Sec. II. Two numerical models are derived and
implemented to compute the elastic wave field in the circulator: The first
model is an approximate coupled mode model (CMM), presented in Sec. III, which
is based on previous research investigating nonreciprocal vibrations in finite
Euler beams [29]. The CMM is a reduced-order model of the circulator and is
therefore computationally efficient, making the CMM well-suited for efficient
design iteration. However, the assumptions employed for the CMM impose
restrictions on the circulator geometry, limiting its ability to capture the
full design space. The second model, presented in Sec. IV, addresses these
limitations by adapting a finite element approach (FEA) that was previously
derived by the authors to explore nonreciprocal elastic wave propagation in
modulated elastic metamaterials [13]. In contrast to the CMM, the FEA requires
a discretized mesh of the circulator geometry and is therefore computationally
expensive, which limits a thorough exploration of the circulator design space.
We therefore calibrate the CMM with the FEA in Sec. V in order to accurately
utilize the CMM to investigate circulator designs, and we show with both
models that the proposed elastic wave circulator exhibits a large degree of
nonreciprocity and can act as a nonreciprocal Lamb wave mode converter. We
then summarize the findings of this work in Sec. VI.

##  II Model Definition

Figure 1: Elastic wave circulator geometry.

Consider the elastic circulator configuration shown in Fig. 1, which consists
of an elastic ring attached to three semi-infinite, external elastic
waveguides acting as input and output ports. The ring has constant radial
thickness $h$ and constant radius of curvature $R$ at its mid-surface, such
that the location of any point on the ring may be determined by the angle
$\theta$. The semi-infinite waveguides each have thickness $h_{\mathrm{w}}$
and are attached at locations $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$,
separated by $120^{\circ}$, creating a three-port network with 3-fold axial
symmetry.

The ring material is considered to be a linear, isotropic solid with modulated
Young’s modulus $E(\theta,t)$, constant Poisson’s ratio $\nu$, and constant
mass density $\rho$. The corresponding material properties of the port
waveguides are constants denoted by $E_{\mathrm{w}}$, $\nu_{\mathrm{w}}$, and
$\rho_{\mathrm{w}}$, respectively, and may generally differ from the ring. The
Young’s modulus of the ring has the assumed form

$E(\theta,t)=E_{0}+E_{\mathrm{m}}E^{\prime}(\theta,t),$ (1)

where $E_{0}$ is the mean (static) Young’s modulus, $E_{\mathrm{m}}$ is the
modulation amplitude, and $E^{\prime}(\theta,t)$ is the variation, or
modulation, of the Young’s modulus about the mean value. The modulation is
restricted to be a periodic function of time,
$E^{\prime}(\theta,t+T_{\mathrm{m}})=E^{\prime}(\theta,t)$, where
$T_{\mathrm{m}}$ is the fundamental period of the modulation. The fundamental
angular frequency of the modulation then follows as
$\omega_{\mathrm{m}}=2\pi/T_{\mathrm{m}}$.

##  III Coupled mode model

###  III.1 Equations of motion

Here, we derive an approximate model for the circulator based on coupled mode
theory. We begin by assuming sufficiently small drive frequencies for which
the wavelength is much larger than the thickness $h$, such that the elastic
ring can be well approximated as a curved Euler-Bernoulli beam with plane
strain assumptions for the out-of-plane direction. Therefore, the complete in-
plane displacement field of the ring can be described by the radial and
tangential displacements, denoted by $W(\theta,t)$ and $V(\theta,t)$,
respectively. The coupled partial differential equations of motion for the
displacements are written in dimensionless form as [30]

$\frac{R^{2}}{c_{0}^{2}}\frac{\partial^{2}w}{\partial
t^{2}}=\gamma^{2}\frac{\partial^{2}}{\partial\theta^{2}}\left[\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}-\frac{\partial^{2}w}{\partial\theta^{2}}\right)\right]\\\
-\bar{E}(\theta,t)\left(w+\frac{\partial v}{\partial\theta}\right)+f_{r},$ (2)
$\frac{R^{2}}{c_{0}^{2}}\frac{\partial^{2}v}{\partial
t^{2}}=\gamma^{2}\frac{\partial}{\partial\theta}\left[\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}-\frac{\partial^{2}w}{\partial\theta^{2}}\right)\right]\\\
+\frac{\partial}{\partial\theta}\left[\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}+w\right)\right]+f_{\theta},$ (3)

where $w(\theta,t)=W(\theta,t)/U$ and $v(\theta,t)=V(\theta,t)/U$ are the
dimensionless radial and tangential displacements, respectively, and are
normalized by an arbitrary reference displacement $U$. The relevant
elastodynamic parameters in Eqns. (2)-(3) are the dimensionless Young’s
modulus $\bar{E}(\theta,t)=E(x,t)/E_{0}$, the longitudinal sound speed
$c_{0}=\sqrt{E_{0}/\rho}$, and the dimensionless parameter

$\gamma=R^{-1}\sqrt{D/(E_{0}h)},$ (4)

which represents the flexural rigidity of the curved plate, where
$D=E_{0}h^{3}/[12(1-\nu^{2})]$ is the flexural rigidity of a planar plate with
the same thickness. The external loading is defined by the forces $f_{r}$ and
$f_{\theta}$, which represent normalized radial and tangential contributions,
respectively, and include sources and radiation losses due to the attached
elastic waveguide ports. The external forces from the waveguide ports can be
idealized as point forces if the external port waveguide thicknesses are much
less than the ring thickness ( $h_{\mathrm{w}}\ll h$). In this case, the
external waveguide ports are assumed to be stress free in the tangential
direction, i.e., $f_{\theta}=0$. Therefore, the only input wave considered in
this work is a time-harmonic radial input force at port 1, which is caused by
an incident longitudinal mode in the port waveguide. The radiation losses due
to the waveguide ports are approximated using the mechanical longitudinal
impedance of a plate [31], such that the normalized exterior forcing term is
written as

$f_{r}=\frac{2RF_{0}}{E_{0}hU}\delta(\theta-\theta_{1})e^{-i\omega_{0}t}+\sum\limits_{q=1}^{3}\frac{R\bar{Z}}{c_{0}}\frac{\partial\bar{w}}{\partial
t}\delta(\theta-\theta_{q}),$ (5)

where the first term on the right-hand side of Eq. (5) is the force due to the
incident longitudinal wave in port 1 ($\theta_{1}$ in Fig. 1), and the second
term on the right-hand side of Eq. (5) is the radiation loss from all three
port waveguides ($\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ in Fig. 1).
Furthermore, $\delta(\theta)$ is the Dirac delta function,
$\bar{Z}=\rho_{\mathrm{w}}c_{\mathrm{w}}h_{\mathrm{w}}/(\rho c_{0}h)$ is the
dimensionless radiation impedance, where
$c_{\mathrm{w}}=\sqrt{E_{\mathrm{w}}/[\rho_{\mathrm{w}}(1-\nu_{\mathrm{w}}^{2})]}$
is the low-frequency limit of the longitudinal wave speed of the port
waveguide, and $F_{0}$ is the amplitude of the time-harmonic input force with
frequency $\omega_{0}$.

###  III.2 Coupled mode formulation

Due to the periodicity of $E^{\prime}(\theta,t)$ in time, we can expand
$E^{\prime}(\theta,t)$ as a truncated Fourier series with $2P+1$ terms, such
that the dimensionless Young’s modulus is written as

$\bar{E}(\theta,t)=1+\alpha\sum\limits_{\begin{subarray}{c}p=-P\\\ p\neq
0\end{subarray}}^{P}\hat{E}^{p}(\theta)e^{-ip\omega_{\mathrm{m}}t},$ (6)

where $\hat{E}^{p}(\theta)$ are the Fourier components of
$E^{\prime}(\theta,t)$, and $\alpha=E_{\mathrm{m}}/E_{0}$ is the dimensionless
modulation amplitude. As a result of the time-periodic form of the Young’s
modulus and the incident wave, Eqns. (2)-(3) possess steady-state solutions of
the form

$\displaystyle w(\theta,t)$
$\displaystyle=\sum\limits_{p=-P}^{P}\hat{w}^{p}(\theta)e^{-i(\omega_{0}+p\omega_{\mathrm{m}})t},$
(7) $\displaystyle v(\theta,t)$
$\displaystyle=\sum\limits_{p=-P}^{P}\hat{v}^{p}(\theta)e^{-i(\omega_{0}+p\omega_{\mathrm{m}})t},$
(8)

where $\hat{w}^{p}(\theta)$ and $\hat{v}^{p}(\theta)$ are the angle-dependent
radial and tangential displacements for each frequency component
$(\omega_{0}+p\omega_{\mathrm{m}})$. Since the displacements of the ring must
be continuous for all $\theta$, including the endpoints $w(0,t)=w(2\pi,t)$ and
$v(0,t)=v(2\pi,t)$, the spatial profiles can be expanded as a truncated
Fourier series with $2N+1$ terms and basis $\psi_{n}(\theta)=e^{in\theta}$,

$\displaystyle\hat{w}^{p}$
$\displaystyle=\sum\limits_{n=-N}^{N}\hat{w}^{p}_{n}\psi_{n}(\theta),$ (9)
$\displaystyle\hat{v}^{p}$
$\displaystyle=\sum\limits_{n=-N}^{N}\hat{v}^{p}_{n}\psi_{n}(\theta),$ (10)

which, when substituted in Eqns. (7)-(8), yields the total displacement
solution

$\displaystyle w(\theta,t)$
$\displaystyle=\sum\limits_{p=-P}^{P}\sum\limits_{n=-N}^{N}\hat{w}^{p}_{n}e^{i\left[n\theta-\left(\omega_{0}+p\omega_{\mathrm{m}}\right)t\right]},$
(11) $\displaystyle v(\theta,t)$
$\displaystyle=\sum\limits_{p=-P}^{P}\sum\limits_{n=-N}^{N}\hat{v}^{p}_{n}e^{i\left[n\theta-\left(\omega_{0}+p\omega_{\mathrm{m}}\right)t\right]}.$
(12)

Equations (11)-(12) show that the total displacement solutions are a sum of
clockwise ($n<0$) and counter-clockwise ($n>0$) propagating plane waves at the
frequencies $(\omega_{0}+p\omega_{\mathrm{m}})$ with amplitudes
$\hat{w}^{p}_{n}$ and $\hat{v}^{p}_{n}$.

Equations for the plane wave amplitudes are derived using a weak formulation
of the equations of motion, which results in a symmetric system of equations
for the plane wave amplitude values, and also allows for arbitrary and
possibly discontinuous stiffness modulation profiles, such as the traveling
modulation form considered in [8]. The weak formulation of the equations of
motion is derived by multiplying Eqns. (2)-(3) by
$\psi^{*}_{m}(\theta)=e^{-im\theta}$ and integrating over $\theta$:

$\frac{R^{2}}{c_{0}^{2}}\int\limits_{0}^{2\pi}\frac{\partial^{2}w}{\partial
t^{2}}\psi^{*}_{m}\,d\theta=\\\
\gamma^{2}\int\limits_{0}^{2\pi}\frac{\partial^{2}}{\partial\theta^{2}}\left[\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}-\frac{\partial^{2}w}{\partial\theta^{2}}\right)\right]\psi^{*}_{m}\,d\theta\\\
-\int\limits_{0}^{2\pi}\bar{E}(\theta,t)\left(w+\frac{\partial
v}{\partial\theta}\right)\psi^{*}_{m}\,d\theta+\int\limits_{0}^{2\pi}f_{r}\psi^{*}_{m}\,d\theta,$
(13)
$\frac{R^{2}}{c_{0}^{2}}\int\limits_{0}^{2\pi}\frac{\partial^{2}v}{\partial
t^{2}}\psi^{*}_{m}\,d\theta=\\\
\gamma^{2}\int\limits_{0}^{2\pi}\frac{\partial}{\partial\theta}\left[\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}-\frac{\partial^{2}w}{\partial\theta^{2}}\right)\right]\psi^{*}_{m}\,d\theta\\\
+\int\limits_{0}^{2\pi}\frac{\partial}{\partial\theta}\left[\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}+w\right)\right]\psi^{*}_{m}\,d\theta.$ (14)

Integration by parts is utilized twice on the first integral on the right hand
side of Eq. (13), and once on the first and second integrals on the right hand
side of Eq. (14), to yield

$\frac{R^{2}}{c_{0}^{2}}\int\limits_{0}^{2\pi}\frac{\partial^{2}w}{\partial
t^{2}}\psi^{*}_{m}\,d\theta=\gamma^{2}\int\limits_{0}^{2\pi}\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}-\frac{\partial^{2}w}{\partial\theta^{2}}\right)\frac{\partial^{2}\psi^{*}_{m}}{\partial\theta^{2}}\,d\theta\\\
-\int\limits_{0}^{2\pi}\bar{E}(\theta,t)\left(w+\frac{\partial
v}{\partial\theta}\right)\psi^{*}_{m}\,d\theta+\int\limits_{0}^{2\pi}f_{r}\psi^{*}_{m}\,d\theta,$
(15)
$\frac{R^{2}}{c_{0}^{2}}\int\limits_{0}^{2\pi}\frac{\partial^{2}v}{\partial
t^{2}}\psi^{*}_{m}\,d\theta=-\gamma^{2}\int\limits_{0}^{2\pi}\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}-\frac{\partial^{2}w}{\partial\theta^{2}}\right)\frac{\partial\psi^{*}_{m}}{\partial\theta}\,d\theta\\\
-\int\limits_{0}^{2\pi}\bar{E}(\theta,t)\left(\frac{\partial
v}{\partial\theta}+w\right)\frac{\partial\psi^{*}_{m}}{\partial\theta}\,d\theta.$
(16)

Note that the boundary terms that arise from integration by parts vanish due
to the spatial periodicity of the displacement solution and Young’s modulus.
Finally, the equations for the plane wave amplitudes are found by substituting
Eqns. (6), (11), and (12) into Eqns. (15)-(16) and utilizing the orthogonality
of the Fourier series in time and space. The resulting expressions can be
written compactly, in matrix-vector form, as

$\left(\left(\bar{\omega}_{0}+p\bar{\omega}_{\mathrm{m}}\right)^{2}\begin{bmatrix}\mathbf{I}&0\\\
0&\mathbf{I}\end{bmatrix}-\begin{bmatrix}\mathbf{K}^{0}_{\mathrm{W,W}}&\mathbf{K}^{0}_{\mathrm{W,V}}\\\
\mathbf{K}^{0}_{\mathrm{V,W}}&\mathbf{K}^{0}_{\mathrm{V,V}}\end{bmatrix}+i\bar{Z}\left(\bar{\omega}_{0}+p\bar{\omega}_{\mathrm{m}}\right)\begin{bmatrix}\mathbf{C}&0\\\
0&0\end{bmatrix}\right)\begin{bmatrix}\bm{w}^{p}\\\
\bm{v}^{p}\end{bmatrix}-\alpha\sum\limits_{\begin{smallmatrix}s=-P\\\ s\neq
p\end{smallmatrix}}^{P}\begin{bmatrix}\mathbf{K}^{p-s}_{\mathrm{W,W}}&\mathbf{K}^{p-s}_{\mathrm{W,V}}\\\
\mathbf{K}^{p-s}_{\mathrm{V,W}}&\mathbf{K}^{p-s}_{\mathrm{V,V}}\end{bmatrix}\begin{bmatrix}\bm{w}^{s}\\\
\bm{v}^{s}\end{bmatrix}=\delta_{p0}\begin{bmatrix}\bm{F}\\\ 0\end{bmatrix},\\\
\quad p\in[-P,P],$ (17)

where
$\bm{w}^{p}=[\hat{w}^{p}_{-N},\hat{w}^{p}_{-N+1},...,\hat{w}^{p}_{N}]^{\mathrm{T}}$
and
$\bm{v}^{p}=[\hat{v}^{p}_{-N},\hat{v}^{p}_{-N+1},...,\hat{v}^{p}_{N}]^{\mathrm{T}}$
are the vectors of plane wave amplitudes; $\bar{\omega}_{0}=R\omega_{0}/c_{0}$
and $\bar{\omega}_{\mathrm{m}}=R\omega_{\mathrm{m}}/c_{0}$ are the
dimensionless drive and modulation frequencies, respectively; and the matrix
$\mathbf{C}$, which represents the radiation loss in the ring due to the port
waveguides, has the components
$C_{mn}=(2\pi)^{-1}\sum\limits_{q=1}^{3}e^{in\theta_{q}}e^{-im\theta_{q}}.$
The forcing vector $\bm{F}$ in Eq. (17) has the values
$F_{n}=(2\pi)^{-1}e^{-in\theta_{\mathrm{s}}}$, and the coupling matrices are

$\left(K^{p-s}_{W,W}\right)_{mn}=-\gamma^{2}\int\limits_{0}^{2\pi}\hat{E}^{p-s}(\theta)\frac{\partial^{2}\psi_{n}}{\partial\theta^{2}}\frac{\partial^{2}\psi_{m}^{*}}{\partial\theta^{2}}\,d\theta\\\
-\int\limits_{0}^{2\pi}\hat{E}^{p-s}(\theta)\psi_{n}\psi_{m}^{*}\,d\theta,$
(18a) $\left(K^{p-s}_{V,W}\right)_{mn}=\left(K^{p-s}_{W,V}\right)_{mn}=\\\
\gamma^{2}\int\limits_{0}^{2\pi}\hat{E}^{p-s}(\theta)\frac{\partial\psi_{n}}{\partial\theta}\frac{\partial^{2}\psi_{m}^{*}}{\partial\theta^{2}}\,d\theta-\int\limits_{0}^{2\pi}\hat{E}^{p-s}(\theta)\frac{\partial\psi_{n}}{\partial\theta}\psi_{m}^{*}\,d\theta,$
(18b)
$\left(K^{p-s}_{V,V}\right)_{mn}=-(\gamma^{2}+1)\int\limits_{0}^{2\pi}\hat{E}^{p-s}(\theta)\frac{\partial\psi_{n}}{\partial\theta}\frac{\partial\psi_{m}^{*}}{\partial\theta}.$
(18c)

The matrices inside the parentheses in Eq. (17) are the so-called mass,
stiffness, and loss matrices, respectively, and are identical to the dynamic
matrices for the unmodulated system. The terms in Eq. (17) proportional to
$\alpha$ arise due to the stiffness modulation, and represent the coupling
between tangential and longitudinal wave modes within the ring at each
frequency $(\bar{\omega}_{0}+p\bar{\omega}_{\mathrm{m}})$. Equation (17) is
numerically solved given the number of retained spatial plane waves $N$,
frequencies $P$, and dimensionless parameters
$(\bar{\omega}_{0},\bar{\omega}_{\mathrm{m}},\bar{Z},\alpha)$. The results are
shown in Sec. V, and are compared with a finite element approach derived in
Sec. IV.

##  IV Finite element approach

While the coupled mode model derived in Sec. III is computationally efficient,
it has a number of assumptions, namely that the curvature $(h/R)$ and drive
frequency must be sufficiently small such that Euler-Bernoulli theory holds.
In addition, the port waveguide thickness must be much smaller than the ring
thickness (i.e., $h_{\mathrm{w}}/h\ll 1$) such that the forcing and radiation
impedance approximations in Eq. (5) are valid. To alleviate these assumptions
and obtain a model with applicability to more general geometries, we also
develop a finite element approach (FEA), which numerically solves the
continuum elastodynamic equations and takes into account the exact coupling
between the ring and exterior waveguide ports. To this end, we extend the FEA
previously derived in Ref. [13] to the present case. The primary modification
that must be made is to enforce outgoing elastic waves for the reflected and
transmitted fields that propagate in the waveguide ports, which is
accomplished by adapting a mode-matching technique [32].

###  IV.1 Derivation

The weak formulation for the elastodynamic equations with spatiotemporally-
modulated elastic moduli is written in indicial notation, with the Einstein
summation convention applied to subscript indices only, as [13]

$\sum\limits_{p,s=-P}^{P}-\int\limits_{\Omega}\hat{v}^{p}_{i,j}\hat{L}^{p-s}_{ijkl}\hat{u}_{k,l}\,d\Omega+\left(\omega_{0}+p\omega_{\mathrm{m}}\right)^{2}\int\limits_{\Omega}\rho\hat{u}^{p}_{i}\hat{v}^{p}_{i}\,d\Omega\\\
+\sum\limits_{q=1}^{3}\int\limits_{\Gamma_{q}}\hat{L}^{0}_{ijkl}\hat{u}^{p,\mathrm{port}}_{k,l}\hat{v}^{p}_{i}n_{j}\,dS=-\int\limits_{\Gamma_{\mathrm{port1}}}\hat{L}^{0}_{ijkl}\hat{u}^{0,\mathrm{inc}}_{k,l}\hat{v}^{0}_{i}n_{j}\,dS,$
(19)

where $\hat{u}^{p}_{i}$ and $\hat{v}^{p}_{i}$ are the trial and test
displacements, respectively, at frequency $(\omega_{0}+p\omega_{\mathrm{m}})$,
$\Omega$ represents the union of the ring and waveguide domains in Fig. 1, and
$\Gamma_{q}$, $q=1,2,3$ are the three exterior waveguide boundaries caused
from the truncation of the semi-infinite waveguides to create a finite
computational domain. The material properties in Eq. (19) include the mass
density $\rho(\bm{x})$ and the fourth-order tangent modulus tensor

$\hat{L}^{s}_{ijkl}(\bm{x})=\hat{\lambda}^{s}(\bm{x})\delta_{ij}\delta_{kl}+\hat{\mu}^{s}(\bm{x})\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right),$
(20)

where
$\hat{\lambda}^{s}(\bm{x})=\hat{E}^{s}(\bm{x})\nu/\left[(1+\nu)(1-2\nu)\right]$
and $\hat{\mu}^{s}(\bm{x})=\hat{E}^{s}(\bm{x})/\left[2(1+\nu)\right]$ are the
Fourier-transformed first and second Lamé parameters, respectively. We note
that the $\hat{E}^{s}(\bm{x})$ functions are defined in Eq. (6) with
$\theta=\mathrm{atan2}(y,x)$ and $\bm{x}=[x,y]^{\mathrm{T}}$. The unknown
displacements at the port waveguide boundaries at each frequency,
$\hat{\bm{u}}^{p,\mathrm{port}}$, must be determined such that the waves
propagating away from the ring in the port waveguides remain purely outgoing
and do not reflect from the computational boundary. Given the plane-strain
assumptions of the model, we assume that the displacement field in the port
waveguides can be completely represented as a linear combination of Rayleigh-
Lamb waves [30]. The drive frequencies considered in this work are limited to
be below the cut-on frequencies of the higher-order Lamb modes, such that only
the zero order symmetric ($S_{0}$) and antisymmetric ($A_{0}$) Lamb waves
propagate in the port waveguides. The components of the displacement field for
the symmetric Lamb wave, $\bm{\psi}^{\mathrm{sym}}$, are [30]

$\displaystyle\psi^{\mathrm{sym}}_{x}$ $\displaystyle=i\left(\beta\cos(\beta
y)+\frac{B}{C}\xi\cos(\eta y)\right)e^{i\xi x},$ (21)
$\displaystyle\psi^{\mathrm{sym}}_{y}$ $\displaystyle=\left(\xi\sin(\beta
y)-\eta\frac{B}{C}\sin(\eta y)\right)e^{i\xi x},$ (22)
$\displaystyle\frac{B}{C}$ $\displaystyle=\frac{(\xi^{2}-\beta^{2})\sin(\beta
h/2)}{2\xi\eta\sin(\eta h/2)},$ (23)

where $\xi$ is the wavenumber, $\eta^{2}=\omega^{2}/c_{1}^{2}-\xi^{2}$,
$\beta^{2}=\omega^{2}/c_{2}^{2}-\xi^{2}$, $c_{1}=\sqrt{(\lambda+2\mu)/\rho}$
is the bulk longitudinal wave speed of the waveguide material, and
$c_{2}=\sqrt{\mu/\rho}$ is the bulk transverse wave speed in the same
material. Similarly, the components of the displacement field for the
antisymmetric Lamb wave, $\bm{\psi}^{\mathrm{asym}}$, are

$\displaystyle\psi^{\mathrm{asym}}_{x}$ $\displaystyle=-i\left(\beta\sin(\beta
y)-\frac{A}{D}\xi\sin(\eta y)\right)e^{i\xi x},$ (24)
$\displaystyle\psi^{\mathrm{asym}}_{y}$ $\displaystyle=\left(\xi\cos(\beta
y)+\eta\frac{A}{D}\cos(\eta y)\right)e^{i\xi x},$ (25)
$\displaystyle\frac{A}{D}$ $\displaystyle=-\frac{(\xi^{2}-\beta^{2})\cos(\beta
h/2)}{2\xi\eta\cos(\eta h/2)}.$ (26)

The wavenumber $\xi$ satisfies the well-known Rayleigh-Lamb dispersion
relation

$\frac{\tan(\beta h/2)}{\tan(\eta
h/2)}+\left[\frac{4\eta\beta\xi^{2}}{(\xi^{2}-\beta^{2})^{2}}\right]^{\pm
1}=0,$ (27)

where the power of the second term on the left-hand side of Eq. (27) is $+1$
for the symmetric Lamb wave and $-1$ for the antisymmetric Lamb wave. The
total displacement at each output port is therefore
$\bm{u}^{p,\mathrm{port}}=A^{p,\mathrm{sym}}\bm{\psi}^{p,\mathrm{sym}}+A^{p,\mathrm{asym}}\bm{\psi}^{p,\mathrm{asym}}$,
where $A^{p,\mathrm{sym}}$ and $A^{p,\mathrm{asym}}$ are the amplitude
coefficients for the symmetric and antisymmetric Lamb waves, respectively, at
the frequency $(\omega_{0}+p\omega_{\mathrm{m}})$.

The Lamb wave amplitudes are determined by creating an auxiliary equation to
Eq. (19) that weakly enforces the continuity of displacement at each
frequency. This is accomplished by forming a weighted inner product equation
with operator $W$,
$\left<W\bm{\psi}^{p},\bm{u}^{p}-\bm{u}^{p,\mathrm{port}}\right>=0$, where the
weighting operator is chosen to be the integrand of the integrals over the
port faces in Eq. (19) [32]. Expanding $\bm{u}^{p,\mathrm{port}}$ in terms of
the symmetric and antisymmetric Lamb wave shapes yields the auxiliary
equations

$\displaystyle\int\limits_{\Gamma_{q}}L^{0}_{ijkl}\psi_{k,l}^{p,\mathrm{sym}}n_{j}\left(\hat{u}^{p}-A^{p,\mathrm{sym}}\psi^{p,\mathrm{sym}}\right)_{i}\,dS+\int\limits_{\Gamma_{q}}L^{0}_{ijkl}\psi_{k,l}^{p,\mathrm{sym}}n_{j}\left(\hat{u}^{p}-A^{p,\mathrm{asym}}\psi^{p,\mathrm{asym}}\right)_{i}\,dS=0,$
(28)
$\displaystyle\int\limits_{\Gamma_{q}}L^{0}_{ijkl}\psi_{k,l}^{p,\mathrm{asym}}n_{j}\left(\hat{u}^{p}-A^{p,\mathrm{sym}}\psi^{p,\mathrm{sym}}\right)_{i}\,dS+\int\limits_{\Gamma_{q}}L^{0}_{ijkl}\psi_{k,l}^{p,\mathrm{asym}}n_{j}\left(\hat{u}^{p}-A^{p,\mathrm{asym}}\psi^{p,\mathrm{asym}}\right)_{i}\,dS=0,$
(29)

for each waveguide boundary $\Gamma_{q}$. Finally, an augmented system system
of equations is solved for the finite element displacement field and the port
Lamb wave amplitudes at each frequency,

$\begin{bmatrix}\omega_{0}^{2}\mathbf{M}+\omega_{0}\mathbf{C}+\mathbf{K}&\mathbf{\Gamma}\\\
\mathbf{\Gamma}^{\mathrm{H}}&\mathbf{\Lambda}\end{bmatrix}\begin{bmatrix}\bm{u}\\\
\bm{A}\end{bmatrix}=\begin{bmatrix}\tilde{\bm{F}}\\\ \bm{S}\end{bmatrix},$
(30)

where $\bm{u}=[\bm{u}^{-P},\bm{u}^{-P+1},...,\bm{u}^{P}]^{\mathrm{T}}$ is an
augmented vector containing the finite element discretized displacements at
each frequency, and
$\bm{A}=[\bm{A}^{-P},\bm{A}^{-P+1},...,\bm{A}^{P}]^{\mathrm{T}}$ is an
augmented vector of the symmetric and antisymmetric Lamb wave amplitudes at
each port waveguide and frequency. The matrices $\mathbf{M}$, $\mathbf{C}$,
and $\mathbf{K}$ in Eq. (30) are from the discretization of the integrals over
$\Omega$ in Eq. (19) and are given in detail in [13]. The matrix
$\mathbf{\Gamma}$ in Eq. (30) is from the integrals over the waveguide port
surfaces on the left-hand side in Eq. (19), where the superscript $\mathrm{H}$
denotes Hermitian transposition, and the matrix $\mathbf{\Lambda}$ is formed
from the integrals proportional to the Lamb wave amplitude coefficients in
Eqns. (28)-(29). Finally, the source terms on the right-hand side of Eq. (30)
due to the incident wave are given by

$\displaystyle\tilde{\bm{F}}$
$\displaystyle=\int\limits_{\Gamma_{1}}L^{0}_{ijkl}u^{\mathrm{inc}}_{k,l}v^{0}_{i}n_{j}\,dS,$
(31) $\displaystyle\bm{S}$
$\displaystyle=\int\limits_{\Gamma_{1}}L^{0}_{ijkl}u^{p,\mathrm{sym}}_{k,l}u_{i}^{\mathrm{inc}}n_{j}\,dS.$
(32)

It is assumed in this work that the incident field is a symmetric Lamb wave,
$\bm{u}^{\mathrm{inc}}=\bm{\psi}^{\mathrm{sym}}$, which is consistent with the
input force for the coupled mode model discussed in Sec. III.1.

##  V Results

Property | Value
---|---
$h/R$ | $0.05$
$h_{\mathrm{w}}/h$ | $0.2$
$E_{\mathrm{0}}$ | $195\,\text{GPa}$
$E_{\mathrm{w}}$ | $12\,\text{MPa}$
$\nu$ | $0.28$
$\nu_{\mathrm{w}}$ | $0.4$
$\rho$ | $7700\,\text{kg}/\text{m}^{3}$
$\rho_{\mathrm{w}}$ | $1000\,\text{kg}/\text{m}^{3}$
Table 1: Material and geometric properties for the CMM and FEA. Figure 2:
Magnitude of the (a) radial displacement, and (b) tangential displacement of
the ring at the waveguide port locations for spatiotemporal parameters
$(\alpha,\bar{\omega}_{\mathrm{m}})=(0.2,0.016)$. The blue, orange, and green
curves denote the displacement magnitude at ports 1, 2, and 3, respectively.
Furthermore, the solid lines and open circles are the results from the CMM and
FEA, respectively. Figure 3: Absolute square of the harmonic modal amplitudes
for (a) $|\hat{w}^{p}_{m}|^{2}$ and (b) $|\hat{v}^{p}_{m}|^{2}$ normalized by
the absolute square of the maximum harmonic modal amplitude in decibels for
spatiotemporal parameters $(\alpha,\bar{\omega}_{\mathrm{m}})=(0.2,0.016)$ and
drive frequency $\bar{\omega}_{0}=0.21604$. Figure 4: Finite element
displacement solution for spatiotemporal parameters
$(\alpha,\bar{\omega}_{\mathrm{m}})=(0.2,0.016)$ at the drive frequency
$\bar{\omega}_{0}=0.21604$. The blue and green arrows denote the symmetric
($S_{0}$) and antisymmetric ($A_{0}$) Lamb wave, respectively, where the
magnitude of the arrows qualitatively designate the relative contributions of
each Lamb wave within the port waveguides. (a) Top-down view of the
circulator. An incident symmetric lamb wave is incident to the circulator in
port 1. (b) Zoom-in plot of port 3, which contains an outgoing antisymmetric
Lamb wave. (c) Zoom-in of port 2, which contains a dominant symmetric and low-
amplitude antisymmetric Lamb wave.

We now utilize both the CMM and FEA of the elastic wave circulator derived in
Sec. III and Sec. IV, respectively, to seek circulator and modulation designs
that yield a large degree of nonreciprocity at the output port waveguides. The
geometric and material parameters used for both models are shown in Table 1.
The ring is modeled using nominal material property values of steel, and the
material properties of the port waveguides are chosen to be representative of
rubber. This choice in the material property contrast results in a
$\bar{Z}=1\mathrm{e}{-4}$, which induces light damping due to the port
waveguides and therefore confines the nonreciprocal behavior to be close to
the resonance frequencies of the unmodulated circulator [25]. A unit input
force is prescribed for both models, i.e., $F_{0}=1$ in Eq. (5), and the
reference displacement is chosen to be $U=8\mathrm{e}{-7}$.

The parameter $\gamma$, which depends on ring geometry, can strongly influence
the accuracy of the predicted resonance frequency for the CMM. The definition
of $\gamma$ in Eq. (4) is only accurate for long-wavelength modes compared to
the circulator thickness, $\lambda/h\ll 1$. A better approximation to $\gamma$
can be sought by calibrating the CMM with FEA. This is accomplished by first
solving for the eigenvalues of the unmodulated circulator without port
waveguides and spatiotemporal modulation, i.e.
$(\alpha,\omega_{\mathrm{m}},\bar{Z})=(0,0,0)$, using FEA and fitting the
parameter $\gamma$ using the equation [30]

$\gamma=\sqrt{\frac{(\bar{\omega}_{n}^{2}-1)(\bar{\omega}_{n}^{2}-n^{2})-n^{2}}{n^{2}(\bar{\omega}_{n}^{2}-1)+n^{4}(\bar{\omega}_{n}^{2}-n^{2}+2)}},$
(33)

which is derived by inverting the dispersion relation of the ring for the
parameter $\gamma$. The parameters in Eq. (33) are the mode number $n$, and
$\bar{\omega}_{n}=\omega_{\mathrm{FEA,n}}R/c_{0}$, where
$\omega_{\mathrm{FEA,n}}$ is the resonance frequency for mode $n$ computed
with FEA.

In this work, we focus only on the mode $n=4$, which we found to be the mode
with the smallest resonance frequency displaying strong nonreciprocity for
small modulation parameters. This results in $\gamma=0.0149467$, which is used
for the remainder of this work. The normalized Young’s modulus of the ring is
modulated with the traveling wave form

$\bar{E}(\theta,t)=1+\alpha\cos(\theta-\omega_{\mathrm{m}}t),$ (34)

which results in $\hat{E}^{\pm 1}(\theta)=e^{\pm i\theta}/2$,
$\hat{E}^{|p|>1}=0$ in Eq. (6).

Figure 2 shows the drive frequency component ($p=0$ in Eqns. (11)-(12)) of the
radial and tangential displacement magnitudes of the ring at the port
waveguide locations computed with the CMM and FEA with prescribed modulation
parameters $(\alpha,\bar{\omega}_{\mathrm{m}})=(0.2,0.016)$ as a function of
the drive frequency. We first note that excellent agreement is obtained
between the two models. The degree of nonreciprocity increases as the drive
frequency approaches the value $\bar{\omega}_{0}=0.21604$, where the radial
displacement magnitude at port 3 exhibits a null while the radial displacement
magnitude at port 2 exhibits a maximum. The tangential displacement also
displays similar behavior, except a maximum in the displacement magnitude is
present at port 3 and a minimum in the displacement magnitude is present at
port 2. Note, however, that the degree of nonreciprocity for the tangential
displacement is not as large as the radial displacement. The effects of this
nonreciprocal ring displacement solution on the waveguide modes will be
discussed further below. The large degree of nonreciprocity shown in Fig. 2 is
further elucidated by investigating the magnitude of the modal amplitudes from
the CMM normalized by the magnitude of the maximum modal amplitude in
decibels, which is shown in Fig. 3. We find that the large degree of
nonreciprocity at the drive frequency ($p=0$) stems from the interference
pattern of the clockwise and counter-clockwise mode ($n=-4,4$), which is
qualitatively similar to previous works [25]. In addition, harmonic
frequencies are generated but are more than $18$ dB lower than the modes
generated at the drive frequency.

Figure 5: Magnitude of the radial and transverse displacements at the port
waveguide locations for various optimal modulation parameters. The modulation
parameters are given in parenthesis above each set of curves. The vertical
dashed line is the resonance frequency of the unmodulated circulator. The
displacements $|w|$ and $|v|$ are identical at ports 2 and 3 for the
unmodulated case. The fading of the curves is utilized to minimize the clutter
due to overlap of multiple curves on a single plot. Figure 6: Magnitude of the
radial and transverse displacements at the port waveguide locations for
various modulation frequencies with fixed $\alpha=0.1$. The vertical dashed
line is the resonance frequency of the unmodulated circulator.

Figure 4 shows the displacement solution of the finite element model at the
drive frequency $\bar{\omega}_{0}=0.21604$ for spatiotemporal modulation
parameters $(\alpha,\bar{\omega}_{\mathrm{m}})=(0.2,0.016)$. Note that
transverse displacements of the ring corresponds to an excitation of the
antisymmetric ($A_{0}$) Lamb wave in the waveguide ports, while radial
displacements of the ring induce the symmetric ($S_{0}$) Lamb wave in the
waveguide ports. Therefore, the nonreciprocity of the ring displayed in Fig. 2
indicates that the present circulator design acts as a nonreciprocal Lamb wave
converter, where the incident symmetric Lamb wave in port 1 is converted to an
antisymmetric Lamb wave in port 3 and a dominant symmetric Lamb wave in port
2, as shown in Fig. 4. While the presence of the null in the radial
displacement magnitude at port 3 in Fig. 2 indicates efficient nonreciprocal
conversion of the incident symmetric Lamb wave to an antisymmetric Lamb wave,
the same cannot be said about the tangential displacement magnitude. As a
consequence, port 2 will contain a combination of symmetric and antisymmetric
Lamb waves, although the magnitude of the antisymmetric Lamb wave is 17 dB
lower than the symmetric Lamb wave for this case. The further reduction of the
antisymmetric Lamb wave in port 2 could potentially be optimized by
considering different material properties of the ring and waveguides, as well
as three-dimensional effects such as port waveguide vertical thicknesses that
are smaller than the ring vertical thickness, which is left for future work.

We now seek to optimize the nonreciprocal response of the proposed elastic
wave circulator by fixing the dimensionless modulation amplitude $\alpha$ and
sweeping over the drive and modulation frequencies. The objective function for
the minimization problem is the ratio of the absolute square of the radial
displacement at port 3 over the radial displacement at port 2, i.e.
$|w(\theta_{3})/w(\theta_{2})|^{2}$. The optimization is carried out using the
minimize function in the scientific library SciPy for Python [33]. Figure 5
shows the results from the optimization for the modulation depth parameter
values of $\alpha=0.1$, $0.2$, and $0.3$. The unmodulated response is also
included as a reference. We find that as $\alpha$ is increased, smaller values
of $\omega_{\mathrm{m}}$ are needed to obtain a large nonreciprocal response.
In addition, increasing $\alpha$ corresponds to a downward shift in the drive
frequency at which the large nonreciprocal response appears. We also note that
the magnitude of the radial displacement at port 2 increases compared to the
unmodulated circulator when modulation is active, while the amplitude of the
tangential displacement at port 3 stays at nearly the same level.

Finally, the effect of the modulation frequency on the degree of
nonreciprocity is further studied by investigating the variation of the ring
displacement magnitude as a function of $\omega_{\mathrm{m}}$ for fixed
$\alpha=0.1$ in Fig. 6. We find that the drive frequency at which the large
nonreciprocal response appears is not sensitive to changes in the modulation
frequency, and is therefore only influenced by $\alpha$. As the modulation
frequency increases the degree of nonreciprocity between the port displacement
increases up to a critical value (bottom left panel in Fig. 6) before the
nonreciprocal performance decreases, suggesting that a unique, optimal value
of $\omega_{\mathrm{m}}$ may exist for a given $\alpha$.

##  VI Conclusion

We have developed a semi-analytical CMM and a FEA to study nonreciprocal
propagation in an elastic wave circulator in which the Young’s modulus is
modulated in a traveling wave fashion around a ring-shaped junction with three
semi-infinite elastic waveguide ports. The CMM is based on curved Euler-
Bernoulli beam theory and is valid for small curvatures and low frequencies of
excitation. In contrast, the FEA solves the exact elastodynamic equations and
incorporates a mode-matching technique that exactly accounts for outward-
propagating Lamb waves in the port waveguides. Both models can be applied to
nonreciprocal elastic systems that contain coupling between multiple resonant
frequencies as well as multiple wave polarizations.

For the cases considered in this work, we have demonstrated that this system
displays a large degree of nonreciprocity near the resonance frequencies of
the unmodulated system. In particular, the thin ring acts as a nonreciprocal
Lamb wave conversion device, in which an incident symmetric Lamb wave is
converted into outgoing symmetric and antisymmetric Lamb waves at the two
output port waveguides. The maximum degree of nonreciprocity and the
corresponding drive frequency are primarily influenced by the modulation
frequency and amplitude, respectively. We have found excellent agreement
between the CMM and FEA for the present case, although we expect the CMM to
become less accurate for cases of larger ring and port thicknesses, as well as
high frequency modes. Therefore, the FEA is needed to explore the entire
parameter space of elastic wave circulator designs. Further design
optimization, as well as finite thickness effects that remove the plane strain
assumptions in this work, should be further explored to find circulator
designs that yield a large degree of nonreciprocity. In addition, experimental
efforts to produce the desired stiffness modulation and measure the
nonreciprocal response of the circulator should be conducted to validate the
numerical results in this study.

###### Acknowledgements.

This work supported by National Science Foundation EFRI award No. 1641078 and
the Postdoctoral Fellowship Program at Applied Research Laboratories at The
University of Texas at Austin.

## References

  * Achenbach and Achenbach [2003] J. Achenbach and J. Achenbach,  _Reciprocity in elastodynamics_ (Cambridge University Press, 2003).
  * Strutt [1871] J. W. Strutt, Some general theorems relating to vibrations, Proceedings of the London Mathematical Society 1, 357 (1871).
  * Nassar  _et al._ [2020] H. Nassar, B. Yousefzadeh, R. Fleury, M. Ruzzene, A. Alù, C. Daraio, A. N. Norris, G. Huang, and M. R. Haberman, Nonreciprocity in acoustic and elastic materials, Nature Reviews Materials 5, 667 (2020).
  * Liang  _et al._ [2009] B. Liang, B. Yuan, and J.-c. Cheng, Acoustic diode: Rectification of acoustic energy flux in one-dimensional systems, Physical review letters 103, 104301 (2009).
  * Mojahed  _et al._ [2019] A. Mojahed, J. Bunyan, S. Tawfick, and A. F. Vakakis, Tunable acoustic nonreciprocity in strongly nonlinear waveguides with asymmetry, Physical Review Applied 12, 034033 (2019).
  * Moore  _et al._ [2018] K. J. Moore, J. Bunyan, S. Tawfick, O. V. Gendelman, S. Li, M. Leamy, and A. F. Vakakis, Nonreciprocity in the dynamics of coupled oscillators with nonlinearity, asymmetry, and scale hierarchy, Physical Review E 97, 012219 (2018).
  * Trainiti and Ruzzene [2016] G. Trainiti and M. Ruzzene, Non-reciprocal elastic wave propagation in spatiotemporal periodic structures, New Journal of Physics 18, 083047 (2016).
  * Nassar  _et al._ [2017] H. Nassar, X. Xu, A. Norris, and G. Huang, Modulated phononic crystals: Non-reciprocal wave propagation and willis materials, Journal of the Mechanics and Physics of Solids 101, 10 (2017).
  * Sugino  _et al._ [2020] C. Sugino, M. Ruzzene, and A. Erturk, Nonreciprocal piezoelectric metamaterial framework and circuit strategies, Physical Review B 102, 014304 (2020).
  * Shen  _et al._ [2019] C. Shen, X. Zhu, J. Li, and S. A. Cummer, Nonreciprocal acoustic transmission in space-time modulated coupled resonators, Physical Review B 100, 054302 (2019).
  * Riva  _et al._ [2020] E. Riva, M. Di Ronco, A. Elabd, G. Cazzulani, and F. Braghin, Non-reciprocal wave propagation in discretely modulated spatiotemporal plates, Journal of Sound and Vibration 471, 115186 (2020).
  * Oudich  _et al._ [2019] M. Oudich, Y. Deng, M. Tao, and Y. Jing, Space-time phononic crystals with anomalous topological edge states, Physical Review Research 1, 033069 (2019).
  * Goldsberry  _et al._ [2019] B. M. Goldsberry, S. P. Wallen, and M. R. Haberman, Non-reciprocal wave propagation in mechanically-modulated continuous elastic metamaterials, The Journal of the Acoustical Society of America 146, 782 (2019).
  * Wallen and Haberman [2019] S. P. Wallen and M. R. Haberman, Nonreciprocal wave phenomena in spring-mass chains with effective stiffness modulation induced by geometric nonlinearity, Physical Review E 99, 013001 (2019).
  * Popa and Cummer [2014] B.-I. Popa and S. A. Cummer, Non-reciprocal and highly nonlinear active acoustic metamaterials, Nature communications 5, 1 (2014).
  * Zhai  _et al._ [2019] Y. Zhai, H.-S. Kwon, and B.-I. Popa, Active willis metamaterials for ultracompact nonreciprocal linear acoustic devices, Physical Review B 99, 220301 (2019).
  * Chen  _et al._ [2020] Y. Chen, X. Li, G. Hu, M. R. Haberman, and G. Huang, An active mechanical willis meta-layer with asymmetric polarizabilities, Nature communications 11, 1 (2020).
  * Koutserimpas and Fleury [2018] T. T. Koutserimpas and R. Fleury, Nonreciprocal gain in non-hermitian time-floquet systems, Physical review letters 120, 087401 (2018).
  * Rosa and Ruzzene [2020] M. I. Rosa and M. Ruzzene, Dynamics and topology of non-hermitian elastic lattices with non-local feedback control interactions, New Journal of Physics 22, 053004 (2020).
  * Caloz  _et al._ [2018] C. Caloz, A. Alu, S. Tretyakov, D. Sounas, K. Achouri, and Z.-L. Deck-Léger, Electromagnetic nonreciprocity, Physical Review Applied 10, 047001 (2018).
  * Wiederhold  _et al._ [2019] C. P. Wiederhold, D. L. Sounas, and A. Alù, Nonreciprocal acoustic propagation and leaky-wave radiation in a waveguide with flow, The Journal of the Acoustical Society of America 146, 802 (2019).
  * Darabi  _et al._ [2019] A. Darabi, L. Fang, A. Mojahed, M. D. Fronk, A. F. Vakakis, and M. J. Leamy, Broadband passive nonlinear acoustic diode, Physical Review B 99, 214305 (2019).
  * Zangeneh-Nejad and Fleury [2018] F. Zangeneh-Nejad and R. Fleury, Doppler-based acoustic gyrator, Applied Sciences 8, 1083 (2018).
  * Fleury  _et al._ [2014] R. Fleury, D. L. Sounas, C. F. Sieck, M. R. Haberman, and A. Alù, Sound isolation and giant linear nonreciprocity in a compact acoustic circulator, Science 343, 516 (2014).
  * Fleury  _et al._ [2015a] R. Fleury, D. L. Sounas, and A. Alù, Subwavelength ultrasonic circulator based on spatiotemporal modulation, Physical Review B 91, 174306 (2015a).
  * Fleury  _et al._ [2015b] R. Fleury, D. Sounas, M. R. Haberman, and A. Alu, Nonreciprocal acoustics, Acoustics Today 11, 14 (2015b).
  * Vila  _et al._ [2017] J. Vila, R. K. Pal, M. Ruzzene, and G. Trainiti, A Bloch-based procedure for dispersion analysis of lattices with periodic time-varying properties, Journal of Sound and Vibration 406, 363 (2017).
  * Beli  _et al._ [2018] D. Beli, P. B. Silva, and J. R. de França Arruda, Mechanical circulator for elastic waves by using the nonreciprocity of flexible rotating rings, Mechanical Systems and Signal Processing 98, 1077 (2018).
  * Goldsberry  _et al._ [2020] B. M. Goldsberry, S. P. Wallen, and M. R. Haberman, Nonreciprocal vibrations of finite elastic structures with spatiotemporally modulated material properties, Physical Review B 102, 014312 (2020).
  * Graff [2012] K. F. Graff,  _Wave motion in elastic solids_ (Courier Corporation, 2012).
  * Heckl [1961] M. A. Heckl,  _Compendium of impedance formulas_ , Tech. Rep. (BOLT BERANEK AND NEWMAN INC CAMBRIDGE MA, 1961).
  * Astley [1996] R. Astley, FE mode-matching schemes for the exterior Helmholtz problem and their relationship to the FE-DtN approach, Communications in numerical methods in engineering 12, 257 (1996).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Scipy 1.0: fundamental algorithms for scientific computing in Python, Nature Methods 17, 261 (2020).