Non-reciprocity in discrete and continuous lattice systems via mechanicalmodulation
Benjamin M. Goldsberry(1) , Samuel P. Wallen(1) , Michael R. Haberman(1)(1) Walker Department of Mechanical Engineering and Applied Research Laboratories, The University of Texas at Austin, Austin,

TX, USA

Abstract
Materials with time- and space-dependent effective properties are of significant interest as a means to introduce
non-reciprocal wave propagation in engineered material systems. Materials displaying this unique behavior havesignificant potential as direction- dependent devices for mechanical wave sensing and transmitting capabilities.
They also open up the possibility for the construction of novel materials or structures that isolate unwanted
vibration. Recent analytical models of spring-mass chains have shown that external application of a nonlinear
mechanical deformation, when applied on time scales that are slow compared to the characteristic times of
propagating linear elastic waves, may induce non-reciprocity via changes in the apparent elastic modulus for
perturbations around that deformation [Wallen et al., Phys.Rev. E., 99, 013001, (2019)]. This work investigates
direct numerical simulation of discrete lattices and a finite element approach for continuous elastic systems
to simulate elastic wave propagation in mechanically-modulated metamaterials. Analysis of the degree of nonreciprocity when varying sub-wavelength geometry or geometric modulations will be shown. Of specific interestis the use of FEM to investigate non- reciprocity in elastic lattices consisting of unit cells with varied geometricasymmetry and more general forms of mechanical modulation.
Keywords: Elastic waves, Metamaterials, Non-reciprocity

1

INTRODUCTION

Acoustic systems that break reciprocity have recently received great attention due to their potential to increase
control over wave propagation and to create devices that make use of acoustic or elastic waves. Tools forthe modeling, design, and analysis of media that generates non-reciprocal wave propagation have strong potential for far-reaching impacts on engineering applications, including more efficient acoustic communicationand imaging devices, analogue signal processing, and vibration isolation. One means to break reciprocity isthrough spatiotemporal modulation of material properties [1]. Recent analytical models of spring-mass chains
have shown that external application of a nonlinear mechanical deformation, when applied on time scales thatare slow compared to the characteristic times of propagating linear elastic waves, may induce non-reciprocity
via changes in the linearized local elastic modulus for perturbations around that deformation [2].
In a recent work, we investigated a continuous tunable mechanical metamaterial structure as a platform toachieve non-reciprocity [3]. A computational approach based on the finite element method was derived andimplemented due to the complex unit cell geometry used to obtain large effective dynamic property modulations.
In this work, we revisit the derivation of our computational approach by characterizing wave propagation in a
medium with space-only and time-only modulation of the material modulus, and determining how these concepts
generalize to the case of spatiotemporal modulation of the material modulus. Our analysis is implemented in a
one-dimensional setting to compare with known exact solutions. Longitduinal wave propagation in a bar is used
as an example case for each study.

5580

2

WAVE PROPAGATION IN A MODULATED MEDIUM

The partial differential equation that describes wave propagation in a modulated medium is the one-dimensional
wave equation of the form


∂
∂u
∂ 2u
a(x,t)− b 2 = 0, −∞ < x < ∞, t ≥ 0,
(1)
∂x
∂x
∂t
with coefficient a(x,t) > 0 that is a non-constant function of time and space, and coefficient b > 0 that is assumed to be constant. Equation (1) can be used to represent various well-studied physical systems, such as
longitudinal wave propagation in a bar. The dependent variable u(x,t) is then the longitudinal particle displacement of the bar, a(x,t) = E(x,t) is the Young’s modulus, and b = ρ0 is the mass per unit volume of thebar.
It is assumed that there is a reference modulation length, λm , in which a(x,t) is periodic in space for all time,a(x + λm ,t) = a(x,t),(2)

and a reference time Tm with angular frequency ωm = 2π/Tm in which a(x,t) is periodic in time for all locationsin space, i.e.
a(x,t + Tm ) = a(x,t).
(3)
Therefore, Bloch-Floquet wave theory can be utilized to derive differential equations, which can be solved forthe allowable traveling-wave modes and their respective dispersion behavior, or frequency-wavenumber spectrum[4]. Solution procedures for Eq. (1) with space-only, time-only, and space-time modulation of a(x,t) are now
derived that are computationally feasible and yield accurate calculations of the Bloch wave modes and dispersion
relations.
2.1 Space-dependent modulus
Consider a modulus that is a function of space only, a(x,t) → a(x). This case physically represents a periodic
heterogeneous medium, which is a common acoustic metamaterial structure and is frequently referred to as aphononic crystal [5, 6]. The traveling wave solution has the following formu(x,t) = Cû(x)ei(kx−ωt) ,(4)

where C is the amplitude constant, û(x) is the Bloch wave mode that is a periodic function on the unit celldomain of the phononic crystal, k is the Bloch wavenumber, and ω is the Bloch wave angular frequency. The
finite element method has been used in previous works to successfully solve Bloch wave problems of this type,which leads to a generalized eigenvalue problem of the form [7]K(k) − ω 2 M û = 0,(5)where K is the stiffness matrix which is a function of the Bloch wavenumber, M is the mass matrix, and û
is the Bloch wave mode that is discretized in space on a computational mesh. The dispersion relation ω(k) issolved by sweeping k in the first Brillouin zone and solving Eq. (5) for ω.
As an example case, consider longitudinal wave propagation in a bar composed of two periodically-alternating
material layers with identical density b = ρ0 . The first layer has Young’s modulus E1 and thickness L1 , and thesecond layer has Young’s modulus E2 and thickness L2 . The periodic interval has length L = L1 + L2 = λm andunit cell defined at −L/2 ≤ x ≤ L/2. Assume, without loss of generality of the method, that the two layers havethe same thickness, L1 = L2 . The modulus a(x) = E(x) is therefore defined as the piecewise function(
E1 , if − L2 ≤ x < 0,
E(x) =(6)
E2 , if 0 ≤ x < L2 .

5581

(a) 10

(b)

8

ω

6
4
2
0

3

2

1

0

kL/π

1

2

3

1.0

0.5

0.0

kL/π

0.5

1.0

Figure 1. (a) Finite element computation of the dispersion relation for propagating longitudinal waves in a bar
with periodic alternating layers. Since the dispersion relation is periodic, the complete dispersion relation canbe obtained by restricting k to lie within the first Brillouin zone (black dashed lines). (b) Comparison of thefinite element computation in the first Brillouin zone from (a) with the exact dispersion relation (orange lines),
Eq. (7).
The exact dispersion relation in this case is the equation [4]





 
 

ω
ω
ω
ω
1 z1 z2cos (kL) = cos
L1 cosL2 −+sin
L1 sin
L2 ,
(7)
c1
c2
2 z2 z1c1c2pwhere c1,2 = E1,2 /ρ0 is the sound speed in each layer, and z1,2 = ρ0 c1,2 is the impedance of each layer. Figure1(a) shows the resulting eigenfrequencies from the finite element calculation using the properties E1 = 1, E2 = 2,
ρ0 = 1, L = 1, and Fig. 1(b) compares the finite element results with the exact dispersion relation provided by
Eq. (7). Excellent agreement is obtained between the exact dispersion relation and the finite element simulation.
As a consequence of the Bloch-Floquet theory, the dispersion relation is periodic with wavenumber periodicity
ω(k + kL ) = ω(k), where kL = 2π/L. The complete dispersion curve can therefore be obtained by limiting k
to the first Brillouin zone, which is the range k ∈ [−π/L, π/L] and is denoted by the vertical dashed lines inFig. 1(a) [4]. In addition, the periodicity of the layers gives rise to Bragg scattering, in which the scatteredwaves caused by the impedance mismatch between the layers destructively interferes with the traveling wave
in certain frequency bands [8, 9]. This interference creates bandgaps in the dispersion relation, which prohibits
traveling waves with frequencies within the gap. The widths of the bandgaps are a function of the contrast of
the material moduli ∆E = E2 − E1 .
2.2 Time-dependent modulus
Now consider the case where the modulus is a periodic function of time, a(x,t) → a(t). For this case, Eq. (1)reduces to
∂ 2u
∂ 2u
a(t) 2 − b 2 = 0.
(8)
∂x
∂t
The traveling wave solution to Eq. (8) has the following formu(x,t) = CΦ(t)ei(kx−ωt) .

5582

(9)

(b)

0
5

dB

10
15

ω

(a) 2.00
1.75
1.50
1.25
1.00
0.75
0.50
0.25
0.00

20
25

1.0

0.5

0.0

k

0.5

1.0

1.0

0.5

0.0

k

0.5

1.0

30

Figure 2. (a) Comparison of the frequency-wavenumber spectrum using Eq. (13) (open circles) with the exact
dispersion solution, Eq. (17). (b) Dispersion relation from (a) weighted by the magnitude of the fundamental
plane wave component in decibels, 20 log10 (||Φ̂0 ||/||Φ̂||).
Note that Eq. (9) has the same form as Eq. (4), except the Bloch function in Eq. (9) is now a periodic functionof time. Since a(t) is assumed to be periodic with frequency ωm = 2π/Tm , a(t) can be written as a Fourier
series in time,Na(t) =∑ ân e−inωmt ,

(10)

n=−N

where it is assumed that the function a(t) can be well-represented with 2N + 1 Fourier terms. The coefficientsân are determined by the Fourier integral
ân =

ωm
2π

Z π/ωma(t)einωm t .

(11)

−π/ωm

Likewise, Φ(t) can also be expanded as a Fourier series in time,NΦ(t) =∑ Φ̂m e−imωmt .

(12)

m=−N

Substitution of Eqns. (10) and (12) into Eq. (8) and utilizing orthogonality of the Fourier series yields aquadratic eigenvalue problem,ω 2 MΦ̂ + ωCΦ̂ + KΦ̂ = 0,(13)where Φ̂ = [Φ̂−N , Φ̂−N+1 , ..., Φ̂N ]T , M is the identity matrix, C depends on ωm , and K depends on â, ωm , andk. For example, if N = 1 then the matrix elements in Eq. (13) are

 

 



Φ̂−1
â0 − ωm2 â−10
1 0 0
−2ωm 002kω 2 0 1 0 + ω  0â0â−1   Φ̂0  = 0.
00  −  â1
(14)
b
0 0 1
0
0 2ωm
0
â1 â0 − ωm2
Φ̂1
As an example, consider the propagation of longitudinal waves in a bar with constant density b = ρ0 and atime-varying Young’s modulus a(t) → E(t), where E(t) varies as a square wave with minimum value E1 andmaximum value E2 , which is written as
E(t) = E1 + ∆EH[cos(ωmt)],

(15)

5583

where H(·) is the Heaviside function, and ∆E = E2 − E1 . The Fourier coefficients in Eq. (10) are thenÊn = E2 δn0 −∆Enπsinc( ),
2
2

(16)

where sinc(x) = sin(x)/x and δn0 is the Kronecker delta. The exact dispersion relation in this case is theequation [11]
1 c1 T1 c2 T2+sin(kc1 T1 ) sin(kc2 T2 ),(17)cos(ωT ) = cos(kc1 T1 ) cos(kc2 T2 ) −2 c2 T2 c1 T1pwhere T1 = T2 = π/ωm , T = T1 + T2 , and c1,2 = E1,2 /ρ0 [4]. This particular case has a similar dispersion
form to the space-only dispersion relation, Eq. (7), except it is a function of time. Figure 2(a) compares thedispersion relation computed using Eq. (13) with parameters N = 5, ωm = 1, E1 = 1, E2 = 3, ρ0 = 1 and two
branches of the exact dispersion relation, Eq. (17). Note that the dispersion curves in Fig. 2(a) are periodic infrequency with respect to the modulation frequency ωm , which is the analogue to wavenumber periodicity inthe case of a spatially-varying modulus. The computed dispersion relation using Eq. (13) agrees well with theexact dispersion curves, and better agreement is obtained, particularly in the higher frequencies, by increasing
the number of Fourier terms N. It is observed from Fig. 2(a) that this system exhibits wavenumber gaps, which
are the analog of bandgaps in frequency for the the case of a spatially-varying modulus. However, propagatingwaves within the wavenumber bandgaps are unstable and exhibit infinite growth in amplitude in the losslesscase [10].
The elements from the computed eigenvector Φ̂ are the amplitudes of the plane waves in Eq. (9). One also notesthat the total energy of the propagating waves, which is related to the squared magnitude of the plane wave amplitudes, are not evenly distributed across each frequency. This can be represented on a frequency-wavenumber
spectrum by weighting each curve by the magnitude of the fundamental plane wave amplitude ||Φ̂0 ||2 . Figure2(b) depicts the same results as Fig. 2(a), but each curve is now assigned a color that represents the normalized
magnitude of the fundamental (m = 0) plane wave component in decibels, specifically 20 log10 (||Φ̂0 ||/||Φ̂||). The
main dispersion branch is now easily identifiable, as well as the branches that interact with the main branchnear the wavenumber gaps.
2.3 Space- and time-dependent modulus
Finally, consider the case where now a(x,t) is a periodic function of time and space. The solution techniques
discussed for the time-only and space-only modulus cases can be leveraged to solve this general case. First,
a(x,t) is expanded as a Fourier series in time,Na(x,t) =∑ ân (x)e−inωmt ,

(18)

n=−N

where ân (x) are the harmonic amplitudes and are periodic functions of space. As in the previous sections, a
Bloch wave solution is assumed, where the Bloch mode U(x,t) is now a periodic function of time and space,u(x,t) = ei(kx−ωt)U(x,t).

(19)

Therefore, U(x,t) can also be expanded as a Fourier series in time with harmonic amplitudes that are a periodic
function of spaceN

U(x,t) =∑ ûm (x)e−imωmt .

(20)

m=−N

Equation (1) is solved by substituting Eqns. (18)-(20) into Eq. (1), utilizing orthogonality of the Fourier series,
and using the finite element method to discretize the resulting coupled ordinary differential equations that are afunction of space into a quadratic eigenvalue problem
ω 2 M̄u + ω C̄u + K̄u = 0,

5584

(21)

(a) 5

(b)

0
5

4

10

1
-1

dB

p=
0

-1

p=

0

1
0

p=

20

p=

15

p=

2

1

ω

p=

3

25

1.0

0.5

0.0

kL/π

0.5

1.0

1.0

0.5

0.0

kL/π

0.5

1.0

30

Figure 3. (a) Finite element computation of the dispersion relation for the case of translating material layerswith modulation speed cm = 0.159 (ωm = 1). Due to the bias in the +x direction, the dispersion curves are notsymmetric about k = 0, i.e. reciprocity is broken. (b) The dispersion curves shown in (a) but with each pointcolored by the magnitude of the fundamental component in decibels, 20 log10 (||u0 ||/||u||). The dashed line is
the first branch of the exact dispersion relation for this case, which is obtained using Eq. (26) and Eq. (7).
where u = [û−N , ..., ûN ] and û is discretized in space on a computational mesh, M̄ is the global assembly of thefinite element mass matrices, C̄ is the global assembly of the finite element matrices proportional to ω, and K̄
is the global assembly of the finite element stiffness matrices.
As an example, consider the case of longitudinal wave propagation in a bar with constant density b = ρ0composed of periodic alternating layers from Sec. 2.1, now translating in the +x direction at speed cm = ωm /km ,where km = 2π/λm is the modulation wavenumber. Therefore, the Young’s modulus a(x,t) → E(x,t) has thefollowing functional form, which was recently studied using a different approach by Trainiti and Ruzzene [1]E(x,t) = E1 + ∆EH [cos (km x − ωmt)] ,(22)

which is a generalization of the time-varying Young’s modulus, Eq. (15). The Fourier coefficients in Eq. (18)
are thus given by [1]
h πn i∆EÊn = E2 δn0 −sinceinkm x .
(23)
2
2
The exact dispersion relation in this case can be obtained by shearing the dispersion relation for the case of
spatial variation of the Young’s modulus, Eq. (7). This is accomplished by replacing k with k̂, ω with ω̂, andc1 , c2 with ĉ1 , ĉ2 , where
c2m
,
c1
c2
ĉ2 = c2 − m ,c2ĉ1 = c1 −

(24)
(25)

and solving Eq. (7) with these transformed variables [11]. The variables (k̂, ω̂) can be transformed back to(k, ω) with the following transformation,k = k̂ +ω̂cm L1ω̂cm L2+,
L(c21 − c2m ) L(c22 − c2m )ω = ω̂ + cm k̂ +ω̂c2m L1ω̂c2m L2+.
L(c21 − c2m ) L(c22 − c2m )

(26)
(27)

5585

(a) 5

(b)

0

4
3

ω

6

2
1

12

0
(c) 5

(d)
18

4

ω

3
24

2
1
0

1.0

0.5

0.0

kL/π

0.5

1.0 1.0

0.5

0.0

kL/π

0.5

1.0

30

Figure 4. Dispersion relations for the translating layers case as a function of the modulation speed cm . (a) No
translation, cm = 0, (b) cm = 0.05, (c) cm = 0.1, (d) cm = 0.3.
Figure 3(a) shows the dispersion calculation from the finite element model using the parameters E1 = 1, E2 =1.5, ρ0 = 1, L = 1, P = 2, and ωm = 1. As in the time-only modulation case, the entire frequency spectrum
can be generated by translating the fundamental (p = 0) branch by ±pωm . A bias in space-time has now
been introduced in the +x direction, which breaks time-reversal symmetry, and the dispersion relation is nolonger symmetric about the k = 0 axis. Therefore, reciprocity is broken. Figure 3(b) shows the same dispersion
curves as in Fig. 3(a), but with each point assigned a color determined by the magnitude of the fundamental
component in decibels, specifically 20 log10 (||û0 ||/||u||). As in the time-only modulation case, it is evident fromthis plot that the harmonic amplitudes for each plane wave in Eq. (20) are not equal. This representation now
shows what appears to be directional bandgaps in the −x direction that do not exist in the +x direction. Thesebandgaps are not identical to bandgaps in a system with only a spatially-varying modulus, since stable modes
do exist within the bandgap. It can be seen from Fig. 3(b) that the branch with the largest magnitude agrees
with the exact dispersion relation computed using the transformed variables, Eq. (26) and Eq. (7).
The frequency ranges of the directional bandgaps can be tuned with the speed of the modulation, cm . Figure
4 depicts the changes in the dispersion relation as cm is increased. For low modulation speed, mode splitting
occurs and the spatial-only dispersion solution splits into the ±ωm curves [12]. For a modulation speed of cm =0.1, the band gap in the +x direction closes, and the band gap in the −x direction splits into two directionalband gaps, where the lowest directional band gap translates down as the modulation speed is increased, and thehigher directional band gap translates up in frequency as the modulation speed is increased. The width of theband gap can be tuned in the same way as the spatial-only modulation case, by increasing the contrast (∆E)
between layers.

5586

3

SUMMARY

In this work, we have derived a computational approach based on the finite element method to study wave
propagation in modulated media. The effects of space-only and time-only modulation of the material modulus
on the propagation of Bloch waves in an elastic bar were explored, and the computational techniques used forthese cases were combined to study the case of a space- and time-varying material modulus. It is straightforward
to generalize this approach to the general elastodynamic case to study non-reciprocal wave propagation in amechanically-modulated medium [3].

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[1] Trainiti, G.; Ruzzene, M. Non-reciprocal elastic wave propagation in spatiotemporal periodic structures.
New Journal of Physics, Vol 18 (8), 2016, 083047.
[2] Wallen, S. P.; Haberman, M. R. Nonreciprocal wave propagation in spring-mass chains with effective stiffness modulation induced by geometric nonlinearity. Physical Review E, Vol 99, 2019, 013001.
[3] Goldsberry, B. M.; Wallen, S. P.; Haberman, M. R. Non-reciprocal wave propagation in mechanicallymodulated continuous elastic metamaterials. Journal of the Acoustical Society of America, In press.
[4] Brillouin, L. Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices. Courier Corporation, 2003.
[5] Hussein, M. I.; Leamy, M. J.; Ruzzene, M. Dynamics of phononic materials and structures: Historicalorigins, recent progress, and future outlook. Applied Mechanics Reviews, Vol 66(4), 2014, 040802.
[6] Adams, S. D.; Craster, R. V.; Guenneau, S. Bloch waves in periodic multi-layered acoustic waveguides.
Proceedings of the Royal Society A, Vol 464 (2098), 2008, 2669–2692.
[7] Langlet, P.; Hladky-Hennion, A.; Decarpigny, J. Analysis of the propagation of plane acoustic waves inpassive periodic materials using the finite element method. Journal of the Acoustical Society of America,
Vol 98 (5), 1995, 2792–2800.
[8] Degraeve, S.; Granger, C.; Dubus, B.; Vasseur, J. O.; Pham Thi, M.; Hladky-Hennion, A. C. Bragg bandgaps tunability in an homogeneous piezoelectric rod with periodic electrical boundary conditions. Journal ofApplied Physics, Vol 115 (19), 2014, 194508.
[9] Deymier, P. A. Acoustic Metamaterials and Phononic Crystals. Vol 173, Springer Science & BusinessMedia, 2013.
[10] Torrent, D.; Parnell, W. J.; Norris, A. N. Loss compensation in time-dependent elastic metamaterials.
Physial Review B, Vol 97 (1), 2018, 014105.
[11] Nassar, H.; Xu, X.; Norris, A. N.; Huang, G. Modulated phononic crystals: Non-reciprocal wave propagation and Willis materials. Journal of the Mechanics and Physics of Solids, Vol 101, 2017, 10–29.
[12] Fleury, R.; Sounas, D. L.; Sieck, C. F.; Haberman, M. R.; Alù, A. Sound Isolation and Giant Linear
Nonreciprocity in a Compact Acoustic Circulator. Science. Vol 343 (6170), 2014, 516–519.

5587