Wave propagation in one-dimensional nonlinear acoustic metamaterials

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


Introduction
Like the phononic crystals [1,2], acoustic metamaterials [3][4][5][6] (AMs) are artificial medias 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, they have become a hot topic during the last decade [7]. By now, most studies [8][9][10][11][12] in this field are dealing with linear AMs (LAMs) which are based on a locally resonant (LR) mechanism [1] to yield negative indexes, super-lensing effects [13], or wave guiding. Different from the Bragg mechanism, LR mechanism paves the way to manipulate low-frequency waves with small "meta-atoms" [13]. However, the low-frequency LR bandgap is generally narrow [7] and heavy "atoms" are necessary to enlarge this bandgap.
The concept of nonlinear metamaterial was first introduced in electromagnetism to investigate the photonic metamaterials in 2003 [14]. According to particular needs, a nonlinear response can be constructed on purpose and even enhanced as compared to conventional nonlinearities found in natural materials, 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.
Compared with NAMs, studies on nonlinear periodic structures [17] (NPSs) have a longer history that dates back to 1983 when Nesterenko highlighted the soliton in a granular crystal [18]. Since then, most of the investigations in the field are dealing with granular crystals interacting in a nonlinear way through Hertzian contact [19,20]. Bistable lattices (another NPS) are studied recently [21]. Actually, NPSs still attract much attention and even they got an intense development since 2006 [19,[21][22][23]. The presence of nonlinear media in periodic structure can be utilized to realize interesting devices, such as acoustic diodes [24][25][26] and acoustic lenses [27]. New physical properties that are different from the linear phononic crystals are found, such as unidirectional transition [21], waves coupling [33], subharmonic frequencies [34,35] , discrete breathers [23], soliton waves [36,37] and surface waves [38]. Both simulations and experiments demonstrate that bifurcations may be observed in granular crystals [22,26,39]. Using a perturbation approach and the harmonic balance method, Narisetti et al. [40,41] recently studied the amplitude-dependent dispersions, stop band properties, and wave beaming in nonlinear periodic granular media. On their side, Manktelow et al. [42] investigated the wave propagation in layered NPS based on a perturbation approach. Furthermore, experimental works highlighted the role played by the critical amplitude in energy transmission [43] and bifurcation-induced bandgap reconfiguration [44] in one-dimensional (1D) NPSs. Recently, a 1D strong NPS has been proved to increase the velocity of sound and therefore the acoustic impedance [45]. However, while the granular crystals are suitable for ultrasonic applications, it is hard to consider them at 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 currently much attention. However, the mechanisms for both properties simultaneously happen are difficult to realize. LAMs consist of linear "meta-atom", but when the linear "meta-atom" gets nonlinear in NAM, properties of the wave propagation get different patterns. In the recent works [46] [47], 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 forbidden bands. This finding demonstrates that chaos is a novel and promising mechanism to achieve simultaneously low-frequency and broad band in both mono-bandgap NAMs and multi-bandgap NAMs and still that a strong nonlinearity is beneficial to expand the bandwidth by several times.
However, there are lots of phenomena arising in NAMs that have not yet been fully explained nor demonstrated. For example, why are the responses in the first pass band similar to ones 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 the nonlinearity has larger influence on the nonlinear LR bandgaps than it has for the Bragg gaps? Furthermore, some problems remain still unexplored, for example, the influence of the nonlinearity on the structure of the bands, the features of chaos and their differences, the jumps, the nonlinear resonances and their stabilities and the influence of the damping on chaos.
In this paper we attempt to answer these questions with the help of the frequency response analysis, the bifurcation theories, the Lyapunov exponents (LEs) and the fractal dimensions.

Diatomic model
The 1D damped diatomic nonlinear metamaterial 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 r with cubic nonlinear stiffness k 1 Δ+k 2 Δ 3 , where Δ denotes the relative displacement and k 1 (k 2 ) symbolizes the linear (nonlinear) stiffness. In addition, a linear viscous damping term c  is taken into consideration. Defining x n and y n as displacements of the linear and nonlinear oscillators in the nth cell, respectively; the differential equation for the nth cell reads where the dot denotes the derivation with respect to time t. The definitions of the generalized parameters are: mass ratio η=m r /m, stiffness ratio β 0 =k 0 /m= 2 s  , β 1 =k 1 /m r = 2 r  , β 2 =k 2 /m r . The generalized frequency is Ω=ω/ω s and the damping coefficient is ζ=c/m r . The strength factor of nonlinearity is defined as σ=3β 2  harmonic balance method or homotopy analysis method. The algorithms for these approaches are described in [47]. The dispersion solutions are eigen waves in the infinite metamaterial without external input.
We have considered a diatomic NAM with a finite length: 8 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.
Beside the diatomic model, tetratomic model is also discussed to figure out the influences on multiple bandgaps.

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 , r n denote the displacements of corresponding oscillators in the nth cell, respectively. The tetratomic NAM has more complicated dynamics.
The definitions of other parameters in Eq. (2) are identical with those in Eq. (1) The algorithms to analyze the bifurcation and characteristics of chaos in the tetratomic model are also identical with the diatomic model. The difference is the chain we considered here has only four cells, so the state space is also 32D. The parameters of the model are: m=1, m r =2, η=2, β 0 =225= 2 and Ω=ω/ω s . ω s is 3 times of ω rs , which is similar with the continuous metamaterials which requires soft resonators. A larger but still weak damping ζ=0.02 is adopted. Assuming that the driving amplitude is A 0 =0.01 and the nonlinear strength is σ=0.6, so it is strongly nonlinear system in this case.

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 ] , i.e., The element product of two vectors can be abbreviated as xy . Therefore, the notation for the mth power of x  A x y A y x . Moreover, it can be easily shown that the partial derivative of the element product is: These algebraic rules allow for a convenient analysis of the high-dimensional nonlinear equations and lead to compact formulas.
(2) Frequency responses Achieving the analytical solutions of Eqs. (1) and (2) for a high-dimensional system is not an easy task. However, the frequency responses can be approximated with the help of the numerical integral method, the averaging method, or the harmonic balance method (HBM). The latter can be implemented to find the approximate steady frequency responses. To this end, we define the coordinate transformation ˆn where M, C, K and N are mass, damping, stiffness and nonlinear coefficient matrices of the whole transformed system; f stands for the node force vector applied on every masses. Let us assume that the steady response has the form cos sin tt , where a and b are constant vectors. The first-order HBM procedure leads to a system of algebraic equations: The solutions of Eq. (7) would accurately describe the responses 22  Y = a b to the driving force, of the metamaterial. However, HBM method does not directly allow investigating the stability of the solutions. This is the reason why we adopted in this work the harmonic average approach (HAA) [48] [49] instead, to study the frequency responses. Within this approach, the solution is assumed to have the form: where t   . The derivatives with respect to time t of the formula in (8) are Comparing the expressions of y in Eqs. (8) and (9), we obtain ( cos sin ) The further substitution of (8) and (9) into (6) gives another form of the equation of motion: Then, assuming that u and v are constant, by calculating (11)×sinθ-(10)×cosθ and integrating the result from 0 to 2π, we obtain Similarly, calculating (11)×cosθ+(10)×sinθ and integrating over the interval [0, 2π] gives: The steady solutions correspond to the condition  u v 0 and their expressions are therefore: 6 The amplitude of the response is 22  The derivation of these solutions using HAA turns the steady response problem into an equilibrium problem of differential equations. Actually, the solutions in Eqs. (14) and (7) are equal. However, HAA allows for analyzing the stability of the solutions of both the systems of differential equations (12) and (13) via the Jacobian matrix. Omitting the terms in 2ω≠0 in the left sides of (12) and (13) since they will not influence the properties of the solutions, the Jacobian matrix J is derived by performing the derivations of (12) and (13) with respect to the vectors u and v, namely Therefore, with the solutions coming from Eq. (15), their stability can be determined: if the real part of an eigenvalue of 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 the NAMs and differential equations with cubic nonlinearity.
The pseudo-arclength continuation algorithm is further adopted to solve Eq. (14). Within this frame, Newton method (MATLAB fsolve) is used to find the solutions. Let us notice that if C=0, this formula reveals that v must read v=0.

Bifurcation analysis method in time domain
From now 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, the integrations of the equations of motion are made in the state space. Furthermore, the method allowing to derive the 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 i.e. it becomes a boundary value problem. Therefore, it becomes a 33D system in this situation.
When calculating the bifurcation diagrams, a differential system is added.
This system has a unique asymptotically stable solution, sin , cos tt . The boundary value at a specified time t 0 is defined as 0 . Therefore, 00 uA   and the whole system becomes 34D.
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 μ i . Let us recall here that, for a continuous differential system, there are three types of codimension-1 bifurcations: saddle-node (SN), period doubling (PD) and torus (TR) bifurcations [54]. The bifurcation associated with the appearance of μ i =1 (μ i =-1) is called a SN (PD) bifurcation; the bifurcation corresponding to the presence of conjugate multipliers μ 1,2 =exp(±iθ 0 ) , 0<θ 0 <π, is called a torus bifurcation. An invariant torus is a quasiperiodic 7 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: λ 1 , λ 2 ,…λ 33 (in descending order). λ 1 is named as the largest LE (LLE). All non-autonomous systems have at least one zero LE that corresponds to the t-component. Therefore, if λ 1 =0, we discard it and make λ n * =λ n+1 be the new spectrum. If λ 1 >0, the motion is Lyapunov chaotic, and if there are more than two positive λ i , the motion is described as Lyapunov hyperchaotic.
Let the quantity

Diatomic model
The bandgap properties of the diatomic NAM model have been laid out in Ref. [46]. It is shown that for this system, the generalized frequency range of the LR bandgap is  resonances, the stable steady solutions become unstable when the frequency exceeds a threshold. Their amplitudes increase then along a branch until a stable solution appears at a certain frequency, where the amplitudes of the stable solutions switch to a lower branch. Subsequently, the amplitude increases along this branch and a "neighbor peak" with lower amplitude germinates.
When the driving amplitude increases to A 0 =0.005, the stable and unstable solutions appear alternately. As a consequence the high amplitudes are unstable and the total frequency width of the unstable regions is broader than the width of the stable regions. However, the "neighbor peak" disappear so a curve of nonlinear resonance becomes a single trajectory bending to the right. These results illustrate the fact that, near a resonance, a high dimensional NAM has a behavior different from that of a single DoF or even of a 2DoF system [55]. 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 ζ=0.01. Therefore, a strong damping effect could allow a chaotic system to turn into a periodic system.
Hereafter, we consider a weak damping ζ=0.01 in the diatomic system.
We have chosen three illustrative frequencies Ω=0.678, 1.124, Ω=2.2, whose locations are labeled in figure   3a. Bifurcation diagrams and chaotic characteristics are illustrated in figure 5 for each of these frequencies.
For the Therefore, the damping ζ=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 5c. In the main, the damping reduces the number of quasi-periodic solutions in the system.
The unstable periodic solution is an essential condition for the chaos to appear. The bifurcation diagram, mean amplitude that directly comes out from numerical integration, LEs and d LD in figure 5 are well consistent with each other.
In the case Ω=0.678 the first unstable solution appears at A 0 =1.867×10 -3 marked by the vertical black line in figure 5. Within the subsequent domain, the stable and unstable periodic solutions appear alternately. When A 0 >6.686×10 -3 , the periodic solution enters into 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×10 -3 <A 0 <1.867×10 -3 , 1   →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 still increases. Actually, both 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, this is a general rule that the system is hyperchaotic when chaos occurs. Moreover, the variation laws of d LD are identical with λ 1 in the chaotic domain, and when strong chaos occurs d LD →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 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 period-doubling bifurcations. The waterfall plot of power spectra illustrates this process, as shown in figure 6. When A 0 →0, we are dealing with a 1:1 resonance. In this situation, the energy gets localized at the driving frequency Ω e . Further increasing A 0 leads to odd harmonic waves 3Ω e and then 5Ω e to appear, so 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.  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 keep 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 the bandgaps via amplitudes, the system can switch suddenly between the stop state (i.e. the bounded state) and the propagation states (i.e. the excited quantum states). These characteristics may be helpful to realize acoustic devices with small dimensions, as for example acoustic diodes or switches. In figure 5b, 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 the jump back to the highest orbit.
The case Ω=2.2 is similar to the case Ω=1.124 except that the motion turns from periodic into chaotic for only one jump at the first unstable periodic solution. The trend in bifurcation diagram agrees well with the mean amplitude variations. After the jump, the amplitude at first slowly increases and then keeps constant as do both LEs and d LD . This suggests a mechanism to suppress elastic waves in strongly NAM: strengthening the nonlinearity decreases the wave transmissibility. Based on all these results, the band structure of a strongly diatomic NAM is sketched in figure 9.
13 Strong hyperchaotic pass band ω Linear (resonant) pass band Transmissibility LR Linear Quasi-linear chaos Figure 9. The sketch of the band structure of diatomic NAM model.

Tetratomic model
The dispersion relations of this model were elaborated in Ref. [47].   In the case Ω=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 further increases. However, 1   <0.03 and d LD <<33, so that 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 → weakly chaotic → strongly hyper-chaotic". However, in these three passbands, the LLEs λ 1 for the strong chaos are 10 times larger than that calculated for the weak chaos in the first passband and d LD →33. The mechanisms that suppress resonances in the passbands are identical to the ones observed with the diatomic NAM: chaos induced by period-doubling bifurcation; the wave amplitudes along the chaotic branches slowly increase, remain constant or even decrease with A 0 . Moreover, as illustrated in figure 11(b,c), periodic windows in the chaotic region are observed, which means that LEs and LDs do not monotonously vary.
Close to Ω=1.42, T A is equal to the transmission coefficient for the corresponding LAM. In fact, there are other specific frequency domains. The comparison with the results obtained with the 10-cells model [47] indicates that these domains relate to the density of the linear modes that is determined by the length of chain.
In the 4-cells model, this domain is a non-resonant region. The periodic solutions for Ω=1. 42  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 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 homogeneously amplification in a broad passband. Another question we have addressed is: why is the nonlinear LR bandgap less efficient than other bandgaps to suppress the elastic waves in NAM? To answer this question, we have set the frequencies to Ω=0.35 and Ω=1. As shown in figure 12, as for the diatomic model, there are multiple bifurcations and branches. However, the stable periodic solutions disappear for an amplitude as low as A 0 =0.003 because of the strong nonlinearity and therefore this bandgap closes at a certain driving amplitude and is replaced by a high energy chaotic orbit, as it is the case in the diatomic model. In contrast, for the BG bandgaps away from this nonlinear bandgap, much higher amplitudes are needed to get the unstable periodic solution. In the example illustrated in figure   11d, when Ω=1 in the first BG bandgap, a periodic motion is observed and the chaos is weak ( 1   <0.04) if the amplitude A 0 <0.01, as for the quasi-linear state in the interval A 0 =0.01~0.015. It is only as A 0 >0.015, that the motion is attracted to the high energy hyper-chaotic state. Combing the mechanisms in diatomic model, one concludes that the nonlinear LR bandgap has weaker ability to reflect the incident waves when it jumps to the high-energy orbits.

Conclusions
The propagation of waves in the nonlinear acoustic metamaterials (NAMs) is fundamentally different from that in the conventional LAMs. However, those 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 broad frequency ranges. The nonlinear wave behaviors, band structures, the bifurcations and the chaos are studied to demonstrate the novel mechanisms that can manipulate wave propagations.
HAA combined with the pseudo-arclength 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 LD ). Our results show that the nonlinear resonances feature multiple branches with unstable peaks and a "neighbor peak". We further have demonstrated that the period doubling process leads the motion from being periodic to chaos and damping has a significant influence on the characteristics of chaotic attractors. Moreover, due to dispersion process in 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 the amplitude so 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, so the responses in these regions are similar to the ones of LAMs. Actually, the hardened NAM influences all the passbands that are higher than the nonlinear LR bandgap, including a small region below this bandgap.
In the nonlinear LR bandgap, we put into evidence "quantum" behaviors because of jumping bifurcations between low-energy and high-energy orbits, whose propagation states (excited states) and stop state (bounded state) have discrete characteristics and switch suddenly. This behavior also explained why the nonlinear LR gets weaker ability to reflect the incident waves. Moreover, jumps in passbands would amplify non-resonant waves.
This work provides an important theoretical base for the understandings and applications of NAMs.