<figure><img src="image_1.png" /></figure>

HAL
open science

# Wave propagation in one-dimensional nonlinear acoustic metamaterials

Xin Fang, Jihong Wen, Bernard Bonello, Jianfei Yin, Dianlong Yu

► To cite this version:

Xin Fang, Jihong Wen, Bernard Bonello, Jianfei Yin, Dianlong Yu. Wave propagation in one-dimensional nonlinear acoustic metamaterials. New Journal of Physics, 2017, 19, pp.053007. 10.1088/1367-2630/aa6d49. hal-01534470

HAL Id: hal-01534470

https://hal.sorbonne-universite.fr/hal-01534470

Submitted on 7 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Wave propagation in one-dimensional nonlinear acoustic metamaterials

This content has been downloaded from IOPscience. Please scroll down to see the full text.

2017 New J. Phys. 19 053007

(http://iopscience.iop.org/1367-2630/19/5/053007)

View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 134.157.80.157

This content was downloaded on 07/06/2017 at 14:54

Please note that terms and conditions apply.

You may also be interested in:

Edge waves in plates with resonators: an elastic analogue of the quantum valley Hall effect

Raj Kumar Pal and Massimo Ruzzene

Electron-phonon metamaterial featuring nonlinear tri-interleaved piezoelectric topologies and its
application in low-frequency vibration control

Bin Bao, Daniel Guyomar and Mickaël Lallart

Dynamics of a hyperbolic system that applies at the onset of the oscillatory instability

Carlos Martel and José M Vega

The Lorenz-84 climate model

Henk Broer, Carles Simó and Renato Vitolo

Hysteresis loops and multi-stability: From periodic orbits to chaotic dynamics (and back) in
diatomic granular crystals

C. Hoogeboom, Y. Man, N. Boechler et al.

Non-reciprocal elastic wave propagation in spatiotemporal periodic structures

G Trainiti and M Ruzzene

Hopf bifurcation with the symmetry of the square

J W Swift

Spatio-temporal phenomena in complex systems with time delays

Serhiy Yanchuk and Giovanni Giacomelli

Noise, chaos, and the Josephson voltage standard

R L Kautz
**New Journal of Physics**

The open access journal at the forefront of physics

Deutsche Physikalische Gesellschaft DPG

IOP Institute of Physics

Published in partnership
with: Deutsche Physikalische
Gesellschaft and the Institute
of Physics

<figure><img src="image_3.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

8 February 2017

REVISED

7 April 2017

ACCEPTED FOR PUBLICATION

13 April 2017

PUBLISHED

16 May 2017

Original content from this work may be used under the terms of the [Creative Commons Attribution 3.0 licence](https://creativecommons.org/licenses/by/3.0/).

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<figure><img src="image_4.png" /></figure>

PAPER

# Wave propagation in one-dimensional nonlinear acoustic metamaterials

Xin Fang<fnref n="1" />, Jihong Wen<fnref n="1" />, Bernard Bonello<fnref n="2" />, Jianfei Yin<fnref n="1" /> and Dianlong Yu<fnref n="1" />

E-mail: wenjihong@vip.sina.com and nmhsyjf@hotmail.com

**Keywords:** nonlinear acoustic metamaterial, wave propagation, band structure, nonlinear dynamics, chaos

Abstract

The propagation of waves in nonlinear acoustic metamaterial (NAM) is fundamentally different from that in conventional linear ones. In this article we consider two one-dimensional (1D) NAM systems featuring respectively a diatomic and a tetratomic meta unit-cell. We investigate the attenuation of waves, band structures, and bifurcations to demonstrate novel nonlinear effects, which can significantly expand the bandwidth for elastic wave suppression and cause nonlinear wave phenomena. The harmonic averaging approach, continuation algorithm, and Lyapunov exponents (LEs) are combined to study the frequency responses, nonlinear modes, bifurcations of periodic solutions, and chaos. The nonlinear resonances are studied, and the influence of damping on hyperchaotic attractors is evaluated. Moreover, a ‘quantum’ behavior is found between the low-energy and high-energy orbits. This work provides a theoretical base for furthering understandings and applications of NAMs.

**1. Introduction**

Acoustic metamaterials [1–6] (AMs) are artificial media that gain their properties from structure rather than composition. They generally feature built-in resonators with subwavelength dimensions, rendering the concepts of effective mass density and effective elastic constants relevant for characterizing them. Because of their ability to manipulate the propagation of sound, AMs have become a hot topic during the last decade [6]. By now, most studies [7–12] are dealing with linear AMs (LAMs), which are based on a locally resonant (LR) mechanism [3] to yield negative indexes, superlensing effects [12], or wave guiding. Different from the Bragg mechanism [13], LR paves the way to manipulate low-frequency waves with small ‘meta-atoms’ [12]. However, the LR bandgap is generally narrow [6], and heavy ‘atoms’ are necessary to enlarge this bandgap.

The concept of nonlinear metamaterial was first introduced in electromagnetism to investigate photonic metamaterials in 2003 [14]. Nonlinear responses can be constructed on purpose, such as tenability, electromagnetically induced transparency, plasmonics, active media, etc [15, 16]. In contrast, it is only recently that nonlinear acoustic metamaterials (NAMs) have emerged.

The presence of nonlinear media in linear phononic crystals can be utilized to realize interesting devices, such as acoustic diodes [17, 18]. Compared with NAMs, studies on one-dimensional (1D) nonlinear periodic structures [19] (NPSs) have a longer history that can date back to 1955 studies on the Fermi–Pasta–Ulam (FPU) models [20, 21] that promoted nonlinear physics. These nonlinear periodic chains are suitable for modeling a number of physical systems arising in different scientific contexts [22, 23]. Based on the FPU model, Nesterenko [24] experimentally highlighted solitons in the granular crystal (a kind of NPS) interacting in a nonlinear way through Hertzian contact [25, 26]. Moreover, a bistable lattice (another NPS) was recently studied [27]. Currently, many new physical properties that are different from linear phononic crystals are found in NPSs, such as unidirectional transition [27], discrete breathers [23, 28], waves coupling [29], subharmonic frequencies [30, 31], soliton waves [32, 33], and surface waves [34]. Both simulations [22, 23] and experiments [35, 36] demonstrate that bifurcations

<footnote n="1">¹ Laboratory of Science and Technology on Integrated Logistics Support, National University of Defense Technology, Changsha 410073, Hunan, People's Republic of China</footnote>

<footnote n="2">² Institut des NanoSciences de Paris (INSP-UMR CNRS 7588), Université Pierre et Marie Curie (box 840) 4, place Jussieu, F-75252 Paris Cedex 05, France</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1. Model of the diatomic NAM.</figcaption></figure>

in NPSs would be relevant to those properties. Acoustic devices such as diodes [36] and lenses [37] can be built upon them. The strong NPS has been proved to increase the velocity of sound and therefore the acoustic impedance [38]. Moreover, the bandgap properties in NPSs attract much attention [39–44]. The perturbation approach and the harmonic balance method (HBM) are adopted [40–42] to study the amplitude-dependent dispersions, stop band properties, and wave beaming in granular crystals; and experimental works highlighted the role played by the critical amplitude in energy transmission [43] and bifurcation-induced bandgap reconfiguration [44] in NPSs. Actually, the granular crystals are suitable for ultrasonic applications; it is hard to consider them at the low-frequency regime because of the high contact stiffness they inherently feature.

Because of their promising applications, AMs with both low-frequency and broadband properties attract much attention. However, the mechanisms for both properties to occur simultaneously are difficult to realize. LAMs consist of linear ‘meta-atoms’, but when this meta-atom becomes nonlinear in NAMs, wave propagation properties show different patterns. In our recent works [45, 46], the wave propagation in diatomic and tetratomic NAMs are analyzed using the homotopy analysis method, and we found that the chaotic bands resulting from bifurcations can significantly enlarge the width of the wave-suppressing bands. This finding demonstrates that chaos is a novel and promising mechanism to simultaneously achieve low-frequency and broadband in both mono-bandgap NAMs and multi-bandgap NAMs; this finding also reiterates still that a strong nonlinearity is beneficial to expand the bandwidth by several times.

However, there are many phenomena arising in NAMs that have not yet been fully explained nor demonstrated. For example, why are the responses in the first passband similar to those observed with LAMs? Under which conditions can the elastic energy propagate in the bandgap? When will the wave be amplified by chaos? In the tetratomic system, why does the nonlinearity have a larger influence on the nonlinear LR bandgaps than it has for Bragg gaps? Furthermore, some problems remain unexplored: for example, the influence of nonlinearity on the structure of bands, the features of chaos and their differences, the jumps, the nonlinear resonances and their stabilities, and the influence of damping on chaos.

In this paper we attempt to answer these questions with the help of frequency response analysis, the bifurcation theories, the Lyapunov exponents (LEs), and the fractal dimensions.

## 2. Models

### 2.1. Diatomic model

The 1D damped diatomic NAM model is illustrated in figure 1. The nonlinear meta-atom of this periodic model consists of a linear oscillator *m* and a damped nonlinear Duffing oscillator *m*<sub>r</sub> with cubic nonlinear stiffness *k*<sub>1</sub>Δ + *k*<sub>2</sub>Δ<sup>3</sup>, where Δ denotes the relative displacement, and *k*<sub>1</sub> (*k*<sub>2</sub>) symbolizes the linear (nonlinear) stiffness. In addition, a linear viscous damping term *c* · Δ̇ is taken into consideration.

Defining *x*<sub>n</sub> and *y*<sub>n</sub> as displacements of the linear and nonlinear oscillators in the *n*th cell, respectively; the differential equation for the *n*th cell reads

$$ \begin{cases} \ddot{x}_n = \beta_0(x_{n+1} + x_{n-1} - 2x_n) + \eta[\zeta(\dot{y}_n - \dot{x}_n) + \beta_1(y_n - x_n) + \beta_2(y_n - x_n)^3] \\ \ddot{y}_n = -\zeta(\dot{y}_n - \dot{x}_n) - \beta_1(y_n - x_n) - \beta_2(y_n - x_n)^3 \end{cases} \quad (1) $$

where the dot denotes the derivation with respect to time *t*. The definitions of the generalized parameters are as follows: mass ratio *η* = *m*<sub>r</sub>/*m*, stiffness ratio *β*<sub>0</sub> = *k*<sub>0</sub>/*m* = *ω*<sub>s</sub><sup>2</sup>, *β*<sub>1</sub> = *k*<sub>1</sub>/*m*<sub>r</sub> = *ω*<sub>r</sub><sup>2</sup>, *β*<sub>2</sub> = *k*<sub>2</sub>/*m*<sub>r</sub>. The generalized frequency is Ω = *ω*/*ω*<sub>s</sub>, and the damping coefficient is *ζ* = *c*/*m*<sub>r</sub>. The strength factor of nonlinearity is defined as *σ* = 3*β*<sub>2</sub>*A*<sub>0</sub><sup>2</sup>/*β*<sub>1</sub>, where *A*<sub>0</sub> stands for the amplitude of the incident wave *u*<sub>0</sub> = *A*<sub>0</sub>sin(*ω* *t*) at the left end. It is a weak nonlinearity if *σ* ≪ 1. The values of the parameters used in this model are as follows: *m* = 1, *m*<sub>r</sub> = 0.5, *β*<sub>0</sub> = 10π, *β*<sub>1</sub> = 15π, and *β*<sub>2</sub> = 10<sup>5</sup>. If *A*<sub>0</sub> = 0.005 and *σ* = 0.16, it is moderate nonlinearity. To increase *σ*, one generally needs to increase the drive because it is difficult to increase *k*<sub>2</sub> in the macroscopic mechanical systems. Just recently, Huang *et al* [47] used chemical bonding interactions to generate a huge *k*<sub>2</sub> on the size scale of about 1 µm.
<figure><img src="image_6.png" /><figcaption>Figure 2. Model of the tetratomic NAM.</figcaption></figure>

Applying the Bloch theorem for periodic structure, the displacement reads $\mathbf{u}(r + a) = \mathbf{u}(r)\exp(i\kappa a)$, where $\kappa$ is the wave vector, and $a$ is the lattice constant. Dispersion solutions are eigen waves in the infinite metamaterial without external input. We can solve the dispersion relations with the perturbation approach, HBM, or homotopy analysis method. These approaches are described in [46].

We have considered a diatomic NAM with a finite length: eight unit cells are involved in the model. Actually, further increasing the number of cells hardly changes the intrinsic characteristics of the NAM but would decrease the accuracy of the bifurcation diagram. Therefore, there are 16 degrees-of-freedom (DoF) and the state space has 32 dimensions (32D, nD symbolizes n-dimensional). In practice, applications would imply the dynamic of waves in a finite system with external excitations. The monochromatic sinusoidal wave $u_0(t)$ drives the system from the left end. The right end is free, and the terminal linear oscillator is analyzed.

There is only one LR bandgap in the diatomic model. In practice, besides the LR bandgap, there are different Bragg bandgaps in metamaterials. The nonlinearity would give rise to different behaviors in these multiple bandgaps. Though the tri-atomic model presents the case of having a LR bandgap and a Bragg bandgap, here we discuss a more complex tetratomic model to determine the influences on multiple bandgaps.

## 2.2. Tetratomic model

The model of the tetratomic NAM is illustrated in figure 2. A ‘molecule’ of this model is composed of three identical linear ‘atoms’ and one nonlinear ‘atom’. Their parameters are labeled in the figure, where $x_n, y_n, z_n$, and $r_n$ denote the displacements of corresponding oscillators in the n-th cell, respectively. The tetratomic NAM has more complicated dynamics.

The generalized motion equation of the n-th cell can be simplified as

$$ 
\begin{cases}
\ddot{x}_n = \beta_0(y_n + z_{n-1} - 2x_n) + \eta[\zeta(\dot{r}_n - \dot{x}_n) + \beta_1(r_n - x_n) + \beta_2(r_n - x_n)^3] \\
\ddot{y}_n = \beta_0(x_n + z_n - 2y_n) \\
\ddot{z}_n = \beta_0(x_{n+1} + y_n - 2z_n) \\
\ddot{r}_n = -\zeta(\dot{r}_n - \dot{x}_n) - \beta_1(r_n - x_n) - \beta_2(r_n - x_n)^3
\end{cases}
\quad (2) $$

The definitions of other parameters in equation (2) are identical to those in equation (1).

The algorithms to analyze the bifurcation and characteristics of chaos in the tetratomic model are also identical to those in the diatomic model. The difference is the chain we considered here has only 4 cells, so the state space is also 32D. The parameters of the model are as follows: $m = 1, m_r = 2, \eta = 2, \beta_0 = 225 = \omega_s^2, \beta_1 = 25 = \omega_{sr}^2, \beta_2 = 5 \times 10^4$, and $\Omega = \omega/\omega_s$. $\omega_s$ is 3 times larger than $\omega_{rs}$, which is similar with the continuous metamaterials that require soft resonators. A larger but still weak damping $\zeta = 0.02$ is adopted. Assuming that the driving amplitude is $A_0 = 0.01$ and the nonlinear strength is $\sigma = 0.6$, it is a strongly nonlinear system in this case.

# 3. Theories

## 3.1. Method in frequency domain

(1) Some algebraic laws.

The analysis of the nonlinear dynamics of high-dimensional systems generally involves high-order operations on matrices and vectors. For convenience, let us define at first some algebraic operations.

We define the ‘element product’ of two nD columns vectors $\mathbf{x} = [x_1 \ x_2 \ \dots \ x_n]^T$ and $\mathbf{y} = [y_1 \ y_2 \ \dots \ y_n]^T$ to be an nD vector $\mathbf{z}$, and the components of which are

$$ \mathbf{z} = \mathbf{x} \circ \mathbf{y} = [\dots x_i y_i \dots]^T, \quad \text{i.e. } z_i = x_i y_i. \quad (3) $$

The element product of two vectors can be abbreviated as $xy$. Therefore, the notation for the m-th power of $\mathbf{x}$ is $\mathbf{x}^m = [\dots x_i^m \dots]^T$. It can be easily established that the element product is commutative $xy = yx$, distributive $(\mathbf{x} + \mathbf{y})\mathbf{z} = \mathbf{x}\mathbf{z} + \mathbf{y}\mathbf{z}$, and associative $xyz = zyx$.
Similarly, we define the element product of two matrices as being $\mathbf{A} \circ \mathbf{B} = [a_{ij}b_{ij}]$, where the dot '$\circ$' cannot be omitted. It is straightforward to show that the element product of two matrices obeys the commutative, distributive, and associative laws as well. However, the product between a matrix and an element product of two vectors does generally not observe the associative law, i.e. $\mathbf{A}(\mathbf{x} \circ \mathbf{y}) \neq (\mathbf{A}\mathbf{x}) \circ \mathbf{y}$, which could become an issue in practice. To overcome this obstacle and ensure this operation becomes associative, we set the following transformation

$$ \mathbf{A}(\mathbf{x} \circ \mathbf{y}) = (\mathbf{A} \circ [[\mathbf{x}]])\mathbf{y}, \quad (4) $$

where the notation [[x]] stands for the square matrix built upon the vector x, namely

$$ [[\mathbf{x}]] = \left[ \frac{\mathbf{x} \quad \mathbf{x} \quad \cdots \quad \mathbf{x}}{n} \right]^T. $$

With this definition, together with the commutativity, it comes $(\mathbf{A} \circ [[\mathbf{x}]])\mathbf{y} = (\mathbf{A} \circ [[\mathbf{y}]])\mathbf{x}$. Moreover, it can be easily shown that the partial derivative of the element product is

$$ \frac{\partial \mathbf{A}(\mathbf{x} \circ \mathbf{y})}{\partial \mathbf{y}} = \mathbf{A} \circ [[\mathbf{x}]], \quad \frac{\partial \mathbf{A}(\mathbf{x}^m)}{\partial \mathbf{x}} = m \mathbf{A} \circ [[\mathbf{x}^{m-1}]]. \quad (5) $$

These algebraic rules allow for a convenient analysis of the high-dimensional nonlinear equations and lead to compact formulas.

## (2) Frequency responses.

It is not an easy task to achieve the analytical solutions of equations (1) and (2) for a high-dimensional system. However, the frequency responses can be approximated with the help of the numerical integral method, the averaging method, or the HBM. The latter can be implemented to find the approximate steady frequency responses. To this end, we define the coordinate transformation $\hat{\mathbf{y}}_n = \mathbf{y}_n - \mathbf{x}_n$ in equation (1), or $\hat{\mathbf{y}}_n = \mathbf{r}_n - \mathbf{x}_n$ in equation (2). With this transformation, the equation of motion for the finite nonlinear metamaterial model reads

$$ \mathbf{M}\ddot{\mathbf{y}} + \mathbf{C}\dot{\mathbf{y}} + \mathbf{K}\mathbf{y} + \mathbf{N}\mathbf{y}^3 = \mathbf{f} \cos \omega t, \quad (6) $$

where $\mathbf{M}, \mathbf{C}, \mathbf{K}$, and $\mathbf{N}$ are mass, damping, stiffness, and nonlinear coefficient matrices, respectively, of the whole transformed system; $\mathbf{f}$ stands for the node force vector applied on every mass.

Let us assume that the steady response has the form $\mathbf{y} = \mathbf{a} \cos \omega t + \mathbf{b} \sin \omega t$, where $\mathbf{a}$ and $\mathbf{b}$ are constant vectors. The first-order HBM leads to a system of algebraic equations:

$$ \begin{cases} [\mathbf{K} - \omega^2 \mathbf{M}] \mathbf{a} + \omega \mathbf{C} \mathbf{b} + 3 \mathbf{N} ((\mathbf{a}^2 + \mathbf{b}^2) \mathbf{a}) / 4 = \mathbf{f} \\ [\mathbf{K} - \omega^2 \mathbf{M}] \mathbf{b} - \omega \mathbf{C} \mathbf{a} + 3 \mathbf{N} ((\mathbf{a}^2 + \mathbf{b}^2) \mathbf{b}) / 4 = \mathbf{0} \end{cases} \quad (7) $$

The solutions of equation (7) would accurately describe the responses $\mathbf{Y} = \sqrt{\mathbf{a}^2 + \mathbf{b}^2}$ to the driving force. However, the HBM does not directly allow investigating the stability of solutions. This is the reason why we adopted the harmonic average approach (HAA) [48, 49] instead, to study the frequency responses. Within this approach, the solution is assumed to have the following form:

$$ \begin{cases} \mathbf{y} = \mathbf{u}(t) \cos \theta + \mathbf{v}(t) \sin \theta \\ \dot{\mathbf{y}} = -\omega \mathbf{u}(t) \sin \theta + \omega \mathbf{v}(t) \cos \theta \end{cases} \quad (8) $$

where $\theta = \omega t$. The derivatives with respect to time $t$ of the formula in (8) are

$$ \begin{cases} \dot{\mathbf{y}} = (\dot{\mathbf{u}} + \omega \mathbf{v}) \cos \theta + (\dot{\mathbf{v}} - \omega \mathbf{u}) \sin \theta \\ \ddot{\mathbf{y}} = (\omega \dot{\mathbf{v}} - \omega^2 \mathbf{u}) \cos \theta - (\omega \dot{\mathbf{u}} + \omega^2 \mathbf{v}) \sin \theta \end{cases} \quad (9) $$

Comparing the expressions of $\dot{\mathbf{y}}$ in equations (8) and (9), we obtain

$$ \omega \mathbf{M}(\dot{\mathbf{u}} \cos \theta + \dot{\mathbf{v}} \sin \theta) = \mathbf{0}. \quad (10) $$

The further substitution of (8) and (9) into (6) gives another form of the equation of motion:

$$ \omega \mathbf{M}(\dot{\mathbf{v}} \cos \theta - \dot{\mathbf{u}} \sin \theta) + [\mathbf{K} - \omega^2 \mathbf{M}](\mathbf{u} \cos \theta + \mathbf{v} \sin \theta) + \omega \mathbf{C}(\mathbf{v} \cos \theta - \mathbf{u} \sin \theta) \\ + \mathbf{N}\{ (\mathbf{u}^3 - 3\mathbf{u}\mathbf{v}^2) \cos^3 \theta + (\mathbf{v}^3 - 3\mathbf{u}^2\mathbf{v}) \sin^3 \theta + 3\mathbf{u}^2\mathbf{v} \sin \theta + 3\mathbf{u}\mathbf{v}^2 \cos \theta \} = \mathbf{f} \cos \theta \quad (11) $$

Then, assuming that $\mathbf{u}$ and $\mathbf{v}$ are constant, by calculating (11) × $\sin \theta$ - (10) × $\cos \theta$ and integrating the result from 0 to $2\pi$, we obtain
$$2\omega\dot{\mathbf{u}} = \mathbf{M}^{-1}[(\mathbf{K} - \omega^2\mathbf{M})\mathbf{v} - \omega\mathbf{C}\mathbf{u} + 3\mathbf{N}((\mathbf{v}^2 + \mathbf{u}^2)\mathbf{v})/4]. \quad (12)$$

Similarly, calculating (11) × cos θ + (10) × sin θ and integrating over the interval [0, 2π] gives

$$2\omega\dot{\mathbf{v}} = -\mathbf{M}^{-1}[(\mathbf{K} - \omega^2\mathbf{M})\mathbf{u} + \omega\mathbf{C}\mathbf{v} + 3\mathbf{N}((\mathbf{v}^2 + \mathbf{u}^2)\mathbf{u})/4 - \mathbf{f}]. \quad (13)$$

The steady solutions correspond to the condition $\dot{\mathbf{u}} = \dot{\mathbf{v}} = \mathbf{0}$, and their expressions are therefore as follows:

$$\left\{ \begin{array}{l} (\mathbf{K} - \omega^2\mathbf{M})\mathbf{v} - \omega\mathbf{C}\mathbf{u} + 3\mathbf{N}((\mathbf{v}^2 + \mathbf{u}^2)\mathbf{v})/4 = 0 \\ (\mathbf{K} - \omega^2\mathbf{M})\mathbf{u} + \omega\mathbf{C}\mathbf{v} + 3\mathbf{N}((\mathbf{v}^2 + \mathbf{u}^2)\mathbf{u})/4 = \mathbf{f}. \end{array} \right. \quad (14)$$

The amplitude of the response is $\mathbf{Y} = \sqrt{\mathbf{u}^2 + \mathbf{v}^2}$.

The derivation of these solutions using HAA turns the steady response problem into an equilibrium problem of differential equations. Actually, solutions in equations (14) and (7) are equal. However, HAA allows for analyzing the stability of the solutions of both systems of differential equations (12) and (13) via the Jacobian matrix. Omitting the terms in $2\omega \neq 0$ in the left sides of (12) and (13) since they will not influence the properties of the solutions, the Jacobian matrix $\mathbf{J}$ is derived by performing the derivations of (12) and (13) with respect to the vectors $\mathbf{u}$ and $\mathbf{v}$, namely,

$$\mathbf{J} = \begin{bmatrix} \mathbf{M}^{-1} & \mathbf{0} \\ \mathbf{0} & -\mathbf{M}^{-1} \end{bmatrix} \begin{bmatrix} -\omega\mathbf{C} + 3\mathbf{N} \circ [[\mathbf{u}\mathbf{v}}]]/2 & \mathbf{K} - \omega^2\mathbf{M} + 3\mathbf{N} \circ [[3\mathbf{v}^2 + \mathbf{u}^2]]/4 \\ \mathbf{K} - \omega^2\mathbf{M} + 3\mathbf{N} \circ [[3\mathbf{u}^2 + \mathbf{v}^2]]/4 & \omega\mathbf{C} + 3\mathbf{N} \circ [[\mathbf{u}\mathbf{v}}]]/2 \end{bmatrix}. \quad (15)$$

Therefore, with the solutions coming from equation (15), their stability can be determined: if the real part of an eigenvalue of $\mathbf{J}$ is positive, the corresponding steady solution is unstable; if a real eigenvalue goes from negative to positive, saddle-node (SN) bifurcation occurs, and if a pair of conjugate complex eigenvalues crosses the imaginary axis, Hopf bifurcation appears [50]. The succinct expressions derived above are universal for NAMs and differential equations with cubic nonlinearity. Let us note that if $\mathbf{C} = \mathbf{0}$, this formula reveals that $\mathbf{v}$ must read $\mathbf{v} = \mathbf{0}$.

Continuation algorithm is further adopted to solve equation (14). Within this frame, the Newton method (MATLAB fsolve) is used to find the solutions. The concept of continuation consists in deriving the solution ($\lambda_n, x_n$) at step $n$ from the solution ($\lambda_{n-1}, x_{n-1}$) at step $n-1$; here, $\lambda$ is the bifurcation parameter, and $x$ is a vector of solution. In this paper, we propose a 'perturbation continuation concept' in which $\lambda$ is specified and $x$ is unknown. If one has found the solution ($\lambda_{n-1}, x_{n-1}$) with an initial guess $X_{n-1}$, the next solution ($\lambda_n, x_n$) is derived from a perturbed initial guess $X_n = (1 + \varepsilon)x_{n-1}$, where $\varepsilon$ is the perturbation parameter. This concept has been very efficient for artificially finding and starting new branches at a specified point.

## 3.2. Bifurcation analysis method in time domain

From here on, the periodic solutions will be mostly derived using time domain methods, which are parts of software such as AUTO [51]. Within this kind of approach, integrations of equations of motion are made in the state space. Furthermore, the method deriving LEs [52, 53] is also based on numerical integrations.

The incident elastic wave transforms the model in a non-autonomous system. However, when calculating the spectra of LEs and bifurcation diagrams, one must transform the non-autonomous system to an autonomous one. To calculate the former, the input wave is expressed as

$$u_0 = A_0 \sin \theta, \quad \dot{\theta} = \omega, \quad \theta(0) = 0 \quad (16)$$

i.e. it becomes a boundary value problem. Therefore, it becomes a 3D system in this situation.

When calculating the bifurcation diagrams, a differential system is added.

$$\left\{ \begin{array}{l} \dot{\phi} = \phi + \omega\phi - \phi(\phi^2 + \phi^2) \\ \dot{\varphi} = -\omega\phi + \varphi - \varphi(\phi^2 + \varphi^2). \end{array} \right. \quad (17)$$

This system has a unique asymptotically stable solution, $\phi = \sin \omega t$, $\varphi = \cos \omega t$. The boundary value at a specified time $t_0$ is defined as $\phi_0 = \sin \omega t_0$, $\varphi_0 = \cos \omega t_0$. Therefore, $u_0 = A_0\phi$ and the whole system becomes 3D.

We used the program AUTO [51] to analyze the bifurcations of periodic solutions. AUTO is based on a Newton iterative scheme to find the solutions. The continuation method is then used to find the branches of solutions. The bifurcations and stabilities of periodic solutions are identified by Floquet multipliers $\mu_i$. Let us recall here that, for a continuous differential system, there are three types of codimension-1 bifurcations: SN, period-doubling (PD), and torus (TR) bifurcations [54]. The bifurcation associated with the appearance of $\mu_i = 1$ ($\mu_i = -1$) is called an SN (PD) bifurcation; the bifurcation corresponding to the presence of conjugate multipliers $\mu_{1,2} = \exp(\pm i\theta_0)$, $0 < \theta_0 < \pi$, is called a TR bifurcation. An invariant TR is a quasiperiodic solution. Both PD bifurcations and breakdowns of invariant tori may induce chaos [50]. At present, most studies
on the bifurcation are focusing on the low-dimensional systems. In our 34D model, there are multiple pairs of conjugate multipliers.

Moreover, the QR decomposition algorithm proposed by Wolf *et al* [52] is adopted to compute the spectra of LEs of the 33D system. There are 33 LEs: $\lambda_1, \lambda_2, \dots, \lambda_{33}$ (in descending order). $\lambda_1$ is named the largest LE (LLE). All non-autonomous systems have at least one zero LE that corresponds to the *t*-component. Therefore, if $\lambda_1 = 0$, we discard it and make $\lambda_n^* = \lambda_{n+1}$ be the new spectrum. If $\lambda_1 > 0$, the motion is Lyapunov chaotic, and if there are more than two positive $\lambda_i$, the motion is described as Lyapunov hyperchaotic.

Let the quantity $s_\lambda = \sum_i \lambda_i^+$ be the sum of all positive LEs. The strength of chaos is determined by the average ratio of the exponential divergence of neighbor orbits in a chaotic attractor. Therefore, $s_\lambda$ may be used to describe the strength of chaos, and the larger $s_\lambda$, the stronger the chaos. Actually, noting that a larger $\lambda_1$ corresponds to a larger $s_\lambda$, it can be more conveniently used to characterize the chaos. Accordingly, if $\lambda_1^+ \to 0$, it is weak chaos and would behave as a quasiperiodic motion.

Finally, the Lyapunov dimension $d_{LD}$ is defined as $d_{LD} = j - \sum_{i=1}^{j} \lambda_i / \lambda_{j+1}$, where *j* is the maximum value of *j* that makes $\sum_{i=1}^{j} \lambda_i > 0$, i.e. $\sum_{i=1}^{j} \lambda_i > 0$, $\sum_{i=1}^{j+1} \lambda_i < 0$. For a high-dimensional system (>2D), $d_{LD}$ can accurately describe the fractal nature of the chaotic attractor. For a chaotic attractor, $0 < d_{LD} < N_d$, where $N_d$ is the dimension of the system. For periodic and quasiperiodic motions, $d_{LD} = 0$. A larger $d_{LD}$ manifests a more complex chaos.

# 4. Properties of wave propagation in NAM models

## 4.1. Diatomic model

Bandgap properties of the diatomic NAM model have been laid out in [45]. It is shown that for this system, the generalized frequency range of the LR bandgap is $\Omega \in [1.088, 1.5]$, and the cutoff frequency of the passband is $\Omega_c = 2.251$. The transmission coefficient is defined as $T_A = A_{max}/A_0$, where $A_{max}$ denotes the maximum amplitude of the output wave, i.e. the maximum amplitude of the last linear oscillator in the chain. $A_{max}$ represents the extreme condition in the long-term chaotic responses. $T_A$ can be derived via three different approaches as shown in figure 3. The hardened nonlinear stiffness leads the lower bound frequency of the LR bandgap to shift upward. On the other hand, when the driving amplitude is set to $A_0 = 0.005$, the bandgap calculated by the homotopy analysis method is $\Omega \in [1.2, 1.758]$ and $\Omega_c = 2.277$, as shown in figure 3(a) [45].

Actually, the passbands of the finite LAM are composed of discrete resonances that become denser when increasing the periodic numbers. As shown in figures 3(a) and (b), the appearance of nonlinearities significantly reduces and even suppresses resonances. However, there are some differences between the numerical integral results and HAA results, especially in the second passband (i.e. the optical band). There are two reasons for this: (1) there is an intrinsic randomness in the chaotic responses, so $A_{max}$ manifests the extreme case, but this does not mean that the system gets a steady amplitude $A_{max}$; and (2) the solution of HAA is an approximately averaged result at a specified frequency, and a more accurate result should require high-order harmonic components. Despite all this, HAA well depicts the nonlinear responses as well as the role played by nonlinearities.

Different from the linear resonances with infinite amplitudes, the nonlinear modes have finite phase volumes so that the amplitudes of the 16 resonances remain bound [45]. Increasing the amplitude reduces the wave transmissibility; meanwhile, the nonlinear upshot of suppressing resonances extends to the first passband (i.e. the acoustic band). This result is further demonstrated by numerical integrations.

Effects on the responses of a weak damping are depicted in figure 3(b). A weak damping only a little attenuates linear resonances. Further comparison shows that it also has only little impact on the amplitudes in NAMs. As shown in figures 3(a) and (b), the behaviors are complex near the resonances of the NAM, especially when the nonlinearity is strong ($A_0 = 0.005$), leading to the occurrence of jumps. In figure 3(b), only one solution (obtained from zero initial value) is shown for a specified frequency. To understand the phenomena near the resonances, let us take the responses in the interval $\Omega \in [1.65, 2.1]$ as an example. The perturbation continuation is combined with the optimized Newton method to find multiple solutions in this frequency range, as illustrated in figures 3(c) and (d).

When the nonlinearity is weak ($A_0 = 0.001$), curves of frequency responses of the hardened NAM bends toward the high-frequency region, as what happens in a single DoF nonlinear system [55]. However, both branch-1 and branch-2 become unstable near resonances then. The solution is stable along branch-3, and a tangential bifurcation (the SN bifurcation) connects branch-2 and branch-3. When the driving amplitude further increases to $A_0 = 0.005$, multiple complex branches bending to the right are found, and all these solutions are unstable. The low-amplitude parts of these branches are discrete, but in the high-amplitude region, as a result of the strong nonlinearity and dense resonances, these branches interweave or converge into a peak. Moreover, the transmission peaks for $A_0 = 0.005$ are much lower than those obtained with $A_0 = 0.001$. These results illustrate the fact that, near a resonance, a high-dimensional NAM has a behavior different from that of a
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Frequency responses of diatomic LAMs and NAMs. (a) Maximum responses solved directly with the numerical integral method; here $A_0 = 0.005$, and $\zeta = 0$ [45]. (b) Solutions found with the Newton iterative method based on equation (14) using 0 as the initial value. In (a), (b), ‘L’ and ‘N’ denote the ‘linear’ and ‘nonlinear’ models, respectively. In (c), (d) the detailed and enlarged local frequency responses are the sets of solutions near a resonance found by the perturbation continuation algorithm; (c) $A_0 = 0.001$ and $\zeta = 0.01$, (d) $A_0 = 0.005$ and $\zeta = 0.01$. The stabilities are determined using the Jacobian matrix in equation (15). The three numbers in (c) denote the three branches of a resonance, respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Influence of damping on the LLEs ($\Omega = 0.678$).</figcaption></figure>

single DoF or even of a 2DoF system [55], especially for strongly nonlinear ones. For both weak and strongly nonlinear AMs, although it is difficult to determine whether the amplitude may generate a ‘jump’ at the catastrophe frequency, the broadband unstable solutions provide an opportunity for the chaos to appear.

In contrast to the steady amplitude, the damping has a significant impact on the characteristics of the chaotic attractors. As shown in figure 4, in the chaotic region the LLE remarkably decreases upon damping as weak as $\zeta = 0.01$. Therefore, a strong damping effect could allow a chaotic solution to turn into a periodic solution. Hereafter, we consider a weak damping $\zeta = 0.01$ in the diatomic system.

We have chosen three illustrative frequencies $\Omega = 0.678, 1.124$, and $2.2$, whose locations are labeled in figure 3(a). Bifurcation diagrams and chaotic characteristics are illustrated in figure 5 for each of these frequencies.

For the conservative 34D system, there are dense TR points along the branches of periodic solutions because of the multiple pairs of conjugate multipliers, which indicates that there are extensive quasiperiodic solutions in the system. The damping does not alter the branches and stabilities of periodic solutions in the bifurcation diagram because the linear damping $\zeta\dot{x}$ does not change the properties of the equilibriums in the time-domain algorithm (the Jacobian matrix from state space). However, phase volumes of the damped systems shrink in the long-term responses, which results in the Floquet multipliers on the unit circle entering the circle. Therefore, the damping $\zeta = 0.01$ reduces the number of TR points along the branches, and the PD points found in the conservative system disappear. In some frequency range, the TR point disappears too, as illustrated by figure 5(c). Results indicate that the damping reduces the number of quasi-periodic solutions in the system.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Bifurcation diagrams and chaos of diatomic model. (a) $\Omega = 0.678$, (b) $\Omega = 1.124$, and (c) $\Omega = 2.2$. $\zeta = 0.01$. In each figure, the four diagrams are corresponding to the bifurcation diagram, the mean amplitude of the last linear oscillator, the LEs ($\lambda_1, \lambda_3, \lambda_5$), and the Lyapunov dimension $d_{LD}$, respectively. In bifurcation diagrams, the solid black (dashed green) lines represent the stable (unstable) periodic solutions, and the red dots (blue asterisks) denote the SN (TR) bifurcation point. The iconographies refer to the detail with enlarged scale. The value $L_2$-norm is defined as $\sqrt{\int_0^1 U_8(x)^2 dx}$, where $U_8(x)$ is the generalized displacement of $u_8$, and $u_8$ is the displacement of the eighth linear oscillator away from the excitation.</figcaption></figure>

The unstable periodic solution is an essential condition for the chaos to appear. The bifurcation diagram, mean amplitude that directly results from numerical integration, LEs, and $d_{LD}$ in figure 5 are well consistent with each other.

In the case of $\Omega = 0.678$, the first unstable solution appears at $A_0 = 1.867 \times 10^{-3}$ marked by the vertical black line in figure 5(a). Within the subsequent domain, the stable and unstable periodic solutions appear alternately. When $A_0 > 6.686 \times 10^{-3}$, the periodic solution enters a constantly unstable region, in which multiple branches are generated with many TR points. The terminal points of the new branches are SN points. In the interval $1 \times 10^{-3} < A_0 < 1.867 \times 10^{-3}$, $\lambda_1^+ \to 0$ and therefore the response is close to the quasiperiodic responses. Both behaviors of LEs and $d_{LD}$ attest that the periodic motion turns into chaos when the driving amplitude continues to increase. Actually, both $\lambda_1^+$ and $d_{LD}$ tend to increase as $A_0$ increases, meaning that the weak chaos becomes a strong chaos. However, LEs and $d_{LD}$ have local fluctuations. For the amplitude at which chaos appears, only one LE is positive, but there are quickly more than two positive LEs, i.e. hyperchaos occurs. Therefore, for the 33D NAM model, there is a general rule that the system is hyperchaotic when chaos occurs. Moreover, the variation laws of $d_{LD}$ are identical with $\lambda_1$ in the chaotic domain, and when strong chaos occurs $d_{LD} \to N_d = 33$.

Combing the bifurcation diagram and LE spectrum, it is known that chaos at TR points would be induced by the breakdowns of invariant tori. A rigorous demonstration of this mechanism is based on the Kolmogorov-Arnold-Moser (KAM) theorem [50]. Such a demonstration falls out of the scope of this article, which rather focuses on the influence of chaos on the responses. At other unstable points, the chaos arises from PD bifurcations. The waterfall plot of power spectra illustrates this process (figure 6). When $A_0 \to 0$, we are dealing with a 1:1 resonance. In this situation, the energy gets localized at the driving frequency $\Omega_e$. Further increasing $A_0$ leads to odd harmonic waves $3\Omega_e$ and then $5\Omega_e$ to appear; thus, the propagating waves become quasiperiodic. At the critical amplitude $A_{0c}$, the propagating wave cascades into chaos. $A_{0c}$ has the same behavior as LEs. In the route toward chaos, energy localization switches to energy dispersion, and therefore the energy is pumped [56] and spread within a broad high-frequency passband and even within the stop bands.

If the NAM is excited by the wave at the bottom of the LR bandgap, the periodic solution exhibits multiple folded branches that induce multiple jumps, as shown by the case $\Omega = 1.124$ in figure 5(b). There are still some discrete TR points on the unstable branches. The SN bifurcation points act as the terminals of branches. The first SN bifurcation occurs at $A_0 = 3.807 \times 10^{-3}$. When $A_0 < 3.807 \times 10^{-3}$, the periodic solutions remain stable and their amplitude tends to zero; both LEs and $d_{LD}$ reveal that the long-term motion of the system is periodic. The amplitudes issued from numerical integration also approaches zero. This result means that the LR bandgap remains a complete stop band and that the bandgap is linearly stable in this interval. After this SN point, the periodic solutions become constantly unstable and a jump occurs. Both the bifurcation diagram and the mean amplitude show that the motion jumps to a high-energy chaotic branch at this SN point. Then, as illustrated both by LEs and by $d_{LD}$, the propagation through the NAM is along a chaotic way. Meanwhile, the hardened nonlinear stiffness shifts up the lower boundary of the LR bandgap.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** 3D waterfall plot of power spectrum density (PSD). The excitation frequency is $\Omega_c = 0.678$. The different colors correspond to the value $\log_{10}(\text{PSD})$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Sketch plot of quantum states in NAMs. The solid circles represent the excited quantum states and the dashed curves denote the fluctuations in the neighborhood of these states. The shade circle and the black point in the center represent the bound quantum state.</figcaption></figure>

Along this chaotic branch, the mean amplitude has almost constant value with LEs and $d_{LD}$ fluctuating in a neighborhood of the corresponding constants. Another SN bifurcation point appears at $A_0 = 7.676 \times 10^{-3}$, where the motion jumps to a higher chaotic orbit. However, this point does not fit well the long-term motion, as attested by the jumps of the mean amplitude to a higher-energy chaotic orbit at $A_0 = 0.013$. Like the mean amplitude, LEs and $d_{LD}$ remain almost constant. These results show that the system has multiple approximately constant states in this frequency range, and that only discrete values of these states can be achieved: the system behaves as a ‘quantum object’, as sketched in figure 7. Therefore, when manipulating waves in bandgaps via amplitudes, the system can switch suddenly between the stop state (i.e. the bound state) and the propagation states (i.e. the excited quantum states). These characteristics may be helpful to realize acoustic devices with small dimensions, such as acoustic diodes or switches. As shown in figure 5(b), it should be noted that, at some value of $A_0$, the amplitude may jump from the highest obit (the third one) to a lower one. However, these jumps are unstable and a small perturbation may stimulate a jump back to the highest orbit.

The case of $\Omega = 2.2$ in figure 5(c) is similar to the case of $\Omega = 1.124$ except that the motion turns from periodic into chaotic for only one jump at the first unstable periodic solution. The trend shown in the bifurcation diagram agrees well with the mean amplitude variations. After the jump, the amplitude first slowly increases and then remains constant as do both LEs and $d_{LD}$. This suggests a mechanism to suppress elastic waves in strongly NAMs: strengthening the nonlinearity decreases the wave transmissibility.

More bifurcation diagrams at other frequencies are depicted in figure 8. According to the statements above, it can be anticipated that the chaos would occur at the first unstable periodic solutions. In the examined amplitude range, neither bifurcations nor unstable solutions occur on the low-frequency periodic orbits ($\Omega \le 0.3$), and the amplitude linearly increases with driving amplitude, i.e. the transmissibility is constant and the linear regime is preserved. At higher frequency the motion turns to weak chaos that behaves as a quasiperiodic wave, which explains why the responses of NAMs in the low-frequency parts of the acoustic branch are similar to those of LAMs. The transmission shows this feature (see above and figure 3). In the acoustic branch, strong chaos can be observed at frequencies near the LR bandgap.

Bifurcation diagrams in the optical branch have similar characteristics. The amplitudes increase along an orbit where both periodic and chaotic behaviors may occur. However, when the motions jump to a higher orbit at a SN point, the amplitudes of the periodic solutions do not increase with the input wave. Instead, they remain constant or even decrease.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Bifurcation diagrams of the diatomic model under different frequencies. The frequencies are labeled in the corresponding plots. The locations of these frequencies are marked in figure 3(a). Other specifications are identical to those in figure 5.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Sketch of the band structure of the diatomic NAM model. The green lines represent the transmission of the NAM.</figcaption></figure>

In the bandgap, the periodic orbits at $\Omega = 1.25$ feature folding and multiple jumps, as is the case for $\Omega = 1.124$. However, a larger critical amplitude is required to force the wave into chaos. This is because the influences of the nonlinearity on the bandgap spread from low to high frequencies. This is further demonstrated by the bifurcation diagram at $\Omega = 1.4$: although there are unstable solutions, the amplitudes are so small that only subharmonic and superharmonic waves can propagate, and this frequency interval keeps its ability to reflect incident waves. In the high-frequency stop band, as shown by the bifurcation diagram at $\Omega = 2.4$, much stronger nonlinearity is needed to generate unstable solutions.

Based on these results, the band structure of a strongly diatomic NAM is sketched in figure 9.

## 4.2. Tetratomic model

The dispersion relations of this model were elaborated in [46]. In the linear regime, there are three bandgaps for this metamaterial: an LR bandgap and two Bragg (BG) bandgaps. The frequency ranges of the linear bandgaps are as follows: LR, [0.3065, 0.4195]; BG1, [1, 1.078]; and BG2, [1.732, 1.777]. It is shown in [46] that, in the nonlinear regime, BG1 and BG2 are not altered if the amplitude $A_0 = 0.01$, but the LR bandgap is shifted to [0.42, 0.5807]. The cutoff frequency of the passband is $\Omega = 2$.

We have considered a weak damping $\zeta = 0.02$ in the tetratomic NAM comprised of four cells. $T_A$ computed with different algorithms, and the transfer function $H(\omega)$ upon sine-sweep excitations [46], are shown in figure 10. For the LAM model, thanks to its strong sensibility to localization properties, only the LR bandgap is clearly evidenced by the four-cell chain. In contrast, more cells are needed for the BG bandgaps to open up. With regard to the nonlinear model, this allows for a fair position of the BG bandgaps. As is the case with the diatomic model, the tetratomic NAM also proves that broadband elastic waves are suppressed. The results from HAA establish that a strong nonlinearity causes the wave transmissibility to decrease. However, the elastic waves can
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Frequency responses of tetatomic LAMs and NAMs. $\zeta = 0.02$. (a) $T_A$ and $H(\omega)$ solved directly using the numerical integral method, where $A_0 = 0.01$, ‘L’ and ‘N’ denote the ‘linear’ and ‘nonlinear’ models, respectively; (b) the solutions found with the Newton iterative method using 0 as the initial value; (c) multiple solutions found by the perturbation continuation algorithm, where blue and brown lines denote the case $A_0 = 0.002$, red and black lines represent the case $A_0 = 0.01$, and the stabilities are distinguished with solid and dashed lines.</figcaption></figure>

still propagate in the nonlinear LR bandgap. The continuation method emphasizes the nonlinear modes that compose the unstable peaks and multiple branches. Overall, these behaviors are similar to that of the diatomic NAM. As shown in figure 10(c), the first resonance in the case $A_0 = 0.01$ germinates multiple bifurcations that start unstable branches, which is different from other unstable nonlinear resonant peaks with only one branch. Moreover, stable and unstable solutions appear alternatively along all the curves in this case. However, in the amplitude range $A_0 \le 0.01$, the nonlinear modes composing of complex unstable branches are not found.

With the discrete models, in conjunction with the results of the diatomic NAM, we can conclude that the hardened NAM actually influences all passbands higher than the nonlinear LR bandgap, including a small region below and near this bandgap.

We have further studied bifurcation and chaos to better understand the wave propagation properties in the system. To this end, we have investigated the wave propagation in the passbands. We have considered three frequencies $\Omega = 0.2463, 0.672$, and $1.42$, located in the first, second, and third passbands, respectively, (see figure 10). The properties of the fourth passband are similar to those of the third passband. The bifurcation diagrams, mean amplitudes, LEs, and LDs are shown in figures 11(a)–(c). Comparing these results, it can be seen that increasing the driving amplitude stimulates the motion into chaos. Moreover, there are differences between the chaotic critical points in bifurcation diagrams and LEs, which are caused by (1) the damping that greatly influences the long-term motion, and (2) the different principles implemented by the two algorithms that lead to different critical points.

In the case of $\Omega = 0.2463$, the bifurcation diagram is consistent with the mean amplitude: between $0.006 < A_0 < 0.009$ the general chaos (only one positive LE) arises, but the motion turns into hyperchaos when $A_0$ increases further. However, $\lambda_i^+ < 0.03$ and $d_{LD} \ll 33$, and thus the wave undergoes a weak hyperchaos whose behavior is similar to the quasiperiodic one. Therefore, it can be deduced that the low-frequency waves in the first passband are quasiperiodic or weakly chaotic.

In fact, the propagation of the waves in the subsequent three passbands are similar except for some local differences. Along with the increasing driving amplitude, the variation law of the wave state in whole is ‘periodic $\rightarrow$ weakly chaotic $\rightarrow$ strongly hyperchaotic’. However, in these three passbands, the LEs $\lambda_1$ for the strong chaos are 10 times larger than those calculated for the weak chaos in the first passband, and $d_{LD} \rightarrow 33$. The mechanisms that suppress resonances in the passbands are identical to those observed with the diatomic NAM. Chaos is induced by PD bifurcation, and the wave amplitudes along the chaotic branches slowly increase, remain constant, or even decrease with $A_0$. Moreover, as illustrated in figures 11(b) and (c), periodic windows in the chaotic region are observed, which means that LEs and LDs do not monotonously vary.

Close to $\Omega = 1.42$, $T_A$ is equal to the transmission for the corresponding LAM. In fact, there are other specific frequency domains. Comparison with the results obtained with the ten-cell model [46] indicates that these domains relate to the density of the linear modes that is determined by the length of the chain. In the four-cell model, this domain is a non-resonant region. The periodic solutions for $\Omega = 1.42$ have two branches: a low-energy one and a high-energy one. Along the low-energy branch, there is an unstable domain in between two stable domains. If the motion had been periodic for $A_0 \le 0.01$, the mean displacement of the last oscillator
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Bifurcation diagrams and chaos of tetratomic model. (a) $\Omega = 0.2463$, (b) $\Omega = 0.672$, (c) $\Omega = 1.42$, and (d) $\Omega = 1$. Other specifications of these figures are the same as those in figure 5. $z_4$ denotes the displacement of the last linear oscillator.</figcaption></figure>

should have been linear against $A_0$. Instead, both the mean amplitude and LEs indicate that the system enters a weak chaos regime in the range $0.006 < A_0 < 0.009$. However, this weakly chaotic motion is different from that described above. Here, the oscillators jump to a high-energy unstable orbit. In the subsequent range $0.009 < A_0 \le 0.011$, the system jumps back to the low-energy weakly chaotic state, featuring quasiperiodic behavior. Therefore, when the amplitude is $A_0 = 0.01$, the linear and nonlinear transmissibilities near $\Omega = 1.42$ are equal. Further increasing the driving amplitude makes the system to jump to the strongly hyperchaotic state, featuring constant amplitude and high energy, except in a narrow window around $A_0 = 0.014$ where the motion becomes periodic. However, in this small interval, amplitudes still are in the high-energy orbit.

The jumps between low-energy and high-energy orbits explains why the waves are amplified under non-resonant conditions. This amplification scheme can be applied to design a broadband wave amplifier. For the finite LAM, the waves can be amplified at the discrete resonant frequencies in the passbands. In contrast, a four-cell NAM can allow for a homogeneous amplification in a broad passband.

Another question we have addressed is why the nonlinear LR bandgap is less efficient than other bandgaps to suppress the elastic waves in NAM. To answer this question, we have set frequencies to $\Omega = 0.35$ and $\Omega = 1$. As shown in figure 12, as for the diatomic model, there are multiple bifurcations and branches. However, stable periodic solutions disappear for an amplitude as low as $A_0 = 0.003$ because of the strong nonlinearity; therefore, this bandgap closes at a certain driving amplitude and is replaced by a high-energy chaotic orbit, as is the case in the diatomic model. In contrast, for the BG bandgaps away from this nonlinear bandgap, much higher amplitudes are needed to obtain the unstable periodic solution. In the example illustrated in figure 11(d), when $\Omega = 1$ in the first BG bandgap, a periodic motion is observed and the chaos is weak ($\lambda_1^+ < 0.04$) if the amplitude $A_0 < 0.01$, as for the quasi-linear state in the interval $A_0 = 0.01 \sim 0.015$. It is only when $A_0 > 0.015$ that the
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Bifurcation diagrams at frequency Ω = 0.35.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Sketch of the band structure of the tetratomic NAM model. The green lines represent the transmission of the strongly nonlinear AM model. The dashed blue lines correspond to the much stronger nonlinearity.</figcaption></figure>

motion is attracted to the high-energy hyperchaotic state. Combing the mechanisms in the diatomic model, one concludes that the nonlinear LR bandgap has weaker ability to reflect incident waves when it jumps to high-energy orbits.

Based on the results from this paper and [46], and combing them with some predications, we sketch the band structure of the tetratomic NAM in figure 13. For the LAM here, the low-frequency LR bandgap is narrow. The strong nonlinearity shifts this nonlinear LR upward; meanwhile, the responses in this gap jumps to the excited state. Below the nonlinear gap, the band structures are similar to the diatomic ones. However, the second and the third linear resonant passbands are replaced by strong hyperchaotic passbands with lower transmissions. When the nonlinearity is weak, BG1, BG2, and the fourth resonant passband would remain in the linear or quasi-linear state; however, when the nonlinearity becomes much stronger, BG1 and BG2 would also obtain multi-state behaviors, while the fourth passband becomes a chaotic passband and its upper boundary spreads to the high-frequency stop band. However, a detailed statement about the band structures should be based on the structure, frequency, and amplitude.

**5. Conclusions**

The propagation of waves in NAMs is fundamentally different from that in conventional LAMs. However, these features are still not fully understood. In this work we investigate the elastic wave propagations in the 1D diatomic and tetratomic NAM models. We further demonstrate that nonlinear effects can greatly suppress elastic waves in a broad band. This band consists of bandgaps and chaotic passbands. The nonlinear wave behaviors, band structures, bifurcations, and chaos are studied to demonstrate the novel mechanisms that can manipulate wave propagations.

HAA combined with the continuation are employed to calculate the frequency responses and nonlinear modes. Bifurcations and chaos are analyzed using both the continuation algorithm and the spectra of LEs (and *d*<sub>LD</sub>). Our results show that the nonlinear resonances feature multiple branches with unstable peaks. We have further demonstrated that the PD process leads the motion from being periodic to chaotic, and the damping has a significant influence on the characteristics of chaotic attractors. Moreover, due to the dispersion process in the chaotic regime, the localized energy spreads in the broadband high-frequency passbands and even in the stop
bands. The trend for the positive LEs and $d_{LD}$ is to increase with amplitude so that the weak chaos (or hyperchaos) turns into strong hyperchaos.

Band structures of both diatomic and tetratomic NAMs are studied. The bifurcations and unstable solutions do not occur on the low-frequency periodic orbits, and thus the responses in these regions are similar to those of LAMs. Actually, the hardened NAM influences all passbands that are higher than the nonlinear LR bandgap, including a small region below this bandgap.

In the nonlinear LR bandgap, we observed 'quantum' behaviors because of jumping bifurcations between low-energy and high-energy orbits, whose propagation states (excited states) and stop states (bound states) have discrete characteristics and switch suddenly. This behavior also explained why the nonlinear LR has a weaker ability to reflect incident waves. Moreover, jumps in passbands would amplify non-resonant waves.

This work provides a theoretical base for further understanding and application of NAMs.

## Acknowledgments

This research was funded by the National Nature Science Foundation of China (Project Nos. 51405502 and 51275519).

## References

[1] Cummer S A, Christensen J and Alù A A 2016 *Nat. Rev. Mater.* **1** 16001

[2] Wang G, Wen X, Wen J, Shao L and Liu Y 2004 *Phys. Rev. Lett.* **93** 154302

[3] Liu Z, Zhang X, Mao Y, Zhu Y Y, Yang Z, Chan C T and Sheng P 2000 *Science* **289** 1734

[4] Xiao Y, Wen J and Wen X 2012 *New J. Phys.* **14** 402380

[5] Molerón M, Serra-Garcia M and Daraio C 2016 *New J. Phys.* **18** 033003

[6] Ma G and Sheng P 2016 *Sci. Adv.* **2** e1501595

[7] Rupin M, Lemoult F, Lerosey G and Roux P 2014 *Phys. Rev. Lett.* **112** 234301

[8] Moiseyenko R P, Pennec Y, Marchal R, Bonello B and Djafari-Rouhani B 2014 *Phys. Rev.* B **90** 134307

[9] Huang H H and Sun C T 2009 *New J. Phys.* **11** 013003

[10] Ma G, Min Y, Xiao S, Yang Z and Sheng P 2014 *Nat. Mater.* **13** 873

[11] Zhang H, Xiao Y, Wen J, Yu D and Wen X 2016 *Appl. Phys. Lett.* **108** 141902

[12] Zhao J, Bonello B and Boyko O 2016 *Phys. Rev.* B **93** 174306

[13] Kushwaha M S, Halevi P, Dobrzynski L and Djafari-Rouhani B 1993 *Phys. Rev. Lett.* **71** 2022–5

[14] Lapine M, Shadrivov I V and Kivshar Y S 2014 *Rev. Mod. Phys.* **86** 1093

[15] Leitenstorfer A, Nelson K A, Reimann K and Tanaka K 2014 *New J. Phys.* **16** 045016

[16] Grady N K *et al* 2013 *New J. Phys.* **15** 105016

[17] Liang L, Yuan B and Cheng J C 2009 *Phys. Rev. Lett.* **103** 104301

[18] Liang B, Guo X S, Tu J, Zhang D and Cheng J C 2010 *Nat. Mater.* **9** 989–92

[19] Vakakis A F 1992 *Acta Mech.* **95** 197

[20] Fermi E, Pasta J and Ulam S 1955 *Collected Papers II* (Chicago, IL: University of Chicago Press) pp 977–88

[21] Gallavotti G (ed) 2008 The fermi–pasta–ulam problem: a status report *Lecture Notes in Physics* 728 (New York: Springer)

[22] Romeo F and Rega G 2006 *Chaos, Solitons Fractals* **27** 606–17

[23] Romeo F and Rega G 2015 *Meccanica* **50** 721–30

[24] Nesterenko V F 2001 *Dynamics of Heterogeneous Materials* (New York: Springer)

[25] Daraio C, Nesterenko V F, Herbold E B and Jin S 2006 *Phys. Rev.* E **73** 026610

[26] Fang X, Zhang C H, Chen X, Wang Y S and Tan Y Y 2015 *Acta Mech.* **226** 1657

[27] Nadkarni N, Arrieta A F, Chong C, Kochmann D M and Daraio C 2016 *Phys. Rev. Lett.* **116** 244501

[28] Boechler N, Theocharis G, Job S, Kevrekidis P G, Porter M A and Daraio C 2010 *Phys. Rev. Lett.* **104** 244302

[29] Manktelow K L, Leamy M J and Ruzzene M 2011 *Nonlinear Dyn.* **63** 193

[30] Scalora M, Bloemer M J, Manka A S, Dowling J P, Bowden C M, Viswanathan R and Haus J W 1997 *Phys. Rev.* A **56** 3166

[31] Meurer T, Qu J and Jacobs L 2002 *Int. J. Solids Struct.* **39** 5585

[32] Kim E, Li F, Chong C, Theocharis G, Yang J and Kevrekidis P G 2015 *Phys. Rev. Lett.* **114** 118002

[33] Wang S Y and Nesterenko V F 2015 *Phys. Rev.* E **91** 062211

[34] Pichard H, Duclos A, Groby J-P, Tournat V, Zheng L and Gusev V E 2016 *Phys. Rev.* E **93** 023008

[35] Lydon J, Theocharis G and Daraio C 2015 *Phys. Rev.* E **91** 023208

[36] Boechler N, Theocharis G and Daraio C 2011 *Nat. Mater.* **10** 665

[37] Donahue C M, Anzel P W J, Bonanomi L, Keller T A and Daraio C 2014 *Appl. Phys. Lett.* **104** 014103

[38] Xu Y and Nesterenko V F 2014 *Phil. Trans. R. Soc.* A **372** 20130186

[39] Herbold E B, Kim J, Nesterenko V F, Wang S Y and Daraio C 2009 *Acta Mech.* **205** 85–103

[40] Narisetti R K, Ruzzene M and Leamy M J 2011 *J. Vib. Acoust.* **133** 061020

[41] Narisetti R K, Ruzzene M and Leamy M J 2012 *Wave Motion* **49** 394

[42] Manktelow K, Leamy M J and Ruzzene M 2013 *J. Mech. Phys. Solids* **61** 2433

[43] Yousefzadeh B and Phani A S 2015 *J. Sound Vib.* **354** 180

[44] Bernard B P, Mazzoleni M J, Garraud N, Arnold D P and Mann B P 2014 *J. Appl. Phys.* **116** 084904

[45] Fang X, Wen J, Yin J, Yu D and Xiao Y 2016 *Phys. Rev.* E **94** 052206

[46] Fang X, Wen J, Yin J and Yu D 2016 *AIP Adv.* **6** 121706

[47] Huang P *et al* 2016 *Nat. Commun.* **7** 11517

[48] Cheung Y K, Chen S H and Lau S L 1990 *J. Sound Vib.* **140** 273–86

[49] Xu D L, Zhang H C, Lu C, Qi E R, Tian C and Wu Y S 2014 *Phys. Rev.* E **89** 042906
[50] Kuznetsov Yuri A 1998 *Elements of Applied Bifurcation Theory* 2nd edn (New York: Springer)

[51] Doedel E J and Oldeman B E AUTO-07P 2012 Continuation and bifurcation software for ordinary differential equations (http://indy.cs.concordia.ca/auto)

[52] Wolf A, Swift J B, Swinney H L and Vastano J A 1985 *Physica D* **16** 285-317

[53] Eckmann J P and Ruelle D 1985 *Rev. Mod. Phys.* **57** 617-56

[54] Strogatz S H 2015 *Nonlinear Dynamics and Chaos* (Boulder, CO: Westview Press)

[55] Nayfeh A H and Mook D T 1979 *Nonlinear Oscillations* (New York: Wiley)

[56] Fang X, Wen J, Yin J and Yu D 2017 *Nonlinear Dyn.* **87** 2677-95